analysis.R

library(R2OpenBUGS)
library(coda) # Load the coda package

# Get the data from Python
Y <- data$Y
pvID <- data$pvID
cropID <- data$cropID
blockID <- data$blockID
mainplotID <- data$mainplotID
N <- length(Y)
nPV <- length(unique(pvID))
nCrops <- length(unique(cropID))
nBlocks <- length(unique(blockID))
nMainPlots <- length(unique(mainplotID))

# BUGS model data
bugs_data <- list("N", "Y", "pvID", "cropID", "blockID", "mainplotID",
                  "nPV", "nCrops", "nBlocks", "nMainPlots")

# BUGS model parameters to monitor
bugs_params <- c("alpha", "beta.pv", "beta.crop", "beta.inter",
                 "sigma.within", "sigma.block", "sigma.mainplot")

# BUGS model initial values
inits <- function() {
  list(
    alpha = rnorm(1, 0, 1),
    beta.pv = c(NA, rnorm(nPV - 1, 0, 1)),
    beta.crop = c(NA, rnorm(nCrops - 1, 0, 1)),
    beta.inter = matrix(rnorm(nPV * nCrops, 0, 1), nrow = nPV, ncol = nCrops),
    tau.within = rgamma(1, 1, 1),
    tau.block = rgamma(1, 1, 1),
    tau.mainplot = rgamma(1, 1, 1)
  )
}

# Create a directory for CODA output if it doesn't exist
if (!dir.exists("coda_output")) {
  dir.create("coda_output")
}

# Run OpenBUGS, requesting CODA-compatible output
# The result 'bugs_coda' will be an mcmc.list object
bugs_coda <- bugs(data = bugs_data,
                    inits = inits,
                    parameters.to.save = bugs_params,
                    model.file = "agrivoltaic_lognormal.txt",
                    n.chains = 3,
                    n.iter = 2000,
                    n.burnin = 1000,
                    n.thin = 1,
                    codaPkg = TRUE, # Return an mcmc.list object
                    debug = TRUE,
                    OpenBUGS.pgm = "C:\\OpenBUGS323\\OpenBUGS.exe")

# --- CODA Output Analysis ---

# 1. Save the summary statistics and quantiles to text files
summary_stats <- summary(bugs_coda)
write.table(summary_stats$statistics, file = "coda_summary_stats.txt", sep = "\t", col.names = NA)
write.table(summary_stats$quantiles, file = "coda_summary_quantiles.txt", sep = "\t", col.names = NA)

# 2. Save the raw MCMC chains to CODA-formatted text files
# This creates coda_output/chain.INDEX and coda_output/chain1.txt, chain2.txt, etc.
write.coda(bugs_coda, stem = "coda_output/chain")

# 3. Generate and save diagnostic plots to a PDF
pdf("coda_diagnostic_plots.pdf")
plot(bugs_coda)
dev.off()

# 4. Print the summary to the R console
print(summary_stats)