Chapter 11 - God Spiked the Integers

This chapter described some of the most common generalized linear models, those used to model counts. It is important to never convert counts to proportions before analysis, because doing so destroys information about sample size. A fundamental difficulty with these models is that parameters are on a different scale, typically log-odds (for binomial) or log-rate (for Poisson), than the outcome variable they describe. Therefore computing implied predictions is even more important than before.

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

11E1. If an event has probability 0.35, what are the log-odds of this event?

log_odds <- log(0.35/(1-0.35))
log_odds
## [1] -0.6190392
#log odd of this event is -0.6190392

11E2. If an event has log-odds 3.2, what is the probability of this event?

prob <- 1/(1+exp(-3.2))
prob
## [1] 0.9608343
#probability of this event is 0.9608343

11E3. Suppose that a coefficient in a logistic regression has value 1.7. What does this imply about the proportional change in odds of the outcome?

coeff <- exp(1.7)
coeff
## [1] 5.473947
#The coefficient in a logistic regression with a value 1.7 implies that log odds will change the outcome by ~5.47 times

11E4. Why do Poisson regressions sometimes require the use of an offset? Provide an example.

#Poisson regressions are measured as events per unit time so in order to make sure they're on scale, an offset is required sometimes. An example would be weather predictions.

11M1. As explained in the chapter, binomial data can be organized in aggregated and disaggregated forms, without any impact on inference. But the likelihood of the data does change when the data are converted between the two formats. Can you explain why?

#The reason why the likelihood of the data does change when the data are converted between the two formats is because during this conversion the multiplier variable c(n,m) is also changed into a constant.

11M2. If a coefficient in a Poisson regression has value 1.7, what does this imply about the change in the outcome?

#Poisson regression value of 1.7 implies that change in the outcome will increase by 5.47 times

11M3. Explain why the logit link is appropriate for a binomial generalized linear model.

#A binomial generalized linear model has a binary outcome, and probability can only lie between 0 and 1 (inclusive), logit link satisfies this criteria so it's appropriate for a binomial generalized linear model.

11M4. Explain why the log link is appropriate for a Poisson generalized linear model.

#A Poisson generalized linear model has a positive outcome, and logit link will be between 0 and 1 which satisfies this criteria so it's appropriate for a Poisson generalized linear model.

11M5. What would it imply to use a logit link for the mean of a Poisson generalized linear model? Can you think of a real research problem for which this would make sense?

#It implies that the outcome of a Poisson generalized model will be positive and between 0 and 1 (inclusive).

11M6. State the constraints for which the binomial and Poisson distributions have maximum entropy. Are the constraints different at all for binomial and Poisson? Why or why not?

#The constraints for binomial and Poisson distributions to have maximum entropy are: for binomial, binary events and fixed probability; for Poisson, the variance must be equal to expected values. The constraints are different because the distributions are different.

11M7. Use quap to construct a quadratic approximate posterior distribution for the chimpanzee model that includes a unique intercept for each actor, m11.4 (page 330). Compare the quadratic approximation to the posterior distribution produced instead from MCMC. Can you explain both the differences and the similarities between the approximate and the MCMC distributions? Relax the prior on the actor intercepts to Normal(0,10). Re-estimate the posterior using both ulam and quap. Do the differences increase or decrease? Why?

data("chimpanzees")
df <-  chimpanzees
df$recipient <- NULL

#Re-estimating using quap
m <- quap(
          alist(
            pulled_left ~ dbinom(1,p),
            logit(p) <- a[actor] + (xp + xpc * condition)*prosoc_left,
            a[actor] ~ dnorm(0,10),
            xp ~ dnorm(0,10),
            xpc ~ dnorm(0,10)
          ),
          data = df
)
pairs(m)

#There seems to be not much of a difference.

11M8. Revisit the data(Kline) islands example. This time drop Hawaii from the sample and refit the models. What changes do you observe?

data("Kline")

df2 <- Kline

df2$P <- scale(log(df2$population))

df2$contact_id <- ifelse(df2$contact == "high", 2,1)

df2
##       culture population contact total_tools mean_TU            P contact_id
## 1    Malekula       1100     low          13     3.2 -1.291473310          1
## 2     Tikopia       1500     low          22     4.7 -1.088550750          1
## 3  Santa Cruz       3600     low          24     4.0 -0.515764892          1
## 4         Yap       4791    high          43     5.0 -0.328773359          2
## 5    Lau Fiji       7400    high          33     5.0 -0.044338980          2
## 6   Trobriand       8000    high          19     4.0  0.006668287          2
## 7       Chuuk       9200    high          40     3.8  0.098109204          2
## 8       Manus      13000     low          28     6.6  0.324317564          1
## 9       Tonga      17500    high          55     5.4  0.518797917          2
## 10     Hawaii     275000     low          71     6.6  2.321008320          1

11H1. Use WAIC or PSIS to compare the chimpanzee model that includes a unique intercept for each actor, m11.4 (page 330), to the simpler models fit in the same section. Interpret the results.

m2 <- map(
          alist(
            pulled_left ~ dbinom(1,p),
            logit(p) <- x,
            x ~ dnorm(0,10)
          ),
          data = df
)

m3 <- map(
          alist(
            pulled_left ~ dbinom(1,p),
            logit(p) <- x + xp*prosoc_left,
            x ~ dnorm(0, 10),
            xp ~ dnorm(0,10)
          ),
          data = df
)

m4 <- map(
          alist(
            pulled_left ~ dbinom(1,p),
            logit(p) <- x + (xp + xpc*condition)*prosoc_left,
            x ~ dnorm(0, 10),
            xp ~ dnorm(0,10),
            xpc ~ dnorm(0,10)
          ),
          data = df
)

m5 <- map(
          alist(
            pulled_left ~ dbinom(1,p),
            logit(p) <- x[actor] + (xp + xpc*condition)*prosoc_left,
            x[actor] ~ dnorm(0, 10),
            xp ~ dnorm(0,10),
            xpc ~ dnorm(0,10)
          ),
          data = df
)
#using compare function to compare the models
compare(m2, m3, m4, m5)
##        WAIC        SE    dWAIC      dSE      pWAIC       weight
## m5 544.4368 18.902931   0.0000       NA 12.9131064 1.000000e+00
## m3 680.5346  9.289938 136.0978 18.34593  2.0188128 2.797277e-30
## m4 682.3067  9.337720 137.8699 18.27002  2.9818336 1.153274e-30
## m2 687.8470  7.092063 143.4103 19.12292  0.9533607 7.225223e-32
#The results indicate that m5 with unique intercept for each actor has the least WAIC value with a weight of 1 whereas other models have a much closer WAIC values and weight of 0. Based on pWAIC values, m5 seems to be more flexible to fit the data. The standard error value is also higher for m5 when compared with other models.