SI for submission

DEBKiss/DEBKiss_llik.R

@@ -0,0 +1,51 @@
+source("DEBKiss_model.R")
+
+DEBKiss_obs_model <- function(data, sim.data, samp){
+
+    w.obs <- simdat <- NA
+    test <- FALSE
+
+    ## choose at which times to evaluate the likelihood
+    w.obs <- which(sim.data[,"time"] %in% data$time)
+
+    ## The simulated data set at only the correct tims
+    simdat <- sim.data[w.obs,]
+
+    ## this is a small correction used to keep the arguments of the
+    ## log normal distribution from being exactly zero
+    ec <- 1e-6
+
+    ## The calculate of the cummulative number of eggs, calculated
+    ## from the total amount of energy in the reproductive buffer,
+    ## times the efficiency of the conversion of the reproductive
+    ## buffer to eggs, divided by the energy per egg
+    cumEggs <- (simdat[,"WR"]*as.numeric(samp[["yBA"]])/as.numeric(samp[["WB0"]])) + ec
+
+    ## this is code to produce a plot of the data + simulation for
+    ## testing purposes
+    if(test){
+        par(mfrow=c(1,2), bty="n")
+        plot(simdat[,"time"], simdat[,"Lw"], type="l", ylab="structural length",
+             xlab="time")
+        points(data$time, data$L, col=2)
+        plot(simdat[,"time"], cumEggs, type="l", ylab="total eggs", xlab="time")
+        points(data$time, data$Egg, col=3)
+    }
+
+    ## log likelihood for the length portion of the data 
+    llik.L <- sum(dlnorm(data$L, meanlog = log(simdat[,"Lw"] + ec),
+                         sdlog = samp[["sdlog.L"]], log = TRUE))
+
+    ## log likelihood for the cumulative eggs
+    llik.E <- sum(dlnorm(data$Egg+ec, meanlog = log(cumEggs),
+                         sdlog = samp[["sdlog.E"]], log = TRUE))
+
+    ## sum the egg and length components together
+    llik<- llik.L + llik.E
+
+    if(test) print(c(llik.L, llik.E, llik)) ## check the likelihood
+
+    return(llik)
+}
+
+

DEBKiss/DEBKiss_model.R

@@ -0,0 +1,58 @@
+## DEBKiss Process Model
+library(deSolve)
+
+DEBKiss1<-function(time,y,params){
+
+    with(as.list(c(y, params)),{
+
+        ec<-0.000001 ## value of WB below which hatching occurs, this
+                     ## is >0 for numerical reasons
+
+        ## physical length, Lw, is the input
+        L<-Lw*deltaM ## structural length        
+        L3<-L^3 ## structural volume         
+
+        ## calculating the length at puberty from the weight at
+        ## puberty
+        if(!is.null(Wp)) Lp<-(Wp/dV)^(1/3)
+
+        JA<-f*JAM*L^2 ## assimilation
+        JX<-JA/yAX ## feeding
+
+        ## first check if still egg stage and modify if needed
+        if(WB>ec){
+            JA <- JA*fb/f ## this sets f to fb, if needed
+            JX <- 0
+            dWB <- - JA
+        }else dWB<-0
+
+        JM<-exp(logJMv)*L3 ## volume maintance costs
+        JV<-yVA*(kappa*JA-JM) ## growth
+        JR<- (1-kappa)*JA ## repo buffer
+
+
+        ## Starvation conditions, after hatching
+        if(WB<=ec){
+            if(JA>JM){
+                if( JV<0 ){
+                    JV<-0
+                    JR<-JA-JM
+                }
+            }else{
+                JV<-(JA-JM)/yAV
+                JR<-0
+            }
+        }
+
+        if(L<Lp) JR<-0 ## check for being mature
+
+
+        dLw <- JV/(3*dV*L^2*deltaM) ## change in length
+        dWR <- JR ## change in reproductive buffer
+
+
+        list(c(dWB,dLw,dWR))
+    })
+
+}
+

DEBKiss/DEBKiss_snail.R

@@ -0,0 +1,207 @@
+### These provide the code from the package, the implemented model,
+### and likelihood, respectively
+library("deBInfer")
+source("DEBKiss_model.R")
+source("DEBKiss_llik.R")
+
+## First read in the simplified dataset
+dat<-read.csv("snail_1food_clean.csv", header=TRUE)
+
+
+## these are the initial parameters based on the fits from the
+## original paper. Note that because we're not beginning at egg laying
+## that we set the initial embryo buffer to 0
+
+y<-c(WB=0.000001, Lw=12.8, WR=0)
+parms<-c(deltaM=0.401, f=1, fb=0.8, JAM=0.11, yAX=0.8,
+         kappa=0.89, logJMv=log(0.008), yVA=0.8, yAV=0.8,
+         dV=0.1, Wp=70, yBA=0.95, WB0=0.15)
+
+
+## Plot the data together with the model output based on the
+## parameters from the paper
+times<-seq(0, 160, by=0.1)
+out <- ode(y, times, DEBKiss1, parms, method="rk4")
+
+par(mfrow=c(1,2),mai=c(.8,.8,.2,.1))# bty="n")
+plot(out[,1], out[,3], type="l", ylab="observed length", xlab="time (days)", lwd=2, col=2, lty=1, pch=23, ylim=c(12, 31), xlim=c(0,148))
+points(dat$time, dat$L, lwd=2)
+plot(out[,1], out[,4]*parms["yBA"]/parms["WB0"], type="l", ylab="total eggs", xlab="time (days)", ylim=c(0, 1100), lwd=2, col=4, lty=1, xlim=c(0,148))
+points(dat$time, dat$Egg, pch=25, lwd=2)
+
+
+##### Below is code necessary to perform inference for this model
+
+## setting up DEB parameters. We choose a subset to estimate based on
+## which parameters were estimated in the original paper. For this
+## example we do not estimate all of the parameters they did as to do
+## so would require fitting all 3 food levels simultaneous.
+
+## DEB parameters to estimate
+kappa <-debinfer_par(name = "kappa", var.type = "de", fixed = FALSE,
+                     value = parms["kappa"], prior = "unif",
+                     hypers = list(min=0, max=1),
+                     prop.var = 0.001, samp.type="rw")
+
+
+## because JMv is small and positive I sample it in log space
+logJMv <-debinfer_par(name = "logJMv", var.type = "de", fixed = FALSE,
+                   value=parms["logJMv"], prior = "norm",
+                   hypers = list(mean = 0, sd = 10 ),
+                   prop.var = 0.01, samp.type = "rw")
+
+
+## Observation parameters to estimate
+sdlog.L <-debinfer_par(name = "sdlog.L", var.type = "obs", fixed = FALSE,
+                       value = 1, prior = "lnorm",
+                       hypers = list(meanlog = 0, sdlog = 1 ),
+                       prop.var = c(4,5), samp.type = "rw-unif")
+
+sdlog.E <-debinfer_par(name = "sdlog.E", var.type = "obs", fixed = FALSE,
+                       value = 1, prior = "lnorm",
+                       hypers = list(meanlog = 0, sdlog = 1),
+                       prop.var = c(4,5), samp.type = "rw-unif")
+
+
+## These were estimated in the original paper, but we fix them here
+Wp <-debinfer_par(name = "Wp", var.type = "de", fixed = TRUE,
+                  value=parms["Wp"], prior = "lnorm",
+                  hypers = list(meanlog = 1, sdlog = 0.1 ),
+                  prop.var = c(4,5), samp.type = "rw-unif")
+
+JAM <-debinfer_par(name = "JAM", var.type = "de", fixed = TRUE,
+                   value=parms["JAM"], prior = "lnorm",
+                   hypers = list(meanlog = 0, sdlog = 0.1 ),
+                   prop.var = c(4,5), samp.type = "rw-unif")
+
+
+## These parameters were assumed and fixed in the original paper
+deltaM <-debinfer_par(name = "deltaM", var.type = "de", fixed = TRUE, value=parms["deltaM"])
+f <-debinfer_par(name = "f", var.type = "de", fixed = TRUE, value=parms["f"])
+fb <-debinfer_par(name = "fb", var.type = "de", fixed = TRUE, value=parms["fb"])
+yAX <-debinfer_par(name = "yAX", var.type = "de", fixed = TRUE, value=parms["yAX"])
+
+yVA <-debinfer_par(name = "yVA", var.type = "de", fixed = TRUE, value=parms["yVA"])
+yAV <-debinfer_par(name = "yAV", var.type = "de", fixed = TRUE, value=parms["yAV"])
+dV <-debinfer_par(name = "dV", var.type = "de", fixed = TRUE, value=parms["dV"])
+WB0 <- debinfer_par(name = "WB0", var.type = "obs", fixed = TRUE, value=parms["WB0"])
+
+## two values of yBA were considered in the paper, we focus on one
+yBA <-debinfer_par(name = "yBA", var.type = "obs", fixed = TRUE, value=parms["yBA"])
+
+
+
+## Initial conditions
+WB <- debinfer_par(name = "WB", var.type = "init", fixed = TRUE, value = y["WB"])
+L <- debinfer_par(name = "Lw", var.type = "init", fixed = TRUE, value = y["Lw"])
+WR <- debinfer_par(name = "WR", var.type = "init", fixed = TRUE, value = y["WR"])
+
+
+## Once all parameters have been individually specified, we combine
+## them using a setup function
+mcmc.pars<-setup_debinfer(kappa, sdlog.L, sdlog.E, deltaM, f, fb, JAM, yAX, logJMv,
+                          yVA, yAV, dV, Wp, yBA, WB0, WB, L, WR)
+
+## do inference with deBInfer
+## MCMC iterations
+iter <- 20000
+
+
+## inference call
+mcmc_samples <- de_mcmc(N = iter, data = dat, de.model = DEBKiss1,
+                        obs.model = DEBKiss_obs_model, all.params = mcmc.pars,
+                        Tmax = max(dat$time), data.times = dat$time, cnt = 500,
+                        plot = FALSE, verbose.mcmc = TRUE,
+                        solver = "ode", method="rk4")
+
+
+## plotting samples using built in functions, based on the CODA package in R
+plot(mcmc_samples, burnin=1000)
+
+pairs(mcmc_samples, burnin=1000)
+
+post_prior_densplot(mcmc_samples, burnin=1000)
+
+
+## creating posterior trajectories
+
+y<-c(WB=0.000001, Lw=12.8, WR=0)
+parms<-c(deltaM=0.401, f=1, fb=0.8, JAM=0.11, yAX=0.8,
+         kappa=0.89, logJMv=log(0.008), yVA=0.8, yAV=0.8,
+         dV=0.1, Wp=70, yBA=0.95, WB0=0.15)
+
+
+## Plot the data together with the model output based on the
+## parameters from the paper
+times<-seq(0, 160, by=0.1)
+
+thin<-seq(1001, 20000, by=10)
+samps.thin<-mcmc_samples$samples[thin,c(1,4)]
+
+ptemp<-parms
+## creating objects to hold the length and reproduction trajectories
+Lws<-WRs<-matrix(NA, nrow=length(times), ncol=length(thin))
+
+## this for loop will call the solver with the thinned parameter samples
+for(i in 1:length(thin)){
+    ## replace kappa and logJMv with the posterior samples
+    ptemp[6:7]<-samps.thin[i,]
+
+    ## solve the ODEs
+    out <- ode(y, times, DEBKiss1, ptemp, method="rk4")
+
+    ## save just the Length and reproduction
+    Lws[,i]<-out[,3]
+    WRs[,i]<-out[,4]
+}
+
+## calculates the posterior mean trajectories
+Lwmean<-apply(Lws, 1, mean)
+WRmean<-apply(WRs, 1, mean)
+
+## calculates the posterior credible intervals of the MEAN trajectory
+## (not the full prediction intervals that would include the
+## observation noise)
+Lwqs<-apply(Lws, 1, quantile, probs=c(0.025, 0.975))
+WRqs<-apply(WRs, 1, quantile, probs=c(0.025, 0.975))
+
+
+## plots the posterior means and CIs of the trajectories with the data
+par(mfrow=c(1,2),mai=c(.8,.8,.2,.1))# bty="n")
+plot(out[,1], Lwmean, type="l", ylab="observed length", xlab="time (days)", lwd=2, col=2, lty=1, pch=23, ylim=c(12, 31), xlim=c(0,148))
+lines(out[,1], Lwqs[1,], lty=3, lwd=2, col=2)
+lines(out[,1], Lwqs[2,], lty=3, lwd=2, col=2)
+points(dat$time, dat$L, lwd=2)
+
+plot(out[,1], WRmean*parms["yBA"]/parms["WB0"], type="l", ylab="total eggs", xlab="time (days)", ylim=c(0, 1100), lwd=2, col=4, lty=1, xlim=c(0,148))
+lines(out[,1], WRqs[1,]*parms["yBA"]/parms["WB0"], lty=3, lwd=2, col=4)
+lines(out[,1], WRqs[2,]*parms["yBA"]/parms["WB0"], lty=3, lwd=2, col=4)
+
+points(dat$time, dat$Egg, pch=25, lwd=2)
+
+
+
+## Extra plots
+source("plotting_extras.R")
+
+samps<-mcmc(mcmc_samples$samples[1001:20000,])
+ss<-summary(samps)
+
+par(mfrow = c(1,3))
+ps<-c("kappa", "logJMv", "sdlog.L", "sdlog.E")
+for (pp in ps[1:2]){
+    pretty_posterior_plot(samps, ss, ref.params = parms,
+                          param = pp, legend = FALSE)
+}
+
+plot.new()
+#for (pp in ps[3:4]){
+#    pretty_posterior_plot(samps, ss, ref.params = parms,
+#                          param = pp, legend = FALSE)
+#}
+
+legend("topleft", legend = c("true parameter value", "posterior mean value", "95% HPDI"), lty = c(2,1, NA), pch = c(NA,NA,15), col = c("black", "black", rethinking::col.alpha("black",0.15)))
+
+
+
+

DEBKiss/plotting_extras.R

@@ -0,0 +1,47 @@
+#HDPI figure code
+# rethinking is not yet on CRAN but can be installed using devtools:
+
+## install.packages(c("devtools","mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
+## library(devtools)
+## install_github("rmcelreath/rethinking")
+
+## The code for the figure itself is 
+
+library(coda)
+library(rethinking)
+
+pretty_posterior_plot <- function(samples, summary, ref.params, param, HDPI = 0.95, legend = TRUE){
+  rethinking::dens(unlist(samples[,param]), show.HPDI = 0.95, main = param)
+  abline(v = ref.params[param], lwd = 2, lty = 2)
+  abline(v = summary$statistics[param, "Mean"])
+  if(legend){
+      legend("topleft",
+             legend = c("true parameter value", "posterior mean value", "95% HPDI"),
+             lty = c(2,1, NA), pch = c(NA,NA,15),
+             col = c("black", "black", rethinking::col.alpha("black",0.15)))
+  }
+}
+
+### where samples is a coda object, i.e. the samples slot of a debinfer_result object (i.e. debinfer_result$samples)
+
+### summary is created using the summary method for the coda object (i.e. summary(debinfer_result$samples))
+
+### ref.params is a named vector with the true values
+
+### param is a string of the parameter to be plotted
+
+#pdf("figs/earthworm_fixedEG_freeobs_estimate_vs_true.pdf")
+#par(mfrow = c(3,3))
+#for (pp in names(freeparams(fixedEG_freeobs1))[1:8]){
+#    pretty_posterior_plot(fixedEG_freeobs_samples,
+#                          fixedEG_freeobs_summary,
+#                          ref.params = c(params_Lumter,
+#                              c(sdlog.EWw = 0.05, sdlog.R = 0.1)),
+#                          param = pp, legend = FALSE)
+#}
+#plot.new()
+#legend("topleft",
+#       legend = c("true parameter value", "posterior mean value", "95% HPDI"),
+#       lty = c(2,1, NA), pch = c(NA,NA,15),
+#       col = c("black", "black", rethinking::col.alpha("black",0.15)))
+#dev.off()

DEBKiss/snail_1food_clean.csv

@@ -0,0 +1,12 @@
+"time","L","Egg","food"
+0,12.92,0,100
+14,18.65,0,100
+28,21.55,0,100
+42,23.76,25.1,100
+56,25.67,142,100
+70,27.15,296,100
+84,27.95,455.3,100
+98,28.19,615.2,100
+112,28.27,766.3,100
+128,29.24,980,100
+140,29.53,1088.3,100

StandardDEB/data/Lumbricus_terrestris_mass_Butt1993.txt

@@ -0,0 +1,7 @@
+days_since_hatching mass_g se_up se_lw
+0.464 0.067 NA NA
+31.520 1.231 1.429 1.084
+60.722 3.429 3.803 3.056
+91.778 4.930 6.066 3.865
+121.444 6.094 7.409 4.790
+152.500 6.337 7.666 5.024