## R code 7.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)
data(rugged)
d <- rugged
head(d)
## isocode isonum country rugged rugged_popw rugged_slope rugged_lsd
## 1 ABW 533 Aruba 0.462 0.380 1.226 0.144
## 2 AFG 4 Afghanistan 2.518 1.469 7.414 0.720
## 3 AGO 24 Angola 0.858 0.714 2.274 0.228
## 4 AIA 660 Anguilla 0.013 0.010 0.026 0.006
## 5 ALB 8 Albania 3.427 1.597 10.451 1.006
## 6 AND 20 Andorra 5.717 6.722 17.774 1.616
## rugged_pc land_area lat lon soil desert tropical dist_coast
## 1 0.000 18 12.508 -69.970 21.324 0.000 100.000 0.001
## 2 39.004 65209 33.833 66.026 27.849 4.356 0.000 0.922
## 3 4.906 124670 -12.299 17.551 26.676 0.425 44.346 0.428
## 4 0.000 9 18.231 -63.064 100.000 0.000 100.000 0.000
## 5 62.133 2740 41.143 20.070 68.088 0.000 0.000 0.048
## 6 99.064 47 42.551 1.576 0.000 0.000 0.000 0.134
## near_coast gemstones rgdppc_2000 rgdppc_1950_m rgdppc_1975_m
## 1 100.0000 0 NA NA NA
## 2 0.0000 0 NA 644.756 720.633
## 3 13.1587 47756 1794.729 1051.822 1073.036
## 4 100.0000 0 NA NA NA
## 5 94.6919 0 3703.113 1001.339 2289.472
## 6 0.0000 0 NA NA NA
## rgdppc_2000_m rgdppc_1950_2000_m q_rule_law cont_africa cont_asia
## 1 NA NA NA 0 0
## 2 565.231 679.791 -1.687 0 1
## 3 765.215 1106.763 -1.567 1 0
## 4 NA NA NA 0 0
## 5 2741.420 1931.784 -0.820 0 0
## 6 NA NA 1.515 0 0
## cont_europe cont_oceania cont_north_america cont_south_america legor_gbr
## 1 0 0 1 0 0
## 2 0 0 0 0 0
## 3 0 0 0 0 0
## 4 0 0 1 0 NA
## 5 1 0 0 0 0
## 6 1 0 0 0 0
## legor_fra legor_soc legor_deu legor_sca colony_esp colony_gbr colony_fra
## 1 1 0 0 0 0 0 0
## 2 1 0 0 0 0 0 0
## 3 1 0 0 0 0 0 0
## 4 NA NA NA NA 0 0 0
## 5 0 1 0 0 0 0 0
## 6 1 0 0 0 0 0 0
## colony_prt colony_oeu africa_region_n africa_region_s africa_region_w
## 1 0 0 0 0 0
## 2 0 0 0 0 0
## 3 1 0 0 0 0
## 4 0 0 0 0 0
## 5 0 0 0 0 0
## 6 0 0 0 0 0
## africa_region_e africa_region_c slave_exports dist_slavemkt_atlantic
## 1 0 0 0 NA
## 2 0 0 0 NA
## 3 0 1 3610000 5.669
## 4 0 0 0 NA
## 5 0 0 0 NA
## 6 0 0 0 NA
## dist_slavemkt_indian dist_slavemkt_saharan dist_slavemkt_redsea pop_1400
## 1 NA NA NA 614
## 2 NA NA NA 1870829
## 3 6.981 4.926 3.872 1223208
## 4 NA NA NA NA
## 5 NA NA NA 200000
## 6 NA NA NA NA
## european_descent
## 1 NA
## 2 0
## 3 2
## 4 NA
## 5 100
## 6 NA
# make log version of outcome
d$log_gdp <- log( d$rgdppc_2000 )
# extract countries with GDP data
dd <- d[ complete.cases(d$rgdppc_2000) , ]
# split countries into Africa and not-Africa
d.A1 <- dd[ dd$cont_africa==1 , ] # Africa
d.A0 <- dd[ dd$cont_africa==0 , ] # not Africa
## R code 7.2
# African nations
m7.1 <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged ,
a ~ dnorm( 8 , 100 ) ,
bR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data=d.A1 )
# non-African nations
m7.2 <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged ,
a ~ dnorm( 8 , 100 ) ,
bR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data=d.A0 )
## R code 7.3
m7.3 <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged ,
a ~ dnorm( 8 , 100 ) ,
bR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data=dd )
## R code 7.4
m7.4 <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged + bA*cont_africa ,
a ~ dnorm( 8 , 100 ) ,
bR ~ dnorm( 0 , 1 ) ,
bA ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data=dd )
## R code 7.5
compare( m7.3 , m7.4 )
## WAIC pWAIC dWAIC weight SE dSE
## m7.4 476.6 4.6 0 1 15.36 NA
## m7.3 539.6 2.7 63 0 13.32 15.17
## R code 7.6
rugged.seq <- seq(from=-1,to=8,by=0.25)
# compute mu over samples, fixing cont_africa=0
mu.NotAfrica <- link( m7.4 ,
data=data.frame(cont_africa=0,
rugged=rugged.seq) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
# compute mu over samples, fixing cont_africa=1
mu.Africa <- link( m7.4 ,
data=data.frame(cont_africa=1,
rugged=rugged.seq) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
# summarize to means and intervals
mu.NotAfrica.mean <- apply( mu.NotAfrica , 2 , mean )
mu.NotAfrica.PI <- apply( mu.NotAfrica , 2 , PI ,
prob=0.97 )
mu.Africa.mean <- apply( mu.Africa , 2 , mean )
mu.Africa.PI <- apply( mu.Africa , 2 , PI , prob=0.97 )
## R code 7.7
m7.5 <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + gamma*rugged + bA*cont_africa ,
gamma <- bR + bAR*cont_africa ,
a ~ dnorm( 8 , 100 ) ,
bA ~ dnorm( 0 , 1 ) ,
bR ~ dnorm( 0 , 1 ) ,
bAR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data=dd )
## R code 7.8
compare( m7.3 , m7.4 , m7.5 )
## WAIC pWAIC dWAIC weight SE dSE
## m7.5 469.5 5.2 0.0 0.97 15.07 NA
## m7.4 476.4 4.4 7.0 0.03 15.38 6.36
## m7.3 539.7 2.8 70.3 0.00 13.32 15.17
## R code 7.9
m7.5b <- map(
alist(
log_gdp ~ dnorm( mu , sigma ) ,
mu <- a + bR*rugged + bAR*rugged*cont_africa +
bA*cont_africa,
a ~ dnorm( 8 , 100 ) ,
bA ~ dnorm( 0 , 1 ) ,
bR ~ dnorm( 0 , 1 ) ,
bAR ~ dnorm( 0 , 1 ) ,
sigma ~ dunif( 0 , 10 )
) ,
data=dd )
## R code 7.10
rugged.seq <- seq(from=-1,to=8,by=0.25)
mu.Africa <- link( m7.5 , data=data.frame(cont_africa=1,rugged=rugged.seq) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.Africa.mean <- apply( mu.Africa , 2 , mean )
mu.Africa.PI <- apply( mu.Africa , 2 , PI , prob=0.97 )
mu.NotAfrica <- link( m7.5 , data=data.frame(cont_africa=0,rugged=rugged.seq) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.NotAfrica.mean <- apply( mu.NotAfrica , 2 , mean )
mu.NotAfrica.PI <- apply( mu.NotAfrica , 2 , PI , prob=0.97 )
## R code 7.11
# plot African nations with regression
d.A1 <- dd[dd$cont_africa==1,]
plot( log(rgdppc_2000) ~ rugged , data=d.A1 ,
col=col.alpha(rangi2,0.7) , ylab="log GDP year 2000" ,
xlab="Terrain Ruggedness Index" )
mtext( "African nations" , 3 )
lines( rugged.seq , mu.Africa.mean , col=rangi2 )
shade( mu.Africa.PI , rugged.seq , col=col.alpha(rangi2,0.15) )
# plot non-African nations with regression
d.A0 <- dd[dd$cont_africa==0,]
plot( log(rgdppc_2000) ~ rugged , data=d.A0 ,
col=col.alpha("black",0.3) , ylab="log GDP year 2000" ,
xlab="Terrain Ruggedness Index" )
mtext( "Non-African nations" , 3 )
lines( rugged.seq , mu.NotAfrica.mean )
shade( mu.NotAfrica.PI , rugged.seq )
## R code 7.12
precis(m7.5)
## Mean StdDev 5.5% 94.5%
## a 9.18 0.14 8.97 9.40
## bA -1.85 0.22 -2.20 -1.50
## bR -0.18 0.08 -0.31 -0.06
## bAR 0.35 0.13 0.14 0.55
## sigma 0.93 0.05 0.85 1.01
## R code 7.13
post <- extract.samples( m7.5 )
gamma.Africa <- post$bR + post$bAR*1
gamma.notAfrica <- post$bR + post$bAR*0
## R code 7.14
mean( gamma.Africa)
## [1] 0.1632458
mean( gamma.notAfrica )
## [1] -0.1837559
## R code 7.15
dens( gamma.Africa , xlim=c(-0.5,0.6) , ylim=c(0,5.5) ,
xlab="gamma" , col=rangi2 )
dens( gamma.notAfrica , add=TRUE )
## R code 7.16
diff <- gamma.Africa - gamma.notAfrica
sum( diff < 0 ) / length( diff )
## [1] 0.0024
## R code 7.17
# get minimum and maximum rugged values
q.rugged <- range(dd$rugged)
# compute lines and confidence intervals
mu.ruggedlo <- link( m7.5 ,
data=data.frame(rugged=q.rugged[1],cont_africa=0:1) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.ruggedlo.mean <- apply( mu.ruggedlo , 2 , mean )
mu.ruggedlo.PI <- apply( mu.ruggedlo , 2 , PI )
mu.ruggedhi <- link( m7.5 ,
data=data.frame(rugged=q.rugged[2],cont_africa=0:1) )
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.ruggedhi.mean <- apply( mu.ruggedhi , 2 , mean )
mu.ruggedhi.PI <- apply( mu.ruggedhi , 2 , PI )
# plot it all, splitting points at median
med.r <- median(dd$rugged)
ox <- ifelse( dd$rugged > med.r , 0.05 , -0.05 )
plot( dd$cont_africa + ox , log(dd$rgdppc_2000) ,
col=ifelse(dd$rugged>med.r,rangi2,"black") ,
xlim=c(-0.25,1.25) , xaxt="n" , ylab="log GDP year 2000" ,
xlab="Continent" )
axis( 1 , at=c(0,1) , labels=c("other","Africa") )
lines( 0:1 , mu.ruggedlo.mean , lty=2 )
shade( mu.ruggedlo.PI , 0:1 )
lines( 0:1 , mu.ruggedhi.mean , col=rangi2 )
shade( mu.ruggedhi.PI , 0:1 , col=col.alpha(rangi2,0.25) )
## R code 7.18
library(rethinking)
data(tulips)
d <- tulips
head(d)
## bed water shade blooms
## 1 a 1 1 0.00
## 2 a 1 2 0.00
## 3 a 1 3 111.04
## 4 a 2 1 183.47
## 5 a 2 2 59.16
## 6 a 2 3 76.75
str(d)
## 'data.frame': 27 obs. of 4 variables:
## $ bed : Factor w/ 3 levels "a","b","c": 1 1 1 1 1 1 1 1 1 2 ...
## $ water : int 1 1 1 2 2 2 3 3 3 1 ...
## $ shade : int 1 2 3 1 2 3 1 2 3 1 ...
## $ blooms: num 0 0 111 183.5 59.2 ...
## R code 7.19
m7.6 <- map(
alist(
blooms ~ dnorm( mu , sigma ) ,
mu <- a + bW*water + bS*shade ,
a ~ dnorm( 0 , 100 ) ,
bW ~ dnorm( 0 , 100 ) ,
bS ~ dnorm( 0 , 100 ) ,
sigma ~ dunif( 0 , 100 )
) ,
data=d )
m7.7 <- map(
alist(
blooms ~ dnorm( mu , sigma ) ,
mu <- a + bW*water + bS*shade + bWS*water*shade ,
a ~ dnorm( 0 , 100 ) ,
bW ~ dnorm( 0 , 100 ) ,
bS ~ dnorm( 0 , 100 ) ,
bWS ~ dnorm( 0 , 100 ) ,
sigma ~ dunif( 0 , 100 )
) ,
data=d )
## R code 7.20
m7.6 <- map(
alist(
blooms ~ dnorm( mu , sigma ) ,
mu <- a + bW*water + bS*shade ,
a ~ dnorm( 0 , 100 ) ,
bW ~ dnorm( 0 , 100 ) ,
bS ~ dnorm( 0 , 100 ) ,
sigma ~ dunif( 0 , 100 )
) ,
data=d ,
method="Nelder-Mead" ,
control=list(maxit=1e4) )
m7.7 <- map(
alist(
blooms ~ dnorm( mu , sigma ) ,
mu <- a + bW*water + bS*shade + bWS*water*shade ,
a ~ dnorm( 0 , 100 ) ,
bW ~ dnorm( 0 , 100 ) ,
bS ~ dnorm( 0 , 100 ) ,
bWS ~ dnorm( 0 , 100 ) ,
sigma ~ dunif( 0 , 100 )
) ,
data=d ,
method="Nelder-Mead" ,
control=list(maxit=1e4) )
## R code 7.21
coeftab(m7.6,m7.7)
## m7.6 m7.7
## a 53.46 -37.13
## bW 76.33 134.56
## bS -38.88 13.01
## sigma 57.40 76.96
## bWS NA -31.4
## nobs 27 27
## R code 7.22
compare( m7.6 , m7.7 )
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## Warning in dnorm(blooms, mu, sigma, log = TRUE): 計算結果が NaN になりまし
## た
## WAIC pWAIC dWAIC weight SE
## m7.6 305.5 5.1 NaN NaN 8.77
## m7.7 NaN NaN NaN NaN NaN
## R code 7.23
d$shade.c <- d$shade - mean(d$shade)
d$water.c <- d$water - mean(d$water)
## R code 7.24
m7.8 <- map(
alist(
blooms ~ dnorm( mu , sigma ) ,
mu <- a + bW*water.c + bS*shade.c ,
a ~ dnorm( 130 , 100 ) ,
bW ~ dnorm( 0 , 100 ) ,
bS ~ dnorm( 0 , 100 ) ,
sigma ~ dunif( 0 , 100 )
) ,
data=d ,
start=list(a=mean(d$blooms),bW=0,bS=0,
sigma=sd(d$blooms)) )
m7.9 <- map(
alist(
blooms ~ dnorm( mu , sigma ) ,
mu <- a + bW*water.c + bS*shade.c +
bWS*water.c*shade.c ,
a ~ dnorm( 130 , 100 ) ,
bW ~ dnorm( 0 , 100 ) ,
bS ~ dnorm( 0 , 100 ) ,
bWS ~ dnorm( 0 , 100 ) ,
sigma ~ dunif( 0 , 100 )
) ,
data=d ,
start=list(a=mean(d$blooms),bW=0,bS=0,bWS=0,
sigma=sd(d$blooms)) )
coeftab(m7.8,m7.9)
## m7.8 m7.9
## a 129.00 129.01
## bW 74.22 74.96
## bS -40.74 -41.14
## sigma 57.35 45.22
## bWS NA -51.87
## nobs 27 27
## R code 7.25
k <- coef(m7.7)
k[1] + k[2]*2 + k[3]*2 + k[4]*2*2
## a
## 132.4133
## R code 7.26
k <- coef(m7.9)
k[1] + k[2]*0 + k[3]*0 + k[4]*0*0
## a
## 129.008
## R code 7.27
precis(m7.9)
## Mean StdDev 5.5% 94.5%
## a 129.01 8.67 115.15 142.87
## bW 74.96 10.60 58.02 91.90
## bS -41.14 10.60 -58.08 -24.20
## bWS -51.87 12.95 -72.57 -31.18
## sigma 45.22 6.15 35.39 55.06
## R code 7.28
# make a plot window with three panels in a single row
par(mfrow=c(1,3)) # 1 row, 3 columns
# loop over values of water.c and plot predictions
shade.seq <- -1:1
for ( w in -1:1 ) {
dt <- d[d$water.c==w,]
plot( blooms ~ shade.c , data=dt , col=rangi2 ,
main=paste("water.c =",w) , xaxp=c(-1,1,2) , ylim=c(0,362) ,
xlab="shade (centered)" )
mu <- link( m7.9 , data=data.frame(water.c=w,shade.c=shade.seq) )
mu.mean <- apply( mu , 2 , mean )
mu.PI <- apply( mu , 2 , PI , prob=0.97 )
lines( shade.seq , mu.mean )
lines( shade.seq , mu.PI[1,] , lty=2 )
lines( shade.seq , mu.PI[2,] , lty=2 )
}
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
## R code 7.29
m7.x <- lm( y ~ x + z + x*z , data=d )
summary(m7.x)
## R code 7.30
m7.x <- lm( y ~ x*z , data=d )
summary(m7.x)
## R code 7.31
m7.x <- lm( y ~ x + x*z - z , data=d )
summary(m7.x)
## R code 7.32
m7.x <- lm( y ~ x*z*w , data=d )
summary(m7.x)
## R code 7.33
x <- z <- w <- 1
colnames( model.matrix(~x*z*w) )
## R code 7.34
d$lang.per.cap <- d$num.lang / d$k.pop