# Psilocybin Research Database - Linear publication growth since 2020
# Source: https://psilocybin-research.com/msca-impact/
# Complete years through 2025; visible, non-false-positive catalogue records.
# Snapshot SHA-256: 756f970ec11384e2b82ac95f6f7a042c56b3eab80cd8b042a7edf2f0299f1415
# Includes the annual counts: no account, database download or credentials needed.
# Run this entire file in RStudio, source("psilocybin-linear-growth.R"), or Rscript psilocybin-linear-growth.R.
# The default snapshot reproduces this download. To use current database counts,
# change "snapshot" to "live" below (internet required). Results may then change.

data_source <- Sys.getenv("PSILOCYBIN_DATA_SOURCE", "snapshot")
snapshot_data <- read.csv(text = paste(c(
  "year,count", "1958,5", "1959,11", "1960,11", "1961,14", "1962,15",
  "1963,20", "1964,10", "1965,10", "1966,9", "1967,12", "1968,15",
  "1969,13", "1970,12", "1971,9", "1972,3", "1973,15", "1974,15",
  "1975,9", "1976,7", "1977,9", "1978,10", "1979,7", "1980,9",
  "1981,16", "1982,8", "1983,4", "1984,7", "1985,6", "1986,8",
  "1987,6", "1988,6", "1989,6", "1990,7", "1991,1", "1992,4",
  "1993,6", "1994,6", "1995,1", "1996,5", "1997,4", "1998,15",
  "1999,18", "2000,17", "2001,8", "2002,17", "2003,18", "2004,12",
  "2005,20", "2006,25", "2007,16", "2008,24", "2009,24", "2010,16",
  "2011,22", "2012,24", "2013,37", "2014,39", "2015,41", "2016,65",
  "2017,73", "2018,106", "2019,114", "2020,170", "2021,304", "2022,423",
  "2023,550", "2024,653", "2025,885"
), collapse = "\n"), colClasses = c("integer", "integer"))
last_complete_year <- 2025L
if (identical(data_source, "snapshot")) {
  annual_counts <- snapshot_data
} else if (identical(data_source, "live")) {
  csv_file <- tempfile(fileext = ".csv")
  tryCatch({
    download.file("https://psilocybin-research.com/msca-impact/analysis-code.php?format=csv", csv_file, method = "libcurl", quiet = TRUE,
                  headers = c("User-Agent" = "Mozilla/5.0 (compatible; PsilocybinResearch-R/1.0)"))
    annual_counts <- read.csv(csv_file, colClasses = c("integer", "integer"))
  }, error = function(e) {
    stop("Could not download annual counts. Use data_source = 'snapshot' to run offline.\n", conditionMessage(e))
  }, finally = unlink(csv_file))
  last_complete_year <- as.integer(format(Sys.time(), "%Y", tz = "UTC")) - 1L
} else {
  stop("data_source must be 'snapshot' or 'live'.")
}
stopifnot(all(c("year", "count") %in% names(annual_counts)))
annual_counts <- annual_counts[complete.cases(annual_counts) &
                               annual_counts$year <= last_complete_year & annual_counts$count > 0L, ]
annual_counts <- annual_counts[order(annual_counts$year), ]
if (anyDuplicated(annual_counts$year)) stop("Annual counts must have unique years.")
cat("Source: Psilocybin Research Database\n")
cat("Data:", data_source, "| complete years through", last_complete_year, "\n")
if (identical(data_source, "snapshot")) cat("Snapshot SHA-256: 756f970ec11384e2b82ac95f6f7a042c56b3eab80cd8b042a7edf2f0299f1415\n")

# Descriptive linear growth, complete years since 2020.
# This reports the same coefficients and raw-count R-squared as the website.
# It does not treat ordinary regression coefficient p-values as time-series tests.

d <- subset(annual_counts, year >= 2020L)
if (nrow(d) < 3L) stop("At least three complete years since 2020 are required.")
d$t <- d$year - 2020L
model <- lm(count ~ t, data = d)
analysis_result <- list(
  first_year = min(d$year), last_year = max(d$year), observations = nrow(d),
  intercept = unname(coef(model)[1L]), annual_increase = unname(coef(model)[2L]),
  r_squared = summary(model)$r.squared
)
print(analysis_result)
if (interactive()) {
  plot(d$year, d$count, xlab = "Publication year", ylab = "Catalogue records",
       main = "Descriptive linear growth", pch = 19, col = "#476f91")
  lines(d$year, fitted(model), col = "#c8943f", lwd = 2)
}