Chapter 8 Conditional manatees

8.1 Building an interaction

$R$ = ruggedness, $G$ = GDP, $C$ = continent, $U$ = unobserved variables.

library(dagitty)
library(rethinking)
dag_8.1 <- dagitty("dag{
                   R -> G
                   C -> G
                   U -> G
                   U -> R
                   }")
coordinates(dag_8.1) <- list(y = c(R = 0, G = 0, C = 0, U = 1),
                             x = c(R = 0, G = 1, U = 1, C = 2))
drawdag(dag_8.1)

$G = f(R,C)$

8.1.1 Making a rugged model

library(rethinking)
data(rugged)
d <- rugged

#log transform GDP
d$log_gdp <- log(d$rgdppc_2000)

#only include countries with GDP data
dd <- d[complete.cases(d$rgdppc_2000),]

#rescale variables
dd$log_gdp_std <- (dd$log_gdp) / mean(dd$log_gdp) # values of 1 is average
dd$rugged_std <- (dd$rugged) / max(dd$rugged) # values range from 0 to max ruggedness (1)

Basic model $$\text{log}(y_{i}) \sim \text{Normal}(\mu_{i}, \sigma)\\ \mu_{i} = \alpha + \beta(r_{i} - \overline{r})\\ \alpha \sim \text{Normal}|(1,1)\\ \beta \sim \text{Normal}(0, 1)\\ \sigma \sim \text{Exponential}(1)$$

In R:

m8.1 <- quap(
  alist(
    log_gdp_std ~ dnorm(mu, sigma),
    mu <- a + b*(rugged_std - 0.215),
    a ~ dnorm(1, 1),
    b ~ dnorm(0, 1),
    sigma ~ dexp(1)
  ), data = dd
)

Sample priors:

set.seed(11)
prior <- extract.prior(m8.1)

#set plot
plot(NULL, xlim=c(0,1), ylim=c(0.5, 1.5), xlab = "ruggedness", ylab = "log GDP")
abline(h=min(dd$log_gdp_std), lty = 2)
abline(h=max(dd$log_gdp_std), lty = 2)

#draw lines from prior
rugged_seq <- seq(from = -0.1, to = 1.1, length.out=30)
mu <- link(m8.1, post = prior, data = data.frame(rugged_std=rugged_seq))
for(i in 1:50){
  lines(rugged_seq, mu[i,], col=col.alpha('black',0.3))
  }

$\alpha$ is too wild. intercept should be somewhere around where the mean of ruggedness hits 1 on the log GDP scale so adjust to Normal(1, 0.1).

$\beta$ is also out of control. we need something (positive or negative) that spans the difference between the dashed lines

Slope should be $\pm 0.6$ which is the differece between the maximum and minimum values of GDP

max(dd$log_gdp_std) - min(dd$log_gdp_std)

## [1] 0.5658058

#proportion of slopes greater than 0.6
sum(abs(prior$b) > 0.6) / length(prior$b)

## [1] 0.54

Let’s fix the model

m8.1 <- quap(
  alist(
    log_gdp_std ~dnorm(mu,sigma),
    mu <-a+b*(rugged_std-0.215),
    a ~dnorm(1,0.1),
    b ~dnorm(0,0.3),
    sigma ~dexp(1)
  ), data = dd
)

precis(m8.1)

##              mean          sd        5.5%      94.5%
## a     0.999998578 0.010412457  0.98335746 1.01663970
## b     0.001994904 0.054795958 -0.08557962 0.08956943
## sigma 0.136503830 0.007397023  0.12468196 0.14832570

No association seen yet

8.1.2 Adding an indicator isn’t enough

Update $\mu$

$$\mu_{i} = \alpha_{CID[i]} + \beta(r_{i} - \overline{r})$$

#make an index variable for Africa (1) and other continents (2)
dd$cid <- ifelse(dd$cont_africa == 1, 1, 2)

Now update the model

m8.2 <- quap(
  alist(
    log_gdp_std ~ dnorm(mu, sigma),
    mu <-a[cid] + b * (rugged_std - 0.215),
    a[cid] ~ dnorm(1, 0.1),
    b ~ dnorm(0, 0.3),
    sigma ~dexp(1)
  ), data = dd
)

compare(m8.1, m8.2)

##           WAIC       SE    dWAIC      dSE    pWAIC       weight
## m8.2 -252.1508 15.25508  0.00000       NA 4.291725 1.000000e+00
## m8.1 -188.9664 13.28913 63.18439 15.11018 2.574284 1.904075e-14

precis(m8.2, depth = 2)

##              mean          sd       5.5%      94.5%
## a[1]   0.88040848 0.015938408  0.8549358 0.90588113
## a[2]   1.04916209 0.010186473  1.0328821 1.06544204
## b     -0.04651791 0.045690786 -0.1195406 0.02650479
## sigma  0.11239762 0.006092464  0.1026607 0.12213455

post <- extract.samples(m8.2)
diff_a1_a2 <- post$a[,1] - post$a[,2]
PI(diff_a1_a2)

##         5%        94% 
## -0.1988472 -0.1384092

rugged.seq <- seq(from = -0.1, to = 1.1, length.out = 30)
mu.NotAfrica <- link(m8.2, data = data.frame(cid=2, rugged_std=rugged.seq))
mu.Africa <- link(m8.2, data = data.frame(cid = 1, rugged_std = rugged.seq))
mu.NotAfrica_mu <- apply(mu.NotAfrica, 2, mean)
mu.NotAfrica_ci <- apply(mu.NotAfrica, 2, PI, prob = 0.97)
mu.Africa_mu <- apply(mu.Africa, 2, mean)
mu.Africa_ci <- apply(mu.Africa, 2, PI)

plot(NULL, xlim=c(0,1), ylim=c(0.5, 1.5), xlab = "ruggedness", ylab = "log GDP")
points(dd$rugged_std, dd$log_gdp_std, col = dd$cid, pch = 16)
lines(rugged.seq, mu.Africa_mu, lwd = 2, col = 1)
shade(mu.Africa_ci, rugged.seq)
lines(rugged.seq, mu.NotAfrica_mu, lwd = 2, col = 2)
shade(mu.NotAfrica_ci, rugged.seq, col = col.alpha(2, 0.3))

8.1.3 Adding an interaction does work

$$\mu_{i} = \alpha_{CID[i]} + \beta_{CID[i]}(r_{i} - \overline{r})$$

m8.3 <- quap(
  alist(
     log_gdp_std ~ dnorm(mu, sigma),
    mu <-a[cid] + b[cid] * (rugged_std - 0.215),
    a[cid] ~ dnorm(1, 0.1),
    b[cid] ~ dnorm(0, 0.3),
    sigma ~dexp(1)
  ), data = dd
)

precis(m8.3, depth = 2)

##             mean          sd        5.5%       94.5%
## a[1]   0.8865632 0.015675727  0.86151037  0.91161605
## a[2]   1.0505709 0.009936627  1.03469025  1.06645155
## b[1]   0.1325019 0.074204597  0.01390861  0.25109517
## b[2]  -0.1425818 0.054749512 -0.23008206 -0.05508147
## sigma  0.1094944 0.005935331  0.10000855  0.11898016

compare(m8.1, m8.2, m8.3, func=PSIS)

## Some Pareto k values are high (>0.5). Set pointwise=TRUE to inspect individual points.

##           PSIS       SE     dPSIS       dSE    pPSIS       weight
## m8.3 -258.8749 15.35435  0.000000        NA 5.314384 9.688047e-01
## m8.2 -252.0033 15.30167  6.871589  6.936621 4.354277 3.119533e-02
## m8.1 -188.4136 13.40307 70.461243 15.674795 2.856543 4.850334e-16

plot(PSIS(m8.3, pointwise = TRUE)$k)

## Some Pareto k values are high (>0.5). Set pointwise=TRUE to inspect individual points.

8.1.4 Plotting the interaction

par(mfrow=c(1,2))
# plot Africa - cid = 1
d.A1 <-dd[dd$cid == 1,]
plot(d.A1$rugged_std, d.A1$log_gdp_std, pch=16, col=rangi2,
xlab="ruggedness (standardized)",ylab="log GDP (as proportion of mean)",
xlim=c(0,1) )
mu <-link(m8.3,data=data.frame(cid=1,rugged_std=rugged_seq))
mu_mean <-apply(mu,2,mean)
mu_ci <-apply(mu,2,PI,prob=0.97)
lines( rugged_seq,mu_mean,lwd=2)
shade( mu_ci,rugged_seq,col=col.alpha(rangi2,0.3))
mtext("African nations")
# plotnon-Africa-cid=2
d.A0 <-dd[dd$cid==2,]
plot( d.A0$rugged_std,d.A0$log_gdp_std,pch=1,col="black",
xlab="ruggedness (standardized)",ylab="logGDP(asproportionofmean)",
xlim=c(0,1) )
mu <-link(m8.3,data=data.frame(cid=2,rugged_std=rugged_seq))
mu_mean <-apply(mu,2,mean)
mu_ci <-apply(mu,2,PI,prob=0.97)
lines( rugged_seq,mu_mean,lwd=2)
shade( mu_ci,rugged_seq)
mtext("Non-African nations")

8.2 Symmetry of interactions

You can break an interaction into 2 identical phrasings
1. GDP ~ ruggedness depending on Africa
2. Africa ~ GDP depending on rugedness

$$\mu_{i} = (2 - CID_{i})(\alpha_{1} + \beta_{1}(r_{i} - \overline{r})) + (CID_{i} - 1)(\alpha_{2} + \beta_{2}(r_{i} - \overline{r}))$$

rugged_seq <- seq(from = -0.2, to = 1.2, length.out = 30)
muA <- link(m8.3, data=data.frame(cid=1, rugged_std=rugged_seq))
muN <- link(m8.3, data=data.frame(cid=2, rugged_std=rugged_seq))
delta <- muA - muN


mu.delta <- apply(delta, 2, mean)
PI.delta <- apply(delta, 2, PI)

plot(x=rugged_seq, type = 'n', xlim = c(0,1), ylim = c(-0.3, 0.2),
     xlab = 'ruggedness (std)', ylab = 'expected difference log GDP')
shade(PI.delta, rugged_seq, col='grey')
abline(h = 0, lty = 2)
text(x = 0.2, y = 0, label = "Africa higher GDP\nAfrica lower GDP")
lines(rugged_seq, mu.delta)

At high ruggedness, being in Africa gives higher than expected GDP.

8.3 Continuous interactions

8.3.1 A winter flower

data(tulips)
d <- tulips
str(d)

## 'data.frame':    27 obs. of  4 variables:
##  $ bed   : Factor w/ 3 levels "a","b","c": 1 1 1 1 1 1 1 1 1 2 ...
##  $ water : int  1 1 1 2 2 2 3 3 3 1 ...
##  $ shade : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ blooms: num  0 0 111 183.5 59.2 ...

8.3.2 the models

Water and Shade work together to create Blooms; $W \rightarrow B \leftarrow S ; B = f(W,S)$

1. water
   $$\beta_{i} \sim \text{Normal}(\mu_{i}, \sigma)\\ \mu_{i} = \alpha + \beta_{W}(W_{i} - \overline{W})\\ \alpha \sim \text{Normal}(0.5,1)\\ \beta_{W} \sim \text{Normal}(0,1)\\ \sigma \sim \text{Exponential}(1)$$
2. shade
   $$\beta_{i} \sim \text{Normal}(\mu_{i}, \sigma)\\ \mu_{i} = \alpha + \beta_{S}(S_{i} - \overline{S})\\ \alpha \sim \text{Normal}(0.5,1)\\ \beta_{S} \sim \text{Normal}(0,1)\\ \sigma \sim \text{Exponential}(1)$$
3. water + shade
   $$\beta_{i} \sim \text{Normal}(\mu_{i}, \sigma)\\ \mu_{i} = \alpha + \beta_{W}(W_{i} - \overline{W}) + \beta_{S}(S_{i} - \overline{S})\\ \alpha \sim \text{Normal}(0.5,1)\\ \beta_{W} \sim \text{Normal}(0,1)\\ \beta_{S} \sim \text{Normal}(0,1)\\ \sigma \sim \text{Exponential}(1)$$
4. water * shade
   

#center predictors and scale outcome
d$blooms_std <- d$blooms / max(d$blooms)
d$water_cent <- d$water - mean(d$water)
d$shade_cent <- d$shade - mean(d$shade)

The $\alpha$ prior is likely too broad. We need it to be between 0 and 1. How much is outside that?

a <- rnorm(1e4, 0.5, 1); sum(a < 0 | a > 1) / length(a)

## [1] 0.623

Let’s tighten it

a <- rnorm(1e4, 0.5, 0.25); sum(a < 0 | a > 1) / length(a)

## [1] 0.047

range of water and shade are each 2 units. range of blooms is one unit. max slopes = 2/1 (0.5) so we can set the prior to 0 with 0.25 sd to get values ranging from -0.5 to 0.5.

#water
m8.4a <- quap(
  alist(
    blooms_std ~ dnorm(mu, sigma),
    mu <- a + bw*water_cent,
    a ~ dnorm(0.5,0.25),
    bw ~ dnorm(0, 0.25),
    sigma ~ dexp(1)
  ), data = d
)

#shade
m8.4b <- quap(
  alist(
    blooms_std ~ dnorm(mu, sigma),
    mu <- a + bs*shade_cent,
    a ~ dnorm(0.5,0.25),
    bs ~ dnorm(0, 0.25),
    sigma ~ dexp(1)
  ), data = d
)

# water + shade
m8.4c <- quap(
  alist(
    blooms_std ~ dnorm(mu, sigma),
    mu <- a + bw*water_cent + bs*shade_cent,
    a ~ dnorm(0.5,0.25),
    bw ~ dnorm(0, 0.25),
    bs ~ dnorm(0, 0.25),
    sigma ~ dexp(1)
  ), data = d
)

#water * shade
m8.4d <- quap(
  alist(
    blooms_std ~ dnorm(mu, sigma),
    mu <- a + bw*water_cent + bs*shade_cent + bws*water_cent*shade_cent,
    a ~ dnorm(0.5,0.25),
    bw ~ dnorm(0, 0.25),
    bs ~ dnorm(0, 0.25),
    bws ~ dnorm(0, 0.25),
    sigma ~ dexp(1)
  ), data = d
)

prior simulations

set.seed(11)
prior_a <- extract.prior(m8.4a) #water
prior_b <- extract.prior(m8.4b) #shade
prior_c <- extract.prior(m8.4c) #water + shade
prior_d <- extract.prior(m8.4d) #water * shade


#set plot
plot(NULL, xlim=c(0,2), ylim=c(0, 1.25), xlab = "Water / shade", ylab = "blooms")
abline(h=min(d$blooms_std), lty = 2)
abline(h=max(d$blooms_std), lty = 2)

#draw lines from prior
water_seq <- seq(from = -1.1, to = 1.1, length.out=30)
shade_seq <- seq(from = -1.1, to = 1.1, length.out=30)

mu_a <- link(m8.4a, post = prior_a, data = data.frame(water_cent=water_seq))
mu_b <- link(m8.4b, post = prior_b, data = data.frame(shade_cent=shade_seq))
mu_c <- link(m8.4c, post = prior_c, data = data.frame(water_cent=water_seq, shade_cent=shade_seq))
mu_d <- link(m8.4d, post = prior_d, data = data.frame(water_cent=water_seq, shade_cent=shade_seq))

#set plot
plot(NULL, xlim=c(-1,1), ylim=c(0, 1.25), xlab = "Water / shade (centered)", ylab = "blooms (std)")
abline(h=min(d$blooms_std), lty = 2)
abline(h=max(d$blooms_std), lty = 2)


for(i in 1:50){
  lines(water_seq, mu_a[i,], col=col.alpha('blue',0.3))
  }

for(i in 1:50){
  lines(water_seq, mu_b[i,], col=col.alpha('black',0.3))
}
for(i in 1:50){
  lines(water_seq, mu_c[i,], col=col.alpha('green',0.3))
}

for(i in 1:50){
  lines(water_seq, mu_d[i,], col=col.alpha('red',0.3))
  }

text(x = -0.75, y = 1.2, label = "Water", col = "blue")
text(x = -0.45, y = 1.2, label = "Shade", col = "black")
text(x = 0, y = 1.2, label = "Water + Shade", col = "green")
text(x = 0.5, y = 1.2, label = "Water * Shade", col = 'red')

precis(m8.4a)

##            mean         sd      5.5%     94.5%
## a     0.3594883 0.03502089 0.3035181 0.4154584
## bw    0.2034854 0.04270338 0.1352371 0.2717336
## sigma 0.1837433 0.02489828 0.1439510 0.2235355

precis(m8.4b)

##             mean         sd       5.5%       94.5%
## a      0.3611171 0.04403635  0.2907385  0.43149566
## bs    -0.1097643 0.05351857 -0.1952973 -0.02423127
## sigma  0.2323701 0.03143809  0.1821260  0.28261427

precis(m8.4c)

##             mean         sd       5.5%       94.5%
## a      0.3587452 0.03022116  0.3104459  0.40704441
## bw     0.2050352 0.03689240  0.1460741  0.26399641
## bs    -0.1125324 0.03687853 -0.1714714 -0.05359337
## sigma  0.1581668 0.02144796  0.1238888  0.19244481

precis(m8.4d)

##             mean         sd        5.5%       94.5%
## a      0.3579980 0.02391747  0.31977332  0.39622278
## bw     0.2067288 0.02923282  0.16000910  0.25344849
## bs    -0.1134595 0.02922579 -0.16016795 -0.06675103
## bws   -0.1431791 0.03567746 -0.20019860 -0.08615967
## sigma  0.1248375 0.01693790  0.09776742  0.15190750

8.3.3 Plotting posterior predictions

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$shade_cent == s,]
  plot(x = idx$water_cent, y = idx$blooms_std, xlim = c(-1,1), ylim = c(0,1),
       xlab = "water", ylab = "blooms", pch = 16, col = rangi2)
  mu <- link(m8.4c, data = data.frame(shade_cent=s, water_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
   mtext(concat("m8.4c post: shade = ", s))
}

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$shade_cent == s,]
  plot(x = idx$water_cent, y = idx$blooms_std, xlim = c(-1,1), ylim = c(0,1),
       xlab = "water", ylab = "blooms", pch = 16, col = rangi2)
  mu <- link(m8.4d, data = data.frame(shade_cent=s, water_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
   mtext(concat("m8.4d post: shade = ", s))
}

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$water_cent == s,]
  plot(x = idx$shade_cent, y = idx$blooms_std, xlim = c(-1,1), ylim = c(0,1),
       xlab = "shade", ylab = "blooms", pch = 16, col = rangi2)
  mu <- link(m8.4c, data = data.frame(water_cent=s, shade_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
  mtext(concat("m8.4c post: water = ", s))
}

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$water_cent == s,]
  plot(x = idx$shade_cent, y = idx$blooms_std, xlim = c(-1,1), ylim = c(0,1),
       xlab = "shade", ylab = "blooms", pch = 16, col = rangi2)
  mu <- link(m8.4d, data = data.frame(water_cent=s, shade_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
  mtext(concat("m8.4c post: water = ", s))
}

8.3.4 Plotting prior predictions

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$shade_cent == s,]
  plot(x = idx$water_cent, y = idx$blooms_std, type = 'n', xlim = c(-1,1), ylim = c(-0.5,1.5),
       xlab = "water", ylab = "blooms")
  mu <- link(m8.4c, post = prior_c, data = data.frame(shade_cent=s, water_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
  lines(-1:1, mu[11,], lwd = 2, col = rangi2)
  abline(h = 0, lty = 2)
  abline(h = 1, lty = 2)
  mtext(concat("m8.4c post: shade = ", s))
}

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$shade_cent == s,]
  plot(x = idx$water_cent, y = idx$blooms_std, type = 'n', xlim = c(-1,1), ylim = c(-0.5,1.5),
       xlab = "water", ylab = "blooms")
  mu <- link(m8.4d, post = prior_d, data = data.frame(shade_cent=s, water_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
  lines(-1:1, mu[11,], lwd = 2, col = rangi2)
  abline(h = 0, lty = 2)
  abline(h = 1, lty = 2)
  mtext(concat("m8.4c post: shade = ", s))
}

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$water_cent == s,]
  plot(x = idx$shade_cent, y = idx$blooms_std, type = 'n', xlim = c(-1,1), ylim = c(-0.5,1.5),
       xlab = "shade", ylab = "blooms")
  mu <- link(m8.4c, post = prior_c, data = data.frame(water_cent=s, shade_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
  lines(-1:1, mu[11,], lwd = 2, col = rangi2)
  abline(h = 0, lty = 2)
  abline(h = 1, lty = 2)
  mtext(concat("m8.4c post: water = ", s))
}

par(mfrow = c(1,3))
for (s in -1:1){
  idx <- d[d$water_cent == s,]
  plot(x = idx$shade_cent, y = idx$blooms_std, type = 'n', xlim = c(-1,1), ylim = c(-0.5,1.5),
       xlab = "shade", ylab = "blooms")
  mu <- link(m8.4d, post = prior_d, data = data.frame(water_cent=s, shade_cent = -1:1))
  for(i in 1:20){
    lines(-1:1, mu[i,], col = col.alpha('black',0.3))
  }
  lines(-1:1, mu[11,], lwd = 2, col = rangi2)
  abline(h = 0, lty = 2)
  abline(h = 1, lty = 2)
  mtext(concat("m8.4c post: water = ", s))
}