# Load R package
library(eDNAoccupancy)

# Input variables
eDNA_detection_file="Goby-Data.csv"
metadata_file="Goby-Metadata.csv"
site_covariates = "~ veg"
sample_covariates = "~ sal+fish+turb"
site_sample_covariates ="~ sal+twg"

# Output files
posterior_summary_file = "posterior_summary.txt"
posterior_median_estimates_file = "posterior_median_estimates.txt"
pdf_file ="plot.pdf"

# Function to plot the relationship between a covariate and the probability of eDNA detection in samples.
plot_covariate= function(metadata,covariate)
{
  vals = metadata[, covariate]
  plot(vals, p$median[,1], ylim=ylimits, cex=3, main=covariate, xlab=covariate, cex.lab=1.5)
  segments(vals, p$lower[,1], vals, p$upper[,1], lwd=2)
}

# Read input data and meta data.
detection_data = read.table(eDNA_detection_file, header=TRUE, sep="," )
metadata = read.table(metadata_file, header=TRUE, sep="," )

# Compute occupancy data matrices
species_detections = occData(detection_data, siteColName = 'site',
                             sampleColName = 'sample')

# Center and scale numeric-valued covariate measurements
metadata.sc = scaleData(metadata)

# Fit the model.
fit = occModel(formulaSite = site_covariates, formulaSiteAndSample = site_sample_covariates, formulaReplicate = sample_covariates ,
               detectionMats = species_detections, siteData = metadata.sc, siteColName = 'site')

# Print Posterior Summaries
posterior.summary=as.data.frame(posteriorSummary(fit, burnin=1000, mcError=TRUE,outputSummary = TRUE))
posterior.summary.df = data.frame("params"=rownames(posterior.summary),posterior.summary)

# Write posterior summaries to file.
write.table(posterior.summary.df,file = posterior_summary_file, sep="\t", row.names=FALSE, quote=FALSE)

# Estimate the posterior summaries of derived parameters of the models - psi, theta and p
# psi = Probability of occurrence of a species in a location 
# theta = Conditional probability of occurrence of a species in a sample given that species is present in the location where the sample is collected. 
# p = Conditional probability of species detection in a subsample of a sample given that species is present in the sample. 

psi = posteriorSummaryOfSiteOccupancy(fit, burnin=1000)
theta = posteriorSummaryOfSampleOccupancy(fit, burnin=1000)
p = posteriorSummaryOfDetection(fit, burnin=1000)

# output estimates of posterior medians
post.medians=cbind(psi=psi$median, theta=theta$median[,1], p=p$median[,1])

# Write estimates of posterior medians to a file.
post.medians.df = data.frame("site"=rownames(post.medians), post.medians)
write.table(post.medians.df,file = posterior_median_estimates_file, sep="\t", row.names=FALSE, quote=FALSE)

# Extract numeric ovariates to plot
covariates_to_plot =metadata[sapply(metadata,is.double)]
col_names = colnames(covariates_to_plot)
ncols= ncol(covariates_to_plot)

# Plot the estimated relationship between each covariate (eg: sal, fish, or turb) 
# and the probability of detection of eDNA in PCR replicates 

pdf(pdf_file, width=10, height=10)
ylimits = c(0.2,1)
par(mfrow=c(1,ncols), mar = c(5, 5, 5, 5), oma=c(0,0,3,0))

for (col in col_names) {
  plot_covariate(metadata,col)
}

mtext("Estimated relationship between covariate and the probability of eDNA detection", outer=TRUE, cex=1.5)
dev.off()
unlink("mc.*")