Practice: Chp. 15

The answers to the questions in this chapter were originally developed with rstan, but running all the cells at once would consistently produce a crash (greater than half the time). The error was similar to:

 *** caught segfault ***
address 0x7fd4069fa9d0, cause 'memory not mapped'
An irrecoverable exception occurred. R is aborting now ...

Stan developers advise moving to CmdStan to avoid these crashes:

• Segfault when installing rstan 2.21.2 in R-4.0.3 on RHEL-7 · Issue #888 · stan-dev/rstan

• Getting started with CmdStanR • cmdstanr

Installation instructions for cmdstan and cmdstanr:

• 1 CmdStan Installation | CmdStan User’s Guide

• R Interface to CmdStan • cmdstanr

This environment includes cmdstan in conda:

library(cmdstanr)
set_cmdstan_path("/opt/conda/bin/cmdstan")

This is cmdstanr version 0.4.0

- Online documentation and vignettes at mc-stan.org/cmdstanr

- Use set_cmdstan_path() to set the path to CmdStan

- Use install_cmdstan() to install CmdStan

CmdStan path set to: /opt/conda/bin/cmdstan

After switching to CmdStan, the crashes completely went away. For more details on the process, see the commits near this change. The warning messages you see in this chapter may look otherwise unfamiliar.

source("iplot.R")
suppressPackageStartupMessages(library(rethinking))

15E1. Rewrite the Oceanic tools model (from Chapter 11) below so that it assumes measured error on the log population sizes of each society. You don’t need to fit the model to data. Just modify the mathematical formula below.

\[\begin{split} \begin{align} T_i & \sim Poisson(\mu_i) \\ log(\mu_i) & = \alpha + \beta log(P_i) \\ \alpha & \sim Normal(0, 1.5) \\ \beta & \sim Normal(0, 1) \end{align} \end{split}\]

Answer. Assuming a \(P_{SE,i}\) is available, though it is not in the data:

\[\begin{split} \begin{align} T_i & \sim Poisson(\mu_i) \\ log(\mu_i) & = \alpha + \beta log(P_{TRUE,i}) \\ P_{OBS,i} & \sim Normal(P_{TRUE,i}, P_{SE,i}) \\ \alpha & \sim Normal(0, 1.5) \\ \beta & \sim Normal(0, 1) \end{align} \end{split}\]

15E2. Rewrite the same model so that it allows imputation of missing values for log population. There aren’t any missing values in the variable, but you can still write down a model formula that would imply imputation, if any values were missing.

Answer. As a challenge we’ll add this to model that also considers measurement error, assuming either both the observed value and the measurement are present, or neither. The priors could be much better:

\[\begin{split} \begin{align} T_i & \sim Poisson(\mu_i) \\ log(\mu_i) & = \alpha + \beta log(P_{TRUE,i}) \\ P_{OBS,i} & \sim Normal(P_{TRUE,i}, P_{SE,i}) \\ P_{OBS,i} & \sim Normal(\nu, \sigma_B) \\ \alpha & \sim Normal(0, 1.5) \\ \beta & \sim Normal(0, 1) \\ \nu & \sim Normal(1e6, 1e6) \\ \sigma_B & \sim Exponential(0.0001) \end{align} \end{split}\]

15M1. Using the mathematical form of the imputation model in the chapter, explain what is being assumed about how the missing values were generated.

ERROR. It’s not clear what the author means by ‘the imputation model’; there are several models in the text.

Answer. In section 15.2.2. the first imputation model assumes a normal distribution for brain sizes. This is somewhat unreasonable, as explained in the text, because these values should be bounded between zero and one.

In the second model for this primate example we still assume a normal distribution, but now correlated with body mass.

15M2. Reconsider the primate milk missing data example from the chapter. This time, assign \(B\) a distribution that is properly bounded between zero and 1. A beta distribution, for example, is a good choice.

Answer. This question requires a start= list to run at all, which isn’t suggested for these imputation models until a later question. We use the dbeta2 parameterization provided by ulam to make the parameters a bit more understandable (as opposed to dbeta). Fitting the model:

data(milk)
d <- milk
d$neocortex.prop <- d$neocortex.perc / 100
d$logmass <- log(d$mass)
stdB = standardize(d$neocortex.prop)
dat_list <- list(
  K = standardize(d$kcal.per.g),
  B = d$neocortex.prop,
  M = standardize(d$logmass),
  center = attr(stdB, "scaled:center"),
  scale = attr(stdB, "scaled:scale")
)

m_milk_beta <- ulam(
  alist(
    K ~ dnorm(mu, sigma),
    mu <- a + bB*B + bM*M,
    B ~ dbeta2(p, theta),
    a ~ dnorm(0, 0.5),
    c(bB, bM) ~ dnorm(0, 0.5),
    p ~ dunif(0, 1),
    theta ~ dunif(0, 1000),
    sigma ~ dexp(1)
  ), cmdstan=TRUE, data=dat_list, chains=4, cores=4, start=list(B_impute=rep(0.7,12))
)

Found 12 NA values in B and attempting imputation.

Without bounds on the parameters to the beta distribution, this model produces a significant number of divergent transitions. To avoid these errors, we are using uniform priors on the two parameters to the beta distribution. Checking the pairs and traceplot:

sel_pars=c("a", "bM", "bB", "sigma", "p", "theta", "B_impute[1]", "B_impute[5]")
iplot(function() {
  pairs(m_milk_beta@stanfit, pars=sel_pars)
})
iplot(function() {
  traceplot(m_milk_beta, pars=sel_pars)
})

[1] 1000
[1] 1
[1] 1000

The precis results:

display(precis(m_milk_beta, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m_milk_beta, depth=3), main="m_milk_beta")
}, ar=2)

             mean        sd          5.5%       94.5%        n_eff    Rhat4    
a            -0.05116851  0.29606393 -0.5052037   0.41696665 1869.260 0.9998132
bM           -0.29726985  0.17142750 -0.5573353  -0.02801925 2800.385 0.9981806
bB            0.08818411  0.40676501 -0.5564030   0.75273561 1911.830 1.0000544
p             0.67494139  0.01423122  0.6519743   0.69658273 1675.502 1.0004415
theta        70.11188688 23.14271794 37.6938177 109.83115043 1310.476 1.0023321
sigma         0.98346024  0.13719629  0.7946166   1.22593809 2584.627 0.9993949
B_impute[1]   0.67572783  0.06029584  0.5795165   0.76839078 2647.676 0.9996259
B_impute[2]   0.67345591  0.06255013  0.5744566   0.76941098 2464.468 0.9982780
B_impute[3]   0.67397771  0.06306257  0.5728357   0.76951075 2406.441 1.0003588
B_impute[4]   0.67462829  0.05883192  0.5788101   0.76225394 2566.681 1.0000882
B_impute[5]   0.67610284  0.05910520  0.5767991   0.76550498 2376.502 1.0004694
B_impute[6]   0.67340057  0.05903189  0.5767701   0.76385963 2268.317 1.0000376
B_impute[7]   0.67489505  0.06102367  0.5802512   0.76690901 2430.159 1.0000165
B_impute[8]   0.67578277  0.06216776  0.5742970   0.77365148 2625.762 1.0022308
B_impute[9]   0.67349779  0.05960898  0.5709987   0.76403209 2393.433 0.9984447
B_impute[10]  0.67375332  0.05872196  0.5730508   0.76118914 2501.458 0.9995646
B_impute[11]  0.67475707  0.06160050  0.5691120   0.77090342 2628.574 0.9985411
B_impute[12]  0.67560322  0.06132239  0.5741296   0.76897318 2721.813 0.9997512

The theta parameter to the beta distribution is highly uncertain because there are simply too many values that are compatible with these observations. The imputed values are essentially nothing more than the mean of the non-NA observations of brain size. That is, any gentle tilt of the imputed values we’d expect from using the observed mass seems to have disappeared.

Let’s zoom in on some parameters of interest:

iplot(function() {
  plot(precis(m_milk_beta, depth=3), main="m_milk_beta", pars=c("a", "bM", "bB", "sigma"))
}, ar=4.5)

Compare these results to the figure produced from R code 15.20 in the chapter. The inference for bM has shrunk even further. It’s likely bB has shrunk as well, though in this new model bB is multiplied by a probability rather than a standardized value so we can expect it to be smaller.

15M3. Repeat the divorce data measurement error models, but this time double the standard errors. Can you explain how doubling the standard errors impacts inference?

Answer. Let’s start by reproducing the results from the chapter, for the sake of comparison:

data(WaffleDivorce)
d <- WaffleDivorce

dlist <- list(
  D_obs = standardize(d$Divorce),
  D_sd = d$Divorce.SE / sd(d$Divorce),
  M = standardize(d$Marriage),
  A = standardize(d$MedianAgeMarriage),
  N = nrow(d)
)
m15.1 <- ulam(
  alist(
    D_obs ~ dnorm(D_true, D_sd),
    vector[N]:D_true ~ dnorm(mu, sigma),
    mu <- a + bA * A + bM * M,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  cmdstan = TRUE, data = dlist, chains = 4, cores = 4
)
display(precis(m15.1, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m15.1, depth=3), main="m15.1")
}, ar=1)

dlist <- list(
  D_obs = standardize(d$Divorce),
  D_sd = d$Divorce.SE / sd(d$Divorce),
  M_obs = standardize(d$Marriage),
  M_sd = d$Marriage.SE / sd(d$Marriage),
  A = standardize(d$MedianAgeMarriage),
  N = nrow(d)
)
m15.2 <- ulam(
  alist(
    D_obs ~ dnorm(D_true, D_sd),
    vector[N]:D_true ~ dnorm(mu, sigma),
    mu <- a + bA * A + bM * M_true[i],
    M_obs ~ dnorm(M_true, M_sd),
    vector[N]:M_true ~ dnorm(0, 1),
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  cmdstan = TRUE, data = dlist, chains = 4, cores = 4
)
display(precis(m15.2, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m15.2, depth=3), main="m15.2")
}, ar=0.5)

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-1367da2ab5.stan', line 24, column 4 to column 34)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.2 seconds.
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 0.2 seconds.
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.2 seconds.
Total execution time: 0.5 seconds.

           mean        sd         5.5%         94.5%       n_eff     Rhat4    
D_true[1]   1.15316574 0.37012833  0.563892475  1.74989110 2203.5811 0.9994987
D_true[2]   0.65546808 0.56020094 -0.236100900  1.52888410 2240.3475 0.9984558
D_true[3]   0.41089610 0.31577544 -0.073016246  0.90617627 2435.4814 0.9993676
D_true[4]   1.40658957 0.46939771  0.659178010  2.14012165 2182.6306 0.9999575
D_true[5]  -0.90023786 0.12403659 -1.092046050 -0.70316698 2380.4359 0.9983087
D_true[6]   0.65046162 0.41016744 -0.002664295  1.30623595 2420.6372 0.9991410
D_true[7]  -1.37022663 0.34764760 -1.925934750 -0.81699102 2471.9079 0.9989182
D_true[8]  -0.34879128 0.48941110 -1.136557650  0.41691948 2333.6308 0.9995918
D_true[9]  -1.89043685 0.59592192 -2.822572700 -0.94111641 1311.5441 1.0016078
D_true[10] -0.61307923 0.17066512 -0.889232460 -0.34026281 2688.5664 0.9985730
D_true[11]  0.76335219 0.29959721  0.285579825  1.24426495 1553.9773 1.0006682
D_true[12] -0.55757996 0.46366725 -1.297839900  0.16153468 1923.6825 0.9996397
D_true[13]  0.17165800 0.49115790 -0.653381430  0.92907183 1253.2028 1.0014375
D_true[14] -0.86955366 0.22710730 -1.245080150 -0.51252277 2193.1867 1.0003510
D_true[15]  0.54938972 0.29721613  0.077860388  1.02241940 2270.1418 0.9993618
D_true[16]  0.28200710 0.37641566 -0.304512715  0.89413573 2167.3267 1.0002984
D_true[17]  0.50583186 0.42135858 -0.185592965  1.17446785 2322.7380 1.0016654
D_true[18]  1.25987628 0.35105252  0.710302160  1.83764880 1963.9749 0.9994051
D_true[19]  0.43917238 0.38441929 -0.175468310  1.06542200 2588.3619 1.0002618
D_true[20]  0.40658838 0.53443690 -0.466622880  1.28072990 1453.2413 1.0001826
D_true[21] -0.55940793 0.32762304 -1.072680550 -0.02175856 2670.8494 0.9987141
D_true[22] -1.10021536 0.25896681 -1.509204500 -0.67054736 2200.2264 0.9994890
D_true[23] -0.27700879 0.25770450 -0.688211900  0.14914577 1805.6023 1.0000475
D_true[24] -1.00828902 0.29680944 -1.476085900 -0.53257590 2036.4173 0.9998370
D_true[25]  0.42644139 0.41907982 -0.211422840  1.08312165 2185.6700 0.9981407
D_true[26] -0.03560633 0.31544041 -0.545365355  0.46128014 2721.7211 1.0001653
D_true[27] -0.04365522 0.49721603 -0.856300380  0.72195147 2272.9411 1.0000432
D_true[28] -0.15249989 0.40306779 -0.804953285  0.49303272 2249.3939 0.9999010
D_true[29] -0.26773353 0.52792144 -1.078230700  0.59383324 1930.6701 0.9989877
D_true[30] -1.80215750 0.23039666 -2.166020450 -1.43927790 2348.1470 0.9996754
D_true[31]  0.17117231 0.41369859 -0.492660205  0.81272506 2847.4546 0.9998587
D_true[32] -1.65761140 0.16596683 -1.924196400 -1.39210395 2585.1019 0.9988227
D_true[33]  0.11941390 0.24567926 -0.293404410  0.50143680 2409.6102 0.9996998
D_true[34] -0.07930125 0.50971164 -0.880127020  0.72891709 1505.7317 1.0007523
D_true[35] -0.11825308 0.22628999 -0.475089525  0.23721810 2610.6556 0.9987296
D_true[36]  1.27163602 0.41808670  0.608068470  1.93706190 1977.8994 0.9995361
D_true[37]  0.23718802 0.35413735 -0.337260090  0.79951773 2256.8464 0.9993212
D_true[38] -1.04186103 0.22394803 -1.383967500 -0.68586310 2273.2483 0.9990624
D_true[39] -0.92725481 0.55033599 -1.788967950 -0.03705807 1940.8838 1.0008087
D_true[40] -0.67815066 0.32289060 -1.216340800 -0.16710560 3102.1919 0.9989761
D_true[41]  0.24909346 0.55565165 -0.652739925  1.12766820 2131.9232 0.9993212
D_true[42]  0.74622692 0.35144307  0.190627455  1.32082000 1965.8029 0.9985666
D_true[43]  0.19387862 0.18360872 -0.103154845  0.48404413 2380.4923 0.9996526
D_true[44]  0.78671897 0.42430177  0.120739375  1.43915070 1982.9018 1.0017695
D_true[45] -0.39407562 0.49855133 -1.178090150  0.42149137 2270.7442 0.9995726
D_true[46] -0.37486130 0.26414860 -0.797083485  0.04861196 2382.6522 1.0019703
D_true[47]  0.13298003 0.30549031 -0.351561445  0.62543107 3017.6014 1.0000734
D_true[48]  0.56240853 0.45162042 -0.133742120  1.27482755 2582.6418 0.9988020
D_true[49] -0.63920292 0.27712076 -1.086768250 -0.20110425 2267.7809 0.9997658
D_true[50]  0.83384483 0.58493775 -0.099355149  1.75366760 1771.8106 0.9999210
a          -0.05526773 0.09543933 -0.202567165  0.09742136 1246.1718 1.0013101
bA         -0.61378817 0.16146817 -0.866782565 -0.35621728  972.5183 1.0030096
bM          0.05121898 0.17144336 -0.229405235  0.32937462  897.8376 1.0006954
sigma       0.58698448 0.10438335  0.433398140  0.76308559  612.9772 1.0054102

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-131e5c493f.stan', line 28, column 4 to column 34)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 finished in 0.4 seconds.
Chain 2 finished in 0.4 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 0.4 seconds.
Chain 4 finished in 0.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.4 seconds.
Total execution time: 0.6 seconds.

           mean        sd         5.5%        94.5%        n_eff     Rhat4    
D_true[1]   1.13388681 0.3700038   0.56357853  1.75856995  2278.124  0.9992381
D_true[2]   0.73414347 0.5577625  -0.17421989  1.63166145  2252.998  1.0010008
D_true[3]   0.42734069 0.3231727  -0.07188098  0.95950344  2668.640  0.9997068
D_true[4]   1.45481411 0.4702182   0.71245354  2.21338270  2868.213  0.9987640
D_true[5]  -0.89763664 0.1295431  -1.10387870 -0.68533850  3201.665  0.9998565
D_true[6]   0.68921785 0.4021074   0.05222334  1.32158590  2784.080  0.9995133
D_true[7]  -1.35656769 0.3304125  -1.89445015 -0.84332208  2641.390  1.0000333
D_true[8]  -0.27556671 0.4762452  -1.07545405  0.47576607  2705.554  0.9986518
D_true[9]  -1.75012361 0.5978167  -2.68661005 -0.79774212  2054.575  0.9996723
D_true[10] -0.62617242 0.1671956  -0.89201643 -0.34875431  3283.469  0.9985267
D_true[11]  0.76679802 0.2793681   0.33288921  1.21969100  2959.183  0.9987223
D_true[12] -0.44602907 0.4979900  -1.23853010  0.35861287  2241.260  1.0003170
D_true[13]  0.21302547 0.4914136  -0.57745285  0.98469074  1315.314  0.9985382
D_true[14] -0.86304926 0.2388686  -1.23765330 -0.49391420  3253.555  0.9994466
D_true[15]  0.53331720 0.2949552   0.07124441  1.02205145  2562.618  0.9995080
D_true[16]  0.28705500 0.3681874  -0.31008845  0.87384413  3002.179  0.9988337
D_true[17]  0.51949443 0.4096071  -0.11601631  1.17522655  3048.076  0.9986472
D_true[18]  1.24400665 0.3456773   0.71768575  1.80893190  2452.218  0.9997295
D_true[19]  0.42481563 0.3630439  -0.16594681  1.02176795  3277.172  0.9990667
D_true[20]  0.24748045 0.5447844  -0.58244826  1.15329695  1194.414  0.9999912
D_true[21] -0.54007075 0.3182891  -1.04094690 -0.03319127  2750.883  0.9993155
D_true[22] -1.10845255 0.2712931  -1.53989645 -0.67339837  3019.962  0.9991770
D_true[23] -0.29044344 0.2511989  -0.68492705  0.12705940  2955.664  0.9984704
D_true[24] -1.02037136 0.2875369  -1.47242130 -0.55460197  2869.834  0.9999754
D_true[25]  0.39554853 0.4133142  -0.26609771  1.08142880  2773.046  0.9988272
D_true[26] -0.03949349 0.3096815  -0.53393152  0.45476229  3102.018  0.9994643
D_true[27] -0.04940019 0.4957922  -0.84797165  0.72680382  3538.592  0.9993781
D_true[28] -0.16155873 0.3894544  -0.78204267  0.45288288  3318.422  0.9985426
D_true[29] -0.31034618 0.4865836  -1.07818190  0.47469244  3085.004  0.9988969
D_true[30] -1.80288174 0.2398179  -2.18279870 -1.42264205  2738.877  0.9997140
⋮          ⋮           ⋮          ⋮           ⋮            ⋮         ⋮        
M_true[25] -0.15304036 0.37391687 -0.74601636  0.446288155 3481.9502 0.9989487
M_true[26] -0.38122032 0.19802365 -0.70368617 -0.048822922 2516.0505 0.9991729
M_true[27] -0.32486394 0.51750346 -1.14167905  0.508964565 3832.8506 0.9989923
M_true[28] -0.15284528 0.35210218 -0.70107837  0.424813170 2722.2272 0.9987227
M_true[29] -0.71017473 0.40327094 -1.35959245 -0.082660369 3404.6344 1.0007958
M_true[30] -1.38435092 0.15037667 -1.61385100 -1.147789750 3500.4376 0.9988590
M_true[31]  0.07443541 0.43386295 -0.61503308  0.755907205 3406.6466 0.9992482
M_true[32] -0.86641755 0.12214125 -1.06473940 -0.677123520 3172.2296 0.9987598
M_true[33]  0.06795647 0.24499380 -0.31880562  0.458021360 3794.1590 0.9993659
M_true[34]  0.95875734 0.60983021 -0.02476590  1.929871950 2558.0188 0.9990821
M_true[35] -0.82016864 0.15483225 -1.05949840 -0.568463975 3476.5566 0.9990808
M_true[36]  0.89210789 0.30801735  0.40521868  1.364540700 2932.2859 0.9987709
M_true[37] -0.28271176 0.27705706 -0.73222911  0.151692345 3522.0454 0.9987309
M_true[38] -1.19894758 0.12378671 -1.39117110 -1.001401400 3116.4656 0.9994999
M_true[39] -1.00450476 0.48378902 -1.78827740 -0.257685960 2789.2188 1.0005413
M_true[40] -0.50229769 0.31186812 -0.99824275  0.005136526 3681.6630 0.9985832
M_true[41]  0.02800274 0.57648471 -0.87603796  0.962350380 3267.1616 0.9987050
M_true[42] -0.15444747 0.20894814 -0.48388113  0.181049900 3259.7368 0.9997349
M_true[43]  0.35034691 0.15872207  0.09996562  0.606964080 2154.1477 1.0011451
M_true[44]  1.95740181 0.43951431  1.25093600  2.669159350 2219.4832 0.9999532
M_true[45] -0.67562754 0.54243173 -1.54826475  0.178041165 3338.7545 0.9999538
M_true[46]  0.07985319 0.21977724 -0.26712863  0.434721930 3107.4275 0.9994109
M_true[47]  0.31756674 0.25685015 -0.09188301  0.737931210 2500.9707 1.0018663
M_true[48]  0.45830875 0.40750127 -0.20002161  1.113281800 4485.8638 0.9984663
M_true[49] -0.73976426 0.19966476 -1.06401035 -0.418241215 3363.9107 0.9993646
M_true[50]  1.28172338 0.71550865  0.14640638  2.397851400 2960.1400 1.0003424
a          -0.04050490 0.09677744 -0.19527547  0.109647935 2248.7523 0.9999561
bA         -0.54366668 0.16430879 -0.80177939 -0.280791395 1125.2827 0.9997204
bM          0.19685452 0.21206789 -0.13557063  0.540942560  769.6873 1.0010293
sigma       0.55713776 0.11399910  0.38579479  0.735343595  600.7142 1.0001022

Doubling the standard error:

d$doubleSE = 2*d$Divorce.SE

dlist <- list(
  D_obs = standardize(d$Divorce),
  D_sd = d$doubleSE / sd(d$Divorce),
  M = standardize(d$Marriage),
  A = standardize(d$MedianAgeMarriage),
  N = nrow(d)
)
m.double.se.1 <- ulam(
  alist(
    D_obs ~ dnorm(D_true, D_sd),
    vector[N]:D_true ~ dnorm(mu, sigma),
    mu <- a + bA * A + bM * M,
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  cmdstan = TRUE, data = dlist, chains = 4, cores = 4
)
display(precis(m.double.se.1, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m.double.se.1, depth=3), main="m.double.se.1")
}, ar=1)

dlist <- list(
  D_obs = standardize(d$Divorce),
  D_sd = d$doubleSE / sd(d$Divorce),
  M_obs = standardize(d$Marriage),
  M_sd = d$Marriage.SE / sd(d$Marriage),
  A = standardize(d$MedianAgeMarriage),
  N = nrow(d)
)
m.double.se.2 <- ulam(
  alist(
    D_obs ~ dnorm(D_true, D_sd),
    vector[N]:D_true ~ dnorm(mu, sigma),
    mu <- a + bA * A + bM * M_true[i],
    M_obs ~ dnorm(M_true, M_sd),
    vector[N]:M_true ~ dnorm(0, 1),
    a ~ dnorm(0, 0.2),
    bA ~ dnorm(0, 0.5),
    bM ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ),
  cmdstan = TRUE, data = dlist, chains = 4, cores = 4
)
display(precis(m.double.se.2, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m.double.se.2, depth=3), main="m.double.se.2")
}, ar=0.5)

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 finished in 0.6 seconds.
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 0.7 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 0.8 seconds.
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 1.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.8 seconds.
Total execution time: 1.1 seconds.

Warning: 57 of 2000 (3.0%) transitions ended with a divergence.
This may indicate insufficient exploration of the posterior distribution.
Possible remedies include: 
  * Increasing adapt_delta closer to 1 (default is 0.8) 
  * Reparameterizing the model (e.g. using a non-centered parameterization)
  * Using informative or weakly informative prior distributions 

           mean        sd         5.5%         94.5%       n_eff     Rhat4   
D_true[1]   0.33304365 0.23489498  0.008330816  0.74066968 394.14731 1.010892
D_true[2]   0.65304745 0.32626561  0.146140005  1.16528795 221.82476 1.029617
D_true[3]   0.04694734 0.20837050 -0.251312265  0.39317153 674.61701 1.007555
D_true[4]   1.14917817 0.32433455  0.647242180  1.64668585 222.51949 1.027255
D_true[5]  -0.66167981 0.16612281 -0.937859155 -0.42221868 227.96174 1.007016
D_true[6]   0.27501919 0.25967921 -0.125722565  0.68801150 272.16273 1.026367
D_true[7]  -1.09468264 0.23998806 -1.457664950 -0.69431590 659.72490 1.006402
D_true[8]  -0.23012897 0.29566653 -0.680935820  0.21537801 285.39505 1.019242
D_true[9]  -2.10517061 0.42893525 -2.779048750 -1.42299100 329.24289 1.013437
D_true[10] -0.51095521 0.19138718 -0.811019810 -0.20767308 447.98941 1.005129
D_true[11]  0.15060838 0.22453924 -0.166040080  0.54561506 186.21052 1.019607
D_true[12] -0.29425871 0.38525591 -0.883182085  0.29825835 231.41162 1.021230
D_true[13]  1.57670940 0.35194272  1.003891800  2.10983240 228.29256 1.026781
D_true[14] -0.75132001 0.19992396 -1.074373650 -0.44377359 830.85457 1.003682
D_true[15]  0.09389170 0.21641041 -0.207696275  0.45381828 511.55607 1.009215
D_true[16]  0.28965894 0.21788624 -0.058580957  0.63657137 475.52376 1.013998
D_true[17]  0.52615361 0.24269913  0.154371500  0.91296724 513.25191 1.013408
D_true[18]  0.68449729 0.25084706  0.310224045  1.08874430 334.58605 1.017593
D_true[19]  0.01225383 0.20744579 -0.306061865  0.36463732 690.28613 1.008829
D_true[20] -0.62358209 0.35533088 -1.152001150 -0.04779766 262.35402 1.008304
D_true[21] -0.83272716 0.23499291 -1.177065100 -0.44448614 659.23075 1.006954
D_true[22] -1.55373343 0.28068948 -1.968546050 -1.06227885 304.13965 1.016618
D_true[23] -0.47808708 0.22764808 -0.825362060 -0.09514581 357.73096 1.006915
D_true[24] -0.57208941 0.26956389 -1.008740950 -0.14090562 292.57534 1.005571
D_true[25] -0.01361509 0.21609970 -0.334490615  0.34207745 454.36716 1.010226
D_true[26]  0.02016651 0.22739609 -0.333545725  0.38074079 584.20173 1.004521
D_true[27] -0.03078444 0.23486930 -0.384448400  0.33763212 521.32886 1.001816
D_true[28]  0.16667402 0.23210542 -0.198637320  0.53399456 646.09569 1.002031
D_true[29] -0.67464394 0.23688893 -1.032188150 -0.29211116 464.93727 1.006057
D_true[30] -1.33643272 0.23838853 -1.718854950 -0.96071635 277.96761 1.009491
D_true[31]  0.02982638 0.21774779 -0.286542135  0.37637218 887.72792 1.009750
D_true[32] -1.54633596 0.23964126 -1.925245050 -1.16406505 460.50985 1.012125
D_true[33]  0.07615650 0.20495008 -0.247307220  0.40731202 656.91387 1.005718
D_true[34]  0.61459842 0.35815178  0.063259779  1.16932830 207.01047 1.027687
D_true[35] -0.38325011 0.21809922 -0.702262005 -0.01397761 325.31119 1.007210
D_true[36]  0.95548072 0.27676839  0.518414340  1.40084490 436.74823 1.016594
D_true[37] -0.14173873 0.21440571 -0.458163710  0.21833585 742.30388 1.004060
D_true[38] -0.93654047 0.22480911 -1.289034400 -0.58345277 340.92831 1.007880
D_true[39] -1.48443219 0.29457772 -1.915089450 -1.01102590 397.90610 1.010878
D_true[40] -0.43822708 0.20381085 -0.775971790 -0.13276657 607.05383 1.002208
D_true[41]  0.10808187 0.22693556 -0.215319455  0.46684571 630.85316 1.009131
D_true[42]  0.30209039 0.23768831 -0.056737318  0.68326158 495.03733 1.005753
D_true[43]  0.35420508 0.20181229  0.028765980  0.67940601 400.20964 1.014392
D_true[44]  1.76058394 0.41574203  1.073125000  2.38668940 207.70426 1.033342
D_true[45] -0.74913597 0.24394977 -1.119639900 -0.34368112 428.55690 1.008449
D_true[46] -0.29159181 0.21400086 -0.631211550  0.05912485 630.70026 1.006663
D_true[47]  0.04092238 0.20723819 -0.288842365  0.37857958 408.04500 1.016986
D_true[48]  0.53672435 0.24881215  0.156925495  0.92921050 554.88683 1.015150
D_true[49] -0.44236666 0.21797338 -0.782920240 -0.11626839 429.86132 1.003001
D_true[50]  1.40768571 0.47761906  0.663213180  2.18029140 192.52649 1.030684
a          -0.11890192 0.09835123 -0.277704255  0.03502495 174.19172 1.036387
bA         -0.63274604 0.15566461 -0.883606615 -0.38579555 276.10146 1.010102
bM          0.21302959 0.18732013 -0.085652424  0.50289488 186.78463 1.021277
sigma       0.16954026 0.09426135  0.047733917  0.33943000  58.93479 1.028302

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 finished in 1.0 seconds.
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 finished in 1.0 seconds.
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 1.1 seconds.
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 7.9 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.7 seconds.
Total execution time: 7.9 seconds.

Warning: 139 of 2000 (7.0%) transitions ended with a divergence.
This may indicate insufficient exploration of the posterior distribution.
Possible remedies include: 
  * Increasing adapt_delta closer to 1 (default is 0.8) 
  * Reparameterizing the model (e.g. using a non-centered parameterization)
  * Using informative or weakly informative prior distributions 

495 of 2000 (25.0%) transitions hit the maximum treedepth limit of 10 or 2^10-1 leapfrog steps.
Trajectories that are prematurely terminated due to this limit will result in slow exploration.
Increasing the max_treedepth limit can avoid this at the expense of more computation.
If increasing max_treedepth does not remove warnings, try to reparameterize the model.

           mean         sd        5.5%        94.5%       n_eff      Rhat4    
D_true[1]   0.345851584 0.2748155 -0.05967765  0.7738881   62.47068  1.106930 
D_true[2]   0.567215294 0.3126899  0.10402547  1.1150698  249.72860  1.022491 
D_true[3]   0.050642343 0.2308589 -0.27524388  0.4543758  130.80984  1.029680 
D_true[4]   1.143340566 0.3380874  0.61033224  1.6881258  170.84009  1.035832 
D_true[5]  -0.655438297 0.1582227 -0.92313702 -0.4279900  101.74412  1.035331 
D_true[6]   0.294572463 0.2695177 -0.09275772  0.7442568  145.50760  1.029644 
D_true[7]  -1.113544847 0.2368149 -1.45630885 -0.7019197  227.07747  1.004008 
D_true[8]  -0.283449256 0.3089853 -0.73119906  0.2687901  195.20008  1.011288 
D_true[9]  -2.087127722 0.4747847 -2.79692805 -1.3077067  155.51209  1.023207 
D_true[10] -0.516124996 0.1988460 -0.83279148 -0.2331847  195.71820  1.020310 
D_true[11]  0.152044710 0.2291893 -0.16334631  0.5591090   85.35617  1.050246 
D_true[12] -0.345921864 0.3594281 -0.86197326  0.2869306  165.11168  1.005692 
D_true[13]  1.521350669 0.3937463  0.84781970  2.0992540  153.46867  1.031202 
D_true[14] -0.772283591 0.1866238 -1.07144530 -0.4882894  364.21776  1.014608 
D_true[15]  0.076649667 0.2070942 -0.22319993  0.4078297  166.28805  1.040489 
D_true[16]  0.299415073 0.2467922 -0.08273458  0.6845722  219.57100  1.030695 
D_true[17]  0.535831064 0.2539483  0.13588757  0.9424322  293.56599  1.037025 
D_true[18]  0.693699611 0.2435804  0.32661323  1.0871142  157.45400  1.049705 
D_true[19]  0.004794195 0.2307726 -0.32514359  0.3644341  200.30640  1.026036 
D_true[20] -0.635858657 0.3284276 -1.15963450 -0.1296354  121.93749  1.021399 
D_true[21] -0.840298873 0.2339080 -1.14985915 -0.4452980  132.71304  1.016325 
D_true[22] -1.593743931 0.2791787 -2.00661680 -1.1379765   88.35963  1.038408 
D_true[23] -0.508799117 0.2199743 -0.84137688 -0.1835776  163.77758  1.019792 
D_true[24] -0.591015865 0.2639985 -0.99833471 -0.1840042  164.34679  1.009509 
D_true[25]  0.016558541 0.2684125 -0.34323010  0.4590949   97.31032  1.046815 
D_true[26]  0.005612387 0.2280891 -0.34199392  0.3516850  175.94131  1.023529 
D_true[27] -0.026065003 0.3011759 -0.51793208  0.4063014  218.09803  1.015612 
D_true[28]  0.167931514 0.2394650 -0.22123523  0.5111828  238.11170  1.028259 
D_true[29] -0.667685858 0.2734505 -1.09641000 -0.2436933  229.47739  1.021412 
D_true[30] -1.364520590 0.2320125 -1.75454670 -1.0078210  249.65691  1.020234 
⋮          ⋮            ⋮         ⋮           ⋮           ⋮          ⋮        
M_true[25] -0.144055111 0.3536391 -0.70385405  0.43772495 1428.57634 1.0004901
M_true[26] -0.377158055 0.2102519 -0.71009822 -0.05402321  859.22827 1.0037720
M_true[27] -0.318411108 0.5244602 -1.20998695  0.47945743  592.91366 1.0017100
M_true[28] -0.136168440 0.3454640 -0.71255225  0.43302051 1011.96712 1.0030505
M_true[29] -0.716917356 0.4311626 -1.41365225 -0.02423970  574.15940 0.9993623
M_true[30] -1.375477236 0.1530471 -1.62525720 -1.12877620 1380.88793 1.0027514
M_true[31]  0.100921943 0.4255508 -0.59248861  0.76512839  650.46351 1.0127802
M_true[32] -0.864850837 0.1157580 -1.05238485 -0.67389823 1170.02918 1.0039008
M_true[33]  0.059856457 0.2474775 -0.31906609  0.48346513 1093.93463 0.9998014
M_true[34]  1.071215626 0.5997121  0.08319105  2.00382520  488.27059 1.0155316
M_true[35] -0.814176695 0.1536622 -1.06817280 -0.57351755  933.46042 1.0053110
M_true[36]  0.895438854 0.3086485  0.39800138  1.39885400  994.25593 1.0013404
M_true[37] -0.273296607 0.2713259 -0.69527436  0.16028290  941.94827 1.0039850
M_true[38] -1.192966414 0.1196683 -1.38808110 -1.00170680 1158.60937 0.9991165
M_true[39] -0.977587905 0.4977788 -1.82537870 -0.21564241  534.79346 1.0110180
M_true[40] -0.486332594 0.2914224 -0.93840780 -0.01272371  666.74721 1.0043081
M_true[41] -0.004221738 0.5396072 -0.88501488  0.87148547  787.52989 1.0052966
M_true[42] -0.158619489 0.2184540 -0.50894779  0.18794969  652.32065 1.0009158
M_true[43]  0.352639372 0.1572148  0.09860271  0.59983577  405.02189 1.0131173
M_true[44]  2.000767040 0.4219890  1.32774450  2.66212120 1060.05936 0.9991249
M_true[45] -0.711090357 0.5634777 -1.62351000  0.18241800  149.20356 1.0273174
M_true[46]  0.073085956 0.2137949 -0.25854329  0.41872518 1049.69761 1.0026678
M_true[47]  0.315990612 0.2511606 -0.09549357  0.71938581  714.52464 1.0031824
M_true[48]  0.438510841 0.4001224 -0.21584149  1.07845795 1010.47596 1.0073289
M_true[49] -0.740397593 0.1966482 -1.05315990 -0.42266740 1394.92812 1.0026860
M_true[50]  1.320617835 0.7040487  0.16763559  2.45795775  586.80804 1.0129785
a          -0.116986600 0.1047352 -0.28653033  0.04375597   53.92701 1.0845520
bA         -0.631187719 0.1639479 -0.88285200 -0.36273862  113.92094 1.0286757
bM          0.252313542 0.1986911 -0.05474765  0.57506156   96.64002 1.0121306
sigma       0.139227967 0.1126487  0.01123033  0.36524222   20.44014 1.2570123

Comparing the most important parameters:

iplot(function() {
  plot(
    coeftab(m.double.se.2, m15.2, m.double.se.1, m15.1), pars=c("a", "bA", "bM", "sigma"),
    main="All models"
  )
}, ar=1.7)

Notice we are getting some divergent transitions. It’s likely we’re running into these issues because there is more parameter space to explore; we’re less certain about all our inferences and so we need to cover more possibilities. Like vague priors, this can lead to divergent transitions.

The a, bA, and bM parameters have increased. With less certainty in the measurements, it’s more plausible that a stronger relationship holds (similar to what we saw in Figure 15.3). In contrast, the sigma parameter has decreased, because we can now explain more variation in terms of measurement error than general noise.

We also see, in general, more uncertainty in parameter estimates across the board (wider HDPI bars).

15M4. Simulate data from this DAG: \(X \rightarrow Y \rightarrow Z\). Now fit a model that predicts \(Y\) using both \(X\) and \(Z\). What kind of confound arises, in terms of inferring the causal influence of \(X\) on \(Y\)?

Answer. First lets simulate, with normal distributions:

N <- 1000
X <- rnorm(N)
Y <- rnorm(N, mean=X)
Z <- rnorm(N, mean=Y)

dat_list <- list(
  X = X,
  Y = Y,
  Z = Z
)

m_confound <- ulam(
  alist(
    Y ~ dnorm(mu, sigma),
    mu <- a + bX*X + bZ*Z,
    c(a, bX, bZ) ~ dnorm(0, 1),
    sigma ~ dexp(1)
  ), cmdstan = TRUE, data=dat_list, chains=4, cores=4
)

display(precis(m_confound, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m_confound, depth=3), main="m_confound")
}, ar=4.5)

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-1353560fc6.stan', line 21, column 4 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-1353560fc6.stan', line 21, column 4 to column 29)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 finished in 0.7 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 0.7 seconds.
Chain 3 finished in 0.7 seconds.
Chain 4 finished in 0.7 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.7 seconds.
Total execution time: 0.9 seconds.

      mean        sd         5.5%        94.5%      n_eff    Rhat4    
bZ    0.499327671 0.01588340  0.47408494 0.52493208 1671.529 1.0006997
bX    0.522257956 0.02570047  0.48127435 0.56275022 1614.158 1.0001688
a     0.001029397 0.02265268 -0.03482089 0.03796807 2213.710 0.9999208
sigma 0.704318764 0.01553018  0.67979880 0.72990238 2063.621 1.0007823

The confound we’re running into is post-treatment bias (conditioning on the outcome).

In this case \(Z\) is a better predictor of \(Y\) than \(X\) (notice the credible intervals). This is because the function between \(X\) and \(Y\) is the same as the function between \(Y\) and \(Z\), specifically the normal distribution with standard deviation equal one. Since \(Y\) will take on more dispersed values than \(X\), and \(Z\) more dispersed values than \(Y\), the more widely dispersed values of \(Z\) will predict \(Y\) better.

15M5. Return to the singing bird model, m15.9, and compare the posterior estimates of cat presence (PrC1) to the true simulated values. How good is the model at inferring the missing data? Can you think of a way to change the simulation so that the precision of the inference is stronger?

set.seed(9)
N_houses <- 100L
alpha <- 5
beta <- (-3)
k <- 0.5
r <- 0.2
cat <- rbern(N_houses, k)
notes <- rpois(N_houses, alpha + beta * cat)
R_C <- rbern(N_houses, r)
cat_obs <- cat
cat_obs[R_C == 1] <- (-9L)
dat <- list(
  notes = notes,
  cat = cat_obs,
  RC = R_C,
  N = as.integer(N_houses)
)

m15.9 <- ulam(
  alist(
    # singing bird model
    notes | RC == 0 ~ poisson(lambda),
    notes | RC == 1 ~ custom(log_sum_exp(
      log(k) + poisson_lpmf(notes | exp(a + b)),
      log(1 - k) + poisson_lpmf(notes | exp(a))
    )),
    log(lambda) <- a + b * cat,
    a ~ normal(5, 3),
    b ~ normal(-3, 2),

    # sneaking cat model
    cat | RC == 0 ~ bernoulli(k),
    k ~ beta(2, 2),

    # imputed values
    gq > vector[N]:PrC1 <- exp(lpC1) / (exp(lpC1) + exp(lpC0)),
    gq > vector[N]:lpC1 <- log(k) + poisson_lpmf(notes[i] | exp(a + b)),
    gq > vector[N]:lpC0 <- log(1 - k) + poisson_lpmf(notes[i] | exp(a))
  ),
  cmdstan = TRUE, data = dat, chains = 4, cores = 4
)
display(precis(m15.9, depth=3), mimetypes="text/plain")
iplot(function() {
  plot(precis(m15.9, depth=3), main="m15.9")
}, ar=0.2)

Warning in '/tmp/RtmpUakS9C/model-133612bcfe.stan', line 3, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

Warning in '/tmp/RtmpUakS9C/model-133612bcfe.stan', line 4, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

Warning in '/tmp/RtmpUakS9C/model-133612bcfe.stan', line 5, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.4 seconds.
Chain 2 finished in 0.4 seconds.
Chain 3 finished in 0.4 seconds.
Chain 4 finished in 0.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.4 seconds.
Total execution time: 0.5 seconds.

          mean        sd          5.5%         94.5%      n_eff    Rhat4    
a          1.6305739  0.06487360   1.5270624    1.7341393 1246.262 1.0007234
b         -0.8686414  0.12404909  -1.0655309   -0.6745582 1370.308 1.0002001
k          0.4613826  0.05293322   0.3792623    0.5478094 1578.813 0.9998864
lpC0[1]   -3.8678825  0.25926699  -4.2997987   -3.4731018 1307.226 1.0008298
lpC0[2]   -2.5370281  0.10940782  -2.7196871   -2.3717162 1604.795 1.0003821
lpC0[3]   -2.3969789  0.13123678  -2.6149412   -2.2036879 1472.008 1.0000513
lpC0[4]   -2.6412583  0.18020927  -2.9358888   -2.3658676 1357.050 1.0002515
lpC0[5]   -2.6412583  0.18020927  -2.9358888   -2.3658676 1357.050 1.0002515
lpC0[6]   -2.5370281  0.10940782  -2.7196871   -2.3717162 1604.795 1.0003821
lpC0[7]   -2.3758426  0.10192148  -2.5441358   -2.2288279 1564.462 0.9999698
lpC0[8]   -2.3969789  0.13123678  -2.6149412   -2.2036879 1472.008 1.0000513
lpC0[9]   -3.1732201  0.23694897  -3.5641388   -2.8074180 1306.062 1.0003836
lpC0[10]  -4.1106469  0.29703792  -4.6066845   -3.6474712 1282.459 1.0004647
lpC0[11]  -2.6412583  0.18020927  -2.9358888   -2.3658676 1357.050 1.0002515
lpC0[12]  -2.3969789  0.13123678  -2.6149412   -2.2036879 1472.008 1.0000513
lpC0[13]  -3.1732201  0.23694897  -3.5641388   -2.8074180 1306.062 1.0003836
lpC0[14]  -2.3758426  0.10192148  -2.5441358   -2.2288279 1564.462 0.9999698
lpC0[15]  -2.5370281  0.10940782  -2.7196871   -2.3717162 1604.795 1.0003821
lpC0[16]  -4.1106469  0.29703792  -4.6066845   -3.6474712 1282.459 1.0004647
lpC0[17]  -5.3072147  0.38228011  -5.9249241   -4.7161928 1275.016 1.0008237
lpC0[18]  -5.7412210  0.35879703  -6.3358815   -5.1745468 1270.034 1.0005170
lpC0[19]  -3.3012319  0.20096205  -3.6413410   -2.9936611 1348.218 1.0008078
lpC0[20]  -2.6412583  0.18020927  -2.9358888   -2.3658676 1357.050 1.0002515
lpC0[21]  -3.1732201  0.23694897  -3.5641388   -2.8074180 1306.062 1.0003836
lpC0[22]  -2.8523643  0.14822054  -3.0976313   -2.6325850 1441.347 1.0007080
lpC0[23]  -3.1732201  0.23694897  -3.5641388   -2.8074180 1306.062 1.0003836
lpC0[24]  -3.3012319  0.20096205  -3.6413410   -2.9936611 1348.218 1.0008078
lpC0[25]  -2.6412583  0.18020927  -2.9358888   -2.3658676 1357.050 1.0002515
lpC0[26]  -3.1732201  0.23694897  -3.5641388   -2.8074180 1306.062 1.0003836
lpC0[27]  -2.5370281  0.10940782  -2.7196871   -2.3717162 1604.795 1.0003821
⋮         ⋮           ⋮           ⋮            ⋮          ⋮        ⋮        
PrC1[71]  0.184080542 0.057826721 0.1004002350 0.28357880 1680.333 0.9999298
PrC1[72]  0.342496963 0.071722656 0.2327824950 0.45816177 1658.033 0.9997594
PrC1[73]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[74]  0.739902231 0.058983217 0.6397687550 0.82782831 1374.728 1.0002978
PrC1[75]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[76]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[77]  0.937878873 0.027476536 0.8864772500 0.97174892 1068.440 1.0006343
PrC1[78]  0.184080542 0.057826721 0.1004002350 0.28357880 1680.333 0.9999298
PrC1[79]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[80]  0.342496963 0.071722656 0.2327824950 0.45816177 1658.033 0.9997594
PrC1[81]  0.342496963 0.071722656 0.2327824950 0.45816177 1658.033 0.9997594
PrC1[82]  0.184080542 0.057826721 0.1004002350 0.28357880 1680.333 0.9999298
PrC1[83]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[84]  0.739902231 0.058983217 0.6397687550 0.82782831 1374.728 1.0002978
PrC1[85]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[86]  0.004186114 0.004196283 0.0006075278 0.01178349 1601.748 1.0012753
PrC1[87]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[88]  0.184080542 0.057826721 0.1004002350 0.28357880 1680.333 0.9999298
PrC1[89]  0.868392604 0.042487911 0.7907626650 0.92642500 1154.959 1.0005724
PrC1[90]  0.868392604 0.042487911 0.7907626650 0.92642500 1154.959 1.0005724
PrC1[91]  0.342496963 0.071722656 0.2327824950 0.45816177 1658.033 0.9997594
PrC1[92]  0.342496963 0.071722656 0.2327824950 0.45816177 1658.033 0.9997594
PrC1[93]  0.548759528 0.071358327 0.4348512500 0.65723994 1568.274 0.9998947
PrC1[94]  0.739902231 0.058983217 0.6397687550 0.82782831 1374.728 1.0002978
PrC1[95]  0.042253297 0.023729624 0.0141282785 0.08730985 1678.285 1.0005301
PrC1[96]  0.739902231 0.058983217 0.6397687550 0.82782831 1374.728 1.0002978
PrC1[97]  0.739902231 0.058983217 0.6397687550 0.82782831 1374.728 1.0002978
PrC1[98]  0.739902231 0.058983217 0.6397687550 0.82782831 1374.728 1.0002978
PrC1[99]  0.090158074 0.039075072 0.0389671790 0.16067187 1679.383 1.0002151
PrC1[100] 0.090158074 0.039075072 0.0389671790 0.16067187 1679.383 1.0002151

Indexes of the missing values:

display(which(R_C != 0))
post <- extract.samples(m15.9)

1. 1
2. 4
3. 14
4. 22
5. 26
6. 28
7. 29
8. 30
9. 32
10. 41
11. 46
12. 47
13. 48
14. 53
15. 55
16. 59
17. 60
18. 65
19. 66
20. 81
21. 83
22. 84
23. 92
24. 100

The mean values of PrC1 at those indexes:

mean_prc1 <- apply(post$PrC1, 2, mean)
display(mean_prc1[R_C != 0])

1. 0.00900593709
2. 0.5487595275
3. 0.1840805422
4. 0.04225329709
5. 0.7399022305
6. 0.3424969625
7. 0.8683926035
8. 0.9378788725
9. 0.5487595275
10. 0.7399022305
11. 0.3424969625
12. 0.8683926035
13. 0.3424969625
14. 0.3424969625
15. 0.04225329709
16. 0.7399022305
17. 0.5487595275
18. 0.7399022305
19. 0.01950867139
20. 0.3424969625
21. 0.5487595275
22. 0.7399022305
23. 0.3424969625
24. 0.0901580741

The true simulated values at those indexes:

display(cat[R_C != 0])

1. 0
2. 0
3. 0
4. 0
5. 1
6. 0
7. 1
8. 1
9. 1
10. 1
11. 0
12. 1
13. 1
14. 1
15. 0
16. 1
17. 0
18. 1
19. 0
20. 0
21. 0
22. 1
23. 0
24. 0

If we take PrC1 > 0.5 to mean a cat is present, then these missing elements were predicted correctly:

cat_prc1 <- as.integer(mean_prc1 > 0.5)
is_correct <- cat_prc1[R_C != 0] == cat[R_C != 0]
display(is_correct)

1. TRUE
2. FALSE
3. TRUE
4. TRUE
5. TRUE
6. TRUE
7. TRUE
8. TRUE
9. TRUE
10. TRUE
11. TRUE
12. TRUE
13. FALSE
14. FALSE
15. TRUE
16. TRUE
17. FALSE
18. TRUE
19. TRUE
20. TRUE
21. FALSE
22. TRUE
23. TRUE
24. TRUE

For a total accuracy of:

display(sum(is_correct) / length(is_correct))

0.791666666666667

15M6. Return to the four dog-eats-homework missing data examples. Simulate each and then fit one or more models to try to recover valid estimates for \(S \rightarrow H\).

Answer. A model fit to completely observed (co) data and a model fit to complete cases (cc) for the first scenario, where dogs eat homework at random:

library(data.table)

N <- 100
S <- rnorm(N)
H <- rbinom(N, size = 10, inv_logit(S))

D <- rbern(N) # dogs completely random
Hm <- H
Hm[D == 1] <- NA

d <- data.frame(S = S, H = H, Hm = Hm)
check_co <- function(d, name) {
  set.seed(27)
  dat.co <- list(
    S = d$S,
    H = d$H
  )
  m.dog.co <- ulam(
    alist(
      H ~ dbinom(10, p),
      logit(p) <- a + bS*S,
      a ~ dnorm(0, 1),
      bS ~ dnorm(0, 1)
    ), cmdstan = TRUE, data=dat.co, cores=4, chains=4
  )
  display(precis(m.dog.co, depth=2), mimetypes="text/plain")
  iplot(function() {
    plot(precis(m.dog.co, depth=2), main=name)
  }, ar=4.0)
}
check_co(d, "m.dog.1.co")

check_cc <- function(d, name) {
  # browser()
  set.seed(27)
  d_cc <- d[complete.cases(d$Hm), ]
  # t_d <- transpose(d_cc)
  # colnames(t_d) <- rownames(d_cc)
  # rownames(t_d) <- colnames(d_cc)
  # display(t_d)
  dat.cc <- list(
    S = d_cc$S,
    Hm = d_cc$Hm
  )
  display_markdown("Sample of data remaining for 'complete case' analysis: ")
  display(head(d_cc))
  m.dog.cc <- ulam(
    alist(
      Hm ~ dbinom(10, p),
      logit(p) <- a + bS*S,
      a ~ dnorm(0, 1),
      bS ~ dnorm(0, 1)
    ), cmdstan = TRUE, data=dat.cc, cores=4, chains=4
  )
  display(precis(m.dog.cc, depth=2), mimetypes="text/plain")
  iplot(function() {
    plot(precis(m.dog.cc, depth=2), main=name)
  }, ar=4.0)
}
check_cc(d, "m.dog.1.cc")

Warning in '/tmp/RtmpUakS9C/model-13625692d9.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean        sd         5.5%       94.5%     n_eff    Rhat4    
a  -0.00386812 0.07272701 -0.1179536 0.1137859 1440.388 0.9994012
bS  1.01038932 0.08952928  0.8723994 1.1575217 1035.581 1.0002041

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

Sample of data remaining for 'complete case' analysis:

A data.frame: 6 × 3

	S	H	Hm
	<dbl>	<int>	<int>
3	-0.3483160	6	6
5	0.2841523	5	5
8	0.3812747	5	5
10	0.3922908	8	8
13	0.1487904	8	8
16	1.1348391	9	9

Warning in '/tmp/RtmpUakS9C/model-134311df7a.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean        sd         5.5%       94.5%     n_eff    Rhat4   
a  -0.03654073 0.09791599 -0.1900768 0.1254365 1339.579 1.001207
bS  0.92498400 0.12369701  0.7329099 1.1265327 1116.149 1.000748

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

As expected (explained in the chapter) we are able to infer bS both with and without all the data.

In the second scenario, dogs of hard-working students eat homework. We don’t need to fit a new completely observed (co) model because it shouldn’t be any different from m.dog.1.co. A model fit to complete cases (cc):

D2 <- ifelse(S > 0, 1, 0)
Hm2 <- H
Hm2[D2 == 1] <- NA
d <- data.frame(S = S, H = H, Hm = Hm2)
check_cc(d, "m.dog.2.cc")

Sample of data remaining for 'complete case' analysis:

A data.frame: 6 × 3

	S	H	Hm
	<dbl>	<int>	<int>
2	-1.29318475	2	2
3	-0.34831602	6	6
4	-0.35097548	2	2
20	-1.44522679	2	2
21	-0.68131726	2	2
24	-0.04805382	6	6

Warning in '/tmp/RtmpUakS9C/model-13719a49e2.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean       sd        5.5%       94.5%      n_eff    Rhat4   
a  -0.2211787 0.1625756 -0.4830884 0.03894572 567.4133 1.007038
bS  0.8135978 0.2028660  0.4879648 1.14727990 421.4382 1.008410

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

As expected (explained in the chapter) we are able to infer bS both with and without all the data.

The text suggests a variation on scenario 2 where the function \(S \rightarrow H\) is nonlinear and unobserved only in the domain of the function where it is non-linear. Let’s simulate this scenario and attempt to infer the coefficient of the (wrong) linear function with and without all the data:

H3 <- ifelse(S > 0, H, 0)
Hm3 <- H3
Hm3[D2 == 1] <- NA
d <- data.frame(S = S, H = H3, Hm = Hm3)
check_co(d, "m.dog.2.var.co")
check_cc(d, "m.dog.2.var.cc")

Warning in '/tmp/RtmpUakS9C/model-1325283428.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean      sd         5.5%      94.5%      n_eff    Rhat4   
a  -1.065430 0.09883981 -1.225324 -0.9143453 983.0190 1.001233
bS  2.289743 0.14768992  2.061293  2.5359345 849.6928 1.001852

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

Sample of data remaining for 'complete case' analysis:

A data.frame: 6 × 3

	S	H	Hm
	<dbl>	<dbl>	<dbl>
2	-1.29318475	0	0
3	-0.34831602	0	0
4	-0.35097548	0	0
20	-1.44522679	0	0
21	-0.68131726	0	0
24	-0.04805382	0	0

Warning in '/tmp/RtmpUakS9C/model-1342371bb4.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean      sd        5.5%       94.5%     n_eff    Rhat4   
a  -3.937813 0.4984334 -4.7425330 -3.171735 937.6243 1.003867
bS  1.580369 0.7318090  0.4746245  2.783566 943.3718 1.004776

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

As expected (explained in the chapter) both these inferences are wrong. The cc inference is wrong because despite complete observability our model was linear when the function was non-linear. The co inference is extra wrong because our sampling is also biased. The data going into this second model is completely uninformative; all the S are random negative numbers and the H are 0. Notice that we confusingly still get a positive inference for bS. The model wants to push everything in the log link as negative as possible, but a is already as negative as it can be given our thin tailed (normal) prior on it. We can push everything farther negative by making bS positive, since all the S values are negative. Be skeptical of models that have any parameter inferences that are extreme relative to the parameter’s priors (like a)! These mean that either your prior is wrong, your data, or your model.

The third scenario, where some unobserved variable like the noise level in the student’s house affects both the homework quality and whether dogs eat homework, is already partially covered in the chapter. Let’s consider the variation suggested in R code 15.14, however:

N2 <- 1000
X <- rnorm(N2)
S2 <- rnorm(N2)
H4 <- rbinom(N2, size = 10, inv_logit(2 + S2 - 2 * X))
# D <- ifelse(X > 1, 1, 0)
D <- ifelse(abs(X) < 1, 1, 0)
Hm4 <- H4
Hm4[D == 1] <- NA
d <- data.frame(S = S2, H = H4, Hm = Hm4)
check_co(d, "m.dog.3.var.co")
check_cc(d, "m.dog.3.var.cc")

Warning in '/tmp/RtmpUakS9C/model-1357ace0f7.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean      sd         5.5%      94.5%     n_eff     Rhat4    
a  1.2221746 0.02611138 1.1807394 1.2637615 1019.2430 0.9993861
bS 0.6320652 0.02500255 0.5934007 0.6724795  944.8489 1.0001442

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 finished in 0.7 seconds.
Chain 3 finished in 0.8 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 0.8 seconds.
Chain 4 finished in 0.8 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.8 seconds.
Total execution time: 0.9 seconds.

Sample of data remaining for 'complete case' analysis:

A data.frame: 6 × 3

	S	H	Hm
	<dbl>	<int>	<int>
1	0.91583404	3	3
3	-0.75727948	9	9
8	0.06453134	10	10
9	1.22364898	10	10
11	-0.85843194	10	10
14	0.12535586	10	10

Warning in '/tmp/RtmpUakS9C/model-137399b0b3.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean      sd         5.5%      94.5%     n_eff     Rhat4   
a  0.6269566 0.03576682 0.5709384 0.6832952 1348.5960 1.000598
bS 0.3857848 0.03922876 0.3241633 0.4471760  732.3931 1.005366

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 finished in 0.2 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 0.3 seconds.
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 0.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.4 seconds.

As expected (see the text) missingness makes the inference for bS worse than with complete observability. We’ve changed the missingness function to remove small values rather than large values. When we remove large values, the remaining observations are farther from a ceiling and therefore the influence of predictors seems larger, which is what we want when we have an unobserved variable. When we remove small values, the remaining observations are closer to a ceiling or a floor and the influence of predictors seems smaller.

In the fourth scenario dogs eat bad homework. Little can be done to get a good inference:

S.4 <- rnorm(N)
H.4 <- rbinom(N, size = 10, inv_logit(S.4))
D.4 <- ifelse(H.4 < 5, 1, 0)
Hm.4 <- H.4
Hm.4[D.4 == 1] <- NA
d <- data.frame(S = S.4, H = H.4, Hm = Hm.4)
check_cc(d, "m.dog.4.cc")

Sample of data remaining for 'complete case' analysis:

A data.frame: 6 × 3

	S	H	Hm
	<dbl>	<int>	<int>
1	1.67531074	8	8
2	0.80350151	7	7
5	0.25107932	7	7
7	0.35701984	6	6
12	0.05818869	5	5
16	0.18096825	5	5

Warning in '/tmp/RtmpUakS9C/model-132a73d6c8.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

   mean      sd        5.5%      94.5%     n_eff    Rhat4    
a  0.3559647 0.1110200 0.1814343 0.5382437 702.4761 0.9993733
bS 0.6936059 0.1352675 0.4730826 0.9099878 747.4805 1.0009603

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

15H1. The data in data(elephants) are counts of matings observed for bull elephants of differing ages. There is a strong positive relationship between age and matings. However, age is not always assessed accurately. First, fit a Poisson model predicting MATINGS with AGE as a predictor. Second, assume that the observed AGE values are uncertain and have a standard error of ±5 years. Re-estimate the relationship between MATINGS and AGE, incorporating this measurement error. Compare the inferences of the two models.

Answer. The full elephants data.frame:

data(elephants)
display(elephants)

q15.h1.a <- function(d) {
  dat <- list(
    M = as.integer(d$MATINGS),
    A = as.numeric(standardize(d$AGE))
  )
  m15.h1.a <- ulam(
    alist(
      M ~ dpois(lambda),
      lambda <- exp(a + bA*A),
      a ~ dnorm(0, 1),
      bA ~ dnorm(0, 1)
    ), cmdstan = TRUE, data=dat, cores = 4, chains = 4, log_lik = TRUE
  )
  display_precis(m15.h1.a, "m15.h1.a", 4.0)
  return(m15.h1.a)
}
q15.h1.b <- function(d, error=5.0) {
  set.seed(31)
  stdA = standardize(d$AGE)
  dat <- list(
    M = as.integer(d$MATINGS),
    A_obs = as.numeric(stdA),
    A_sd = error / attr(stdA, "scaled:scale"),
    N = nrow(d)
  )
  m15.h1.b <- ulam(
    alist(
      M ~ dpois(lambda),
      A_obs ~ dnorm(A_true, A_sd),
      lambda <- exp(a + bA*A_true[i]),
      vector[N]:A_true ~ dnorm(0, 1),
      a ~ dnorm(0, 1),
      bA ~ dnorm(0, 1)
    ), cmdstan = TRUE, data=dat, cores = 4, chains = 4, log_lik = TRUE
  )
  display_precis(m15.h1.b, "m15.h1.b", 0.8)
  return(m15.h1.b)
}

# m15.h1.a <- q15.h1.a(elephants)
# m15.h1.b <- q15.h1.b(elephants)
#
# iplot(function() {
#   plot(coeftab(m15.h1.b, m15.h1.a), pars=c("a", "bA"), main="Parameter comparison for elephant models")
# }, ar=3.0)

A data.frame: 41 × 2

AGE	MATINGS
<int>	<int>
27	0
28	1
28	1
28	1
28	3
29	0
29	0
29	0
29	2
29	2
29	2
30	1
32	2
33	4
33	3
33	3
33	3
33	2
34	1
34	1
34	2
34	3
36	5
36	6
37	1
37	1
37	6
38	2
39	1
41	3
42	4
43	0
43	2
43	3
43	4
43	9
44	3
45	5
47	7
48	2
52	9

The second model infers both a larger and more uncertain estimate (wider credible intervals) for bA.

15H2. Repeat the model fitting problem above, now increasing the assumed standard error on AGE. How large does the standard error have to get before the posterior mean for the coefficient on AGE reaches zero?

Answer. The coefficient’s uncertainty gets large enough to cover zero when we get to roughly 20-25 years of age uncertainty:

q15.h1.b(elephants, error=25)

Warning in '/tmp/RtmpUakS9C/model-135a0b78b7.stan', line 3, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.4 seconds.

Raw data (following plot):

           mean        sd        5.5%       94.5%     n_eff       Rhat4   
A_true[1]  -0.60117469 1.0091844 -2.0722995 1.0873623   14.344321 1.120206
A_true[2]  -0.36787650 0.8513137 -1.7027519 0.9999452   28.260337 1.056953
A_true[3]  -0.36570746 0.8614791 -1.7286127 0.9951490   33.241592 1.057124
A_true[4]  -0.37982283 0.8580098 -1.7140640 0.9987858   35.989888 1.050788
A_true[5]   0.07493666 0.6996594 -1.0779699 1.1663960  595.483364 1.013497
A_true[6]  -0.60846831 1.0250175 -2.1094692 1.1398493   11.788646 1.131486
A_true[7]  -0.61409542 1.0013380 -2.0913275 1.0568451   12.110301 1.132563
A_true[8]  -0.62425165 1.0104447 -2.2182763 1.1239754   11.248228 1.138278
A_true[9]  -0.13277755 0.7341117 -1.3427536 1.0298548 1617.384900 1.006307
A_true[10] -0.13387963 0.7824652 -1.4082887 1.0681719 1429.816665 1.004412
A_true[11] -0.14135320 0.7575667 -1.3215731 1.0423720 1685.314777 1.007532
A_true[12] -0.36961043 0.8986039 -1.8168895 1.0650409   37.249962 1.053382
A_true[13] -0.11639938 0.8248625 -1.4488624 1.1664826 2303.501889 1.005051
A_true[14]  0.28239312 0.7969648 -1.0607951 1.4582005   37.533339 1.057184
A_true[15]  0.07194796 0.7655739 -1.1629314 1.2292131 1136.502776 1.009605
A_true[16]  0.09594969 0.7445255 -1.1186045 1.2374490 1617.463497 1.003458
A_true[17]  0.09693792 0.7447062 -1.0956143 1.2555275 1004.926624 1.009668
A_true[18] -0.09630907 0.7648403 -1.3584912 1.1129924 1837.140523 1.005444
A_true[19] -0.32319797 0.8443859 -1.6560516 1.0874029   27.516203 1.061404
A_true[20] -0.32582686 0.8627404 -1.7263214 1.0351113   63.919385 1.037648
A_true[21] -0.08356380 0.7170041 -1.2493470 1.0399564 3002.595613 1.001265
A_true[22]  0.10043619 0.7392507 -1.1157200 1.2598396 1280.651269 1.006613
A_true[23]  0.44967315 0.8409441 -1.0400442 1.6324739   11.980876 1.146094
A_true[24]  0.62802098 0.9662659 -1.2159055 1.9125331    9.630480 1.190744
A_true[25] -0.29544474 0.8628739 -1.6330445 1.0469239   46.111393 1.048206
A_true[26] -0.29177282 0.8602039 -1.6205196 1.1353206   49.601660 1.047328
A_true[27]  0.67260803 0.9505727 -1.1801111 1.9040425    9.303567 1.193147
A_true[28] -0.08947878 0.7602801 -1.3539159 1.0541851 2115.792679 1.003276
A_true[29] -0.30100218 0.8948009 -1.7469175 1.1253371   33.215895 1.057725
A_true[30]  0.13275653 0.7520887 -1.1455423 1.2932678 1170.307847 1.008382
A_true[31]  0.34172283 0.7643616 -0.9776645 1.5050318   31.367545 1.056077
A_true[32] -0.48979044 0.9963426 -2.0226607 1.2049862   12.988025 1.126033
A_true[33] -0.05266336 0.7399035 -1.2451754 1.0873137 2050.562082 1.003732
A_true[34]  0.14068062 0.7345973 -1.0965141 1.3081742  642.349157 1.009799
A_true[35]  0.33953829 0.7762237 -0.9566545 1.5095188   29.663239 1.058023
A_true[36]  1.09639730 1.2996878 -1.7307875 2.5033904    5.481716 1.371427
A_true[37]  0.17087846 0.7288982 -1.0229196 1.2856901 1310.801033 1.007596
A_true[38]  0.53593754 0.8557067 -0.9830917 1.7093324   10.907019 1.148637
A_true[39]  0.83153739 1.0640228 -1.3882139 2.1190046    6.843605 1.274743
A_true[40] -0.02768258 0.7587293 -1.2566854 1.1798185 1898.656795 1.004493
A_true[41]  1.12761587 1.2898501 -1.7181836 2.5366333    5.789249 1.355103
a           0.81522363 0.1435477  0.5869336 1.0458022 1203.674478 1.001330
bA          0.39535337 0.4501382 -0.6236069 0.8170445    4.838490 1.464206

Hamiltonian Monte Carlo approximation
2000 samples from 4 chains

Sampling durations (seconds):
        warmup sample total
chain:1   0.18   0.10  0.28
chain:2   0.18   0.09  0.27
chain:3   0.18   0.11  0.29
chain:4   0.17   0.10  0.28

Formula:
M ~ dpois(lambda)
A_obs ~ dnorm(A_true, A_sd)
lambda <- exp(a + bA * A_true[i])
vector[N]:A_true ~ dnorm(0, 1)
a ~ dnorm(0, 1)
bA ~ dnorm(0, 1)

15H3. The fact that information flows in all directions among parameters sometimes leads to rather unintuitive conclusions. Here’s an example from missing data imputation, in which imputation of a single datum reverses the direction of an inferred relationship. Use these data:

R code 15.32

set.seed(100)
x <- c( rnorm(10) , NA )
y <- c( rnorm(10,x) , 100 )
d <- list(x=x,y=y)

These data comprise 11 cases, one of which has a missing predictor value. You can quickly confirm that a regression of y on x for only the complete cases indicates a strong positive relationship between the two variables. But now fit this model, imputing the one missing value for x:

\[\begin{split} \begin{align} y_i & \sim Normal(\mu_i, \sigma) \\ \mu_i & = \alpha + \beta x_i \\ x_i & \sim Normal(0, 1) \\ \alpha & \sim Normal(0, 100) \\ \beta & \sim Normal(0, 100) \\ \sigma & \sim Normal(1) \end{align} \end{split}\]

Be sure to run multiple chains. What happens to the posterior distribution of \(\beta\)? Be sure to inspect the full density. Can you explain the change in inference?

ERROR: The last line of the model should be: $\( \sigma \sim Exponential(1) \)$

Answer. The list d produced by R code 15.32:

r.code.15.32 <- function() {
  set.seed(100)
  x <- c(rnorm(10), NA)
  y <- c(rnorm(10, x), 100)
  d <- list(x = x, y = y)
  display(d)
  return(data.frame(d))
}
d.15.h3 <- r.code.15.32()

$x

1. -0.502192350531457
2. 0.131531165327303
3. -0.07891708981887
4. 0.886784809417845
5. 0.116971270510841
6. 0.318630087617032
7. -0.58179068471591
8. 0.714532710891568
9. -0.825259425862769
10. -0.359862131395465
11. <NA>

$y

1. -0.412306206753926
2. 0.227805625612433
3. -0.280551042002224
4. 1.62662530929628
5. 0.240350771599711
6. 0.289313378387685
7. -0.97064493161941
8. 1.22538896826149
9. -1.73907361123195
10. 1.95043469138359
11. 100

Fitting a model to complete cases, using the same model in the question minus the \(x_i\) prior:

q.15.h3.a <- function(d) {
  d_cc <- d[complete.cases(d$x),]
  dat <- list(
    Y = d_cc$y,
    X = d_cc$x
  )
  m15.h3.a <- ulam(
    alist(
      Y ~ dnorm(mu, sigma),
      mu <- a + b*X,
      c(a, b) ~ dnorm(0, 1),
      sigma ~ dnorm(0, 1)
    ), cmdstan = TRUE, data=dat, cores=4, chains=4, log_lik=TRUE
  )
  display_precis(m15.h3.a, "m15.h3.a", 4.5)
  return(m15.h3.a)
}
m15.h3.a <- q.15.h3.a(d.15.h3)

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Scale parameter is -4.90991, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 2 Rejecting initial value:

Chain 2   Error evaluating the log probability at the initial value.

Chain 2 Exception: normal_lpdf: Scale parameter is -1.99224, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Exception: normal_lpdf: Scale parameter is -1.99224, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Rejecting initial value:

Chain 2   Error evaluating the log probability at the initial value.

Chain 2 Exception: normal_lpdf: Scale parameter is -0.908125, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Exception: normal_lpdf: Scale parameter is -0.908125, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Rejecting initial value:

Chain 2   Error evaluating the log probability at the initial value.

Chain 2 Exception: normal_lpdf: Scale parameter is -0.850954, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Exception: normal_lpdf: Scale parameter is -0.850954, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Rejecting initial value:

Chain 2   Error evaluating the log probability at the initial value.

Chain 2 Exception: normal_lpdf: Scale parameter is -0.806948, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Exception: normal_lpdf: Scale parameter is -0.806948, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Rejecting initial value:

Chain 2   Error evaluating the log probability at the initial value.

Chain 2 Exception: normal_lpdf: Scale parameter is -0.966028, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Exception: normal_lpdf: Scale parameter is -0.966028, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -1.72548, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -0.637511, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -19.6802, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -1.36893, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -0.58501, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -3.18212, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is -3.81429, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Rejecting initial value:

Chain 3   Error evaluating the log probability at the initial value.

Chain 3 Exception: normal_lpdf: Scale parameter is -1.38994, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Exception: normal_lpdf: Scale parameter is -1.38994, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Rejecting initial value:

Chain 3   Error evaluating the log probability at the initial value.

Chain 3 Exception: normal_lpdf: Scale parameter is -0.394764, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Exception: normal_lpdf: Scale parameter is -0.394764, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Rejecting initial value:

Chain 3   Error evaluating the log probability at the initial value.

Chain 3 Exception: normal_lpdf: Scale parameter is -1.13481, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Exception: normal_lpdf: Scale parameter is -1.13481, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Rejecting initial value:

Chain 3   Error evaluating the log probability at the initial value.

Chain 3 Exception: normal_lpdf: Scale parameter is -0.175318, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Exception: normal_lpdf: Scale parameter is -0.175318, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: normal_lpdf: Scale parameter is -0.333645, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: normal_lpdf: Scale parameter is -0.819312, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: normal_lpdf: Scale parameter is -1.57384, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: normal_lpdf: Scale parameter is -0.695607, but must be positive! (in '/tmp/RtmpUakS9C/model-134c4692d5.stan', line 18, column 4 to column 29)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.0 seconds.
Chain 4 finished in 0.0 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.4 seconds.

Raw data (following plot):

      mean      sd        5.5%       94.5%     n_eff     Rhat4    
b     1.0796106 0.4921676  0.2869450 1.8460359 1080.7259 1.0021330
a     0.2105812 0.2997821 -0.2759974 0.6750308 1180.1113 0.9997719
sigma 0.9642299 0.2458819  0.6491778 1.3962750  877.1299 1.0038247

Fitting a model to all the data, with the imputation model suggested in the question:

q.15.h3.b <- function(d) {
  dat <- list(
    Y = d$y,
    X = d$x
  )
  m15.h3.b <- ulam(
    alist(
      Y ~ dnorm(mu, sigma),
      mu <- a + b*X,
      X ~ dnorm(0, 1),
      c(a, b) ~ dnorm(0, 1),
      sigma ~ dexp(1)
    ), cmdstan = TRUE, data=dat, cores=4, chains=4, log_lik=TRUE
  )
  display_precis(m15.h3.b, "m15.h3.b", 4.0)
  return(m15.h3.b)
}
m15.h3.b <- q.15.h3.b(d.15.h3)

Found 1 NA values in X and attempting imputation.

Warning in '/tmp/RtmpUakS9C/model-132c5f8f42.stan', line 4, column 26: Declaration

    of arrays by placing brackets after a type is deprecated and will be

    removed in Stan 2.32.0. Instead use the array keyword before the type.

    This can be changed automatically using the auto-format flag to stanc

Warning in '/tmp/RtmpUakS9C/model-132c5f8f42.stan', line 17, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.2 seconds.

Raw data (following plot):

            mean        sd        5.5%      94.5%     n_eff    Rhat4    
b            0.02013376 1.0382012 -1.684889  1.686190 2102.224 0.9996513
a            0.29524993 0.9979384 -1.314590  1.880945 2317.508 0.9989613
sigma       18.87964020 2.2106872 15.649452 22.568184 2517.860 1.0004584
X_impute[1] -0.01639253 1.0044745 -1.625428  1.616405 2165.268 0.9994991

Notice the model has imputed a value near zero for the single missing value. This point is extremely influential because of the large value of y associated with it (the 100) so it dominates the inference results. A value of zero associated with 100 implies almost no relationship between x and y so the model infers that sigma is large instead.

15H4. Using data(Primates301), consider the relationship between brain volume (brain) and body mass (body). These variables are presented as single values for each species. However, there is always a range of sizes in a species, and some of these measurements are taken from very small samples. So these values are measured with some unknown error.

We don’t have the raw measurements to work with - that would be best. But we can imagine what might happen if we had them. Suppose error is proportional to the measurement. This makes sense, because larger animals have larger variation. As a consequence, the uncertainty is not uniform across the values and this could mean trouble.

Let’s make up some standard errors for these measurements, to see what might happen. Load the data and scale the the measurements so the maximum is 1 in both cases:

R code 15.33

library(rethinking)
data(Primates301)
d <- Primates301
cc <- complete.cases( d$brain , d$body )
B <- d$brain[cc]
M <- d$body[cc]
B <- B / max(B)
M <- M / max(M)

Now I’ll make up some standard errors for \(B\) and \(M\), assuming error is 10% of the measurement.

R code 15.34

Bse <- B*0.1
Mse <- M*0.1

Let’s model these variables with this relationship:

\[\begin{split} \begin{align} B_i & \sim LogNormal(\mu_i, \sigma) \\ \mu_i & = \alpha + \beta log M_i \end{align} \end{split}\]

This says that brain volume is a log-normal variable, and the mean on the log scale is given by \(\mu\). What this model implies is that the expected value of B is:

\[ E(B_i | M_i) = exp(\alpha)M_i^{\beta} \]

So this is a standard allometric scaling relationship - incredibly common in biology.

Ignoring measurement error, the corresponding ulam model is:

R code 15.35

dat_list <- list( B = B , M = M )
m15H4 <- ulam(
 alist(
  B ~ dlnorm( mu , sigma ),
  mu <- a + b*log(M),
  a ~ normal(0,1),
  b ~ normal(0,1),
  sigma ~ exponential(1)
 ) , cmdstan = TRUE, data=dat_list )

Your job is to add the measurement errors to this model. Use the divorce/marriage example in the chapter as a guide. It might help to initialize the unobserved true values of \(B\) and \(M\) using the observed values, by adding a list like this to ulam:

R code 15.36

start=list( M_true=dat_list$M , B_true=dat_list$B )

Compare the inference of the measurement error model to those of m1.1 above. Has anything changed? Why or why not?

ERROR:

Compare the inference of the measurement error model to those of m1.1 above.

The author meant m15H4.

Solution. The math is a bit fuzzy for the allometric relationship. The mean of a LogNormal random variable is:

\[ E[X] = exp(\mu + \sigma^2/2) \]

The author seems to be ignoring the \(\sigma^2/2\) term. If you accept that ‘simplification’ you can get the second equation he presents with:

\[\begin{split} \begin{align} E(B_i|M_i) & = exp(\mu_i) \\ & = exp(\alpha + \beta \cdot log(M_i)) \\ & = e^{\alpha} e^{log(M_i) \cdot \beta} \\ & = e^{\alpha} M_i^{\beta} \end{align} \end{split}\]

load.d.15.h4 <- function() {
  data(Primates301)
  d <- Primates301
  cc <- complete.cases( d$brain , d$body )
  B <- d$brain[cc]
  M <- d$body[cc]
  B <- B / max(B)
  M <- M / max(M)

  Bse <- B*0.1
  Mse <- M*0.1

  d.15.h4 <- data.frame(B=B, M=M, Bse=Bse, Mse=Mse)
  return(d.15.h4)
}
d.15.h4 <- load.d.15.h4()

q.15.h4.a <- function(d) {
  dat_list <- list(B = d$B, M = d$M)
  m15H4 <- ulam(
    alist(
      B ~ dlnorm(mu, sigma),
      mu <- a + b * log(M),
      a ~ normal(0, 1),
      b ~ normal(0, 1),
      sigma ~ exponential(1)
    ), cmdstan = TRUE, data=dat_list, chains=4, cores=4
  )
  display_precis(m15H4, "m15H4", 4.0)
  return(m15H4)
}
m15H4 <- q.15.h4.a(d.15.h4)

q.15.h4.b <- function(d) {
  dat_list <- list(
    B_obs = d$B,
    M_obs = d$M,
    B_sd = d$Bse,
    M_sd = d$Mse,
    N_obs = nrow(d)
  )
  m.15.h4.b <- ulam(
    alist(
      B_obs ~ dnorm(B_true, B_sd),
      vector[N_obs]:B_true ~ dlnorm(mu, sigma),
      mu <- a + b * log(M_true[i]),
      M_obs ~ dnorm(M_true, M_sd),
      vector[N_obs]:M_true ~ dexp(1),
      a ~ normal(0, 1),
      b ~ normal(0, 1),
      sigma ~ exponential(1)
    ), cmdstan = TRUE, data=dat_list, chains=4, cores=4, log_lik=TRUE
  )
  display_precis(m.15.h4.b, "m.15.h4.b", 0.2)
  return(m.15.h4.b)
}
m.15.h4.b <- q.15.h4.b(d.15.h4)

# iplot_coeftab <- function(..., pars, ar, main="coeftab") {
#   iplot(function() {
#     plot(coeftab(...), pars=pars, main=main)
#   }, ar=ar)
# }
# iplot_coeftab(m15H4, m.15.h4.b, pars=c("a", "b", "sigma"), ar=2.0)
iplot(function() {
  plot(coeftab(m15H4, m.15.h4.b), pars=c("a", "b", "sigma"), main="coeftab")
}, ar=2)

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Scale parameter is 0, but must be positive finite! (in '/tmp/RtmpUakS9C/model-1371183445.stan', line 18, column 4 to column 32)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Scale parameter is 0, but must be positive finite! (in '/tmp/RtmpUakS9C/model-1371183445.stan', line 18, column 4 to column 32)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Scale parameter is inf, but must be positive finite! (in '/tmp/RtmpUakS9C/model-1371183445.stan', line 18, column 4 to column 32)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.5 seconds.

Raw data (following plot):

      mean      sd         5.5%      94.5%     n_eff     Rhat4   
a     0.4264338 0.05988520 0.3312353 0.5224070  700.3517 1.003542
b     0.7830678 0.01441620 0.7598551 0.8058107  711.1302 1.003888
sigma 0.2936612 0.01548044 0.2704216 0.3199446 1012.5430 1.000912

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[2] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[9] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[9] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[9] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[25] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[74] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[79] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Location parameter[81] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 21, column 4 to column 36)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[2] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Location parameter[25] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[10] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[10] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[10] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[10] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[27] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[27] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Location parameter[81] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 21, column 4 to column 36)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Location parameter[10] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[1] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[5] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[8] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[11] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[27] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[79] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[79] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: normal_lpdf: Location parameter[81] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 21, column 4 to column 36)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Location parameter[5] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[3] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[3] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[4] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[4] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[9] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[9] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[13] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[13] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[54] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: normal_lpdf: Location parameter[81] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 21, column 4 to column 36)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[1] is -inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[3] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Location parameter[13] is inf, but must be finite! (in '/tmp/RtmpUakS9C/model-1311d775da.stan', line 25, column 4 to column 37)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 finished in 2.1 seconds.
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 2.2 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 finished in 2.4 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 2.6 seconds.

All 4 chains finished successfully.
Mean chain execution time: 2.3 seconds.
Total execution time: 2.7 seconds.

Raw data (following plot):

            mean         sd           5.5%         94.5%       n_eff   
B_true[1]   0.116501049  0.0108402991 0.098894523  0.133888925 4151.292
B_true[2]   0.110231006  0.0101802896 0.094368209  0.126886765 4166.197
B_true[3]   0.108207786  0.0102709931 0.092084153  0.124479100 5120.396
B_true[4]   0.106226831  0.0095893161 0.090907526  0.121450265 4162.346
B_true[5]   0.104855479  0.0098984460 0.089063354  0.121058310 4432.129
B_true[6]   0.110209395  0.0094850710 0.094820756  0.125500275 4889.156
B_true[7]   0.122010135  0.0120237453 0.102857645  0.141140965 4963.691
B_true[8]   0.114012183  0.0100984946 0.097958978  0.130092905 5001.726
B_true[9]   0.041210848  0.0039537034 0.034655030  0.047464372 4923.059
B_true[10]  0.031967147  0.0030520460 0.027005560  0.036873166 4120.390
B_true[11]  0.033990621  0.0032671643 0.028946007  0.039204179 5147.285
B_true[12]  0.011615569  0.0011530438 0.009795893  0.013472539 3759.581
B_true[13]  0.013879664  0.0013195969 0.011810195  0.016042645 5508.006
B_true[14]  0.226184862  0.0229068566 0.189930505  0.262839200 5025.154
B_true[15]  0.224777534  0.0220302396 0.189755270  0.259387860 4569.248
B_true[16]  0.205150567  0.0197638475 0.174193665  0.236584235 4188.408
B_true[17]  0.205278207  0.0203633079 0.171790135  0.238753335 5450.778
B_true[18]  0.021432650  0.0018572777 0.018437834  0.024353876 4434.788
B_true[19]  0.017127153  0.0015081086 0.014662933  0.019543994 4186.990
B_true[20]  0.211951499  0.0227170391 0.174935600  0.247565000 5339.681
B_true[21]  0.140854892  0.0151150537 0.116608590  0.164992495 4918.582
B_true[22]  0.128264715  0.0137268858 0.106459525  0.150669660 4512.236
B_true[23]  0.022477453  0.0021018929 0.019144861  0.025896720 4993.448
B_true[24]  0.015836114  0.0014732390 0.013535502  0.018186942 5526.985
B_true[25]  0.014493791  0.0014178371 0.012279597  0.016777236 6602.060
B_true[26]  0.014668672  0.0013989743 0.012521752  0.016804271 4540.035
B_true[27]  0.008086627  0.0007714216 0.006856515  0.009321702 6237.143
B_true[28]  0.121937390  0.0125656239 0.102376120  0.141602715 4603.340
B_true[29]  0.124567333  0.0134158195 0.103564470  0.146810950 4432.053
B_true[30]  0.134406476  0.0142944320 0.111640525  0.157261400 3971.116
⋮           ⋮            ⋮            ⋮            ⋮           ⋮       
M_true[156] 0.0218681341 2.198526e-03 0.0184078130 0.025338225 4455.942
M_true[157] 0.0740924054 7.433089e-03 0.0621375165 0.086132004 4547.262
M_true[158] 0.1119229651 1.180245e-02 0.0931663405 0.130695530 5311.109
M_true[159] 0.0030852139 2.975357e-04 0.0026144887 0.003569972 4734.040
M_true[160] 0.0039918902 3.786749e-04 0.0033638231 0.004593379 4221.715
M_true[161] 0.0040281731 4.056766e-04 0.0033881246 0.004692708 5922.510
M_true[162] 0.0043007993 4.155236e-04 0.0036304378 0.004958710 3708.773
M_true[163] 0.0045064638 4.245070e-04 0.0038333670 0.005177733 4897.942
M_true[164] 0.0029434441 2.700648e-04 0.0025215720 0.003363746 5842.937
M_true[165] 0.0033507184 3.271177e-04 0.0028392695 0.003860236 4454.267
M_true[166] 0.0063850683 5.773970e-04 0.0054606781 0.007314350 4608.478
M_true[167] 0.0064372238 5.972085e-04 0.0054688097 0.007427873 5125.998
M_true[168] 0.1110327053 1.096326e-02 0.0936650190 0.128780385 5002.494
M_true[169] 0.0878963590 8.510442e-03 0.0745290905 0.101560850 5301.929
M_true[170] 0.0009667255 9.704217e-05 0.0008089893 0.001120210 4625.569
M_true[171] 0.0008680478 8.382024e-05 0.0007257090 0.001002774 3178.937
M_true[172] 0.0009747715 9.358022e-05 0.0008231044 0.001120076 4790.759
M_true[173] 0.1173411332 1.179082e-02 0.0990117515 0.136820390 4913.305
M_true[174] 0.0481464885 5.088016e-03 0.0399478245 0.056162464 4459.045
M_true[175] 0.0764081500 7.633404e-03 0.0643051320 0.088423083 5720.210
M_true[176] 0.0866757924 8.740503e-03 0.0729003225 0.100575090 5683.205
M_true[177] 0.0530617235 5.169180e-03 0.0448247865 0.061317172 3356.499
M_true[178] 0.0569924843 5.634844e-03 0.0477599575 0.065910313 4385.082
M_true[179] 0.0898655397 8.890999e-03 0.0755746165 0.104114510 5428.557
M_true[180] 0.0474043347 4.704109e-03 0.0398315630 0.054710805 4639.041
M_true[181] 0.0258553452 2.651107e-03 0.0215936340 0.030070616 4128.754
M_true[182] 0.0265882314 2.706218e-03 0.0222069710 0.030818392 5536.600
a           0.4191816405 5.932353e-02 0.3262413600 0.513002280 2767.671
b           0.7838324045 1.417906e-02 0.7607243800 0.806090665 2770.061
sigma       0.2627569540 1.680301e-02 0.2372583950 0.290605660 2309.916
            Rhat4    
B_true[1]   0.9985442
B_true[2]   0.9997402
B_true[3]   0.9983380
B_true[4]   0.9990252
B_true[5]   0.9981828
B_true[6]   0.9987771
B_true[7]   0.9984331
B_true[8]   0.9983743
B_true[9]   0.9992680
B_true[10]  0.9983893
B_true[11]  0.9986528
B_true[12]  1.0000221
B_true[13]  0.9983626
B_true[14]  0.9982209
B_true[15]  0.9983817
B_true[16]  0.9990960
B_true[17]  0.9988796
B_true[18]  0.9992431
B_true[19]  0.9995469
B_true[20]  0.9984028
B_true[21]  0.9990997
B_true[22]  0.9989766
B_true[23]  0.9993852
B_true[24]  0.9986114
B_true[25]  0.9985109
B_true[26]  0.9993952
B_true[27]  0.9986008
B_true[28]  0.9988681
B_true[29]  0.9993563
B_true[30]  0.9993881
⋮           ⋮        
M_true[156] 0.9992022
M_true[157] 0.9989419
M_true[158] 0.9991944
M_true[159] 0.9983307
M_true[160] 0.9986911
M_true[161] 0.9986974
M_true[162] 0.9994798
M_true[163] 0.9998614
M_true[164] 0.9984022
M_true[165] 0.9991358
M_true[166] 0.9987511
M_true[167] 0.9982739
M_true[168] 0.9986564
M_true[169] 0.9987730
M_true[170] 0.9986466
M_true[171] 1.0001176
M_true[172] 0.9992279
M_true[173] 0.9987013
M_true[174] 0.9988869
M_true[175] 0.9985858
M_true[176] 0.9982319
M_true[177] 0.9990161
M_true[178] 0.9993060
M_true[179] 0.9996466
M_true[180] 0.9993738
M_true[181] 0.9995555
M_true[182] 0.9989652
a           0.9989791
b           0.9991454
sigma       1.0007652

Little has changed with the update. It looks like sigma is slightly smaller, perhaps because more uncertainty in the data has naturally led to more confidence in the model.

15H5. Now consider missing values - this data set is lousy with them. You can ignore measurement error in this problem. Let’s get a quick idea of the missing values by counting them in each variable: R code

R code 15.37

library(rethinking)
data(Primates301)
d <- Primates301
colSums( is.na(d) )

We’ll continue to focus on just brain and body, to stave off insanity. Consider only those species with measured body masses:

R code 15.38

cc <- complete.cases( d$body )
M <- d$body[cc]
M <- M / max(M)
B <- d$brain[cc]
B <- B / max( B , na.rm=TRUE )

You should end up with 238 species and 56 missing brain values among them.

First, consider whether there is a pattern to the missing values. Does it look like missing values are associated with particular values of body mass? Draw a DAG that represents how missingness works in this case. Which type (MCAR, MAR, MNAR) is this?

Second, impute missing values for brain size. It might help to initialize the 56 imputed variables to a valid value:

R code 15.39

start=list( B_impute=rep(0.5,56) )

This just helps the chain get started.

Compare the inferences to an analysis that drops all the missing values. Has anything changed? Why or why not? Hint: Consider the density of data in the ranges where there are missing values. You might want to plot the imputed brain sizes together with the observed values.

ERROR:

by counting them in each variable: R code

load.d15h5 <- function() {
  data(Primates301)
  d <- Primates301
  display_markdown("**Answer.** The `colSums` output:")
  display(colSums(is.na(d)), mimetypes="text/plain")

  cc <- complete.cases( d$body )
  M <- d$body[cc]
  M <- M / max(M)
  B <- d$brain[cc]
  B <- B / max( B , na.rm=TRUE )
  d15h5 <- data.frame(B = B, M = M)
  display_markdown("<br/> To confirm there are 238 species, preview the data.frame:")
  display(d15h5)
  display_markdown("<br/> To confirm there are 56 missing brain values, check the `summary`:")
  display(summary(d15h5))
  return(d15h5)
}

e15h5a <- r"(
The negative intercept `a` merely indicates that most values are 0 rather than 1 in D (not missing).
More interesting is `bM`, which is reliably negative. A negative `bM` indicates that larger values
of `M` are associated with being missing less. That is, brain sizes are missing more often for
species with small masses.

Clearly this is not MCAR. See section **15.2.2** for potential DAGs. If you believe that brain
volumes are only missing because researchers prefer to study species with a larger mass, then this
may be a salvageable MAR situation. If you believe brain volumes are missing because researchers
prefer to study species with larger brains, this is MNAR.
)"
q15h5a <- function(d) {
  dat <- list(
    D = as.integer(is.na(d$B)),
    M = as.numeric(standardize(d$M))
  )
  display_markdown("<br/> Let's check if we can predict `D` (deleted values) from `M`:")
  m15h5a <- ulam(
    alist(
      D ~ dbinom(1, p),
      logit(p) <- a + bM*M,
      c(a, bM) ~ dnorm(0, 1)
    ), cmdstan = TRUE, data=dat, cores=4, chains=4
  )
  display_precis(m15h5a, "m15h5a", ar=4.0)
  display_markdown(e15h5a)
}

q15h5b <- function(d) {
  dat <- list(
    B = d$B,
    M = d$M
  )
  # browser()
  m15h5b <- ulam(
    alist(
      B ~ dlnorm(mu, sigma),
      mu <- a + b * log(M),
      a ~ normal(0, 1),
      b ~ normal(0, 1),
      sigma ~ exponential(1)
    ), cmdstan=TRUE, data=dat, cores=4, chains=4, start=list( B_impute=rep(0.5,56) )
  )
  display_precis(m15h5b, "m15h5b", ar=0.7)
  return(m15h5b)
}

q15h5cc <- function(d) {
  d_cc <- d[complete.cases(d$B),]
  dat <- list(B = d_cc$B, M = d_cc$M)
  m15H4cc <- ulam(
    alist(
      B ~ dlnorm(mu, sigma),
      mu <- a + b * log(M),
      a ~ normal(0, 1),
      b ~ normal(0, 1),
      sigma ~ exponential(1)
    ), cmdstan=TRUE, data=dat, chains=4, cores=4
  )
  display_precis(m15H4cc, "m15H4cc", 4.0)
  return(m15H4cc)
}

e15h5 <- r"(
There doesn't seem to be any change in the inferences for `a`, `b`, and `sigma`. When you impute
output values, the model isn't really doing anything but posterior predictive inference. We haven't
supplied any new assumptions about what the distribution of `B` should look like, other than what we
know from `M` (which we only learned from complete cases).
)"
q15h5 <- function() {
  d15h5 <- load.d15h5()
  q15h5a(d15h5)
  m15h5b <- q15h5b(d15h5)
  m15H4cc <- q15h5cc(d15h5)

  iplot(function() {
    plot(coeftab(m15H4cc, m15h5b), pars=c("a", "b", "sigma"), main="coeftab")
  }, ar=2)
  display_markdown(e15h5)
}
q15h5()

Answer. The colSums output:

               name               genus             species          subspecies 
                  0                   0                   0                 267 
             spp_id            genus_id     social_learning     research_effort 
                  0                   0                  98                 115 
              brain                body          group_size           gestation 
                117                  63                 114                 161 
            weaning           longevity        sex_maturity maternal_investment 
                185                 181                 194                 197 

To confirm there are 238 species, preview the data.frame:

A data.frame: 238 × 2

B	M
<dbl>	<dbl>
0.11810206	0.0358076923
NA	0.0006006923
0.10755796	0.0491923077
0.10713050	0.0414076923
0.10523745	0.0398076923
0.10153276	0.0480769231
0.10407719	0.0685769231
0.12025973	0.0508541538
0.11240255	0.0457692308
0.04207462	0.0092692308
0.03317931	0.0056461538
NA	0.0060848462
NA	0.0073692308
0.03429886	0.0076076923
NA	0.0054076923
0.01196898	0.0016153846
0.01408594	0.0023769231
0.23819895	0.0628230769
0.23254015	0.0694230769
0.21391496	0.0579615385
0.21139088	0.0636923077
0.02007043	0.0092846154
0.01618255	0.0061615385
NA	0.0810562308
0.22529363	0.0517538462
0.15470108	0.0243461538
0.13998412	0.0225769231
NA	0.0069051538
NA	0.0082123846
NA	0.0073702308
⋮	⋮
0.016162192	0.0030846154
0.020640381	0.0039769231
NA	0.0031377692
0.019744743	0.0040384615
0.019907586	0.0043307692
0.022574145	0.0044923077
0.019296924	0.0028846154
0.019866876	0.0033153846
NA	0.0029619231
NA	0.0061496154
0.051031001	0.0060692308
0.049137949	0.0061461538
NA	0.0068190000
0.225802512	0.1134000000
0.251389256	0.0868846154
0.006432308	0.0009692308
0.006106622	0.0008692308
0.006839416	0.0009692308
0.271398620	0.1180769231
NA	0.0747661538
0.117776376	0.0491846154
NA	0.0626148462
0.165489446	0.0780769231
0.172206729	0.0892307692
0.126447778	0.0542769231
0.148268773	0.0575000000
0.210963421	0.0907230769
0.124758280	0.0479769231
0.063264600	0.0266923077
0.065381562	0.0275000000

To confirm there are 56 missing brain values, check the summary:

       B                 M            
 Min.   :0.00332   Min.   :0.0002402  
 1st Qu.:0.02440   1st Qu.:0.0056880  
 Median :0.12075   Median :0.0273346  
 Mean   :0.14075   Mean   :0.0522706  
 3rd Qu.:0.17789   3rd Qu.:0.0574231  
 Max.   :1.00000   Max.   :1.0000000  
 NA's   :56                           

Let's check if we can predict D (deleted values) from M:

Warning in '/tmp/RtmpUakS9C/model-132bf8a4d7.stan', line 2, column 4: Declaration

    of arrays by placing brackets after a variable name is deprecated and

    will be removed in Stan 2.32.0. Instead use the array keyword before the

    type. This can be changed automatically using the auto-format flag to

    stanc

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.2 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 0.2 seconds.
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 0.2 seconds.
Chain 4 finished in 0.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.2 seconds.
Total execution time: 0.4 seconds.

Raw data (following plot):

   mean       sd        5.5%      94.5%       n_eff    Rhat4    
bM -0.5072257 0.3028944 -1.048434 -0.06574995 1014.347 0.9996125
a  -1.2086073 0.1610078 -1.474740 -0.95873821 1127.009 0.9996780

The negative intercept a merely indicates that most values are 0 rather than 1 in D (not missing). More interesting is bM, which is reliably negative. A negative bM indicates that larger values of M are associated with being missing less. That is, brain sizes are missing more often for species with small masses.

Clearly this is not MCAR. See section 15.2.2 for potential DAGs. If you believe that brain volumes are only missing because researchers prefer to study species with a larger mass, then this may be a salvageable MAR situation. If you believe brain volumes are missing because researchers prefer to study species with larger brains, this is MNAR.

Found 56 NA values in B and attempting imputation.

Raw data (following plot):

             mean        sd           5.5%        94.5%       n_eff    
a            0.429209747 0.0570991572 0.339948276 0.520651112  889.8638
b            0.783647438 0.0138930038 0.762394288 0.805935970  948.3123
sigma        0.293430363 0.0156388062 0.269724497 0.318633929 1072.5169
B_impute[1]  0.004786200 0.0014914046 0.002857822 0.007525111 1137.7001
B_impute[2]  0.029086884 0.0083617833 0.018039909 0.043735841 1575.0531
B_impute[3]  0.034563101 0.0103652766 0.021077614 0.053294099 1017.0016
B_impute[4]  0.026491233 0.0082289481 0.015560821 0.040971634 1489.8989
B_impute[5]  0.225915045 0.0726967881 0.130827723 0.349635062  836.4177
B_impute[6]  0.032250991 0.0094624808 0.019980671 0.048312020 1516.9083
B_impute[7]  0.037880390 0.0118272347 0.022293129 0.058300099 1398.9808
B_impute[8]  0.034554255 0.0108554389 0.020622710 0.053117521 1520.0070
B_impute[9]  0.045598981 0.0131587578 0.027158449 0.068392413 1863.5973
B_impute[10] 0.041774107 0.0128845153 0.024896337 0.064855272 1125.7264
B_impute[11] 0.018121363 0.0054244604 0.010799199 0.027650251 2035.6780
B_impute[12] 0.014231789 0.0044881272 0.008344933 0.022448822 1586.1794
B_impute[13] 0.015329548 0.0046909879 0.009262177 0.023614649 1571.3961
B_impute[14] 0.016122983 0.0047035829 0.009709436 0.024484244 2297.1859
B_impute[15] 0.016484752 0.0049648826 0.009806028 0.025120559 1808.9191
B_impute[16] 0.018505499 0.0053459654 0.011731652 0.027859439 1797.5330
B_impute[17] 0.093398955 0.0280276883 0.055897370 0.143725021 1457.8011
B_impute[18] 0.130130424 0.0408690371 0.076146237 0.203106584  999.3712
B_impute[19] 0.129359063 0.0383762879 0.078740431 0.199986927 1201.2827
B_impute[20] 0.070147251 0.0211930144 0.042256754 0.106495169  866.6901
B_impute[21] 0.012061878 0.0037404611 0.007220431 0.018747894 1621.0253
B_impute[22] 0.050166179 0.0152995489 0.030274479 0.077493567 1254.9425
B_impute[23] 0.860432847 0.2620480707 0.513400500 1.325823939 1116.8501
B_impute[24] 0.142664219 0.0460889820 0.082237482 0.215528521  812.1062
B_impute[25] 0.025412435 0.0075556927 0.015203344 0.039167057 1689.2951
B_impute[26] 0.028134694 0.0081404623 0.017316218 0.042867656 1299.7887
B_impute[27] 0.169314478 0.0546392671 0.099653470 0.260259193  819.2582
B_impute[28] 0.093156644 0.0287961199 0.056129526 0.143993327  982.3735
B_impute[29] 0.116278165 0.0360459097 0.068648519 0.179953866 1137.0522
B_impute[30] 0.221431434 0.0632861315 0.138755290 0.328337691 1196.3263
B_impute[31] 0.216952919 0.0652025838 0.127738585 0.333300396 1212.2576
B_impute[32] 0.002483832 0.0007899518 0.001431282 0.003865538 1293.6024
B_impute[33] 0.004399372 0.0013400503 0.002625177 0.006772746 1195.9556
B_impute[34] 0.002301103 0.0006781731 0.001391275 0.003527719 1508.4254
B_impute[35] 0.003809495 0.0011719956 0.002227190 0.006039362 1356.6416
B_impute[36] 0.003345134 0.0010734592 0.001956075 0.005216269 1472.9661
B_impute[37] 0.004320865 0.0012847241 0.002571459 0.006529500 1509.1870
B_impute[38] 0.041955488 0.0121875359 0.025458276 0.063306769 1582.9235
B_impute[39] 0.152073431 0.0461289557 0.090522707 0.234132139 1215.9228
B_impute[40] 0.167228264 0.0503866096 0.100458543 0.253082815 1251.5582
B_impute[41] 0.017628984 0.0054501432 0.010653640 0.026924360 1646.1347
B_impute[42] 0.229739458 0.0675851325 0.139328718 0.349393552 1091.0886
B_impute[43] 0.193920914 0.0597911106 0.116504261 0.302735805 1090.9389
B_impute[44] 0.180186612 0.0512559379 0.109414988 0.272320041 1329.3292
B_impute[45] 0.068089177 0.0210701081 0.040722215 0.105434886  853.7422
B_impute[46] 0.094615481 0.0285466364 0.056778052 0.145485725 1245.0110
B_impute[47] 0.205792596 0.0675916585 0.119281179 0.324160656  803.9435
B_impute[48] 0.233534784 0.0740091353 0.136150060 0.369226377 1175.6142
B_impute[49] 0.252462732 0.0768465961 0.147657038 0.381118539 1226.3441
B_impute[50] 0.019421691 0.0058181517 0.011540242 0.029780839 2163.5233
B_impute[51] 0.017648020 0.0053207744 0.010650871 0.027172169 1294.9259
B_impute[52] 0.016768843 0.0050340439 0.009744324 0.025868531 2547.8019
B_impute[53] 0.029890798 0.0091762387 0.017554664 0.046480812 1613.0035
B_impute[54] 0.032199437 0.0098086945 0.018890907 0.049411995 1640.9744
B_impute[55] 0.209300633 0.0602691992 0.129023668 0.313635790 1205.6605
B_impute[56] 0.183497458 0.0561860323 0.109550257 0.283267329  976.8538
             Rhat4    
a            1.0012201
b            1.0007573
sigma        1.0025655
B_impute[1]  1.0003645
B_impute[2]  1.0007252
B_impute[3]  1.0054650
B_impute[4]  1.0000994
B_impute[5]  0.9991215
B_impute[6]  1.0013527
B_impute[7]  1.0005014
B_impute[8]  1.0022631
B_impute[9]  0.9998690
B_impute[10] 1.0006931
B_impute[11] 1.0001202
B_impute[12] 1.0011720
B_impute[13] 0.9994300
B_impute[14] 0.9999795
B_impute[15] 0.9991311
B_impute[16] 0.9996744
B_impute[17] 1.0008510
B_impute[18] 1.0034034
B_impute[19] 1.0015601
B_impute[20] 1.0010288
B_impute[21] 1.0001575
B_impute[22] 0.9998241
B_impute[23] 1.0017472
B_impute[24] 1.0001907
B_impute[25] 1.0026519
B_impute[26] 1.0000780
B_impute[27] 1.0001801
B_impute[28] 1.0011353
B_impute[29] 1.0003357
B_impute[30] 1.0053738
B_impute[31] 1.0007850
B_impute[32] 1.0033039
B_impute[33] 1.0031479
B_impute[34] 0.9997142
B_impute[35] 1.0052562
B_impute[36] 1.0018235
B_impute[37] 1.0016042
B_impute[38] 0.9997651
B_impute[39] 0.9999042
B_impute[40] 0.9998815
B_impute[41] 0.9989642
B_impute[42] 1.0030449
B_impute[43] 1.0011753
B_impute[44] 1.0005502
B_impute[45] 1.0019154
B_impute[46] 0.9996109
B_impute[47] 1.0002261
B_impute[48] 1.0019765
B_impute[49] 0.9993826
B_impute[50] 0.9998903
B_impute[51] 1.0008995
B_impute[52] 0.9989648
B_impute[53] 0.9990433
B_impute[54] 0.9997337
B_impute[55] 0.9995189
B_impute[56] 1.0000349

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Scale parameter is inf, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Scale parameter is inf, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Scale parameter is inf, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: lognormal_lpdf: Scale parameter is 0, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: lognormal_lpdf: Scale parameter is 0, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 

Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 3 Exception: lognormal_lpdf: Scale parameter is inf, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 3 

Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: lognormal_lpdf: Scale parameter is 0, but must be positive finite! (in '/tmp/RtmpUakS9C/model-135ebb0538.stan', line 18, column 4 to column 32)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.3 seconds.
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 0.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.5 seconds.

Raw data (following plot):

      mean      sd         5.5%      94.5%     n_eff    Rhat4   
a     0.4267492 0.05525369 0.3397942 0.5146442 923.0049 1.003886
b     0.7830870 0.01328885 0.7622681 0.8047012 899.6270 1.003091
sigma 0.2929341 0.01539626 0.2687894 0.3185917 705.9164 1.001046

There doesn't seem to be any change in the inferences for a, b, and sigma. When you impute output values, the model isn't really doing anything but posterior predictive inference. We haven't supplied any new assumptions about what the distribution of B should look like, other than what we know from M (which we only learned from complete cases).

15H6. Return to the divorce rate measurement error model. This time try to incorporate the full generative system: \(A \rightarrow M \rightarrow D\), \(A \rightarrow D\). What this means is that the prior for \(M\) should include \(A\) somehow, because it is influenced by \(A\).

load.d15h6 <- function() {
  data(WaffleDivorce)
  d <- WaffleDivorce

  dat <- list(
    D_obs = standardize(d$Divorce),
    D_sd = d$Divorce.SE / sd(d$Divorce),
    M_obs = standardize(d$Marriage),
    M_sd = d$Marriage.SE / sd(d$Marriage),
    A = standardize(d$MedianAgeMarriage),
    N = nrow(d)
  )
  return(dat)
}

a15h6a <- function(dat) {
  m15.2 <- ulam(
    alist(
      D_obs ~ dnorm(D_true, D_sd),
      vector[N]:D_true ~ dnorm(mu, sigma),
      mu <- a + bA * A + bM * M_true[i],
      M_obs ~ dnorm(M_true, M_sd),
      vector[N]:M_true ~ dnorm(0, 1),
      a ~ dnorm(0, 0.2),
      bA ~ dnorm(0, 0.5),
      bM ~ dnorm(0, 0.5),
      sigma ~ dexp(1)
    ),
    cmdstan = TRUE, data = dat, chains = 4, cores = 4
  )
  display_precis(m15.2, "m15.2", ar=0.4)
  return(m15.2)
}

q15h6b <- function(dat) {
  m15h6 <- ulam(
    alist(
      D_obs ~ dnorm(D_true, D_sd),
      vector[N]:D_true ~ dnorm(mu, sigma),
      mu <- a + bA * A + bM * M_true[i],
      M_obs ~ dnorm(M_true, M_sd),
      vector[N]:M_true ~ dnorm(mu_M, 1),
      mu_M <- aM + bAM * A,
      aM ~ dnorm(0, 0.2),
      a ~ dnorm(0, 0.2),
      bA ~ dnorm(0, 0.5),
      bM ~ dnorm(0, 0.5),
      bAM ~ dnorm(0, 0.5),
      sigma ~ dexp(1)
    ),
    cmdstan = TRUE, data = dat, chains = 4, cores = 4
  )
  display_precis(m15h6, "m15h6", ar=0.4)
  return(m15h6)
}

e15h6b <- r"(
Notice `bAM` is nonzero. However, most of the original parameters are unchanged, as expected based
on the discussion in section **5.1.4** (`M` adds little predictive power once you know `A`).
)"
q15h6 <- function() {
  dat <- load.d15h6()
  display_markdown("**Answer.** Let's start by reproducing results from the chapter:")
  m15.2 <- a15h6a(dat)
  display_markdown("Fitting the new model:")
  m15h6 <- q15h6b(dat)

  iplot(function() {
    plot(coeftab(m15.2, m15h6), pars=c("a", "bA", "bM", "sigma"), main="coeftab")
  }, ar=2)
  display_markdown(e15h6b)
}
q15h6()

Answer. Let's start by reproducing results from the chapter:

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-137bbf99b3.stan', line 28, column 4 to column 34)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 4 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-137bbf99b3.stan', line 28, column 4 to column 34)

Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 4 

Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 finished in 0.3 seconds.
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 finished in 0.4 seconds.
Chain 3 finished in 0.4 seconds.
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 0.4 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.4 seconds.
Total execution time: 0.6 seconds.

Raw data (following plot):

           mean          sd         5.5%        94.5%       n_eff     Rhat4    
D_true[1]   1.12385798   0.3554066   0.54481399  1.68852375 1914.401  0.9998953
D_true[2]   0.74050110   0.5266029  -0.08015318  1.59452610 2444.684  0.9998991
D_true[3]   0.42053533   0.3161150  -0.07594199  0.93404446 2597.937  0.9991635
D_true[4]   1.46071123   0.4480914   0.74302428  2.18795620 2510.655  0.9986221
D_true[5]  -0.89449025   0.1259901  -1.09231825 -0.68827442 2127.438  0.9997797
D_true[6]   0.68512870   0.3867411   0.07769300  1.30872455 2292.443  0.9995023
D_true[7]  -1.35476021   0.3430152  -1.90592495 -0.80999384 2481.670  0.9996085
D_true[8]  -0.28698893   0.4548214  -1.01598670  0.41175814 2475.836  0.9990076
D_true[9]  -1.78497443   0.6018430  -2.71190090 -0.82044548 1997.348  0.9997020
D_true[10] -0.62441650   0.1650684  -0.89387593 -0.36438402 2698.861  0.9987466
D_true[11]  0.77751538   0.2935488   0.31978693  1.25578840 2016.897  0.9987995
D_true[12] -0.45209473   0.4729781  -1.20371080  0.29501435 1803.887  1.0021392
D_true[13]  0.24294624   0.4715200  -0.52307459  0.96934553 1402.331  1.0011768
D_true[14] -0.86766586   0.2254476  -1.22821000 -0.49401419 2411.601  0.9997595
D_true[15]  0.54340814   0.3058431   0.05794533  1.03787930 2339.175  0.9983663
D_true[16]  0.28457517   0.3785022  -0.32169396  0.88385075 2560.544  1.0009084
D_true[17]  0.51067483   0.4163774  -0.15068532  1.16342025 3012.344  0.9995116
D_true[18]  1.23537672   0.3389570   0.70068147  1.78937960 2185.490  1.0013140
D_true[19]  0.43054590   0.3816789  -0.16673172  1.05539540 2935.632  0.9996269
D_true[20]  0.24361644   0.5390114  -0.60930384  1.11169005 1259.548  1.0034180
D_true[21] -0.55886426   0.3156862  -1.06725495 -0.04662007 2947.657  0.9994248
D_true[22] -1.10886942   0.2473762  -1.50504785 -0.71464697 2560.147  1.0022444
D_true[23] -0.30507948   0.2487208  -0.71469129  0.08518569 2185.574  1.0019136
D_true[24] -1.02378336   0.2922501  -1.48872680 -0.57042309 2238.040  0.9990170
D_true[25]  0.40349368   0.3942681  -0.20790610  1.02558770 1876.685  0.9994361
D_true[26] -0.04755081   0.3053162  -0.54331310  0.43034702 2843.097  0.9993076
D_true[27] -0.04985283   0.5001203  -0.85829017  0.74741533 2852.683  1.0006701
D_true[28] -0.14513949   0.4056516  -0.81104229  0.48514383 2881.513  0.9993679
D_true[29] -0.30315273   0.4723143  -1.03334420  0.47123007 1925.885  0.9998698
D_true[30] -1.80815678   0.2307198  -2.17006860 -1.43145570 2293.208  0.9984504
⋮          ⋮             ⋮          ⋮           ⋮           ⋮         ⋮        
M_true[25] -0.1448684243 0.37293504 -0.74455669  0.43709568 1892.4066 0.9989351
M_true[26] -0.3811124354 0.20494118 -0.71277557 -0.04695257 2586.3208 0.9997921
M_true[27] -0.3235572561 0.52736772 -1.15066815  0.53674010 2911.4103 0.9982081
M_true[28] -0.1433853533 0.35678454 -0.71350194  0.44161150 2615.8537 0.9991915
M_true[29] -0.7045548847 0.40187257 -1.33865510 -0.03454238 2315.0763 0.9989793
M_true[30] -1.3776183825 0.15088344 -1.62501045 -1.13820930 2282.7380 0.9999063
M_true[31]  0.0660791122 0.43749275 -0.61072549  0.75275760 2219.3813 1.0007310
M_true[32] -0.8663235600 0.12741769 -1.06660815 -0.66400894 2411.1577 0.9998176
M_true[33]  0.0717017755 0.24752616 -0.32552679  0.45669929 3181.0730 0.9986132
M_true[34]  0.9689018099 0.62724536 -0.03917871  1.99773575 2405.8041 0.9990945
M_true[35] -0.8213976310 0.16609261 -1.09209135 -0.55897318 2767.5670 0.9994034
M_true[36]  0.9071350519 0.33482186  0.35913469  1.45461885 2784.3482 0.9983553
M_true[37] -0.2667911624 0.27913679 -0.71303866  0.16792863 2635.3455 1.0004072
M_true[38] -1.1997689575 0.12348601 -1.39374705 -1.00127745 3088.8690 0.9997407
M_true[39] -1.0074613184 0.46519207 -1.76555605 -0.27849233 2504.4498 1.0007309
M_true[40] -0.5011537417 0.29221312 -0.95379836 -0.02462128 2572.5901 1.0001853
M_true[41] -0.0007725454 0.57542887 -0.92589931  0.88127204 2416.7046 0.9987366
M_true[42] -0.1637778936 0.22125186 -0.51520496  0.18682741 2693.1536 0.9985345
M_true[43]  0.3513243068 0.15965578  0.10107341  0.60185929 2658.6843 0.9982000
M_true[44]  1.9525186970 0.43109855  1.29961745  2.64893650 2469.4158 0.9997747
M_true[45] -0.6865810001 0.52788819 -1.50918005  0.17976398 2595.0752 0.9986795
M_true[46]  0.0930301573 0.20905387 -0.24043334  0.43013554 2097.9095 0.9993468
M_true[47]  0.3251173919 0.24122606 -0.05922683  0.71370753 2507.4179 0.9990207
M_true[48]  0.4699208470 0.39142002 -0.15237856  1.09293820 2236.3509 1.0014293
M_true[49] -0.7466404448 0.19900366 -1.06804200 -0.43654851 2829.3651 0.9995158
M_true[50]  1.2858845407 0.70408722  0.17539639  2.45705750 2117.5642 1.0017119
a          -0.0425067809 0.09662695 -0.19308193  0.11207957 1432.8261 0.9984967
bA         -0.5483560592 0.16292837 -0.81169573 -0.28882952 1081.4427 1.0017287
bM          0.2020118926 0.21267107 -0.12967567  0.53664554  758.0593 1.0051126
sigma       0.5510849735 0.10375394  0.38863759  0.72575790  619.3232 1.0077815

Fitting the new model:

Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpUakS9C/model-13186b2286.stan', line 36, column 4 to column 34)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 finished in 0.5 seconds.
Chain 2 finished in 0.4 seconds.
Chain 3 finished in 0.4 seconds.
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 0.5 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.5 seconds.
Total execution time: 0.7 seconds.

Raw data (following plot):

           mean        sd         5.5%        94.5%       n_eff     Rhat4    
D_true[1]   1.15741184 0.3595291   0.59412252  1.74972665 2074.339  1.0006261
D_true[2]   0.73003959 0.5469587  -0.12951053  1.60639965 2306.957  0.9987878
D_true[3]   0.43048803 0.3482953  -0.11874783  0.99947320 2745.864  0.9996444
D_true[4]   1.45810499 0.4762043   0.70619835  2.21900915 2612.755  0.9990083
D_true[5]  -0.89844902 0.1289282  -1.10140970 -0.69021863 2891.040  0.9991381
D_true[6]   0.68066719 0.3973069   0.06299353  1.30106485 3260.033  0.9987280
D_true[7]  -1.36457105 0.3540565  -1.94241340 -0.80054167 2338.352  0.9995367
D_true[8]  -0.32822943 0.4806232  -1.12201045  0.41759172 1921.289  1.0005649
D_true[9]  -1.86883291 0.6021902  -2.81836470 -0.88254134 1784.200  1.0021782
D_true[10] -0.62529896 0.1579800  -0.87109117 -0.37367789 3233.745  0.9991621
D_true[11]  0.76868559 0.2796956   0.33666835  1.22805045 1956.473  0.9994012
D_true[12] -0.52527207 0.4899950  -1.29960710  0.24368810 2053.003  1.0006026
D_true[13]  0.16622776 0.4815465  -0.61262437  0.90119968 1303.675  0.9997353
D_true[14] -0.86828884 0.2237228  -1.23626130 -0.50700680 2547.783  0.9999712
D_true[15]  0.53343476 0.2985346   0.06041430  1.02060690 2433.651  0.9997813
D_true[16]  0.28778024 0.3777882  -0.31387357  0.88401106 2227.041  0.9993436
D_true[17]  0.51062931 0.4274562  -0.13550758  1.19631340 2896.135  0.9992179
D_true[18]  1.25345616 0.3455893   0.71119767  1.81143085 2220.506  0.9995943
D_true[19]  0.43898114 0.3791973  -0.15190421  1.04743770 2311.989  0.9985908
D_true[20]  0.35956008 0.5607253  -0.47115521  1.29438130 1256.289  1.0016227
D_true[21] -0.55297449 0.3135544  -1.06210115 -0.04361986 2774.016  1.0009602
D_true[22] -1.09450045 0.2523780  -1.48913400 -0.69456875 2808.485  0.9988513
D_true[23] -0.28616501 0.2611433  -0.71201247  0.13880295 2067.197  0.9988679
D_true[24] -1.00610083 0.2906452  -1.47546390 -0.54541531 1919.609  0.9992381
D_true[25]  0.42796329 0.4121669  -0.23009722  1.09007990 2678.018  1.0002640
D_true[26] -0.04345925 0.3177689  -0.54807182  0.43438509 2494.603  1.0001580
D_true[27] -0.03687699 0.4940858  -0.84317314  0.75180598 2415.242  0.9996449
D_true[28] -0.16204792 0.3966675  -0.79008256  0.48044281 2991.713  0.9985323
D_true[29] -0.27487396 0.5107981  -1.05462010  0.56081215 2084.902  1.0010248
D_true[30] -1.81449135 0.2447264  -2.20513550 -1.42679180 2331.857  0.9990249
⋮          ⋮           ⋮          ⋮           ⋮           ⋮         ⋮        
M_true[27] -0.28923060 0.54690892 -1.18023410  0.58945606 2698.8430 0.9987006
M_true[28] -0.09182479 0.35302763 -0.65251846  0.45166202 2785.4455 0.9987070
M_true[29] -0.80201348 0.41152929 -1.44995700 -0.14778865 2602.0373 0.9989405
M_true[30] -1.39175353 0.15530592 -1.63198830 -1.14317020 2697.6732 0.9998723
M_true[31]  0.08217810 0.42530037 -0.58524281  0.76108709 2623.9126 1.0007984
M_true[32] -0.87636522 0.12236965 -1.06934135 -0.68060106 2646.8222 0.9995971
M_true[33]  0.07919226 0.25055348 -0.33089932  0.48249644 2726.8911 0.9987757
M_true[34]  1.15444743 0.60967329  0.15492086  2.12917970 3131.7459 0.9990234
M_true[35] -0.83005297 0.15572244 -1.07748365 -0.58163307 3003.4485 0.9986447
M_true[36]  0.97891253 0.32159342  0.43837574  1.49501505 2601.8884 1.0001305
M_true[37] -0.28830185 0.27470838 -0.73073531  0.15900299 3000.0338 0.9985948
M_true[38] -1.20582236 0.12291540 -1.39910685 -1.00237735 3103.4038 0.9984800
M_true[39] -1.27723121 0.47164140 -2.01509910 -0.52995943 2484.1396 0.9983936
M_true[40] -0.51197200 0.29538098 -0.99158171 -0.04300853 3171.2424 1.0000859
M_true[41]  0.06659838 0.58581477 -0.86348428  1.00052695 2868.5150 0.9995583
M_true[42] -0.14762759 0.21755571 -0.49363252  0.18497697 2665.3211 1.0004558
M_true[43]  0.36461604 0.15730744  0.11083406  0.61445859 2989.6685 0.9998724
M_true[44]  2.23985825 0.43964666  1.54402570  2.93478195 2075.8127 0.9995402
M_true[45] -0.82197830 0.56628402 -1.75761095  0.08632466 2600.6192 0.9986284
M_true[46]  0.07905651 0.21048304 -0.25326136  0.40616174 2929.1503 0.9989778
M_true[47]  0.32459071 0.24795602 -0.07911612  0.72100663 3165.4918 0.9992154
M_true[48]  0.55095543 0.38791439 -0.08567644  1.15856975 2700.9058 1.0004629
M_true[49] -0.75362239 0.20906057 -1.07695540 -0.41069346 2919.1697 0.9992028
M_true[50]  1.79173403 0.71098989  0.65998150  2.93193095 2709.0857 0.9992667
aM         -0.03838471 0.12467626 -0.23042831  0.16585141 2454.9284 0.9996420
a          -0.04708609 0.09969729 -0.20777371  0.10764634 1872.9229 0.9986087
bA         -0.56984227 0.17168467 -0.83871434 -0.29010205  868.6785 1.0064443
bM          0.12587028 0.20792856 -0.19685901  0.46149217  696.3589 1.0062410
bAM        -0.63835680 0.14670222 -0.87692010 -0.40137013 1877.6531 0.9994998
sigma       0.58195392 0.11223721  0.41170270  0.76920186  613.5693 1.0022078

Notice bAM is nonzero. However, most of the original parameters are unchanged, as expected based on the discussion in section 5.1.4 (M adds little predictive power once you know A).

15H7. Some lad named Andrew made an eight-sided spinner. He wanted to know if it is fair. So he spun it a bunch of times, recording the counts of each value. Then he accidentally spilled coffee over the 4s and 5s. The surviving data are summarized below.

library(data.table)

load.d15h7 <- function() {
  V = seq(8)
  Freq = c(18, 19, 22, NA, NA, 19, 20, 22)
  d <- data.frame(Value = V, Frequency = Freq)
  t_d <- transpose(d)
  colnames(t_d) <- rownames(d)
  rownames(t_d) <- colnames(d)
  display(t_d)
  return(d)
}
e15h7b <- r"(
Your job is to impute the two missing values in the table above. Andrew doesn’t remember how many
times he spun the spinner. So you will have to assign a prior distribution for the total number of
spins and then marginalize over the unknown total. Andrew is not sure the spinner is fair (every
value is equally likely), but he’s confident that none of the values is twice as likely as any
other. Use a Dirichlet distribution to capture this prior belief. Plot the joint posterior
distribution of 4s and 5s.
)"

e15h7c <- r"(
**Answer.** We could model this as a categorical distribution where we predict the category i.e. the
spinner value from no data; we don't need a predictor to spin. To do this, we'd have to expand the
given table into many trials. How would we model a variable number of observations, though?

To model a variable number of observations, it seems like we'd need to use the multinomial
distribution. It doesn't look like this distribution is available through `ulam`. It's also not
clear how we'd deal with what is essentially a single observation in this case.

What if you took the $Value$ as a categorical predictor? The assumption that no value is twice as
likely as any other would become the variance on a prior. But the Dirichlet distribution strongly
implies either a categorical or multinomial.

To plot the joint posterior distribution, use `pairs` on the model.

Skipping this question for now.
)"
q15h7 <- function() {
  d <- load.d15h7()
  display_markdown(e15h7b)
  display_markdown(e15h7c)
}
q15h7()

A data.frame: 2 × 8

	1	2	3	4	5	6	7	8
	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
Value	1	2	3	4	5	6	7	8
Frequency	18	19	22	NA	NA	19	20	22

Your job is to impute the two missing values in the table above. Andrew doesn’t remember how many times he spun the spinner. So you will have to assign a prior distribution for the total number of spins and then marginalize over the unknown total. Andrew is not sure the spinner is fair (every value is equally likely), but he’s confident that none of the values is twice as likely as any other. Use a Dirichlet distribution to capture this prior belief. Plot the joint posterior distribution of 4s and 5s.

Answer. We could model this as a categorical distribution where we predict the category i.e. the spinner value from no data; we don't need a predictor to spin. To do this, we'd have to expand the given table into many trials. How would we model a variable number of observations, though?

To model a variable number of observations, it seems like we'd need to use the multinomial distribution. It doesn't look like this distribution is available through ulam. It's also not clear how we'd deal with what is essentially a single observation in this case.

What if you took the $Value$ as a categorical predictor? The assumption that no value is twice as likely as any other would become the variance on a prior. But the Dirichlet distribution strongly implies either a categorical or multinomial.

To plot the joint posterior distribution, use pairs on the model.

Skipping this question for now.