knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.width = 8, fig.height = 6)
rm(list = ls()) library(eliciteg)
Here we show how to sample from the Mixture Dirichlet posterior. This is done to sample from the multiple experts' prior beliefs about movement of right whales in and around the mid-Atlantic region. The protocol for sampling is as follows:
In the model from Schick et al. (2013), monthly movement probabilities are sampled from a Dirichlet posterior. Informed priors come from one expert—Philip Hamilton, New England Aquarium. The sightings data come from the North Atlantic Right Whale Consortium; missing data are imputed with multiple imputation. Because we are using mulitple imputation, we need to sample to approximate the posterior.
In this vignette, we will work through a toy example using simulated data first, and then proceed to use real data. Here is a toy dataset, representing movements of whales in 5 regions.
maleExT <- matrix(c(1, 0, 2, 1, 1, 6, 750, 0, 1, 3, 3, 0, 0, 3, 1, 9, 250, 13, 5, 2, 1, 300, 0, 0, 3), nrow = 5, ncol = 5, byrow = TRUE) colnames(maleExT) <- c("BOF", "GOM", "GSC", "JL", "MIDA")
For this toy example, we'll assume two experts. Here are their priors:
priorT_1 <- matrix(c(7.502, 0.39, 0.39, 0.39, 3.952, 0.144, 0.047, 0.047, 0.047, 0.003, 10.144, 0.048, 0.041, 0.047, 0.003, 0.144, 0.049, 0.047, 0.047, 0.003, 5, 0.003, 0.003, 0.003, 2.52), nrow = 5, ncol = 5, byrow = TRUE) priorT_2 <- matrix(c(5, 2.394, 1.394, 1.3948, 3.952, 5, 1.047, 1.0479, 1.0479, 0.003, 5, 0.047, 1.0479, 2.0479, 0.003, 5, 0.047, 1.0479, 1.0479, 0.003, 5, 1.003, 1.003, 1.003, 2.052), nrow = 5, ncol = 5, byrow = TRUE) priorT <- array(1, dim = c(5, 5, 2)) priorT[, , 1] <- priorT_1 priorT[, , 2] <- priorT_2 nreg <- ncol(priorT_1)
With those in place, now we proceed to calculate the $C$ normalizing constants. First, we set up a placeholder matrix:
cmatOut <- matrix(NA, nrow = dim(priorT)[3], ncol = nreg, dimnames = list(c('expert1', 'expert2'), colnames(maleExT)))
Note that the rows in cmatOut
correspond to the number of experts, while the columns correspond to the number of regions. With that set up, we can populate it with a nested set of loops:
for(i in seq_along(1:dim(maleExT)[2])){ for(j in seq_along(1:dim(priorT)[3])){ prior <- priorT[, , j] cmatOut[j, i] <- calcC(maleExT[, i], prior[, i]) } } cmatOut
This looks fine, but you'll note quite large (negative) values for GOM. In situtations where the values in cmatOut
get very big and negative, taking the exponent yields 0 values, which we can't have in this algorithm.
The solution Michail came up with is to add some numerical padding to all experts' values in the instances where one expert has a very negative value. This preserves the relationships among experts, while keeping things tractable numerically.
Here's the solution in the form of a new function:
transformC
The processing logic in the function is to:
sweep
to add those column-wise to the matrixIn use, the function produces this output:
cmatOutexp <- transformC(cmatOut, tval = -600) cmatOutexp
With the $C$ matrix assembled, we can now assemble the $K$ weights.
Mathematically, the overall posterior is represented as:
$$ f^{(1)} (\pi) = \sum^J_{j=1} k_j^{(1)} f^{(1)}_j (\pi),$$
with
$$k_j^{(1)} = \frac{k_j^{(0)} C_j}{\sum_{i=1}^J k_j^{(0)} C_j} .$$
To sample, you choose a component distribution $j$with probability $k_j^{(0)}$. We then sample from $f_j^{(1)} (\pi) = Dir \left( d^{(0)}{j1} + n_1, \ldots, d^{(0)}{jM} + n_M\right )$, assuming our data are $(n_1, \ldots, n_M)$.
Our function to assemble this matrix is calck
. A priori the prior weights for $K$ are equal to 1 / (# of experts)
, i.e. we give them equal weight. calcK
returns the posterior weights for the experts.
Note that above we have the prior represented as $k_j^{(0)}$, which in our algorithm will always be equal to 1 / (# of experts)
. That is, we'll be giving each expert equal weight. In the calcK
function $k_j^{(0)}$ is represented as kprior
. What we return from this function is $k_j^{(1)}$, i.e. the posterior weights.
calcK
In use:
kmat <- calcK(cmatOutexp) kmat colSums(kmat)
# this block is to save these objects for later comparison to the outputs from Michail's code cmatOutexpSim <- cmatOutexp cmatSim <- cmatOut kmatSim <- kmat priorTsim <- priorT
That output, then, is a matrix where each column is the probability vector (that sums to 1), which we will use to select the expert within the Gibbs loop. For this transition that we're exploring here, Expert 2 is much more likely to be sampled during each step, i.e. their weight is much higher. You can also see in the from MIDA to all regions transition, where the priors are almost identical, we'll be sampling from these two experts almost equally.
With the toy example, we'll sample directly from a Dirichlet making use of the rdirichlet()
function in the gtools
library. First we can view the prior for each expert:
Let's see how our experts stack up for the prior for the first region: Bay of Fundy (BOF). First we sample them to create the mixture prior
i <- 1 # the region where we assume a whale is present library(gtools) prdirsample <- matrix(NA, nrow = 1000, ncol = nreg) for (ig in 1:1000) { numexp <- nrow(kmat) k <- rep(1 / numexp, numexp) idx <- which(rmultinom(1, 1, k) == 1) prior <- priorT[, i, idx] prdirsample[ig, ] <- rdirichlet(1, prior) } colnames(prdirsample) <- c(colnames(kmat))
And then we can plot those to see the mixture prior distribution:
library(ggplot2) library(reshape2) xlp <- melt(as.data.frame(prdirsample)) # draw the histograms ggplot(data = xlp, aes(x = value))+ geom_density()+ facet_wrap(~variable, nrow = 1)+ xlim(0, 1)+ ylim(0, 10)
And then we can move on to the posterior
postdirsample <- rdirichlet(1000, rep(1, nreg)) # placeholder for sampled values for (ig in 1:1000) { idx <- which(rmultinom(1, 1, kmat[, i]) == 1) prior <- priorT[, i, idx] postdirsample[ig, ] <- rdirichlet(1, maleExT[, i] + prior) }
And then we can plot the posterior (black line) and include the prior (blue lines):
colnames(postdirsample) <- colnames(kmat) xl <- melt(postdirsample) colnames(xlp) <- c('Var2', 'value') # draw the histograms ggplot(data = xl, aes(x = value, group = Var2))+ geom_density(lwd = 1)+ geom_density(data = xlp, aes(x = value), colour = 'blue')+ facet_wrap(~Var2, nrow = 1)+ xlim(0, 1)+ ylim(0, 10)
Recall that the data for this region was: r maleExT[, 1]
. So we're essentially just getting the data back since we're adding the same prior for each region. We say essentially, because while the means would be the same the variance for $D(1, 1, 1, 1, 1)$ is different than $D(5, 5, 5, 5, 5)$.
Since that seems to work with a toy data example for one region, let's try it with a full 9 * 9 data matrix. Again, we won't use multiple imputation here, though we will in the paper; This is just to see it in action. Note that where above we sampled directly from the Dirichlet, here we'll use a slightly different algorithm (though n.b. if you look at the internals of rdirichlet
, you'll see it's the same sampling logic).
We'll remake the two supporting matrices (actually I'll print 3 out: 1) cmat
with logc values; 2) cmat
after transforming it; and 3) kmat
). I won't show the code, as it's identical to above (save for the dimensions).
data(moveData) cmatOut <- matrix(NA, nrow = length(priorList$males), ncol = dim(maleEx)[2], dimnames = list(c('exp1', 'exp2', 'exp3', 'exp4', 'exp5', 'exp6', 'exp7', 'exp8'), colnames(maleEx))) for(i in seq_along(1:dim(maleEx)[2])){ for(j in seq_along(1:length(priorList$males))){ priorT <- priorList$males[[j]] cmatOut[j, i] <- calcC(maleEx[, i], priorT[, i]) } } round(cmatOut, 2) cmatOutexp <- transformC(cmatOut) cmatOutexp kmat <- calcK(cmatOutexp) round(kmat, 2)
Here I'll simulate a Gibbs loop for males.
ng <- 100 nregions <- ncol(kmat) post <- matrix(0, nregions, nregions) summh <- post for(i in seq_along(1:ng)){ for(j in seq_along(1:nregions)){ idx <- which(rmultinom(1, 1, kmat[, j]) == 1) # chooses the expert prior <- priorList$males[[idx]][, j] # gets prior for that expert di <- matrix(rgamma(nregions * 1, shape = maleEx[, j] + prior, scale = 1), nregions, 1) post[, j] <- di / matrix(colSums(di), nregions, 1, byrow = T) } summh <- summh + post } malemovePostprob <- summh / ng colnames(malemovePostprob) <- row.names(malemovePostprob) <- colnames(maleEx) round(malemovePostprob, 3)
As with males, we'll remake the two (err...3) supporting matrices as above.
data(moveData) cmatOut <- matrix(NA, nrow = length(priorList$females), ncol = dim(femEx)[2], dimnames = list(c('exp1', 'exp2', 'exp3', 'exp4', 'exp5', 'exp6', 'exp7', 'exp8'), colnames(femEx))) for(i in seq_along(1:dim(femEx)[2])){ for(j in seq_along(1:length(priorList$females))){ priorT <- priorList$females[[j]] cmatOut[j, i] <- calcC(femEx[, i], priorT[, i]) } } round(cmatOut, 2) cmatOutexp <- transformC(cmatOut) cmatOutexp kmat <- calcK(cmatOutexp) round(kmat, 2)
Here I'll simulate a Gibbs loop for females.
ng <- 100 nregions <- ncol(kmat) post <- matrix(0, nregions, nregions) sumfh <- post for(i in seq_along(1:ng)){ for(j in seq_along(1:nregions)){ idx <- which(rmultinom(1, 1, kmat[, j]) == 1) # chooses the expert prior <- priorList$females[[idx]][, j] # gets prior for that expert di <- matrix(rgamma(nregions * 1, shape = femEx[, j] + prior, scale = 1), nregions, 1) post[, j] <- di / matrix(colSums(di), nregions, 1, byrow = T) } sumfh <- sumfh + post } femalemovePostprob <- sumfh / ng colnames(femalemovePostprob) <- row.names(femalemovePostprob) <- colnames(femEx) round(femalemovePostprob, 3)
The code above is mostly written by me (Rob), and this section will compare it against code written by Michail. If kmat
agree using the same input data, then we should be confident about the methods. I saved the initial data, priors, and matrices to objects with new names for comparison:
cmatOutexpSim
cmatSim
kmatSim
priorTsim
Now we'll use Michail's code to re-create the matrices.
k <- matrix(NA, nrow = nreg, ncol = 2) # nreg for the number of regions where whales are assumed present in turn # 2 for the 2 experts matnames <- list(c(colnames(kmatSim)), c(row.names(kmatSim))) logc_pr <- matrix(NA, nrow = nreg, ncol = 2, dimnames = matnames) logc_pos <- matrix(NA, nrow = nreg, ncol = 2, dimnames = matnames) logc <- matrix(NA, nrow = nreg, ncol = 2, dimnames = matnames) cmat <- matrix(NA, nrow = nreg, ncol = 2, dimnames = matnames) kpost <- matrix(NA, nrow = nreg, ncol = 2, dimnames = matnames) for (i in 1:nreg){# nreg for the # of regions where we will assume a whale is currently present for (ex in 1:2){ # as many experts as you have k[i, ex] <- 0.5 logc_pr[i, ex] <- lgamma(sum(priorTsim[, i, ex])) - sum(lgamma(priorTsim[, i, ex])) logc_pos[i, ex] <- lgamma(sum(priorTsim[, i, ex] + maleExT[, i])) - sum(lgamma(priorTsim[, i, ex] + maleExT[, i])) logc[i, ex] <- logc_pr[i, ex] - logc_pos[i, ex] cmat[i, ex] <- exp(logc[i, ex]) } if (logc[i, 1] < -600 & logc[i, 2] < -600){ addforc <- -max(logc[i, 1], logc[i, 2]) kpost[i, 1] <- k[i, 1] * exp(logc[i, 1] + addforc) / (k[i, 1] * exp(logc[i, 1] + addforc) + k[i, 2] * exp(logc[i, 2] + addforc))# sum for 2 exp. in denom. kpost[i, 2] <- k[i, 2] * exp(logc[i, 2] + addforc) / (k[i, 1] * exp(logc[i, 1] + addforc) + k[i, 2] * exp(logc[i, 2] + addforc)) } else { kpost[i, 1] <- k[i, 1] * cmat[i, 1] / (k[i, 1] * cmat[i, 1] + k[i, 2] * cmat[i, 2]) # sum for 2 experts in denominator kpost[i, 2] <- k[i, 2] * cmat[i, 2] / (k[i, 1] * cmat[i, 1] + k[i, 2] * cmat[i, 2]) } }
Let's view the output(s) side by side:
logc t(cmatSim) all.equal(logc, t(cmatSim))
cmat t(cmatOutexpSim) all.equal(cmat, t(cmatOutexpSim))
Note that these aren't equal for the one region that has to be adjusted (in this case it's GOM). This is because in Michail's code, he does the adjusting in the loop to create kpost
, whereas in my code I send the whole transformed matrix (cmatOutexpSim
) to the calcK
function.
kpost t(kmatSim) all.equal(kpost, t(kmatSim))
So we can rest assured that the algorithm is working, and can proceed with the main sampling.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.