I'm interested in a simple random-intercepts model: $$Y_{ij} = \alpha_0 + \gamma_i + \epsilon_{ij}$$ where $\gamma_i \sim N(0, \sigma_{\gamma}^2)$ independently of $\epsilon_{ij} \sim N(0, \sigma_{\epsilon}^2)$.
In nlme
:
library(nlme)
data(ergoStool)
l1 <- lme(effort ~ 1, data = ergoStool, random = ~ 1|Subject, method="ML")
and gives estimates of (what I assume are) $\sigma_\gamma$ and $\sigma_\epsilon$:
summary(l1)
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 0.9090556 2.020727
In gamlss
, the same model can be fit, with the random intercepts included in a couple of ways:
library(gamlss)
t1 <- gamlss(effort ~ 1 + random(Subject), data = ergoStool)
t2 <- gamlss(effort ~ 1 + re(random = ~ 1|Subject), data = ergoStool)
The gamlss
models give the same fit, and each has three variance 'components', which agree with one another:
c(sigb = getSmo(t1)$sigb, sige = getSmo(t1)$sige, exp(t1$sigma.coefficients))
sigb sige (Intercept)
0.9090593 1.0610962 1.9043752
summary(t2)
Random effects:
Formula: ~1 | Subject
(Intercept) Residual
StdDev: 0.9090594 1.061096
exp(t2$sigma.coefficients)
(Intercept)
1.904375
The first is obviously the same $\hat{\sigma}_\gamma$, but I can't work out how the other two (the residual standard deviation from the random intercept part, and the estimate of the standard deviation of the outcome) relate to $\hat{\sigma}_\epsilon$. Is there a relationship between these estimates?