## R code 9.1
n <- c( 12, 36 , 7 , 41 )
q <- n / sum(n)
q
## [1] 0.12500000 0.37500000 0.07291667 0.42708333
## R code 9.2
sum(q)
## [1] 1
## R code 9.3
p <- cumsum(q)
p
## [1] 0.1250000 0.5000000 0.5729167 1.0000000
## R code 9.4
log(p/(1-p))
## [1] -1.9459101 0.0000000 0.2937611 Inf
## R code 9.5
plot( 1:4 , p , xlab="rating" , ylab="cumulative proportion" ,
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" )
## R code 9.6
library(rethinking)
data(Hurricanes)
d <- Hurricanes
d$fmnnty_std <- ( d$femininity - mean(d$femininity) )/sd(d$femininity)
head(d)
## name year deaths category min_pressure damage_norm female femininity
## 1 Easy 1950 2 3 960 1590 1 6.77778
## 2 King 1950 4 3 955 5350 0 1.38889
## 3 Able 1952 3 1 985 150 0 3.83333
## 4 Barbara 1953 1 1 987 58 1 9.83333
## 5 Florence 1953 0 1 985 15 1 8.33333
## 6 Carol 1954 60 3 960 19321 1 8.11111
## fmnnty_std
## 1 -0.0009347453
## 2 -1.6707580424
## 3 -0.9133139565
## 4 0.9458703621
## 5 0.4810742824
## 6 0.4122162926
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 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.11 seconds (Warm-up)
## # 0.113 seconds (Sampling)
## # 0.223 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.105 seconds (Warm-up)
## # 0.122 seconds (Sampling)
## # 0.227 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.105 seconds (Warm-up)
## # 0.102 seconds (Sampling)
## # 0.207 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.105 seconds (Warm-up)
## # 0.107 seconds (Sampling)
## # 0.212 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.7
precis(m1)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 3.00 0.02 2.96 3.04 2297 1
## bf 0.24 0.03 0.20 0.28 2423 1
## R code 9.8
m0 <- map2stan(
alist(
deaths ~ dpois(lambda),
log(lambda) <- a,
a ~ dnorm(0,10)
),
data=list(deaths=d$deaths) , chains=4 )
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.06 seconds (Warm-up)
## # 0.058 seconds (Sampling)
## # 0.118 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.06 seconds (Warm-up)
## # 0.058 seconds (Sampling)
## # 0.118 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.06 seconds (Warm-up)
## # 0.055 seconds (Sampling)
## # 0.115 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.06 seconds (Warm-up)
## # 0.055 seconds (Sampling)
## # 0.115 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.9
compare(m0,m1)
## WAIC pWAIC dWAIC weight SE dSE
## m1 4436.4 140.1 0.0 1 1010.13 NA
## m0 4448.1 78.7 11.8 0 1075.11 136.03
## R code 9.10
# 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)
# superimpose 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 )
## R code 9.11
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 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.343 seconds (Warm-up)
## # 0.335 seconds (Sampling)
## # 0.678 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 1"
## count
## Exception thrown at line 20: neg_binomial_2_log: Location parameter[1] is 1.#INF, but must be finite 2
## Exception thrown at line 20: neg_binomial_2_log: Precision parameter is 1.#INF, but must be finite! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.337 seconds (Warm-up)
## # 0.343 seconds (Sampling)
## # 0.68 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 2"
## count
## Exception thrown at line 20: neg_binomial_2_log: Location parameter[1] is 1.#INF, but must be finite 2
## Exception thrown at line 20: neg_binomial_2_log: Location parameter[4] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.351 seconds (Warm-up)
## # 0.362 seconds (Sampling)
## # 0.713 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 3"
## count
## Exception thrown at line 20: neg_binomial_2_log: Location parameter[1] is 1.#INF, but must be finite 2
## Exception thrown at line 20: neg_binomial_2_log: Location parameter[4] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.353 seconds (Warm-up)
## # 0.359 seconds (Sampling)
## # 0.712 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 4"
## count
## Exception thrown at line 20: neg_binomial_2_log: Precision parameter is 1.#INF, but must be finite! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.12
precis(m3)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 3.02 0.16 2.76 3.27 2591 1
## bf 0.21 0.16 -0.03 0.46 2011 1
## scale 0.45 0.06 0.36 0.55 2821 1
## R code 9.13
plot(coeftab(m1,m3))
## R code 9.14
# 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)
# superimpose 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)
# 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 )
## R code 9.15
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) )
## R code 9.16
# 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 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.288 seconds (Warm-up)
## # 0.273 seconds (Sampling)
## # 0.561 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.28 seconds (Warm-up)
## # 0.275 seconds (Sampling)
## # 0.555 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.275 seconds (Warm-up)
## # 0.304 seconds (Sampling)
## # 0.579 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.277 seconds (Warm-up)
## # 0.286 seconds (Sampling)
## # 0.563 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.17
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 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.196 seconds (Warm-up)
## # 0.177 seconds (Sampling)
## # 0.373 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.186 seconds (Warm-up)
## # 0.194 seconds (Sampling)
## # 0.38 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.188 seconds (Warm-up)
## # 0.18 seconds (Sampling)
## # 0.368 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.192 seconds (Warm-up)
## # 0.183 seconds (Sampling)
## # 0.375 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.18
compare(m_f_p,m_f_p_fp)
## WAIC pWAIC dWAIC weight SE dSE
## m_f_p 3489.0 197.4 0.0 1 1117.08 NA
## m_f_p_fp 3569.2 260.1 80.3 0 1143.80 45.6
## R code 9.19
precis(m_f_p_fp)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 2.72 0.03 2.68 2.77 1689 1.00
## bf 0.26 0.03 0.21 0.31 1729 1.01
## bp -0.74 0.02 -0.77 -0.70 1756 1.00
## bfp 0.07 0.03 0.03 0.11 1788 1.00
## R code 9.20
minpress_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)
# now try plotting together
# will use sqrt scale for deaths,
# to make differences easier to see
# cannot use log scale, bc of zeros in data
# note uses of sqrt() throughout code
plot( d$min_pressure_std , sqrt(d$deaths) ,
pch=ifelse(d$fmnnty_std>0,16,1) , col=rangi2 ,
xlab="minimum pressure" , 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 )
## R code 9.21
# standardize predictor
d$damage_std <- (d$damage_norm - mean(d$damage_norm))/sd(d$damage_norm)
# interaction model
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 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.275 seconds (Warm-up)
## # 0.28 seconds (Sampling)
## # 0.555 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.265 seconds (Warm-up)
## # 0.261 seconds (Sampling)
## # 0.526 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.273 seconds (Warm-up)
## # 0.233 seconds (Sampling)
## # 0.506 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.277 seconds (Warm-up)
## # 0.316 seconds (Sampling)
## # 0.593 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
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 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.194 seconds (Warm-up)
## # 0.204 seconds (Sampling)
## # 0.398 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.187 seconds (Warm-up)
## # 0.186 seconds (Sampling)
## # 0.373 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.182 seconds (Warm-up)
## # 0.185 seconds (Sampling)
## # 0.367 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.2 seconds (Warm-up)
## # 0.201 seconds (Sampling)
## # 0.401 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.22
compare( m_f_d , m_f_d_fd )
## WAIC pWAIC dWAIC weight SE dSE
## m_f_d 3212.3 128.2 0.0 1 790.04 NA
## m_f_d_fd 3225.2 142.1 12.9 0 794.09 27.16
## R code 9.23
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 )
## R code 9.24
precis(m_f_d_fd)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 2.81 0.03 2.77 2.85 1750 1
## bf 0.22 0.03 0.18 0.27 1805 1
## bd 0.46 0.01 0.44 0.48 1985 1
## bfd 0.03 0.01 0.01 0.05 2050 1
## R code 9.25
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 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.649 seconds (Warm-up)
## # 0.667 seconds (Sampling)
## # 1.316 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 1"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[59] is 1.#INF, but must be finit 2
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[1] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.677 seconds (Warm-up)
## # 0.633 seconds (Sampling)
## # 1.31 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 2"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[6] is 1.#INF, but must be finite 2
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[16] is 1.#INF, but must be finit 1
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[4] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.637 seconds (Warm-up)
## # 0.613 seconds (Sampling)
## # 1.25 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 3"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[1] is 1.#INF, but must be finite 2
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[6] is 1.#INF, but must be finite 2
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[16] is 1.#INF, but must be finit 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.649 seconds (Warm-up)
## # 0.664 seconds (Sampling)
## # 1.313 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 4"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[6] is 1.#INF, but must be finite 3
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[1] is 0, but must be > 0! 1
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[16] is 1.#INF, but must be finit 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
precis(mgP_f_d_fd)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 2.60 0.13 2.39 2.81 2517 1
## bf 0.09 0.14 -0.12 0.31 2471 1
## bd 1.25 0.22 0.92 1.59 2333 1
## bfd 0.32 0.20 0.00 0.63 2310 1
## scale 0.68 0.11 0.52 0.86 1997 1
## R code 9.26
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 )
## R code 9.27
d$log_damage <- log(d$damage_norm)
d$log_damage_std <- (d$log_damage - mean(d$log_damage))/sd(d$log_damage)
dat <- list(
deaths=d$deaths,
fmnnty=d$fmnnty_std,
log_damage=d$log_damage_std )
# capital 'D' here denotes log scale
m_f_D_fD <- map2stan(
alist(
deaths ~ dpois(lambda),
log(lambda) <- a + bf*fmnnty + bD*log_damage +
bfD*fmnnty*log_damage,
a ~ dnorm(0,10),
c(bf,bD,bfD) ~ dnorm(0,1)
),
data=dat , chains=4 )
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.411 seconds (Warm-up)
## # 0.416 seconds (Sampling)
## # 0.827 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.414 seconds (Warm-up)
## # 0.452 seconds (Sampling)
## # 0.866 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.409 seconds (Warm-up)
## # 0.4 seconds (Sampling)
## # 0.809 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.392 seconds (Warm-up)
## # 0.408 seconds (Sampling)
## # 0.8 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
m_f_D <- map2stan(
alist(
deaths ~ dpois(lambda),
log(lambda) <- a + bf*fmnnty + bD*log_damage,
a ~ dnorm(0,10),
c(bf,bD) ~ dnorm(0,1)
),
data=dat , chains=4 )
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.275 seconds (Warm-up)
## # 0.28 seconds (Sampling)
## # 0.555 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.265 seconds (Warm-up)
## # 0.237 seconds (Sampling)
## # 0.502 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.264 seconds (Warm-up)
## # 0.253 seconds (Sampling)
## # 0.517 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (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)
## # Elapsed Time: 0.253 seconds (Warm-up)
## # 0.293 seconds (Sampling)
## # 0.546 seconds (Total)
##
##
## SAMPLING FOR MODEL 'deaths ~ dpois(lambda)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## R code 9.28
compare( m_f_d , m_f_d_fd , m_f_D , m_f_D_fD )
## WAIC pWAIC dWAIC weight SE dSE
## m_f_D_fD 2090.1 102.2 0.0 1 439.04 NA
## m_f_D 2127.4 95.2 37.3 0 447.87 47.61
## m_f_d 3212.3 128.2 1122.2 0 790.04 430.27
## m_f_d_fd 3225.2 142.1 1135.1 0 794.09 434.27
## R code 9.29
precis(m_f_D_fD)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 2.18 0.04 2.11 2.25 1528 1
## bf 0.01 0.04 -0.06 0.08 1313 1
## bD 1.51 0.04 1.45 1.57 1578 1
## bfD 0.30 0.04 0.23 0.36 1312 1
## R code 9.30
log_damage_seq <- seq(from=-3.2,to=2,length.out=30)
# 'masculine' storms
d_pred <- data.frame(
fmnnty=-1,
log_damage=log_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,median)
lambda_m.PI <- apply(lambda_m,2,PI)
# 'feminine' storms
d_pred <- data.frame(
fmnnty=1,
log_damage=log_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,median)
lambda_f.PI <- apply(lambda_f,2,PI)
# plot
plot( d$log_damage_std , sqrt(d$deaths) ,
pch=ifelse(d$fmnnty_std>0,16,1) , col=rangi2 ,
xlab="log damage (std)" , ylab="sqrt(deaths)" )
lines( log_damage_seq , sqrt(lambda_m.mu) , lty=2 )
shade( sqrt(lambda_m.PI) , log_damage_seq )
lines( log_damage_seq , sqrt(lambda_f.mu) , lty=1 )
shade( sqrt(lambda_f.PI) , log_damage_seq )
## R code 9.31
mgP_f_D_fD <- map2stan(
alist(
deaths ~ dgampois(lambda,scale),
log(lambda) <- a + bf*fmnnty + bD*log_damage +
bfD*fmnnty*log_damage,
a ~ dnorm(0,10),
c(bf,bD,bfD) ~ dnorm(0,1),
scale ~ dcauchy(0,1)
),
data=dat , chains=4 )
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.646 seconds (Warm-up)
## # 0.633 seconds (Sampling)
## # 1.279 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 1"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[6] is 1.#INF, but must be finite 2
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[2] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.659 seconds (Warm-up)
## # 0.572 seconds (Sampling)
## # 1.231 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 2"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[39] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.656 seconds (Warm-up)
## # 0.6 seconds (Sampling)
## # 1.256 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 3"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[1] is 0, but must be > 0! 1
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[16] is 0, but must be > 0! 1
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[8] is 0, but must be > 0! 1
## Exception thrown at line 25: neg_binomial_2_log: Precision parameter is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (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)
## # Elapsed Time: 0.622 seconds (Warm-up)
## # 0.683 seconds (Sampling)
## # 1.305 seconds (Total)
##
## [1] "The following numerical problems occured the indicated number of times on chain 4"
## count
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[5] is 0, but must be > 0! 2
## Exception thrown at line 25: neg_binomial_2_log: Location parameter[6] is 0, but must be > 0! 1
## [1] "If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,"
## [2] "but if this warning occurs often then your model may be either severely ill-conditioned or misspecified."
##
## SAMPLING FOR MODEL 'deaths ~ dgampois(lambda, scale)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0 seconds (Sampling)
## # 0 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
precis(mgP_f_D_fD)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 2.25 0.12 2.05 2.42 2308 1
## bf 0.04 0.12 -0.15 0.22 2338 1
## bD 1.37 0.13 1.16 1.57 2429 1
## bfD 0.17 0.12 -0.03 0.36 2864 1
## scale 1.06 0.18 0.78 1.33 2569 1
## R code 9.32
# 'masculine' storms
d_pred <- data.frame(
fmnnty=-1,
log_damage=log_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,
log_damage=log_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$log_damage_std , sqrt(d$deaths) ,
pch=ifelse(d$fmnnty_std>0,16,1) , col=rangi2 ,
xlab="log damage (std)" , ylab="sqrt(deaths)" )
lines( log_damage_seq , sqrt(lambda_m.mu) , lty=2 )
shade( sqrt(lambda_m.PI) , log_damage_seq )
lines( log_damage_seq , sqrt(lambda_f.mu) , lty=1 )
shade( sqrt(lambda_f.PI) , log_damage_seq )
## R code 9.33
library(rethinking)
data(Trolley)
d <- Trolley
d$female <- 1 - d$male
## R code 9.34
m1a <- map(
alist(
response ~ dordlogit( phi , c(a1,a2,a3,a4,a5,a6) ),
phi <- bA*action + bI*intention + bC*contact +
bAI*action*intention + bCI*contact*intention +
bf*female + bfC*female*contact,
c(a1,a2,a3,a4,a5,a6) ~ dnorm(0,10),
c(bA,bI,bC,bAI,bCI,bf,bfC) ~ dnorm(0,10)
) ,
data=d ,
start=list(a1=-2,a2=-1.5,a3=-0.5,a4=0,a5=1,a6=1.5) )
## R code 9.35
dat <- list(
response=d$response,
action=d$action,
intention=d$intention,
contact=d$contact,
female=d$female
)
m1a_stan <- map2stan(
alist(
response ~ dordlogit( phi , cutpoints ),
phi <- bA*action + bI*intention + bC*contact +
bAI*action*intention + bCI*contact*intention +
bf*female + bfC*female*contact,
cutpoints ~ dnorm(0,10),
c(bA,bI,bC,bAI,bCI,bf,bfC) ~ dnorm(0,10)
) ,
data=dat , chains=3 , cores=3 ,
start=list(cutpoints=c(-2,-1.5,-0.5,0,1,1.5)) )
## Precompiling model...
## the number of chains is less than 1; sampling not done
## Dispatching 3 chains to 3 cores.
##
## SAMPLING FOR MODEL 'response ~ dordlogit(phi, cutpoints)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0.047 seconds (Sampling)
## # 0.047 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## Warning in array(0, dim = c(n, n_cases, K)): Reached total allocation of
## 1535Mb: see help(memory.size)
## Warning in array(0, dim = c(n, n_cases, K)): Reached total allocation of
## 1535Mb: see help(memory.size)
## Warning in array(0, dim = c(n, n_cases, K)): Reached total allocation of
## 1535Mb: see help(memory.size)
## Warning in array(0, dim = c(n, n_cases, K)): Reached total allocation of
## 1535Mb: see help(memory.size)
## R code 9.36
precis( m1a )
## Mean StdDev 5.5% 94.5%
## a1 -2.95 0.06 -3.04 -2.86
## a2 -2.25 0.05 -2.34 -2.17
## a3 -1.65 0.05 -1.73 -1.57
## a4 -0.59 0.05 -0.67 -0.52
## a5 0.10 0.05 0.02 0.17
## a6 1.02 0.05 0.95 1.10
## bA -0.48 0.05 -0.56 -0.39
## bI -0.28 0.06 -0.38 -0.19
## bC -0.43 0.08 -0.56 -0.30
## bAI -0.45 0.08 -0.58 -0.32
## bCI -1.30 0.10 -1.46 -1.14
## bf -0.62 0.04 -0.68 -0.55
## bfC 0.21 0.09 0.07 0.35
## R code 9.37
post <- extract.samples(m1a)
phi.link <- function(A,I,C,f) {
post$bA*A + post$bI*I + post$bC*C + post$bAI*A*I +
post$bCI*C*I + post$bf*f + post$bfC*f*C
}
## R code 9.38
female.noC <- phi.link(0,0,0,1)
female.C <- phi.link(0,0,1,1)
male.noC <- phi.link(0,0,0,0)
male.C <- phi.link(0,0,1,0)
## R code 9.39
female.diff <- female.C - female.noC
male.diff <- male.C - male.noC
## R code 9.40
precis( data.frame( female.diff , male.diff ) )
## Mean StdDev |0.89 0.89|
## female.diff -0.22 0.08 -0.36 -0.09
## male.diff -0.43 0.08 -0.56 -0.29
## R code 9.41
post <- extract.samples(m1a)
kA <- 0
kC <- 1
kI <- 1
kf <- 0:1
plot( 1 , 1 , type="n" , xlab="female" , ylab="probability" ,
xlim=c(0,1) , ylim=c(0,1) , xaxp=c(0,1,1) , yaxp=c(0,1,2) )
for ( s in 1:100 ) {
p <- post[s,]
ak <- as.numeric( c(p$a1,p$a2,p$a3,p$a4,p$a5,p$a6) )
phi <- p$bf*kf + p$bfC*kf*kC + p$bA*kA + p$bI*kI +
p$bC*kC + p$bAI*kA*kI + p$bCI*kC*kI
pk <- pordlogit( 1:6 , phi=phi , a=ak )
for ( i in 1:6 )
lines( 0:1 , pk[,i] ,
col=col.alpha("slateblue",0.1) , lwd=0.5 )
}
abline( h=0 , lty=2 , col="lightgray" )
abline( h=1 , lty=2 , col="lightgray" )
mtext( paste("A =",kA,", I =",kI,", C =",kC ) , 3 )
## R code 9.42
library(rethinking)
data(Fish)
d <- Fish
str(d)
## 'data.frame': 250 obs. of 6 variables:
## $ fish_caught: int 0 0 0 0 1 0 0 0 0 1 ...
## $ livebait : int 0 1 1 1 1 1 1 1 0 1 ...
## $ camper : int 0 1 0 1 0 1 0 0 1 1 ...
## $ persons : int 1 1 1 2 1 4 3 4 3 1 ...
## $ child : int 0 0 0 1 0 2 1 3 2 0 ...
## $ hours : num 21.124 5.732 1.323 0.548 1.695 ...
## R code 9.43
d$loghours <- log(d$hours)
m2a <- map2stan(
alist(
fish_caught ~ dzipois(p,mu),
logit(p) <- a0 + bp0*persons + bc0*child,
log(mu) <- a + bp*persons + bc*child + loghours,
c(a0,a) ~ dnorm(0,10),
c(bp0,bc0,bp,bc) ~ dnorm(0,1)
) ,
data=d , iter=1e4 , chains=2 )
## Warning: 使われていないコネクション 84 (<-kuma:11485) を閉じます
## Warning: 使われていないコネクション 83 (<-kuma:11485) を閉じます
## Warning: 使われていないコネクション 82 (<-kuma:11485) を閉じます
##
## SAMPLING FOR MODEL 'fish_caught ~ dzipois(p, mu)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 1, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 1, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 1, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 1, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 1, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 1, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 1, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 1, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 1, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 1, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 1, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 41.364 seconds (Warm-up)
## # 51.536 seconds (Sampling)
## # 92.9 seconds (Total)
##
##
## SAMPLING FOR MODEL 'fish_caught ~ dzipois(p, mu)' NOW (CHAIN 2).
##
## Chain 2, Iteration: 1 / 10000 [ 0%] (Warmup)
## Chain 2, Iteration: 1000 / 10000 [ 10%] (Warmup)
## Chain 2, Iteration: 2000 / 10000 [ 20%] (Warmup)
## Chain 2, Iteration: 3000 / 10000 [ 30%] (Warmup)
## Chain 2, Iteration: 4000 / 10000 [ 40%] (Warmup)
## Chain 2, Iteration: 5000 / 10000 [ 50%] (Warmup)
## Chain 2, Iteration: 5001 / 10000 [ 50%] (Sampling)
## Chain 2, Iteration: 6000 / 10000 [ 60%] (Sampling)
## Chain 2, Iteration: 7000 / 10000 [ 70%] (Sampling)
## Chain 2, Iteration: 8000 / 10000 [ 80%] (Sampling)
## Chain 2, Iteration: 9000 / 10000 [ 90%] (Sampling)
## Chain 2, Iteration: 10000 / 10000 [100%] (Sampling)
## # Elapsed Time: 39.953 seconds (Warm-up)
## # 45.788 seconds (Sampling)
## # 85.741 seconds (Total)
##
##
## SAMPLING FOR MODEL 'fish_caught ~ dzipois(p, mu)' NOW (CHAIN 1).
##
## Chain 1, Iteration: 1 / 1 [100%] (Sampling)
## # Elapsed Time: 0 seconds (Warm-up)
## # 0.002 seconds (Sampling)
## # 0.002 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]
## R code 9.44
precis(m2a)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a0 0.76 0.56 -0.12 1.64 4076 1
## a -2.30 0.15 -2.52 -2.06 3862 1
## bp0 -1.00 0.27 -1.44 -0.58 3807 1
## bc0 1.00 0.57 0.09 1.84 3932 1
## bp 0.67 0.04 0.61 0.74 3847 1
## bc 0.56 0.09 0.42 0.71 5141 1
## R code 9.45
zip_link <- link(m2a)
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
str(zip_link)
## List of 2
## $ p : num [1:1000, 1:250] 0.526 0.511 0.438 0.427 0.186 ...
## $ mu: num [1:1000, 1:250] 4.3 4.96 4.11 4.6 4.21 ...
## R code 9.46
zeros <- rbinom(1e4,1,0.5)
obs_fish <- (1-zeros)*rpois(1e4,1)
simplehist(obs_fish)
## R code 9.47
fish_sim <- sim(m2a)
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
str(fish_sim)
## num [1:1000, 1:250] 0 0 1 0 4 9 0 2 3 0 ...
## R code 9.48
# new data
pred_dat <- list(
loghours=log(1), # note that this is zero, the baseline rate
persons=1,
child=0 )
# sim predictions - want expected number of fish, but must use both processes
fish_link <- link( m2a , data=pred_dat )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
# summarize
p <- fish_link$p
mu <- fish_link$mu
( expected_fish_mean <- mean( (1-p)*mu ) )
## [1] 0.1104318
( expected_fish_PI <- PI( (1-p)*mu ) )
## 5% 94%
## 0.08262274 0.14021904
