EM step (single iteration) for random-effects variance

Formula: for the random-intercept model y_ij = x_ij' b + u_j + e_ij with u_j ~ N(0, sigma2_u) and e_ij ~ N(0, sigma2_e), the E-step conditional moments are closed form,


  uhat_j     = sigma2_u sum_i r_ij / (sigma2_e + n_j sigma2_u)
  var(u_j|y) = sigma2_u sigma2_e / (sigma2_e + n_j sigma2_u)

with r_ij = y_ij - x_ij' b, and the M-step is sigma2_u' = mean_j (uhat_j^2 + var(u_j|y)) and sigma2_e' = (sum_ij (r_ij - uhat_j)^2 + sum_j n_j var(u_j|y)) / N.

Usage

Emaxr(y, X, cluster, sigma2_u, sigma2_e, beta = NULL)

Arguments

y: Response, length N.
X: N x p fixed-effects design.
cluster: Length-N cluster labels.
sigma2_u: Current between-cluster variance (>= 0).
sigma2_e: Current residual variance (> 0, or >= 0 with beta).
beta: Optional fixed-effect vector held fixed.

Value

List with estimate, sigma2_u, sigma2_e, beta, u_hat, var_u, J, n, method.

Details

b is the GLS estimate at the current variances, formed cluster-wise from the Sherman-Morrison inverse V_j^-1 = (I - sigma2_u J_j / (sigma2_e + n_j sigma2_u)) / sigma2_e, or supplied fixed via beta, in which case sigma2_e = 0 is admissible and the update reduces to its exact limit.

References

Dempster, Laird & Rubin (1977), JRSS-B 39(1):1-38; Laird & Ware (1982), Biometrics 38(4):963-974, doi:10.2307/2529876.

Examples

set.seed(1)
r <- Emaxr(y = rnorm(10), X = rnorm(10), cluster = rnorm(10), sigma2_u = 0.5, sigma2_e = 0.5); TRUE
#> [1] TRUE