Practice: Chp. 15
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:
Installation instructions for cmdstan
and 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 mcstan.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.
Answer. Assuming a \(P_{SE,i}\) is available, though it is not in the data:
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:
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 nonNA 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/model1367da2ab5.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 illconditioned 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/model131e5c493f.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 illconditioned 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 noncentered 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 noncentered parameterization)
* Using informative or weakly informative prior distributions
495 of 2000 (25.0%) transitions hit the maximum treedepth limit of 10 or 2^101 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/model1353560fc6.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 illconditioned 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/model1353560fc6.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 illconditioned 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 posttreatment 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/model133612bcfe.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 autoformat flag to
stanc
Warning in '/tmp/RtmpUakS9C/model133612bcfe.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 autoformat flag to
stanc
Warning in '/tmp/RtmpUakS9C/model133612bcfe.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 autoformat 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
 4
 14
 22
 26
 28
 29
 30
 32
 41
 46
 47
 48
 53
 55
 59
 60
 65
 66
 81
 83
 84
 92
 100
The mean values of PrC1
at those indexes:
mean_prc1 < apply(post$PrC1, 2, mean)
display(mean_prc1[R_C != 0])
 0.00900593709
 0.5487595275
 0.1840805422
 0.04225329709
 0.7399022305
 0.3424969625
 0.8683926035
 0.9378788725
 0.5487595275
 0.7399022305
 0.3424969625
 0.8683926035
 0.3424969625
 0.3424969625
 0.04225329709
 0.7399022305
 0.5487595275
 0.7399022305
 0.01950867139
 0.3424969625
 0.5487595275
 0.7399022305
 0.3424969625
 0.0901580741
The true simulated values at those indexes:
display(cat[R_C != 0])
 0
 0
 0
 0
 1
 0
 1
 1
 1
 1
 0
 1
 1
 1
 0
 1
 0
 1
 0
 0
 0
 1
 0
 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)
 TRUE
 FALSE
 TRUE
 TRUE
 TRUE
 TRUE
 TRUE
 TRUE
 TRUE
 TRUE
 TRUE
 TRUE
 FALSE
 FALSE
 TRUE
 TRUE
 FALSE
 TRUE
 TRUE
 TRUE
 FALSE
 TRUE
 TRUE
 TRUE
For a total accuracy of:
display(sum(is_correct) / length(is_correct))
15M6. Return to the four dogeatshomework 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/model13625692d9.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 autoformat 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:
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/model134311df7a.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 autoformat 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 hardworking 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:
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/model13719a49e2.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 autoformat 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 nonlinear. 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/model1325283428.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 autoformat 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:
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/model1342371bb4.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 autoformat 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 nonlinear. 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/model1357ace0f7.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 autoformat 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:
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/model137399b0b3.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 autoformat 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:
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/model132a73d6c8.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 autoformat 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. Reestimate 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)
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 2025 years of age uncertainty:
q15.h1.b(elephants, error=25)
Warning in '/tmp/RtmpUakS9C/model135a0b78b7.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 autoformat 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
:
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

 0.502192350531457
 0.131531165327303
 0.07891708981887
 0.886784809417845
 0.116971270510841
 0.318630087617032
 0.58179068471591
 0.714532710891568
 0.825259425862769
 0.359862131395465
 <NA>
 $y

 0.412306206753926
 0.227805625612433
 0.280551042002224
 1.62662530929628
 0.240350771599711
 0.289313378387685
 0.97064493161941
 1.22538896826149
 1.73907361123195
 1.95043469138359
 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model134c4692d5.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 illconditioned 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/model132c5f8f42.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 autoformat flag to stanc
Warning in '/tmp/RtmpUakS9C/model132c5f8f42.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 autoformat 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:
This says that brain volume is a lognormal variable, and the mean on the log scale is given by \(\mu\). What this model implies is that the expected value of B is:
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:
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:
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/model1371183445.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 illconditioned 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/model1371183445.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 illconditioned 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/model1371183445.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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/model1311d775da.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 illconditioned 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.198526e03 0.0184078130 0.025338225 4455.942
M_true[157] 0.0740924054 7.433089e03 0.0621375165 0.086132004 4547.262
M_true[158] 0.1119229651 1.180245e02 0.0931663405 0.130695530 5311.109
M_true[159] 0.0030852139 2.975357e04 0.0026144887 0.003569972 4734.040
M_true[160] 0.0039918902 3.786749e04 0.0033638231 0.004593379 4221.715
M_true[161] 0.0040281731 4.056766e04 0.0033881246 0.004692708 5922.510
M_true[162] 0.0043007993 4.155236e04 0.0036304378 0.004958710 3708.773
M_true[163] 0.0045064638 4.245070e04 0.0038333670 0.005177733 4897.942
M_true[164] 0.0029434441 2.700648e04 0.0025215720 0.003363746 5842.937
M_true[165] 0.0033507184 3.271177e04 0.0028392695 0.003860236 4454.267
M_true[166] 0.0063850683 5.773970e04 0.0054606781 0.007314350 4608.478
M_true[167] 0.0064372238 5.972085e04 0.0054688097 0.007427873 5125.998
M_true[168] 0.1110327053 1.096326e02 0.0936650190 0.128780385 5002.494
M_true[169] 0.0878963590 8.510442e03 0.0745290905 0.101560850 5301.929
M_true[170] 0.0009667255 9.704217e05 0.0008089893 0.001120210 4625.569
M_true[171] 0.0008680478 8.382024e05 0.0007257090 0.001002774 3178.937
M_true[172] 0.0009747715 9.358022e05 0.0008231044 0.001120076 4790.759
M_true[173] 0.1173411332 1.179082e02 0.0990117515 0.136820390 4913.305
M_true[174] 0.0481464885 5.088016e03 0.0399478245 0.056162464 4459.045
M_true[175] 0.0764081500 7.633404e03 0.0643051320 0.088423083 5720.210
M_true[176] 0.0866757924 8.740503e03 0.0729003225 0.100575090 5683.205
M_true[177] 0.0530617235 5.169180e03 0.0448247865 0.061317172 3356.499
M_true[178] 0.0569924843 5.634844e03 0.0477599575 0.065910313 4385.082
M_true[179] 0.0898655397 8.890999e03 0.0755746165 0.104114510 5428.557
M_true[180] 0.0474043347 4.704109e03 0.0398315630 0.054710805 4639.041
M_true[181] 0.0258553452 2.651107e03 0.0215936340 0.030070616 4128.754
M_true[182] 0.0265882314 2.706218e03 0.0222069710 0.030818392 5536.600
a 0.4191816405 5.932353e02 0.3262413600 0.513002280 2767.671
b 0.7838324045 1.417906e02 0.7607243800 0.806090665 2770.061
sigma 0.2627569540 1.680301e02 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:
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/model132bf8a4d7.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 autoformat 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/model135ebb0538.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 illconditioned 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/model135ebb0538.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 illconditioned 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/model135ebb0538.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 illconditioned 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/model135ebb0538.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 illconditioned 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/model135ebb0538.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 illconditioned 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/model135ebb0538.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 illconditioned 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/model135ebb0538.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 illconditioned 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/model137bbf99b3.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 illconditioned 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/model137bbf99b3.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 illconditioned 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/model13186b2286.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 illconditioned 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 eightsided 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()
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.