Chapter4

Medium

4M1

trials <- 1e4
mu.prior.samples <- rnorm(n = trials, mean = 0, sd = 10)
sigma.prior.samples <- runif(n = trials, min = 0, max = 10)
simulated.heights.from.prior <- rnorm(n = trials, mean = mu.prior.samples, sd = sigma.prior.samples)

dens(simulated.heights.from.prior)

4M2

ff <- alist(
  y = dnorm(mu, sigma),
  mu = dnorm(0,10),
  sigma = dunif(0,10)
)

4M4

ff <- alist(
  heigth = dnorm(mu, sigma),
  mu <- a + exp(log_b)*year,
  
  a = dnorm(175,20),
  log_b = dnorm(0, 20), # use log so that b is always > 0, as students are unlikely to shrink.
  sigma = dunif(0,20) # year is unlikely to account for much of the variation, so expect a large SD
  
  
)

4M5

ff <- alist(
  heigth = dnorm(mu, sigma),
  mu <- a + exp(log_b)*year,
  
  a = dnorm(120,25),
  log_b = dnorm(0, 20), # use log so that b is always > 0, as students are unlikely to shrink.
  sigma = dunif(0,sqrt(64) ) # year is unlikely to account for much of the variation, so expect a large SD
)

4H1

# load data
data(Howell1)
d <- Howell1
d$weight.centered <- scale(d$weight)

# build model
model <- map(
  alist(
    height ~ dnorm(mean = mu, sd = sigma),
    mu <- alpha + beta*weight.centered,
    alpha ~ dnorm(mean = 0, sd = 100),
    beta ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 64)
  ),
  data = d
)

individual.weights <- c(46.95, 43.72, 64.78, 32.59, 54.63)
individual.weights.centered <- (individual.weights - mean(d$weight)) / sd(d$weight)
simulated.heights <- sim(model, data = list(weight.centered = individual.weights.centered))

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

sim.mean<-apply(simulated.heights, 2, mean)
sim.89pi<- apply(simulated.heights,2, HPDI,0.89)

4H3

d3 <- subset(d, age<18)

d3$weight.centered <- scale(d3$weight)

# build model
model <- map(
  alist(
    height ~ dnorm(mean = mu, sd = sigma),
    mu <- alpha + beta*weight.centered,
    alpha ~ dnorm(mean = 100, sd = 30),
    beta ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 64)
  ),
  data = d3
)

precis(model)

##         Mean StdDev   5.5%  94.5%
## alpha 108.32   0.61 107.34 109.29
## beta   24.23   0.61  23.25  25.20
## sigma   8.44   0.43   7.75   9.13

plot(y=d3$height,x=d3$weight.centered, col=ifelse(d3$age>14,"red","black"))

cf=coef(model)
abline(a=cf["alpha"], b=cf["beta"])

heights <- sim(model, data=list(weight.centered=seq(-3,3,0.01)))

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

heights.mean <- apply(heights,2,mean)
lines(heights.mean~seq(-3,3,0.01))
heights.89PI <- apply(heights,2,HPDI)
shade(heights.89PI, seq(-3,3,0.01),col=col.alpha("blue",0.3))

4H4

d$weight.log <- log(d$weight)

# build model
model <- map(
  alist(
    height ~ dnorm(mean = mu, sd = sigma),
    mu <- alpha + beta*weight.log,
    alpha ~ dnorm(mean = 150, sd = 30),
    beta ~ dnorm(mean = 0, sd = 10),
    sigma ~ dunif(min = 0, max = 64)
  ),
  data = d
)

precis(model)

##         Mean StdDev   5.5%  94.5%
## alpha -23.24   1.33 -25.37 -21.11
## beta   46.92   0.38  46.31  47.53
## sigma   5.14   0.16   4.89   5.38

plot(y=d$height,x=d$weight.log)

cf=coef(model)
abline(a=cf["alpha"], b=cf["beta"])

sim.weights = seq(1,5,0.01)

sim.mu <- link(model, data=list(weight.log=sim.weights))

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

sim.mu.89PI <- apply(sim.mu,2,HPDI)
shade(sim.mu.89PI, sim.weights,col=col.alpha("red",0.3))

heights <- sim(model, data=list(weight.log=sim.weights))

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

heights.mean <- apply(heights,2,mean)
lines(heights.mean~sim.weights)
heights.89PI <- apply(heights,2,HPDI)
shade(heights.89PI, sim.weights,col=col.alpha("blue",0.3))