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)

$log_cum_odds =

-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:

\[\begin{split} \begin{align} Pr(0 | p_z,n,p_b) & = p_z + (1-p_z)Pr(0|n,p_b) \\ & = p_z + (1-p_z)\binom{n}{0}p_b^0(1-p_b)^{n-0} \\ & = p_z + (1-p_z)(1-p_b)^n \end{align} \end{split}\]

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
)

A data.frame: 92 × 10
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))

A data.frame: 6 × 12
	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

[Package rethinking version 2.13 ]

The head of Fish:

display(head(Fish))

A data.frame: 6 × 6
	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.

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

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

Executable Notes

Practice: Chp. 12

Practice: Chp. 12#

Fishing data

Description

Usage

Format

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