# This code can be used to generate the simulated datasets and BLR, KMR, and LM parameter estimates, as well as PSIS-LOO statistics to compare BLR and LM for the simulation study performed using parameter values in Table S1 in Vitense et al (in press).

library(R2jags)
library(R2WinBUGS)
library(mcmcplots)
library(ggplot2)
library(gridExtra)
library(KernSmooth)
library(loo)


source("SimDatAdd_2sig.R") # source function that creates simulated datasets


# The parameters below set the number of simulations, sample size per simulation, Chla perturbation strength, SAV measurement error, and are used to name the directory where output is saved

nsims <- 1 # number of simulated datasets to create and analyze (at each noise level if doing multiple)
n <- 100 # number of lakes (sample size) in each simulation (change as needed)
SAV.sigma <- .5 # SAV error (.5 was used in the paper, .2 (low SAV error) was used for comparison in Figure S10 in Appendix S1)


sigma.mat <- matrix(c(rep(.05, 2), rep(.125, 2), rep(.2, 2)), ncol=2, byrow=T) # standard deviation for perturbations to Chla. Rows are different noise levels for each set of simulations (low, medium, and high in the paper). First column is noise parameter for lakes in the clear state, second column is noise parameter for lakes in the turbid state (which are set to be equal in the paper).


set.seed(13)
seeds <- sample.int(1000000000,size=1000, replace=FALSE) # sample 1000 seeds to loop through for reproducibility. The same seeds are always sampled by setting the seed to 13 (or another number) first.






######### BLR MODEL #########

latent_reg_blr <-function(){
  
  # Priors 
  a0 ~  dnorm(0,0.01) # N(0,10^2)
  a1adj ~  dnorm(0,0.1)  
  b0 ~ dunif(0,10)   # slopes should be positive
  b1 ~ dunif(0,10)  # N(0,3.16^2)
  alpha ~ dnorm(0,.01) # N(0,10^2) -- intercept for logistic regression (with SAV)
  negbeta ~ dlnorm( .5, 1/1^2 ) # slope for logistic regression (with SAV)
  sigma1 ~ dunif(0, 20)  # Residual standard deviation
  sigma2 ~ dunif(0, 20)  # Residual standard deviation
  tau1 <- 1 / ( sigma1 * sigma1)
  tau2 <- 1 / ( sigma2 * sigma2)
  
  beta <- -1*negbeta  # make beta negative -- gives a mean of -e (-2.72)
  
  b2 <- (a1+b1*logp1-a0-b0*logp2) /(logp1-logp2) # derived slope of unstable line
  constraint <- step(b2)
  
  
  logp1 ~ dunif(0,6.5)   # left bifurcation point
  logp2 ~ dunif(0,6.5)   # right bifurcation point
  p1 <- exp(logp1)
  p2 <- exp(logp2)
  b1adj <- b1-b0 # difference in slopes
  a1 <- a1adj+a0 # difference in intercepts
  
  
  
  for(i in 1:nsamp){
    lnC[i] ~ dnorm(mu[i], tau1*(1-S[i]) + tau2*S[i])
    mu[i] <- a0+a1adj*S[i] + b0*lnP[i]*(1-S[i])+b1*lnP[i]*S[i] # linear relationship between log(TP) and log(Chla). a1 and b1 are adjustments to intercept and slope for turbid state
    S[i] ~ dbern(p[i])
    logit(ptmp[i])<-alpha+beta*SAVc[i]         
    p[i] <- ptmp[i]*step(lnP[i]-logp1)*step(logp2-lnP[i])+step(lnP[i]-logp2) # sets p to 0 if TP less than p1, p depends on SAV if between p1 and p2, and p is 1 if TP greater than p2
    ForceNeg[i] ~ dbern(constraint)
    
    log_lik[i] <- log( dnorm(lnC[i], mu[i], tau1*(1-S[i]) + tau2*S[i]) )
    
  }
}

# Initial values
inits.blr <- function(){ list(a0=-2.17, a1adj=3.274869, logp1=4.23, logp2=5.71, b0=.803, b1=.535, alpha=runif(1,-4,-1), negbeta=runif(1,5,10), sigma1=runif(1, .4,.6), sigma2=runif(1, .4,.6) )}


# Parameters to estimate
params.blr <- c("a0","a1","a1adj", "b0", "b1","b1adj","b2","alpha","beta","p","p1","p2","logp1","logp2", "sigma1", "sigma2", "S", "log_lik") 

# MCMC settings
nc.blr <- 4 # number of chains
ni.blr <- 50000 # total number of iterations  
nb.blr <- 20000 # burn-in   
nt.blr <- 12     # thinning rate

###### END OF BLR MODEL #####






########### LINEAR MODEL #############
latent_reg_lm <-function(){
  
  # Priors 
  a0 ~  dnorm(0,0.01) # N(0,10^2)
  b0 ~ dunif(0,10)   # slopes should be positive
  sigma1 ~ dunif(0, 20)  # Residual standard deviation
  tau1 <- 1 / ( sigma1 * sigma1)
  
  for(i in 1:nsamp){
    lnC[i] ~ dnorm(mu[i], tau1)
    mu[i] <- a0 + b0*lnP[i] # linear relationship between log(TP) and log(Chla). 
    
    log_lik[i] <- log( dnorm(lnC[i], mu[i], tau1) )
  }
}


# Inits function
inits.lm <- function(){ list(a0=rnorm(1, 0, 1),  b0=runif(1, 0, 2),  sigma1=runif(1, .2, 1))} 

# Parameters to estimate
params.lm <- c("a0", "b0",  "sigma1", "log_lik") 

# MCMC settings - NOT NEEDED WHEN USING SAME SETTINGS AS BLR
# nc.lm <- 4
# ni.lm <- 40000  
# nb.lm <- 20000    
# nt.lm <- 8     

################## END OF LINEAR MODEL ############






### Do simulation

for(j in 1:nrow(sigma.mat)){
  
sigma <- sigma.mat[j,] # pull out one row (one set of noise parameters)
  
  
dir.create(paste0("N",n,"_nsims",nsims, "_SAV",SAV.sigma,"_sigma",sigma[1]), recursive=TRUE, showWarnings = FALSE) # Name directory using sample size, number of simulations, perturbation strength, and SAV measurement error

setwd(paste0("N",n,"_nsims",nsims, "_SAV",SAV.sigma,"_sigma",sigma[1]))


##### Create data frames for storing output

# BLR data
sigma1.dat <- sigma2.dat <- a0.dat <- a1.dat <- a1adj.dat <- b0.dat <- b1.dat <- b1adj.dat <- b2.dat <- alpha.dat <- beta.dat <- logp1.dat <- logp2.dat <- p1.dat <- p2.dat <- data.frame(mean=rep(NA,nsims),sd=rep(NA,nsims),Q2.5=rep(NA,nsims),Q25=rep(NA,nsims),Q50=rep(NA,nsims),Q75=rep(NA,nsims),Q97.5=rep(NA,nsims),Rhat=rep(NA,nsims),n.eff=rep(NA,nsims),min=rep(NA,nsims),max=rep(NA,nsims),mode.d=rep(NA,nsims),mode.dSJ=rep(NA,nsims),mode.dKS=rep(NA,nsims))

# LM data
a0LM.dat <- b0LM.dat <- sigma1LM.dat <- data.frame(mean=rep(NA,nsims),sd=rep(NA,nsims),Q2.5=rep(NA,nsims),Q25=rep(NA,nsims),Q50=rep(NA,nsims),Q75=rep(NA,nsims),Q97.5=rep(NA,nsims),Rhat=rep(NA,nsims),n.eff=rep(NA,nsims),min=rep(NA,nsims),max=rep(NA,nsims) )

# Alpha estimate for non-centered SAV + similarity between BLR and KMR, BLR and True states, KMR and true states
alpham.adj <- Percent.Match.LatentKmeans <- Percent.Match.LatentTrue <- Percent.Match.KmeansTrue <- rep(NA,nsims)

# KMR data
a0.kmeans <- a1.kmeans <- b0.kmeans <- b1.kmeans <- data.frame(est=rep(NA,nsims),se=rep(NA,nsims))
#pvalues.kmeans <- data.frame(pInt=rep(NA,nsims), pSlope=rep(NA,nsims))
sigma1.kmeans <- sigma2.kmeans <- rep(NA,nsims)
p1.kmeans <- data.frame(min=rep(NA,nsims),Q2.5=rep(NA,nsims),Q5=rep(NA,nsims))
p2.kmeans <- data.frame(max=rep(NA,nsims),Q97.5=rep(NA,nsims),Q95=rep(NA,nsims))

# PSIS-LOO data
loo.dat <- data.frame(elpd_loo_BLR=rep(NA,nsims),se_elpd_loo_BLR=rep(NA,nsims),p_loo_BLR=rep(NA,nsims),se_p_loo_BLR=rep(NA,nsims),looic_BLR=rep(NA,nsims),se_looic_BLR=rep(NA,nsims),NumBad_loo_BLR=rep(NA,nsims),PropBad_loo_BLR=rep(NA,nsims), elpd_loo_LM=rep(NA,nsims), se_elpd_loo_LM=rep(NA,nsims), p_loo_LM=rep(NA,nsims), se_p_loo_LM=rep(NA,nsims), looic_LM=rep(NA,nsims), se_looic_LM=rep(NA,nsims), NumBad_loo_LM=rep(NA,nsims), PropBad_loo_LM=rep(NA,nsims), elpd_diff_loo=rep(NA,nsims), se_elpd_diff_loo=rep(NA,nsims))


######################




for (i in 1:nsims){
  print(i)
  set.seed(seeds[i])
  
  # Simulate data
  mulogTP <- 4.5
  sigmalogTP <- 1
  # Simulate data
  LakeSim <- run.lake.sim.ADD(alpha=14, p=8, r=.63, s=.17, pi=270, c=1/200,mulogTP=mulogTP,sigmalogTP=sigmalogTP,n=n, sigma=sigma,dt=.001,numdays=30,plotpoints=F, TPdist="lnorm")
  
  #### Add sampling error to SAV (not process error)
  SAV1.err <- LakeSim$SAV1+rnorm(length(LakeSim$SAV1),0,SAV.sigma)
  SAV1.err <- SAV1.err/max(SAV1.err); SAV1.err[SAV1.err<0] <- 0
  
  
  # Pull out state variables
  TPug <- LakeSim$TP  # ug/L
  Chla <- LakeSim$Chla1  # ug/L
  lnChla <- log(Chla)
  SAVms <- SAV1.err      
  lnSAVms <- log(SAVms+.02)
  lnTPug <- log(TPug)
  
  
  #### Rename variables to match old code:
  lnP <- lnTPug
  P <- TPug
  lnC <- lnChla
  Chl <- Chla
  SAV <- SAVms
  SAV.mean <- mean(SAV)
  SAVc <- SAV-SAV.mean # centered SAV
  lnSAV <- lnSAVms
  
  
######## Do K-Means Cluster with Regression Analysis #########
  clust.dat <- data.frame(lnSAV,lnC)
  fit <- kmeans(clust.dat, 2, algorithm="Hartigan-Wong", nstart=50) # 2 cluster solution
  
  #### Append new cluster assignment:
  # Ensure coding is 0=clear, 1=turbid
  if(fit$centers[2,2] > fit$centers[1,2]){ fit.cluster <- as.factor(fit$cluster-1) }
  if(fit$centers[2,2] < fit$centers[1,2]){
    fit.cluster <- fit$cluster; fit.cluster[fit.cluster==2] <- 0
    fit.cluster <- as.factor(fit.cluster)
  }
  fit.dat <- data.frame(clust.dat, fit.cluster=fit.cluster, lnP=lnP)
  fit.dat0 <- subset(fit.dat, fit.cluster=="0")
  fit.dat1 <- subset(fit.dat, fit.cluster=="1")
  
  #### Regress log(Chla) on log(TP) in each cluster
  mod0 <- lm(lnC~lnP, data=fit.dat0)
  mod1 <- lm(lnC~lnP, data=fit.dat1)
  
  a0.kmeans[i,] <- coef(summary(mod0))[1,1:2]
  a1.kmeans[i,] <- coef(summary(mod1))[1,1:2]
  
  b0.kmeans[i,] <- coef(summary(mod0))[2,1:2]
  b1.kmeans[i,] <- coef(summary(mod1))[2,1:2]
  
  sigma1.kmeans[i] <- summary(mod0)$sigma
  sigma2.kmeans[i] <- summary(mod1)$sigma
  
  p1.kmeans[i,1] <- min(fit.dat1$lnP)
  p1.kmeans[i,2] <- quantile(fit.dat1$lnP, probs=.025)[[1]]
  p1.kmeans[i,3] <- quantile(fit.dat1$lnP, probs=.05)[[1]]
  
  p2.kmeans[i,1] <- max(fit.dat0$lnP)
  p2.kmeans[i,2] <- quantile(fit.dat0$lnP, probs=.975)[[1]]
  p2.kmeans[i,3] <- quantile(fit.dat0$lnP, probs=.95)[[1]]
  
######### End of k-means ###########
  
  
  
  
  
#### Run BLR and LM JAGS routines:
ForceNeg <- rep(0,length(lnC))
nsamp <- length(lnC)
  
win.data.blr <- list("lnC", "lnP", "nsamp", "SAVc", "ForceNeg")
win.data.lm <- list("lnC", "lnP", "nsamp")
  
  
##### Start MCMC chains ####
# BLR
out.lat.centerSIM <- do.call(jags.parallel, list(win.data.blr, inits.blr, params.blr, latent_reg_blr, nc.blr, ni.blr, nb.blr, nt.blr))

# LM  
out.lat.centerSIM.lin <- do.call(jags.parallel, list(win.data.lm, inits.lm, params.lm, latent_reg_lm, nc.blr, ni.blr, nb.blr, nt.blr))
#out.lat.centerSIM.lin <- do.call(jags.parallel, list(win.data.lm, inits.lm, params.lm, latent_reg_lm, nc.lm, ni.lm, nb.lm, nt.lm)) # use this one for setting different MCMC parameters for LM
  



##### PSIS-LOO model selection statistics for BLR vs. LM ####

llik <- out.lat.centerSIM$BUGSoutput$sims.list$log_lik
llik.lin <- out.lat.centerSIM.lin$BUGSoutput$sims.list$log_lik

loo <- loo(llik) # PSIS-LOO for BLR

loo.dat[i,1] <- loo$elpd_loo    # elpd_loo_BLR = Expected log pointwise predictive density
loo.dat[i,2] <- loo$se_elpd_loo # se_elpd_loo_BLR = Standard error of Expected log pointwise predictive density
loo.dat[i,3] <- loo$p_loo     # p_loo_BLR = Estimated effective number of parameters
loo.dat[i,4] <- loo$se_p_loo  # se_p_loo_BLR = Standard error of effective number of parameters
loo.dat[i,5] <- loo$looic     # looic_BLR = The LOO information criterion (-2*elpd_loo, i.e., converted to deviance scale)
loo.dat[i,6] <- loo$se_looic  # se_looic_BLR = Standard error of LOO information criterion
loo.dat[i,7] <- sum(loo$pareto_k > .7) # NumBad_loo_BLR = number points approximated poorly (pareto k values greater than 0.7)
loo.dat[i,8] <- sum(loo$pareto_k > .7)/n # PropBad_loo_BLR = proportion points approximated poorly


loo.lin <- loo(llik.lin) # PSIS-LOO for LM

loo.dat[i,9] <- loo.lin$elpd_loo    # elpd_loo_LM = Expected log pointwise predictive density
loo.dat[i,10] <- loo.lin$se_elpd_loo # se_elpd_loo_LM = Standard error of Expected log pointwise predictive density
loo.dat[i,11] <- loo.lin$p_loo     # p_loo_LM = Estimated effective number of parameters
loo.dat[i,12] <- loo.lin$se_p_loo  # se_p_loo_LM = Standard error of effective number of parameters
loo.dat[i,13] <- loo.lin$looic     # looic_LM = The LOO information criterion (-2*elpd_loo, i.e., converted to deviance scale)
loo.dat[i,14] <- loo.lin$se_looic  # se_looic_LM = Standard error of LOO information criterion
loo.dat[i,15] <- sum(loo.lin$pareto_k > .7) # NumBad_loo_LM = number points approximated poorly  (pareto k values greater than 0.7)
loo.dat[i,16] <- sum(loo.lin$pareto_k > .7)/n # PropBad_loo_LM = proportion points approximated poorly

loo_diff <- compare(loo, loo.lin) # second minus first (higher one is better)
loo.dat[i,17] <- loo_diff[[1]]  # elpd_diff_loo = difference in expected predictive accuracy (difference in elpd_loo)
loo.dat[i,18] <- loo_diff[[2]]  # se_elpd_diff_loo = standard error of elpd_diff
# The difference will be positive if the expected predictive accuracy for the second model is higher. So NEGATIVE values mean BLR is better

# CAN USE POINT ESTIMATE AND SE OF DIFFERENCE TO DETERMINE SIGNIFICANCE (see how this was done in BLR_applied_to_MDNRdata.Rmd)





#### sAVE COEFS FROM LM ####

a0LM.row <- which(row.names(out.lat.centerSIM.lin$BUGSoutput$summary)=="a0")
a0LM.dat[i,1:9] <- out.lat.centerSIM.lin$BUGSoutput$summary[a0LM.row,]
a0LM.dat[i,10] <- min(out.lat.centerSIM.lin$BUGSoutput$sims.list$a0)
a0LM.dat[i,11] <- max(out.lat.centerSIM.lin$BUGSoutput$sims.list$a0)

b0LM.row <- which(row.names(out.lat.centerSIM.lin$BUGSoutput$summary)=="b0")
b0LM.dat[i,1:9] <- out.lat.centerSIM.lin$BUGSoutput$summary[b0LM.row,]
b0LM.dat[i,10] <- min(out.lat.centerSIM.lin$BUGSoutput$sims.list$b0)
b0LM.dat[i,11] <- max(out.lat.centerSIM.lin$BUGSoutput$sims.list$b0)

sigma1LM.row <- which(row.names(out.lat.centerSIM.lin$BUGSoutput$summary)=="sigma1")
sigma1LM.dat[i,1:9] <- out.lat.centerSIM.lin$BUGSoutput$summary[sigma1LM.row,]
sigma1LM.dat[i,10] <- min(out.lat.centerSIM.lin$BUGSoutput$sims.list$sigma1)
sigma1LM.dat[i,11] <- max(out.lat.centerSIM.lin$BUGSoutput$sims.list$sigma1)








#### Save summary statistics for BLR posterior distributions ####
  
  # n.eff = an estimate of the number of independent samples from the posterior that would hold the same information
  a0.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="a0")
  a1.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="a1")
  a1adj.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="a1adj")
  b0.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="b0")
  b1.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="b1")
  b1adj.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="b1adj")
  b2.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="b2")
  alpha.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="alpha")
  beta.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="beta")
  logp1.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="logp1")
  logp2.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="logp2")
  p1.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="p1")
  p2.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="p2")
  sigma1.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="sigma1")
  sigma2.row <- which(row.names(out.lat.centerSIM$BUGSoutput$summary)=="sigma2")
  
  sigma1.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[sigma1.row,]
  sigma1.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$sigma1)
  sigma1.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$sigma1)
  
  sigma2.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[sigma2.row,]
  sigma2.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$sigma2)
  sigma2.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$sigma2)

  a0.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[a0.row,]
  a0.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$a0)
  a0.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$a0)
  
  a1.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[a1.row,]
  a1.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$a1)
  a1.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$a1)
  
  a1adj.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[a1adj.row,]
  a1adj.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$a1adj)
  a1adj.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$a1adj)
  
  b0.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[b0.row,]
  b0.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$b0)
  b0.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$b0)
  
  b1.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[b1.row,]
  b1.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$b1)
  b1.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$b1)
  
  b1adj.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[b1adj.row,]
  b1adj.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$b1adj)
  b1adj.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$b1adj)
  
  b2.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[b2.row,]
  b2.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$b2)
  b2.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$b2)
  
  alpha.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[alpha.row,]
  alpha.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$alpha)
  alpha.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$alpha)
  
  beta.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[beta.row,]
  beta.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$beta)
  beta.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$beta)
  
  logp1.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[logp1.row,]
  logp1.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$logp1)
  logp1.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$logp1)
  
  logp2.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[logp2.row,]
  logp2.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$logp2)
  logp2.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$logp2)
  
  p1.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[p1.row,]
  p1.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$p1)
  p1.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$p1)
  
  p2.dat[i,1:9] <- out.lat.centerSIM$BUGSoutput$summary[p2.row,]
  p2.dat[i,10] <- min(out.lat.centerSIM$BUGSoutput$sims.list$p2)
  p2.dat[i,11] <- max(out.lat.centerSIM$BUGSoutput$sims.list$p2)
  
  alpham.adj[i] <- out.lat.centerSIM$BUGSoutput$mean$alpha-SAV.mean*out.lat.centerSIM$BUGSoutput$mean$beta
  
  
########### KDE estimate for modes columns 12-14 ##########

### Note: there are three different methods used here. We found the 'bkde' estimate to be best.

# sigma1
d.sigma1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$sigma1)
sigma1.dat[i,12] <- d.sigma1$x[which.max(d.sigma1$y)]
dSJ.sigma1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$sigma1,bw = "sj")
sigma1.dat[i,13] <- dSJ.sigma1$x[which.max(dSJ.sigma1$y)]
dKS.sigma1 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$sigma1)
sigma1.dat[i,14] <- dKS.sigma1$x[which.max(dKS.sigma1$y)] 

# sigma2
d.sigma2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$sigma2)
sigma2.dat[i,12] <- d.sigma2$x[which.max(d.sigma2$y)]
dSJ.sigma2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$sigma2,bw = "sj")
sigma2.dat[i,13] <- dSJ.sigma2$x[which.max(dSJ.sigma2$y)]
dKS.sigma2 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$sigma2)
sigma2.dat[i,14] <- dKS.sigma2$x[which.max(dKS.sigma2$y)]

# a0
d.a0 <- density(out.lat.centerSIM$BUGSoutput$sims.list$a0)
a0.dat[i,12] <- d.a0$x[which.max(d.a0$y)]
dSJ.a0 <- density(out.lat.centerSIM$BUGSoutput$sims.list$a0,bw = "sj")
a0.dat[i,13] <- dSJ.a0$x[which.max(dSJ.a0$y)]
dKS.a0 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$a0)
a0.dat[i,14] <- dKS.a0$x[which.max(dKS.a0$y)] 

# a1
d.a1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$a1)
a1.dat[i,12] <- d.a1$x[which.max(d.a1$y)]
dSJ.a1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$a1,bw = "sj")
a1.dat[i,13] <- dSJ.a1$x[which.max(dSJ.a1$y)]
dKS.a1 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$a1)
a1.dat[i,14] <- dKS.a1$x[which.max(dKS.a1$y)] 

# a1adj
d.a1adj <- density(out.lat.centerSIM$BUGSoutput$sims.list$a1adj)
a1adj.dat[i,12] <- d.a1adj$x[which.max(d.a1adj$y)]
dSJ.a1adj <- density(out.lat.centerSIM$BUGSoutput$sims.list$a1adj,bw = "sj")
a1adj.dat[i,13] <- dSJ.a1adj$x[which.max(dSJ.a1adj$y)]
dKS.a1adj <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$a1adj)
a1adj.dat[i,14] <- dKS.a1adj$x[which.max(dKS.a1adj$y)] 

# b0
d.b0 <- density(out.lat.centerSIM$BUGSoutput$sims.list$b0)
b0.dat[i,12] <- d.b0$x[which.max(d.b0$y)]
dSJ.b0 <- density(out.lat.centerSIM$BUGSoutput$sims.list$b0,bw = "sj")
b0.dat[i,13] <- dSJ.b0$x[which.max(dSJ.b0$y)]
dKS.b0 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$b0)
b0.dat[i,14] <- dKS.b0$x[which.max(dKS.b0$y)] 

# b1
d.b1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$b1)
b1.dat[i,12] <- d.b1$x[which.max(d.b1$y)]
dSJ.b1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$b1,bw = "sj")
b1.dat[i,13] <- dSJ.b1$x[which.max(dSJ.b1$y)]
dKS.b1 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$b1)
b1.dat[i,14] <- dKS.b1$x[which.max(dKS.b1$y)] 

# b1adj
d.b1adj <- density(out.lat.centerSIM$BUGSoutput$sims.list$b1adj)
b1adj.dat[i,12] <- d.b1adj$x[which.max(d.b1adj$y)]
dSJ.b1adj <- density(out.lat.centerSIM$BUGSoutput$sims.list$b1adj,bw = "sj")
b1adj.dat[i,13] <- dSJ.b1adj$x[which.max(dSJ.b1adj$y)]
dKS.b1adj <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$b1adj)
b1adj.dat[i,14] <- dKS.b1adj$x[which.max(dKS.b1adj$y)] 

# b2
d.b2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$b2)
b2.dat[i,12] <- d.b2$x[which.max(d.b2$y)]
dSJ.b2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$b2,bw = "sj")
b2.dat[i,13] <- dSJ.b2$x[which.max(dSJ.b2$y)]
dKS.b2 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$b2)
b2.dat[i,14] <- dKS.b2$x[which.max(dKS.b2$y)] 

# alpha
d.alpha <- density(out.lat.centerSIM$BUGSoutput$sims.list$alpha)
alpha.dat[i,12] <- d.alpha$x[which.max(d.alpha$y)]
dSJ.alpha <- density(out.lat.centerSIM$BUGSoutput$sims.list$alpha,bw = "sj")
alpha.dat[i,13] <- dSJ.alpha$x[which.max(dSJ.alpha$y)]
dKS.alpha <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$alpha)
alpha.dat[i,14] <- dKS.alpha$x[which.max(dKS.alpha$y)] 

# beta
d.beta <- density(out.lat.centerSIM$BUGSoutput$sims.list$beta)
beta.dat[i,12] <- d.beta$x[which.max(d.beta$y)]
dSJ.beta <- density(out.lat.centerSIM$BUGSoutput$sims.list$beta,bw = "sj")
beta.dat[i,13] <- dSJ.beta$x[which.max(dSJ.beta$y)]
dKS.beta <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$beta)
beta.dat[i,14] <- dKS.beta$x[which.max(dKS.beta$y)] 

# logp1
d.logp1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$logp1)
logp1.dat[i,12] <- d.logp1$x[which.max(d.logp1$y)]
dSJ.logp1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$logp1,bw = "sj")
logp1.dat[i,13] <- dSJ.logp1$x[which.max(dSJ.logp1$y)]
dKS.logp1 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$logp1)
logp1.dat[i,14] <- dKS.logp1$x[which.max(dKS.logp1$y)] 

# logp2
d.logp2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$logp2)
logp2.dat[i,12] <- d.logp2$x[which.max(d.logp2$y)]
dSJ.logp2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$logp2,bw = "sj")
logp2.dat[i,13] <- dSJ.logp2$x[which.max(dSJ.logp2$y)]
dKS.logp2 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$logp2)
logp2.dat[i,14] <- dKS.logp2$x[which.max(dKS.logp2$y)] 

# p1
d.p1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$p1)
p1.dat[i,12] <- d.p1$x[which.max(d.p1$y)]
dSJ.p1 <- density(out.lat.centerSIM$BUGSoutput$sims.list$p1,bw = "sj")
p1.dat[i,13] <- dSJ.p1$x[which.max(dSJ.p1$y)]
dKS.p1 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$p1)
p1.dat[i,14] <- dKS.p1$x[which.max(dKS.p1$y)] 

# p2
d.p2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$p2)
p2.dat[i,12] <- d.p2$x[which.max(d.p2$y)]
dSJ.p2 <- density(out.lat.centerSIM$BUGSoutput$sims.list$p2,bw = "sj")
p2.dat[i,13] <- dSJ.p2$x[which.max(dSJ.p2$y)]
dKS.p2 <- bkde(out.lat.centerSIM$BUGSoutput$sims.list$p2)
p2.dat[i,14] <- dKS.p2$x[which.max(dKS.p2$y)] 




#################
  
  # Classify Discrete State: 0 if mean(State)<.5 ; 1 if mean(State) > .5
  Latent.State <- as.factor(as.numeric(out.lat.centerSIM$BUGSoutput$mean$S > .5))
  Latent.State <- factor(Latent.State, levels=levels(fit.cluster))
  
  #### Percent agreement of latent classification with k-means classification:
  Percent.Match.LatentKmeans[i] <- mean(fit.cluster==Latent.State)
  
  #### Percent agreement of latent classification with true basin:
  Percent.Match.LatentTrue[i] <- mean(LakeSim$TrueBasin1==Latent.State)
  
  #### Percent agreement of k=means classification with true basin:
  Percent.Match.KmeansTrue[i] <- mean(fit.cluster==LakeSim$TrueBasin1)
  
}




# Save output
write.csv(sigma1.dat,"sigma1dat.csv", row.names=FALSE)
write.csv(sigma2.dat,"sigma2dat.csv", row.names=FALSE)
write.csv(a0.dat,"a0dat.csv", row.names=FALSE)
write.csv(a1.dat,"a1dat.csv", row.names=FALSE)
write.csv(a1adj.dat,"a1adjdat.csv", row.names=FALSE)
write.csv(b0.dat,"b0dat.csv", row.names=FALSE)
write.csv(b1.dat,"b1dat.csv", row.names=FALSE)
write.csv(b1adj.dat,"b1adjdat.csv", row.names=FALSE)
write.csv(b2.dat,"b2dat.csv", row.names=FALSE)
write.csv(alpha.dat,"alphadat.csv", row.names=FALSE)
write.csv(beta.dat,"betadat.csv", row.names=FALSE)
write.csv(logp1.dat,"logp1dat.csv", row.names=FALSE)
write.csv(logp2.dat,"logp2dat.csv", row.names=FALSE)
write.csv(p1.dat,"p1dat.csv", row.names=FALSE)
write.csv(p2.dat,"p2dat.csv", row.names=FALSE)

write.csv(Percent.Match.LatentKmeans,"LatentKmeans.csv", row.names=FALSE)
write.csv(Percent.Match.LatentTrue,"LatentTrue.csv", row.names=FALSE)
write.csv(Percent.Match.KmeansTrue,"KmeansTrue.csv", row.names=FALSE)


write.csv(a0.kmeans,"a0Kmeans.csv", row.names=FALSE)
write.csv(a1.kmeans,"a1Kmeans.csv", row.names=FALSE)
write.csv(b0.kmeans,"b0Kmeans.csv", row.names=FALSE)
write.csv(b1.kmeans,"b1Kmeans.csv", row.names=FALSE)
write.csv(sigma1.kmeans,"sigma1Kmeans.csv", row.names=FALSE)
write.csv(sigma2.kmeans,"sigma2Kmeans.csv", row.names=FALSE)

write.csv(p1.kmeans,"p1Kmeans.csv", row.names=FALSE)
write.csv(p2.kmeans,"p2Kmeans.csv", row.names=FALSE)

write.csv(loo.dat,"loodat.csv", row.names=FALSE)

write.csv(a0LM.dat,"a0LMdat.csv", row.names=FALSE)
write.csv(b0LM.dat,"b0LMdat.csv", row.names=FALSE)
write.csv(sigma1LM.dat,"sigma1LMdat.csv", row.names=FALSE)




# This stuff saved for plotting last sim if I want to:
LakeSim$SAVerr <- SAV1.err
Latent.State <- as.factor(as.numeric(out.lat.centerSIM$BUGSoutput$mean$S > .5))
State.BLR <- as.numeric(Latent.State)-1
State.BLR[State.BLR==0] <- "Clear" ; State.BLR[State.BLR==1] <- "Turbid"
LakeSim$State.BLR <- State.BLR

Kmeans.state <- fit.dat$fit.cluster
State.KM <- as.numeric(fit.dat$fit.cluster)-1
State.KM[State.KM==0] <- "Clear"; State.KM[State.KM==1] <- "Turbid"
LakeSim$State.KM <- State.KM

LakeSim$log.TP <- lnP
LakeSim$log.Chla <- lnC
LakeSim$log.SAV <- lnSAV

write.csv(LakeSim,"LastSim.csv", row.names=FALSE)

setwd("..")
}