In robschick/eliciteg:
The goal of this vignette is to take the data gathered from the expert elicitation, turn it into individual priors for each expert, and then sample from it with the algorithm written by Professor Jim Clark, at Duke University. The methods for splitting the individual beta distributions into a Dirichlet are detailed in the vignette: Split-Beta.Rmd. The methods for sampling from the mixture distribution are given in: test.Rmd.
Sampling From the Dirichlet
To set up the sampler, we need to assemble the priors and the data. First, the data. We take these from a previously completed run of the Gibbs sampler. What I did was simply load in the data, run updateStates() once, and save it out to an rdata file:
femEx <- moveFem[, , k]
maleEx <- moveMale[, , k]
save(femEx, maleEx, file = 'data/moveData.rda')
Testing the Algorithm
In this section we will explain the sampling algorithm we are using within the Gibbs sampler, highlight how we can extend it to work with the mixture of Dirichlet's, and then test it with just the data.
Structure of the Algorithm
Jim wrote the algorithm as part of the Schick et al. 2013a PLoS-ONE paper. We used Gibbs sampling to estimate the parameters as opposed to drawing directly from the posterior because we impute the missing data for location. Therefore each time through the Gibbs loop, missing states can change. So the algorithm (based on Gelman et al. 2014) looks like this:
- Draw from $x \sim Gamma(a, 1)$
- Set $\theta = \frac{x}{\sum_{i=1}^n x}$
This is the implementation of that in R:
di <- matrix(rgamma(nreg * nreg,
shape = malePrior[, , k] + maleDat[, , k],
scale = 1),
nreg, nreg)
malePost[, , k] <- di / matrix(colSums(di), nreg, nreg, byrow = T)
where:
nreg equals the number of regions (9)malePrior[, , k] is the prior for the $k^{th}$ monthmaleDat[, , k] is the data for the $k^{th}$ monthmalePost[, , k] is the Dirichlet posterior
These components look like this. n.b. that in our original code, we chose different names for these data frames. Note also that the data is reinitialised and imputed each time through the loop, so we don't have a glimpse of it.
load('../data/eg_202_ng_25000_BIG_2500.rdata')
k <- 12
round(maleMove[, , k], 2) # prior
round(maleMoveProb[, , k], 2) # posterior
Sampling with the Mixture Approach
Ok, with the Dirichlet algorithm outlined, we need to use the sampling details from test.Rmd to simulate the posterior. Specifically, we need to set up a structure where by each time we go through the Gibbs loop, we sample from the mixture distribution. Before we had just one prior (as shown above); now the prior is a mixture of all experts' answers to the elicitation exercise.
The process flow will be as in test.Rmd. This is all within the gibbs loop (because the data update each time):
- Split the betas and assemble one Dirichlet matrix containing all priors (Prior to the Gibbs loop)
- Calculate
c for all experts using the estimate data and the individual prior
- Use that to assemble the $K$ probability vector
- At each iteration draw from a Multinomial distribution with probability $K$.
- Use the $K^{th}$ expert's Dirichlet distribution
- Merge it into the existing prior
- Sample from the Posterior
First thing that would be useful is to summarise the text of each question so we have it for reference. The reason I think this is necessary, is that we're going to have split this differently based on the direction of movement we ask about. I'll explain each question and data object that they refer to:
Question 9 - Consider a population of 100 adult female right whales that are present in the northern regions during December and move south to the mid-Atlantic during January.
In this scenario, there are two possible transitions. The animals are in one of 7 possible Northern regions, and can either move south to MIDA or remain in their Northern region. The answers are in dfList[[1]] (n.b. all the answers are different elements in this list):
library(eliciteg)
dfList[[1]]
The moves to the south are represented with the moveNortherntoMIDA column. The animals remaining in the Northern regions are in the moveNortherntoNorthern. In this case all of these transitions need to be split, because the first one represents 7 transitions from North to MIDA, and the second part of the question refers to animals remaining in 7 different Northern regions. So we'll have a matrix with 7 transitions from North to MIDA populated, as well as 7 North[i] to North[i] transitions.
Question 10 - Consider a population of 100 adult female right whales that are present in the mid-Atlantic during December and remain in the mid-Atlantic for any duration during January.
In this scenario, the animals are in the MIDA, and can make three transitions: 1) they can remain in MIDA; 2) move south to MIDA; or 3) move North. The answers are in dfList[[2]]:
dfList[[2]]
This one is a bit more straightforward. We have the MIDA to MIDA transition - moveMIDAtoMIDA - that does not need to be split. Then we have the MIDA to SEUS - moveMIDAtoSEUS - which also does not need to be split. Finally, we have the MIDA to Northern regions - moveMIDAtoNorthern - which does need to be split.
Question 11 - Consider a population of 100 adult female right whales that are present in the southeastern United States during December and move north to the mid-Atlantic during January.
In this scenario, the animals are in the SEUS and can make one of two transitions. The answers are in dfList[[3]].
dfList[[3]]
So far this is the most straightforward, as neither of the questions need to be split. We have animals starting in the SEUS, and can either remain - moveSEUStoSEUS - or move north to MIDA - moveSEUStoMIDA
Now we move on to the males. Question 12 - Consider a population of 100 adult male right whales that are present in the northern regions during December and move south to the mid-Atlantic during January.
As before with females, in this scenario, there are two possible transitions. The animals are in one of 7 possible Northern regions, and can either move south to MIDA or remain in their Northern region. The answers are in dfList[[4]]:
dfList[[4]]
The moves to the south are represented with the moveNortherntoMIDA column. The animals remaining in the Northern regions are in the moveNortherntoNorthern. In this case all of these transitions need to be split, because the first one represents 7 transitions from North to MIDA, and the second part of the question refers to animals remaining in 7 different Northern regions. So we'll have a matrix with 7 transitions from North to MIDA populated, as well as 7 North[i] to North[i] transitions.
Question 13 - Consider a population of 100 adult male right whales that are present in the mid-Atlantic during December and remain in the mid-Atlantic for any duration during January.
In this scenario, the animals are in the MIDA, and can make three transitions: 1) they can remain in MIDA; 2) move south to MIDA; or 3) move North. The answers are in dfList[[5]]:
dfList[[5]]
This one is a bit more straightforward. We have the MIDA to MIDA transition - moveMIDAtoMIDA - that does not need to be split. Then we have the MIDA to SEUS - moveMIDAtoSEUS - which also does not need to be split. Finally, we have the MIDA to Northern regions - moveMIDAtoNorthern - which does need to be split.
Question 14 - Consider a population of 100 adult male right whales that are present in the southeastern United States during December and move north to the mid-Atlantic during January.
In this scenario, the animals are in the SEUS and can make one of two transitions. The answers are in dfList[[6]].
dfList[[6]]
So far this is the most straightforward, as neither of the questions need to be split. We have animals starting in the SEUS, and can either remain - moveSEUStoSEUS - or move north to MIDA - moveSEUStoMIDA
Summary ok then, what all do we have for each of the 6 questions? For questions 9 and 12, all the answers need to be split. For example, we have for Question 9, expert 6: Beta(35, 4) and that will need to be split into Beta(5, 0.5714) for each of the 7 southern transitions. Similarly, we have Beta(65, 4), which will split to Beta(9.285714, 0.5714) for each of the 7 remaining transitions.
For questions 10 and 13, we need to split just one of the transitions - the moveMIDAtoNorthern answer. For question 10, expert 6, we'd go from Beta(20, 4) to Beta(2.8571429, 0.5714286) for transition from MIDA to each of the 7 Northern regions.
For questions 11 and 14, we do not need to split any transitions. They can be assembled as is.
Now that we've gone through that, let's start in on assembling the data into a form that we can use within the sampler. That means taking these individual Beta distributions, and assembling one Dirichlet.
Step 1 - Split the Betas
The data are stored in the dfList.rda data object. What we want, then, is a data structure that will allow us to easily extract from this and apply it to the matrix prior. Let's assume that we have once slice of a (9 * 9 * 12) array. If we ignore all of Philip Hamilton's priors, that will leave us with something like this:
regNorth <- regID[!regID %in% c('MIDA', 'SEUS')] # this vector contains the Northern regions
nNorthreg <- length(regNorth) # This will be used to divide down the Betas
movePrior <- matrix(1, nrow = length(regID), ncol = length(regID))
rownames(movePrior) <- colnames(movePrior) <- regID
movePrior2 <- movePrior
movePrior
However this raises a question - namely what do we want the flat prior values to be? I've started at 1, but perhaps it should be lower? One reason I ask this, is that as you'll see in the Northern to Northern transitions, you can end up with prior values << 1.
Back to our reality. Each time through the loop, we'll sample an expert, and populate movePrior. Let's show this for one expert for each of these two questions. Note that in the loop we'll sample the expert with the $c$ and $K$ machinery previously covered; here I'll just choose one for illustration purposes.
Recall that we used the following parameterisation for each Beta distribution: $\alpha = n \times m_1$ and $\beta = n - n \times m_1$, where $n$ = inputConfRaw, and $m_1$ = moveNortherntoMIDA / 100. Obviously, the column corresponding to $m_1$ changes for each question.
Here are two functions to return these values per each of the nreg transitions:
betaA <- function(n, m, nreg) {
expBeta1a <- (n * m) / nreg
expBeta1a
}
betaB <- function(n, m, nreg) {
expBeta1b <- (n - (n * m)) / nreg
expBeta1b
}
First we'll update elements of movePrior with data for the transitions south (North to MIDA). As noted above, the idea is that animals could be in any of the 7 Northern regions and could make a transition to MIDA:
expert <- 4
dfList[[1]][which(dfList[[1]]$expert == expert),]
n <- dfList[[1]]$inputConfRaw[which(dfList[[1]]$expert == expert)]
m <- dfList[[1]]$moveNortherntoMIDA[which(dfList[[1]]$expert == expert)] / 100
expBeta1a <- betaA(n, m, nreg = nNorthreg)
expBeta1b <- betaB(n, m, nreg = nNorthreg)
movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) %in% regNorth] <- expBeta1a
movePrior
And then adding in the North to North transitions:
n <- dfList[[1]]$inputConfRaw[which(dfList[[1]]$expert == expert)]
m <- dfList[[1]]$moveNortherntoNorthern[which(dfList[[1]]$expert == expert)] / 100
expBeta1a <- betaA(n, m, nreg = nNorthreg)
expBeta1b <- betaB(n, m, nreg = nNorthreg)
diag(movePrior)[which(rownames(movePrior) %in% regNorth)] <- expBeta1a
movePrior
OK, that was dfList[[1]], i.e. Question 9, what does the process and result look like for dfList[[2]], i.e. Question 10? Recall from above that we need to split the MIDA to Northern transitions moveMIDAtoNorthern, because it could go from MIDA to each of 7 areas in the North.
First for the remaining in MIDA transition (i.e., MIDA to MIDA):
expert2 <- 1
dfList[[2]][which(dfList[[2]]$expert == expert2),]
n <- dfList[[2]]$inputConfRaw[which(dfList[[2]]$expert == expert2)]
m <- dfList[[2]]$moveMIDAtoMIDA[which(dfList[[2]]$expert == expert2)] / 100
expBeta1a <- betaA(n, m, nreg = 1)
expBeta1b <- betaB(n, m, nreg = 1)
movePrior2[rownames(movePrior2) == 'MIDA', rownames(movePrior2) == 'MIDA'] <- expBeta1a
movePrior2
And then adding in the MIDA to SEUS transition:
n <- dfList[[2]]$inputConfRaw[which(dfList[[2]]$expert == expert2)]
m <- dfList[[2]]$moveMIDAtoSEUS[which(dfList[[2]]$expert == expert2)] / 100
expBeta1a <-betaA(n, m, nreg = 1)
expBeta1b <- betaB(n, m, nreg = 1)
movePrior2[rownames(movePrior2) == 'SEUS', colnames(movePrior2) == 'MIDA'] <- expBeta1a
movePrior2
And finally, the ones that need to be split:
n <- dfList[[2]]$inputConfRaw[which(dfList[[2]]$expert == expert2)]
m <- dfList[[2]]$moveMIDAtoNorthern[which(dfList[[2]]$expert == expert2)] / 100
expBeta1a <- betaA(n, m, nreg = nNorthreg)
expBeta1b <- betaB(n, m, nreg = nNorthreg)
movePrior2[rownames(movePrior2) %in% regNorth, colnames(movePrior2) == 'MIDA'] <- expBeta1a
movePrior2
Note above that the splitting we assign the split value to all of the regions in regNorth
Let's examine them together to see if the differences make sense given the structure of the questions.
dfList[[1]][dfList[[1]]$expert == expert, ]
movePrior
dfList[[2]][dfList[[2]]$expert == expert2, ]
movePrior2
So this seems to make sense to me. For expert 4 in the first answer, they had a majority of whales staying in the northern regions (30 going south, 70 remaining). So all of the diagonals are 0.7, and indeed higher than the transitions to MIDA. The matrix is set up as from columns to rows. So the MIDA row has been populated with 0.3 for the from Northern to MIDA transitions. Their confidence in their answers was low (7).
In contrast for expert 1 in the next question, they had much higher confidence - so we expect the values in the prior to be higher, and indeed they are. Also, the individual has equal probability of MIDA to MIDA as they do with MIDA to SEUS; those values are identical at 4.29. Since the number of whales must sum to 100, the value moving North is smaller (20), and that's reflected in the matrix.
Ok, so it at least passes the intuition test. Now we need to assemble all of these into one matrix/array/list. This will be static because these priors never change. However, they should be in the same form as the data, which is an array, because I'm going to have to add them together. I'm not 100% sure it matters all that much, because I'm just going to be sending them both to calcC().
Step 2 - Assemble the Prior Matrix
If we use an array, which is fine, than what are the dimensions that we need? We have:
- one month's transition (December to January),
- 9 regions
- 8 experts
- 3 questions
- 2 sexes (M/F)
We need to add a 9x9 slice to a 9x9 slice for each question, where above each array slice is a month. So here we'll just have one slice for the December to January transition. So for each question (1-6) we'll have 8 slices, where each slice represents one experts' prior that has been split based on their answers. I could then just have a recursive list, where each question has a sub-list containing each expert's answers, or each questions' answers. If the latter, we'd have list[[1]] correspond to question 9, and then within that list component, we'd have 8 sub-lists, each containing a 9*9 slice that corresponded to each expert.
What would that look like? Let's try it with fake-ish data:
nexp <- 8
nreg <- 9
nq <- 6
allqs <- list(q1 = vector('list', nexp),
q2 = vector('list', nexp),
q3 = vector('list', nexp),
q4 = vector('list', nexp),
q5 = vector('list', nexp),
q6 = vector('list', nexp))
for(i in seq_along(1:nq)){
for(j in seq_along(1:nexp)){
allqs[[i]][[j]] <- matrix(data = round(runif(nreg * nreg), 0) + j, ncol = nreg, nrow = nreg)
}
}
allqs[[2]][[2]]
Not the most elegant way of doing it, I don't think, but it works and is fast enough for this application. Ok, let's populate it with real data now:
# This chunk of code will iterate over the question/expert combinations and populate a list element with the new prior
nexp <- 8
nreg <- 9
nq <- 6
flatpriorVal <- 0.001
allqs <- list(q1 = vector('list', nexp),
q2 = vector('list', nexp),
q3 = vector('list', nexp),
q4 = vector('list', nexp),
q5 = vector('list', nexp),
q6 = vector('list', nexp))
for(i in seq_along(1:nq)){
for(j in seq_along(1:nexp)){
allqs[[i]][[j]] <- matrix(data = flatpriorVal, ncol = nreg, nrow = nreg)
rownames(allqs[[i]][[j]]) <- colnames(allqs[[i]][[j]]) <- regID
}
}
Ok, the list now contains the flat prior, and we can populate it. The goal of this will be to iterate over each question, peel off the answers, and the split them up and put them into the correct list element. To be anal for now, I'm going to split them into two code blocks for females (q9-11) and males (q12-14).
Step 2 - Assemble Prior Matrix - 2 element list
I'm going to redo this now as a two element list, instead of a 6 element list. Let's start with the data as before. The logic will be that the priors all start at 0, and then we will fill these prior values up as we did before, one transition type at a time. Then we'll add the three slices together. That's the idea anyway.
# Question 9 - requires splitting: 1) northern to mida; 2) northern to northern
listElement <- 1
expert <- dfList[[listElement]]$expert
nexp <- length(expert)
for(j in seq_along(1:nexp)){
movePrior <- allqs[[listElement]][[j]]
dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),]
n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])]
m <- dfList[[listElement]]$moveNortherntoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - Northern to MIDA
m2 <- dfList[[listElement]]$moveNortherntoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- Northern to Northern
expBeta1a <- betaA(n, m, nreg = nNorthreg)
expBeta1a2 <- betaA(n, m2, nreg = nNorthreg * nNorthreg)
movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) %in% regNorth] <- expBeta1a
movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) %in% regNorth] <- expBeta1a2
allqs[[listElement]][[j]] <- movePrior
}
# Question 10 - requires splitting MIDA to Northern transitions `moveMIDAtoNorthern`
listElement <- 2
expert <- dfList[[listElement]]$expert
nexp <- length(expert)
for(j in seq_along(1:nexp)){
movePrior <- allqs[[listElement]][[j]]
n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])]
m <- dfList[[listElement]]$moveMIDAtoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100
expBeta1a <- betaA(n, m, nreg = 1)
movePrior[rownames(movePrior) == 'MIDA', rownames(movePrior) == 'MIDA'] <- expBeta1a
m2 <- dfList[[listElement]]$moveMIDAtoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100
expBeta1a2 <-betaA(n, m2, nreg = 1)
movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'MIDA'] <- expBeta1a2
m3 <- dfList[[listElement]]$moveMIDAtoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100
expBeta1a3 <- betaA(n, m3, nreg = nNorthreg)
movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) == 'MIDA'] <- expBeta1a3
allqs[[listElement]][[j]] <- movePrior
}
# Question 11 - doesn't require splitting
listElement <- 3
expert <- dfList[[listElement]]$expert
nexp <- length(expert)
for(j in seq_along(1:nexp)){
movePrior <- allqs[[listElement]][[j]]
dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),]
n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])]
m <- dfList[[listElement]]$moveSEUStoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - SEUS to MIDA
m2 <- dfList[[listElement]]$moveSEUStoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- SEUS to SEUS
expBeta1a <- betaA(n, m, nreg = 1)
expBeta1a2 <- betaA(n, m2, nreg = 1)
movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) == 'SEUS'] <- expBeta1a
movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'SEUS'] <- expBeta1a2
allqs[[listElement]][[j]] <- movePrior
}
# Question 12 - requires splitting: 1) northern to mida; 2) northern to northern
listElement <- 4
expert <- dfList[[listElement]]$expert
nexp <- length(expert)
for(j in seq_along(1:nexp)){
movePrior <- allqs[[listElement]][[j]]
dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),]
n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])]
m <- dfList[[listElement]]$moveNortherntoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - Northern to MIDA
m2 <- dfList[[listElement]]$moveNortherntoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- Northern to Northern
expBeta1a <- betaA(n, m, nreg = nNorthreg)
expBeta1a2 <- betaA(n, m2, nreg = nNorthreg * nNorthreg)
movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) %in% regNorth] <- expBeta1a
movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) %in% regNorth] <- expBeta1a2
allqs[[listElement]][[j]] <- movePrior
}
# Question 13 - requires splitting MIDA to Northern transitions `moveMIDAtoNorthern`
listElement <- 5
expert <- dfList[[listElement]]$expert
nexp <- length(expert)
for(j in seq_along(1:nexp)){
movePrior <- allqs[[listElement]][[j]]
n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])]
m <- dfList[[listElement]]$moveMIDAtoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100
expBeta1a <- betaA(n, m, nreg = 1)
movePrior[rownames(movePrior) == 'MIDA', rownames(movePrior) == 'MIDA'] <- expBeta1a
m2 <- dfList[[listElement]]$moveMIDAtoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100
expBeta1a2 <-betaA(n, m2, nreg = 1)
movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'MIDA'] <- expBeta1a2
m3 <- dfList[[listElement]]$moveMIDAtoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100
expBeta1a3 <- betaA(n, m3, nreg = nNorthreg)
movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) == 'MIDA'] <- expBeta1a3
allqs[[listElement]][[j]] <- movePrior
}
# Question 14 - doesn't require splitting
listElement <- 6
expert <- dfList[[listElement]]$expert
nexp <- length(expert)
for(j in seq_along(1:nexp)){
movePrior <- allqs[[listElement]][[j]]
dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),]
n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])]
m <- dfList[[listElement]]$moveSEUStoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - SEUS to MIDA
m2 <- dfList[[listElement]]$moveSEUStoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- SEUS to SEUS
expBeta1a <- betaA(n, m, nreg = 1)
expBeta1a2 <- betaA(n, m2, nreg = 1)
movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) == 'SEUS'] <- expBeta1a
movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'SEUS'] <- expBeta1a2
allqs[[listElement]][[j]] <- movePrior
}
No we take these split-out elements and add them together to form a 2 element recursive list.
priorList <- list()
priorList$females$exp1 <- allqs[[1]][[1]] + allqs[[2]][[1]] + allqs[[3]][[1]]
priorList$females$exp2 <- allqs[[1]][[2]] + allqs[[2]][[2]] + allqs[[3]][[2]]
priorList$females$exp3 <- allqs[[1]][[3]] + allqs[[2]][[3]] + allqs[[3]][[3]]
priorList$females$exp4 <- allqs[[1]][[4]] + allqs[[2]][[4]] + allqs[[3]][[4]]
priorList$females$exp5 <- allqs[[1]][[5]] + allqs[[2]][[5]] + allqs[[3]][[5]]
priorList$females$exp6 <- allqs[[1]][[6]] + allqs[[2]][[6]] + allqs[[3]][[6]]
priorList$females$exp7 <- allqs[[1]][[7]] + allqs[[2]][[7]] + allqs[[3]][[7]]
priorList$females$exp8 <- allqs[[1]][[8]] + allqs[[2]][[8]] + allqs[[3]][[8]]
priorList$males$exp1 <- allqs[[4]][[1]] + allqs[[5]][[1]] + allqs[[6]][[1]]
priorList$males$exp2 <- allqs[[4]][[2]] + allqs[[5]][[2]] + allqs[[6]][[2]]
priorList$males$exp3 <- allqs[[4]][[3]] + allqs[[5]][[3]] + allqs[[6]][[3]]
priorList$males$exp4 <- allqs[[4]][[4]] + allqs[[5]][[4]] + allqs[[6]][[4]]
priorList$males$exp5 <- allqs[[4]][[5]] + allqs[[5]][[5]] + allqs[[6]][[5]]
priorList$males$exp6 <- allqs[[4]][[6]] + allqs[[5]][[6]] + allqs[[6]][[6]]
priorList$males$exp7 <- allqs[[4]][[7]] + allqs[[5]][[7]] + allqs[[6]][[7]]
priorList$males$exp8 <- allqs[[4]][[8]] + allqs[[5]][[8]] + allqs[[6]][[8]]
devtools::use_data(priorList, overwrite = TRUE)
robschick/eliciteg documentation built on May 27, 2019, 11:58 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
The goal of this vignette is to take the data gathered from the expert elicitation, turn it into individual priors for each expert, and then sample from it with the algorithm written by Professor Jim Clark, at Duke University. The methods for splitting the individual beta distributions into a Dirichlet are detailed in the vignette: Split-Beta.Rmd. The methods for sampling from the mixture distribution are given in: test.Rmd.
Sampling From the Dirichlet
To set up the sampler, we need to assemble the priors and the data. First, the data. We take these from a previously completed run of the Gibbs sampler. What I did was simply load in the data, run updateStates() once, and save it out to an rdata file:
femEx <- moveFem[, , k] maleEx <- moveMale[, , k] save(femEx, maleEx, file = 'data/moveData.rda')
Testing the Algorithm
In this section we will explain the sampling algorithm we are using within the Gibbs sampler, highlight how we can extend it to work with the mixture of Dirichlet's, and then test it with just the data.
Structure of the Algorithm
Jim wrote the algorithm as part of the Schick et al. 2013a PLoS-ONE paper. We used Gibbs sampling to estimate the parameters as opposed to drawing directly from the posterior because we impute the missing data for location. Therefore each time through the Gibbs loop, missing states can change. So the algorithm (based on Gelman et al. 2014) looks like this:
- Draw from $x \sim Gamma(a, 1)$
- Set $\theta = \frac{x}{\sum_{i=1}^n x}$
This is the implementation of that in R:
di <- matrix(rgamma(nreg * nreg, shape = malePrior[, , k] + maleDat[, , k], scale = 1), nreg, nreg) malePost[, , k] <- di / matrix(colSums(di), nreg, nreg, byrow = T)
where:
nregequals the number of regions (9)malePrior[, , k]is the prior for the $k^{th}$ monthmaleDat[, , k]is the data for the $k^{th}$ monthmalePost[, , k]is the Dirichlet posterior
These components look like this. n.b. that in our original code, we chose different names for these data frames. Note also that the data is reinitialised and imputed each time through the loop, so we don't have a glimpse of it.
load('../data/eg_202_ng_25000_BIG_2500.rdata') k <- 12 round(maleMove[, , k], 2) # prior round(maleMoveProb[, , k], 2) # posterior
Sampling with the Mixture Approach
Ok, with the Dirichlet algorithm outlined, we need to use the sampling details from test.Rmd to simulate the posterior. Specifically, we need to set up a structure where by each time we go through the Gibbs loop, we sample from the mixture distribution. Before we had just one prior (as shown above); now the prior is a mixture of all experts' answers to the elicitation exercise.
The process flow will be as in test.Rmd. This is all within the gibbs loop (because the data update each time):
- Split the betas and assemble one Dirichlet matrix containing all priors (Prior to the Gibbs loop)
- Calculate
cfor all experts using the estimate data and the individual prior - Use that to assemble the $K$ probability vector
- At each iteration draw from a Multinomial distribution with probability $K$.
- Use the $K^{th}$ expert's Dirichlet distribution
- Merge it into the existing prior
- Sample from the Posterior
First thing that would be useful is to summarise the text of each question so we have it for reference. The reason I think this is necessary, is that we're going to have split this differently based on the direction of movement we ask about. I'll explain each question and data object that they refer to:
Question 9 - Consider a population of 100 adult female right whales that are present in the northern regions during December and move south to the mid-Atlantic during January.
In this scenario, there are two possible transitions. The animals are in one of 7 possible Northern regions, and can either move south to MIDA or remain in their Northern region. The answers are in dfList[[1]] (n.b. all the answers are different elements in this list):
library(eliciteg)
dfList[[1]]
The moves to the south are represented with the moveNortherntoMIDA column. The animals remaining in the Northern regions are in the moveNortherntoNorthern. In this case all of these transitions need to be split, because the first one represents 7 transitions from North to MIDA, and the second part of the question refers to animals remaining in 7 different Northern regions. So we'll have a matrix with 7 transitions from North to MIDA populated, as well as 7 North[i] to North[i] transitions.
Question 10 - Consider a population of 100 adult female right whales that are present in the mid-Atlantic during December and remain in the mid-Atlantic for any duration during January.
In this scenario, the animals are in the MIDA, and can make three transitions: 1) they can remain in MIDA; 2) move south to MIDA; or 3) move North. The answers are in dfList[[2]]:
dfList[[2]]
This one is a bit more straightforward. We have the MIDA to MIDA transition - moveMIDAtoMIDA - that does not need to be split. Then we have the MIDA to SEUS - moveMIDAtoSEUS - which also does not need to be split. Finally, we have the MIDA to Northern regions - moveMIDAtoNorthern - which does need to be split.
Question 11 - Consider a population of 100 adult female right whales that are present in the southeastern United States during December and move north to the mid-Atlantic during January.
In this scenario, the animals are in the SEUS and can make one of two transitions. The answers are in dfList[[3]].
dfList[[3]]
So far this is the most straightforward, as neither of the questions need to be split. We have animals starting in the SEUS, and can either remain - moveSEUStoSEUS - or move north to MIDA - moveSEUStoMIDA
Now we move on to the males. Question 12 - Consider a population of 100 adult male right whales that are present in the northern regions during December and move south to the mid-Atlantic during January.
As before with females, in this scenario, there are two possible transitions. The animals are in one of 7 possible Northern regions, and can either move south to MIDA or remain in their Northern region. The answers are in dfList[[4]]:
dfList[[4]]
The moves to the south are represented with the moveNortherntoMIDA column. The animals remaining in the Northern regions are in the moveNortherntoNorthern. In this case all of these transitions need to be split, because the first one represents 7 transitions from North to MIDA, and the second part of the question refers to animals remaining in 7 different Northern regions. So we'll have a matrix with 7 transitions from North to MIDA populated, as well as 7 North[i] to North[i] transitions.
Question 13 - Consider a population of 100 adult male right whales that are present in the mid-Atlantic during December and remain in the mid-Atlantic for any duration during January.
In this scenario, the animals are in the MIDA, and can make three transitions: 1) they can remain in MIDA; 2) move south to MIDA; or 3) move North. The answers are in dfList[[5]]:
dfList[[5]]
This one is a bit more straightforward. We have the MIDA to MIDA transition - moveMIDAtoMIDA - that does not need to be split. Then we have the MIDA to SEUS - moveMIDAtoSEUS - which also does not need to be split. Finally, we have the MIDA to Northern regions - moveMIDAtoNorthern - which does need to be split.
Question 14 - Consider a population of 100 adult male right whales that are present in the southeastern United States during December and move north to the mid-Atlantic during January.
In this scenario, the animals are in the SEUS and can make one of two transitions. The answers are in dfList[[6]].
dfList[[6]]
So far this is the most straightforward, as neither of the questions need to be split. We have animals starting in the SEUS, and can either remain - moveSEUStoSEUS - or move north to MIDA - moveSEUStoMIDA
Summary ok then, what all do we have for each of the 6 questions? For questions 9 and 12, all the answers need to be split. For example, we have for Question 9, expert 6: Beta(35, 4) and that will need to be split into Beta(5, 0.5714) for each of the 7 southern transitions. Similarly, we have Beta(65, 4), which will split to Beta(9.285714, 0.5714) for each of the 7 remaining transitions.
For questions 10 and 13, we need to split just one of the transitions - the moveMIDAtoNorthern answer. For question 10, expert 6, we'd go from Beta(20, 4) to Beta(2.8571429, 0.5714286) for transition from MIDA to each of the 7 Northern regions.
For questions 11 and 14, we do not need to split any transitions. They can be assembled as is.
Now that we've gone through that, let's start in on assembling the data into a form that we can use within the sampler. That means taking these individual Beta distributions, and assembling one Dirichlet.
Step 1 - Split the Betas
The data are stored in the dfList.rda data object. What we want, then, is a data structure that will allow us to easily extract from this and apply it to the matrix prior. Let's assume that we have once slice of a (9 * 9 * 12) array. If we ignore all of Philip Hamilton's priors, that will leave us with something like this:
regNorth <- regID[!regID %in% c('MIDA', 'SEUS')] # this vector contains the Northern regions nNorthreg <- length(regNorth) # This will be used to divide down the Betas movePrior <- matrix(1, nrow = length(regID), ncol = length(regID)) rownames(movePrior) <- colnames(movePrior) <- regID movePrior2 <- movePrior movePrior
However this raises a question - namely what do we want the flat prior values to be? I've started at 1, but perhaps it should be lower? One reason I ask this, is that as you'll see in the Northern to Northern transitions, you can end up with prior values << 1.
Back to our reality. Each time through the loop, we'll sample an expert, and populate movePrior. Let's show this for one expert for each of these two questions. Note that in the loop we'll sample the expert with the $c$ and $K$ machinery previously covered; here I'll just choose one for illustration purposes.
Recall that we used the following parameterisation for each Beta distribution: $\alpha = n \times m_1$ and $\beta = n - n \times m_1$, where $n$ = inputConfRaw, and $m_1$ = moveNortherntoMIDA / 100. Obviously, the column corresponding to $m_1$ changes for each question.
Here are two functions to return these values per each of the nreg transitions:
betaA <- function(n, m, nreg) { expBeta1a <- (n * m) / nreg expBeta1a } betaB <- function(n, m, nreg) { expBeta1b <- (n - (n * m)) / nreg expBeta1b }
First we'll update elements of movePrior with data for the transitions south (North to MIDA). As noted above, the idea is that animals could be in any of the 7 Northern regions and could make a transition to MIDA:
expert <- 4 dfList[[1]][which(dfList[[1]]$expert == expert),] n <- dfList[[1]]$inputConfRaw[which(dfList[[1]]$expert == expert)] m <- dfList[[1]]$moveNortherntoMIDA[which(dfList[[1]]$expert == expert)] / 100 expBeta1a <- betaA(n, m, nreg = nNorthreg) expBeta1b <- betaB(n, m, nreg = nNorthreg) movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) %in% regNorth] <- expBeta1a movePrior
And then adding in the North to North transitions:
n <- dfList[[1]]$inputConfRaw[which(dfList[[1]]$expert == expert)] m <- dfList[[1]]$moveNortherntoNorthern[which(dfList[[1]]$expert == expert)] / 100 expBeta1a <- betaA(n, m, nreg = nNorthreg) expBeta1b <- betaB(n, m, nreg = nNorthreg) diag(movePrior)[which(rownames(movePrior) %in% regNorth)] <- expBeta1a movePrior
OK, that was dfList[[1]], i.e. Question 9, what does the process and result look like for dfList[[2]], i.e. Question 10? Recall from above that we need to split the MIDA to Northern transitions moveMIDAtoNorthern, because it could go from MIDA to each of 7 areas in the North.
First for the remaining in MIDA transition (i.e., MIDA to MIDA):
expert2 <- 1 dfList[[2]][which(dfList[[2]]$expert == expert2),] n <- dfList[[2]]$inputConfRaw[which(dfList[[2]]$expert == expert2)] m <- dfList[[2]]$moveMIDAtoMIDA[which(dfList[[2]]$expert == expert2)] / 100 expBeta1a <- betaA(n, m, nreg = 1) expBeta1b <- betaB(n, m, nreg = 1) movePrior2[rownames(movePrior2) == 'MIDA', rownames(movePrior2) == 'MIDA'] <- expBeta1a movePrior2
And then adding in the MIDA to SEUS transition:
n <- dfList[[2]]$inputConfRaw[which(dfList[[2]]$expert == expert2)] m <- dfList[[2]]$moveMIDAtoSEUS[which(dfList[[2]]$expert == expert2)] / 100 expBeta1a <-betaA(n, m, nreg = 1) expBeta1b <- betaB(n, m, nreg = 1) movePrior2[rownames(movePrior2) == 'SEUS', colnames(movePrior2) == 'MIDA'] <- expBeta1a movePrior2
And finally, the ones that need to be split:
n <- dfList[[2]]$inputConfRaw[which(dfList[[2]]$expert == expert2)] m <- dfList[[2]]$moveMIDAtoNorthern[which(dfList[[2]]$expert == expert2)] / 100 expBeta1a <- betaA(n, m, nreg = nNorthreg) expBeta1b <- betaB(n, m, nreg = nNorthreg) movePrior2[rownames(movePrior2) %in% regNorth, colnames(movePrior2) == 'MIDA'] <- expBeta1a movePrior2
Note above that the splitting we assign the split value to all of the regions in regNorth
Let's examine them together to see if the differences make sense given the structure of the questions.
dfList[[1]][dfList[[1]]$expert == expert, ] movePrior dfList[[2]][dfList[[2]]$expert == expert2, ] movePrior2
So this seems to make sense to me. For expert 4 in the first answer, they had a majority of whales staying in the northern regions (30 going south, 70 remaining). So all of the diagonals are 0.7, and indeed higher than the transitions to MIDA. The matrix is set up as from columns to rows. So the MIDA row has been populated with 0.3 for the from Northern to MIDA transitions. Their confidence in their answers was low (7).
In contrast for expert 1 in the next question, they had much higher confidence - so we expect the values in the prior to be higher, and indeed they are. Also, the individual has equal probability of MIDA to MIDA as they do with MIDA to SEUS; those values are identical at 4.29. Since the number of whales must sum to 100, the value moving North is smaller (20), and that's reflected in the matrix.
Ok, so it at least passes the intuition test. Now we need to assemble all of these into one matrix/array/list. This will be static because these priors never change. However, they should be in the same form as the data, which is an array, because I'm going to have to add them together. I'm not 100% sure it matters all that much, because I'm just going to be sending them both to calcC().
Step 2 - Assemble the Prior Matrix
If we use an array, which is fine, than what are the dimensions that we need? We have:
- one month's transition (December to January),
- 9 regions
- 8 experts
- 3 questions
- 2 sexes (M/F)
We need to add a 9x9 slice to a 9x9 slice for each question, where above each array slice is a month. So here we'll just have one slice for the December to January transition. So for each question (1-6) we'll have 8 slices, where each slice represents one experts' prior that has been split based on their answers. I could then just have a recursive list, where each question has a sub-list containing each expert's answers, or each questions' answers. If the latter, we'd have list[[1]] correspond to question 9, and then within that list component, we'd have 8 sub-lists, each containing a 9*9 slice that corresponded to each expert.
What would that look like? Let's try it with fake-ish data:
nexp <- 8 nreg <- 9 nq <- 6 allqs <- list(q1 = vector('list', nexp), q2 = vector('list', nexp), q3 = vector('list', nexp), q4 = vector('list', nexp), q5 = vector('list', nexp), q6 = vector('list', nexp)) for(i in seq_along(1:nq)){ for(j in seq_along(1:nexp)){ allqs[[i]][[j]] <- matrix(data = round(runif(nreg * nreg), 0) + j, ncol = nreg, nrow = nreg) } } allqs[[2]][[2]]
Not the most elegant way of doing it, I don't think, but it works and is fast enough for this application. Ok, let's populate it with real data now:
# This chunk of code will iterate over the question/expert combinations and populate a list element with the new prior nexp <- 8 nreg <- 9 nq <- 6 flatpriorVal <- 0.001 allqs <- list(q1 = vector('list', nexp), q2 = vector('list', nexp), q3 = vector('list', nexp), q4 = vector('list', nexp), q5 = vector('list', nexp), q6 = vector('list', nexp)) for(i in seq_along(1:nq)){ for(j in seq_along(1:nexp)){ allqs[[i]][[j]] <- matrix(data = flatpriorVal, ncol = nreg, nrow = nreg) rownames(allqs[[i]][[j]]) <- colnames(allqs[[i]][[j]]) <- regID } }
Ok, the list now contains the flat prior, and we can populate it. The goal of this will be to iterate over each question, peel off the answers, and the split them up and put them into the correct list element. To be anal for now, I'm going to split them into two code blocks for females (q9-11) and males (q12-14).
Step 2 - Assemble Prior Matrix - 2 element list
I'm going to redo this now as a two element list, instead of a 6 element list. Let's start with the data as before. The logic will be that the priors all start at 0, and then we will fill these prior values up as we did before, one transition type at a time. Then we'll add the three slices together. That's the idea anyway.
# Question 9 - requires splitting: 1) northern to mida; 2) northern to northern listElement <- 1 expert <- dfList[[listElement]]$expert nexp <- length(expert) for(j in seq_along(1:nexp)){ movePrior <- allqs[[listElement]][[j]] dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),] n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])] m <- dfList[[listElement]]$moveNortherntoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - Northern to MIDA m2 <- dfList[[listElement]]$moveNortherntoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- Northern to Northern expBeta1a <- betaA(n, m, nreg = nNorthreg) expBeta1a2 <- betaA(n, m2, nreg = nNorthreg * nNorthreg) movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) %in% regNorth] <- expBeta1a movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) %in% regNorth] <- expBeta1a2 allqs[[listElement]][[j]] <- movePrior } # Question 10 - requires splitting MIDA to Northern transitions `moveMIDAtoNorthern` listElement <- 2 expert <- dfList[[listElement]]$expert nexp <- length(expert) for(j in seq_along(1:nexp)){ movePrior <- allqs[[listElement]][[j]] n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])] m <- dfList[[listElement]]$moveMIDAtoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 expBeta1a <- betaA(n, m, nreg = 1) movePrior[rownames(movePrior) == 'MIDA', rownames(movePrior) == 'MIDA'] <- expBeta1a m2 <- dfList[[listElement]]$moveMIDAtoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100 expBeta1a2 <-betaA(n, m2, nreg = 1) movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'MIDA'] <- expBeta1a2 m3 <- dfList[[listElement]]$moveMIDAtoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100 expBeta1a3 <- betaA(n, m3, nreg = nNorthreg) movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) == 'MIDA'] <- expBeta1a3 allqs[[listElement]][[j]] <- movePrior } # Question 11 - doesn't require splitting listElement <- 3 expert <- dfList[[listElement]]$expert nexp <- length(expert) for(j in seq_along(1:nexp)){ movePrior <- allqs[[listElement]][[j]] dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),] n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])] m <- dfList[[listElement]]$moveSEUStoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - SEUS to MIDA m2 <- dfList[[listElement]]$moveSEUStoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- SEUS to SEUS expBeta1a <- betaA(n, m, nreg = 1) expBeta1a2 <- betaA(n, m2, nreg = 1) movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) == 'SEUS'] <- expBeta1a movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'SEUS'] <- expBeta1a2 allqs[[listElement]][[j]] <- movePrior }
# Question 12 - requires splitting: 1) northern to mida; 2) northern to northern listElement <- 4 expert <- dfList[[listElement]]$expert nexp <- length(expert) for(j in seq_along(1:nexp)){ movePrior <- allqs[[listElement]][[j]] dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),] n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])] m <- dfList[[listElement]]$moveNortherntoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - Northern to MIDA m2 <- dfList[[listElement]]$moveNortherntoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- Northern to Northern expBeta1a <- betaA(n, m, nreg = nNorthreg) expBeta1a2 <- betaA(n, m2, nreg = nNorthreg * nNorthreg) movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) %in% regNorth] <- expBeta1a movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) %in% regNorth] <- expBeta1a2 allqs[[listElement]][[j]] <- movePrior } # Question 13 - requires splitting MIDA to Northern transitions `moveMIDAtoNorthern` listElement <- 5 expert <- dfList[[listElement]]$expert nexp <- length(expert) for(j in seq_along(1:nexp)){ movePrior <- allqs[[listElement]][[j]] n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])] m <- dfList[[listElement]]$moveMIDAtoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 expBeta1a <- betaA(n, m, nreg = 1) movePrior[rownames(movePrior) == 'MIDA', rownames(movePrior) == 'MIDA'] <- expBeta1a m2 <- dfList[[listElement]]$moveMIDAtoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100 expBeta1a2 <-betaA(n, m2, nreg = 1) movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'MIDA'] <- expBeta1a2 m3 <- dfList[[listElement]]$moveMIDAtoNorthern[which(dfList[[listElement]]$expert == expert[j])] / 100 expBeta1a3 <- betaA(n, m3, nreg = nNorthreg) movePrior[rownames(movePrior) %in% regNorth, colnames(movePrior) == 'MIDA'] <- expBeta1a3 allqs[[listElement]][[j]] <- movePrior } # Question 14 - doesn't require splitting listElement <- 6 expert <- dfList[[listElement]]$expert nexp <- length(expert) for(j in seq_along(1:nexp)){ movePrior <- allqs[[listElement]][[j]] dfList[[listElement]][which(dfList[[listElement]]$expert == expert[j]),] n <- dfList[[listElement]]$inputConfRaw[which(dfList[[listElement]]$expert == expert[j])] m <- dfList[[listElement]]$moveSEUStoMIDA[which(dfList[[listElement]]$expert == expert[j])] / 100 # part 1 - SEUS to MIDA m2 <- dfList[[listElement]]$moveSEUStoSEUS[which(dfList[[listElement]]$expert == expert[j])] / 100 # Part 2- SEUS to SEUS expBeta1a <- betaA(n, m, nreg = 1) expBeta1a2 <- betaA(n, m2, nreg = 1) movePrior[rownames(movePrior) == 'MIDA', colnames(movePrior) == 'SEUS'] <- expBeta1a movePrior[rownames(movePrior) == 'SEUS', colnames(movePrior) == 'SEUS'] <- expBeta1a2 allqs[[listElement]][[j]] <- movePrior }
No we take these split-out elements and add them together to form a 2 element recursive list.
priorList <- list() priorList$females$exp1 <- allqs[[1]][[1]] + allqs[[2]][[1]] + allqs[[3]][[1]] priorList$females$exp2 <- allqs[[1]][[2]] + allqs[[2]][[2]] + allqs[[3]][[2]] priorList$females$exp3 <- allqs[[1]][[3]] + allqs[[2]][[3]] + allqs[[3]][[3]] priorList$females$exp4 <- allqs[[1]][[4]] + allqs[[2]][[4]] + allqs[[3]][[4]] priorList$females$exp5 <- allqs[[1]][[5]] + allqs[[2]][[5]] + allqs[[3]][[5]] priorList$females$exp6 <- allqs[[1]][[6]] + allqs[[2]][[6]] + allqs[[3]][[6]] priorList$females$exp7 <- allqs[[1]][[7]] + allqs[[2]][[7]] + allqs[[3]][[7]] priorList$females$exp8 <- allqs[[1]][[8]] + allqs[[2]][[8]] + allqs[[3]][[8]] priorList$males$exp1 <- allqs[[4]][[1]] + allqs[[5]][[1]] + allqs[[6]][[1]] priorList$males$exp2 <- allqs[[4]][[2]] + allqs[[5]][[2]] + allqs[[6]][[2]] priorList$males$exp3 <- allqs[[4]][[3]] + allqs[[5]][[3]] + allqs[[6]][[3]] priorList$males$exp4 <- allqs[[4]][[4]] + allqs[[5]][[4]] + allqs[[6]][[4]] priorList$males$exp5 <- allqs[[4]][[5]] + allqs[[5]][[5]] + allqs[[6]][[5]] priorList$males$exp6 <- allqs[[4]][[6]] + allqs[[5]][[6]] + allqs[[6]][[6]] priorList$males$exp7 <- allqs[[4]][[7]] + allqs[[5]][[7]] + allqs[[6]][[7]] priorList$males$exp8 <- allqs[[4]][[8]] + allqs[[5]][[8]] + allqs[[6]][[8]] devtools::use_data(priorList, overwrite = TRUE)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.