From 1d657c098a7135b1cabf3bb14a4d80a797abdeea Mon Sep 17 00:00:00 2001
From: Philipp Boersch-Supan <philipp.boerschsupan@bto.org>
Date: Mon, 25 Jun 2018 17:29:17 +0100
Subject: [PATCH] figure edits, added prior to hdpi plot function
---
DEBKiss/DEBKiss_snail.R | 16 ++++++++++++----
DEBKiss/plotting_extras.R | 15 +++++++++------
2 files changed, 21 insertions(+), 10 deletions(-)
diff --git a/DEBKiss/DEBKiss_snail.R b/DEBKiss/DEBKiss_snail.R
index 7702c75..16758be 100644
--- a/DEBKiss/DEBKiss_snail.R
+++ b/DEBKiss/DEBKiss_snail.R
@@ -161,9 +161,9 @@ parms<-c(deltaM=0.401, f=1, fb=0.8, JAM=0.11, yAX=0.8,
## Plot the data together with the model output based on the
## parameters from the paper
-times<-seq(0, 160, by=0.1)
+times<-seq(0, 160, by=0.5)
-samps.thin<-window(mcmc_samples$samples[,c('kappa','logJMv')], 1001, 20000, thin = 50)
+samps.thin<-window(mcmc_samples$samples[,c('kappa','logJMv')], 1001, 20000, thin = 10)
ptemp<-parms
## creating objects to hold the length and reproduction trajectories
@@ -171,6 +171,7 @@ Lws<-WRs<-matrix(NA, nrow=length(times), ncol=nrow(samps.thin))
## this for loop will call the solver with the thinned parameter samples
## NB: simulating the trajectories will take 3-5 minutes on a modern desktop computer
+Sys.time()
for(i in 1:nrow(samps.thin)){
## replace kappa and logJMv with the posterior samples
ptemp[6:7]<-samps.thin[i,]
@@ -182,6 +183,7 @@ for(i in 1:nrow(samps.thin)){
Lws[,i]<-out[,3]
WRs[,i]<-out[,4]
}
+Sys.time()
## calculates the posterior mean trajectories
@@ -197,12 +199,12 @@ 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))
+plot(out[,1], Lwmean, type="l", ylab="Observed shell length (mm)", 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))
+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)
@@ -211,6 +213,12 @@ points(dat$time, dat$Egg, pch=25, lwd=2)
## Extra plots
+## The plots showing shaded HDPIs use functionality from the rethinking package, which 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")
+
source("plotting_extras.R")
samps<-window(mcmc_samples$samples, 1001, 20000, thin = 10)
diff --git a/DEBKiss/plotting_extras.R b/DEBKiss/plotting_extras.R
index 8e0ee35..5089a07 100644
--- a/DEBKiss/plotting_extras.R
+++ b/DEBKiss/plotting_extras.R
@@ -10,16 +10,19 @@
library(coda)
library(rethinking)
-pretty_posterior_plot <- function(samples, summary, ref.params, param, HDPI = 0.95, legend = TRUE){
+pretty_posterior_plot <- function(samples, summary, ref.params, param, HDPI = 0.95, legend = TRUE, deBInfer_result = NULL){
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)))
+ if(!is.null(deBInfer_result)){
+ dprior <- paste("d", deBInfer_result$all.params[[param]]$prior,
+ sep = "")
+ plot.range <- seq(par("usr")[1], par("usr")[2], length.out = 100)
+ prior.dens <- do.call(dprior, c(list(x = plot.range),
+ deBInfer_result$all.params[[param]]$hypers))
+ lines(plot.range, prior.dens, col = 'red', lty = 3, lwd = 1.5)
}
+ 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)