Chapter5

Hard

5H1

library(rethinking)
data("foxes")
colnames(foxes)

## [1] "group"     "avgfood"   "groupsize" "area"      "weight"

foxes.sc = cbind(foxes[,1,drop=F],apply(foxes[,-1],2,scale) )

H5.1_area <- map(
  alist(
    weight ~ dnorm( mu , sigma ) ,
    mu <- a + ba*area ,
    a ~ dnorm( 0 , 10 ) ,
    ba ~ dnorm( 0 , 1 ) ,
    sigma ~ dunif( 0 , 10 )
  ) ,
  data = foxes.sc )
plot(precis(H5.1_area))

# prepare new counterfactual data
# W.avg <- mean( foxes.sc$weight )
A.seq <- seq( from=-3 , to=3 , length.out=30 )
pred.data <- data.frame(
  area=A.seq#,
#  MedianWeigth=W.avg
)

# compute counterfactual mean divorce (mu)
mu <- link( H5.1_area , data=pred.data )

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

mu.mean <- apply( mu , 2 , mean )
mu.PI <- apply( mu , 2 , PI )

# simulate counterfactual divorce outcomes
R.sim <- sim( H5.1_area , data=pred.data , n=1e4 )

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

R.PI <- apply( R.sim , 2 , PI )

# display predictions, hiding raw data with type="n"
plot( weight ~ area , data=foxes.sc , type="n" )
lines( A.seq , mu.mean )
shade( mu.PI , A.seq )
shade( R.PI , A.seq )

H5.1_groupsize <- map(
  alist(
    weight ~ dnorm( mu , sigma ) ,
    mu <- a + bs*groupsize ,
    a ~ dnorm( 0 , 10 ) ,
    bs ~ dnorm( 0 , 1 ) ,
    sigma ~ dunif( 0 , 10 )
  ) ,
  data = foxes.sc )
plot(precis(H5.1_groupsize))

# prepare new counterfactual data
# W.avg <- mean( foxes.sc$weight )
GS.seq <- seq( from=-3 , to=3 , length.out=30 )
pred.data <- data.frame(
  area=GS.seq
)

# compute counterfactual mean divorce (mu)
mu <- link( H5.1_area , data=pred.data )

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

mu.mean <- apply( mu , 2 , mean )
mu.PI <- apply( mu , 2 , PI )

# simulate counterfactual divorce outcomes
R.sim <- sim( H5.1_area , data=pred.data , n=1e4 )

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

R.PI <- apply( R.sim , 2 , PI )

# display predictions, hiding raw data with type="n"
plot( weight ~ groupsize , data=foxes.sc , type="n" )
lines( GS.seq , mu.mean )
shade( mu.PI , GS.seq )
shade( R.PI , GS.seq )

5H2

H5.2 <- map(
  alist(
    weight ~ dnorm( mu , sigma ) ,
    mu <- a + bs*groupsize + ba*area ,
    a ~ dnorm( 0 , 10 ) ,
    c(bs,ba) ~ dnorm( 0 , 1 ) ,
    sigma ~ dunif( 0 , 10 )
  ) ,
  data = foxes.sc )
plot(precis(H5.2))

Counterfactual plot area

# prepare new counterfactual data
GS.avg <- mean( foxes.sc$groupsize )
A.seq <- seq( from=-3 , to=3 , length.out=30 )
pred.data <- data.frame(
  area=A.seq,
  groupsize=GS.avg
)

# compute counterfactual mean divorce (mu)
mu <- link( H5.2 , data=pred.data )

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

mu.mean <- apply( mu , 2 , mean )
mu.PI <- apply( mu , 2 , PI )

# simulate counterfactual divorce outcomes
R.sim <- sim( H5.2 , data=pred.data , n=1e4 )

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

R.PI <- apply( R.sim , 2 , PI )

# display predictions, hiding raw data with type="n"
plot( weight ~ area , data=foxes.sc , type="n" )
lines( A.seq , mu.mean )
shade( mu.PI , A.seq )
shade( R.PI , A.seq )

Counterfactual plot groupsize

# prepare new counterfactual data
GS.seq <- seq( from=-3 , to=3 , length.out=30 )
A.avg <- mean( foxes.sc$area )
pred.data <- data.frame(
  area=A.avg,
  groupsize=GS.seq
)

# compute counterfactual mean divorce (mu)
mu <- link( H5.2 , data=pred.data )

## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]

mu.mean <- apply( mu , 2 , mean )
mu.PI <- apply( mu , 2 , PI )

# simulate counterfactual divorce outcomes
R.sim <- sim( H5.2 , data=pred.data , n=1e4 )

## [ 1000 / 10000 ]
[ 2000 / 10000 ]
[ 3000 / 10000 ]
[ 4000 / 10000 ]
[ 5000 / 10000 ]
[ 6000 / 10000 ]
[ 7000 / 10000 ]
[ 8000 / 10000 ]
[ 9000 / 10000 ]
[ 10000 / 10000 ]

R.PI <- apply( R.sim , 2 , PI )

# display predictions, hiding raw data with type="n"
plot( weight ~ groupsize , data=foxes.sc , type="n" )
lines( A.seq , mu.mean )
shade( mu.PI , A.seq )
shade( R.PI , A.seq )

5H3

H5.3a <- map(
  alist(
    weight ~ dnorm( mu , sigma ) ,
    mu <- a + bs*groupsize + bf*avgfood ,
    a ~ dnorm( 0 , 10 ) ,
    c(bs,bf) ~ dnorm( 0 , 1 ) ,
    sigma ~ dunif( 0 , 10 )
  ) ,
  data = foxes.sc )
plot(precis(H5.3a))

H5.3b <- map(
  alist(
    weight ~ dnorm( mu , sigma ) ,
    mu <- a + bs*groupsize + bf*avgfood + ba*area ,
    a ~ dnorm( 0 , 10 ) ,
    c(bs,bf,ba) ~ dnorm( 0 , 1 ) ,
    sigma ~ dunif( 0 , 10 )
  ) ,
  data = foxes.sc )
plot(precis(H5.3b))

pairs( ~ groupsize + avgfood + area + weight, data=foxes.sc , col=rangi2 )

H5.3c <- map(
  alist(
    weight ~ dnorm( mu , sigma ) ,
    mu <- a + bs*groupsize + ba*area ,
    a ~ dnorm( 0 , 10 ) ,
    c(bs,ba) ~ dnorm( 0 , 1 ) ,
    sigma ~ dunif( 0 , 10 )
  ) ,
  data = foxes.sc )
plot(precis(H5.3c))