Chapter 12 - Monsters and Mixtures

This chapter introduced several new types of regression, all of which are generalizations of generalized linear models (GLMs). Ordered logistic models are useful for categorical outcomes with a strict ordering. They are built by attaching a cumulative link function to a categorical outcome distribution. Zero-inflated models mix together two different outcome distributions, allowing us to model outcomes with an excess of zeros. Models for overdispersion, such as beta-binomial and gamma-Poisson, draw the expected value of each observation from a distribution that changes shape as a function of a linear model.

Place each answer inside the code chunk (grey box). The code chunks should contain a text response or a code that completes/answers the question or activity requested. Make sure to include plots if the question requests them.

Finally, upon completion, name your final output .html file as: YourName_ANLY505-Year-Semester.html and publish the assignment to your R Pubs account and submit the link to Canvas. Each question is worth 5 points.

Questions

12-1. 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.

n <- c(12, 36, 7, 41)
q <- n / sum(n) 
q
## [1] 0.12500000 0.37500000 0.07291667 0.42708333
sum(q)
## [1] 1
p <- cumsum(q) 
p
## [1] 0.1250000 0.5000000 0.5729167 1.0000000
log(p/(1-p))
## [1] -1.9459101  0.0000000  0.2937611        Inf

12-2. Make a version of Figure 12.5 for the employee ratings data given just above.

n = c(12, 36, 7, 41)
p = n / sum(n)
cum_p = cumsum(p)
plot(
  y = cum_p,
  x = 1:4,
  type = "b",
  ylim = c(0, 1)
)
segments(1:4, 0, 1:4, cum_p)
for (i in 1:4) {
    segments(i + 0.05, c(0, cum_p)[i], i + 0.05, cum_p[i], col = "green")
}

12-3. In 2014, a paper was published that was entitled “Female hurricanes are deadlier than male hurricanes.”191 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.

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 quap 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?

data(Hurricanes)
data1 <- Hurricanes
m1 <- map(
  alist(
    deaths ~ dpois( lambda ),
    log(lambda) <- a + bF*femininity,
    a ~ dnorm(0,10),
    bF ~ dnorm(0,5)
  ) ,
  data=data1)

m2 <- map(
  alist(
    deaths ~ dpois( lambda ),
    log(lambda) <- a ,
    a ~ dnorm(0,10)
  ) ,
  data=data1)

compare(m1,m2)
##        WAIC       SE    dWAIC      dSE     pWAIC       weight
## m1 4427.365 1008.028  0.00000       NA 134.04635 9.999999e-01
## m2 4459.241 1081.337 31.87607 139.5389  87.84111 1.197289e-07
y <- sim(m1)

y.mean <- colMeans(y)
y.PI <- apply(y, 2, PI)

plot(y=data1$deaths, x=data1$femininity, col=rangi2, ylab="deaths", xlab="femininity", pch=16)
points(y=y.mean, x=data1$femininity, pch=1)
segments(x0=data1$femininity, x1= data1$femininity, y0=y.PI[1,], y1=y.PI[2,])

lines(y= y.mean[order(data1$femininity)], x=sort(data1$femininity))
lines( y.PI[1,order(data1$femininity)], x=sort(data1$femininity), lty=2 )
lines( y.PI[2,order(data1$femininity)], x=sort(data1$femininity), lty=2 )

# We can tell that femininity is not strong related to number of deaths, especially at the high end. 
# Poisson models exhibits high over-dispersion, thus poorly fitted.

12-4. 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?

library(rethinking)
library(cmdstanr) 
 
data(Hurricanes)

d <- Hurricanes # load data on object called d
d$fem_std <- (d$femininity - mean(d$femininity)) / sd(d$femininity) # standardised femininity
dat <- list(D = d$deaths, F = d$fem_std)

# model formula
f <- alist(
  D ~ dpois(lambda), # poisson outcome distribution
  log(lambda) <- a + bF * F, # log-link for lambda with linear model
  # priors in log-space, 0 corresponds to outcome of 1
  a ~ dnorm(1, 1),
  bF ~ dnorm(0, 1)
)

N <- 1e3
a <- rnorm(N, 1, 1)
bF <- rnorm(N, 0, 1)
F_seq <- seq(from = -2, to = 2, length.out = 30) # sequence from -2 to 2 because femininity data is standardised
plot(NULL,
  xlim = c(-2, 2), ylim = c(0, 500),
  xlab = "name femininity (std)", ylab = "deaths"
)
for (i in 1:N) {
  lines(F_seq,
    exp(a[i] + bF[i] * F_seq), # inverse link to get outcome scale
    col = grau(), lwd = 1.5
  )
}

12-5. 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?

library(rethinking)
data(Hurricanes)
data=Hurricanes
data$fem_std = (data$femininity - mean(data$femininity)) / sd(data$femininity)
dat = list(D = data$deaths, F = data$fem_std)
dat$P = standardize(data$min_pressure)
dat$S = standardize(data$damage_norm)

##Unable to instal CmdStan

mH3a = ulam( alist( D ~ dgampois(lambda, scale), log(lambda) <- a + bF * F + bP * P + bFP * F * P, a ~ dnorm(1, 1), c(bF, bP, bFP) ~ dnorm(0, 1), scale ~ dexp(1) ), data = dat, cores = 4, chains = 4, log_lik = TRUE ) precis(mH3a)

P_seq = seq(from = -3, to = 2, length.out = 1e2) d_pred = data.frame(F = -1, P = P_seq) lambda_m = link(mH3a, data = d_pred) lambda_m.mu = apply(lambda_m, 2, mean) lambda_m.PI = apply(lambda_m, 2, PI)

d_pred = data.frame(F = 1, P = P_seq) lambda_f = link(mH3a, data = d_pred) lambda_f.mu = apply(lambda_f, 2, mean) lambda_f.PI = apply(lambda_f, 2, PI)

Plotting, sqrt() to make differences easier to spot

plot(dat\(P, sqrt(dat\)D), pch = 1, lwd = 2, col = ifelse(dat$F > 0, “red”, “dark gray”), xlab = “minimum pressure (std)”, ylab = “sqrt(deaths)” ) lines(P_seq, sqrt(lambda_m.mu), lty = 2) shade(sqrt(lambda_m.PI), P_seq) lines(P_seq, sqrt(lambda_f.mu), lty = 1, col = “blue”) shade(sqrt(lambda_f.PI), P_seq, col = col.alpha(“blue”, 0.2))

ggplot(as.data.frame(PSISk(mH3a)), aes(x = PSISk(mH3a)))+ theme_bw() +labs(title = “Paraeto-K values”, subtitle = “Values > 1 indicate highly influential data”)