## R code 6.1
library(rethinking)
## Loading required package: rstan
## Loading required package: Rcpp
## Loading required package: ggplot2
## rstan (Version 2.8.0, packaged: 2015-09-19 14:48:38 UTC, GitRev: 05c3d0058b6a)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())
## Loading required package: parallel
## rethinking (Version 1.57)
sppnames <- c( "afarensis","africanus","habilis",
"boisei",
"rudolfensis","ergaster","sapiens")
brainvolcc <- c( 438 , 452 , 612, 521, 752, 871, 1350 )
masskg <- c( 37.0 , 35.5 , 34.5 , 41.5 , 55.5 , 61.0 ,
53.5 )
d <- data.frame( species=sppnames , brain=brainvolcc ,
mass=masskg )
d
## species brain mass
## 1 afarensis 438 37.0
## 2 africanus 452 35.5
## 3 habilis 612 34.5
## 4 boisei 521 41.5
## 5 rudolfensis 752 55.5
## 6 ergaster 871 61.0
## 7 sapiens 1350 53.5
## R code 6.2
m6.1 <- lm( brain ~ mass , data=d )
summary(m6.1)
##
## Call:
## lm(formula = brain ~ mass, data = d)
##
## Residuals:
## 1 2 3 4 5 6 7
## -99.86 -54.83 125.86 -109.96 -168.60 -163.39 470.77
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -227.629 439.794 -0.518 0.6268
## mass 20.689 9.436 2.192 0.0798 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 252.1 on 5 degrees of freedom
## Multiple R-squared: 0.4902, Adjusted R-squared: 0.3882
## F-statistic: 4.807 on 1 and 5 DF, p-value: 0.07985
## R code 6.3
1 - var(resid(m6.1))/var(d$brain)
## [1] 0.490158
## R code 6.4
m6.2 <- lm( brain ~ mass + I(mass^2) , data=d )
summary(m6.2)
##
## Call:
## lm(formula = brain ~ mass + I(mass^2), data = d)
##
## Residuals:
## 1 2 3 4 5 6 7
## -103.43 -21.66 186.35 -193.16 -206.14 -61.50 399.53
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2618.137 3831.633 -0.683 0.532
## mass 127.306 169.905 0.749 0.495
## I(mass^2) -1.133 1.802 -0.629 0.564
##
## Residual standard error: 268.8 on 4 degrees of freedom
## Multiple R-squared: 0.536, Adjusted R-squared: 0.304
## F-statistic: 2.31 on 2 and 4 DF, p-value: 0.2153
## R code 6.5
m6.3 <- lm( brain ~ mass + I(mass^2) + I(mass^3) , data=d )
summary(m6.3)
##
## Call:
## lm(formula = brain ~ mass + I(mass^2) + I(mass^3), data = d)
##
## Residuals:
## 1 2 3 4 5 6 7
## -31.26 -50.07 68.92 -12.22 -320.11 51.46 293.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21990.3741 21519.7671 1.022 0.382
## mass -1473.8357 1389.1938 -1.061 0.367
## I(mass^2) 32.8102 29.2975 1.120 0.344
## I(mass^3) -0.2351 0.2025 -1.161 0.330
##
## Residual standard error: 257.9 on 3 degrees of freedom
## Multiple R-squared: 0.6798, Adjusted R-squared: 0.3595
## F-statistic: 2.123 on 3 and 3 DF, p-value: 0.2761
m6.4 <- lm( brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) ,
data=d )
summary(m6.4)
##
## Call:
## lm(formula = brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4),
## data = d)
##
## Residuals:
## 1 2 3 4 5 6 7
## 139.072 1.863 -63.331 -111.818 -193.166 22.594 204.785
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.229e+05 2.506e+05 1.289 0.326
## mass -2.795e+04 2.201e+04 -1.270 0.332
## I(mass^2) 8.922e+02 7.139e+02 1.250 0.338
## I(mass^3) -1.243e+01 1.013e+01 -1.228 0.344
## I(mass^4) 6.393e-02 5.307e-02 1.205 0.352
##
## Residual standard error: 240.4 on 2 degrees of freedom
## Multiple R-squared: 0.8144, Adjusted R-squared: 0.4433
## F-statistic: 2.194 on 4 and 2 DF, p-value: 0.3367
m6.5 <- lm( brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) +
I(mass^5) , data=d )
summary(m6.5)
##
## Call:
## lm(formula = brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) +
## I(mass^5), data = d)
##
## Residuals:
## 1 2 3 4 5 6 7
## 43.0291 -64.3854 28.7485 -8.5915 -4.4456 0.4078 5.2371
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.535e+06 4.777e+05 -3.214 0.192
## mass 1.800e+05 5.313e+04 3.389 0.183
## I(mass^2) -8.325e+03 2.343e+03 -3.553 0.175
## I(mass^3) 1.896e+02 5.120e+01 3.704 0.168
## I(mass^4) -2.127e+00 5.541e-01 -3.838 0.162
## I(mass^5) 9.397e-03 2.376e-03 3.956 0.158
##
## Residual standard error: 83.33 on 1 degrees of freedom
## Multiple R-squared: 0.9889, Adjusted R-squared: 0.9331
## F-statistic: 17.74 on 5 and 1 DF, p-value: 0.1782
m6.6 <- lm( brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) +
I(mass^5) + I(mass^6) , data=d )
summary(m6.5)
##
## Call:
## lm(formula = brain ~ mass + I(mass^2) + I(mass^3) + I(mass^4) +
## I(mass^5), data = d)
##
## Residuals:
## 1 2 3 4 5 6 7
## 43.0291 -64.3854 28.7485 -8.5915 -4.4456 0.4078 5.2371
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.535e+06 4.777e+05 -3.214 0.192
## mass 1.800e+05 5.313e+04 3.389 0.183
## I(mass^2) -8.325e+03 2.343e+03 -3.553 0.175
## I(mass^3) 1.896e+02 5.120e+01 3.704 0.168
## I(mass^4) -2.127e+00 5.541e-01 -3.838 0.162
## I(mass^5) 9.397e-03 2.376e-03 3.956 0.158
##
## Residual standard error: 83.33 on 1 degrees of freedom
## Multiple R-squared: 0.9889, Adjusted R-squared: 0.9331
## F-statistic: 17.74 on 5 and 1 DF, p-value: 0.1782
## R code 6.6
m6.7 <- lm( brain ~ 1 , data=d )
summary(m6.7)
##
## Call:
## lm(formula = brain ~ 1, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -275.71 -227.21 -101.71 97.79 636.29
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 713.7 121.8 5.86 0.00109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 322.2 on 6 degrees of freedom
## R code 6.7
d.new <- d[ -1 , ]
## R code 6.8
plot( brain ~ mass , d , col="slateblue" )
for ( i in 1:nrow(d) ) {
d.new <- d[ -i , ]
m0 <- lm( brain ~ mass, d.new )
abline( m0 , col=col.alpha("black",0.5) )
}
## R code 6.9
p <- c( 0.3 , 0.7 )
-sum( p*log(p) )
## [1] 0.6108643
## R code 6.10
# fit model with lm
m6.1 <- lm( brain ~ mass , d )
# compute deviance by cheating
(-2) * logLik(m6.1)
## 'log Lik.' 94.92499 (df=3)
## R code 6.11
# standardize the mass before fitting
d$mass.s <- (d$mass-mean(d$mass))/sd(d$mass)
m6.8 <- map(
alist(
brain ~ dnorm( mu , sigma ) ,
mu <- a + b*mass.s
) ,
data=d ,
start=list(a=mean(d$brain),b=0,sigma=sd(d$brain)) ,
method="Nelder-Mead" )
# extract MAP estimates
theta <- coef(m6.8)
# compute deviance
dev <- (-2)*sum( dnorm(
d$brain ,
mean=theta[1]+theta[2]*d$mass.s ,
sd=theta[3] ,
log=TRUE ) )
dev
## [1] 94.92499
## R code 6.12
N <- 20
kseq <- 1:5
dev <- sapply( kseq , function(k) {
print(k);
r <- replicate( 1e4 , sim.train.test( N=N, k=k ) );
c( mean(r[1,]) , mean(r[2,]) , sd(r[1,]) , sd(r[2,]) )
} )
## R code 6.13
r <- mcreplicate( 1e4 , sim.train.test( N=N, k=k ) ,
mc.cores=4 )
## R code 6.14
plot( 1:5 , dev[1,] , ylim=c( min(dev[1:2,])-5 ,
max(dev[1:2,])+10 ) ,
xlim=c(1,5.1) , xlab="number of parameters" ,
ylab="deviance" ,
pch=16 , col=rangi2 )
mtext( concat( "N = ",N ) )
points( (1:5)+0.1 , dev[2,] )
for ( i in kseq ) {
pts_in <- dev[1,i] + c(-1,+1)*dev[3,i]
pts_out <- dev[2,i] + c(-1,+1)*dev[4,i]
lines( c(i,i) , pts_in , col=rangi2 )
lines( c(i,i)+0.1 , pts_out )
}
## R code 6.15
data(cars)
m <- map(
alist(
dist ~ dnorm(mu,sigma),
mu <- a + b*speed,
a ~ dnorm(0,100),
b ~ dnorm(0,10),
sigma ~ dunif(0,30)
) , data=cars )
post <- extract.samples(m,n=1000)
## R code 6.16
n_samples <- 1000
ll <- sapply( 1:n_samples ,
function(s) {
mu <- post$a[s] + post$b[s]*cars$speed
dnorm( cars$dist , mu , post$sigma[s] , log=TRUE )
} )
## R code 6.17
n_cases <- nrow(cars)
lppd <- sapply( 1:n_cases , function(i)
log_sum_exp(ll[i,]) - log(n_samples) )
## R code 6.18
pWAIC <- sapply( 1:n_cases , function(i) var(ll[i,]) )
## R code 6.19
-2*( sum(lppd) - sum(pWAIC) )
## [1] 420.362
## R code 6.20
waic_vec <- -2*( lppd - pWAIC )
sqrt( n_cases*var(waic_vec) )
## [1] 13.99977
## R code 6.21
data(milk)
d <- milk[ complete.cases(milk) , ]
d$neocortex <- d$neocortex.perc / 100
dim(d)
## [1] 17 9
## R code 6.22
a.start <- mean(d$kcal.per.g)
sigma.start <- log(sd(d$kcal.per.g))
m6.11 <- map(
alist(
kcal.per.g ~ dnorm( a , exp(log.sigma) )
) ,
data=d , start=list(a=a.start,log.sigma=sigma.start) )
m6.12 <- map(
alist(
kcal.per.g ~ dnorm( mu , exp(log.sigma) ) ,
mu <- a + bn*neocortex
) ,
data=d , start=list(a=a.start,bn=0,log.sigma=sigma.start) )
m6.13 <- map(
alist(
kcal.per.g ~ dnorm( mu , exp(log.sigma) ) ,
mu <- a + bm*log(mass)
) ,
data=d , start=list(a=a.start,bm=0,log.sigma=sigma.start) )
m6.14 <- map(
alist(
kcal.per.g ~ dnorm( mu , exp(log.sigma) ) ,
mu <- a + bn*neocortex + bm*log(mass)
) ,
data=d , start=list(a=a.start,bn=0,bm=0,log.sigma=sigma.start) )
## R code 6.23
WAIC( m6.14 )
## Constructing posterior predictions
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## [1] -15.45523
## attr(,"lppd")
## [1] 12.3435
## attr(,"pWAIC")
## [1] 4.61588
## attr(,"se")
## [1] 7.278786
## R code 6.24
( milk.models <- compare( m6.11 , m6.12 , m6.13 , m6.14 ) )
## WAIC pWAIC dWAIC weight SE dSE
## m6.14 -15.5 4.6 0.0 0.95 7.66 NA
## m6.11 -8.3 1.8 7.2 0.03 4.44 7.30
## m6.13 -7.6 3.2 8.0 0.02 5.86 5.35
## m6.12 -6.3 2.9 9.2 0.01 4.31 7.68
## R code 6.25
plot( milk.models , SE=TRUE , dSE=TRUE )
## R code 6.26
diff <- rnorm( 1e5 , 6.7 , 7.26 )
sum(diff<0)/1e5
## [1] 0.17844
## R code 6.27
coeftab(m6.11,m6.12,m6.13,m6.14)
## m6.11 m6.12 m6.13 m6.14
## a 0.66 0.35 0.71 -1.09
## log.sigma -1.79 -1.80 -1.85 -2.16
## bn NA 0.45 NA 2.79
## bm NA NA -0.03 -0.10
## nobs 17 17 17 17
## R code 6.28
plot( coeftab(m6.11,m6.12,m6.13,m6.14) )
## R code 6.29
# compute counterfactual predictions
# neocortex from 0.5 to 0.8
nc.seq <- seq(from=0.5,to=0.8,length.out=30)
d.predict <- list(
kcal.per.g = rep(0,30), # empty outcome
neocortex = nc.seq, # sequence of neocortex
mass = rep(4.5,30) # average mass
)
pred.m6.14 <- link( m6.14 , data=d.predict )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu <- apply( pred.m6.14 , 2 , mean )
mu.PI <- apply( pred.m6.14 , 2 , PI )
# plot it all
plot( kcal.per.g ~ neocortex , d , col=rangi2 )
lines( nc.seq , mu , lty=2 )
lines( nc.seq , mu.PI[1,] , lty=2 )
lines( nc.seq , mu.PI[2,] , lty=2 )
## R code 6.30
milk.ensemble <- ensemble( m6.11 , m6.12 , m6.13 , m6.14 , data=d.predict )
## Constructing posterior predictions
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## Constructing posterior predictions
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## Constructing posterior predictions
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## Constructing posterior predictions
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu <- apply( milk.ensemble$link , 2 , mean )
mu.PI <- apply( milk.ensemble$link , 2 , PI )
lines( nc.seq , mu )
shade( mu.PI , nc.seq )
## R code 6.31
library(rethinking)
data(Howell1)
d <- Howell1
d$age <- (d$age - mean(d$age))/sd(d$age)
head(d)
## height weight age male
## 1 151.765 47.82561 1.6222002 1
## 2 139.700 36.48581 1.6222002 0
## 3 136.525 31.86484 1.7186002 0
## 4 156.845 53.04191 0.5618002 1
## 5 145.415 41.27687 1.0438002 0
## 6 163.830 62.99259 0.2726002 1
set.seed( 1000 )
i <- sample(1:nrow(d),size=nrow(d)/2)
d1 <- d[ i , ]
d2 <- d[ -i , ]
## R code 6.32
#sum( dnorm( d2$height , mu , sigma , log=TRUE ) )