4H1

I assume that the object here is to use the linear regression fit to the data from adults for the predictions. Copying (and slightly editing) McElreath’s code for fitting the model (centered version – it would be a little more straightforward to use the uncentered version):

> set.seed(19827) # for reproducibility
> 
> library(rethinking)
Loading required package: rstan
Loading required package: ggplot2
Loading required package: StanHeaders
rstan (Version 2.12.1, packaged: 2016-09-11 13:07:50 UTC, GitRev: 85f7a56811da)
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.59)
> data(Howell1)
> d <- Howell1
> d2 <- d[d$age >= 18, ]
> 
> d2$weight.c <- d2$weight - mean(d2$weight)
> 
> # fit model
> m4.4 <- map(
+     alist(
+         height ~ dnorm(mu, sigma),
+         mu <- a + b*weight.c ,
+         a ~ dnorm(178 ,100),
+         b ~ dnorm(0, 10),
+         sigma ~ dunif(0, 50)
+     ),
+     data=d2)
> 
> # summarize fit
> precis(m4.4, corr=TRUE)
        Mean StdDev   5.5%  94.5% a b sigma
a     154.60   0.27 154.17 155.03 1 0     0
b       0.91   0.04   0.84   0.97 0 1     0
sigma   5.07   0.19   4.77   5.38 0 0     1

To generate predictions using this model, it’s necessary to subtract the mean from each supplied weight:

> weight <- c(46.95, 43.72, 64.78, 32.59, 54.63)
> (weight.c <- weight - mean(d2$weight))
[1]   1.959514  -1.270486  19.789514 -12.400486   9.639514

Then, getting both PI and HPDI intervals for the predictions of height:

> sim.height <- sim(m4.4, data=list(weight.c=weight.c))
[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
> colMeans(sim.height)
[1] 156.4575 153.2647 172.6661 143.5328 163.3402
> (height.PI <- apply(sim.height, 2, PI, prob=0.89))
        [,1]     [,2]     [,3]     [,4]     [,5]
5%  148.2225 145.1733 164.4094 135.6866 155.5564
94% 165.0105 161.3084 180.6698 152.1499 171.4783
> (height.HPDI <- apply(sim.height, 2, HPDI, prob=0.89))
          [,1]     [,2]     [,3]     [,4]     [,5]
|0.89 147.5293 145.3687 164.4165 134.9333 155.6011
0.89| 164.2229 161.4611 180.7248 151.1782 171.5142

4H2

Selecting those below 18:

> d3 <- d[d$age < 18, ]
> d3$weight.c <- with(d3, weight - mean(weight))
> dim(d3)
[1] 192   5
> head(d3)
   height   weight  age male  weight.c
19 121.92 19.61785 12.0    1  1.203661
20 105.41 13.94795  8.0    0 -4.466239
21  86.36 10.48931  6.5    0 -7.924878
24 129.54 23.58678 13.0    1  5.172591
25 109.22 15.98912  7.0    0 -2.425075
29 137.16 27.32892 17.0    1  8.914725
> zapsmall(colMeans(d3))
   height    weight       age      male  weight.c 
108.31885  18.41419   7.72187   0.47917   0.00000

[a]

Specifying and fitting a linear regression to the kids’ data:

> mod.kids <- map(
+     alist(
+         height ~ dnorm(mu, sigma),
+         mu <- a + b*weight.c ,
+         a ~ dnorm(108, 100),
+         b ~ dnorm(0, 10),
+         sigma ~ dunif(0, 50)
+     ),
+     data=d3)
> 
> # summarize fit
> precis(mod.kids)
        Mean StdDev   5.5%  94.5%
a     108.32   0.61 107.35 109.29
b       2.72   0.07   2.61   2.83
sigma   8.44   0.43   7.75   9.13

Interpretation:

\(a = 108.32\): average height at average weight (essentially average marginal height).

\(b = 2.72\): each increment of 1 kg in weight is associated on average with an additional 2.72 cm of height. Thus a 10 kg increment in weight would be associated on average with a 27.2 cm increase in height.

\(\sigma = 8.44\): points are typically about 8.5 cm from the regression line.

[b]

Graphing the (centered) data, showing the fitted line:

> plot(height ~ weight.c, data=d3, ylim=c(50, 200))  # plot data
> 
> weight.c.seq <- with(d3, seq(min(weight.c), max(weight.c), length=100))
> 
> # pointwise HPDI for the (conditional) mean height
> mu <- link(mod.kids, data=data.frame(weight.c=weight.c.seq))
[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
> mu.HPDI <- apply(mu, 2, HPDI, prob=0.89)
> shade(mu.HPDI, weight.c.seq, col="magenta")
> 
> # pointwise PI for height
> sim.height <- sim(mod.kids, data=list(weight.c=weight.c.seq))
[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
> height.PI <- apply(sim.height, 2, PI, prob=0.89)
> shade(height.PI, weight.c.seq)
> 
> abline(coef(mod.kids)[c("a", "b")]) # plot fitted line

[c]

The relationship is obviously nonlinear, though monotone. Since the bulge points to the left, I’d try a transformation of weight down the ladder of powers; I’d start with the log transformation. Alternatively, thinking about the dimensions of the variables, height is a linear measure of size but weight is probably related to volume That suggest the cube-root transformation of weight (1/3 power), probably not too different from log (“0” power). Finally, I’m aware that BMI is defined as weight divided by height\(^2\), which suggests the square-root transformation of weight, though I’m not sure what the rationale is for the definition.

What do Profs. Box and Tidwell say?

> library(car)

Attaching package: 'car'
The following object is masked from 'package:rethinking':

    logit
> boxTidwell(height ~ weight, data=d)
 Score Statistic p-value MLE of lambda
       -35.38973       0    -0.1128969

iterations =  4

So, logs look pretty good (\(-0.11 \approx 0\)).

4H3

[a]

Using the log transformation as suggested (but logs base-2):

> mod.log <- map(
+     alist(
+         height ~ dnorm(mu, sigma),
+         mu <- a + b*log2(weight) ,
+         a ~ dnorm(178, 100),
+         b ~ dnorm(0, 140), # compensating for log2 (but where did 100 come from?)
+         sigma ~ dunif(0, 50)
+     ),
+     data=d)
> 
> # summarize fit
> precis(mod.log)
        Mean StdDev   5.5%  94.5%
a     -23.79   1.34 -25.92 -21.65
b      32.63   0.27  32.21  33.05
sigma   5.13   0.16   4.89   5.38

Interpretation:

\(a = -23.79\) isn’t directly interpretable since 0 on the log\(_2\) scale correspond to a weight of 1 kg (and, of course, height can’t be negative).

\(b = 32.63\): Doubling weight is associated on average with an increase of about 33 cm of height.

\(\sigma = 5.13\): On average the points are about 5 cm from the fitted regression.

[b]

Here is the plot with the fitted line, the 97% HPDI for the conditional mean height, and the 97% HPDI for predicted height:

> plot(height ~ weight, data=d, col=col.alpha(rangi2, 0.4), ylim=c(40, 180))
> 
> weight.seq <- with(d, seq(min(weight), max(weight), length=100))
> 
> # pointwise HPDI for the (conditional) mean height
> mu <- link(mod.log, data=data.frame(weight=weight.seq))
[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
> mu.HPDI <- apply(mu, 2, HPDI, prob=0.97)
> shade(mu.HPDI, weight.seq, col="magenta")
> 
> # pointwise HPDI for height
> sim.height <- sim(mod.log, data=list(weight=weight.seq))
[ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
> height.PI <- apply(sim.height, 2, HPDI, prob=0.97)
> shade(height.PI, weight.seq)
> 
> # plot fitted line
> mu.mean <- colMeans(mu)
> lines(weight.seq, mu.mean)

That seems like a reasonable fit to the data.