From b789d086bf216bedc8f8d3730dda60ad64bc77db Mon Sep 17 00:00:00 2001 From: Gregory Britten Date: Tue, 9 Nov 2021 18:50:41 -0500 Subject: [PATCH] transferring to public --- R/events.R | 114 ++++++++++++++++++++++++ R/lab.r | 111 +++++++++++++++++++++++ R/plots.R | 216 +++++++++++++++++++++++++++++++++++++++++++++ README.md | 4 +- src/events_h2.stan | 79 +++++++++++++++++ src/lab_h2.stan | 84 ++++++++++++++++++ 6 files changed, 607 insertions(+), 1 deletion(-) create mode 100644 R/events.R create mode 100644 R/lab.r create mode 100644 R/plots.R create mode 100644 src/events_h2.stan create mode 100644 src/lab_h2.stan diff --git a/R/events.R b/R/events.R new file mode 100644 index 0000000..f78e6e9 --- /dev/null +++ b/R/events.R @@ -0,0 +1,114 @@ +rm(list=ls()) + +library(tidyverse) +library(rstan) +library(RColorBrewer) +options(mc.cores = parallel::detectCores()) +rstan_options(auto_write = TRUE) + +setwd('~/dropbox/working/taylor_host_parasite/github/') +################################################# +## DATA ######################################### +################################################# +d <- read_csv('../data/data.csv') +dd <- read.csv('../data/replicates.csv') %>% group_by(date) %>% summarise(sdI=sd(I)/sqrt(2),sdP=sd(P)/sqrt(2)) +d <- merge(d,dd,by='date') + +d$sdI[is.na(d$sdI)] <- median(d$sdI,na.rm=TRUE) +d$sdP[is.na(d$sdP)] <- median(d$sdP,na.rm=TRUE) +d$sdP[d$sdP==0] <- min(d$sdP[d$sdP!=0]) + +d <- d[order(d$days),] + +X1 <- d %>% filter(days < 120) +X2 <- d %>% filter(days > 120 & days < 170) +X3 <- d %>% filter(days > 170) + +XX <- list() +XX[[1]] <- X1 +XX[[2]] <- X2 +XX[[3]] <- X3 + +##################################################### +## STAN ############################################# +##################################################### +mod_events <- stan_model(file='src/events_h2.stan') + +dt <- 0.1 +ndays <- 25 +sig_perc <- 0.1 +rs <- c(0.21,0.24,0.35) +as <- c(4.40E-7, 5.58E-8, 1.34E-8) +hs <- c(2.96,2.79,2.46) +es <- c(600,550,40) +ms <- c(0.48,1.02,0.26) + +initf <- function(){list( + H0 = runif(1,0,2.0*12080.0), + I0 = fit_dat$I[1], + P0 = fit_dat$P[1], + r = mean(rs), + h1 = mean(hs), + h2 = 0.5, + a = mean(as), + m = mean(ms), + e = mean(es))} + +tmod <- seq(0,ndays,length.out=ndays/dt) + +POST_events <- list() + +for(i in 1:3){ + fit_dat <- XX[[i]] + fit_time <- fit_dat$days - fit_dat$days[1] + + sigma_I <- rep(mean(fit_dat$sdI),nrow(fit_dat)) + sigma_P <- rep(mean(fit_dat$sdP),nrow(fit_dat)) + + dat <- list(N=nrow(fit_dat), + I=fit_dat$I, + P=fit_dat$P, + dt=dt, + time_points=fit_time/dt + 1, + sigma_I=sigma_I, + sigma_P=sigma_P, + maxt = ndays/dt, + K = 1E7) + + mcmc <- sampling(mod_events,data=dat,init=initf) + + post <- extract(mcmc) + POST_events[[i]] <- post +} + +post <- POST_events[[1]] +flux_df_event1 <- data.frame(time_days=seq(0,ndays,length.out=dat$maxt), + growH=colMeans(post$growH*dt),growH_sd=apply(post$growH*dt,2,sd), + lossH=colMeans(post$lossH*dt),lossH_sd=apply(post$lossH*dt,2,sd), + hand =colMeans(post$hand*dt), hand_sd=apply(post$hand*dt, 2,sd), + growP=colMeans(post$growP*dt),growP_sd=apply(post$growP*dt,2,sd), + lossP=colMeans(post$lossP*dt),lossP_sd=apply(post$lossP*dt,2,sd)) +post <- POST_events[[2]] +flux_df_event2 <- data.frame(time_days=seq(0,ndays,length.out=dat$maxt), + growH=colMeans(post$growH*dt),growH_sd=apply(post$growH*dt,2,sd), + lossH=colMeans(post$lossH*dt),lossH_sd=apply(post$lossH*dt,2,sd), + hand =colMeans(post$hand*dt), hand_sd=apply(post$hand*dt, 2,sd), + growP=colMeans(post$growP*dt),growP_sd=apply(post$growP*dt,2,sd), + lossP=colMeans(post$lossP*dt),lossP_sd=apply(post$lossP*dt,2,sd)) +post <- POST_events[[3]] +flux_df_event3 <- data.frame(time_days=seq(0,ndays,length.out=dat$maxt), + growH=colMeans(post$growH*dt),growH_sd=apply(post$growH*dt,2,sd), + lossH=colMeans(post$lossH*dt),lossH_sd=apply(post$lossH*dt,2,sd), + hand =colMeans(post$hand*dt), hand_sd=apply(post$hand*dt, 2,sd), + growP=colMeans(post$growP*dt),growP_sd=apply(post$growP*dt,2,sd), + lossP=colMeans(post$lossP*dt),lossP_sd=apply(post$lossP*dt,2,sd)) + +write.csv(file='~/dropbox/working/taylor_host_parasite/results/flux_df_event1_h2.csv',flux_df_event1,row.names=FALSE) +write.csv(file='~/dropbox/working/taylor_host_parasite/results/flux_df_event2_h2.csv',flux_df_event2,row.names=FALSE) +write.csv(file='~/dropbox/working/taylor_host_parasite/results/flux_df_event3_h2.csv',flux_df_event3,row.names=FALSE) + +##--SAVE RESULTS--####################################### +save(file='../results/POST_events.rdata',POST_events) + + + diff --git a/R/lab.r b/R/lab.r new file mode 100644 index 0000000..92f16c0 --- /dev/null +++ b/R/lab.r @@ -0,0 +1,111 @@ +rm(list=ls()) + +library(tidyverse) +library(rstan) +library(RColorBrewer) +options(mc.cores = parallel::detectCores()) +rstan_options(auto_write = TRUE) + +setwd('~/dropbox/working/taylor_host_parasite/github/') +################################################# +## DATA ######################################### +################################################# +dh <- read.csv('../data/lab_heterocapsa.csv') +ds <- read.csv('../data/lab_scrippsiella.csv') +ddh <- dh %>% group_by(days) %>% summarise(Hmean=mean(H), Imean=mean(I), Pmean=mean(P), sdH=sd(H)/sqrt(3), sdI=sd(I)/sqrt(3),sdP=sd(P)/sqrt(3)) +dds <- ds %>% group_by(days) %>% summarise(Hmean=mean(H), Imean=mean(I), Pmean=mean(P), sdH=sd(H)/sqrt(3), sdI=sd(I)/sqrt(3),sdP=sd(P)/sqrt(3)) + +DAT <- list() +DAT[[1]] <- ddh +DAT[[2]] <- dds +##################################################### +## STAN ############################################# +##################################################### +mod <- stan_model(file='src/lab_h2.stan') + +dt <- 0.01 #divide by 100 to match time variable resolution in lab data +ndays <- 8 + +rs <- c(0.21,0.24,0.35) +as <- c(4.40E-7, 5.58E-8, 1.34E-8) +hs <- c(2.96,2.79,2.46) +es <- c(600,550,40) +ms <- c(0.48,1.02,0.26) + +initf <- function(){list( + H0 = fit_dat$Hmean[1], + I0 = fit_dat$Imean[1], + P0 = fit_dat$Pmean[1], + r = mean(rs), + h1 = mean(hs), + h2 = 0.5, + a = mean(as), + m = mean(ms), + e = mean(es))} + +tmod <- seq(0,ndays,length.out=ndays/dt) + +##--SET DATASET--########################### +POST_lab <- list() +for(i in 1:2){ + fit_dat <- DAT[[i]] #units of 1/100 day to match lab time resolution + fit_time <- (fit_dat$days - fit_dat$days[1])*100 + + sigma_H <- rep(mean(fit_dat$sdH),nrow(fit_dat)) + sigma_I <- rep(mean(fit_dat$sdI),nrow(fit_dat)) + sigma_P <- rep(mean(fit_dat$sdP),nrow(fit_dat)) + + dat <- list(N=nrow(fit_dat), + H=fit_dat$Hmean, + I=fit_dat$Imean, + P=fit_dat$Pmean, + dt=dt, + time_points=as.integer(fit_time + 1), + sigma_H=sigma_H, + sigma_I=sigma_I, + sigma_P=sigma_P, + maxt = ndays/dt, + K = 1E7) + + mcmc <- sampling(mod,data=dat,init=initf) + post <- extract(mcmc) + + POST_lab[[i]] <- post +} + +save(file='../results/POST_lab.rdata',POST_lab) + + +##--TABLES--############################################# +df1 <- data.frame(mean = c(mean(POST_lab[[1]]$r), mean(POST_lab[[1]]$a/(1+POST_lab[[1]]$h1*POST_lab[[1]]$a)), mean(POST_lab[[1]]$h2), mean(POST_lab[[1]]$m), mean(POST_lab[[1]]$e)), + sd = c(sd(POST_lab[[1]]$r), sd(POST_lab[[1]]$a/(1+POST_lab[[1]]$h1*POST_lab[[1]]$a)), sd(POST_lab[[1]]$h2), sd(POST_lab[[1]]$m), sd(POST_lab[[1]]$e))) +df1$cv <- df1$sd/df1$mean +df1 + +df2 <- data.frame(mean = c(mean(POST_lab[[2]]$r), mean(POST_lab[[2]]$a/(1+POST_lab[[2]]$h1*POST_lab[[2]]$a)), mean(POST_lab[[2]]$h2), mean(POST_lab[[2]]$m), mean(POST_lab[[2]]$e)), + sd = c(sd(POST_lab[[2]]$r), sd(POST_lab[[2]]$a/(1+POST_lab[[2]]$h1*POST_lab[[2]]$a)), sd(POST_lab[[2]]$h2), sd(POST_lab[[2]]$m), sd(POST_lab[[2]]$e))) +df2$cv <- df2$sd/df2$mean +df2 + +##--FLUXES--############################################# +#post <- POST[[1]] +post <- POST_lab[[1]] +flux_df_heterocapsa <- data.frame(time_days=seq(0,ndays,length.out=dat$maxt), + growH=colMeans(post$growH*dt),growH_sd=apply(post$growH*dt,2,sd), + lossH=colMeans(post$lossH*dt),lossH_sd=apply(post$lossH*dt,2,sd), + hand =colMeans(post$hand*dt), hand_sd=apply(post$hand*dt,2,sd), + growP=colMeans(post$growP*dt),growP_sd=apply(post$growP*dt,2,sd), + lossP=colMeans(post$lossP*dt),lossP_sd=apply(post$lossP*dt,2,sd)) + +#post <- POST[[2]] +post <- POST_lab[[2]] +flux_df_scrippsiella <- data.frame(time_days=seq(0,ndays,length.out=dat$maxt), + growH=colMeans(post$growH*dt),growH_sd=apply(post$growH*dt,2,sd), + lossH=colMeans(post$lossH*dt),lossH_sd=apply(post$lossH*dt,2,sd), + hand =colMeans(post$hand*dt), hand_sd=apply(post$hand*dt,2,sd), + growP=colMeans(post$growP*dt),growP_sd=apply(post$growP*dt,2,sd), + lossP=colMeans(post$lossP*dt),lossP_sd=apply(post$lossP*dt,2,sd)) + +write.csv(file='~/dropbox/working/taylor_host_parasite/results/flux_df_heterocapsa_h2.csv',flux_df_heterocapsa,row.names=FALSE) +write.csv(file='~/dropbox/working/taylor_host_parasite/results/flux_df_scrippsiella_h2.csv',flux_df_scrippsiella,row.names=FALSE) + diff --git a/R/plots.R b/R/plots.R new file mode 100644 index 0000000..be75547 --- /dev/null +++ b/R/plots.R @@ -0,0 +1,216 @@ +rm(list=ls()) +library(tidyverse) +library(RColorBrewer) + +############################################################################### +##--FIELD EVENTS--############################################################# +############################################################################### +##--LOAD DATA--############################## +setwd('~/dropbox/working/taylor_host_parasite/github/') +d <- read_csv('../data/data.csv') +dd <- read.csv('../data/replicates.csv') %>% group_by(date) %>% summarise(sdI=sd(I)/sqrt(2),sdP=sd(P)/sqrt(2)) +d <- merge(d,dd,by='date') + +##--ESTIMATE OF OBSERVATION ERROR--#################### +d$sdI[is.na(d$sdI)] <- median(d$sdI,na.rm=TRUE) +d$sdP[is.na(d$sdP)] <- median(d$sdP,na.rm=TRUE) +d$sdP[d$sdP==0] <- min(d$sdP[d$sdP!=0]) + +d <- d[order(d$days),] + +##--SPLIT UP INTO EVENTS--################ +XX <- list() +XX[[1]] <- d %>% filter(days < 120) +XX[[2]] <- d %>% filter(days > 120 & days < 170) +XX[[3]] <- d %>% filter(days > 170) + +##--TIME-STEPPING PARAMETERS FOR MODEL--######################## +dt <- 0.1 +ndays <- 25 +tmod <- seq(0,ndays,length.out=ndays/dt) + +##--LOAD EVENTS FITTING RESULTS FOR PLOTTING--################### +load('../results/POST_events.rdata') + +############################################################################### +##--LAB EVENTS--############################################################### +############################################################################### +##--LOAD LAB DATA--######## +setwd('~/dropbox/working/taylor_host_parasite/github/') +dh <- read.csv('../data/lab_heterocapsa.csv') +ds <- read.csv('../data/lab_scrippsiella.csv') + +DAT <- list() +DAT[[1]] <- dh %>% group_by(days) %>% summarise(Hmean=mean(H), Imean=mean(I), Pmean=mean(P), sdH=sd(H)/sqrt(3), sdI=sd(I)/sqrt(3),sdP=sd(P)/sqrt(3)) +DAT[[2]] <- ds %>% group_by(days) %>% summarise(Hmean=mean(H), Imean=mean(I), Pmean=mean(P), sdH=sd(H)/sqrt(3), sdI=sd(I)/sqrt(3),sdP=sd(P)/sqrt(3)) + +dt <- 0.01 #divide by 100 to match time variable resolution in lab data +ndays <- 8 +tmod_lab <- seq(0,ndays,length.out=ndays/dt) + +##--LOAD FITTING RESULTS FOR LAB EVENTS--########################## +load('../results/POST_lab.rdata') + +################################################################################# +##--PLOTS--###################################################################### +################################################################################# + +##--FIELD EVENT TIME SERIES--######################## +pdf('../plots/time_series_fits_events.pdf',height=6,width=8.5) +par(mfcol=c(3,3),mar=c(2,2,2,2),oma=c(2,2,2,2)) +for(i in 1:3){ + post <- POST_events[[i]] + X <- XX[[i]] + + Hmod <- colMeans(post$x[,,1]) + Hmodsd <- apply(post$x[,,1],2,sd) + Imod <- colMeans(post$x[,,2]) + Imodsd <- apply(post$x[,,2],2,sd) + Pmod <- colMeans(post$x[,,3]) + Pmodsd <- apply(post$x[,,3],2,sd) + + plot(tmod,Hmod,type='l',lty=1,xlab='',ylab='',ylim=c(0,max(Hmod + 2*Hmodsd)),lwd=2) + lines(tmod,Hmod - 2*Hmodsd,type='l',lwd=1,lty=2) + lines(tmod,Hmod + 2*Hmodsd,type='l',lwd=1,lty=2) + mtext(paste('Event #',i,sep=''),adj=0) + if(i==1) mtext(side=2,'Host Abundance',line=2.5) + plot(tmod,Imod,type='l',lty=1,xlab='',ylab='',ylim=c(0,max(Imod + 2*Imodsd, X$I)),lwd=2) + lines(tmod,Imod - 2*Imodsd,type='l',lwd=1,lty=2) + lines(tmod,Imod + 2*Imodsd,type='l',lwd=1,lty=2) + points(X$days-X$days[1], X$I,pch=8,col='red') + if(i==1) mtext(side=2,'Infected Host Abundance',line=2.5) + plot(tmod,Pmod,type='l',lty=1,xlab='',ylab='',ylim=c(0,max(Pmod + 2*Pmodsd, X$P)),lwd=2) + lines(tmod,Pmod - 2*Pmodsd,type='l',lwd=1,lty=2) + lines(tmod,Pmod + 2*Pmodsd,type='l',lwd=1,lty=2) + points(X$days-X$days[1], X$P,pch=8,col='red') + if(i==1) mtext(side=2,'Spore Abundance',line=2.5) +} +mtext(side=1,outer=TRUE,'Days',line=0.5) +dev.off() + +##--LAB EVENT TIME SERIES--############################ +pdf('../plots/time_series_fits_lab.pdf',height=6,width=8.5*(2/3)) +par(mfcol=c(3,2),mar=c(2,2,2,2),oma=c(2,2,2,2)) +for(i in 1:2){ + post <- POST_lab[[i]] + dat <- DAT[[i]] + + Hmod <- colMeans(post$x[,,1]) + Hmodsd <- apply(post$x[,,1],2,sd) + Imod <- colMeans(post$x[,,2]) + Imodsd <- apply(post$x[,,2],2,sd) + Pmod <- colMeans(post$x[,,3]) + Pmodsd <- apply(post$x[,,3],2,sd) + + plot(tmod_lab,Hmod,type='l',lty=1,xlab='',ylab='',ylim=c(0,max(Hmod + 2*Hmodsd, dat$Hmean)),lwd=2) + lines(tmod_lab,Hmod - 2*Hmodsd,type='l',lwd=1,lty=2) + lines(tmod_lab,Hmod + 2*Hmodsd,type='l',lwd=1,lty=2) + if(i==1) mtext(side=2,'Host Abundance',line=2.5,cex=0.8) + points(dat$days-dat$days[1], dat$Hmean,pch=8,col='red') + if(i==1) mtext(expression(italic('heterocapsa')),adj=0) + if(i==2) mtext(expression(italic('scrippsiella')),adj=0) + + plot(tmod_lab,Imod,type='l',lty=1,xlab='',ylab='',ylim=c(0,max(Imod + 2*Imodsd, dat$Imean)),lwd=2) + lines(tmod_lab,Imod - 2*Imodsd,type='l',lwd=1,lty=2) + lines(tmod_lab,Imod + 2*Imodsd,type='l',lwd=1,lty=2) + points(dat$days-dat$days[1], dat$Imean,pch=8,col='red') + if(i==1) mtext(side=2,'Infected Host Abundance',line=2.5,cex=0.8) + + plot(tmod_lab,Pmod,type='l',lty=1,xlab='',ylab='',ylim=c(0,max(Pmod + 2*Pmodsd, dat$Pmean)),lwd=2) + lines(tmod_lab,Pmod - 2*Pmodsd,type='l',lwd=1,lty=2) + lines(tmod_lab,Pmod + 2*Pmodsd,type='l',lwd=1,lty=2) + points(dat$days-dat$days[1], dat$Pmean,pch=8,col='red') + if(i==1) mtext(side=2,'Spore Abundance',line=2.5,cex=0.8) +} +mtext(side=1,outer=TRUE,'Days',line=0.5) +dev.off() + + +########################################################################################### +## POSTERIOR PARAMETER DISTRIBUTIONS ACROSS EVENTS ######################################## +########################################################################################### + +post1_events <- POST_events[[1]] +post2_events <- POST_events[[2]] +post3_events <- POST_events[[3]] +post_scrip <- POST_lab[[2]] + +cols <- brewer.pal(n = 8, name = "Dark2") +nbrks <- 60 + +hist2 <- function(x,xlim,col,breaks,ylim,add){ + hist(x,xlim=xlim,col=col,breaks=breaks,main='',ylim=ylim,add=add,border=TRUE) + abline(v=quantile(x,probs=c(0.025,0.5,0.975)),lty=c(2,1,2),col=col,lwd=0.8) +} + +##--PARAMETER DISTRIBUTIONS--############################### +pdf('../plots/parameter_distributions_11_05_2021.pdf',height=5,width=8) +par(mfrow=c(2,3),mar=c(2,2,2,2),oma=c(2,2,2,2)) +uplim <- max(c(post1_events$r,post3_events$r,post_scrip$r)) + +hist2(x=post1_events$r, xlim=c(0,uplim), col=cols[1], breaks=seq(-0.1,uplim,length.out=nbrks),ylim=c(0,1500),add=FALSE) +hist2(x=post3_events$r, xlim=c(0,1), col=cols[3], breaks=seq(-0.1,uplim,length.out=nbrks),ylim=c(0,1500),add=TRUE) + mtext(expression(italic('r')),side=1,line=2.5) +hist2(x=post_scrip$r, xlim=c(0,uplim), col=cols[2], breaks=seq(-0.1,uplim,length.out=nbrks),ylim=c(0,1500),add=TRUE) + legend('topright',col=cols[c(1,3,2)],legend=c('Event #1','Event #3',expression(italic('Scrippsiella'))),bty='n',pch=15) + +uplim <- log10(max(c(post1_events$a/(1+post1_events$a*post1_events$h1),post3_events$a/(1+post3_events$a*post3_events$h1),post_scrip$a/(1+post_scrip$a*post_scrip$h1)))) +hist2(x=log10(post1_events$a/(1+post1_events$a*post1_events$h1)), xlim=c(-8.3,uplim),col=cols[1], breaks=seq(-8.3,uplim,length.out=nbrks),ylim=c(0,1500),add=FALSE) +hist2(x=log10(post3_events$a/(1+post3_events$a*post3_events$h1)), xlim=c(-8.3,uplim),col=cols[3], breaks=seq(-8.3,uplim,length.out=nbrks),ylim=c(0,1500),add=TRUE) +hist2(x=log10(post_scrip$a/(1+post_scrip$a*post_scrip$h1)), xlim=c(-8.3,uplim),col=cols[2], breaks=seq(-8.3,uplim,length.out=nbrks),ylim=c(0,1500),add=TRUE) +mtext(expression('log10{'~italic('a/(1+ah'[1]*')')~'}'),side=1,line=2.5) + +uplim <- max(c(post1_events$h2,post3_events$h2,post_scrip$h2)) +hist2(x=post1_events$h2, xlim=c(0,uplim), col=cols[1], breaks=seq(-0.1,uplim,length.out=nbrks), ylim=c(0,1500), add=FALSE) +hist2(x=post3_events$h2, xlim=c(0,uplim), col=cols[3], breaks=seq(-0.1,uplim,length.out=nbrks), ylim=c(0,1500), add=TRUE) +hist2(x=post_scrip$h2, xlim=c(0,uplim), col=cols[2], breaks=seq(-0.1,uplim,length.out=nbrks), ylim=c(0,1500), add=TRUE) + mtext(expression(italic('h'[2])),side=1,line=2.5) + +uplim <- max(c(post1_events$m,post3_events$m,post_scrip$m)) +hist2(x=post1_events$m, xlim=c(0,uplim),col=cols[1], breaks=seq(0,uplim,length.out=nbrks),ylim=c(0,1500),add=FALSE) +hist2(x=post3_events$m, xlim=c(0,uplim),col=cols[3], breaks=seq(0,uplim,length.out=nbrks),ylim=c(0,1500),add=TRUE) +hist2(x=post_scrip$m, xlim=c(0,uplim),col=cols[2], breaks=seq(0,uplim,length.out=nbrks),ylim=c(0,1500),add=TRUE) + mtext(expression(italic('m')),side=1,line=2.5) + +uplim <- max(c(post1_events$e,post3_events$e,post_scrip$e)) +hist2(x=post1_events$e,xlim=c(0,uplim),col=cols[1], breaks=seq(0,uplim,length.out=nbrks*2),ylim=c(0,1500),add=FALSE) +hist2(x=post3_events$e,xlim=c(0,uplim),col=cols[3], breaks=seq(0,uplim,length.out=nbrks*2),ylim=c(0,1500),add=TRUE) +hist2(x=post_scrip$e,xlim=c(0,uplim),col=cols[2], breaks=seq(0,uplim,length.out=nbrks*2),ylim=c(0,1500),add=TRUE) + mtext(expression(italic('e')),side=1,line=2.5) +dev.off() + + +##--PARAMETER DISTRIBUTIONS FOR H1 and H2--########################################## +pdf('../plots/parameter_distributions.pdf',height=4,width=8) +par(mfrow=c(1,2),mar=c(2,2,2,2),oma=c(2,2,2,2),cex.axis=0.8) + uplim <- max(c(post1_events$h1,post3_events$h1,post_scrip$h1)) + hist2(x=post1_events$h1, xlim=c(0,uplim),col=cols[1], breaks=seq(0,uplim,length.out=nbrks),ylim=c(0,600),add=FALSE) + hist2(x=post3_events$h1, xlim=c(0,uplim),col=cols[3], breaks=seq(0,uplim,length.out=nbrks),ylim=c(0,600),add=TRUE) + hist2(x=post_scrip$h1, xlim=c(0,uplim),col=cols[2], breaks=seq(0,uplim,length.out=nbrks),ylim=c(0,600),add=TRUE) + mtext(expression(italic('h'[1])),side=1,line=2.5) + legend('topleft',col=cols[c(1,3,2)],legend=c('Event #1','Event #3',expression(italic('Scrippsiella'))),bty='n',pch=15) + + uplim <- max(c(post1_events$h2,post3_events$h2,post_scrip$h2)) + hist2(x=post1_events$h2, xlim=c(0,uplim), col=cols[1], breaks=seq(-0.1,uplim,length.out=nbrks), ylim=c(0,600), add=FALSE) + hist2(x=post3_events$h2, xlim=c(0,uplim), col=cols[3], breaks=seq(-0.1,uplim,length.out=nbrks), ylim=c(0,600), add=TRUE) + hist2(x=post_scrip$h2, xlim=c(0,uplim), col=cols[2], breaks=seq(-0.1,uplim,length.out=nbrks), ylim=c(0,600), add=TRUE) + mtext(expression(italic('h'[2])),side=1,line=2.5) +dev.off() + +############################################################################# +##--TABLE OF RATES--######################################################### +#############################################################################S +df1 <- data.frame(mean = c(mean(post1_events$r), mean(post1_events$a/(1+post1_events$h1*post1_events$a)), mean(post1_events$h2), mean(post1_events$m), mean(post1_events$e)), + sd = c(sd(post1_events$r), sd(post1_events$a/(1+post1_events$h1*post1_events$a)), sd(post1_events$h2), sd(post1_events$m), sd(post1_events$e))) +df1$cv <- df1$sd/df1$mean +df1 + +df2 <- data.frame(mean = c(mean(post2_events$r), mean(post2_events$a/(1+post2_events$h1*post2_events$a)), mean(post2_events$h2), mean(post2_events$m), mean(post2_events$e)), + sd = c(sd(post2_events$r), sd(post2_events$a/(1+post2_events$h1*post2_events$a)), sd(post2_events$h2), sd(post2_events$m), sd(post2_events$e))) +df2$cv <- df2$sd/df2$mean +df2 + +df3 <- data.frame(mean = c(mean(post3_events$r), mean(post3_events$a/(1+post3_events$h1*post3_events$a)), mean(post3_events$h2), mean(post3_events$m), mean(post3_events$e)), + sd = c(sd(post3_events$r), sd(post3_events$a/(1+post3_events$h1*post3_events$a)), sd(post3_events$h2), sd(post3_events$m), sd(post3_events$e))) +df3$cv <- df3$sd/df3$mean +df3 diff --git a/README.md b/README.md index a2a905c..1bbe616 100644 --- a/README.md +++ b/README.md @@ -1 +1,3 @@ -# host_parasite_public \ No newline at end of file +# Code for submission of Sehein et al. +`/src` constains .stan files +`/R` contains r scripts for fitting models and plotting results diff --git a/src/events_h2.stan b/src/events_h2.stan new file mode 100644 index 0000000..f021535 --- /dev/null +++ b/src/events_h2.stan @@ -0,0 +1,79 @@ +data { + int N; //number of data points + vector[N] I; //read new data + vector[N] P; + real dt; + int time_points[N]; + real sigma_I[N]; + real sigma_P[N]; + int maxt; + real K; +} +parameters { + real H0; //initial conditions are parameters + real I0; + real P0; + real r; //biological parameters + real h1; + real h2; + real a; + real m; + real e; +} +transformed parameters{ + real x[maxt,3]; //matrix to store variables while looping + + real growH[maxt]; //declare intermediaries + real lossH[maxt]; + real dHdt; + real hand[maxt]; + real dIdt; + real growP[maxt]; + real lossP[maxt]; + real dPdt; + + x[1,1] = H0; //specify initial conditions in x + x[1,2] = I0; + x[1,3] = P0; + + growH[1] = r * x[1,1] * (1-x[1,1]/K); + lossH[1] = (a * x[1,1]/(1+(a * h1 * x[1,1]))) * x[1,3]; + hand[1] = x[1,2]/h2; + growP[1] = e * hand[1]; + lossP[1] = m * x[1,3]; + + for(j in 2:maxt){ //time step the model + growH[j] = r * x[j-1,1] * (1-x[j-1,1]/K); + lossH[j] = (a * x[j-1,1]/(1+(a * h1 * x[j-1,1]))) * x[j-1,3]; + dHdt = growH[j] - lossH[j]; + + hand[j] = x[j-1,2]/h2; + dIdt = lossH[j] - hand[j]; + + growP[j] = e * hand[j]; + lossP[j] = m * x[j-1,3]; + dPdt = growP[j] - lossP[j]; + + x[j,1] = fmax(x[j-1,1] + dHdt*dt, 1E-15); //store variable in x; fmax disallows negative values + x[j,2] = fmax(x[j-1,2] + dIdt*dt, 1E-15); + x[j,3] = fmax(x[j-1,3] + dPdt*dt, 1E-15); + } +} +model { + + //--PRIORS-------------------------------------// + r ~ normal(0.2666667, 0.07371115); + a ~ normal(1.697333e-07, 2.350159e-07); + h1 ~ normal(2.736667, 0.2542309); + h2 ~ normal(0.5, 0.5); + m ~ normal(0.5866667, 0.3910669); + e ~ normal(396.6667, 0.5*309.8925); + + I0 ~ uniform(0, 10.0*I[1]); + P0 ~ uniform(0, 10.0*P[1]); + + //--LIKELIHOOD-----------------------------------// + I ~ normal(x[time_points,2], sigma_I); + P ~ normal(x[time_points,3], sigma_P); +} + diff --git a/src/lab_h2.stan b/src/lab_h2.stan new file mode 100644 index 0000000..d3bd557 --- /dev/null +++ b/src/lab_h2.stan @@ -0,0 +1,84 @@ +data { + int N; //number of data points + vector[N] H; //read new data + vector[N] I; //read new data + vector[N] P; + real dt; + int time_points[N]; + real sigma_H[N]; + real sigma_I[N]; + real sigma_P[N]; + int maxt; + real K; +} +parameters { + real H0; //initial conditions are parameters + real I0; + real P0; + real r; //biological parameters + real h1; + real h2; + + real a; + real m; + real e; +} +transformed parameters{ + real x[maxt,3]; //matrix to store variables while looping + + real growH[maxt]; //declare intermediaries + real lossH[maxt]; + real dHdt; + real hand[maxt]; + real dIdt; + real growP[maxt]; + real lossP[maxt]; + real dPdt; + + x[1,1] = H0; //specify initial conditions in x + x[1,2] = I0; + x[1,3] = P0; + + growH[1] = r * x[1,1] * (1-x[1,1]/K); + lossH[1] = (a * x[1,1]/(1+(a * h1 * x[1,1]))) * x[1,3]; + hand[1] = x[1,2]/h2; + growP[1] = e * hand[1]; + lossP[1] = m * x[1,3]; + + + for(j in 2:maxt){ //time step the model + growH[j] = r * x[j-1,1] * (1-x[j-1,1]/K); + lossH[j] = (a * x[j-1,1]/(1+(a * h1 * x[j-1,1]))) * x[j-1,3]; + dHdt = growH[j] - lossH[j]; + + hand[j] = x[j-1,2]/h2; + dIdt = lossH[j] - hand[j]; + + growP[j] = e * hand[j]; + lossP[j] = m * x[j-1,3]; + dPdt = growP[j] - lossP[j]; + + x[j,1] = fmax(x[j-1,1] + dHdt*dt, 1E-15); //store variable in x; fmax disallows negative values + x[j,2] = fmax(x[j-1,2] + dIdt*dt, 1E-15); + x[j,3] = fmax(x[j-1,3] + dPdt*dt, 1E-15); + } +} +model { + //--PRIORS-------------------------------------// + r ~ normal(0.2666667, 0.07371115); + a ~ normal(1.697333e-07, 2.350159e-07); + h1 ~ normal(2.736667, 0.2542309); + h2 ~ normal(0.5, 0.5); + m ~ normal(0.5866667, 0.3910669); + e ~ normal(396.6667, 0.5*309.8925); + + H0 ~ normal(H[1], 0.5*H[1]); + I0 ~ normal(I[1], 0.5*I[1]); + P0 ~ normal(P[1], 0.5*P[1]); + + //--LIKELIHOOD-----------------------------------// + H ~ normal(x[time_points,1], sigma_H); + I ~ normal(x[time_points,2], sigma_I); + P ~ normal(x[time_points,3], sigma_P); +} +