##This code estimates carrying capacity based on uninfected population data ##Load necessary packages library(R2jags) library(mcmcplots) library(reshape2) library(ggplot2) library(deSolve) library(bbmle) library(MASS) library(rmisc) #Delete previous files require(here) ##Load posterior estimates of healthy individuals lifetables data load("./code_output/jagsout_birthrate.Rdata") load("./code_output/jagsout_deathrate.Rdata") load("./code_output/jagsout_maturationrate.Rdata") deathrate.posterior<-drestimate$BUGSoutput$sims.list$trait.mean birthrate.posterior<-brestimate$BUGSoutput$sims.list$trait.mean maturationrate.posterior<-juestimate$BUGSoutput$sims.list$trait.mean ##Ok, now want to fit model to actual population data. ##format spreadsheets epidemic = read.csv("./raw_data/experiment_raw_data_zeros.csv") epidemic$treatment<-as.factor(epidemic$treatment) epidemic_sum=subset(epidemic, select=c(date,treatment,replicate,metschonly,pastonly,coinf,uninf)) epidemic_max=subset(epidemic, select=c(date,treatment,replicate,animal,juveniles)) epidemic_mean=subset(epidemic, select=c(date,treatment,replicate,eggs,metschonly,pastonly,coinf,uninf)) summary=aggregate(. ~ date+replicate+treatment, data = epidemic_sum, FUN = sum ) summary_max=aggregate(. ~ date+replicate+treatment, data = epidemic_max, FUN = max ) summary_mean=aggregate(. ~ date+replicate+treatment, data = epidemic_mean, FUN = mean ) summary$pop_size<-summary_max$animal summary$juveniles<-summary_max$juveniles summary$fertility<-summary_mean$eggs summary$metschprev<-summary_mean$metschonly summary$pastprev<-summary_mean$pastonly summary$coinfprev<-summary_mean$coinf summary$uninfprev<-summary_mean$uninf summary$metschonly<-summary$metschonly summary$pastonly<-summary$pastonly summary$coinf<-summary$coinf summary$uninf<-summary$uninf #Now aggregate along date and treatment for mean and 95% confidence intervals epi_single <- summary[ which(summary$treatment=='2'), ] #################################################################### ######Model for only using modes of distribution #################################################################### model.loc = ("uninf-singlesfit.txt") jagsscript = cat(" model { ## Priors, trying to find carrying capacity #feeding rate distribution, currently drawn from distribution in shocket et al. 2018 K ~ dunif(1,200) #create number for randomized sampling of posterior LTrand ~ dunif(0,length.lifetable) LTrand.round <- round(LTrand) BR <- birthrate.posterior[LTrand.round] DR <- deathrate.posterior[LTrand.round] SR <- 0.98 ##death via sampling V <- 1000 ##Volume J <- 6 ##Vast majority of juvenile posterior falls between 5.5 nd 6.5, and cannot have stochastic indexing, so set this as a steady # #Pre time points so that time delay indexing does not go negative for (z in 1:40){ S[z] <- 0 } S[41] <- 100 ##Starting population size for (j in 41:110){ S[j+1] <- ifelse((S[j - J] * BR * (1 - S[j - J] / K)) > 0, ifelse((S[j] * SR * (1 - DR) + (S[j - J] * BR * (1 - S[j - J] / K)) * (SR ^ J)) < 1, 0, (S[j] * SR * (1 - DR) + (S[j - J] * BR * (1 - S[j - J] / K)) * (SR ^ J))), ifelse((S[j] * SR * (1 - DR)) < 1, 0, (S[j] * SR * (1 - DR)))) } for ( i in 1:X ){ uninfdata[i] ~ dpois(lambdas[i]) log(lambdas[i]) <- S[(date[i]+40)] } } ", file = model.loc) jags.data = list( uninfdata = epi_single$uninf*10, X = length(epi_single$uninf), date = epi_single$date, deathrate.posterior = as.vector(deathrate.posterior), birthrate.posterior = as.vector(birthrate.posterior), length.lifetable = length(deathrate.posterior)) jags.params = c("K") inits<-function(){list(K = 10)} single_K_estimate = jags(jags.data, parameters.to.save = jags.params, model.file = model.loc, n.chains = 3, n.burnin = 1, n.thin = 15, n.iter = 1000000, DIC = TRUE, inits = inits) mcmcplot(single_K_estimate) # check diagnostics - change name as needed single_K_estimate$BUGSoutput$summary save(single_K_estimate, file = "./code_output/uninf_K_est.Rdata")