Assignment #9
Chapter 12
Valerie Gono
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.
# Create data set
ratings <- c(12, 36, 7, 41)
# Compute proportions of each element in ratings
proportions <- ratings / sum(ratings)
# Compute cumulative sum of proportions
cum_proportions <- cumsum(proportions)
# Compute odds
odds <- cum_proportions / (1-cum_proportions)
# Compute log cumulative odds of each rating
log_odds <- log(odds)
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.
# Plot
plot (1:4, cum_proportions, type="b",
xlab="Ratings", ylab="Cumulative Proportion", ylim=c(0,1), xlim=c(0.5,4.5),
main="Cumulative Proportions of Ratings")
# Plot gray bars over each outcome to show cumulative probability
for (x in 1:4) lines(c(x, x), c(0, cum_proportions[x]), col = "gray", lwd = 4)
# Plot blue line segments to show the discrete probability of each individual outcome
for (x in 1:4) lines(c(x, x) + 0.05, c(cum_proportions[x] - proportions[x], cum_proportions[x]), col = "blue", lwd = 4)
# Labels
text(1:4 + 0.2, cum_proportions - proportions / 2, labels=1:4, col="blue", cex=1.5)
12-3. In 2014, a paper was published that was entitled “Female hurricanes are deadlier than male hurricanes”. 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?
# Read-in data
data(Hurricanes)
d <- Hurricanes
# Add transformed femininity variable
d$femininity_std <- (d$femininity - mean(d$femininity)) / sd(d$femininity)
# Build simple model
simple_mod <- ulam(
alist(
deaths ~ dpois(lambda),
log(lambda) <- a + b*fem,
a ~ dnorm(0,10),
b ~ dnorm(0,1)
),
data=list(
deaths=d$deaths,
fem=d$femininity_std
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(simple_mod)## mean sd 5.5% 94.5% n_eff Rhat4
## a 3.0004358 0.02398149 2.961011 3.0395216 1356.140 1.0035001
## b 0.2388104 0.02517820 0.198338 0.2777833 1427.633 0.9983317# The simple model shows a positive association between femininity of hurricane names and deaths.
# Build intercept only model
int_mod <- ulam(
alist(
deaths ~ dpois(lambda),
log(lambda) <- a,
a ~ dnorm(0,10)
),
data=list(deaths=d$deaths
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(int_mod)## mean sd 5.5% 94.5% n_eff Rhat4
## a 3.02812 0.0227608 2.991532 3.064604 728.2449 1.005113# Compare model
compare(int_mod, simple_mod)## WAIC SE dWAIC dSE pWAIC weight
## simple_mod 4416.274 1001.495 0.00000 NA 133.90953 1.000000e+00
## int_mod 4454.969 1078.628 38.69494 142.0289 78.84276 3.958235e-09# Comparison between the simple model and the intercept only model showed that WAIC is lower for the model which include femininity of names.
# Check to see if model retrodict well
# PLot raw data
plot(d$femininity_std, d$deaths, pch=16, lwd=2,
col=rangi2, xlab="Name Femininity" , ylab="Deaths",
main="Check Model Prediction")
# Compute trend from model
prediction_d <- list(fem=seq(from=-2,to=2,length.out=1e2))
prediction <- link(simple_mod, data=prediction_d)
prediction_mean <- apply(prediction, 2, mean)
prediction_interval <- apply(prediction, 2, PI)
# Superimpose trend to plot
lines(prediction_d$fem, prediction_mean)
# Compute sampling distribution
death_post_sim <- sim(simple_mod,data=prediction_d)
death_post_sim.PI <- apply(death_post_sim, 2, PI)
# Superimpose sampling interval - dashed lines
lines(prediction_d$fem, death_post_sim.PI[1,] , lty=2)
lines(prediction_d$fem, death_post_sim.PI[2,] , lty=2)
# From our plot, we can see that the model does not retrodict well. While the model is confident with its prediction, we see that most death tolls are either above or below the band restricted by the dashed lines. This is a sign of over-dispersion, which the model does not take into account. We do see a weak positive trend which is highly influenced by the larger death tolls that occurred on the right side of the chart. 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?
# Read-in data
data(Hurricanes)
d <- Hurricanes
# Add transformed femininity variable
d$femininity_std <- (d$femininity - mean(d$femininity)) / sd(d$femininity)
# Build a Gamma-Poisson model
gam_poisson_mod <- ulam(
alist(
deaths ~ dgampois(lambda, scale),
log(lambda) <- a + b*fem,
a ~ dnorm(0,10),
b ~ dnorm(0,1),
scale ~ dexp(1)
),
data=list(
deaths=d$deaths,
fem=d$femininity_std
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(gam_poisson_mod)## mean sd 5.5% 94.5% n_eff Rhat4
## a 3.0254579 0.16110221 2.77822611 3.2848185 1917.445 0.9999162
## b 0.2094084 0.15459363 -0.04709664 0.4476123 1928.454 0.9998885
## scale 0.4536779 0.06450874 0.36125002 0.5654199 1975.334 0.9984546# The Gamma-Poisson model still shows a positive association between femininity of hurricane names and deaths.
# Check to see if model retrodict well
# PLot raw data
plot(d$femininity_std, d$deaths, pch=16, lwd=2,
col=rangi2, xlab="Name Femininity" , ylab="Deaths",
main="Check Model Prediction")
# Compute trend from model
prediction_d <- list(fem=seq(from=-2,to=2,length.out=1e2))
prediction <- link(gam_poisson_mod, data=prediction_d)
prediction_mean <- apply(prediction, 2, mean)
prediction_interval <- apply(prediction, 2, PI)
# Superimpose trend to plot
lines(prediction_d$fem, prediction_mean)
# Compute sampling distribution
death_post_sim <- sim(gam_poisson_mod, data=prediction_d)
death_post_sim.PI <- apply(death_post_sim, 2, PI)
# Superimpose sampling interval - dashed lines
lines(prediction_d$fem, death_post_sim.PI[1,] , lty=2)
lines(prediction_d$fem, death_post_sim.PI[2,] , lty=2)
# Wider band of dashed lines.
# Show wavering of positive association
plot(coeftab(simple_mod, gam_poisson_mod))
# Plot shows that on the posterior for the femininity coefficient is more uncertain for the Gamma-Poisson model in comparison to the simple model. This association is diminished due the nature of the gamma distribution, which allows for individual computation of death rates. The gamma distribution are informed by a and b in our model, and since it's done repeatedly for many individual computation of death rates, it widens our posterior 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?
# Read-in data
data(Hurricanes)
d <- Hurricanes
# Add transformed femininity variable and damage_norm
d$femininity_std <- (d$femininity - mean(d$femininity)) / sd(d$femininity)
# Standardized min_pressure and damage_norm
d$min_pressure_std <- (d$min_pressure - mean(d$min_pressure)) / sd(d$min_pressure)
d$damage_norm_std <- (d$damage_norm - mean(d$damage_norm)) / sd(d$damage_norm)
# Model with femininity and minimum pressure with interaction
fem_min_press_mod_int <- ulam(
alist(
deaths ~ dgampois(lambda, scale),
log(lambda) <- a + bf*fem + bm*minp + bfm*fem*minp,
a ~ dnorm(1, 1),
c(bf,bm, bfm) ~ dnorm(0,1),
scale ~ dexp(1)
),
data=list(
deaths=d$deaths,
fem=d$femininity_std,
minp=d$min_pressure_std
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(fem_min_press_mod_int)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.7528697 0.14171039 2.52835344 2.9782445 2482.754 1.0000077
## bfm 0.3028388 0.14835607 0.07566799 0.5459886 2783.417 1.0002948
## bm -0.6711998 0.14258343 -0.88868707 -0.4374135 2247.800 0.9992889
## bf 0.3010074 0.13920527 0.07581048 0.5157033 2742.835 1.0003921
## scale 0.5543735 0.07660023 0.44069274 0.6766872 2698.133 0.9990649# Model with femininity and minimum pressure without interaction
fem_min_press_mod <- ulam(
alist(
deaths ~ dgampois(lambda, scale),
log(lambda) <- a + bf*fem + bm*minp,
a ~ dnorm(1, 1),
c(bf,bm) ~ dnorm(0,1),
scale ~ dexp(1)
),
data=list(
deaths=d$deaths,
fem=d$femininity_std,
minp=d$min_pressure_std
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(fem_min_press_mod)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.7757224 0.14053844 2.55271346 2.9981571 1823.656 0.9994707
## bm -0.5895972 0.13817902 -0.80447376 -0.3739658 1922.557 0.9990967
## bf 0.2866462 0.14127591 0.05427609 0.5057955 2104.527 0.9993789
## scale 0.5382735 0.08209751 0.41361912 0.6789402 2316.603 0.9985025# Compare model with and without interactions
compare(fem_min_press_mod, fem_min_press_mod_int)## WAIC SE dWAIC dSE pWAIC weight
## fem_min_press_mod_int 694.3395 38.50040 0.000000 NA 8.542271 0.8821908
## fem_min_press_mod 698.3662 41.44551 4.026684 5.543618 8.753209 0.1178092# The model with interaction has a lower WAIC than the model without interaction. For the model with interaction, as minimum pressure reduce, the storm gets stronger, which causes more death. The femininity coefficient is still positive, but does not overlap 0. This is similar to the bfm coefficient which is the interaction between femininity and minimum pressure.
# Plotting to check on result
press <- seq(from = -3, to = 2, length.out = 1e2)
# 'masculine' storms
prediction_d_m <- data.frame(fem = -1, minp = press)
prediction_m <- link(fem_min_press_mod_int, data = prediction_d_m)
prediction_mean_m <- apply(prediction_m, 2, mean)
prediction_interval_m <- apply(prediction_m, 2, PI)
# 'feminine' storms
prediction_d_f <- data.frame(fem = 1, minp = press)
prediction_f <- link(fem_min_press_mod_int, data = prediction_d_f)
prediction_mean_f <- apply(prediction_f, 2, mean)
prediction_interval_f <- apply(prediction_f, 2, PI)
# Plotting, use sqrt of deaths for scaling purposes
plot(d$min_pressure_std, sqrt(d$deaths),
pch = 1, lwd = 2, col = ifelse(d$femininity_std > 0, "red", "dark gray"),
xlab = "Normalized Minimum Pressure", ylab = "Sqrt of Deaths",
main="Normalized Minimum Pressure and Deaths"
)
lines(press, sqrt(prediction_mean_m), lty = 2, col = "dark gray")
lines(press, sqrt(prediction_mean_f), lty = 2, col = "red")
# On average, the masculine (gray) storms are less deadly than the feminine (red) storms. However, as Minimum Pressure reduces, we see that the masculine storms leads to more deaths. However, at the higher minimum pressure (to the right side), the difference between masculine and feminine deaths are closer together.
# Model with femininity and damage norm with interaction
fem_min_dam_mod_int <- ulam(
alist(
deaths ~ dgampois(lambda, scale),
log(lambda) <- a + bf*fem + bd*damn + bfd*fem*damn,
a ~ dnorm(1, 1),
c(bf, bd, bfd) ~ dnorm(0,1),
scale ~ dexp(1)
),
data=list(
deaths=d$deaths,
fem=d$femininity_std,
damn=d$damage_norm_std
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(fem_min_dam_mod_int)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.57203776 0.12972262 2.3646113 2.7771137 3070.269 0.9992286
## bfd 0.30519098 0.20025233 -0.0173470 0.6234954 2165.766 1.0002979
## bd 1.25269919 0.21776652 0.9066918 1.6102407 2464.601 1.0007694
## bf 0.08863789 0.13029913 -0.1274077 0.2977380 3108.156 0.9994661
## scale 0.68200183 0.09816009 0.5325255 0.8449778 2208.256 0.9982251# Model with femininity and minimum pressure without interaction
fem_min_dam_mod <- ulam(
alist(
deaths ~ dgampois(lambda, scale),
log(lambda) <- a + bf*fem + bd*damn,
a ~ dnorm(1, 1),
c(bf, bd) ~ dnorm(0,1),
scale ~ dexp(1)
),
data=list(
deaths=d$deaths,
fem=d$femininity_std,
damn=d$damage_norm_std
), chains=4, cores=4, log_lik=TRUE)## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
## /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1precis(fem_min_dam_mod)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.58401389 0.13178894 2.3727492 2.7983683 2504.276 0.9991704
## bd 1.22726226 0.23004971 0.8613325 1.5841426 2183.588 1.0013330
## bf 0.09962618 0.13175443 -0.1083490 0.3052264 1943.945 1.0004342
## scale 0.67479940 0.09780692 0.5294284 0.8358938 2089.326 1.0000149# Compare model with and without interactions
compare(fem_min_dam_mod, fem_min_dam_mod_int)## WAIC SE dWAIC dSE pWAIC weight
## fem_min_dam_mod_int 670.4158 33.85500 0.0000000 NA 6.702325 0.5669006
## fem_min_dam_mod 670.9542 33.62997 0.5384337 4.323696 5.643369 0.4330994# The model without the interaction has a lower WAIC than the model with the interaction. However, the difference is minimum. Checking for counterfactual predictions for model with interactions again.
# Plotting to check on result
dam <- seq(from = -1, to = 6, length.out = 1e2)
# 'masculine' storms
prediction_d_m <- data.frame(fem = -1, damn = dam)
prediction_m <- link(fem_min_dam_mod_int, data = prediction_d_m)
prediction_mean_m <- apply(prediction_m, 2, mean)
prediction_interval_m <- apply(prediction_m, 2, PI)
# 'feminine' storms
prediction_d_f <- data.frame(fem = 1, damn = dam)
prediction_f <- link(fem_min_dam_mod_int, data = prediction_d_f)
prediction_mean_f <- apply(prediction_f, 2, mean)
prediction_interval_f <- apply(prediction_f, 2, PI)
# Plotting, use sqrt of deaths for scaling purposes
plot(d$damage_norm_std, sqrt(d$deaths),
pch = 1, lwd = 2, col = ifelse(d$femininity_std > 0, "red", "dark gray"),
xlab = "Normalized Damage", ylab = "Sqrt of Deaths",
main="Normalized Damage and Deaths"
)
lines(dam, sqrt(prediction_mean_m), lty = 2, col = "dark gray")
lines(dam, sqrt(prediction_mean_f), lty = 2, col = "red")
# With this new model, it can be observed that there are less difference between masculine and feminine hurricanes overall. The prediction is still not great however. However, we do see a strong interaction effect, likely due to a few "outlier" feminine hurricane that skew our model. We would need more data to create a better model for prediction.