Setup details are described here, and the meta-post about these solutions is here.

Materials

The summmer course1 is based off of the second edition of Statistical Rethinking by Richard McElreath. This submission covers the following exercise questions:

  • Chapter 9
    • E{3,4,5,6}
    • M{1,2,3}
  • Chapter 11
    • E{1,2,3,4}
    • M{2,3,4,5,6,8}
  • Chapter 12
    • E{4}
    • H{1,2}

Packages

A colophon with details is provided at the end, but the following packages and theme parameters are used throughout.

1libsUsed<-c("tidyverse","tidybayes","orgutils","dagitty",
2            "rethinking","tidybayes.rethinking",
3            "ggplot2","kableExtra","dplyr","glue",
4            "latex2exp","data.table","printr","devtools")
5invisible(lapply(libsUsed, library, character.only = TRUE));
6theme_set(theme_grey(base_size=24))
7set.seed(1995)

Chapter IX: Markov Chain Monte Carlo

Easy Questions (Ch9)

9E3

Which sort of parameters can Hamiltonian Monte Carlo not handle? Can you explain why?

Solution

Hamiltonian Monte Carlo is derived by adding the concept of momentum which requires that the Hessian is non-negative, which in term requires a continuous smooth function. Thus HMC cannot handle discrete parameters by construction. More formally, the HMC requires a transform from the D-dimensional parameter space to a 2D-dimensional phase space cite:betancourtConceptualIntroductionHamiltonian2018.

9E4

Explain the difference between the effective number of samples, n_eff as calculated by Stan, and the actual number of samples.

Solution

We will invoke the precise definition of the effective sample size cite:betancourtConceptualIntroductionHamiltonian2018

\[ ESS = \frac{N}{1+2\sum_{l=1}^{\infty}\rho_{l}} \]

Where we note that \(\rho_{l}\) is the lag-l autocorrelation of \(f\) over the Markov chain (in time). In essence, this is the number of independent samples which have equivalent information of the posterior. This is relevant, because the samples from a Marko chain are sequentially correlated (autocorrelated).

9E5

Which value should Rhat approach, when a chain is sampling the posterior distribution correctly?

Solution

The literature cite:gelmanBayesianDataAnalysis2014 often cites a value of \(1.01\) for convergence. However, newer versions of Stan tend are documented to suggest \(1.05\) since they use newer formulations of the Rhat value cite:vehtariRanknormalizationFoldingLocalization2020. It should also be noted that cite:royConvergenceDiagnosticsMarkov2020 the Rhat value does not necessarily indicate convergence, it is not a necessary and sufficient condition, but a heuristic, and should be understood as such.

HOLD 9E6

Sketch a good trace plot for a Markov chain, one that is effectively sampling from the posterior distribution. What is good about its shape? Then sketch a trace plot for a malfunctioning Markov chain. What about its shape indicates malfunction?

Solution

Recall that the “health” of a chain can be determined by the following qualities in the trace plot.

Stationarity
This ensures that the chain is sampling the high probability portion of the posterior distribution
Mixing
This ensures that the chain explores the full region
Convergence
Convergence implies that independent chains agree on the same region of high probability

We will require a sample model to plot.

 1data(rugged)
 2rugDat<-rugged
 3rugDat<-rugDat %>% dplyr::mutate(logGDP=log(rgdppc_2000)) %>% tidyr::drop_na() %>% dplyr::mutate(logGDP_std=logGDP/mean(logGDP),
 4                                                                                                 rugged_std=rugged/max(rugged),
 5                                                                                                 cid=ifelse(cont_africa==1,1,2))
 6datList<-list(
 7  logGDP_std=rugDat$logGDP_std,
 8  rugged_std=rugDat$rugged_std,
 9  cid=as.integer(rugDat$cid)
10)
1m91unif<-ulam(
2  alist(
3    logGDP_std ~ dnorm(mu,sigma),
4    mu<-a[cid] + b[cid]*(rugged_std-0.215),
5    a[cid]~dnorm(1,0.1),
6    b[cid]~dnorm(0,0.3),
7    sigma~dunif(0,1)
8  ), data=datList, chains=4, cores=4
9)

We would like to check the trace and trace rank plots.

1m91unif %>% traceplot
1m91unif %>% trankplot

Clearly this is a good model, with well mixed chains, as can be seen in the trank and trace plots.

We will now check the plots for the unhealthy chain described in the chapter.

1m9e4un<-ulam(
2  alist(
3    y ~ dnorm(mu,sigma),
4    mu<-alpha,
5    alpha ~ dnorm(0,1000),
6    sigma~dexp(0.0001)
7  ),data=list(y=c(-1,1)),chains=4,cores=4
8)
  1SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 1).
  2Chain 1:
  3Chain 1: Gradient evaluation took 8.1e-05 seconds
  4Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.81 seconds.
  5Chain 1: Adjust your expectations accordingly!
  6Chain 1:
  7Chain 1:
  8Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
  9Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
 10
 11SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 2).
 12Chain 2:
 13Chain 2: Gradient evaluation took 3.6e-05 seconds
 14Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
 15Chain 2: Adjust your expectations accordingly!
 16Chain 2:
 17Chain 2:
 18Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
 19Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
 20Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
 21Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
 22Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
 23Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
 24Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
 25
 26SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 3).
 27Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
 28Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
 29Chain 3:
 30Chain 3: Gradient evaluation took 2.5e-05 seconds
 31Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
 32Chain 3: Adjust your expectations accordingly!
 33Chain 3:
 34Chain 3:
 35Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
 36Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
 37Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
 38Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
 39Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
 40Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
 41Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
 42Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
 43Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
 44Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
 45Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
 46Chain 2:
 47Chain 2:  Elapsed Time: 0.069084 seconds (Warm-up)
 48Chain 2:                0.023083 seconds (Sampling)
 49Chain 2:                0.092167 seconds (Total)
 50Chain 2:
 51Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
 52
 53SAMPLING FOR MODEL '726d002e27cec1633082261fcfedb813' NOW (CHAIN 4).
 54Chain 4:
 55Chain 4: Gradient evaluation took 3.1e-05 seconds
 56Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
 57Chain 4: Adjust your expectations accordingly!
 58Chain 4:
 59Chain 4:
 60Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
 61Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
 62Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
 63Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
 64Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
 65Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
 66Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
 67Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
 68Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
 69Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
 70Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
 71Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
 72Chain 1:
 73Chain 1:  Elapsed Time: 0.082116 seconds (Warm-up)
 74Chain 1:                0.095511 seconds (Sampling)
 75Chain 1:                0.177627 seconds (Total)
 76Chain 1:
 77Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
 78Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
 79Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
 80Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
 81Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
 82Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
 83Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
 84Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
 85Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
 86Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
 87Chain 3:
 88Chain 3:  Elapsed Time: 0.073606 seconds (Warm-up)
 89Chain 3:                0.076388 seconds (Sampling)
 90Chain 3:                0.149994 seconds (Total)
 91Chain 3:
 92Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
 93Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
 94Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
 95Chain 4:
 96Chain 4:  Elapsed Time: 0.064213 seconds (Warm-up)
 97Chain 4:                0.091667 seconds (Sampling)
 98Chain 4:                0.15588 seconds (Total)
 99Chain 4:
100Warning messages:
1011: There were 55 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See
102http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
1032: Examine the pairs() plot to diagnose sampling problems
104
1053: The largest R-hat is 1.08, indicating chains have not mixed.
106Running the chains for more iterations may help. See
107http://mc-stan.org/misc/warnings.html#r-hat
1084: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
109Running the chains for more iterations may help. See
110http://mc-stan.org/misc/warnings.html#bulk-ess
1115: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
112Running the chains for more iterations may help. See
113http://mc-stan.org/misc/warnings.html#tail-ess
1m9e4un %>% traceplot
1m9e4un %>% trankplot

Clearly these plots show a model which is unable to converge.

1m9e4un %>% precis
1        mean     sd    5.5%   94.5% n_eff Rhat4
2alpha -22.01 267.51 -414.95  316.61   154  1.02
3sigma 420.97 983.52    5.39 1894.87   151  1.04

This has clear repercussions on the actual predictions as well.

Questions of Medium Complexity (Ch9)

HOLD 9M1

Re-estimate the terrain ruggedness model from the chapter, but now using a uniform prior for the standard deviation, sigma. The uniform prior should be dunif(0,1). Use ulam to estimate the posterior. Does the different prior have any detectable influence on the posterior distribution of sigma? What or why not?

Solution

Instead of using the complete.cases formulation in the book, we will instead use a more tidyverse friendly approach.

 1data(rugged)
 2rugDat<-rugged
 3rugDat<-rugDat %>% dplyr::mutate(logGDP=log(rgdppc_2000)) %>% tidyr::drop_na() %>% dplyr::mutate(logGDP_std=logGDP/mean(logGDP),
 4                                                                                                 rugged_std=rugged/max(rugged),
 5                                                                                                 cid=ifelse(cont_africa==1,1,2))
 6datList<-list(
 7  logGDP_std=rugDat$logGDP_std,
 8  rugged_std=rugDat$rugged_std,
 9  cid=as.integer(rugDat$cid)
10)

We can now formulate a model with a uniform prior on sigma.

1m91unif<-ulam(
2  alist(
3    logGDP_std ~ dnorm(mu,sigma),
4    mu<-a[cid] + b[cid]*(rugged_std-0.215),
5    a[cid]~dnorm(1,0.1),
6    b[cid]~dnorm(0,0.3),
7    sigma~dunif(0,1)
8  ), data=datList, chains=4, cores=4
9)
1m91exp<-ulam(
2  alist(
3    logGDP_std ~ dnorm(mu,sigma),
4    mu<-a[cid] + b[cid]*(rugged_std-0.215),
5    a[cid]~dnorm(1,0.1),
6    b[cid]~dnorm(0,0.3),
7    sigma~dexp(1)
8  ), data=datList, chains=4, cores=4
9)
 1SAMPLING FOR MODEL '9b462775c5cc2badb2b667c53f2020c8' NOW (CHAIN 1).
 2Chain 1:
 3Chain 1: Gradient evaluation took 2.8e-05 seconds
 4Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
 5Chain 1: Adjust your expectations accordingly!
 6Chain 1:
 7Chain 1:
 8Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
 9
10SAMPLING FOR MODEL '9b462775c5cc2badb2b667c53f2020c8' NOW (CHAIN 2).
11Chain 2:
12Chain 2: Gradient evaluation took 3.3e-05 seconds
13Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
14Chain 2: Adjust your expectations accordingly!
15Chain 2:
16Chain 2:
17Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
18
19SAMPLING FOR MODEL '9b462775c5cc2badb2b667c53f2020c8' NOW (CHAIN 3).
20Chain 3:
21Chain 3: Gradient evaluation took 2.6e-05 seconds
22Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
23Chain 3: Adjust your expectations accordingly!
24Chain 3:
25Chain 3:
26Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
27Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
28Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
29Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
30Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
31
32SAMPLING FOR MODEL '9b462775c5cc2badb2b667c53f2020c8' NOW (CHAIN 4).
33Chain 4:
34Chain 4: Gradient evaluation took 2.5e-05 seconds
35Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
36Chain 4: Adjust your expectations accordingly!
37Chain 4:
38Chain 4:
39Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
40Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
41Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
42Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
43Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
44Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
45Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
46Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
47Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
48Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
49Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
50Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
51Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
52Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
53Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
54Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
55Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
56Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
57Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
58Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
59Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
60Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
61Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
62Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
63Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
64Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
65Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
66Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
67Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
68Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
69Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
70Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
71Chain 1:
72Chain 1:  Elapsed Time: 0.048459 seconds (Warm-up)
73Chain 1:                0.017728 seconds (Sampling)
74Chain 1:                0.066187 seconds (Total)
75Chain 1:
76Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
77Chain 2:
78Chain 2:  Elapsed Time: 0.043574 seconds (Warm-up)
79Chain 2:                0.018569 seconds (Sampling)
80Chain 2:                0.062143 seconds (Total)
81Chain 2:
82Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
83Chain 3:
84Chain 3:  Elapsed Time: 0.03201 seconds (Warm-up)
85Chain 3:                0.024666 seconds (Sampling)
86Chain 3:                0.056676 seconds (Total)
87Chain 3:
88Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
89Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
90Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
91Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
92Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
93Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
94Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
95Chain 4:
96Chain 4:  Elapsed Time: 0.040804 seconds (Warm-up)
97Chain 4:                0.018479 seconds (Sampling)
98Chain 4:                0.059283 seconds (Total)
99Chain 4:

The posterior distributions are simply:

1m91exp %>% extract.samples %>% .$sigma %>% dens(.,xlab="sigma")
2m91unif %>% extract.samples  %>% .$sigma %>% dens(.,add=TRUE,col="blue")
3mtext("posterior")

With the priors being:

1m91exp %>% extract.prior %>% .$sigma %>% dens(.,xlab="sigma")
2m91unif %>% extract.prior  %>% .$sigma %>% dens(.,add=TRUE,col="blue")
3mtext("prior")

This makes sense, since we know that uniform prior is essentially a step function between 0 and 1 with a value of 1, while the exponential function decays normally, but should actually be spiked upwards to 1 as well.

HOLD 9M2

Modify the terrain ruggedness model again. This times, change the prior for b[cid] to dexp(0.3). What does this do to the posterior distribution? Can you explain it?

Solution
1m92exp<-ulam(
2  alist(
3    logGDP_std ~ dnorm(mu,sigma),
4    mu<-a[cid] + b[cid]*(rugged_std-0.215),
5    a[cid]~dnorm(1,0.1),
6    b[cid]~dexp(0.3),
7    sigma~dexp(1)
8  ),data=datList, chains=4, cores=4
9)

Priors:

1m92exp %>% extract.prior %>% .$sigma %>% dens(.,xlab="sigma",col="blue")
2mtext("prior")

Posterior:

1m92exp %>% extract.samples %>% .$sigma %>% dens(.,xlab="sigma",col="blue")
2mtext("posterior")
1m92exp %>% precis(.,depth = 2)
2m91exp %>% precis(.,depth = 2)
3m91unif %>% precis(.,depth = 2)
 1      mean   sd 5.5% 94.5% n_eff Rhat4
 2a[1]  1.01 0.02 0.98  1.03  1161     1
 3b[1]  0.13 0.08 0.02  0.27   821     1
 4sigma 0.12 0.01 0.10  0.15  1097     1
 5
 6      mean   sd  5.5% 94.5% n_eff Rhat4
 7a[1]  1.01 0.02  0.98  1.04  1609     1
 8b[1]  0.09 0.10 -0.06  0.24  1518     1
 9sigma 0.12 0.01  0.10  0.15  1493     1
10
11      mean   sd  5.5% 94.5% n_eff Rhat4
12a[1]  1.01 0.02  0.97  1.04  1744     1
13b[1]  0.10 0.09 -0.05  0.25  1874     1
14sigma 0.12 0.01  0.10  0.15  1824     1

We can see that there isn’t much difference, however, the main difference is in the b parameter, which seems to have fewer samples, and is also no longer takes any negative values.

HOLD 9M3

Re-estimate one of the Stan models from the chapter, but at different numbers of warmup iterations. Be sure to use the same number of sampling iterations in each case. Compare the n_eff values. How much warmup is enough?

Solution

For brevity, we will re-use the same data and model as used in the previous questions.

1warmTrial<-seq.int(10,10000,length.out = 10)
2nSampleEff<-matrix(NA,nrow=length(warmTrial),ncol=3)
3nSampleEffExp<-matrix(NA,nrow=length(warmTrial),ncol=3)
  • Uniform Model

    1for(i in 1:length(warmTrial)){
    2  tmp<-ulam(m91unif,chains=4,cores=4,refresh=-1,warmup=warmTrial[i],iter=1000+warmTrial[i])
    3  nSampleEff[i,]<-precis(tmp,2)$n_eff
    4}
    
      1
      2Chain 1:
      3Chain 1: Gradient evaluation took 9.7e-05 seconds
      4Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.97 seconds.
      5Chain 1: Adjust your expectations accordingly!
      6Chain 1:
      7Chain 1:
      8Chain 1: WARNING: No variance estimation is
      9Chain 1:          performed for num_warmup < 20
     10Chain 1:
     11Chain 2:
     12Chain 2: Gradient evaluation took 9.6e-05 seconds
     13Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
     14Chain 2: Adjust your expectations accordingly!
     15Chain 2:
     16Chain 2:
     17Chain 2: WARNING: No variance estimation is
     18Chain 2:          performed for num_warmup < 20
     19Chain 2:
     20Chain 3:
     21Chain 3: Gradient evaluation took 0.000107 seconds
     22Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.07 seconds.
     23Chain 3: Adjust your expectations accordingly!
     24Chain 3:
     25Chain 3:
     26Chain 3: WARNING: No variance estimation is
     27Chain 3:          performed for num_warmup < 20
     28Chain 3:
     29Chain 1:
     30Chain 1:  Elapsed Time: 0.00162 seconds (Warm-up)
     31Chain 1:                0.1738 seconds (Sampling)
     32Chain 1:                0.17542 seconds (Total)
     33Chain 1:
     34Chain 3:
     35Chain 3:  Elapsed Time: 0.003755 seconds (Warm-up)
     36Chain 3:                0.059869 seconds (Sampling)
     37Chain 3:                0.063624 seconds (Total)
     38Chain 3:
     39Chain 4:
     40Chain 4: Gradient evaluation took 9.2e-05 seconds
     41Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.92 seconds.
     42Chain 4: Adjust your expectations accordingly!
     43Chain 4:
     44Chain 4:
     45Chain 4: WARNING: No variance estimation is
     46Chain 4:          performed for num_warmup < 20
     47Chain 4:
     48Chain 2:
     49Chain 2:  Elapsed Time: 0.003924 seconds (Warm-up)
     50Chain 2:                0.257815 seconds (Sampling)
     51Chain 2:                0.261739 seconds (Total)
     52Chain 2:
     53Chain 4:
     54Chain 4:  Elapsed Time: 0.003735 seconds (Warm-up)
     55Chain 4:                0.127719 seconds (Sampling)
     56Chain 4:                0.131454 seconds (Total)
     57Chain 4:
     58Chain 1:
     59Chain 1: Gradient evaluation took 0.000209 seconds
     60Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.09 seconds.
     61Chain 1: Adjust your expectations accordingly!
     62Chain 1:
     63Chain 1:
     64Chain 2:
     65Chain 2: Gradient evaluation took 7.8e-05 seconds
     66Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
     67Chain 2: Adjust your expectations accordingly!
     68Chain 2:
     69Chain 2:
     70Chain 3:
     71Chain 3: Gradient evaluation took 6.2e-05 seconds
     72Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
     73Chain 3: Adjust your expectations accordingly!
     74Chain 3:
     75Chain 3:
     76Chain 4:
     77Chain 4: Gradient evaluation took 8.2e-05 seconds
     78Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
     79Chain 4: Adjust your expectations accordingly!
     80Chain 4:
     81Chain 4:
     82Chain 1:
     83Chain 1:  Elapsed Time: 0.118609 seconds (Warm-up)
     84Chain 1:                0.075057 seconds (Sampling)
     85Chain 1:                0.193666 seconds (Total)
     86Chain 1:
     87Chain 2:
     88Chain 2:  Elapsed Time: 0.122167 seconds (Warm-up)
     89Chain 2:                0.05402 seconds (Sampling)
     90Chain 2:                0.176187 seconds (Total)
     91Chain 2:
     92Chain 3:
     93Chain 3:  Elapsed Time: 0.118156 seconds (Warm-up)
     94Chain 3:                0.035528 seconds (Sampling)
     95Chain 3:                0.153684 seconds (Total)
     96Chain 3:
     97Chain 4:
     98Chain 4:  Elapsed Time: 0.073505 seconds (Warm-up)
     99Chain 4:                0.040445 seconds (Sampling)
    100Chain 4:                0.11395 seconds (Total)
    101Chain 4:
    102Chain 1:
    103Chain 1: Gradient evaluation took 9.8e-05 seconds
    104Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.98 seconds.
    105Chain 1: Adjust your expectations accordingly!
    106Chain 1:
    107Chain 1:
    108Chain 2:
    109Chain 2: Gradient evaluation took 7.3e-05 seconds
    110Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.73 seconds.
    111Chain 2: Adjust your expectations accordingly!
    112Chain 2:
    113Chain 2:
    114Chain 3:
    115Chain 3: Gradient evaluation took 0.000107 seconds
    116Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.07 seconds.
    117Chain 3: Adjust your expectations accordingly!
    118Chain 3:
    119Chain 3:
    120Chain 4:
    121Chain 4: Gradient evaluation took 0.000109 seconds
    122Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
    123Chain 4: Adjust your expectations accordingly!
    124Chain 4:
    125Chain 4:
    126Chain 1:
    127Chain 1:  Elapsed Time: 0.23286 seconds (Warm-up)
    128Chain 1:                0.032903 seconds (Sampling)
    129Chain 1:                0.265763 seconds (Total)
    130Chain 1:
    131Chain 2:
    132Chain 2:  Elapsed Time: 0.196139 seconds (Warm-up)
    133Chain 2:                0.032946 seconds (Sampling)
    134Chain 2:                0.229085 seconds (Total)
    135Chain 2:
    136Chain 3:
    137Chain 3:  Elapsed Time: 0.15298 seconds (Warm-up)
    138Chain 3:                0.042238 seconds (Sampling)
    139Chain 3:                0.195218 seconds (Total)
    140Chain 3:
    141Chain 4:
    142Chain 4:  Elapsed Time: 0.109015 seconds (Warm-up)
    143Chain 4:                0.041119 seconds (Sampling)
    144Chain 4:                0.150134 seconds (Total)
    145Chain 4:
    146Chain 1:
    147Chain 1: Gradient evaluation took 7.4e-05 seconds
    148Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.74 seconds.
    149Chain 1: Adjust your expectations accordingly!
    150Chain 1:
    151Chain 1:
    152Chain 2:
    153Chain 2: Gradient evaluation took 7.9e-05 seconds
    154Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.79 seconds.
    155Chain 2: Adjust your expectations accordingly!
    156Chain 2:
    157Chain 2:
    158Chain 3:
    159Chain 3: Gradient evaluation took 0.000109 seconds
    160Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.09 seconds.
    161Chain 3: Adjust your expectations accordingly!
    162Chain 3:
    163Chain 3:
    164Chain 4:
    165Chain 4: Gradient evaluation took 0.000101 seconds
    166Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.01 seconds.
    167Chain 4: Adjust your expectations accordingly!
    168Chain 4:
    169Chain 4:
    170Chain 1:
    171Chain 1:  Elapsed Time: 0.304209 seconds (Warm-up)
    172Chain 1:                0.038004 seconds (Sampling)
    173Chain 1:                0.342213 seconds (Total)
    174Chain 1:
    175Chain 2:
    176Chain 2:  Elapsed Time: 0.274087 seconds (Warm-up)
    177Chain 2:                0.034831 seconds (Sampling)
    178Chain 2:                0.308918 seconds (Total)
    179Chain 2:
    180Chain 3:
    181Chain 3:  Elapsed Time: 0.231719 seconds (Warm-up)
    182Chain 3:                0.038131 seconds (Sampling)
    183Chain 3:                0.26985 seconds (Total)
    184Chain 3:
    185Chain 4:
    186Chain 4:  Elapsed Time: 0.177718 seconds (Warm-up)
    187Chain 4:                0.038546 seconds (Sampling)
    188Chain 4:                0.216264 seconds (Total)
    189Chain 4:
    190Chain 1:
    191Chain 1: Gradient evaluation took 0.000117 seconds
    192Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.17 seconds.
    193Chain 1: Adjust your expectations accordingly!
    194Chain 1:
    195Chain 1:
    196Chain 2:
    197Chain 2: Gradient evaluation took 7.3e-05 seconds
    198Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.73 seconds.
    199Chain 2: Adjust your expectations accordingly!
    200Chain 2:
    201Chain 2:
    202Chain 3:
    203Chain 3: Gradient evaluation took 7.2e-05 seconds
    204Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.72 seconds.
    205Chain 3: Adjust your expectations accordingly!
    206Chain 3:
    207Chain 3:
    208Chain 4:
    209Chain 4: Gradient evaluation took 7.3e-05 seconds
    210Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.73 seconds.
    211Chain 4: Adjust your expectations accordingly!
    212Chain 4:
    213Chain 4:
    214Chain 1:
    215Chain 1:  Elapsed Time: 0.322041 seconds (Warm-up)
    216Chain 1:                0.048995 seconds (Sampling)
    217Chain 1:                0.371036 seconds (Total)
    218Chain 1:
    219Chain 2:
    220Chain 2:  Elapsed Time: 0.293325 seconds (Warm-up)
    221Chain 2:                0.032541 seconds (Sampling)
    222Chain 2:                0.325866 seconds (Total)
    223Chain 2:
    224Chain 3:
    225Chain 3:  Elapsed Time: 0.264383 seconds (Warm-up)
    226Chain 3:                0.04051 seconds (Sampling)
    227Chain 3:                0.304893 seconds (Total)
    228Chain 3:
    229Chain 4:
    230Chain 4:  Elapsed Time: 0.220301 seconds (Warm-up)
    231Chain 4:                0.040218 seconds (Sampling)
    232Chain 4:                0.260519 seconds (Total)
    233Chain 4:
    234Chain 1:
    235Chain 1: Gradient evaluation took 4.9e-05 seconds
    236Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.49 seconds.
    237Chain 1: Adjust your expectations accordingly!
    238Chain 1:
    239Chain 1:
    240Chain 2:
    241Chain 2: Gradient evaluation took 5.1e-05 seconds
    242Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
    243Chain 2: Adjust your expectations accordingly!
    244Chain 2:
    245Chain 2:
    246Chain 3:
    247Chain 3: Gradient evaluation took 3.9e-05 seconds
    248Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
    249Chain 3: Adjust your expectations accordingly!
    250Chain 3:
    251Chain 3:
    252Chain 4:
    253Chain 4: Gradient evaluation took 3.9e-05 seconds
    254Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
    255Chain 4: Adjust your expectations accordingly!
    256Chain 4:
    257Chain 4:
    258Chain 1:
    259Chain 1:  Elapsed Time: 0.306918 seconds (Warm-up)
    260Chain 1:                0.048194 seconds (Sampling)
    261Chain 1:                0.355112 seconds (Total)
    262Chain 1:
    263Chain 2:
    264Chain 2:  Elapsed Time: 0.299256 seconds (Warm-up)
    265Chain 2:                0.052776 seconds (Sampling)
    266Chain 2:                0.352032 seconds (Total)
    267Chain 2:
    268Chain 3:
    269Chain 3:  Elapsed Time: 0.300093 seconds (Warm-up)
    270Chain 3:                0.052132 seconds (Sampling)
    271Chain 3:                0.352225 seconds (Total)
    272Chain 3:
    273Chain 4:
    274Chain 4:  Elapsed Time: 0.278897 seconds (Warm-up)
    275Chain 4:                0.053532 seconds (Sampling)
    276Chain 4:                0.332429 seconds (Total)
    277Chain 4:
    278Chain 1:
    279Chain 1: Gradient evaluation took 9.4e-05 seconds
    280Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.94 seconds.
    281Chain 1: Adjust your expectations accordingly!
    282Chain 1:
    283Chain 1:
    284Chain 2:
    285Chain 2: Gradient evaluation took 8.4e-05 seconds
    286Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.84 seconds.
    287Chain 2: Adjust your expectations accordingly!
    288Chain 2:
    289Chain 2:
    290Chain 3:
    291Chain 3: Gradient evaluation took 7.1e-05 seconds
    292Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.71 seconds.
    293Chain 3: Adjust your expectations accordingly!
    294Chain 3:
    295Chain 3:
    296Chain 4:
    297Chain 4: Gradient evaluation took 6.2e-05 seconds
    298Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.62 seconds.
    299Chain 4: Adjust your expectations accordingly!
    300Chain 4:
    301Chain 4:
    302Chain 1:
    303Chain 1:  Elapsed Time: 0.437519 seconds (Warm-up)
    304Chain 1:                0.04052 seconds (Sampling)
    305Chain 1:                0.478039 seconds (Total)
    306Chain 1:
    307Chain 3:
    308Chain 3:  Elapsed Time: 0.3479 seconds (Warm-up)
    309Chain 3:                0.033007 seconds (Sampling)
    310Chain 3:                0.380907 seconds (Total)
    311Chain 3:
    312Chain 2:
    313Chain 2:  Elapsed Time: 0.412403 seconds (Warm-up)
    314Chain 2:                0.052948 seconds (Sampling)
    315Chain 2:                0.465351 seconds (Total)
    316Chain 2:
    317Chain 4:
    318Chain 4:  Elapsed Time: 0.315601 seconds (Warm-up)
    319Chain 4:                0.038006 seconds (Sampling)
    320Chain 4:                0.353607 seconds (Total)
    321Chain 4:
    322Chain 1:
    323Chain 1: Gradient evaluation took 8.8e-05 seconds
    324Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.88 seconds.
    325Chain 1: Adjust your expectations accordingly!
    326Chain 1:
    327Chain 1:
    328Chain 2:
    329Chain 2: Gradient evaluation took 5.8e-05 seconds
    330Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
    331Chain 2: Adjust your expectations accordingly!
    332Chain 2:
    333Chain 2:
    334Chain 3:
    335Chain 3: Gradient evaluation took 5.9e-05 seconds
    336Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.59 seconds.
    337Chain 3: Adjust your expectations accordingly!
    338Chain 3:
    339Chain 3:
    340Chain 4:
    341Chain 4: Gradient evaluation took 4.4e-05 seconds
    342Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds.
    343Chain 4: Adjust your expectations accordingly!
    344Chain 4:
    345Chain 4:
    346Chain 1:
    347Chain 1:  Elapsed Time: 0.387624 seconds (Warm-up)
    348Chain 1:                0.033915 seconds (Sampling)
    349Chain 1:                0.421539 seconds (Total)
    350Chain 1:
    351Chain 2:
    352Chain 2:  Elapsed Time: 0.385989 seconds (Warm-up)
    353Chain 2:                0.044089 seconds (Sampling)
    354Chain 2:                0.430078 seconds (Total)
    355Chain 2:
    356Chain 3:
    357Chain 3:  Elapsed Time: 0.409957 seconds (Warm-up)
    358Chain 3:                0.039642 seconds (Sampling)
    359Chain 3:                0.449599 seconds (Total)
    360Chain 3:
    361Chain 4:
    362Chain 4:  Elapsed Time: 0.363549 seconds (Warm-up)
    363Chain 4:                0.042624 seconds (Sampling)
    364Chain 4:                0.406173 seconds (Total)
    365Chain 4:
    366Chain 1:
    367Chain 1: Gradient evaluation took 7.8e-05 seconds
    368Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.78 seconds.
    369Chain 1: Adjust your expectations accordingly!
    370Chain 1:
    371Chain 1:
    372Chain 2:
    373Chain 2: Gradient evaluation took 5.1e-05 seconds
    374Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
    375Chain 2: Adjust your expectations accordingly!
    376Chain 2:
    377Chain 2:
    378Chain 3:
    379Chain 3: Gradient evaluation took 4.7e-05 seconds
    380Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.47 seconds.
    381Chain 3: Adjust your expectations accordingly!
    382Chain 3:
    383Chain 3:
    384Chain 4:
    385Chain 4: Gradient evaluation took 5.1e-05 seconds
    386Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
    387Chain 4: Adjust your expectations accordingly!
    388Chain 4:
    389Chain 4:
    390Chain 2:
    391Chain 2:  Elapsed Time: 0.355205 seconds (Warm-up)
    392Chain 2:                0.051537 seconds (Sampling)
    393Chain 2:                0.406742 seconds (Total)
    394Chain 2:
    395Chain 1:
    396Chain 1:  Elapsed Time: 0.394264 seconds (Warm-up)
    397Chain 1:                0.037432 seconds (Sampling)
    398Chain 1:                0.431696 seconds (Total)
    399Chain 1:
    400Chain 3:
    401Chain 3:  Elapsed Time: 0.383287 seconds (Warm-up)
    402Chain 3:                0.036322 seconds (Sampling)
    403Chain 3:                0.419609 seconds (Total)
    404Chain 3:
    405Chain 4:
    406Chain 4:  Elapsed Time: 0.335487 seconds (Warm-up)
    407Chain 4:                0.048271 seconds (Sampling)
    408Chain 4:                0.383758 seconds (Total)
    409Chain 4:
    410Chain 1:
    411Chain 1: Gradient evaluation took 5e-05 seconds
    412Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.5 seconds.
    413Chain 1: Adjust your expectations accordingly!
    414Chain 1:
    415Chain 1:
    416Chain 2:
    417Chain 2: Gradient evaluation took 5.5e-05 seconds
    418Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.55 seconds.
    419Chain 2: Adjust your expectations accordingly!
    420Chain 2:
    421Chain 2:
    422Chain 3:
    423Chain 3: Gradient evaluation took 6.3e-05 seconds
    424Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.63 seconds.
    425Chain 3: Adjust your expectations accordingly!
    426Chain 3:
    427Chain 3:
    428Chain 4:
    429Chain 4: Gradient evaluation took 5.1e-05 seconds
    430Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
    431Chain 4: Adjust your expectations accordingly!
    432Chain 4:
    433Chain 4:
    434Chain 1:
    435Chain 1:  Elapsed Time: 0.454863 seconds (Warm-up)
    436Chain 1:                0.037816 seconds (Sampling)
    437Chain 1:                0.492679 seconds (Total)
    438Chain 1:
    439Chain 2:
    440Chain 2:  Elapsed Time: 0.415499 seconds (Warm-up)
    441Chain 2:                0.037093 seconds (Sampling)
    442Chain 2:                0.452592 seconds (Total)
    443Chain 2:
    444Chain 3:
    445Chain 3:  Elapsed Time: 0.363295 seconds (Warm-up)
    446Chain 3:                0.061002 seconds (Sampling)
    447Chain 3:                0.424297 seconds (Total)
    448Chain 3:
    449Chain 4:
    450Chain 4:  Elapsed Time: 0.462541 seconds (Warm-up)
    451Chain 4:                0.043532 seconds (Sampling)
    452Chain 4:                0.506073 seconds (Total)
    453Chain 4:
    454Warning messages:
    4551: There were 470 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See
    456http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
    4572: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
    458http://mc-stan.org/misc/warnings.html#bfmi-low
    4593: Examine the pairs() plot to diagnose sampling problems
    460
    4614: The largest R-hat is 1.05, indicating chains have not mixed.
    462Running the chains for more iterations may help. See
    463http://mc-stan.org/misc/warnings.html#r-hat
    4645: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
    465Running the chains for more iterations may help. See
    466http://mc-stan.org/misc/warnings.html#bulk-ess
    4676: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
    468Running the chains for more iterations may help. See
    469http://mc-stan.org/misc/warnings.html#tail-ess
    
    1nSampleEff %>% tibble(nWarmup=warmTrial)
    
     1# A tibble: 10 x 2
     2   .[,"a[1]"] [,"b[1]"] [,"sigma"] nWarmup
     3        <dbl>     <dbl>      <dbl>   <int>
     4 1      1051.      399.       45.8      10
     5 2      3101.     3085.     3135.     1120
     6 3      3515.     3529.     3214.     2230
     7 4      3122.     3277.     3522.     3340
     8 5      3145.     3382.     3322.     4450
     9 6      3378.     3193.     3701.     5560
    10 7      3299.     3539.     3149.     6670
    11 8      3570.     3050.     3079.     7780
    12 9      3247.     3148.     3340.     8890
    1310      3159.     2929.     2960.    10000
    
  • Exponential Model

    1for(i in 1:length(warmTrial)){
    2  tmp<-ulam(m91exp,chains=4,cores=4,refresh=-1,warmup=warmTrial[i],iter=1000+warmTrial[i])
    3  nSampleEffExp[i,]<-precis(tmp,2)$n_eff
    4}
    
      1
      2Chain 1:
      3Chain 1: Gradient evaluation took 3e-05 seconds
      4Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
      5Chain 1: Adjust your expectations accordingly!
      6Chain 1:
      7Chain 1:
      8Chain 1: WARNING: No variance estimation is
      9Chain 1:          performed for num_warmup < 20
     10Chain 1:
     11Chain 2:
     12Chain 2: Gradient evaluation took 2.8e-05 seconds
     13Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
     14Chain 2: Adjust your expectations accordingly!
     15Chain 2:
     16Chain 2:
     17Chain 2: WARNING: No variance estimation is
     18Chain 2:          performed for num_warmup < 20
     19Chain 2:
     20Chain 3:
     21Chain 3: Gradient evaluation took 2.7e-05 seconds
     22Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
     23Chain 3: Adjust your expectations accordingly!
     24Chain 3:
     25Chain 3:
     26Chain 3: WARNING: No variance estimation is
     27Chain 3:          performed for num_warmup < 20
     28Chain 3:
     29Chain 4:
     30Chain 4: Gradient evaluation took 2.5e-05 seconds
     31Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
     32Chain 4: Adjust your expectations accordingly!
     33Chain 4:
     34Chain 4:
     35Chain 4: WARNING: No variance estimation is
     36Chain 4:          performed for num_warmup < 20
     37Chain 4:
     38Chain 1:
     39Chain 1:  Elapsed Time: 0.000825 seconds (Warm-up)
     40Chain 1:                0.074412 seconds (Sampling)
     41Chain 1:                0.075237 seconds (Total)
     42Chain 1:
     43Chain 3:
     44Chain 3:  Elapsed Time: 0.000433 seconds (Warm-up)
     45Chain 3:                0.05124 seconds (Sampling)
     46Chain 3:                0.051673 seconds (Total)
     47Chain 3:
     48Chain 4:
     49Chain 4:  Elapsed Time: 0.001827 seconds (Warm-up)
     50Chain 4:                0.059028 seconds (Sampling)
     51Chain 4:                0.060855 seconds (Total)
     52Chain 4:
     53Chain 2:
     54Chain 2:  Elapsed Time: 0.005092 seconds (Warm-up)
     55Chain 2:                0.111894 seconds (Sampling)
     56Chain 2:                0.116986 seconds (Total)
     57Chain 2:
     58Chain 1:
     59Chain 1: Gradient evaluation took 2.7e-05 seconds
     60Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
     61Chain 1: Adjust your expectations accordingly!
     62Chain 1:
     63Chain 1:
     64Chain 2:
     65Chain 2: Gradient evaluation took 2.7e-05 seconds
     66Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
     67Chain 2: Adjust your expectations accordingly!
     68Chain 2:
     69Chain 2:
     70Chain 3:
     71Chain 3: Gradient evaluation took 2.7e-05 seconds
     72Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
     73Chain 3: Adjust your expectations accordingly!
     74Chain 3:
     75Chain 3:
     76Chain 4:
     77Chain 4: Gradient evaluation took 3.3e-05 seconds
     78Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
     79Chain 4: Adjust your expectations accordingly!
     80Chain 4:
     81Chain 4:
     82Chain 1:
     83Chain 1:  Elapsed Time: 0.048625 seconds (Warm-up)
     84Chain 1:                0.03552 seconds (Sampling)
     85Chain 1:                0.084145 seconds (Total)
     86Chain 1:
     87Chain 2:
     88Chain 2:  Elapsed Time: 0.045476 seconds (Warm-up)
     89Chain 2:                0.03068 seconds (Sampling)
     90Chain 2:                0.076156 seconds (Total)
     91Chain 2:
     92Chain 3:
     93Chain 3:  Elapsed Time: 0.059723 seconds (Warm-up)
     94Chain 3:                0.034697 seconds (Sampling)
     95Chain 3:                0.09442 seconds (Total)
     96Chain 3:
     97Chain 4:
     98Chain 4:  Elapsed Time: 0.055658 seconds (Warm-up)
     99Chain 4:                0.030508 seconds (Sampling)
    100Chain 4:                0.086166 seconds (Total)
    101Chain 4:
    102Chain 1:
    103Chain 1: Gradient evaluation took 3.7e-05 seconds
    104Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
    105Chain 1: Adjust your expectations accordingly!
    106Chain 1:
    107Chain 1:
    108Chain 2:
    109Chain 2: Gradient evaluation took 2.8e-05 seconds
    110Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    111Chain 2: Adjust your expectations accordingly!
    112Chain 2:
    113Chain 2:
    114Chain 3:
    115Chain 3: Gradient evaluation took 2.6e-05 seconds
    116Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
    117Chain 3: Adjust your expectations accordingly!
    118Chain 3:
    119Chain 3:
    120Chain 4:
    121Chain 4: Gradient evaluation took 2.6e-05 seconds
    122Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
    123Chain 4: Adjust your expectations accordingly!
    124Chain 4:
    125Chain 4:
    126Chain 1:
    127Chain 1:  Elapsed Time: 0.10462 seconds (Warm-up)
    128Chain 1:                0.035886 seconds (Sampling)
    129Chain 1:                0.140506 seconds (Total)
    130Chain 1:
    131Chain 3:
    132Chain 3:  Elapsed Time: 0.08377 seconds (Warm-up)
    133Chain 3:                0.036938 seconds (Sampling)
    134Chain 3:                0.120708 seconds (Total)
    135Chain 3:
    136Chain 2:
    137Chain 2:  Elapsed Time: 0.096442 seconds (Warm-up)
    138Chain 2:                0.0419 seconds (Sampling)
    139Chain 2:                0.138342 seconds (Total)
    140Chain 2:
    141Chain 4:
    142Chain 4:  Elapsed Time: 0.091513 seconds (Warm-up)
    143Chain 4:                0.062068 seconds (Sampling)
    144Chain 4:                0.153581 seconds (Total)
    145Chain 4:
    146Chain 1:
    147Chain 1: Gradient evaluation took 3.1e-05 seconds
    148Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    149Chain 1: Adjust your expectations accordingly!
    150Chain 1:
    151Chain 1:
    152Chain 2:
    153Chain 2: Gradient evaluation took 3.1e-05 seconds
    154Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    155Chain 2: Adjust your expectations accordingly!
    156Chain 2:
    157Chain 2:
    158Chain 3:
    159Chain 3: Gradient evaluation took 2.6e-05 seconds
    160Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
    161Chain 3: Adjust your expectations accordingly!
    162Chain 3:
    163Chain 3:
    164Chain 4:
    165Chain 4: Gradient evaluation took 2.4e-05 seconds
    166Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
    167Chain 4: Adjust your expectations accordingly!
    168Chain 4:
    169Chain 4:
    170Chain 1:
    171Chain 1:  Elapsed Time: 0.12662 seconds (Warm-up)
    172Chain 1:                0.033816 seconds (Sampling)
    173Chain 1:                0.160436 seconds (Total)
    174Chain 1:
    175Chain 2:
    176Chain 2:  Elapsed Time: 0.119929 seconds (Warm-up)
    177Chain 2:                0.048545 seconds (Sampling)
    178Chain 2:                0.168474 seconds (Total)
    179Chain 2:
    180Chain 4:
    181Chain 4:  Elapsed Time: 0.118087 seconds (Warm-up)
    182Chain 4:                0.040127 seconds (Sampling)
    183Chain 4:                0.158214 seconds (Total)
    184Chain 4:
    185Chain 3:
    186Chain 3:  Elapsed Time: 0.131859 seconds (Warm-up)
    187Chain 3:                0.053882 seconds (Sampling)
    188Chain 3:                0.185741 seconds (Total)
    189Chain 3:
    190Chain 1:
    191Chain 1: Gradient evaluation took 2.8e-05 seconds
    192Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    193Chain 1: Adjust your expectations accordingly!
    194Chain 1:
    195Chain 1:
    196Chain 2:
    197Chain 2: Gradient evaluation took 3.1e-05 seconds
    198Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    199Chain 2: Adjust your expectations accordingly!
    200Chain 2:
    201Chain 2:
    202Chain 3:
    203Chain 3: Gradient evaluation took 3.8e-05 seconds
    204Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
    205Chain 3: Adjust your expectations accordingly!
    206Chain 3:
    207Chain 3:
    208Chain 4:
    209Chain 4: Gradient evaluation took 2.8e-05 seconds
    210Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    211Chain 4: Adjust your expectations accordingly!
    212Chain 4:
    213Chain 4:
    214Chain 1:
    215Chain 1:  Elapsed Time: 0.178377 seconds (Warm-up)
    216Chain 1:                0.051505 seconds (Sampling)
    217Chain 1:                0.229882 seconds (Total)
    218Chain 1:
    219Chain 2:
    220Chain 2:  Elapsed Time: 0.19769 seconds (Warm-up)
    221Chain 2:                0.028853 seconds (Sampling)
    222Chain 2:                0.226543 seconds (Total)
    223Chain 2:
    224Chain 3:
    225Chain 3:  Elapsed Time: 0.209976 seconds (Warm-up)
    226Chain 3:                0.049889 seconds (Sampling)
    227Chain 3:                0.259865 seconds (Total)
    228Chain 3:
    229Chain 4:
    230Chain 4:  Elapsed Time: 0.233406 seconds (Warm-up)
    231Chain 4:                0.03065 seconds (Sampling)
    232Chain 4:                0.264056 seconds (Total)
    233Chain 4:
    234Chain 1:
    235Chain 1: Gradient evaluation took 3.1e-05 seconds
    236Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    237Chain 1: Adjust your expectations accordingly!
    238Chain 1:
    239Chain 1:
    240Chain 2:
    241Chain 2: Gradient evaluation took 3e-05 seconds
    242Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
    243Chain 2: Adjust your expectations accordingly!
    244Chain 2:
    245Chain 2:
    246Chain 3:
    247Chain 3: Gradient evaluation took 2.9e-05 seconds
    248Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.29 seconds.
    249Chain 3: Adjust your expectations accordingly!
    250Chain 3:
    251Chain 3:
    252Chain 4:
    253Chain 4: Gradient evaluation took 3.3e-05 seconds
    254Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
    255Chain 4: Adjust your expectations accordingly!
    256Chain 4:
    257Chain 4:
    258Chain 1:
    259Chain 1:  Elapsed Time: 0.226122 seconds (Warm-up)
    260Chain 1:                0.031202 seconds (Sampling)
    261Chain 1:                0.257324 seconds (Total)
    262Chain 1:
    263Chain 2:
    264Chain 2:  Elapsed Time: 0.229103 seconds (Warm-up)
    265Chain 2:                0.028798 seconds (Sampling)
    266Chain 2:                0.257901 seconds (Total)
    267Chain 2:
    268Chain 3:
    269Chain 3:  Elapsed Time: 0.238441 seconds (Warm-up)
    270Chain 3:                0.03459 seconds (Sampling)
    271Chain 3:                0.273031 seconds (Total)
    272Chain 3:
    273Chain 4:
    274Chain 4:  Elapsed Time: 0.228212 seconds (Warm-up)
    275Chain 4:                0.036574 seconds (Sampling)
    276Chain 4:                0.264786 seconds (Total)
    277Chain 4:
    278Chain 1:
    279Chain 1: Gradient evaluation took 3e-05 seconds
    280Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
    281Chain 1: Adjust your expectations accordingly!
    282Chain 1:
    283Chain 1:
    284Chain 2:
    285Chain 2: Gradient evaluation took 2.7e-05 seconds
    286Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
    287Chain 2: Adjust your expectations accordingly!
    288Chain 2:
    289Chain 2:
    290Chain 3:
    291Chain 3: Gradient evaluation took 2.8e-05 seconds
    292Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    293Chain 3: Adjust your expectations accordingly!
    294Chain 3:
    295Chain 3:
    296Chain 4:
    297Chain 4: Gradient evaluation took 2.4e-05 seconds
    298Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
    299Chain 4: Adjust your expectations accordingly!
    300Chain 4:
    301Chain 4:
    302Chain 1:
    303Chain 1:  Elapsed Time: 0.285094 seconds (Warm-up)
    304Chain 1:                0.048246 seconds (Sampling)
    305Chain 1:                0.33334 seconds (Total)
    306Chain 1:
    307Chain 4:
    308Chain 4:  Elapsed Time: 0.263814 seconds (Warm-up)
    309Chain 4:                0.034411 seconds (Sampling)
    310Chain 4:                0.298225 seconds (Total)
    311Chain 4:
    312Chain 2:
    313Chain 2:  Elapsed Time: 0.305613 seconds (Warm-up)
    314Chain 2:                0.039627 seconds (Sampling)
    315Chain 2:                0.34524 seconds (Total)
    316Chain 2:
    317Chain 3:
    318Chain 3:  Elapsed Time: 0.314974 seconds (Warm-up)
    319Chain 3:                0.040995 seconds (Sampling)
    320Chain 3:                0.355969 seconds (Total)
    321Chain 3:
    322Chain 1:
    323Chain 1: Gradient evaluation took 3e-05 seconds
    324Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
    325Chain 1: Adjust your expectations accordingly!
    326Chain 1:
    327Chain 1:
    328Chain 2:
    329Chain 2: Gradient evaluation took 2.8e-05 seconds
    330Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.28 seconds.
    331Chain 2: Adjust your expectations accordingly!
    332Chain 2:
    333Chain 2:
    334Chain 3:
    335Chain 3: Gradient evaluation took 3e-05 seconds
    336Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
    337Chain 3: Adjust your expectations accordingly!
    338Chain 3:
    339Chain 3:
    340Chain 4:
    341Chain 4: Gradient evaluation took 2.7e-05 seconds
    342Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
    343Chain 4: Adjust your expectations accordingly!
    344Chain 4:
    345Chain 4:
    346Chain 1:
    347Chain 1:  Elapsed Time: 0.293873 seconds (Warm-up)
    348Chain 1:                0.028447 seconds (Sampling)
    349Chain 1:                0.32232 seconds (Total)
    350Chain 1:
    351Chain 3:
    352Chain 3:  Elapsed Time: 0.288431 seconds (Warm-up)
    353Chain 3:                0.04518 seconds (Sampling)
    354Chain 3:                0.333611 seconds (Total)
    355Chain 3:
    356Chain 2:
    357Chain 2:  Elapsed Time: 0.36752 seconds (Warm-up)
    358Chain 2:                0.027453 seconds (Sampling)
    359Chain 2:                0.394973 seconds (Total)
    360Chain 2:
    361Chain 4:
    362Chain 4:  Elapsed Time: 0.383786 seconds (Warm-up)
    363Chain 4:                0.032154 seconds (Sampling)
    364Chain 4:                0.41594 seconds (Total)
    365Chain 4:
    366Chain 1:
    367Chain 1: Gradient evaluation took 3.1e-05 seconds
    368Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    369Chain 1: Adjust your expectations accordingly!
    370Chain 1:
    371Chain 1:
    372Chain 2:
    373Chain 2: Gradient evaluation took 2.7e-05 seconds
    374Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
    375Chain 2: Adjust your expectations accordingly!
    376Chain 2:
    377Chain 2:
    378Chain 3:
    379Chain 3: Gradient evaluation took 2.7e-05 seconds
    380Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
    381Chain 3: Adjust your expectations accordingly!
    382Chain 3:
    383Chain 3:
    384Chain 4:
    385Chain 4: Gradient evaluation took 3.1e-05 seconds
    386Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    387Chain 4: Adjust your expectations accordingly!
    388Chain 4:
    389Chain 4:
    390Chain 1:
    391Chain 1:  Elapsed Time: 0.354637 seconds (Warm-up)
    392Chain 1:                0.027996 seconds (Sampling)
    393Chain 1:                0.382633 seconds (Total)
    394Chain 1:
    395Chain 2:
    396Chain 2:  Elapsed Time: 0.339298 seconds (Warm-up)
    397Chain 2:                0.030313 seconds (Sampling)
    398Chain 2:                0.369611 seconds (Total)
    399Chain 2:
    400Chain 4:
    401Chain 4:  Elapsed Time: 0.319312 seconds (Warm-up)
    402Chain 4:                0.027904 seconds (Sampling)
    403Chain 4:                0.347216 seconds (Total)
    404Chain 4:
    405Chain 3:
    406Chain 3:  Elapsed Time: 0.333814 seconds (Warm-up)
    407Chain 3:                0.032747 seconds (Sampling)
    408Chain 3:                0.366561 seconds (Total)
    409Chain 3:
    410Chain 1:
    411Chain 1: Gradient evaluation took 3.8e-05 seconds
    412Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
    413Chain 1: Adjust your expectations accordingly!
    414Chain 1:
    415Chain 1:
    416Chain 2:
    417Chain 2: Gradient evaluation took 3.7e-05 seconds
    418Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.37 seconds.
    419Chain 2: Adjust your expectations accordingly!
    420Chain 2:
    421Chain 2:
    422Chain 3:
    423Chain 3: Gradient evaluation took 2.6e-05 seconds
    424Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
    425Chain 3: Adjust your expectations accordingly!
    426Chain 3:
    427Chain 3:
    428Chain 4:
    429Chain 4: Gradient evaluation took 3.1e-05 seconds
    430Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
    431Chain 4: Adjust your expectations accordingly!
    432Chain 4:
    433Chain 4:
    434Chain 1:
    435Chain 1:  Elapsed Time: 0.343002 seconds (Warm-up)
    436Chain 1:                0.028223 seconds (Sampling)
    437Chain 1:                0.371225 seconds (Total)
    438Chain 1:
    439Chain 3:
    440Chain 3:  Elapsed Time: 0.305501 seconds (Warm-up)
    441Chain 3:                0.027824 seconds (Sampling)
    442Chain 3:                0.333325 seconds (Total)
    443Chain 3:
    444Chain 2:
    445Chain 2:  Elapsed Time: 0.358116 seconds (Warm-up)
    446Chain 2:                0.02975 seconds (Sampling)
    447Chain 2:                0.387866 seconds (Total)
    448Chain 2:
    449Chain 4:
    450Chain 4:  Elapsed Time: 0.330673 seconds (Warm-up)
    451Chain 4:                0.030407 seconds (Sampling)
    452Chain 4:                0.36108 seconds (Total)
    453Chain 4:
    454Warning messages:
    4551: There were 14 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See
    456http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
    4572: Examine the pairs() plot to diagnose sampling problems
    
    1nSampleEffExp %>% tibble(nWarmup=warmTrial)
    
     1# A tibble: 10 x 2
     2   .[,1]  [,2]  [,3] nWarmup
     3   <dbl> <dbl> <dbl>   <int>
     4 1 3063. 1265. 1045.      10
     5 2 3238. 2909. 3425.    1120
     6 3 3268. 2949. 2966.    2230
     7 4 3123. 3257. 3218.    3340
     8 5 3227. 3449. 3554.    4450
     9 6 3494. 3743. 3343.    5560
    10 7 3200. 3012. 2879.    6670
    11 8 3210. 3207. 3125.    7780
    12 9 2919. 3329. 3060.    8890
    1310 2951. 3461. 3078.   10000
    

    It is important to note that the divergent transitions are probably why the number of effective samples decrease in the last two rows.

  • Results

    We can see that the number of effective samples increases almost constantly. This is probably due to correlations in the chain, which are removed during the warmup period.

Chapter XI: God Spiked The Integers

Easy Questions (Ch11)

11E1

If an event has probability \(0.35\), what are the log-odds of this event?

Solution
1log(0.35/(1-0.35))
1[1] -0.6190392

11E2

If an event has log-odds \(3.2\), what is the probability of this event?

Solution
1logistic(3.2)
1[1] 0.9608343

11E3

Suppose that a coefficient in a logistic regression has value \(1.7\). What does this imply about the proportional change in odds of the outcome?

Solution
1exp(1.7)
1[1] 5.473947

Note that this is not really the change in the variable, but the proportional odds.

HOLD 11E4

Why do Poisson regressions sometimes require the use of an offset? Provide an example.

Solution

The Poisson distribution is often understood as a limiting distribution of the Binomial where \(λ=np\) as \(n→∞\) and \(p→0\). The single parameter thus expresses the expected value, but is often used to encode different time-steps as well. Essentially, the distribution assumes a constant rate in time or space, and thus the change in exposure, is expressed by the offset, which is the logarithm of the exposure.

For any case where samples are drawn from populations which have different aggregation time periods but are still within the purview of a Poisson distribution, the offset is a natural way of expressing these.

To leverage the example of the book, when constructing a model to account for the fact that one Monastery calculates their averages on a weekly basis, while the other averages by day, this constraint should be modeled by having differing offsets.

Questions of Medium Complexity (Ch11)

HOLD 11M2

If a coefficient in a Poisson regression has values \(1.7\), what does this imply about the change in the outcome?

Solution
1exp(1.7)
1[1] 5.473947

The coefficient in a Poisson regression implies that the proportional change in the count will be ~5.474 when the predictor variable increase by one unit.

HOLD 11M3

Explain why the logit link is appropriate for a binomial generalized linear model.

Solution

The logit link essentially connects a parameter constrained between zero and one and the real space. The logit function is defined as:

\[\mathrm{logit}(pᵢ)=\log{\frac{pᵢ}{1-pᵢ}}\]

Where \(pᵢ\) is a probability mass. This link makes sense for a GLM since the predicted value is a probability distribution parameter, and we would like to obtain this from a linear model which spans the entire set of real numbers.

1curve(logit,from=-0.5,to=1.5)

The link function maps a parameter onto a linear model.

HOLD 11M4

Explain why the log is appropriate for a Poisson generalized linear model.

Solution

The log function ensures that the parameter cannot take values which are less than zero. This is a natural consequence of the function definition.

1curve(log,from=-0.5,to=100000)

The log link assumes that the parameter value is the exponentiation of the linear model.

This makes sense for a Poisson GLM as the Poisson distribution does not accept negative values.

HOLD 11M5

What would it imply to use a logit link for the mean of a Poisson generalized linear model? Can you think of a real research problem for which this would make sense?

Solution

We should write this out more explicitly.

This implies that the mean μ lies between zero and one. Since the Poisson distribution is defined by a single parameter, this does limit the model outputs. The premise of a Poisson regression problem is that the GLM models a count with an unknown maximum, so it does seem to be a very severe restriction.

To my mind this is feasible for constrained problems, where the Poisson distribution is to be followed but only within a particular range for some reason, and when the Binomial (of which the Poisson is a special case), decreases too slowly.

It was mentioned on the class forums, that the COVID-19 problem was modeled with a two parameter generalized link function, i.e. \(\log{\frac{p}{S-p}}\) which essentially constrains the model to have Poisson dynamics but with an output mean between 0 and S.

HOLD 11M6

State the constraints for which the binomial and Poisson distributions have maximum entropy. Are the constraints different at all for binomial and Poisson? Why or why not?

Solution

The Binomial distribution is defined to be the maximum entropy distributions are:

  1. Discrete binary outcomes
  2. Constant probability (or expectation)

This is defined by the number of outcomes (n) as well as the probability (p). The experiment is essentially reduced to a series of independent and identical Bernoulli trials with only two outcomes. The Poisson distribution is derived as a limiting form of the Binomial, where \(n→∞\) and \(p→0\). Since this does not change the underlying constraints, this is still a maximum entropy distribution.

HOLD 11M8

Revisit the data(Kline) islands example. This time drop Hawaii from the sample and refit the models. What changes do you observe?

Solution
1data(Kline)
2kDat<-Kline
3kDat<-kDat %>% dplyr::mutate(cid=ifelse(contact=="high",2,1),
4                             stdPop=standardize(log(population))) %>% filter(culture!="Hawaii")
5datList<-list(
6  totTools=kDat$total_tools,
7  stdPop=kDat$stdPop,
8  cid=as.integer(kDat$cid)
9)

We can now fit this.

1m11m10res<-ulam(
2  alist(
3    totTools ~ dpois(lambda),
4    log(lambda)<-a[cid]+b[cid]*stdPop,
5    a[cid] ~ dnorm(3,0.5),
6    b[cid] ~ dnorm(0,0.2)
7  ),data=datList, chains=4, cores=4
8)
1m11m10res %>% precis(2)
1     mean   sd  5.5% 94.5% n_eff Rhat4
2a[1] 3.18 0.12  2.99  3.37  1621     1
3a[2] 3.61 0.08  3.48  3.73  1962     1
4b[1] 0.19 0.13 -0.01  0.39  1639     1
5b[2] 0.19 0.16 -0.06  0.44  1830     1

We see that the slopes are now the same, which makes sense since in this data-set Hawaii was the only outlier.

Chapter XII: Monsters and Mixtures

Easy Questions (Ch12)

HOLD 12E4

Over-dispersion is common in count data. Give an example of a natural process that might produce over-dispersed counts. Can you also give an example of a process that might produce /under-/dispersed counts?

Solution
  • Over-dispersion

    Over dispersion is essentially the occurrence of greater variability than accounted for based on the statistical model. The presence of over-dispersion is typically due to heterogeneity in populations. This heterogeneity may arise from simple aggregation issues like in the case considered in the text, of Monasteries which accumulate data weekly or daily, in-spite of following the same Poisson model.

  • Under-dispersion

    Under dispersion is essentially the occurrence of less variability than accounted for based on the statistical model. The clearest example of under-dispersion is from the draws of an MCMC sampler. The number of effective samples is typically lower than the number of samples, as the data is highly correlated (autocorrelated) as the sampler draws sequential samples. For a count model, if a hidden rate limiting variable exists and has not been accounted for, then the variation in counts is lowered, and will show up as under-dispersion.

Hard Questions (Ch12)

12H1

In 2014, a paper was published that was entitle “Female hurricanes are deadlier than male hurricanes.” As the title suggests, the paper claimed that hurricanes with female names have caused greater loss of life, and the explanation given is that people unconsciously rate female hurricanes as less dangerous and so are less likely to evacuate. Statisticians severely criticized the paper after publication. Here, you’ll explore the complete data used in the paper and consider the hypothesis that hurricanes with female names are deadlier. Load the data with:

1library(rethinking)
2data(Hurricanes)

Acquaint yourself with the columns by inspecting the help ?Hurricanes. In this problem, you’ll focus on predicting deaths using feminity as a predictor. You can use quap or ulam. Compare the model to an intercept-only Poisson model of deaths. How strong is the association between feminity of name and deaths? Which storms does the model fit (retrodict) well? Which storms does it fit poorly?

Solution

Since I have no understanding of hurricanes except that it is unlikely to have too much of an effect. I will run through some sample priors. Presumably, most hurricanes do not kill over a thousand people. Furthermore, a-priori, I would not like to assume that femininity is positive or negative, so I will instead encode a belief that it shouldn’t matter much either way, ergo a Gaussian.

1N<-100
2a<-rnorm(N,1,0.5)
3bF<-rnorm(N,0.5,2)
4seqF<-seq(from=-2,to=2,length.out=100)
5plot(NULL,xlim=c(-2,2),ylim=c(0,1000),xlab="Femininity",ylab="deaths")
6for(i in 1:N) lines(seqF,exp(a[i]+bF[i]*seqF),col=grau())

This seems to be reasonable to me. It does have a bit of an unreasonable focus on 0, but it does also seem to mostly hug the x-axis in a way indicating my prior belief that it should not matter all that much. There is enough diversity in the priors to allow for stronger trends, but they are by and large unlikely.

Now we can actually use these in a model.

1data(Hurricanes)
2hurDat<-Hurricanes %>% as.data.frame
3hurDat<-hurDat %>% dplyr::mutate(femStd=standardize(femininity))
4datListH<-list(deaths=hurDat$deaths,femStd=hurDat$femStd)
5m12h1norm<-ulam(alist(deaths ~ dpois(lambda),
6                 log(lambda) <- a+bF*femStd,
7                 a ~ dnorm(1,0.5),
8                 bF ~ dnorm(0.5,2)
9                 ),data=datListH, chains=4, cores=4,log_lik = TRUE)
 1SAMPLING FOR MODEL 'bd16fb771b491de48a3f8ce09fc68301' NOW (CHAIN 1).
 2
 3SAMPLING FOR MODEL 'bd16fb771b491de48a3f8ce09fc68301' NOW (CHAIN 2).
 4Chain 2:
 5Chain 1:
 6Chain Chain 12: : Gradient evaluation took 5.1e-05 secondsGradient evaluation took 3.9e-05 seconds
 7
 8Chain Chain 12: : 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
 9
10Chain Chain 12: : Adjust your expectations accordingly!Adjust your expectations accordingly!
11
12Chain Chain 12: :
13
14Chain Chain 12: :
15
16Chain Chain 21: Iteration:   1 / 1000 [  0%]  (Warmup)
17: Iteration:   1 / 1000 [  0%]  (Warmup)
18Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
19Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
20
21SAMPLING FOR MODEL 'bd16fb771b491de48a3f8ce09fc68301' NOW (CHAIN 3).
22Chain 3:
23Chain 3: Gradient evaluation took 4.3e-05 seconds
24Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds.
25Chain 3: Adjust your expectations accordingly!
26Chain 3:
27Chain 3:
28Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
29Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
30Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
31Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
32Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
33
34SAMPLING FOR MODEL 'bd16fb771b491de48a3f8ce09fc68301' NOW (CHAIN 4).
35Chain 4:
36Chain 4: Gradient evaluation took 3.1e-05 seconds
37Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
38Chain 4: Adjust your expectations accordingly!
39Chain 4:
40Chain 4:
41Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
42Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
43Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
44Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
45Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
46Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
47Chain 4: Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
48Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
49Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
50Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
51Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
52Iteration: 200 / 1000 [ 20%]  (Warmup)
53Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
54Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
55Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
56Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
57Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
58Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
59Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
60Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
61Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
62Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
63Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
64Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
65Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
66Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
67Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
68Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
69Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
70Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
71Chain 1:
72Chain 1:  Elapsed Time: 0.055115 seconds (Warm-up)
73Chain 1:                0.042775 seconds (Sampling)
74Chain 1:                0.09789 seconds (Total)
75Chain 1:
76Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
77Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
78Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
79Chain 2:
80Chain 2:  Elapsed Time: 0.050636 seconds (Warm-up)
81Chain 2:                0.049372 seconds (Sampling)
82Chain 2:                0.100008 seconds (Total)
83Chain 2:
84Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
85Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
86Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
87Chain 3:
88Chain 3:  Elapsed Time: 0.056256 seconds (Warm-up)
89Chain 3:                0.039069 seconds (Sampling)
90Chain 3:                0.095325 seconds (Total)
91Chain 3:
92Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
93Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
94Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
95Chain 4:
96Chain 4:  Elapsed Time: 0.054509 seconds (Warm-up)
97Chain 4:                0.036572 seconds (Sampling)
98Chain 4:                0.091081 seconds (Total)
99Chain 4:
1m12h1norm %>% precis
1   mean   sd 5.5% 94.5% n_eff Rhat4
2a  3.00 0.02 2.96  3.04  1334     1
3bF 0.24 0.02 0.20  0.28  1361     1

There are several points to be noted:

  • The number of effective samples is far lower than the number simulated (4000)
  • The bounds are quite tight for the values

All told, we can see that the model is quite certain of the values, but given the low number of effective samples it would make sense to run this model longer.

Furthermore, the model infers that from our data, there is a positive correlation (quite a high one) between femininity and deaths.

This is best seen by actually visualizing the results.

1m12h1norm %>% pairs

We should check the posterior as well.

 1k<-PSIS(m12h1norm,pointwise=TRUE)$k
 2plot(hurDat$femStd,hurDat$deaths,xlab="Standardized Femininity",ylab="Deaths",col=rangi2,pch=hurDat$female, lwd=2,cex=1+normalize(k))
 3## Axis for predictions
 4ns<-500
 5femininity<-seq(from=min(hurDat$femStd), to=max(hurDat$femStd),length.out = ns)
 6## Female
 7lambda<-link(m12h1norm,data=data.frame(femStd=femininity))
 8lmu<-apply(lambda,2,mean)
 9lci<-apply(lambda,2,PI)
10lines(femininity,lmu,lty=2,lwd=1.5)
11shade(lci,femininity,xpd=TRUE)

We can see that the model prediction does not actually handle the data very well, in that it is evident the model simply cannot account for the high death rate values. We have plotted the 89% interval as well (the default for PI).

We can also inspect the expect PSISk values.

1m12h1norm %>% PSISk %>% summary
1Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
2   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
3-0.1800 -0.0500  0.0600  0.1521  0.1500  2.5800

Clearly, there are some points with high leverage.

Finally we can plot the posterior distribution.

 1posterior<-m12h1norm %>% extract.samples
 2fem<-sample(hurDat$femStd,3)
 3for(i in 1:10){
 4  curve(dgamma(x,exp(posterior$a[i]+posterior$bF[i]*(fem[1]))), from=0,to=100, col="red",ylab="Density",xlab="Average deaths",add=ifelse(i==1,FALSE,TRUE))
 5  curve(dgamma(x,exp(posterior$a[i]+posterior$bF[i]*(fem[2]))), from=0,to=100, col="blue",add=TRUE)
 6  curve(dgamma(x,exp(posterior$a[i]+posterior$bF[i]*(fem[3]))), from=0,to=100, col="black",add=TRUE)}
 7legend("topright",
 8       legend=c(sprintf("Femininity=%.3f",fem[1]), sprintf("Femininity=%.3f",fem[2]), sprintf("Femininity=%.3f",fem[3])),
 9       col=c("red","blue","black"),
10       pch=19
11       )

As is to be expected, for each value of femininity, we have one family of gamma distributions.

HOLD 12H2

Counts are nearly always over-dispersed relative to Poisson. So fit a gamma-Poisson (aka negative-binomial) model to predict deaths using feminity. Show that the over-dispersed model no longer shows as precise a positive association between feminity and deaths, with an $89$% interval that overlaps zero. Can you explain why the association diminished in strength?

Solution

Recall that the gamma-Poisson has two parameters, one for the rate, and the other for the dispersion of rates. Larger values of the dispersion imply that the distribution is more similar to a pure Poisson process. For ensuring meaningful comparisons, we will keep the same priors as before. We will need a scale parameter, but we will postulate a simple exponential prior for that.

 1data(Hurricanes)
 2hurDat<-Hurricanes %>% as.data.frame
 3hurDat<-hurDat %>% dplyr::mutate(femStd=standardize(femininity))
 4datListH<-list(deaths=hurDat$deaths,femStd=hurDat$femStd)
 5m12h2norm<-ulam(alist(deaths ~ dgampois(lambda,scale),
 6                 log(lambda) <- a+bF*femStd,
 7                 a ~ dnorm(1,0.5),
 8                 bF ~ dnorm(0.5,2),
 9                 scale ~ dexp(1)
10                 ),data=datListH, chains=4, cores=4,log_lik = TRUE)
 1SAMPLING FOR MODEL '5dc94bd836781d34a208695cf643c56c' NOW (CHAIN 1).
 2Chain 1:
 3Chain 1: Gradient evaluation took 0.00012 seconds
 4Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
 5Chain 1: Adjust your expectations accordingly!
 6Chain 1:
 7Chain 1:
 8
 9SAMPLING FOR MODEL '5dc94bd836781d34a208695cf643c56c' NOW (CHAIN 2).
10Chain 2:
11Chain 2: Gradient evaluation took 9.6e-05 seconds
12Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.96 seconds.
13Chain 2: Adjust your expectations accordingly!
14Chain 2:
15Chain 2:
16Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
17Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
18
19SAMPLING FOR MODEL '5dc94bd836781d34a208695cf643c56c' NOW (CHAIN 3).
20Chain 3:
21Chain 3: Gradient evaluation took 9.1e-05 seconds
22Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.91 seconds.
23Chain 3: Adjust your expectations accordingly!
24Chain 3:
25Chain 3:
26Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
27
28SAMPLING FOR MODEL '5dc94bd836781d34a208695cf643c56c' NOW (CHAIN 4).
29Chain 4:
30Chain 4: Gradient evaluation took 8.2e-05 seconds
31Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.82 seconds.
32Chain 4: Adjust your expectations accordingly!
33Chain 4:
34Chain 4:
35Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
36Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
37Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
38Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
39Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
40Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
41Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
42Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
43Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
44Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
45Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
46Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
47Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
48Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
49Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
50Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
51Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
52Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
53Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
54Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
55Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
56Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
57Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
58Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
59Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
60Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
61Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
62Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
63Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
64Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
65Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
66Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
67Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
68Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
69Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
70Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
71Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
72Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
73Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
74Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
75Chain 2:
76Chain 2:  Elapsed Time: 0.119773 seconds (Warm-up)
77Chain 2:                0.091705 seconds (Sampling)
78Chain 2:                0.211478 seconds (Total)
79Chain 2:
80Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
81Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
82Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
83Chain 4:
84Chain 4:  Elapsed Time: 0.122953 seconds (Warm-up)
85Chain 4:                0.082244 seconds (Sampling)
86Chain 4:                0.205197 seconds (Total)
87Chain 4:
88Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
89Chain 3:
90Chain 3:  Elapsed Time: 0.115457 seconds (Warm-up)
91Chain 3:                0.114422 seconds (Sampling)
92Chain 3:                0.229879 seconds (Total)
93Chain 3:
94Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
95Chain 1:
96Chain 1:  Elapsed Time: 0.110931 seconds (Warm-up)
97Chain 1:                0.129622 seconds (Sampling)
98Chain 1:                0.240553 seconds (Total)
99Chain 1:
1m12h2norm %>% precis
2m12h1norm %>% precis
1      mean   sd  5.5% 94.5% n_eff Rhat4
2a     2.86 0.14  2.64  3.08  1705     1
3bF    0.21 0.14 -0.01  0.44  1946     1
4scale 0.45 0.06  0.35  0.55  1957     1
5
6   mean   sd 5.5% 94.5% n_eff Rhat4
7a  3.00 0.02 2.96  3.04  1334     1
8bF 0.24 0.02 0.20  0.28  1361     1

We note that the effective number of samples in the second model are greater, which implies that this model is less prone to correlations. We can quantify this with the WAIC as well.

1WAIC(m12h2norm) %>% rbind(WAIC(m12h1norm)) %>% tibble(model=c("Gamma-Poisson","Poisson")) %>% toOrg
1|            WAIC |              lppd |          penalty |          std_err | model         |
2|-----------------+-------------------+------------------+------------------+---------------|
3| 710.78471929582 | -351.457619486557 | 3.93474016135277 | 34.6128979470592 | Gamma-Poisson |
4| 4427.5667952452 | -2080.92835360918 | 132.855044013418 | 1009.13483879188 | Poisson       |

The WAIC values show that the gamma-Poisson model is less likely to over-fit.

We would like to see the models together.

1m12h1norm %>% precis(pars=c("a","bF")) %>% plot(col="blue")
2m12h2norm %>% precis(pars=c("a","bF")) %>% plot(add=TRUE,col="red")

We can see that there is little to no difference in the means, though the intervals seem wider than before. This is more clear in the coeftab plot.

1plot(coeftab(m12h1norm,m12h2norm))

An important consequence of this is that the model is no longer completely sure that there is any effect of femininity on the death count, as can be seen from the wider uncertainty interval, which includes 0.

We can also visualize the model with pairs.

1pairs(m12h2norm)
 1k<-PSIS(m12h2norm,pointwise=TRUE)$k
 2plot(hurDat$femStd,hurDat$deaths,xlab="Standardized Femininity",ylab="Deaths",col=rangi2,pch=hurDat$female, lwd=2,cex=1+normalize(k))
 3## Axis for predictions
 4ns<-500
 5femininity<-seq(from=min(hurDat$femStd),to=max(hurDat$femStd),length.out = ns)
 6## Gamma Poisson
 7lambda<-link(m12h2norm,data=data.frame(femStd=femininity))
 8lmu<-apply(lambda,2,mean)
 9lci<-apply(lambda,2,PI)
10shade(lci,femininity,xpd=TRUE,col="red")
11lines(femininity,lmu,lty=2,lwd=1.5,col="white")
12## Poisson
13lambda<-link(m12h1norm,data=data.frame(femStd=femininity))
14lmu<-apply(lambda,2,mean)
15lci<-apply(lambda,2,PI)
16shade(lci,femininity,xpd=TRUE,col="blue")
17lines(femininity,lmu,lty=2,lwd=1.5,col="white")

Clearly, the uncertainty of the newer model is much greater, even though the predictions do not differ much. Unfortunately, both models fail to account for storms with high death counts.

We would also like to plot the predicted distributions.

 1posterior<-m12h2norm %>% extract.samples
 2fem<-sample(hurDat$femStd,2)
 3for(i in 1:100){
 4  curve(dgamma2(x,exp(posterior$a[i]+posterior$bF[i]*(fem[1])),posterior$scale[i]), from=0,to=100,col="red",ylab="Density",xlab="Average deaths",add=ifelse(i==1,FALSE,TRUE))
 5  curve(dgamma2(x,exp(posterior$a[i]+posterior$bF[i]*(fem[2])),posterior$scale[i]), from=0,to=100,col="blue",add=TRUE)
 6  }
 7legend("topright",
 8       legend=c(sprintf("Femininity=%.3f",fem[1]), sprintf("Femininity=%.3f",fem[2])),
 9       col=c("red","blue"),
10       pch=19
11       )

This clearly has more spread than the previous predictions. By definition, the dispersion term tends to spread the distribution out, with higher values of the dispersion corresponding to a “true” Poisson distribution.

A: Colophon

To ensure that this document is fully reproducible at a later date, we will record the session info.

1devtools::session_info()
  1─ Session info ───────────────────────────────────────────────────────────────
  2 setting  value
  3 version  R version 4.0.0 (2020-04-24)
  4 os       Arch Linux
  5 system   x86_64, linux-gnu
  6 ui       X11
  7 language (EN)
  8 collate  en_US.UTF-8
  9 ctype    en_US.UTF-8
 10 tz       Iceland
 11 date     2020-06-21
 12
 13─ Packages ───────────────────────────────────────────────────────────────────
 14 package              * version    date       lib   source
 15 arrayhelpers           1.1-0      2020-02-04 [167] CRAN (R 4.0.0)
 16 assertthat             0.2.1      2019-03-21 [34]  CRAN (R 4.0.0)
 17 backports              1.1.6      2020-04-05 [68]  CRAN (R 4.0.0)
 18 boot                   1.3-24     2019-12-20 [5]   CRAN (R 4.0.0)
 19 broom                  0.5.6      2020-04-20 [67]  CRAN (R 4.0.0)
 20 callr                  3.4.3      2020-03-28 [87]  CRAN (R 4.0.0)
 21 cellranger             1.1.0      2016-07-27 [55]  CRAN (R 4.0.0)
 22 cli                    2.0.2      2020-02-28 [33]  CRAN (R 4.0.0)
 23 coda                   0.19-3     2019-07-05 [169] CRAN (R 4.0.0)
 24 colorspace             1.4-1      2019-03-18 [97]  CRAN (R 4.0.0)
 25 crayon                 1.3.4      2017-09-16 [35]  CRAN (R 4.0.0)
 26 curl                   4.3        2019-12-02 [26]  CRAN (R 4.0.0)
 27 dagitty              * 0.2-2      2016-08-26 [244] CRAN (R 4.0.0)
 28 data.table           * 1.12.8     2019-12-09 [27]  CRAN (R 4.0.0)
 29 DBI                    1.1.0      2019-12-15 [77]  CRAN (R 4.0.0)
 30 dbplyr                 1.4.3      2020-04-19 [76]  CRAN (R 4.0.0)
 31 desc                   1.2.0      2018-05-01 [84]  CRAN (R 4.0.0)
 32 devtools             * 2.3.0      2020-04-10 [219] CRAN (R 4.0.0)
 33 digest                 0.6.25     2020-02-23 [42]  CRAN (R 4.0.0)
 34 dplyr                * 0.8.5      2020-03-07 [69]  CRAN (R 4.0.0)
 35 ellipsis               0.3.0      2019-09-20 [30]  CRAN (R 4.0.0)
 36 evaluate               0.14       2019-05-28 [82]  CRAN (R 4.0.0)
 37 fansi                  0.4.1      2020-01-08 [36]  CRAN (R 4.0.0)
 38 forcats              * 0.5.0      2020-03-01 [29]  CRAN (R 4.0.0)
 39 fs                     1.4.1      2020-04-04 [109] CRAN (R 4.0.0)
 40 generics               0.0.2      2018-11-29 [71]  CRAN (R 4.0.0)
 41 ggplot2              * 3.3.0      2020-03-05 [78]  CRAN (R 4.0.0)
 42 glue                 * 1.4.0      2020-04-03 [37]  CRAN (R 4.0.0)
 43 gridExtra              2.3        2017-09-09 [123] CRAN (R 4.0.0)
 44 gtable                 0.3.0      2019-03-25 [79]  CRAN (R 4.0.0)
 45 haven                  2.2.0      2019-11-08 [28]  CRAN (R 4.0.0)
 46 hms                    0.5.3      2020-01-08 [44]  CRAN (R 4.0.0)
 47 htmltools              0.4.0      2019-10-04 [112] CRAN (R 4.0.0)
 48 httr                   1.4.1      2019-08-05 [100] CRAN (R 4.0.0)
 49 inline                 0.3.15     2018-05-18 [162] CRAN (R 4.0.0)
 50 jsonlite               1.6.1      2020-02-02 [101] CRAN (R 4.0.0)
 51 kableExtra           * 1.1.0      2019-03-16 [212] CRAN (R 4.0.0)
 52 knitr                  1.28       2020-02-06 [113] CRAN (R 4.0.0)
 53 latex2exp            * 0.4.0      2015-11-30 [211] CRAN (R 4.0.0)
 54 lattice                0.20-41    2020-04-02 [6]   CRAN (R 4.0.0)
 55 lifecycle              0.2.0      2020-03-06 [38]  CRAN (R 4.0.0)
 56 loo                    2.2.0      2019-12-19 [163] CRAN (R 4.0.0)
 57 lubridate              1.7.8      2020-04-06 [106] CRAN (R 4.0.0)
 58 magrittr               1.5        2014-11-22 [21]  CRAN (R 4.0.0)
 59 MASS                   7.3-51.5   2019-12-20 [7]   CRAN (R 4.0.0)
 60 matrixStats            0.56.0     2020-03-13 [164] CRAN (R 4.0.0)
 61 memoise                1.1.0      2017-04-21 [229] CRAN (R 4.0.0)
 62 modelr                 0.1.6      2020-02-22 [107] CRAN (R 4.0.0)
 63 munsell                0.5.0      2018-06-12 [96]  CRAN (R 4.0.0)
 64 mvtnorm                1.1-0      2020-02-24 [243] CRAN (R 4.0.0)
 65 nlme                   3.1-147    2020-04-13 [11]  CRAN (R 4.0.0)
 66 orgutils             * 0.4-1      2017-03-21 [209] CRAN (R 4.0.0)
 67 pillar                 1.4.3      2019-12-20 [39]  CRAN (R 4.0.0)
 68 pkgbuild               1.0.6      2019-10-09 [86]  CRAN (R 4.0.0)
 69 pkgconfig              2.0.3      2019-09-22 [43]  CRAN (R 4.0.0)
 70 pkgload                1.0.2      2018-10-29 [83]  CRAN (R 4.0.0)
 71 plyr                   1.8.6      2020-03-03 [73]  CRAN (R 4.0.0)
 72 prettyunits            1.1.1      2020-01-24 [58]  CRAN (R 4.0.0)
 73 printr               * 0.1        2017-05-19 [214] CRAN (R 4.0.0)
 74 processx               3.4.2      2020-02-09 [88]  CRAN (R 4.0.0)
 75 ps                     1.3.2      2020-02-13 [89]  CRAN (R 4.0.0)
 76 purrr                * 0.3.4      2020-04-17 [50]  CRAN (R 4.0.0)
 77 R6                     2.4.1      2019-11-12 [48]  CRAN (R 4.0.0)
 78 Rcpp                   1.0.4.6    2020-04-09 [10]  CRAN (R 4.0.0)
 79 readr                * 1.3.1      2018-12-21 [45]  CRAN (R 4.0.0)
 80 readxl                 1.3.1      2019-03-13 [54]  CRAN (R 4.0.0)
 81 remotes                2.1.1      2020-02-15 [233] CRAN (R 4.0.0)
 82 reprex                 0.3.0      2019-05-16 [108] CRAN (R 4.0.0)
 83 rethinking           * 2.01       2020-06-06 [242] local
 84 rlang                  0.4.5      2020-03-01 [31]  CRAN (R 4.0.0)
 85 rmarkdown              2.1        2020-01-20 [110] CRAN (R 4.0.0)
 86 rprojroot              1.3-2      2018-01-03 [85]  CRAN (R 4.0.0)
 87 rstan                * 2.19.3     2020-02-11 [161] CRAN (R 4.0.0)
 88 rstudioapi             0.11       2020-02-07 [91]  CRAN (R 4.0.0)
 89 rvest                  0.3.5      2019-11-08 [120] CRAN (R 4.0.0)
 90 scales                 1.1.0      2019-11-18 [93]  CRAN (R 4.0.0)
 91 sessioninfo            1.1.1      2018-11-05 [231] CRAN (R 4.0.0)
 92 shape                  1.4.4      2018-02-07 [193] CRAN (R 4.0.0)
 93 StanHeaders          * 2.19.2     2020-02-11 [165] CRAN (R 4.0.0)
 94 stringi                1.4.6      2020-02-17 [52]  CRAN (R 4.0.0)
 95 stringr              * 1.4.0      2019-02-10 [74]  CRAN (R 4.0.0)
 96 svUnit                 1.0.3      2020-04-20 [168] CRAN (R 4.0.0)
 97 testthat               2.3.2      2020-03-02 [81]  CRAN (R 4.0.0)
 98 textutils              0.2-0      2020-01-07 [210] CRAN (R 4.0.0)
 99 tibble               * 3.0.1      2020-04-20 [32]  CRAN (R 4.0.0)
100 tidybayes            * 2.0.3      2020-04-04 [166] CRAN (R 4.0.0)
101 tidybayes.rethinking * 2.0.3.9000 2020-06-07 [246] local
102 tidyr                * 1.0.2      2020-01-24 [75]  CRAN (R 4.0.0)
103 tidyselect             1.0.0      2020-01-27 [49]  CRAN (R 4.0.0)
104 tidyverse            * 1.3.0      2019-11-21 [66]  CRAN (R 4.0.0)
105 usethis              * 1.6.0      2020-04-09 [238] CRAN (R 4.0.0)
106 V8                     3.0.2      2020-03-14 [245] CRAN (R 4.0.0)
107 vctrs                  0.2.4      2020-03-10 [41]  CRAN (R 4.0.0)
108 viridisLite            0.3.0      2018-02-01 [99]  CRAN (R 4.0.0)
109 webshot                0.5.2      2019-11-22 [213] CRAN (R 4.0.0)
110 withr                  2.2.0      2020-04-20 [90]  CRAN (R 4.0.0)
111 xfun                   0.13       2020-04-13 [116] CRAN (R 4.0.0)
112 xml2                   1.3.2      2020-04-23 [122] CRAN (R 4.0.0)
113
114[1] /nix/store/xzd8h53xkyvfm3kvj5ab6znp685wi04w-r-car-3.0-7/library
115[2] /nix/store/mhr8zw9bmxarc3n821b83i0gz2j9zlrq-r-abind-1.4-5/library
116[3] /nix/store/hp86nhr0787vib3l8mkw0gf9nxwb45im-r-carData-3.0-3/library
117[4] /nix/store/vhw7s2h5ds6sp110z2yvilchv8j9jch5-r-lme4-1.1-23/library
118[5] /nix/store/987n8g0zy9sjvfvnsck1bkkcknw05yvb-r-boot-1.3-24/library
119[6] /nix/store/jxxxxyz4c1k5g3drd35gsrbjdg028d11-r-lattice-0.20-41/library
120[7] /nix/store/q9zfm5h53m8rd08xcsdcwaag31k4z1pf-r-MASS-7.3-51.5/library
121[8] /nix/store/kjkm50sr144yvrhl5axfgykbiy13pbmg-r-Matrix-1.2-18/library
122[9] /nix/store/8786z5lgy8h3akfjgj3yq5yq4s17rhjy-r-minqa-1.2.4/library
123[10] /nix/store/93wv3j0z1nzqp6fjsm9v7v8bf8d1xkm2-r-Rcpp-1.0.4.6/library
124[11] /nix/store/akfw6zsmawmz8lmjkww0rnqrazm4mqp0-r-nlme-3.1-147/library
125[12] /nix/store/rxs0d9bbn8qhw7wmkfb21yk5abp6lpq1-r-nloptr-1.2.2.1/library
126[13] /nix/store/8n0jfiqn4275i58qgld0dv8zdaihdzrk-r-RcppEigen-0.3.3.7.0/library
127[14] /nix/store/8vxrma33rhc96260zsi1jiw7dy3v2mm4-r-statmod-1.4.34/library
128[15] /nix/store/2y46pb5x9lh8m0hdmzajnx7sc1bk9ihl-r-maptools-0.9-9/library
129[16] /nix/store/iwf9nxx1v883wlv0p88q947hpz5lhfh7-r-foreign-0.8-78/library
130[17] /nix/store/rl9sjqply6rjbnz5k792ghm62ybv76px-r-sp-1.4-1/library
131[18] /nix/store/ws4bkzyv2vj5pyn1hgwyy6nlp48arz0n-r-mgcv-1.8-31/library
132[19] /nix/store/307dzxrmnqk4p86560a02r64x1fhhmxb-r-nnet-7.3-13/library
133[20] /nix/store/g2zpzkdb9hzkza1wpcbrk58119v1wyaf-r-pbkrtest-0.4-8.6/library
134[21] /nix/store/p0l503fr8960vld70w6ilmknxs5qwq77-r-magrittr-1.5/library
135[22] /nix/store/rmjpcaw3i446kwnjgcxcaid0yac36cj2-r-quantreg-5.55/library
136[23] /nix/store/10mzmnvc5jjgk2xzasia522pk60a30qz-r-MatrixModels-0.4-1/library
137[24] /nix/store/6qwdzvmnnmhjwdnvg2zmvv6wafd1vf91-r-SparseM-1.78/library
138[25] /nix/store/aa9c39a3yiqkh1h7pbngjlbr7czvc7yi-r-rio-0.5.16/library
139[26] /nix/store/2fx4vqlybgwp5rhhy6pssqx7h1a927fn-r-curl-4.3/library
140[27] /nix/store/k4m3fn1kqvvvn8y33kd57gq49hr3ar8y-r-data.table-1.12.8/library
141[28] /nix/store/651hfjylqzmsf565wyx474vyjny771gy-r-haven-2.2.0/library
142[29] /nix/store/a3rnz28irmqvmj8axj5x5j1am2c3gzs4-r-forcats-0.5.0/library
143[30] /nix/store/j8v4gzib137q2cml31hvvfkrc0f60pp5-r-ellipsis-0.3.0/library
144[31] /nix/store/xaswqlnamf4k8vwx0x3wav3l0x60sag0-r-rlang-0.4.5/library
145[32] /nix/store/dqm3xpix2jwhhhr67s6fgrwbw7hizap7-r-tibble-3.0.1/library
146[33] /nix/store/v7xfsq6d97wpn6m0hjrac78w5xawbr8a-r-cli-2.0.2/library
147[34] /nix/store/fikjasr98klhk9cf44x4lhi57vh3pmkg-r-assertthat-0.2.1/library
148[35] /nix/store/3fya6cd38vsqdj0gjb7bcsy00sirlyw1-r-crayon-1.3.4/library
149[36] /nix/store/payqi9bwh216rwhaq07jgc26l4fv1zsb-r-fansi-0.4.1/library
150[37] /nix/store/h6a61ghws7yrdxlg412xl1im37z5r28i-r-glue-1.4.0/library
151[38] /nix/store/y8mjbia1wbnq26dkigr0p3xxwrbzsc2r-r-lifecycle-0.2.0/library
152[39] /nix/store/kwaghh12cnifgvcbvlv2anx0hd5f4ild-r-pillar-1.4.3/library
153[40] /nix/store/k1phn8j10nni7gzvcgp0vc25dby6bb77-r-utf8-1.1.4/library
154[41] /nix/store/k3b77y8v7zsshpp1ccs8jwk2i2g4rm9a-r-vctrs-0.2.4/library
155[42] /nix/store/iibjmbh7vj0d0bfafz98yn29ymg43gkw-r-digest-0.6.25/library
156[43] /nix/store/aqsj4k3pgm80qk4jjg7sh3ac28n6alv0-r-pkgconfig-2.0.3/library
157[44] /nix/store/i7c5v8s4hd9rlqah3bbvy06yywjqwdgk-r-hms-0.5.3/library
158[45] /nix/store/2fyrk58cmcbrxid66rbwjli7y114lvrm-r-readr-1.3.1/library
159[46] /nix/store/163xq2g5nblqgh7qhvzb6mvgg6qdrirj-r-BH-1.72.0-3/library
160[47] /nix/store/dr27b6k49prwgrjs0v30b6mf5lxa36pk-r-clipr-0.7.0/library
161[48] /nix/store/bghvqg9mcaj2jkbwpy0di6c563v24acz-r-R6-2.4.1/library
162[49] /nix/store/nq8jdq7nlg9xns4xpgyj6sqv8p4ny1wz-r-tidyselect-1.0.0/library
163[50] /nix/store/zlwhf75qld7vmwx3d4bdws057ld4mqbp-r-purrr-0.3.4/library
164[51] /nix/store/0gbmmnbpqlr69l573ymkcx8154fvlaca-r-openxlsx-4.1.4/library
165[52] /nix/store/1m1q4rmwx56dvx9rdzfsfq0jpw3hw0yx-r-stringi-1.4.6/library
166[53] /nix/store/mhy5vnvbsl4q7dcinwx3vqlyywxphbfd-r-zip-2.0.4/library
167[54] /nix/store/88sp7f7q577i6l5jjanqiv5ak6nv5357-r-readxl-1.3.1/library
168[55] /nix/store/6q9zwivzalhmzdracc8ma932wirq8rl5-r-cellranger-1.1.0/library
169[56] /nix/store/jh2n6k2ancdzqych5ix8n4rq9w514qq9-r-rematch-1.0.1/library
170[57] /nix/store/22xjqikqd6q556absb5224sbx6q0kp0c-r-progress-1.2.2/library
171[58] /nix/store/9vp32wa1qvv6lkq6p70qlli5whrxzfbi-r-prettyunits-1.1.1/library
172[59] /nix/store/r9rhqb6fsk75shihmb7nagqb51pqwp0y-r-class-7.3-16/library
173[60] /nix/store/z1kad071y43wij1ml9lpghh7jbimmcli-r-cluster-2.1.0/library
174[61] /nix/store/i8wr965caf6j1rxs2dsvpzhlh4hyyb4y-r-codetools-0.2-16/library
175[62] /nix/store/8iglq3zr68a39hzswvzxqi2ffhpw9p51-r-KernSmooth-2.23-16/library
176[63] /nix/store/n3k50zv40i40drpdf8npbmy2y08gkr6w-r-rpart-4.1-15/library
177[64] /nix/store/b4r6adzcvpm8ivflsmis7ja7q4r5hkjy-r-spatial-7.3-11/library
178[65] /nix/store/zqg6hmrncl8ax3vn7z5drf4csddwnhcx-r-survival-3.1-12/library
179[66] /nix/store/4anrihkx11h8mzb269xdyi84yp5v7grl-r-tidyverse-1.3.0/library
180[67] /nix/store/945haq0w8nfm9ib7r0nfngn5lk2i15ix-r-broom-0.5.6/library
181[68] /nix/store/52viqxzrmxl7dk0zji293g5b0b9grwh8-r-backports-1.1.6/library
182[69] /nix/store/zp1k42sw2glqy51w4hnzsjs8rgi8xzx2-r-dplyr-0.8.5/library
183[70] /nix/store/mkjd98mnshch2pwnj6h31czclqdaph3f-r-plogr-0.2.0/library
184[71] /nix/store/kflrzax6y5pwfqwzgfvqz433a3q3hnhn-r-generics-0.0.2/library
185[72] /nix/store/xi1n5h5w17c33y6ax3dfhg2hgzjl9bxz-r-reshape2-1.4.4/library
186[73] /nix/store/vn63z92zkpbaxmmhzpb6mq2fvg0xa26h-r-plyr-1.8.6/library
187[74] /nix/store/wmpyxss67bj44rin7hlnr9qabx66p5hj-r-stringr-1.4.0/library
188[75] /nix/store/330qbgbvllwz3h0i2qidrlk50y0mbgph-r-tidyr-1.0.2/library
189[76] /nix/store/cx3x4pqb65l1mhss65780hbzv9jdrzl6-r-dbplyr-1.4.3/library
190[77] /nix/store/gsj49bp3hpw9jlli3894c49amddryqsq-r-DBI-1.1.0/library
191[78] /nix/store/kvymhwp4gac0343c2yi1qvdpavx4gdn2-r-ggplot2-3.3.0/library
192[79] /nix/store/knv51jvpairvibrkkq48b6f1l2pa1cv8-r-gtable-0.3.0/library
193[80] /nix/store/158dx0ddv20ikwag2860nlg9p3hbh1zc-r-isoband-0.2.1/library
194[81] /nix/store/fprs9rp1jlhxzj7fp6l79akyf8k3p7zd-r-testthat-2.3.2/library
195[82] /nix/store/0pmlnkyn0ir3k9bvxihi1r06jyl64w3i-r-evaluate-0.14/library
196[83] /nix/store/7210bjjqn5cjndxn5isnd4vip00xhkhy-r-pkgload-1.0.2/library
197[84] /nix/store/9a12ybd74b7dns40gcfs061wv7913qjy-r-desc-1.2.0/library
198[85] /nix/store/na9pb1apa787zp7vvyz1kzym0ywjwbj0-r-rprojroot-1.3-2/library
199[86] /nix/store/pa2n7bh61qxyarn5i2ynd62k6knb1np1-r-pkgbuild-1.0.6/library
200[87] /nix/store/1hxm1m7h4272zxk9bpsaq46mvnl0dbss-r-callr-3.4.3/library
201[88] /nix/store/bigvyk6ipglbiil93zkf442nv4y3xa1x-r-processx-3.4.2/library
202[89] /nix/store/370lr0wf7qlq0m72xnmasg2iahkp2n52-r-ps-1.3.2/library
203[90] /nix/store/rr72q61d8mkd42zc5fhcd2rqjghvc141-r-withr-2.2.0/library
204[91] /nix/store/9gw77p7fmz89fa8wi1d9rvril6hd4sxy-r-rstudioapi-0.11/library
205[92] /nix/store/9x4v4pbrgmykbz2801h77yz2l0nmm5nb-r-praise-1.0.0/library
206[93] /nix/store/pf8ssb0dliw5bzsncl227agc8przb7ic-r-scales-1.1.0/library
207[94] /nix/store/095z4wgjrxn63ixvyzrj1fm1rdv6ci95-r-farver-2.0.3/library
208[95] /nix/store/5aczj4s7i9prf5i32ik5ac5baqvjwdb1-r-labeling-0.3/library
209[96] /nix/store/wch26phipzz9gxd4vbr4fynh7v28349j-r-munsell-0.5.0/library
210[97] /nix/store/3w8fh756mszhsjx5fwgwydcpn8vkwady-r-colorspace-1.4-1/library
211[98] /nix/store/8cmaj81v2vm4f8p59ylbnsby8adkbmhd-r-RColorBrewer-1.1-2/library
212[99] /nix/store/h4x4ygax7gpz6f0c2v0xacr62080qwb8-r-viridisLite-0.3.0/library
213[100] /nix/store/qhx0i2nn5syb6vygdn8fdxgl7k56yj81-r-httr-1.4.1/library
214[101] /nix/store/lxnb4aniv02i4jhdvz02aaql1kznbpxb-r-jsonlite-1.6.1/library
215[102] /nix/store/13dcry4gad3vfwqzqb0ii4n06ybrxybr-r-mime-0.9/library
216[103] /nix/store/2can5l8gscc92a3bqlak8hfcg96v5hvf-r-openssl-1.4.1/library
217[104] /nix/store/piwsgxdz5w2ak8c6fcq0lc978qbxwdp1-r-askpass-1.1/library
218[105] /nix/store/3sj5h6dwa1l27d2hvdchclygk0pgffsr-r-sys-3.3/library
219[106] /nix/store/2z0p88g0c03gigl2ip60dlsfkdv1k30h-r-lubridate-1.7.8/library
220[107] /nix/store/1pkmj8nqjg2iinrkg2w0zkwq0ldc01za-r-modelr-0.1.6/library
221[108] /nix/store/bswkzvn8lczwbyw3y7n0p0qp2q472s0g-r-reprex-0.3.0/library
222[109] /nix/store/yid22gad8z49q52d225vfba2m4cgj2lx-r-fs-1.4.1/library
223[110] /nix/store/d185qiqaplm5br9fk1pf29y0srlabw83-r-rmarkdown-2.1/library
224[111] /nix/store/iszqviydsdj31c3ww095ndqy1ld3cibs-r-base64enc-0.1-3/library
225[112] /nix/store/i89wfw4cr0fz3wbd7cg44fk4dwz8b6h1-r-htmltools-0.4.0/library
226[113] /nix/store/qrl28laqwmhpwg3dpcf4nca8alv0px0g-r-knitr-1.28/library
227[114] /nix/store/jffaxc4a3bbf2g6ip0gdcya73dmg53mb-r-highr-0.8/library
228[115] /nix/store/717srph13qpnbzmgsvhx25q8pl51ivpj-r-markdown-1.1/library
229[116] /nix/store/mxqmyq3ybdfyc6p0anhfy2kfw0iz5k4n-r-xfun-0.13/library
230[117] /nix/store/b8g6hadva0359l6j1aq4dbvxlqf1acxc-r-yaml-2.2.1/library
231[118] /nix/store/rrl05vpv7cw58zi0k9ykm7m4rjb9gjv3-r-tinytex-0.22/library
232[119] /nix/store/2ziq8nzah6xy3dgmxgim9h2wszz1f89f-r-whisker-0.4/library
233[120] /nix/store/540wbw4p1g2qmnmbfk0rhvwvfnf657sj-r-rvest-0.3.5/library
234[121] /nix/store/n3prn77gd9sf3z4whqp86kghr55bf5w8-r-selectr-0.4-2/library
235[122] /nix/store/gv28yjk5isnglq087y7767xw64qa40cw-r-xml2-1.3.2/library
236[123] /nix/store/693czdcvkp6glyir0mi8cqvdc643whvc-r-gridExtra-2.3/library
237[124] /nix/store/3sykinp7lyy70dgzr0fxjb195nw864dv-r-future-1.17.0/library
238[125] /nix/store/bqi2l53jfxncks6diy0hr34bw8f86rvk-r-globals-0.12.5/library
239[126] /nix/store/dydyl209klklzh69w9q89f2dym9xycnp-r-listenv-0.8.0/library
240[127] /nix/store/lni0bi36r4swldkx7g4hql7gfz9b121b-r-gganimate-1.0.5/library
241[128] /nix/store/hh92jxs79kx7vxrxr6j6vin1icscl4k7-r-tweenr-1.0.1/library
242[129] /nix/store/0npx3srjnqgh7bib80xscjqvfyzjvimq-r-GGally-1.5.0/library
243[130] /nix/store/x5nzxklmacj6l162g7kg6ln9p25r3f17-r-reshape-0.8.8/library
244[131] /nix/store/q29z7ckdyhfmg1zlzrrg1nrm36ax756j-r-ggfortify-0.4.9/library
245[132] /nix/store/1rvm1w9iv2c5n22p4drbjq8lr9wa2q2r-r-cowplot-1.0.0/library
246[133] /nix/store/rp8jhnasaw1vbv5ny5zx0mw30zgcp796-r-ggrepel-0.8.2/library
247[134] /nix/store/wb7y931mm8nsj7w9xin83bvbaq8wvi4d-r-corrplot-0.84/library
248[135] /nix/store/gdzcqivfvgdrsz247v5kmnnw1v6p9c1p-r-rpart.plot-3.0.8/library
249[136] /nix/store/6yqg37108r0v22476cm2kv0536wyilki-r-caret-6.0-86/library
250[137] /nix/store/6fjdgcwgisiqz451sg5fszxnn9z8vxg6-r-foreach-1.5.0/library
251[138] /nix/store/c3ph5i341gk7jdinrkkqf6y631xli424-r-iterators-1.0.12/library
252[139] /nix/store/sjm1rxshlpakpxbrynfhsjnnp1sjvc3r-r-ModelMetrics-1.2.2.2/library
253[140] /nix/store/vgk4m131d057xglmrrb9rijhzdr2qhhp-r-pROC-1.16.2/library
254[141] /nix/store/bv1kvy1wc2jx3v55rzn3cg2qjbv7r8zp-r-recipes-0.1.10/library
255[142] /nix/store/001h42q4za01gli7avjxhq7shpv73n9k-r-gower-0.2.1/library
256[143] /nix/store/ssffpl6ydffqyn9phscnccxnj71chnzg-r-ipred-0.9-9/library
257[144] /nix/store/baliqip8m6p0ylqhqcgqak29d8ghral1-r-prodlim-2019.11.13/library
258[145] /nix/store/j4n2wsv98asw83qiffg6a74dymk8r2hl-r-lava-1.6.7/library
259[146] /nix/store/hf5wq5kpsf6p9slglq5iav09s4by0y5i-r-numDeriv-2016.8-1.1/library
260[147] /nix/store/s58hm38078mx4gyqffvv09zn575xn648-r-SQUAREM-2020.2/library
261[148] /nix/store/g63ydzd53586pvr9kdgk8kf5szq5f2bc-r-timeDate-3043.102/library
262[149] /nix/store/0jkarmlf1kjv4g8a3svkc7jfarpp77ny-r-mlr3-0.2.0/library
263[150] /nix/store/g1m0n1w7by213v773iyn7vnxr25pkf56-r-checkmate-2.0.0/library
264[151] /nix/store/fc2ah8cz2sj6j2jk7zldvjmsjn1yakpn-r-lgr-0.3.4/library
265[152] /nix/store/0i2hs088j1s0a6i61124my6vnzq8l27m-r-mlbench-2.1-1/library
266[153] /nix/store/vzcs6k21pqrli3ispqnvj5qwkv14srf5-r-mlr3measures-0.1.3/library
267[154] /nix/store/h2yqqaia46bk3b1d1a7bq35zf09p1b1a-r-mlr3misc-0.2.0/library
268[155] /nix/store/c9mrkc928cmsvvnib50l0jb8lsz59nyk-r-paradox-0.2.0/library
269[156] /nix/store/vqpbdipi4p4advl2vxrn765mmgcrabvk-r-uuid-0.1-4/library
270[157] /nix/store/xpclynxnfq4h9218gk4y62nmgyyga6zl-r-mlr3viz-0.1.1/library
271[158] /nix/store/7w6pld5vir3p9bybay67kq0qwl0gnx17-r-mlr3learners-0.2.0/library
272[159] /nix/store/ca50rp6ha5s51qmhb1gjlj62r19xfzxs-r-mlr3pipelines-0.1.3/library
273[160] /nix/store/9hg0xap4pir64mhbgq8r8cgrfjn8aiz5-r-mlr3filters-0.2.0/library
274[161] /nix/store/jgqcmfix0xxm3y90m8wy3xkgmqf2b996-r-rstan-2.19.3/library
275[162] /nix/store/mvv1gjyrrpvf47fn7a8x722wdwrf5azk-r-inline-0.3.15/library
276[163] /nix/store/zmkw51x4w4d1v1awcws0xihj4hnxfr09-r-loo-2.2.0/library
277[164] /nix/store/30xxalfwzxl05bbfvj5sy8k3ysys6z5y-r-matrixStats-0.56.0/library
278[165] /nix/store/fhkww2l0izx87bjnf0pl9ydl1wprp0xv-r-StanHeaders-2.19.2/library
279[166] /nix/store/aflck5pzxa8ym5q1dxchx5hisfmfghkr-r-tidybayes-2.0.3/library
280[167] /nix/store/jhlbhiv4fg0wsbxwjz8igc4hcg79vw94-r-arrayhelpers-1.1-0/library
281[168] /nix/store/fv089zrnvicnavbi08hnzqpi9g1z4inj-r-svUnit-1.0.3/library
282[169] /nix/store/xci2rgjizx1fyb33818jx5s1bgn8v8k6-r-coda-0.19-3/library
283[170] /nix/store/dch9asd38yldz0sdn8nsgk9ivjrkbhva-r-HDInterval-0.2.0/library
284[171] /nix/store/rs8dri2m5cqdmpiw187rvl4yhjn0jg2v-r-e1071-1.7-3/library
285[172] /nix/store/qs1zyh3sbvccgnqjzas3br6pak399zgc-r-pvclust-2.2-0/library
286[173] /nix/store/sh3zxvdazp7rkjn1iczrag1h2358ifm1-r-forecast-8.12/library
287[174] /nix/store/h67kaxqr2ppdpyj77wg5hm684jypznji-r-fracdiff-1.5-1/library
288[175] /nix/store/fh0z465ligbpqyam5l1fwiijc7334kbk-r-lmtest-0.9-37/library
289[176] /nix/store/0lnsbwfg0axr80h137q52pa50cllbjpf-r-zoo-1.8-7/library
290[177] /nix/store/p7k4s3ivf83dp2kcxr1cr0wlc1rfk6jx-r-RcppArmadillo-0.9.860.2.0/library
291[178] /nix/store/ssnxv5x6zid2w11v8k5yvnyxis6n1qfk-r-tseries-0.10-47/library
292[179] /nix/store/zrbskjwaz0bzz4v76j044d771m24g6h8-r-quadprog-1.5-8/library
293[180] /nix/store/2x3w5sjalrfm6hf1dxd951j8y94nh765-r-quantmod-0.4.17/library
294[181] /nix/store/7g55xshf49s9379ijm1zi1qnh1vbsifq-r-TTR-0.23-6/library
295[182] /nix/store/6ilyzph46q6ijyanq4p7f0ccyni0d7j0-r-xts-0.12-0/library
296[183] /nix/store/17xhqghcnqha7pwbf98dxsq1729slqd5-r-urca-1.3-0/library
297[184] /nix/store/722lyn0k8y27pj1alik56r4vpjnncd9z-r-swdft-1.0.0/library
298[185] /nix/store/36n0zgy10fsqcq76n0qmdwjxrwh7pn9n-r-xgboost-1.0.0.2/library
299[186] /nix/store/ac0ar7lf75qx84xsdjv6j02rkdgnhybz-r-ranger-0.12.1/library
300[187] /nix/store/i1ighkq42x10dirqmzgbx2mhbnz1ynkb-r-DALEX-1.2.0/library
301[188] /nix/store/28fqnhsfng1bkphl0wvr7lg5y3p6va46-r-iBreakDown-1.2.0/library
302[189] /nix/store/dpym77x9qc2ksr4mwjm3pb9ar1kvwhdl-r-ingredients-1.2.0/library
303[190] /nix/store/sp4d281w6dpr31as0xdjqizdx8hhb01q-r-DALEXtra-0.2.1/library
304[191] /nix/store/ckhp9kpmjcs0wxb113pxn25c2wip2d0n-r-ggdendro-0.1-20/library
305[192] /nix/store/f3k7dxj1dsmqri2gn0svq4c9fvvl9g7q-r-glmnet-3.0-2/library
306[193] /nix/store/l6ccj6mwkqybjvh6dr8qzalygp0i7jyb-r-shape-1.4.4/library
307[194] /nix/store/418mqfwlafh6984xld8lzhl7rv29qw68-r-reticulate-1.15/library
308[195] /nix/store/qwh982mgxd2mzrgbjk14irqbasywa1jk-r-rappdirs-0.3.1/library
309[196] /nix/store/6sxs76abll23c6372h6nf101wi8fcr4c-r-FactoMineR-2.3/library
310[197] /nix/store/39d2va10ydgyzddwr07xwdx11fwk191i-r-ellipse-0.4.1/library
311[198] /nix/store/4lxym5nxdn8hb7l8a566n5vg9paqcfi2-r-flashClust-1.01-2/library
312[199] /nix/store/wp161zbjjs41fq4kn4k3m244c7b8l2l2-r-leaps-3.1/library
313[200] /nix/store/irghsaplrpb3hg3y7j831bbklf2cqs6d-r-scatterplot3d-0.3-41/library
314[201] /nix/store/09ahkf50g1q9isxanbdykqgcdrp8mxl1-r-factoextra-1.0.7/library
315[202] /nix/store/zi9bq7amsgc6w2x7fvd62g9qxz69vjfm-r-dendextend-1.13.4/library
316[203] /nix/store/wcywb7ydglzlxg57jf354x31nmy63923-r-viridis-0.5.1/library
317[204] /nix/store/pvnpg4vdvv93pmwrlgmy51ihrb68j55f-r-ggpubr-0.2.5/library
318[205] /nix/store/qpapsc4l9pylzfhc72ha9d82hcbac41z-r-ggsci-2.9/library
319[206] /nix/store/h0zg4x3bmkc82ggx8h4q595ffckcqgx5-r-ggsignif-0.6.0/library
320[207] /nix/store/vn5svgbf8vsgv8iy8fdzlj0izp279q15-r-polynom-1.4-0/library
321[208] /nix/store/mc1mlsjx5h3gc8nkl7jlpd4vg145nk1z-r-lindia-0.9/library
322[209] /nix/store/z1k4c8lhabp9niwfg1xylg58pf99ld9r-r-orgutils-0.4-1/library
323[210] /nix/store/ybj4538v74wx4f1l064m0qn589vyjmzg-r-textutils-0.2-0/library
324[211] /nix/store/hhm5j0wvzjc0bfd53170bw8w7mij2wnh-r-latex2exp-0.4.0/library
325[212] /nix/store/njlv5mkxgjyx3x8p984nr84dwa2v1iqp-r-kableExtra-1.1.0/library
326[213] /nix/store/lf2sb84ylh259m421ljbj731a4prjhsl-r-webshot-0.5.2/library
327[214] /nix/store/n6b8ap54b78h8l70kyx9nvayp44rnfzf-r-printr-0.1/library
328[215] /nix/store/02g1v6d3ly8zylpckigwk6w3l1mx2i9d-r-microbenchmark-1.4-7/library
329[216] /nix/store/ri6qm0fp8cyx2qnysxjv2wsk0nndl1x9-r-webchem-0.5.0/library
330[217] /nix/store/cg95rqc1gmaqxf5kxja3cz8m5w4vl76l-r-RCurl-1.98-1.2/library
331[218] /nix/store/qbpinv148778fzdz8372x8gp34hspvy1-r-bitops-1.0-6/library
332[219] /nix/store/1g0lbrx6si76k282sxr9cj0mgknrw0lx-r-devtools-2.3.0/library
333[220] /nix/store/hnvww0128czlx6w8aipjn0zs7nvmvak9-r-covr-3.5.0/library
334[221] /nix/store/p4nv59przmb14sxi49jwqarkv0l40jsp-r-rex-1.2.0/library
335[222] /nix/store/vnysmc3vkgkligwah1zh9l4sahr533a8-r-lazyeval-0.2.2/library
336[223] /nix/store/d638w33ahybsa3sqr52fafvxs2b7w9x3-r-DT-0.13/library
337[224] /nix/store/35nqc34wy2nhd9bl7lv6wriw0l3cghsw-r-crosstalk-1.1.0.1/library
338[225] /nix/store/03838i63x5irvgmpgwj67ah0wi56k9d7-r-htmlwidgets-1.5.1/library
339[226] /nix/store/l4640jxlsjzqhw63c18fziar5vc0xyhk-r-promises-1.1.0/library
340[227] /nix/store/rxrb8p3dxzsg10v7yqaq5pi3y3gk6nqh-r-later-1.0.0/library
341[228] /nix/store/giprr32bl6k18b9n4qjckpf102flarly-r-git2r-0.26.1/library
342[229] /nix/store/bbkpkf44b13ig1pkz7af32kw5dzp12vb-r-memoise-1.1.0/library
343[230] /nix/store/m31vzssnfzapsapl7f8v4m15003lcc8r-r-rcmdcheck-1.3.3/library
344[231] /nix/store/hbiylknhxsin9hp9zaa6dwc2c9ai1mqx-r-sessioninfo-1.1.1/library
345[232] /nix/store/8vwlbx3s345gjccrkiqa6h1bm9wq4s9q-r-xopen-1.0.0/library
346[233] /nix/store/mjnwnlv60cn56ap0rrzvrkqlh5qisszx-r-remotes-2.1.1/library
347[234] /nix/store/1rq4zyzqymml7cc11q89rl5g514ml9na-r-roxygen2-7.1.0/library
348[235] /nix/store/2658mrn1hpkq0fv629rvags91qg65pbn-r-brew-1.0-6/library
349[236] /nix/store/nvjalws9lzva4pd4nz1z2131xsb9b5p6-r-commonmark-1.7/library
350[237] /nix/store/qx900vivd9s2zjrxc6868s92ljfwj5dv-r-rversions-2.0.1/library
351[238] /nix/store/1drg446wilq5fjnxkglxnnv8pbp1hllg-r-usethis-1.6.0/library
352[239] /nix/store/p3f3wa41d304zbs5cwvw7vy4j17zd6nq-r-gh-1.1.0/library
353[240] /nix/store/769g7jh93da8w15ad0wsbn2aqziwwx56-r-ini-0.3.1/library
354[241] /nix/store/p7kifw1l6z2zg68a71s4sdbfj8gdmnv5-r-rematch2-2.1.1/library
355[242] /nix/store/6zhdqip9ld9vl6pvifqcf4gsqy2f5wix-r-rethinking/library
356[243] /nix/store/496p28klmflihdkc83c8p1cywg85mgk4-r-mvtnorm-1.1-0/library
357[244] /nix/store/xb1zn7ab4nka7h1vm678ginzfwg4w9wf-r-dagitty-0.2-2/library
358[245] /nix/store/3zj4dkjbdwgf3mdsl9nf9jkicpz1nwgc-r-V8-3.0.2/library
359[246] /nix/store/qiqsh62w69b5xgj2i4wjamibzxxji0mf-r-tidybayes.rethinking/library
360[247] /nix/store/4j6byy1klyk4hm2k6g3657682cf3wxcj-R-4.0.0/lib/R/library

  1. Summer of 2020 ↩︎