Assignment #10

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. Problems are labeled Easy (E), Medium (M), and Hard (H).

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

12E1. What is the difference between an ordered categorical variable and an unordered one? Define and then give an example of each.

# An ordered categorical variable has an explicit ranking of levels, such as a 1-to-5 scale on a questionnaire. However, the difference between each level is not equal to other differences. For example, the distance from 1 to 2 does not necessarily equal the difference between 4 and 5.

# An unordered categorical variable is not constrained to any order among the values. The different values merely represent different discrete outcomes, without any implied ordering. Examples include choices among colors, individual identities, and individual words.

12E2. What kind of link function does an ordered logistic regression employ? How does it differ from an ordinary logit link?

# Ordered logistic regression uses a cumulative logit link function. The difference with an ordinary logit link is that the probability is a cumulative probability instead of a discrete probability of a single event.

12E3. When count data are zero-inflated, using a model that ignores zero-inflation will tend to induce which kind of inferential error?

# Ignoring zero-inflation will tend to underestimate the rate of events, because a count distribution with extra zeros added to it will have a lower mean.

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 underdispersed counts?

# Over-dispersion occurs when there is more variation at the min and max than expected, leading to a distribution with thicker tails. A typical example is aggregating count data over a period of time.

# Under-dispersion occurs when there is less variation in the data than predicted. Usually financial data (stock prices) are autocorrelated, and would be causing underdispersed counts.

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.

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

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

plot(1:4, p, xlab = "rating", ylab = "cumulative proportion", xlim = c(0.7, 4.3), ylim = c(0, 1), xaxt = "n")
axis(1, at = 1:4, labels = 1:4)

# plot gray cumulative probability lines
for (x in 1:4) lines(c(x, x), c(0, p[x]), col = "gray", lwd = 2)

# plot blue discrete probability segments
for (x in 1:4) lines(c(x, x) + 0.1, c(p[x] - q[x], p[x]), col = "slateblue", lwd = 2)

# add number labels
text(1:4 + 0.2, p - q/2, labels = 1:4, col = "slateblue")

12M3. Can you modify the derivation of the zero-inflated Poisson distribution (ZIPoisson) from the chapter to construct a zero-inflated binomial distribution?

# The probability of a zero: 
# Pr(0|p_0, q, n) = p_0 + (1 − p_0)(1 − q)^n

# The probability of a non-zero observation (y): 
# Pr(y|p_0, q, n) = (1 − p_0)(n!/(y!(n − y)!))(q^y)((1 − q)^(n−y))

12H1. 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. Load the data with:

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

?Hurricanes

## starting httpd help server ... done

d <- Hurricanes
d$fmnnty_std <- standardize(d$femininity)
dat <- list(D = d$deaths, FM = d$fmnnty_std)

m12.h1.1 <- ulam(alist(
    D ~ dpois(lambda),
    log(lambda) <- a + bF*FM,
    a ~ dnorm(0, 10),
    bF ~ dnorm(0, 1)),
    data = dat, chains = 4, log_lik = TRUE)

## 
## SAMPLING FOR MODEL '10c653a22e9ae6c6229a552021610d3b' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 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.113 seconds (Warm-up)
## Chain 1:                0.124 seconds (Sampling)
## Chain 1:                0.237 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '10c653a22e9ae6c6229a552021610d3b' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 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.128 seconds (Warm-up)
## Chain 2:                0.115 seconds (Sampling)
## Chain 2:                0.243 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '10c653a22e9ae6c6229a552021610d3b' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 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.053 seconds (Warm-up)
## Chain 3:                0.059 seconds (Sampling)
## Chain 3:                0.112 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '10c653a22e9ae6c6229a552021610d3b' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 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.123 seconds (Warm-up)
## Chain 4:                0.125 seconds (Sampling)
## Chain 4:                0.248 seconds (Total)
## Chain 4:

precis(m12.h1.1)

##         mean         sd      5.5%     94.5%    n_eff    Rhat4
## a  3.0002991 0.02421138 2.9609280 3.0400093 1039.288 1.001843
## bF 0.2403679 0.02473479 0.2010309 0.2801747 1301.772 1.003025

dat2 <- list(D = d$deaths)

m12.h1.2 <- ulam(alist(
    D ~ dpois(lambda),
    log(lambda) <- a,
    a ~ dnorm(0, 10)),
    data = dat2, chains = 4, log_lik = TRUE)

## 
## SAMPLING FOR MODEL '44dfa17951690fd17b6e6a2c8851ff60' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 0 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 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.016 seconds (Warm-up)
## Chain 1:                0.016 seconds (Sampling)
## Chain 1:                0.032 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '44dfa17951690fd17b6e6a2c8851ff60' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 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.016 seconds (Warm-up)
## Chain 2:                0.02 seconds (Sampling)
## Chain 2:                0.036 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '44dfa17951690fd17b6e6a2c8851ff60' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 0 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 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.017 seconds (Warm-up)
## Chain 3:                0.016 seconds (Sampling)
## Chain 3:                0.033 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '44dfa17951690fd17b6e6a2c8851ff60' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 0 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 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.02 seconds (Warm-up)
## Chain 4:                0.013 seconds (Sampling)
## Chain 4:                0.033 seconds (Total)
## Chain 4:

precis(m12.h1.2)

##       mean         sd     5.5%    94.5%    n_eff    Rhat4
## a 3.028129 0.02285499 2.991708 3.064325 536.9151 1.003097

compare(m12.h1.2, m12.h1.1, func = PSIS)

## 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.

##              PSIS        SE    dPSIS      dSE     pPSIS      weight
## m12.h1.1 4381.406  985.2235  0.00000       NA 114.66134 1.00000e+00
## m12.h1.2 4426.971 1067.1565 45.56468 145.5334  68.91161 1.27572e-10

# plot raw data
plot(d$fmnnty_std, d$deaths, pch = 16, col = rangi2, xlab = "femininity", ylab = "deaths")

# compute model-based trend
pred_dat <- list(FM = seq(from = -2, to = 1.5, length.out = 30))
lambda <- link(m12.h1.1, data = pred_dat)
lambda.mu <- apply(lambda, 2, mean)
lambda.PI <- apply(lambda, 2, PI)

# superimpose trend
lines(pred_dat$FM, lambda.mu)
shade(lambda.PI, pred_dat$FM)

# compute sampling distribution
deaths_sim <- sim(m12.h1.1, data = pred_dat)
deaths_sim.PI <- apply(deaths_sim, 2, PI)

# superimpose sampling interval as dashed lines
lines(pred_dat$FM, deaths_sim.PI[1, ], lty = 2)
lines(pred_dat$FM, deaths_sim.PI[2, ], lty = 2)

# The model including the femininity rating fits the data better than the intercept-only model, but the plot shows that the reason for this better fit is due to a few influential outliers on the right side of the plot.

# There are four storms that were particularly deadly (over 100 deaths), and all four had names that received high ratings of femininity, leading the model to believe that female hurricanes are deadlier than male hurricanes.