Practice: Chp. 12
Practice: Chp. 12#
source("iplot.R")
suppressPackageStartupMessages(library(rethinking))
12E1. What is the difference between an ordered categorical variable and an unordered one? Define and then give an example of each.
Answer. An ordered categorical variable is one in which the categories are ranked, but the distance between the categories is not necessarily constant. An example of an ordered categorical variable is the Likert scale; an example of an unordered categorical variable is the Blood type of a person.
See also Categorical variable and Ordinal data.
12E2. What kind of link function does an ordered logistic regression employ? How does it differ from an ordinary logit link?
Answer. It employs a cumulative link, which maps the probability of a category or any category of lower rank (the cumulative probability) to the log-cumulative-odds scale. The major difference from the ordinary logit link is the type of probabilties being mapped; the cumulative link maps cumulative probabilities to log-cumulative-odds and the ordinary logit link maps absolute probabilities to log-odds.
12E3. When count data are zero-inflated, using a model that ignores zero-inflation will tend to induce which kind of inferential error?
Answer. A ‘zero’ means that nothing happened; nothing can happen either because the rate of events is low because the process that generates events failed to get started. We may incorrectly infer the rate of events is lower than it is in reality, when in fact it never got started.
12E4. Over-dispersion is common in count data. Give an example of a natural process that might produce over-dispersed counts. Can you also give an example of a process that might produce under-dispersed counts?
Answer. The sex ratios of families actually skew toward either all boys or girls; see Overdispersion.
In contexts where the waiting time for a service is bounded (e.g. a store is only open a short time) but the model still uses a Poisson distribution you may see underdispersed counts because large values are simply not possible. If one queue feeds into another, and the first queue has a limited number of workers, the waiting times for the second queue may be underdispersed because of the homogeneity enforced by the bottleneck in the first queue. “)
display_markdown(” 12M1. At a certain university, employees are annually rated from 1 to 4 on their productivity, with 1 being least productive and 4 most productive. In a certain department at this certain university in a certain year, the numbers of employees receiving each rating were (from 1 to 4): 12, 36, 7, 41. Compute the log cumulative odds of each rating.
Answer. As a list:
freq_prod <- c(12, 36, 7, 41)
pr_k <- freq_prod / sum(freq_prod)
cum_pr_k <- cumsum(pr_k)
lco = list(log_cum_odds = logit(cum_pr_k))
display(lco)
- -1.94591014905531
- 0
- 0.293761118528163
- Inf
12M2. Make a version of Figure 12.5 for the employee ratings data given just above.
iplot(function() {
off <- 0.04
plot(
1:4, cum_pr_k, type="b",
main="Version of Figure 12.5 for Employee Ratings Data", xlab="response", ylab="cumulative proportion",
ylim=c(0,1)
)
# See `simplehist`
for (j in 1:4) lines(c(j,j), c(0,cum_pr_k[j]), lwd=3, col="gray")
for (j in 1:4) lines(c(j,j)+off, c(if (j==1) 0 else cum_pr_k[j-1],cum_pr_k[j]), lwd=3, col=rangi2)
})

12M3. Can you modify the derivation of the zero-inflated Poisson distribution (ZIPoisson) from the chapter to construct a zero-inflated binomial distribution?
Answer. Call \(p_z\) the additional probability of a zero, and \(p_b\) the probability of success associated with the binomial distribution. The likelihood of observing a zero is:
The likelihood of observing a non-zero value \(k\) is: $\( \begin{align} Pr(k | k>0,p_z,n,p_b) & = (1-p_z)Pr(k|n,p_b) \\ & = (1-p_z)\binom{n}{k}p_b^k(1-p_b)^{n-k} \end{align} \)$
To construct a GLM, use a logit link with \(p_z\) and \(p_b\) (in the typical case where \(n\) is known).
12H1. In 2014, a paper was published that was entitled “Female hurricanes are deadlier than male hurricanes.” As the title suggests, the paper claimed that hurricanes with female names have caused greater loss of life, and the explanation given is that people unconsciously rate female hurricanes as less dangerous and so are less likely to evacuate.
Statisticians severely criticized the paper after publication. Here, you’ll explore the complete data used in the paper and consider the hypothesis that hurricanes with female names are deadlier. Load the data with:
R code 12.38
library(rethinking)
data(Hurricanes)
Acquaint yourself with the columns by inspecting the help ?Hurricanes
. In this problem, you’ll
focus on predicting deaths
using femininity
of each hurricane’s name. Fit and interpret the
simplest possible model, a Poisson model of deaths
using femininity
as a predictor. You can use
map
or ulam
. Compare the model to an intercept-only Poisson model of deaths. How strong is the
association between femininity of name and deaths? Which storms does the model fit (retrodict) well?
Which storms does it fit poorly?
Answer. The Hurricanes
data.frame with an index and a standardized femininity:
data(Hurricanes)
hur <- Hurricanes
hur$index = 1:92
hur$StdFem = standardize(hur$femininity)
display(hur)
hur_dat = list(
Deaths = hur$deaths,
Femininity = hur$StdFem
)
m_hur_pois <- map(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a + bFemininity * Femininity,
a ~ dnorm(3, 0.5),
bFemininity ~ dnorm(0, 1.0)
),
data = hur_dat
)
m_hur_intercept <- map(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a,
a ~ dnorm(3, 0.5)
),
data = hur_dat
)
name | year | deaths | category | min_pressure | damage_norm | female | femininity | index | StdFem |
---|---|---|---|---|---|---|---|---|---|
<fct> | <int> | <int> | <int> | <int> | <int> | <int> | <dbl> | <int> | <dbl> |
Easy | 1950 | 2 | 3 | 960 | 1590 | 1 | 6.77778 | 1 | -0.0009347453 |
King | 1950 | 4 | 3 | 955 | 5350 | 0 | 1.38889 | 2 | -1.6707580424 |
Able | 1952 | 3 | 1 | 985 | 150 | 0 | 3.83333 | 3 | -0.9133139565 |
Barbara | 1953 | 1 | 1 | 987 | 58 | 1 | 9.83333 | 4 | 0.9458703621 |
Florence | 1953 | 0 | 1 | 985 | 15 | 1 | 8.33333 | 5 | 0.4810742824 |
Carol | 1954 | 60 | 3 | 960 | 19321 | 1 | 8.11111 | 6 | 0.4122162926 |
Edna | 1954 | 20 | 3 | 954 | 3230 | 1 | 8.55556 | 7 | 0.5499353710 |
Hazel | 1954 | 20 | 4 | 938 | 24260 | 1 | 9.44444 | 8 | 0.8253673305 |
Connie | 1955 | 0 | 3 | 962 | 2030 | 1 | 8.50000 | 9 | 0.5327193242 |
Diane | 1955 | 200 | 1 | 987 | 14730 | 1 | 9.88889 | 10 | 0.9630864089 |
Ione | 1955 | 7 | 3 | 960 | 6200 | 0 | 5.94444 | 11 | -0.2591568553 |
Flossy | 1956 | 15 | 2 | 975 | 1540 | 1 | 7.00000 | 12 | 0.0679232445 |
Helene | 1958 | 1 | 3 | 946 | 540 | 1 | 9.88889 | 13 | 0.9630864089 |
Debra | 1959 | 0 | 1 | 984 | 430 | 1 | 9.88889 | 14 | 0.9630864089 |
Gracie | 1959 | 22 | 3 | 950 | 510 | 1 | 9.77778 | 15 | 0.9286574139 |
Donna | 1960 | 50 | 4 | 930 | 53270 | 1 | 9.27778 | 16 | 0.7737253874 |
Ethel | 1960 | 0 | 1 | 981 | 35 | 1 | 8.72222 | 17 | 0.6015773141 |
Carla | 1961 | 46 | 4 | 931 | 15850 | 1 | 9.50000 | 18 | 0.8425833773 |
Cindy | 1963 | 3 | 1 | 996 | 300 | 1 | 9.94444 | 19 | 0.9802993570 |
Cleo | 1964 | 3 | 2 | 968 | 6450 | 1 | 7.94444 | 20 | 0.3605712508 |
Dora | 1964 | 5 | 2 | 966 | 16260 | 1 | 9.33333 | 21 | 0.7909383355 |
Hilda | 1964 | 37 | 3 | 950 | 2770 | 1 | 8.83333 | 22 | 0.6360063090 |
Isbell | 1964 | 3 | 2 | 974 | 800 | 1 | 9.44444 | 23 | 0.8253673305 |
Betsy | 1965 | 75 | 3 | 948 | 20000 | 1 | 8.33333 | 24 | 0.4810742824 |
Alma | 1966 | 6 | 2 | 982 | 730 | 1 | 8.77778 | 25 | 0.6187933608 |
Inez | 1966 | 3 | 1 | 983 | 99 | 1 | 8.27778 | 26 | 0.4638613343 |
Beulah | 1967 | 15 | 3 | 950 | 5060 | 1 | 7.27778 | 27 | 0.1539972812 |
Gladys | 1968 | 3 | 2 | 977 | 800 | 1 | 8.94444 | 28 | 0.6704353039 |
Camille | 1969 | 256 | 5 | 909 | 23040 | 1 | 9.05556 | 29 | 0.7048673975 |
Celia | 1970 | 22 | 3 | 945 | 6870 | 1 | 9.44444 | 30 | 0.8253673305 |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
Bertha | 1996 | 8 | 2 | 974 | 700 | 1 | 8.50000 | 63 | 0.5327193 |
Fran | 1996 | 26 | 3 | 954 | 8260 | 1 | 7.16667 | 64 | 0.1195683 |
Danny | 1997 | 10 | 1 | 984 | 200 | 0 | 2.22222 | 65 | -1.4125390 |
Bonnie | 1998 | 3 | 2 | 964 | 1650 | 1 | 9.38889 | 66 | 0.8081544 |
Earl | 1998 | 3 | 1 | 987 | 160 | 0 | 1.88889 | 67 | -1.5158260 |
Georges | 1998 | 1 | 2 | 964 | 3870 | 0 | 2.27778 | 68 | -1.3953230 |
Bret | 1999 | 0 | 3 | 951 | 94 | 0 | 2.33333 | 69 | -1.3781100 |
Floyd | 1999 | 56 | 2 | 956 | 8130 | 0 | 1.83333 | 70 | -1.5330421 |
Irene | 1999 | 8 | 1 | 964 | 1430 | 1 | 9.27778 | 71 | 0.7737254 |
Lili | 2002 | 2 | 1 | 963 | 1260 | 1 | 10.33333 | 72 | 1.1008024 |
Claudette | 2003 | 3 | 1 | 979 | 250 | 1 | 9.16667 | 73 | 0.7392964 |
Isabel | 2003 | 51 | 2 | 957 | 4980 | 1 | 9.38889 | 74 | 0.8081544 |
Alex | 2004 | 1 | 1 | 972 | 5 | 0 | 4.16667 | 75 | -0.8100239 |
Charley | 2004 | 10 | 4 | 941 | 20510 | 0 | 2.88889 | 76 | -1.2059620 |
Frances | 2004 | 7 | 2 | 960 | 12620 | 1 | 6.00000 | 77 | -0.2419408 |
Gaston | 2004 | 8 | 1 | 985 | 170 | 0 | 2.66667 | 78 | -1.2748200 |
Ivan | 2004 | 25 | 3 | 946 | 18590 | 0 | 1.05556 | 79 | -1.7740450 |
Jeanne | 2004 | 5 | 3 | 950 | 10210 | 1 | 8.50000 | 80 | 0.5327193 |
Cindy | 2005 | 1 | 1 | 991 | 350 | 1 | 9.94444 | 81 | 0.9802994 |
Dennis | 2005 | 15 | 3 | 946 | 2650 | 0 | 2.44444 | 82 | -1.3436810 |
Ophelia | 2005 | 1 | 1 | 982 | 91 | 1 | 9.16667 | 83 | 0.7392964 |
Rita | 2005 | 62 | 3 | 937 | 10690 | 1 | 9.50000 | 84 | 0.8425834 |
Wilma | 2005 | 5 | 3 | 950 | 25960 | 1 | 8.61111 | 85 | 0.5671483 |
Humberto | 2007 | 1 | 1 | 985 | 51 | 0 | 2.38889 | 86 | -1.3608940 |
Dolly | 2008 | 1 | 1 | 967 | 1110 | 1 | 9.83333 | 87 | 0.9458704 |
Gustav | 2008 | 52 | 2 | 954 | 4360 | 0 | 1.72222 | 88 | -1.5674711 |
Ike | 2008 | 84 | 2 | 950 | 20370 | 0 | 1.88889 | 89 | -1.5158260 |
Irene | 2011 | 41 | 1 | 952 | 7110 | 1 | 9.27778 | 90 | 0.7737254 |
Isaac | 2012 | 5 | 1 | 966 | 24000 | 0 | 1.94444 | 91 | -1.4986131 |
Sandy | 2012 | 159 | 2 | 942 | 75000 | 1 | 9.00000 | 92 | 0.6876514 |
The model with the femininity predictor appears to have inferred a minor but significant association of deaths with femininity:
iplot(function() {
plot(precis(m_hur_pois), main="m_hur_pois")
}, ar=4.5)

The model with the femininity predictor performs only marginally better on WAIC than the intercept-only model:
iplot(function() {
plot(compare(m_hur_pois, m_hur_intercept))
}, ar=4.5)

The model is not accurate when it attempts to retrodict the death toll on storms with many deaths. The predictions it makes are less dispersed than the data, that is, the data are overdispersed relative to the predictions (Overdispersion).
iplot(function() {
result <- postcheck(m_hur_pois, window=92)
display_markdown("The raw data:")
display(result)
})
The raw data:
- $mean
-
- 20.1251659857895
- 13.5272690213414
- 16.1946676460679
- 25.2301287082935
- 22.5782266740179
- 22.2100755222089
- 22.9525690372696
- 24.513735856352
- 22.8583946338462
- 25.3341934023582
- 18.9237585345459
- 20.4583741902435
- 25.3341934023582
- 25.3341934023582
- 25.1265150113655
- 24.2130552851921
- 23.2374194592045
- 24.6148119744012
- 25.43867345702
- 21.9379438722294
- 24.3128593461609
- 23.4293096819226
- 24.513735856352
- 22.5782266740179
- 23.3331736632992
- 22.4856224225171
- 20.8827526901423
- 23.6228029801984
- 23.8179305320726
- 24.513735856352
- 22.8583946338462
- 21.0549844953731
- 25.5436083208327
- 23.142079863412
- 23.2374194592045
- 23.6228029801984
- 26.3989288665572
- 20.2910782055103
- 13.8065131719318
- 13.8630519306067
- 14.6797282101263
- 14.8609723008321
- 25.2301287082935
- 25.43867345702
- 13.8065131719318
- 14.3826130109495
- 25.0233131904701
- 24.6148119744012
- 14.0915799092339
- 24.920558586997
- 24.4130976457744
- 15.1061505280087
- 13.9768434241141
- 22.5782266740179
- 23.8179305320726
- 15.1061505280087
- 14.5007199645343
- 13.8065131719318
- 14.3826130109495
- 25.2301287082935
- 20.7971642453911
- 22.8583946338462
- 22.8583946338462
- 20.7119459426219
- 14.3826130109495
- 24.4130976457744
- 14.0340982113646
- 14.4415496593931
- 14.5007199645343
- 13.9768434241141
- 24.2130552851921
- 26.1824197235577
- 24.0146708789476
- 24.4130976457744
- 16.5975275516719
- 15.1061505280087
- 19.0015518427568
- 14.8609723008321
- 13.1997116297174
- 22.8583946338462
- 25.43867345702
- 14.6198081648425
- 24.0146708789476
- 24.6148119744012
- 23.0471188819599
- 14.5601462541308
- 25.2301287082935
- 13.8630519306067
- 14.0340982113646
- 24.2130552851921
- 14.0915799092339
- 23.7201637909248
- $PI
A matrix: 2 × 92 of type dbl 5% 19.38088 12.35608 15.20263 24.07366 21.69864 21.37226 22.03918 23.46239 21.95723 24.15860 ⋯ 23.00099 23.54719 22.13124 13.45796 24.07366 12.71296 12.89525 23.17868 12.95659 22.72140 94% 20.90360 14.73865 17.17299 26.42291 23.45277 23.05732 23.87516 25.60047 23.76346 26.54061 ⋯ 25.04226 25.71440 23.98079 15.66709 26.42291 15.03912 15.19748 25.25856 15.24904 24.71972 - $outPI
A matrix: 2 × 92 of type dbl 5% 13 8 10 18 15 15 16 17 16 17.000 ⋯ 16 17 15.945 9 17.945 8 8 17 8 16 94% 28 20 23 34 31 30 31 32 31 33.055 ⋯ 32 33 30.000 21 33.000 20 20 32 21 32

12H2. Counts are nearly always over-dispersed relative to Poisson. So fit a gamma-Poisson (aka
negative-binomial) model to predict deaths
using femininity
. Show that the over-dispersed model
no longer shows as precise a positive association between femininity and deaths, with an 89% interval
that overlaps zero. Can you explain why the association diminished in strength?
data(Hurricanes)
hur <- Hurricanes
hur$index = 1:92
hur$StdFem = standardize(hur$femininity)
hur_dat = list(
Deaths = hur$deaths,
Femininity = hur$StdFem
)
m_hur_gamma_pois <- map(
alist(
Deaths ~ dgampois(lambda, phi),
log(lambda) <- a + bFemininity * Femininity,
a ~ dnorm(3, 0.5),
bFemininity ~ dnorm(0, 1.0),
phi ~ dexp(1)
),
data = hur_dat
)
Answer. The model no longer infers an association of deaths with femininity:
iplot(function() {
plot(precis(m_hur_gamma_pois), main="m_hur_gamma_pois")
}, ar=4.5)

The association has diminished in strength because the influence of outliers are diminished. Notice the storms with the greatest death toll are female in the original data: Diane, Camille, and Sandy.
12H3. In the data, there are two measures of a hurricane’s potential to cause death:
damage_norm
and min_pressure
. Consult ?Hurricanes
for their meanings. It makes some sense to
imagine that femininity of a name matters more when the hurricane is itself deadly. This implies an
interaction between femininity and either or both of damage_norm
and min_pressure
. Fit a series
of models evaluating these interactions. Interpret and compare the models. In interpreting the
estimates, it may help to generate counterfactual predictions contrasting hurricanes with masculine
and feminine names. Are the effect sizes plausible?
data(Hurricanes)
hur <- Hurricanes
hur$index = 1:92
hur$StdFem = standardize(hur$femininity)
hur$StdDam = standardize(hur$damage_norm)
hur$StdPres = standardize(hur$min_pressure)
Answer. To get some unreasonable answers, lets go back to the Poisson distribution rather than the Gamma-Poisson. It’s likely this is what the authors of the original paper were using.
hur_dat = list(
Deaths = hur$deaths,
StdFem = hur$StdFem,
StdDam = hur$StdDam
)
m_fem_dam_inter <- quap(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a + bFem * StdFem + bDam * StdDam + bFemDamInter * StdFem * StdDam,
a ~ dnorm(3, 0.5),
bFem ~ dnorm(0, 1.0),
bDam ~ dnorm(0, 1.0),
bFemDamInter ~ dnorm(0, 0.4)
),
data = hur_dat
)
In the first model, lets assume an interaction between femininity
and damage_norm
. Damage is
unsurprisingly associated with more deaths. The model also believes there is a minor interaction
between damage and femininity, where together they increase deaths more than either alone.
iplot(function() {
plot(precis(m_fem_dam_inter), main="m_fem_dam_inter")
}, ar=4.5)
hur_dat = list(
Deaths = hur$deaths,
StdFem = hur$StdFem,
StdPres = hur$StdPres
)
m_fem_pres_inter <- quap(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a + bFem * StdFem + bPres * StdPres + bFemPresInter * StdFem * StdPres,
a ~ dnorm(3, 0.5),
bFem ~ dnorm(0, 1.0),
bPres ~ dnorm(0, 1.0),
bFemPresInter ~ dnorm(0, 0.4)
),
data = hur_dat
)

In the second model, lets assume an interaction between femininity
and min_pressure
. Pressure is
unsurprisingly negatively associated with more deaths; lower pressure indicates a stronger storm.
There also appears to be a minor interaction between pressure and femininity, though now the
association is in the opposite direction. That is, when both pressure and femininity are positive,
which means the storm is weak and feminine, more deaths occur.
iplot(function() {
plot(precis(m_fem_pres_inter), main="m_fem_pres_inter")
}, ar=4.5)
hur_dat = list(
Deaths = hur$deaths,
StdFem = hur$StdFem,
StdPres = hur$StdPres,
StdDam = hur$StdDam
)
m_fem_dam_pres_inter <- quap(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a + bFem * StdFem + bPres * StdPres + bDam * StdDam +
bFemPresInter * StdFem * StdPres +
bFemDamInter * StdFem * StdDam +
bDamPresInter * StdDam * StdPres +
bFDPInter * StdPres * StdDam * StdPres,
a ~ dnorm(3, 0.5),
bFem ~ dnorm(0, 1.0),
bDam ~ dnorm(0, 1.0),
bPres ~ dnorm(0, 1.0),
bFemPresInter ~ dnorm(0, 0.4),
bFemDamInter ~ dnorm(0, 0.4),
bDamPresInter ~ dnorm(0, 0.4),
bFDPInter ~ dnorm(0, 0.4)
),
data = hur_dat
)

In the third model, lets assume interactions between all three predictors. The influence of femininity has dropped significantly. The three-way interaction term is hard to interpret on the parameter scale. Lets try to interpret it on the outcome scale.
iplot(function() {
plot(precis(m_fem_dam_pres_inter), main="m_fem_dam_pres_inter")
}, ar=3.5)

To some extent this is a counterfactual plot rather than just a posterior prediction plot because we’re considering extremes of at least some predictors. Still, many of the predictors vary here only within the range of the data.
In the following Enneaptych, we hold damage and pressure at typical values and vary femininity. We observe the effect on deaths. Many of the effect sizes are inplausible. In the bottom-right plot, notice an implication of the non-negative three-way interaction parameter: when pressure is high and damage is low femininity has a dramatic effect on deaths.
iplot(function() {
par(mfrow=c(3,3))
for (d in -1: 1) {
for (p in -1:1) {
sim_dat = data.frame(StdPres=p, StdDam=d, StdFem=seq(from=-2, to=2, length.out=30))
s <- sim(m_fem_dam_pres_inter, data=sim_dat)
plot(sim_dat$StdFem, colMeans(s),
main=paste("StdPres = ", p, ",", "StdDam = ", d),
type="l", ylab="Deaths", xlab="StdFem")
}
}
})

Ignoring the known outliers, PSIS believes the full interaction model performs best:
iplot(function() {
plot(compare(m_fem_dam_pres_inter, m_fem_pres_inter, m_fem_dam_inter, func=PSIS))
}, ar=4.5)
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.

12H4. In the original hurricanes paper, storm damage (damage_norm
) was used directly. This
assumption implies that mortality increases exponentially with a linear increase in storm strength,
because a Poisson regression uses a log link. So it’s worth exploring an alternative hypothesis:
that the logarithm of storm strength is what matters. Explore this by using the logarithm of
damage_norm
as a predictor. Using the best model structure from the previous problem, compare a
model that uses log(damage_norm)
to a model that uses damage_norm
directly. Compare their
PSIS/WAIC values as well as their implied predictions. What do you conclude?
Answer We’ll use the full interaction model as our baseline, though it probably doesn’t matter which model we start from.
data(Hurricanes)
hur <- Hurricanes
hur$index = 1:92
hur$StdFem = standardize(hur$femininity)
hur$StdDam = standardize(hur$damage_norm)
hur$StdLogDam = standardize(log(hur$damage_norm))
hur$StdPres = standardize(hur$min_pressure)
hur_dat = list(
Deaths = hur$deaths,
StdFem = hur$StdFem,
StdPres = hur$StdPres,
StdDam = hur$StdDam
)
m_fem_dam_pres_inter <- quap(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a + bFem * StdFem + bPres * StdPres + bDam * StdDam +
bFemPresInter * StdFem * StdPres +
bFemDamInter * StdFem * StdDam +
bDamPresInter * StdDam * StdPres +
bFDPInter * StdPres * StdDam * StdPres,
a ~ dnorm(3, 0.5),
bFem ~ dnorm(0, 1.0),
bDam ~ dnorm(0, 1.0),
bPres ~ dnorm(0, 1.0),
bFemPresInter ~ dnorm(0, 0.4),
bFemDamInter ~ dnorm(0, 0.4),
bDamPresInter ~ dnorm(0, 0.4),
bFDPInter ~ dnorm(0, 0.4)
),
data = hur_dat
)
First we’ll run postcheck
on the baseline model. Notice there is more dispersion than in the plain
old Poisson model in 12H1:
iplot(function() {
result <- postcheck(m_fem_dam_pres_inter, window=92)
display_markdown("The raw data:")
display(result)
})
hur_dat = list(
Deaths = hur$deaths,
StdFem = hur$StdFem,
StdPres = hur$StdPres,
StdLogDam = hur$StdLogDam
)
m_fem_logdam_pres_inter <- quap(
alist(
Deaths ~ dpois(lambda),
log(lambda) <- a + bFem * StdFem + bPres * StdPres + bDam * StdLogDam +
bFemPresInter * StdFem * StdPres +
bFemDamInter * StdFem * StdLogDam +
bDamPresInter * StdLogDam * StdPres +
bFDPInter * StdPres * StdLogDam * StdPres,
a ~ dnorm(3, 0.5),
bFem ~ dnorm(0, 1.0),
bDam ~ dnorm(0, 1.0),
bPres ~ dnorm(0, 1.0),
bFemPresInter ~ dnorm(0, 0.4),
bFemDamInter ~ dnorm(0, 0.4),
bDamPresInter ~ dnorm(0, 0.4),
bFDPInter ~ dnorm(0, 0.4)
),
data = hur_dat
)
The raw data:
- $mean
-
- 12.4289629811358
- 18.5555907721241
- 2.1249070586138
- 1.77374475898114
- 2.06701979381536
- 39.5683369755651
- 15.9939793092382
- 37.4206334005963
- 11.3393368856691
- 96.6869545199769
- 16.2813093468006
- 5.70521738475244
- 16.91337602597
- 2.48568771306356
- 15.2596820401918
- 114.838235470021
- 2.86766572841986
- 34.5128187022824
- 0.856829259425056
- 14.1797118534941
- 42.1928980075304
- 17.3290292624905
- 5.31717940816366
- 30.8172826315106
- 3.08573684193247
- 2.48545947274957
- 19.4141586210783
- 4.41744064025897
- 203.106284084467
- 21.5776644407842
- 3.78105321641471
- 3.67383539337091
- 0.908907885314257
- 184.984153821049
- 15.4361111412495
- 18.0008369177879
- 3.43397244991166
- 0.824230228928947
- 1.90433426037738
- 8.89528698983109
- 22.0565207809061
- 25.0936386773702
- 22.8871246559745
- 15.5198373640191
- 0.319483073147672
- 1.76904239289811
- 14.4654254906871
- 20.1824722641769
- 10.122581589849
- 8.30801811862062
- 1.34098254166381
- 1.29640697342803
- 0.933284011415338
- 2.24111365545684
- 2.04205676669216
- 26.1365654250733
- 2.58118066040024
- 13.7981829597164
- 66.615328996312
- 9.78792501462256
- 6.50495949367915
- 23.2557773959279
- 5.28556240933658
- 20.0449904479539
- 2.34938006223169
- 9.92009766804856
- 1.76765633021552
- 12.7572493066998
- 20.724355786377
- 18.5368837752248
- 9.7380836695525
- 9.7864089263809
- 3.4710114383841
- 16.0821194776559
- 5.87258812291837
- 21.9002218721581
- 23.0819441177566
- 2.13353996480914
- 19.1758534832144
- 22.7315653259822
- 1.37215259835675
- 24.4221358194378
- 2.6812724963672
- 26.6361393883323
- 39.9601777458033
- 2.08160460381955
- 8.09843833341641
- 18.9855117772041
- 19.8208111691232
- 19.7288010624798
- 39.4585510163461
- 125.255203562995
- $PI
A matrix: 2 × 92 of type dbl 5% 11.68685 16.65435 1.762939 1.514754 1.797806 37.57049 15.13025 35.14764 10.62842 87.18363 ⋯ 2.347914 24.56441 37.48099 1.658119 7.408450 17.03529 17.91018 18.57943 32.84712 111.1257 94% 13.18583 20.65130 2.545884 2.069954 2.358545 41.74312 16.90335 39.78237 12.05322 106.80668 ⋯ 3.041965 28.82114 42.48130 2.576003 8.807731 21.13936 21.89069 20.89448 46.56996 141.6279 - $outPI
A matrix: 2 × 92 of type dbl 5% 7 12 0 0 0 29 10 27 6 79 ⋯ 1 19.000 30 0 4 12 13 13 27 102 94% 19 26 5 4 5 50 23 48 17 115 ⋯ 6 35.055 51 5 13 26 27 27 52 149

Ignoring the known outliers, comparing the models using PSIS, it appears the log(damage)
model
improves performance:
iplot(function() {
plot(compare(m_fem_dam_pres_inter, m_fem_logdam_pres_inter, func=PSIS))
}, ar=4.5)
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.
Some Pareto k values are very high (>1). Set pointwise=TRUE to inspect individual points.

Run postcheck
to visualize how much the model has changed on the outcome scale:
iplot(function() {
result <- postcheck(m_fem_logdam_pres_inter, window=92)
display_markdown("The raw data:")
display(result)
})
The raw data:
- $mean
-
- 10.9140498819251
- 16.3313710022569
- 2.6709357607139
- 0.322526496273671
- 0.136160362124281
- 40.9683395319938
- 14.2580819168725
- 49.8113315443322
- 11.9564068702583
- 154.596179681551
- 20.5892701945882
- 12.2153330272158
- 4.49122488660067
- 3.18140578723594
- 4.53948056380432
- 113.659356226795
- 0.398583050207092
- 45.2021587791377
- 1.3848939343451
- 26.8510567912674
- 49.2758181095555
- 12.6829339090324
- 6.49005192517415
- 36.1190315021328
- 6.09839407651378
- 0.985349989263909
- 17.2100121657768
- 6.68745222404144
- 226.247559918251
- 21.1517893910113
- 3.06888928147871
- 5.14531160281342
- 0.824883825194171
- 116.523475702912
- 8.0641871709459
- 20.7742255615973
- 4.26613646512088
- 0.336503528493622
- 2.15921940466006
- 17.6042019668685
- 19.471387174271
- 12.4919902588978
- 33.5125543173247
- 3.88133868051414
- 0.980499511312163
- 3.27154926980215
- 17.603607823164
- 12.858401141245
- 21.9762686587169
- 9.03992355758047
- 0.0191863016958919
- 1.05680367616152
- 0.028132573873713
- 0.0431267785317184
- 2.19936585828156
- 27.3093114353624
- 4.90419915633675
- 16.6633629041742
- 82.4226128407474
- 1.60139950190683
- 12.4178519108521
- 21.783740514807
- 6.30210027489317
- 22.3899535257731
- 4.3845208266936
- 10.4265057085354
- 3.44908259779698
- 17.5834546186034
- 6.55064690808527
- 18.2385297667872
- 9.5489541419671
- 8.2831244154469
- 2.34924285804395
- 18.9706169543237
- 0.800492465445044
- 23.6536111949234
- 28.2980353485363
- 3.51696407225748
- 18.9235904552031
- 25.6757596927057
- 2.0637764086095
- 13.2976653987726
- 0.821268377189356
- 29.7819961636452
- 42.9099807127748
- 1.61784913313187
- 7.9762091532931
- 15.6121558301295
- 20.7781808862812
- 22.0385062348842
- 31.5885609198836
- 86.8549157202277
- $PI
A matrix: 2 × 92 of type dbl 5% 10.13220 14.85667 2.089515 0.2280939 0.09321956 38.65071 13.30539 46.44437 11.14502 139.1676 ⋯ 0.6348161 27.66339 40.25265 1.144068 7.153467 14.10175 18.75725 20.62962 26.19048 78.75532 94% 11.70628 17.87311 3.359151 0.4401200 0.19064352 43.36988 15.28017 53.01502 12.81762 170.1517 ⋯ 1.0467692 31.85393 45.78290 2.201723 8.841428 17.14417 22.94061 23.47136 37.87033 95.57261 - $outPI
A matrix: 2 × 92 of type dbl 5% 6 10 0 0 0 31 8 38 7 128 ⋯ 0 21 33 0 4 9 13 15 21 71 94% 16 23 6 1 1 52 21 61 18 180 ⋯ 2 39 54 4 13 22 28 30 44 105

12H5. One hypothesis from developmental psychology, usually attributed to Carol Gilligan,
proposes that women and men have different average tendencies in moral reasoning. Like most
hypotheses in social psychology, it is merely descriptive, not causal. The notion is that women are
more concerned with care (avoiding harm), while men are more concerned with justice and rights.
Evaluate this hypothesis, using the Trolley
data, supposing that contact
provides a proxy for
physical harm. Are women more or less bothered by contact than are men, in these data? Figure out
the model(s) that is needed to address this question.
Answer. The head
of the Trolley
data.frame:
data(Trolley)
d <- Trolley
display(head(Trolley))
case | response | order | id | age | male | edu | action | intention | contact | story | action2 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
<fct> | <int> | <int> | <fct> | <int> | <int> | <fct> | <int> | <int> | <int> | <fct> | <int> | |
1 | cfaqu | 4 | 2 | 96;434 | 14 | 0 | Middle School | 0 | 0 | 1 | aqu | 1 |
2 | cfbur | 3 | 31 | 96;434 | 14 | 0 | Middle School | 0 | 0 | 1 | bur | 1 |
3 | cfrub | 4 | 16 | 96;434 | 14 | 0 | Middle School | 0 | 0 | 1 | rub | 1 |
4 | cibox | 3 | 32 | 96;434 | 14 | 0 | Middle School | 0 | 1 | 1 | box | 1 |
5 | cibur | 3 | 4 | 96;434 | 14 | 0 | Middle School | 0 | 1 | 1 | bur | 1 |
6 | cispe | 3 | 9 | 96;434 | 14 | 0 | Middle School | 0 | 1 | 1 | spe | 1 |
We’ll start from model m12.5
in the text, converting bC
and bIC
to index variables based on
gender. Before doing the conversion, lets reproduce the precis
summary for m12.5
in the text:
dat <- list(
R = d$response,
A = d$action,
I = d$intention,
C = d$contact
)
m12.5 <- ulam(
alist(
R ~ dordlogit(phi, cutpoints),
phi <- bA * A + bC * C + BI * I,
BI <- bI + bIA * A + bIC * C,
c(bA, bI, bC, bIA, bIC) ~ dnorm(0, 0.5),
cutpoints ~ dnorm(0, 1.5)
),
data = dat, chains = 4, cores = min(detectCores(), 4), log_lik=TRUE
)
iplot(function() {
plot(precis(m12.5), main="m12.5")
}, ar=4)
dat <- list(
R = d$response,
A = d$action,
I = d$intention,
C = d$contact,
gid = ifelse(d$male == 1, 2, 1)
)
m_trolley_care <- ulam(
alist(
R ~ dordlogit(phi, cutpoints),
phi <- bA * A + bC[gid] * C + BI * I,
BI <- bI + bIA * A + bIC[gid] * C,
c(bA, bI, bIA) ~ dnorm(0, 0.5),
bC[gid] ~ dnorm(0, 0.5),
bIC[gid] ~ dnorm(0, 0.5),
cutpoints ~ dnorm(0, 1.5)
),
data = dat, chains = 4, cores = min(detectCores(), 4), log_lik=TRUE
)
flush.console()
6 vector or matrix parameters hidden. Use depth=2 to show them.

The following plot shows the same precis
summary after converting bC
and bIC
to index
variables based on gender. It appears women disapprove of contact more than men from this data and
model, lending some support to the theory:
iplot(function() {
plot(precis(m_trolley_care, depth=2), main="m_trolley_care")
}, ar=2.5)

Comparing the models using PSIS, to check for overfitting or outliers:
iplot(function() {
plot(compare(m12.5, m_trolley_care, func=PSIS))
}, ar=4.5)

12H6. The data in data(Fish)
are records of visits to a national park. See ?Fish
for
details. The question of interest is how many fish an average visitor takes per hour, when fishing.
The problem is that not everyone tried to fish, so the fish_caught
numbers are zero-inflated. As
with the monks example in the chapter, there is a process that determines who is fishing (working)
and another process that determines fish per hour (manuscripts per day), conditional on fishing
(working). We want to model both. Otherwise we’ll end up with an underestimate of rate of fish
extraction from the park.
You will model these data using zero-inflated Poisson GLMs. Predict fish_caught
as a function
of any of the other variables you think are relevant. One thing you must do, however, is use a
proper Poisson offset/exposure in the Poisson portion of the zero-inflated model. Then use the
hours
variable to construct the offset. This will adjust the model for the differing amount of
time individuals spent in the park.
Answer. The short documentation on Fish
:
display(help(Fish))
data(Fish)
Fish {rethinking} | R Documentation |
Fishing data
Description
Fishing data from visitors to a national park.
Usage
data(Fish)
Format
fish_caught : Number of fish caught during visit
livebait : Whether or not group used livebait to fish
camper : Whether or not group had a camper
persons : Number of adults in group
child : Number of children in group
hours : Number of hours group spent in park
The head
of Fish
:
display(head(Fish))
fish_caught | livebait | camper | persons | child | hours | |
---|---|---|---|---|---|---|
<int> | <int> | <int> | <int> | <int> | <dbl> | |
1 | 0 | 0 | 0 | 1 | 0 | 21.124 |
2 | 0 | 1 | 1 | 1 | 0 | 5.732 |
3 | 0 | 1 | 0 | 1 | 0 | 1.323 |
4 | 0 | 1 | 1 | 2 | 1 | 0.548 |
5 | 1 | 1 | 0 | 1 | 0 | 1.695 |
6 | 0 | 1 | 1 | 4 | 2 | 0.493 |
Notice several outliers in the full dataframe (not represented in head
):
display(summary(Fish))
fish_caught livebait camper persons
Min. : 0.000 Min. :0.000 Min. :0.000 Min. :1.000
1st Qu.: 0.000 1st Qu.:1.000 1st Qu.:0.000 1st Qu.:2.000
Median : 0.000 Median :1.000 Median :1.000 Median :2.000
Mean : 3.296 Mean :0.864 Mean :0.588 Mean :2.528
3rd Qu.: 2.000 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:4.000
Max. :149.000 Max. :1.000 Max. :1.000 Max. :4.000
child hours
Min. :0.000 Min. : 0.0040
1st Qu.:0.000 1st Qu.: 0.2865
Median :0.000 Median : 1.8315
Mean :0.684 Mean : 5.5260
3rd Qu.:1.000 3rd Qu.: 7.3427
Max. :3.000 Max. :71.0360
One complication introduced by the question and data is how to infer fish per ‘average visitor’ when
every observation includes several visitors, both adults and children. This solution disregards
children, assuming any fishing they do is offset by the time they take away from fishing. Addressing
persons
can be done the same way we address hours
, by treating it as an exposure. That is:
$\(
\begin{align}
log \left(\frac{fish}{person \cdot hour}\right) & = log(fish) - log(person \cdot hour) \\
& = log(fish) - log(person) - log(hour)
\end{align}
\)$
fish_dat <- list(
FishCaught = Fish$fish_caught,
LogPersonsInGroup = log(Fish$persons),
LogHours = log(Fish$hours)
)
m_fish <- ulam(
alist(
FishCaught ~ dzipois(p, fish_per_person_hour),
log(fish_per_person_hour) <- LogPersonsInGroup + LogHours + al,
logit(p) <- ap,
ap ~ dnorm(-0.5, 1),
al ~ dnorm(1, 2)
), data = fish_dat, chains = 4
)
flush.console()
SAMPLING FOR MODEL '1a5d1436509eaa3c57b3600d89299fd9' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000111 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.11 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
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:
Chain 1: Elapsed Time: 0.360352 seconds (Warm-up)
Chain 1: 0.28316 seconds (Sampling)
Chain 1: 0.643512 seconds (Total)
Chain 1:
SAMPLING FOR MODEL '1a5d1436509eaa3c57b3600d89299fd9' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 8.9e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.89 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
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:
Chain 2: Elapsed Time: 0.388114 seconds (Warm-up)
Chain 2: 0.290798 seconds (Sampling)
Chain 2: 0.678912 seconds (Total)
Chain 2:
SAMPLING FOR MODEL '1a5d1436509eaa3c57b3600d89299fd9' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.00012 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
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:
Chain 3: Elapsed Time: 0.353938 seconds (Warm-up)
Chain 3: 0.310407 seconds (Sampling)
Chain 3: 0.664345 seconds (Total)
Chain 3:
SAMPLING FOR MODEL '1a5d1436509eaa3c57b3600d89299fd9' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 8.8e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.88 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
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:
Chain 4: Elapsed Time: 0.422965 seconds (Warm-up)
Chain 4: 0.313815 seconds (Sampling)
Chain 4: 0.73678 seconds (Total)
Chain 4:
The precision of the parameters, which is rather hard to interpret. We can take from ap
that most
visitors who were a part of this survey chose to fish. That is, the probability of not fishing looks
solidly below 0.5, since ap
is solidly negative and any score less than zero implies a probability
of less than 0.5:
iplot(function() {
plot(precis(m_fish))
}, ar=4.5)

Finally, a posterior check. Despite the large uncertainty in ap
the model is still surprised by
outliers. An obvious alternative based on this chapter is a Gamma-Poisson, but this solution doesn’t
pursue that.
iplot(function() {
result <- postcheck(m_fish, window=250)
display_markdown("The raw data:")
display(result)
})
The raw data:
- $mean
-
- 6.76595490376015
- 1.83594269590765
- 0.423752998375056
- 0.351045567814861
- 0.542903501319516
- 0.631625784426009
- 0.0317094080416709
- 0.0051247528148155
- 0.208513380152805
- 1.29111741228258
- 0.185772289537062
- 0.080714856833344
- 4.71797556013951
- 2.20236252216696
- 5.37682659389423
- 2.77281156986611
- 0.0589346573703781
- 0.136446543694462
- 0.531693104537107
- 0.853271343666779
- 1.28375058011128
- 8.50452729618631
- 0.477242605879693
- 4.42138049098206
- 4.10204433120887
- 10.7197017003903
- 6.714707375612
- 0.221004965138918
- 0.34720200320375
- 0.010249505629631
- 0.805226786027884
- 0.494858943680621
- 42.7839988744871
- 1.52845752701872
- 0.693122818203795
- 3.35799428190785
- 0.336952497574119
- 4.44700425505614
- 0.0989717887361242
- 0.293392098648187
- 6.56640984103327
- 13.5101296080573
- 0.0442009930277836
- 18.6131022234099
- 0.0896831742592711
- 2.7036274068661
- 0.868005008009374
- 1.51372386267613
- 1.19919215866683
- 0.0819960450370479
- 0.851349561361224
- 0.315812892213005
- 5.83196870326003
- 5.12539340891734
- 13.3217949421129
- 20.9282093075028
- 2.67544126638461
- 1.23730750772702
- 4.52259435907467
- 1.63319466267151
- 1.89295557097247
- 5.20610826575069
- 0.595111920620449
- 16.3748664315392
- 0.0192178230555581
- 1.58771248144002
- 0.973703034814943
- 0.060536142625008
- 2.39422045567161
- 0.177444566212986
- 9.28893477390401
- 0.00896831742592712
- 0.0858396096481595
- 0.0384356461111162
- 0.798500547958439
- 0.416386166203759
- 0.0413183195694499
- 8.40587580450111
- 0.504788152259326
- 5.16446964913031
- 0.850708967259372
- 0.465071317944506
- 13.6561850632796
- 0.409659928134313
- 1.78661695006505
- 0.555395086305629
- 0.0999326798889021
- 3.07837495644948
- 45.6000505462282
- 0.060536142625008
- 0.650523310430642
- 0.700169353324167
- 2.539635316792
- 0.655968360296383
- 0.226770312055585
- 0.214278727069473
- 0.0166554466481503
- 0.622657467000082
- 12.3993394354461
- 3.02872891355596
- 2.51465214681978
- 6.17724892415822
- 0.161429713666688
- 0.859997581736225
- 0.456103000518579
- 0.0884019860555672
- 0.654687172092679
- 0.140290108305574
- 0.185772289537062
- 4.99791518264881
- 7.68744951927416
- 5.36113203839886
- 0.050927231097229
- 2.66615265190776
- 0.108580700263903
- 2.33624668945401
- 0.497741617138955
- 0.594471326518597
- 0.411261413388943
- 0.373786658430605
- 0.126837632166683
- 3.95919184649589
- 0.0768712922222324
- 0.573972315259335
- 0.0666217865926014
- 0.609845584963044
- 3.08542149156985
- 30.2286747752401
- 2.34393381867624
- 91.0104852383084
- 20.4650597718638
- 2.47781798596329
- 15.8726406556873
- 0.176803972111134
- 2.28051500259289
- 2.92239029264853
- 1.3221862262224
- 1.93811745515303
- 0.545145580675998
- 6.6826776705194
- 0.0140930702407426
- 0.0967297093796424
- 1.34909117850018
- 0.207552489000027
- 1.18221641496775
- 0.334390121166711
- 45.3707178577652
- 14.177948959238
- 0.122033176402794
- 0.737964405333431
- 2.54796304011608
- 0.397168343148201
- 0.373146064328753
- 0.393324778537089
- 1.40290108305574
- 1.62454664229651
- 4.80701814029693
- 0.876332731333449
- 3.12289624652819
- 5.00303993546362
- 4.21350770493111
- 0.38531735226394
- 0.0374747549583383
- 0.26136239355559
- 0.0550910927592665
- 3.02104178433373
- 1.18509908842608
- 11.6613750301127
- 0.0970500064305684
- 0.0243425758703736
- 13.0191142289878
- 2.98004376181521
- 7.96258468601957
- 7.05165987318612
- 0.070465351203713
- 15.3941168616039
- 3.72633589047271
- 1.03936393025477
- 7.56669753107507
- 0.453220327060245
- 9.66079965002906
- 0.27801784020374
- 12.1347740713812
- 0.427916860037094
- 0.199865359777804
- 3.31315269477822
- 0.029787625736115
- 0.171999516347245
- 41.2619472884869
- 1.16588126537052
- 12.66358450246
- 1.63479614792614
- 2.22478331573178
- 0.370583687921345
- 12.5261770676128
- 2.32151302511142
- 0.872168869671411
- 0.347842597305602
- 0.117869314740756
- 2.34329322457438
- 0.0800742627314921
- 6.3828796308527
- 2.90669573715316
- 3.08189822400967
- 0.538099045555627
- 0.30268071312504
- 4.20870324916722
- 0.0128118820370387
- 3.76893539824587
- 0.414144086847277
- 5.99211722872301
- 11.3346720381682
- 1.30777285893073
- 5.90371524266745
- 0.0115306938333349
- 0.207552489000027
- 0.0525287163518588
- 0.0643797072361196
- 2.27410906157437
- 0.144774267018538
- 0.352326756018565
- 0.965695608541794
- 2.78594374895407
- 0.834694114713073
- 0.329905962453747
- 11.2059126236959
- 20.0787815284471
- 0.590307464856559
- 0.091284659513901
- 2.53354967282441
- 2.72508730927814
- 0.0819960450370479
- 0.161429713666688
- 1.74497833344468
- 0.411902007490795
- 8.18935499807516
- 1.65465456508355
- 0.685435688981572
- 8.32420005651499
- 1.62006248358355
- 22.7884945792808
- 0.343358438592638
- 1.84491101333358
- 0.461227753333394
- 3.58124132640325
- 6.54623112682494
- 0.0051247528148155
- 1.26132978654646
- 0.0720668364583429
- 0.149258425731501
- $PI
A matrix: 2 × 250 of type dbl 5% 6.381444 1.731606 0.3996710 0.3310956 0.5120502 0.5957304 0.02990736 0.004833512 0.1966635 1.217743 ⋯ 21.49342 0.3238453 1.740064 0.4350161 3.377719 6.174208 0.004833512 1.189648 0.06797126 0.1407760 94% 7.149916 1.940131 0.4478006 0.3709671 0.5737127 0.6674699 0.03350889 0.005415577 0.2203463 1.364387 ⋯ 24.08172 0.3628437 1.949608 0.4874020 3.784473 6.917723 0.005415577 1.332909 0.07615656 0.1577287 - $outPI
A matrix: 2 × 250 of type dbl 5% 0 0 0 0 0 0 0 0 0 0 ⋯ 0 0 0 0 0 0 0 0 0 0 94% 11 4 1 1 2 2 0 0 1 3 ⋯ 30 1 4 1 7 10 0 3 0 1

12H7. In the trolley data — data(Trolley)
— we saw how education level (modeled as an ordered
category) is associated with responses. But is this association causal? One plausible confound is
that education is also associated with age, through a causal process: People are older when they
finish school than when they begin it. Reconsider the Trolley data in this light. Draw a DAG that
represents hypothetical causal relationships among response, education, and age. Which statical
model or models do you need to evaluate the causal influence of education on responses? Fit these
models to the trolley data. What do you conclude about the causal relationships among these three
variables?
ERROR:
Which statical [sic] model or models do you need
Answer. Consider the following DAGs. Implicit in this list of DAGs is that the response would
never influence age or education. Because we already know that education is associated with the
response, we also only include DAGs where Education
and Response
are d-connected.
A more difficult question is what arrows to allow between education and age. Getting older does not
necessarily lead to more education, but education always leads to greater age. That is, someone
could stop schooling at any point in their lives. This perspective implies Education
should cause
Age
, because the choice to get a particular education is also a choice to get older before
continuing with ‘real’ work.
The issue with allowing this arrow is it contradicts the following definition of Causal dependence from David Lewis:
An event E causally depends on C if, and only if, (i) if C had occurred, then E would have occurred, and (ii) if C had not occurred, then E would not have occurred.
Said another way, causality is almost always associated with the advancement of time. This is
inherent in the definition of causality above, and in lists (see #4) such as the Bradford-Hill
criteria. See also the definition of Causal inference. The Age
covariate is confusingly
a measurement of time. It is perhaps best to think of Age
as a measurement of the potential to get
education; greater Age
indicates an individual is ‘richer’ in terms of time opportunity to get an
education.
Said yet another way, age is not controllable. We will get older regardless of whether we get an education. Education on the other hand is clearly within our control. The issue seems to be that we require two things to get a particular educational degree: a personal commitment to do so and an investment of our time. We talk about investing our time when in the end we are always forced to invest it in something; we can’t choose to stay in a frozen present.
Said yet another way, we should not confuse prediction and causal inference. Someone’s education may help us predict their age, but this does not imply the relationship is causal.
Described as a rule for DAGs, indirect measures of time (like Age) can only have arrows pointing out
of them. A convenient side effect of this rule is that to get the total causal effect of a variable
like Age
there should never be any back door paths to close.
Do these rules contradict general relativity? If there was a variable for how fast or far someone had traveled in their life (e.g. if they were an astronaut) couldn’t that cause a change in their Age? Perhaps in this extra special case you’d need to define a standardized time that all other time measurements were made relative to.
This perspective leads to other questions. Age does not necessarily indicate the level of education we will eventually achieve. Are people who had a healthy upbringing, e.g. with more material resources, more or less likely to provide disapproving responses, and also more or less likely to get an education? There are plenty of possible confounds that could influence both eventual education level and response, such as temperament. Would a 14 year old who will eventually get a PhD answer the same way to these questions now and after they finish their PhD?
Practially speaking, based on the wording in the question and the last paragraph of section 12.4, it
seems the author wanted readers to only consider that Age
influences Education
(not vice versa).
See also the next question, which adds a few more comments on this issue.
Conditional independencies, if any, are listed under each DAG.
In this first model, education is a mediator between age and response:
library(dagitty)
age_dag <- dagitty('
dag {
bb="0,0,1,1"
Age [pos="0.4,0.2"]
Education [pos="0.6,0.2"]
Response [pos="0.5,0.3"]
Age -> Education
Education -> Response
}')
iplot(function() plot(age_dag), scale=10)
display(impliedConditionalIndependencies(age_dag))

Age _||_ Rspn | Edct
In this model, all variables are associated regardless of conditioning:
mediator_dag <- dagitty('
dag {
bb="0,0,1,1"
Age [pos="0.4,0.2"]
Education [pos="0.6,0.2"]
Response [pos="0.5,0.3"]
Age -> Education
Age -> Response
Education -> Response
}')
iplot(function() plot(mediator_dag), scale=10)
display(impliedConditionalIndependencies(mediator_dag))

In this model, education has no direct causal effect on the response:
independent_dag <- dagitty('
dag {
bb="0,0,1,1"
Age [pos="0.4,0.2"]
Education [pos="0.6,0.2"]
Response [pos="0.5,0.3"]
Age -> Education
Age -> Response
}')
iplot(function() plot(independent_dag), scale=10)
display(impliedConditionalIndependencies(independent_dag))

Edct _||_ Rspn | Age
This answer will not consider this model because, as stated above, it seems more reasonable to
assume that age leads to education than vice versa. Notice it has the same conditional
independencies as the last model, and we could construct similar models where Education -> Age
and
with the same conditional independencies as the first and second models. Although that approach
would work for this question, see comments in the next question.
independent_dag <- dagitty('
dag {
bb="0,0,1,1"
Age [pos="0.4,0.2"]
Education [pos="0.6,0.2"]
Response [pos="0.5,0.3"]
Education -> Age
Age -> Response
}')
iplot(function() plot(independent_dag), scale=10)
display(impliedConditionalIndependencies(independent_dag))

Edct _||_ Rspn | Age
To test which model is correct, we only need one more model beyond what we learned from m12.6
in
the text: a model that includes age as a predictor. If we find that the education coefficient bE
is near zero in this model, the third causal model is correct. If we find that the age coefficient
is near zero, the first causal model is correct. If both coefficients are non-zero, we can infer the
second causal model is correct. Of course, all causal inferences are limited by the variables we are
considering.
data(Trolley)
d <- Trolley
edu_levels <- c(6, 1, 8, 4, 7, 2, 5, 3)
d$edu_new <- edu_levels[d$edu]
dat <- list(
R = d$response,
action = d$action,
intention = d$intention,
contact = d$contact,
E = as.integer(d$edu_new),
StdAge = standardize(d$age),
alpha = rep(2, 7)
)
m_trolley_age <- ulam(
alist(
R ~ ordered_logistic(phi, kappa),
phi <- bE * sum(delta_j[1:E]) + bStdAge * StdAge + bA * action + bI * intention + bC * contact,
kappa ~ normal(0, 1.5),
c(bA, bI, bC, bE, bStdAge) ~ normal(0, 1),
vector[8]:delta_j <<- append_row(0, delta),
simplex[7]:delta ~ dirichlet(alpha)
),
data = dat, chains = 4, cores = min(detectCores(), 4), iter=2000
)
iplot(function() {
plot(precis(m_trolley_age, depth=2), main="m_trolley_age")
}, ar=1.8)
display_markdown("Raw data (preceding plot):")
display(precis(m_trolley_age, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff Rhat4
kappa[1] -2.68361325 0.08865314 -2.81809233 -2.55714508 934.8242 1.0045896
kappa[2] -1.99898884 0.08602333 -2.13072845 -1.87498384 956.9940 1.0043763
kappa[3] -1.41489799 0.08459355 -1.54145036 -1.29335548 930.5575 1.0041952
kappa[4] -0.39022666 0.08408004 -0.51812506 -0.27072486 912.5518 1.0048259
kappa[5] 0.27999016 0.08492709 0.14992044 0.40188572 911.4397 1.0043455
kappa[6] 1.18651290 0.08672497 1.05659667 1.31209534 934.4431 1.0038228
bStdAge -0.09906705 0.02135479 -0.13219621 -0.06503196 2383.0371 1.0002423
bE 0.23368574 0.10691089 0.07107577 0.37026784 725.3296 1.0056800
bC -0.95951533 0.04991494 -1.03958491 -0.88054049 4455.0481 1.0006841
bI -0.71919289 0.03580042 -0.77702762 -0.66303366 3907.5633 1.0004625
bA -0.70773068 0.03971195 -0.77110260 -0.64431515 3836.2686 1.0004352
delta[1] 0.10779876 0.07067768 0.02181789 0.23894013 4253.9584 1.0008411
delta[2] 0.12025954 0.07794762 0.02225073 0.26652180 4951.3060 0.9996026
delta[3] 0.08666906 0.06125356 0.01645026 0.19803202 2473.1828 1.0024484
delta[4] 0.06414256 0.04969173 0.01182441 0.15042389 1764.9911 1.0015228
delta[5] 0.44279365 0.13735187 0.18244217 0.63211812 915.1823 1.0039116
delta[6] 0.08272248 0.05728302 0.01653252 0.19350915 3970.9883 1.0004653
delta[7] 0.09561396 0.06337485 0.01884149 0.21152727 4067.4154 1.0001127
As implied in the text, it is not correct to causally infer the associated of education bE
is
negative based on m12.6
. In this model, the association of education with response has reversed,
becoming positive. Since the error bars on neither bStdAge
or bE
cross zero, we’ll infer for now
that the second causal model above is correct.
12H8. Consider one more variable in the trolley data: Gender. Suppose that gender might influence education as well as response directly. Draw the DAG now that includes response, education, age, and gender. Using only the DAG, is it possible that the inferences from 12H7 above are confounded by gender? If so, define any additional models you need to infer the causal influence of education on response. What do you conclude?
Answer. First, lets consider a DAG that doesn’t work. In the following DAG we assume
Education -> Age
rather than that Age -> Gender
:
incorrect_dag <- dagitty('
dag {
bb="0,0,1,1"
Education [pos="0.4,0.2"]
Age [pos="0.6,0.2"]
Response [pos="0.5,0.3"]
Gender [pos="0.3,0.3"]
Education -> Age
Education -> Response
Age -> Response
Gender -> Education
Gender -> Response
}')
iplot(function() plot(incorrect_dag), scale=10)
display(impliedConditionalIndependencies(incorrect_dag))

Age _||_ Gndr | Edct
With the additions implied by this question, this DAG implies that Gender
will be able to predict
Age
without conditioning, which is not correct.
Instead, we’ll add this question’s covariate to the second model from the last question:
gender_dag <- dagitty('
dag {
bb="0,0,1,1"
Age [pos="0.4,0.2"]
Education [pos="0.6,0.2"]
Response [pos="0.5,0.3"]
Gender [pos="0.7,0.3"]
Age -> Education
Education -> Response
Age -> Response
Gender -> Education
Gender -> Response
}')
iplot(function() plot(gender_dag), scale=10)
display(impliedConditionalIndependencies(gender_dag))

Age _||_ Gndr
This DAG implies Age
will always be independent of Gender
, which seems reasonable in this
context, ignoring that females live slightly longer than males. To infer either the direct or total
causal effect of Education
on Response
the required adjustment is to condition on both Age
and
Gender
:
display(adjustmentSets(gender_dag, exposure="Education", outcome="Response", effect="total"))
{ Age, Gender }
We can do so by fitting a model that includes all three variables as predictors of Response
:
data(Trolley)
d <- Trolley
edu_levels <- c(6, 1, 8, 4, 7, 2, 5, 3)
d$edu_new <- edu_levels[d$edu]
dat <- list(
R = d$response,
action = d$action,
intention = d$intention,
contact = d$contact,
E = as.integer(d$edu_new),
StdAge = standardize(d$age),
gid = ifelse(d$male == 1, 2, 1),
alpha = rep(2, 7)
)
# This fit takes more than three hours with 4 chains/cores.
m_trolley_gender <- ulam(
alist(
R ~ ordered_logistic(phi, kappa),
phi <- bE * sum(delta_j[1:E]) + bStdAge * StdAge + bGender[gid]
+ bA * action + bI * intention + bC * contact,
kappa ~ normal(0, 1.5),
bGender[gid] ~ normal(0, 1),
c(bA, bI, bC, bE, bStdAge) ~ normal(0, 1),
vector[8]:delta_j <<- append_row(0, delta),
simplex[7]:delta ~ dirichlet(alpha)
),
data = dat, chains = 2, cores = min(detectCores(), 2), iter=1600
)
flush.console()
iplot(function() {
plot(precis(m_trolley_gender, depth=2), main="m_trolley_gender")
}, ar=1.8)
display_markdown("Raw data (preceding plot):")
display(precis(m_trolley_gender, depth = 2), mimetypes="text/plain")

Raw data (preceding plot):
mean sd 5.5% 94.5% n_eff Rhat4
kappa[1] -2.266657325 0.46288890 -3.00485360 -1.5165712 573.0954 1.0016901
kappa[2] -1.582421818 0.46177762 -2.31693874 -0.8349806 572.3607 1.0018609
kappa[3] -0.993131342 0.46200481 -1.72707767 -0.2474544 570.7043 1.0017850
kappa[4] 0.051192897 0.46220897 -0.67312307 0.7982066 574.1828 1.0017592
kappa[5] 0.738588857 0.46212487 0.01278630 1.4784193 575.5787 1.0016047
kappa[6] 1.664283873 0.46311785 0.93955882 2.4109957 577.9330 1.0016990
bGender[1] 0.301659579 0.45912846 -0.41939969 1.0395143 561.7080 1.0012553
bGender[2] 0.862541251 0.46088869 0.13650775 1.6024925 565.7668 1.0011667
bStdAge -0.066221873 0.02193411 -0.10169805 -0.0318095 690.4008 0.9990740
bE 0.002776721 0.17079243 -0.27589601 0.2518510 334.8590 1.0001277
bC -0.973445079 0.04982326 -1.05475852 -0.8931221 1228.7393 1.0019420
bI -0.725770562 0.03634345 -0.78454521 -0.6683968 1536.4701 0.9997787
bA -0.713728611 0.03985093 -0.77665578 -0.6505416 1458.8686 1.0037167
delta[1] 0.150196008 0.09681367 0.03337167 0.3302251 1387.4576 0.9993277
delta[2] 0.146686397 0.09427857 0.03088713 0.3291819 1614.6417 0.9995624
delta[3] 0.138256497 0.09033356 0.02936921 0.2970326 1220.4606 1.0010909
delta[4] 0.136899029 0.09815857 0.02294144 0.3208046 904.8349 0.9990669
delta[5] 0.186898355 0.15213432 0.01796063 0.4656690 494.4601 1.0018564
delta[6] 0.114728654 0.07435170 0.02224116 0.2496879 1319.6433 1.0002505
delta[7] 0.126335060 0.08040863 0.02693776 0.2792283 1802.2701 0.9989585
At this point, education appears to have no causal effect on response. That is, the mean and HPDI of
bE
are near zero.