Assignment #9
Chapter 12
Saurabh Nehete
2023-02-01
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 valyh mue 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)
p <- cumsum(q)
logodds <- log(p/(1-p))
logodds## [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(1:4,p,xlab="rating",ylab="cumulativeproportion",
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")
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?
#example(stan_model, package = “rstan”, run.dontrun = TRUE)
data(hurricanes)
library(rethinking)
data(Hurricanes)
d <- Hurricanes
d$fmnnty_std <- ( d$femininity- mean(d$femininity) )/sd(d$femininity)
m1 <- map2stan(
alist( deaths ~ dpois(lambda),
log(lambda) <- a + bf*fmnnty,
a ~ dnorm(0,10),
bf ~ dnorm(0,1)
),
data=list(
deaths=d$deaths,
fmnnty=d$fmnnty_std
) , chains=4 )##
## SAMPLING FOR MODEL 'ed5c538b45ceca816d5b9bc85caa1644' 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 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.079 seconds (Warm-up)
## Chain 1: 0.072 seconds (Sampling)
## Chain 1: 0.151 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'ed5c538b45ceca816d5b9bc85caa1644' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.062 seconds (Warm-up)
## Chain 2: 0.063 seconds (Sampling)
## Chain 2: 0.125 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'ed5c538b45ceca816d5b9bc85caa1644' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.068 seconds (Warm-up)
## Chain 3: 0.083 seconds (Sampling)
## Chain 3: 0.151 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'ed5c538b45ceca816d5b9bc85caa1644' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.089 seconds (Warm-up)
## Chain 4: 0.109 seconds (Sampling)
## Chain 4: 0.198 seconds (Total)
## Chain 4:precis(m1)## mean sd 5.5% 94.5% n_eff Rhat4
## a 3.0003155 0.02295428 2.9632239 3.0370862 2542.442 1.0003746
## bf 0.2385263 0.02535505 0.1980824 0.2778208 2345.516 0.9992429m0<-map2stan(
alist( deaths~dpois(lambda),
log(lambda)<-a,
a~dnorm(0,10) ),
data=list(deaths=d$deaths),chains=4)##
## SAMPLING FOR MODEL '17eb1b0b15637d7ebd2937eaec48a17d' 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 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.041 seconds (Warm-up)
## Chain 1: 0.038 seconds (Sampling)
## Chain 1: 0.079 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '17eb1b0b15637d7ebd2937eaec48a17d' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.06 seconds (Warm-up)
## Chain 2: 0.054 seconds (Sampling)
## Chain 2: 0.114 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '17eb1b0b15637d7ebd2937eaec48a17d' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.054 seconds (Warm-up)
## Chain 3: 0.05 seconds (Sampling)
## Chain 3: 0.104 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '17eb1b0b15637d7ebd2937eaec48a17d' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.048 seconds (Warm-up)
## Chain 4: 0.051 seconds (Sampling)
## Chain 4: 0.099 seconds (Total)
## Chain 4:compare(m0,m1)## WAIC SE dWAIC dSE pWAIC weight
## m1 4409.635 997.3152 0.00000 NA 126.90346 1.000000e+00
## m0 4443.394 1072.4423 33.75904 142.1065 84.29563 4.670014e-08#plot raw data
plot(d$fmnnty_std,d$deaths,pch=16,
col=rangi2,xlab="femininity",ylab="deaths")
#compute model-based trend
pred_dat<-list(fmnnty=seq(from=-2,to=1.5,length.out=30))
lambda<-link(m1,data=pred_dat) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda.mu<-apply(lambda,2,mean)
lambda.PI<-apply(lambda,2,PI)
#super impose trend
lines(pred_dat$fmnnty,lambda.mu)
shade(lambda.PI,pred_dat$fmnnty)
#compute sampling distribution
deaths_sim<-sim(m1,data=pred_dat) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]deaths_sim.PI<-apply(deaths_sim,2,PI)
#superimpose sampling interval as dashed lines
lines(pred_dat$fmnnty,deaths_sim.PI[1,],lty=2)
lines(pred_dat$fmnnty,deaths_sim.PI[2,],lty=2)
## The sampling distribution is narrow and so is the 89% interval of expected value. Femininity accounts for very little of the variation in deaths.There is a lot of over-dispersion, which is common in Poisson models. So, this homogeneous Poisson model does a poor job for most of the hurricanes in the sample, as most of them lie outside the prediction envelope (the dashed boundaries).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?
library(rethinking)
data(Hurricanes)
d<-Hurricanes
d$fmnnty_std<-(d$femininity-mean(d$femininity))/sd(d$femininity)
m3<-map2stan(
alist(
deaths~dgampois(lambda,scale),
log(lambda)<-a+bf*fmnnty,
a~dnorm(0,10),
bf~dnorm(0,1),
scale~dcauchy(0,1)
#dexp would also be fine here
),
data=list( deaths=d$deaths,
fmnnty=d$fmnnty_std
),chains=4)##
## SAMPLING FOR MODEL '402c4bbc3d55ec520ca69377f3b47bb3' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.541 seconds (Warm-up)
## Chain 1: 0.451 seconds (Sampling)
## Chain 1: 0.992 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '402c4bbc3d55ec520ca69377f3b47bb3' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.474 seconds (Warm-up)
## Chain 2: 0.378 seconds (Sampling)
## Chain 2: 0.852 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '402c4bbc3d55ec520ca69377f3b47bb3' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.488 seconds (Warm-up)
## Chain 3: 0.396 seconds (Sampling)
## Chain 3: 0.884 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '402c4bbc3d55ec520ca69377f3b47bb3' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.373 seconds (Warm-up)
## Chain 4: 0.38 seconds (Sampling)
## Chain 4: 0.753 seconds (Total)
## Chain 4:precis(m3)## mean sd 5.5% 94.5% n_eff Rhat4
## a 3.0301443 0.16401969 2.77936481 3.3012052 3183.963 1.0006108
## bf 0.2125296 0.15820758 -0.04174512 0.4582841 2972.191 0.9994542
## scale 0.4543599 0.06237849 0.36100037 0.5592017 3175.316 0.9994333library(arulesViz)
library(raster)
#plot raw data
plot(d$fmnnty_std,d$deaths,pch=16,
col=rangi2,xlab="femininity",ylab="deaths")
#compute model-based trend
pred_dat<-list(fmnnty=seq(from=-2,to=1.5,length.out=30))
lambda<-link(m3,data=pred_dat) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda.mu<-apply(lambda,2,mean)
lambda.PI<-apply(lambda,2,PI)
#super impose trend
lines(pred_dat$fmnnty,lambda.mu)
shade(lambda.PI,pred_dat$fmnnty)
#compute sampling distribution
deaths_sim<-sim(m3,data=pred_dat) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]deaths_sim.PI<-apply(deaths_sim,2,PI)
#super impose sampling interval as dashed lines
lines(pred_dat$fmnnty,deaths_sim.PI[1,],lty=2)
lines(pred_dat$fmnnty,deaths_sim.PI[2,],lty=2)
post <- extract.samples(m3)
fem <- (-1)
# 1 stddev below mean
for ( i in 1:100 )
curve( dgamma2(x,exp(post$a[i]+post$bf[i]*fem),post$scale[i]) ,
from=0 , to=70 , xlab="mean deaths" , ylab="Density" ,
ylim=c(0,0.19) , col=col.alpha("black",0.2) ,
add=ifelse(i==1,FALSE,TRUE) )
mtext( concat("femininity = ",fem) )
fem <- (0)
# 1 stddev below mean
for ( i in 1:100 )
curve( dgamma2(x,exp(post$a[i]+post$bf[i]*fem),post$scale[i]) ,
from=0 , to=70 , xlab="mean deaths" , ylab="Density" ,
ylim=c(0,0.19) , col=col.alpha("black",0.2) ,
add=ifelse(i==1,FALSE,TRUE) )
mtext( concat("femininity = ",fem) )
fem <- (1)
# 1 stddev below mean
for ( i in 1:100 )
curve( dgamma2(x,exp(post$a[i]+post$bf[i]*fem),post$scale[i]) ,
from=0 , to=70 , xlab="mean deaths" , ylab="Density" ,
ylim=c(0,0.19) , col=col.alpha("black",0.2) ,
add=ifelse(i==1,FALSE,TRUE) )
mtext( concat("femininity = ",fem) )
## Note that the 89% interval forbf now overlaps zero, because it has more than 5 times the standard deviation as the analogous parameter in model m1. Model m3 has nearly the same posterior means for the intercept and slope bf, but much more uncertainty.
## Each gray curve is a gamma distribution of mean death rates, sampled from the posterior distribution of the model. A gamma distribution is defined by two parameters: a mean and a scale. The mean in this model is controlled by the linear model and its two parameters, a and bf. The code above takes single values of a and bf from the posterior distribution and builds a single linear model. It’s exponential in the code because this model uses a log link. 100 gamma distributions are drawn for each given value of femininity. This visualizes the uncertainty in the posterior about the variation in death rates.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?
# standardize new predictor
d$min_pressure_std <- (d$min_pressure- mean(d$min_pressure))/sd(d$min_pressure)
# interaction model
m_f_p_fp <- map2stan(
alist(
deaths ~ dpois(lambda),
log(lambda) <- a + bf*fmnnty + bp*minpress +
bfp*fmnnty*minpress,
a ~ dnorm(0,10),
c(bf,bp,bfp) ~ dnorm(0,1)
),
data=list(
deaths=d$deaths,
fmnnty=d$fmnnty_std,
minpress=d$min_pressure_std
) , chains=4 )##
## SAMPLING FOR MODEL '3743497a812dbf035f9ad51590a00035' 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 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.315 seconds (Warm-up)
## Chain 1: 0.289 seconds (Sampling)
## Chain 1: 0.604 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '3743497a812dbf035f9ad51590a00035' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.308 seconds (Warm-up)
## Chain 2: 0.288 seconds (Sampling)
## Chain 2: 0.596 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '3743497a812dbf035f9ad51590a00035' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.353 seconds (Warm-up)
## Chain 3: 0.249 seconds (Sampling)
## Chain 3: 0.602 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '3743497a812dbf035f9ad51590a00035' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.319 seconds (Warm-up)
## Chain 4: 0.302 seconds (Sampling)
## Chain 4: 0.621 seconds (Total)
## Chain 4:m_f_p<-map2stan(
alist(
deaths~dpois(lambda),
log(lambda)<-a+bf*fmnnty+bp*minpress,
a~dnorm(0,10),
c(bf,bp)~dnorm(0,1) ),
data=list( deaths=d$deaths,
fmnnty=d$fmnnty_std,
minpress=d$min_pressure_std
),chains=4)##
## SAMPLING FOR MODEL 'e9d54f68d666c9ecfac25c407be230d3' 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 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.154 seconds (Warm-up)
## Chain 1: 0.156 seconds (Sampling)
## Chain 1: 0.31 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'e9d54f68d666c9ecfac25c407be230d3' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.165 seconds (Warm-up)
## Chain 2: 0.171 seconds (Sampling)
## Chain 2: 0.336 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'e9d54f68d666c9ecfac25c407be230d3' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.179 seconds (Warm-up)
## Chain 3: 0.193 seconds (Sampling)
## Chain 3: 0.372 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'e9d54f68d666c9ecfac25c407be230d3' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.115 seconds (Warm-up)
## Chain 4: 0.171 seconds (Sampling)
## Chain 4: 0.286 seconds (Total)
## Chain 4:compare(m_f_p,m_f_p_fp)## WAIC SE dWAIC dSE pWAIC weight
## m_f_p 3485.094 1110.612 0.00000 NA 202.6558 1.000000e+00
## m_f_p_fp 3582.416 1149.235 97.32205 52.60328 264.4032 7.358433e-22precis(m_f_p_fp)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.72181926 0.02961156 2.6742930 2.7682659 2111.031 1.001230
## bf 0.26079178 0.03230637 0.2085626 0.3139309 1944.905 0.999845
## bp -0.73537083 0.02246576 -0.7706629 -0.6988443 1868.975 1.002384
## bfp 0.07028868 0.02591167 0.0293863 0.1121150 1696.698 1.000261minpress_seq<-seq(from=-3,to=2,length.out=30)
#'masculine'storms
d_pred<-data.frame(
fmnnty=-1,
minpress=minpress_seq)
lambda_m<-link(m_f_p_fp,data=d_pred) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda_m.mu<-apply(lambda_m,2,mean)
lambda_m.PI<-apply(lambda_m,2,PI)
#'feminine'storms
d_pred<-data.frame(
fmnnty=1,
minpress=minpress_seq)
lambda_f<-link(m_f_p_fp,data=d_pred) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda_f.mu<-apply(lambda_f,2,mean)
lambda_f.PI<-apply(lambda_f,2,PI)
# using sqrt scale for deaths so that the difference is easier to see
plot(d$min_pressure_std,sqrt(d$deaths),
pch=ifelse(d$fmnnty_std>0,16,1),col=rangi2,
xlab="minimumpressure",ylab="sqrt(deaths)")
lines(minpress_seq,sqrt(lambda_m.mu),lty=2)
shade(sqrt(lambda_m.PI),minpress_seq)
lines(minpress_seq,sqrt(lambda_f.mu),lty=1)
shade(sqrt(lambda_f.PI),minpress_seq)
## The dashed trend is “masculine” storms,and the solid trend is “feminine” storms. It is observed from the interaction model that feminine storms are a little more deadly, but the difference between masculine and feminine storms actually decreases as pressure drops.
#standardizepredictor
d$damage_std<-(d$damage_norm-mean(d$damage_norm))/sd(d$damage_norm)
#interactionmodel
m_f_d_fd<-map2stan(
alist(
deaths~dpois(lambda),
log(lambda)<-a+bf*fmnnty+bd*damage+
bfd*fmnnty*damage,
a~dnorm(0,10), c
(bf,bd,bfd)~dnorm(0,1) ),
data=list( deaths=d$deaths,
fmnnty=d$fmnnty_std,
damage=d$damage_std
),chains=4) ##
## SAMPLING FOR MODEL 'a9ccc9552fe4cbc419be6a4099896014' 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 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.275 seconds (Warm-up)
## Chain 1: 0.315 seconds (Sampling)
## Chain 1: 0.59 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'a9ccc9552fe4cbc419be6a4099896014' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.284 seconds (Warm-up)
## Chain 2: 0.242 seconds (Sampling)
## Chain 2: 0.526 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'a9ccc9552fe4cbc419be6a4099896014' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.288 seconds (Warm-up)
## Chain 3: 0.28 seconds (Sampling)
## Chain 3: 0.568 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'a9ccc9552fe4cbc419be6a4099896014' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.298 seconds (Warm-up)
## Chain 4: 0.291 seconds (Sampling)
## Chain 4: 0.589 seconds (Total)
## Chain 4:m_f_d<-map2stan(
alist(
deaths~dpois(lambda),
log(lambda)<-a+bf*fmnnty+bd*damage,
a~dnorm(0,10),
c(bf,bd)~dnorm(0,1) ),
data=list(
deaths=d$deaths,
fmnnty=d$fmnnty_std,
damage=d$damage_std
),chains=4)##
## SAMPLING FOR MODEL 'daf6193224979c17bf30447bceeb3e5e' 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 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.138 seconds (Warm-up)
## Chain 1: 0.152 seconds (Sampling)
## Chain 1: 0.29 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL 'daf6193224979c17bf30447bceeb3e5e' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.179 seconds (Warm-up)
## Chain 2: 0.118 seconds (Sampling)
## Chain 2: 0.297 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL 'daf6193224979c17bf30447bceeb3e5e' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.117 seconds (Warm-up)
## Chain 3: 0.151 seconds (Sampling)
## Chain 3: 0.268 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL 'daf6193224979c17bf30447bceeb3e5e' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.14 seconds (Warm-up)
## Chain 4: 0.13 seconds (Sampling)
## Chain 4: 0.27 seconds (Total)
## Chain 4:damage_seq<-seq(from=-0.6,to=5.5,length.out=30)
#'masculine'storms
d_pred<-data.frame(
fmnnty=-1, damage=damage_seq)
lambda_m<-link(m_f_d_fd,data=d_pred) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda_m.mu<-apply(lambda_m,2,mean)
lambda_m.PI<-apply(lambda_m,2,PI)
#'feminine'storms
d_pred<-data.frame(
fmnnty=1, damage=damage_seq)
lambda_f<-link(m_f_d_fd,data=d_pred) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda_f.mu<-apply(lambda_f,2,mean)
lambda_f.PI<-apply(lambda_f,2,PI)
#plot
plot(d$damage_std,sqrt(d$deaths),
pch=ifelse(d$fmnnty_std>0,16,1),col=rangi2,
xlab="damage",ylab="sqrt(deaths)")
lines(damage_seq,sqrt(lambda_m.mu),lty=2)
shade(sqrt(lambda_m.PI),damage_seq)
lines(damage_seq,sqrt(lambda_f.mu),lty=1)
shade(sqrt(lambda_f.PI),damage_seq)
## Now we see the relationship: feminine storms (solid trend) are stronger at all damage values and their difference from masculine storms increases with damage. Effect is very small.
precis(m_f_d_fd)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.81344145 0.02718093 2.77076424 2.85730782 2661.128 1.0008823
## bf 0.22112207 0.02816413 0.17645281 0.26703064 2058.891 0.9995264
## bd 0.45941501 0.01171186 0.44057092 0.47779184 2554.237 0.9998514
## bfd 0.03365944 0.01350270 0.01237077 0.05519579 2059.537 0.9994329## The interaction coefficient is quite small. So it doesn’t make much difference, until storm damage grows massive.
dat<-list(
deaths=d$deaths,
fmnnty=d$fmnnty_std,
damage=d$damage_std )
mgP_f_d_fd<-map2stan(
alist(
deaths~dgampois(lambda,scale),
log(lambda)<-a+bf*fmnnty+bd*damage+
bfd*fmnnty*damage,
a~dnorm(0,10),
c(bf,bd,bfd)~dnorm(0,1),
scale~dcauchy(0,1)
),
data=dat,chains=4) ##
## SAMPLING FOR MODEL '9fa6f38ed8149185859cdb9120d42983' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.001 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.735 seconds (Warm-up)
## Chain 1: 0.69 seconds (Sampling)
## Chain 1: 1.425 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '9fa6f38ed8149185859cdb9120d42983' 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 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.601 seconds (Warm-up)
## Chain 2: 0.522 seconds (Sampling)
## Chain 2: 1.123 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '9fa6f38ed8149185859cdb9120d42983' 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 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.624 seconds (Warm-up)
## Chain 3: 0.602 seconds (Sampling)
## Chain 3: 1.226 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '9fa6f38ed8149185859cdb9120d42983' 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 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.641 seconds (Warm-up)
## Chain 4: 0.628 seconds (Sampling)
## Chain 4: 1.269 seconds (Total)
## Chain 4:precis(mgP_f_d_fd)## mean sd 5.5% 94.5% n_eff Rhat4
## a 2.59785446 0.1318940 2.39451470 2.8137585 4613.554 0.9998397
## bf 0.08645912 0.1342637 -0.13558668 0.3012519 4489.632 1.0001499
## bd 1.25009896 0.2222657 0.90614297 1.6088969 4436.792 0.9995311
## bfd 0.30615618 0.2030470 -0.02672891 0.6155897 4459.075 0.9994560
## scale 0.68347467 0.1036640 0.52963090 0.8562883 4434.236 0.9995279## This model finds a stronger relationship with damage but estimates are more variable. This can be expected from a gamma-Poisson model. The gamma-distributed rates allow for a wider range of parameter values to produce similar predictions.
damage_seq<-seq(from=-0.6,to=5.5,length.out=30)
#'masculine'storms
d_pred<-data.frame( fmnnty=-1, damage=damage_seq)
lambda_m<-link(mgP_f_d_fd,data=d_pred) ## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda_m.mu<-apply(lambda_m,2,median)
lambda_m.PI<-apply(lambda_m,2,PI)
#'feminine'storms
d_pred<-data.frame( fmnnty=1, damage=damage_seq)
lambda_f<-link(mgP_f_d_fd,data=d_pred)## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]lambda_f.mu<-apply(lambda_f,2,median)
lambda_f.PI<-apply(lambda_f,2,PI)
#plot
plot(d$damage_std,sqrt(d$deaths),
pch=ifelse(d$fmnnty_std>0,16,1),col=rangi2,
xlab="damage",ylab="sqrt(deaths)")
lines(damage_seq,sqrt(lambda_m.mu),lty=2)
shade(sqrt(lambda_m.PI),damage_seq)
lines(damage_seq,sqrt(lambda_f.mu),lty=1)
shade(sqrt(lambda_f.PI),damage_seq)
## Overlapping prediction regions for masculine and feminine storms create that dark wedge in the middle. Feminine storms are more deadly at all damage amounts, and the difference grows with damage. Gamma-distributed variation is observed because predictions are all over the place, especially for masculine storms.