Assignment #9
Chapter 12
Shiyu Bai
2023-01-30
Chapter 12 - Monsters & 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.
df <- c(12, 36, 7, 41)
q <- df/sum(df)
q## [1] 0.12500000 0.37500000 0.07291667 0.42708333p <- cumsum(q)
p## [1] 0.1250000 0.5000000 0.5729167 1.0000000log_odds <- log(p/(1-p))
log_odds## [1] -1.9459101 0.0000000 0.2937611 Inf12-2. Make a version of Figure 12.5 for the employee ratings data given just above.
c_prop=df/sum(df)
c_prop=sapply(1:4,function(i){sum(c_prop[1:i])})
plot(y=c_prop,x=1:4,type="b", ylim=c(0,1))
segments(1:4,0,1:4,c_prop)
for(i in 1:4){segments(i+0.05,c(0,c_prop)[i],i+0.05,c_prop[i], col = "red")}
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 4414.857 1000.499 0.00000 NA 125.65537 9.999998e-01
## m2 4445.876 1074.092 31.01883 139.0969 77.55112 1.838007e-07y <- 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 )
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?
data(Hurricanes)
d <- Hurricanes
d$fem_std <- (d$femininity - mean(d$femininity)) / sd(d$femininity) # standardised femininity
dat <- list(D = d$deaths, F = d$fem_std)The gamma-Poisson has two parameters, one for the rate and the other for the dispersion of rates. When dispersion is high, the distribution is closer to a pure Poisson process. For consistency, we keep the same priors as before and add a simple exponential prior for the scale parameter, which has a wider range than previous predictions. The second model also has a higher effective number of samples, making it less prone to correlations, which can be quantified using the WAIC. High dispersion values result in a wider distribution, which is closer to a “true” Poisson distribution.
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?
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)
dat$P <- standardize(d$min_pressure)
dat$S <- standardize(d$damage_norm)