In robschick/eliciteg:

knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.width = 8, fig.height = 6)

rm(list = ls())
library(eliciteg)

Goal of the Vignette

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:

1. Assemble movement data
2. Assemble priors from multiple individual experts
3. Calculate normalizing constant $C$
4. Error check the $C$ matrix for numerical stability
5. Exponentiate the $C$ matrix
6. Use the data and $C$ to calculate the $K$ priors
7. Sample from the Dirichlet to simulate the $K$ posteriors

The Data

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:

1. iterate across columns, and test if any values are smaller than a threshold (default is -600)
2. if any are true, then to calculate the numerical value to add; this value is the largest in the column
3. use sweep to add those column-wise to the matrix
4. return the exponent of the matrix

In use, the function produces this output:

cmatOutexp <- transformC(cmatOut, tval = -600)
cmatOutexp

Assemble $K$

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.

Sampling From the Dirichlet

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:

Priors

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)

Posteriors

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)$.

Real Data

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).

Males

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)

Females

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)

Appendix - Comparing Two Methods

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:

1. cmatOutexpSim
2. cmatSim
3. kmatSim
4. 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:

the raw c values

logc
t(cmatSim)
all.equal(logc, t(cmatSim))

The transformed c values

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.

The k values

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.

robschick/eliciteg documentation built on May 27, 2019, 11:58 a.m.