##Here we estimate degradation rate of pathogens from our singly infected populations library(R2jags) library(mcmcplots) library(reshape2) library(ggplot2) library(deSolve) library(bbmle)= library(MASS) library(Rmisc) require(here) load("./code_output/jagsout_maturationrate.Rdata") load("./code_output/jagsout_deathrate.Rdata") load("./code_output/jagsout_birthrate.Rdata") load("./code_output/mass_pastinf.Rdata") load("./code_output/mass_metschinf.Rdata") load("./code_output/jagsout_metschspores.Rdata") load("./code_output/jagsout_pastspores.Rdata") load("./code_output/uninf_K_est.Rdata") load("./code_output/jagsout_metschdeath.Rdata") #################################################################################3 deathrate.posterior<-drestimate$BUGSoutput$sims.list$trait.mean birthrate.posterior<-brestimate$BUGSoutput$sims.list$trait.mean maturationrate.posterior<-juestimate$BUGSoutput$sims.list$trait.mean pastinf.posterior<-pastmassprobestimate$BUGSoutput$sims.list$infectivity metschinf.posterior<-metschmassprobestimate$BUGSoutput$sims.list$infectivity pastsingleposterior<-pastestimate$BUGSoutput$sims.list$intercept metschsingleposterior<-metschestimate$BUGSoutput$sims.list$intercept K.posterior<-single_K_estimate$BUGSoutput$sims.list$K ##Ok, now want to fit model to actual population data. ##Do all the transformations from data visualization, except don't multiply things by ten , and leave each replicate in 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*1 summary$juveniles<-summary_max$juveniles*1 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*1 summary$pastonly<-summary$pastonly*1 summary$coinf<-summary$coinf*1 summary$uninf<-summary$uninf*1 ##############################################################3 #Remove coinfected treatments epi_single <- summary[ which(summary$treatment!='1'), ] epi_single <- epi_single[ which(epi_single$treatment!='5'), ] epi_single <- epi_single[ which(epi_single$treatment!='6'), ] epi_single <- epi_single[ which(epi_single$treatment!='7'), ] #################################################################### ######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 feedtau <- 1 / (feedvar * feedvar) F ~ dnorm(feednorm,feedtau)T(0.1,) MD ~ dunif(0,1) PD ~ dunif(0,1) #randomized sampling of posterior LTrand ~ dunif(0,length.lifetable) IFPrand ~ dunif(0,length.infPfunction) IFMrand ~ dunif(0,length.infMfunction) SRrand ~ dunif(0,length.spores) Krand ~ dunif(0,length.K) LTrand.round <- round(LTrand) IFPrand.round <- round(IFPrand) IFMrand.round <- round(IFPrand) SRrand.round <- round(SRrand) Krand.round <- round(Krand) B <- birthrate.posterior[LTrand.round] DR <- deathrate.posterior[LTrand.round] MI <- metschinf.posterior[IFMrand.round] PI <- pastinf.posterior[IFPrand.round] Bpp <- pastsingleposterior[SRrand.round] * 1000 Bmm <- metschsingleposterior[SRrand.round] * 1000 K <- K.posterior[Krand.round] ##placeholders MM <- 6 MJ <- 7 PM <- 11 PJ <- 14 SR <- 0.98 V <- 1000 J <- 6 #Pre time points for time delay for (z in 1:40){ for (y in 1:3){ S[z,y] <- 0 Jco[z,y] <- 0 IM[z,y] <- 0 IMco[z,y] <- 0 IP[z,y] <- 0 IPco[z,y] <- 0 EM[z,y] <- 0 EMco[z,y] <- 0 EP[z,y] <- 0 EPco[z,y] <- 0 EnvM[z,y] <- 0 EnvP[z,y] <- 0 N[z,y] <- 0 Mprev[z,y] <- 0 Prev[z,y] <- 0 } } S[41,1] <- 100 S[41,2] <- 100 S[41,3] <- 100 Jco[41,1] <- 0 Jco[41,2] <- 0 Jco[41,3] <- 0 IM[41,1] <- 0 IM[41,2] <- 0 IM[41,3] <- 0 IMco[41,1] <- 0 IMco[41,2] <- 0 IMco[41,3] <- 0 IP[41,1] <- 0 IP[41,2] <- 0 IP[41,3] <- 0 IPco[41,1] <- 0 IPco[41,2] <- 0 IPco[41,3] <- 0 EM[41,1] <- 0 EM[41,2] <- 0 EM[41,3] <- 0 EMco[41,1] <- 0 EMco[41,2] <- 0 EMco[41,3] <- 0 EP[41,1] <- 0 EP[41,2] <- 0 EP[41,3] <- 0 EPco[41,1] <- 0 EPco[41,2] <- 0 EPco[41,3] <- 0 EnvM[41,1] <- 0 EnvM[41,2] <- 190000 EnvM[41,3] <- 0 EnvP[41,1] <- 0 EnvP[41,2] <- 0 EnvP[41,3] <- 1000000 N[41,1] <- 100 N[41,2] <- 100 N[41,3] <- 100 Mprev[41,1] <- 0 Mprev[41,2] <- 0 Mprev[41,3] <- 0 Pprev[41,1] <- 0 Pprev[41,2] <- 0 Pprev[41,3] <- 0 for (j in 41:106){ for (m in 1:3){ EMco[j+1,m] <- ifelse( (S[j,m] * (((1 - exp( -(F * N[j,m]) / V)) / (N[j,m]+0.000001)) * EnvM[j,m] * MI)) > S[j,m] * SR * (1 - DR), S[j,m] * SR * (1 - DR), (S[j,m] * (((1 - exp( -(F * N[j,m]) / V)) / (N[j,m]+0.000001)) * EnvM[j,m] * MI))) EM[j+1,m] <- ifelse((EM[j,m] * SR * (1 - DR) + EMco[j+1,m] - EMco[j-MJ,m] * (pow((SR * (1 - DR)), MJ))) > 0, EM[j,m] * SR * (1 - DR) + EMco[j+1,m] - EMco[j-MJ,m] * (pow((SR * (1 - DR)), MJ)), 0) IM[j+1,m] <- ifelse(( IM[j,m] * SR * (1 - 1/MM) + (EMco[j-MJ,m] * (pow((SR * (1 - DR)), MJ)))) > 0, (IM[j,m] * SR * (1 - 1/MM) + (EMco[j-MJ,m] * (pow((SR * (1 - DR)), MJ)))), 0) EnvM[j+1,m] <-ifelse((EnvM[j,m] * MD * exp(-(N[j,m] * F) / V) + IM[j,m] * (1/MM) * (Bmm)) > 0, EnvM[j,m] * MD * exp(-(N[j,m] * F) / V) + IM[j,m] * (1/MM) * (Bmm), 0) EPco[j+1,m] <- ifelse((S[j,m] * (((1 - exp( -(F * N[j,m]) / V)) / (N[j,m]+0.000001)) * EnvP[j,m]) * PI) > S[j,m] * SR * (1 - DR), S[j,m] * SR * (1 - DR), (S[j,m] * (((1 - exp( -(F * N[j,m]) / V)) / (N[j,m]+0.000001)) * EnvP[j,m]) * PI) ) EP[j+1,m] <- ifelse((EP[j,m] * SR * (1 - DR) + EPco[j+1,m] - EPco[j-PJ,m] * (pow((SR * (1 - DR)), PJ))) > 0, EP[j,m] * SR * (1 - DR) + EPco[j+1,m] - EPco[j-PJ,m] * (pow((SR * (1 - DR)), PJ)), 0) IP[j+1,m] <- ifelse(( IP[j,m] * SR * (1 - 1/PM) + (EPco[j-PJ,m] * (pow((SR * (1 - DR)), PJ)))) > 0, (IP[j,m] * SR * (1 - 1/PM) + (EPco[j-PJ,m] * (pow((SR * (1 - DR)), PJ)))), 0) EnvP[j+1,m] <- ifelse((EnvP[j,m] * PD - EnvP[j,m] * (1-exp(-(N[j,m] * F) / V)) * PI + IP[j,m] * (1/PM) * (Bpp)) > 0, EnvP[j,m] * PD - EnvP[j,m] * (1-exp(-(N[j,m] * F) / V)) * PI + IP[j,m] * (1/PM) * (Bpp), 0) S[j+1,m] <- ifelse(ifelse((S[j,m] * SR * (1 - DR) + Jco[j - J,m] * (pow(SR,J)) - EMco[j+1,m] - EPco[j+1,m]) > 0, S[j,m] * SR * (1 - DR) + Jco[j - J,m] * (pow(SR, J)) - EMco[j+1,m] - EPco[j+1,m], 0) < 1, 0, ifelse((S[j,m] * SR * (1 - DR) + Jco[j - J,m] * (pow(SR,J)) - EMco[j+1,m] - EPco[j+1,m]) > 0, S[j,m] * SR * (1 - DR) + Jco[j - J,m] * (pow(SR, J)) - EMco[j+1,m] - EPco[j+1,m], 0)) Jco[j+1,m] <- ifelse( (S[j,m] * B * (1 - N[j,m] / K)) > 0, (S[j,m] * B * (1 - N[j,m]/K)), 0) N[j+1,m] <- S[j+1,m] + EM[j+1,m] + IM[j+1,m] + EP[j+1,m] + IP[j+1,m] Mprev[j+1,m] <- (IM[j+1,m] / (N[j+1,m] + 0.0000001)) + 0.000000001 Pprev[j+1,m] <- (IP[j+1,m] / (N[j+1,m] + 0.0000001)) + 0.000000001 #EnvP[j+1,m] <- ifelse((EnvP[j,m] * PD - EnvP[j,m] * (1-exp(-(N[j,m] * F) / V)) * PI + IP[j,m] * DR * (Bpp)) > 0, EnvP[j,m] * PD - EnvP[j,m] * (1-exp(-(N[j,m] * F) / V)) * PI + IP[j,m] * DR * (Bpp), 0) } } for ( i in 1:X ){ #so dbinom would be fed in the number of trials, and the probability of success, to generate possible number of successes. # So the thing that would match up with is the number of infecteds in the raw data. Is it possible to put the raw data for n into the distribution estimating # raw prevalence? I think so. metschonlydata[i] ~ dbin(Mprev[(date[i]+40),treatment[i]],pop_size[i]) pastonlydata[i] ~ dbin(Pprev[(date[i]+40),treatment[i]],pop_size[i]) } } ", file = model.loc) jags.data = list( feednorm=3.8, feedvar=0.3, pastonlydata = epi_single$pastonly*10, metschonlydata = epi_single$metschonly*10, pop_size = epi_single$pop_size*10, X = length(epi_single$uninf), treatment = as.numeric(epi_single$treatment), date = epi_single$date, deathrate.posterior = as.vector(deathrate.posterior), birthrate.posterior = as.vector(birthrate.posterior), K.posterior = as.vector(K.posterior), length.K = length(K.posterior), pastinf.posterior = as.vector(pastinf.posterior), metschinf.posterior = as.vector(metschinf.posterior), pastsingleposterior = as.vector(pastsingleposterior), metschsingleposterior = as.vector(metschsingleposterior), length.lifetable = length(deathrate.posterior), length.infPfunction = length(pastinf.posterior), length.infMfunction = length(metschinf.posterior),length.spores = length(pastsingleposterior)) jags.params = c("MD","PD") inits<-function(){list( MD = 0.95, PD = 0.95)} single_deg_estimate = jags(jags.data, parameters.to.save = jags.params, model.file = model.loc, n.chains = 3, n.burnin = 10, n.thin = 15, n.iter = 10000, DIC = TRUE, inits = inits) mcmcplot(single_deg_estimate) # check diagnostics - change name as needed single_deg_estimate$BUGSoutput$summary save(single_deg_estimate, file = "./code_output/single_deginf_estimate.Rdata") ##Ok, found mistake. Currently, initial EnvM values set at per ml, but additions are to the total spores counted. Ok, adjust starting number, ##And asjust number added. #spores counted*10,000= spores/ml #*.1 ml = *1000