Rarely do we ever work with one model. Instead, you will typically want to evaluate how the posterior of a parameter changes with changes to data or model. It would be nice to scale up the functionality of ‘postpack’ to allow simple or complex post-processing tasks from posterior samples from similar models using the same consistent workflow. This vignette illustrates how this can be acheived using ‘postpack’.
I have only recently began working with the ‘postpack’ functions in this regard efficiently, so the ideas are not fully fleshed out yet. In future versions of ‘postpack’, more functionality may be added in this regard.
The most significant trick to use is to store multiple
mcmc.list objects as elements of a larger list
object. Suppose you have two mcmc.list objects from two
highly similar models, named cjs and
cjs_no_rho (see vignette("example-mcmclists")
or ?cjs for more details).
Load these objects into your session:
library(postpack)
data(cjs)
data(cjs_no_rho)And create a list object with them, where each element
is an mcmc.list object:
post_list = list(cjs, cjs_no_rho)Be sure to assign element names, which allows tracking which output is from which model later on:
names(post_list) = c("est_rho", "no_rho")From here, the world is wide open to you. Anything you would do with
one mcmc.list object, you can do with two (or any number)
through the use of the base R apply() family of functions.
For example, extract the dimensions of the saved output for each
model:
sapply(post_list, post_dim)## est_rho no_rho
## burn 11000 11000
## post_burn 50000 50000
## thin 200 200
## chains 2 2
## saved 500 500
## params 21 21Notice the two have identical dimensions. You can see that each model has the same parameters saved:
sapply(post_list, get_params, type = "base_index")## est_rho no_rho
## [1,] "B0" "B0"
## [2,] "sig_B0" "sig_B0"
## [3,] "B1" "B1"
## [4,] "sig_B1" "sig_B1"
## [5,] "b0[1]" "b0[1]"
## [6,] "b0[2]" "b0[2]"
## [7,] "b0[3]" "b0[3]"
## [8,] "b0[4]" "b0[4]"
## [9,] "b0[5]" "b0[5]"
## [10,] "b1[1]" "b1[1]"
## [11,] "b1[2]" "b1[2]"
## [12,] "b1[3]" "b1[3]"
## [13,] "b1[4]" "b1[4]"
## [14,] "b1[5]" "b1[5]"
## [15,] "SIG[1,1]" "SIG[1,1]"
## [16,] "SIG[2,1]" "SIG[2,1]"
## [17,] "SIG[1,2]" "SIG[1,2]"
## [18,] "SIG[2,2]" "SIG[2,2]"
## [19,] "p[2]" "p[2]"
## [20,] "p[3]" "p[3]"
## [21,] "p[4]" "p[4]"You could verify that all parameters in both models converged well
according to the Rhat statistic:
sapply(post_list, function(model) post_summ(model, ".", Rhat = TRUE)["Rhat",])## est_rho no_rho
## B0 1.011 1.020
## sig_B0 1.022 1.030
## B1 1.033 1.020
## sig_B1 1.043 1.010
## b0[1] 1.000 1.004
## b0[2] 0.998 0.998
## b0[3] 1.004 1.008
## b0[4] 0.999 1.002
## b0[5] 0.999 1.008
## b1[1] 1.005 0.999
## b1[2] 1.006 1.011
## b1[3] 1.007 1.009
## b1[4] 1.000 1.013
## b1[5] 1.017 0.999
## SIG[1,1] 1.182 1.141
## SIG[2,1] 1.009 NaN
## SIG[1,2] 1.009 NaN
## SIG[2,2] 1.199 1.074
## p[2] 1.000 1.000
## p[3] 0.999 1.004
## p[4] 1.005 1.003(The NaN values for "SIG[1,2]" and
"SIG[2,1]" in the no_rho model are because
those had the same value each MCMC iteration).
You could extract the summaries of the global survival coefficients from each model to see that the qualitative inference does not depend on the feature that differs between these two models:
lapply(post_list, function(model) post_summ(model, "^B", digits = 2))## $est_rho
## B0 B1
## mean 1.59 0.42
## sd 0.49 0.21
## 50% 1.57 0.40
## 2.5% 0.55 0.06
## 97.5% 2.59 0.86
##
## $no_rho
## B0 B1
## mean 1.62 0.38
## sd 0.47 0.19
## 50% 1.60 0.39
## 2.5% 0.81 -0.02
## 97.5% 2.48 0.69Or verify that the detection probabilities are not affected either:
lapply(post_list, function(model) post_summ(model, "p", digits = 2))## $est_rho
## p[2] p[3] p[4]
## mean 0.74 0.70 0.63
## sd 0.02 0.03 0.04
## 50% 0.74 0.70 0.63
## 2.5% 0.69 0.64 0.55
## 97.5% 0.79 0.76 0.72
##
## $no_rho
## p[2] p[3] p[4]
## mean 0.74 0.70 0.64
## sd 0.03 0.03 0.05
## 50% 0.75 0.70 0.64
## 2.5% 0.69 0.64 0.56
## 97.5% 0.79 0.76 0.73You can embed more steps in the function that is applied to each
mcmc.list object. For example, the code below obtains the
posterior mean correlation matrix for each model:
lapply(post_list, function(model) {
SIG_decomp = vcov_decomp(model, "SIG")
rho_mean = post_summ(SIG_decomp, "rho")["mean",]
array_format(rho_mean)
})## $est_rho
## [,1] [,2]
## [1,] 1.0000000 0.1847515
## [2,] 0.1847515 1.0000000
##
## $no_rho
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1Or more complex still, to predict the survival probability between two consecutive detection arrays for fish of different sizes in each year and model:
lapply(post_list, function(model) {
# 2SDs below average length, average length, and 2SDs above average length
# model was fitted with length scaled and centered
pred_length = c(-2,0,2)
# extract posterior mean of random coefficients
b0_mean = post_summ(model, "b0")["mean",]
b1_mean = post_summ(model, "b1")["mean",]
# predict survival each year from coefficients at each length
pred_phi = sapply(1:5, function(y) {
logit_phi = b0_mean[y] + b1_mean[y] * pred_length
phi = exp(logit_phi)/(1 + exp(logit_phi))
round(phi, 2)
})
# give dimension names
dimnames(pred_phi) = list(c("small", "average", "large"), paste0("y", 1:5))
# return the predicted survival
return(pred_phi)
})## $est_rho
## y1 y2 y3 y4 y5
## small 0.66 0.81 0.67 0.70 0.60
## average 0.80 0.91 0.85 0.82 0.75
## large 0.89 0.96 0.94 0.90 0.86
##
## $no_rho
## y1 y2 y3 y4 y5
## small 0.66 0.83 0.66 0.70 0.59
## average 0.80 0.91 0.85 0.82 0.75
## large 0.89 0.96 0.95 0.89 0.86