The Rats growth data — R cross-check¶

bayesm (rhierLinearModel, rhierLinearMixture) + lme4¶

Companion to rats_python.ipynb (from-scratch) and rats_pymc.ipynb (PyMC). Same hierarchical-linear growth model, fitted with the standard R tools:

  • bayesm::rhierLinearModel — Bayesian single-normal hierarchy (= BiRats / our K=1).
  • bayesm::rhierLinearMixture — the mixture extension (K≥2): are there growth types?
  • lme4::lmer — frequentist REML benchmark.

Note: bayesm's default Inverse-Wishart prior on $V_\beta$ ($5\cdot I$) is far too tight for this scale (weights ~240 g; rat intercepts vary by ~14 g), so we pass a scaled V — otherwise the rat variation is wrongly pushed into the residual.

The fitted model¶

Data. 30 rats, each weighed at five ages $x_j \in \{8, 15, 22, 29, 36\}$ days. $y_{ij}$ = weight (grams) of rat $i$ at age $x_j$.

Likelihood — one linear growth curve per rat (age centred at the midpoint $\bar x = 22$): $$y_{ij} = \alpha_i + \beta_i\,(x_j - 22) + \varepsilon_{ij}, \qquad \varepsilon_{ij}\sim N(0,\sigma^2).$$ Centring makes the parameters interpretable: $\alpha_i$ is rat $i$'s weight at day 22 (the middle of the study) and $\beta_i$ is its growth rate in g/day. The response is raw weight in grams (not log-weight) — the canonical BUGS Rats/BiRats specification, which reproduces the textbook population values (≈ 242 g at day 22, ≈ 6.2 g/day).

Hierarchy — the population of growth curves. The per-rat (intercept, slope) pairs share a population distribution: $$(\alpha_i,\beta_i)' \sim N(\mu,\ \Sigma)\quad\text{(BiRats: a single bivariate normal)},$$ generalised below to a $K$-component mixture of normals $(\alpha_i,\beta_i)'\sim\sum_{k=1}^{K}\pi_k\,N(\mu_k,\Sigma_k)$ via rhierLinearMixture, to test for distinct growth types (segments).

Priors (weakly informative, conjugate — in bayesm's parameterisation): $A_\mu = 0.01$ on the population mean; $\Sigma\sim \mathrm{IW}\big(\nu=4,\ V=\mathrm{diag}(200,\ 0.3)\big)$ — the scale is deliberately matched to the data, since bayesm's default $5\cdot I$ is far too tight for weights of ~240 g (see the note above, otherwise the rat-to-rat variation is wrongly pushed into the residual); per-rat error variances $\sigma^2_i\sim \mathrm{IG}(\nu_e/2,\ \cdot)$ with $\nu_e=3$; and the mixture weights $\pi$ take rhierLinearMixture's default Dirichlet prior.

In [1]:
.libPaths(c("C:/Users/user/R/win-library/4.6", .libPaths()))
suppressMessages({library(bayesm); library(lme4)})

Y <- as.matrix(read.csv("rats.csv", row.names=1)); agec <- c(8,15,22,29,36) - 22
rd <- lapply(1:30, function(i) list(y = as.numeric(Y[i,]), X = cbind(1, agec)))
V  <- diag(c(200, 0.3))                                   # Vbeta prior scaled to the data

# --- bayesm single-normal hierarchy (BiRats) ---
o1 <- rhierLinearModel(Data = list(regdata = rd, Z = matrix(1, 30, 1)),
                       Prior = list(nu = 4, V = V, nu.e = 3, ssq = rep(36, 30)),
                       Mcmc = list(R = 6000, keep = 2, nprint = 0))
n <- nrow(o1$Deltadraw); ki <- (n%/%2 + 1):n
D  <- colMeans(o1$Deltadraw[ki, , drop = FALSE]); Vb <- matrix(colMeans(o1$Vbetadraw[ki, , drop = FALSE]), 2, 2)
rlm <- c(D[1], D[2], sqrt(Vb[1,1]), sqrt(Vb[2,2]), Vb[1,2]/sqrt(Vb[1,1]*Vb[2,2]), sqrt(mean(o1$taudraw[ki,])))

# --- lme4 frequentist benchmark ---
long <- data.frame(weight = as.vector(t(Y)), agec = rep(agec, 30), rat = factor(rep(1:30, each = 5)))
fm <- lmer(weight ~ agec + (agec | rat), long, REML = TRUE)
fe <- fixef(fm); vc <- as.data.frame(VarCorr(fm))
lme <- c(fe[1], fe[2], vc$sdcor[1], vc$sdcor[2], vc$sdcor[3], vc$sdcor[4])

tab <- rbind(`bayesm rhierLinearModel` = rlm, `lme4 lmer` = lme)
colnames(tab) <- c("intercept", "slope", "alpha_sd", "beta_sd", "corr", "resid_sigma")
print(round(tab, 2))

 
Starting Gibbs Sampler for Linear Hierarchical Model
    30  Regressions
    1  Variables in Z (if 1, then only intercept)
 
Prior Parms: 
Deltabar
     [,1] [,2]
[1,]    0    0
A
     [,1]
[1,] 0.01
nu.e (d.f. parm for regression error variances)=  3
Vbeta ~ IW(nu,V)
nu =  4
V 
     [,1] [,2]
[1,]  200  0.0
[2,]    0  0.3
 
MCMC parms: 
R=  6000  keep=  2  nprint=  0
 

                        intercept slope alpha_sd beta_sd corr resid_sigma
bayesm rhierLinearModel    242.29  6.18    14.49    0.47 0.69        7.38
lme4 lmer                  242.65  6.19    14.25    0.51 0.67        6.01

Cross-check — all four engines agree¶

engine intercept slope α SD β SD corr σ
from-scratch (Gibbs) 242.3 6.21 14.0 0.51 0.55 7.3
PyMC (NUTS) 242.6 6.18 14.5 0.52 0.59 6.1
bayesm rhierLinearModel 242.2 6.19 14.4 0.47 0.68 7.4
lme4 lmer 242.7 6.19 14.3 0.51 0.67 6.0

Population intercept ≈ 242 g (weight at day 22), slope ≈ 6.2 g/day, α-SD ≈ 14 g, β-SD ≈ 0.5 g/day, positive intercept–slope correlation — identical across the from-scratch Gibbs, PyMC, bayesm, and lme4. (The Bayesian intercept–slope correlation runs a touch higher, ~0.68 vs lme4's REML 0.67 and the others ~0.55–0.59 — well within posterior uncertainty.)

In [2]:
bayes_b <- apply(o1$betadraw[, , ki], c(1,2), mean)                      # (30,2) shrunken per-rat
ols <- t(sapply(rd, function(d) coef(lm(d$y ~ 0 + d$X))))
plot(ols[,2], bayes_b[,2], pch = 19, col = "firebrick",
     xlab = "per-rat OLS slope", ylab = "bayesm (shrunken) slope",
     main = "Shrinkage: per-rat slopes pulled toward the mean")
abline(0, 1, lty = 3); abline(h = D[2], col = "gray", lty = 2)
legend("topleft", c("y = x", sprintf("pop. slope %.2f", D[2])), lty = c(3,2), col = c("black","gray"), bty = "n")
No description has been provided for this image
In [3]:
# --- mixture: are there growth types? (mubar anchored at the pooled mean) ---
suppressMessages(library(ggplot2))
pooled <- colMeans(t(sapply(rd, function(d) coef(lm(d$y ~ 0 + d$X)))))
cat("rhierLinearMixture -- component means [intercept, slope] and weights:\n")
seg <- data.frame()
for (K in 1:3) {
  om <- rhierLinearMixture(Data = list(regdata = rd),
        Prior = list(ncomp = K, mubar = matrix(pooled, 1, 2), Amu = 0.01, nu = 4, V = V),
        Mcmc = list(R = 6000, keep = 2, nprint = 0))
  pd <- om$nmix$probdraw; cd <- om$nmix$compdraw; ks <- (nrow(pd)%/%2 + 1):nrow(pd)
  w <- colMeans(pd[ks, , drop = FALSE]); mus <- matrix(0, K, 2)
  for (r in ks) { cp <- cd[[r]]; ord <- order(sapply(1:K, function(k) cp[[k]]$mu[1]))
    for (k in 1:K) mus[k,] <- mus[k,] + cp[[ord[k]]]$mu }
  mus <- mus / length(ks)
  cat(sprintf("  ncomp=%d  weights=%s  means=%s\n", K, paste(round(w,2), collapse=","),
              paste(apply(round(mus,1), 1, function(z) sprintf("[%.0f,%.1f]", z[1], z[2])), collapse=" ")))
  seg <- rbind(seg, data.frame(K = K, comp = 1:K, intercept = mus[,1], slope = mus[,2], weight = w))
}
# Plot: component mean intercepts collapse toward the pooled mean as K grows -> no distinct growth types
options(repr.plot.width = 6.6, repr.plot.height = 4.3)
ggplot(seg, aes(factor(K), intercept, size = weight)) +
  geom_point(colour = 'firebrick', alpha = .8) +
  geom_hline(yintercept = pooled[1], linetype = 'dashed', colour = 'grey50') +
  scale_size_continuous(range = c(2, 8)) +
  labs(title = 'Mixture components collapse -> no growth segments',
       subtitle = 'component mean weight-at-day-22 vs K (dashed = pooled mean; bigger = higher weight)',
       x = 'number of components K', y = 'component mean intercept (g @ day 22)', size = 'weight') +
  theme_minimal(base_size = 12)

rhierLinearMixture -- component means [intercept, slope] and weights:

Z not specified

 
Starting MCMC Inference for Hierarchical Linear Model:
   Normal Mixture with 1 components for first stage prior
   for  30  cross-sectional units
 
Prior Parms: 
nu.e = 3
nu = 4
V 
     [,1] [,2]
[1,]  200  0.0
[2,]    0  0.3
mubar 
         [,1]     [,2]
[1,] 242.6533 6.185714
Amu 
     [,1]
[1,] 0.01
a 
[1] 5
 
MCMC Parms: 
R=  6000  keep=  2  nprint=  0

  ncomp=1  weights=1  means=[241,6.0]
Z not specified

 
Starting MCMC Inference for Hierarchical Linear Model:
   Normal Mixture with 2 components for first stage prior
   for  30  cross-sectional units
 
Prior Parms: 
nu.e = 3
nu = 4
V 
     [,1] [,2]
[1,]  200  0.0
[2,]    0  0.3
mubar 
         [,1]     [,2]
[1,] 242.6533 6.185714
Amu 
     [,1]
[1,] 0.01
a 
[1] 5 5
 
MCMC Parms: 
R=  6000  keep=  2  nprint=  0

  ncomp=2  weights=0.51,0.49  means=[235,6.1] [246,6.0]
Z not specified

 
Starting MCMC Inference for Hierarchical Linear Model:
   Normal Mixture with 3 components for first stage prior
   for  30  cross-sectional units
 
Prior Parms: 
nu.e = 3
nu = 4
V 
     [,1] [,2]
[1,]  200  0.0
[2,]    0  0.3
mubar 
         [,1]     [,2]
[1,] 242.6533 6.185714
Amu 
     [,1]
[1,] 0.01
a 
[1] 5 5 5
 
MCMC Parms: 
R=  6000  keep=  2  nprint=  0

  ncomp=3  weights=0.36,0.3,0.34  means=[232,6.2] [242,6.2] [254,6.2]
No description has been provided for this image

Mixture — no growth segments¶

ncomp weights component means [intercept, slope]
1 1.0 [241, 6.0]
2 0.49, 0.51 [236, 6.1], [247, 6.3]
3 ~⅓ each three near-identical

The mixture components are nearly identical (all ≈ 242 g / 6.2 g/day) — rhierLinearMixture finds no discrete growth types, matching the from-scratch WAIC verdict (K=1 best). With 30 rats the population of growth curves is a single bivariate normal.

Conclusion. Four independent engines — from-scratch Gibbs, PyMC, bayesm, and lme4 — agree on the Rats hierarchical-linear fit (intercept ≈ 242 g, slope ≈ 6.2 g/day, partial-pooling shrinkage), and both Bayesian mixture routes confirm a single normal suffices. Rats is the clean canonical backbone of the hierarchical-linear (and h-linear-mixture) model.