#added “knitr::opts_chunk$set(error = TRUE)” because I could not knit the document. All the code works fine in R, just would not knit.

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)
d <- df/sum(df)
d
## [1] 0.12500000 0.37500000 0.07291667 0.42708333
d.2 <- cumsum(q)
## Error in eval(expr, envir, enclos): cannot coerce type 'closure' to vector of type 'double'
d.2
## Error in eval(expr, envir, enclos): object 'd.2' not found
log <- log(d.2/(1-d.2))
## Error in eval(expr, envir, enclos): object 'd.2' not found
log
## function (x, base = exp(1))  .Primitive("log")

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

prop=df/sum(df)
prop=sapply(1:4,function(i){sum(prop[1:i])})
plot(y=prop,x=1:4,type="b", ylim=c(0,2))

segments(1:4,0,1:4,c_prop)
## Error in segments(1:4, 0, 1:4, c_prop): object 'c_prop' not found
for(i in 1:4)
{segments(i+0.05,
          c(0,c_prop)[i],
          i+0.05,
          c_prop[i],
          col = "red")}
## Error in segments(i + 0.05, c(0, c_prop)[i], i + 0.05, c_prop[i], col = "red"): object 'c_prop' not found

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")

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

m.2 <- map(
  alist(
    deaths ~ dpois( lambda ),
    log(lambda) <- a ,
    a ~ dnorm(0,10)
  ) ,
  data=data)
compare(m.1,m.2)
##         WAIC        SE   dWAIC      dSE     pWAIC       weight
## m.1 4405.979  996.3649  0.0000       NA 127.69296 1.000000e+00
## m.2 4461.023 1080.0996 55.0448 145.4782  88.09173 1.114738e-12
y <- sim(m.1)

mean <- colMeans(y)

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

plot(y=data$deaths, x=data$femininity, col=rangi2, ylab="deaths", xlab="femininity", pch=20)
points(y=y.mean, x=data$femininity, pch=1)
## Error in xy.coords(x, y): object 'y.mean' not found
segments(x0=data$femininity, x1= data$femininity, y0=PI.y[1,], y1=PI.y[2,])

lines(y= mean[order(data$femininity)],  x=sort(data$femininity))

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

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

#The relationship is fairly weak. This tell us that female hurricanes are have a higher relation to deaths. 

#The second model fits well and it is utilizing femininity as it's only predictor.

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)
data <- Hurricanes 
data$fem_std <- (data$femininity - mean(data$femininity)) / sd(data$femininity)
data.output <- list(D = data$deaths, F = data$fem_std)
data.output
## $D
##  [1]   2   4   3   1   0  60  20  20   0 200   7  15   1   0  22  50   0  46   3
## [20]   3   5  37   3  75   6   3  15   3 256  22   0   2   0 117   1  21   5   0
## [39]   1  15   5   2  21   3   0   1   4   8  12   5   3   5   0   1  13  21   3
## [58]  15  62   3   6   9   8  26  10   3   3   1   0  56   8   2   3  51   1  10
## [77]   7   8  25   5   1  15   1  62   5   1   1  52  84  41   5 159
## 
## $F
##  [1] -0.0009347453 -1.6707580424 -0.9133139565  0.9458703621  0.4810742824
##  [6]  0.4122162926  0.5499353710  0.8253673305  0.5327193242  0.9630864089
## [11] -0.2591568553  0.0679232445  0.9630864089  0.9630864089  0.9286574139
## [16]  0.7737253874  0.6015773141  0.8425833773  0.9802993570  0.3605712508
## [21]  0.7909383355  0.6360063090  0.8253673305  0.4810742824  0.6187933608
## [26]  0.4638613343  0.1539972812  0.6704353039  0.7048673975  0.8253673305
## [31]  0.5327193242  0.1884262762  0.9975154038  0.5843643659  0.6015773141
## [36]  0.6704353039  1.1352313836  0.0334942496 -1.5846840057 -1.5674710576
## [41] -1.3264649944 -1.2748199526  0.9458703621  0.9802993570 -1.5846840057
## [46] -1.4125390310  0.9114413671  0.8425833773 -1.4986130677  0.8942284190
## [51]  0.8081543823 -1.2059619628 -1.5330420626  0.4810742824  0.7048673975
## [56] -1.2059619628 -1.3781100361 -1.5846840057 -1.4125390310  0.9458703621
## [61]  0.1367812344  0.5327193242  0.5327193242  0.1195682863 -1.4125390310
## [66]  0.8081543823 -1.5158260158 -1.3953229842 -1.3781100361 -1.5330420626
## [71]  0.7737253874  1.1008023886  0.7392963925  0.8081543823 -0.8100238730
## [76] -1.2059619628 -0.2419408085 -1.2748199526 -1.7740450272  0.5327193242
## [81]  0.9802993570 -1.3436810411  0.7392963925  0.8425833773  0.5671483191
## [86] -1.3608939893  0.9458703621 -1.5674710576 -1.5158260158  0.7737253874
## [91] -1.4986130677  0.6876513507
# In order to make any significant comparisons, we must keep the same priors. For the scale parameter we can just use an exponential prior. It can be determined that this has more spread than the last predictions. I also notice that there are more samples in the second model meaning this model is less exposed to correlations. This is where a WAIC can be used. The dispersion number usually spreads the distribution out, so that there are larger values of the dispersion with a real 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)
data <- Hurricanes 
data$fem_std <- (data$femininity - mean(data$femininity)) / sd(data$femininity)
df <- list(D = data$deaths, F = data$fem_std)
df$P <- standardize(data$min_pressure)
df$S <- standardize(data$damage_norm)
# I start with a basic model which builds on the previous gamma-Poisson models by adding an interaction between femininity and min_pressure: When the minimum pressure gets lower, the storm grows stronger. The model here makes less of a distinction between feminine and masucline hurricanes.When we get closer to the right portion of the plot, the distances grow rapidly. This is most likely very hard for the model to take into consideration. The interaction effect is strong because there are 3 or 4 highly influential feminine storms at the upper-righthand corner of our plot. This tells us that when feminine storms are very damaging, they create a high death toll. I dont really think the effect sizes are plasuble for many reasons, mostly because this is a very simple model.