MME方程

在求解混合线性模型的时候，往往需要和MME方程打交道，因此这篇博文重点讲述使用REML进行 b 和 u 的估计。
其中X为固定效应的设计矩阵，Z为随机效应的设计矩阵；G为随机效应方差矩阵，R为误差的方差矩阵

MME方程

将MME方差的左边简化为分块矩阵

分块矩阵

随机效应方差估计的更新式子：

随机效应方差估计理论

同理可得随机误差估计的更新式子：

随机误差的估计

R语言版—广义MME方程版本

# 加载必要包
library(MASS)   # 用于ginv()

# ------------------------- 生成模拟数据 -------------------------
set.seed(123)
n <- 100        # 总观测数
q <- 10          # 随机效应水平数（例如10个品种）
p <- 2            # 固定效应数（含截距）

# 设计矩阵
X <- cbind(1, rnorm(n))           # 固定效应设计矩阵
Z <- matrix(0, n, q)
for(i in 1:n){
  Z[i, sample(1:q, 1)] <- 1       # 随机效应分配
}

# 真实方差组分
sigma2_u_true <- 4
sigma2_e_true <- 1

# 随机效应和误差
u_true <- rnorm(q, 0, sqrt(sigma2_u_true))
e_true <- rnorm(n, 0, sqrt(sigma2_e_true))

# 固定效应系数
beta_true <- c(2, -0.5)

# 响应变量
y <- X %*% beta_true + Z %*% u_true + e_true

# ------------------------- 初始化 -------------------------
# 初始方差值
sigma2_u <- 1
sigma2_e <- 1

# 迭代参数
max_iter <- 100
tol <- 1e-6
iter <- 0
converged <- FALSE

# 存储迭代历史
history <- data.frame(iter = 0, sigma2_u = sigma2_u, sigma2_e = sigma2_e)

# ------------------------- EM迭代求解 -------------------------
while(iter < max_iter && !converged){
  
  iter <- iter + 1
  
  # 计算R逆和G逆（对角阵）
  R_inv <- diag(1/sigma2_e, n)
  G_inv <- diag(1/sigma2_u, q)
  
  # 构建MME左侧矩阵, 参考MME方程
  ## X 为固定小于设计矩阵
  ## Z 为随机效应设计矩阵
  ## R_inv 为随机效应方程矩阵的逆矩阵
  XtRinvX <- t(X) %*% R_inv %*% X
  XtRinvZ <- t(X) %*% R_inv %*% Z
  ZtRinvX <- t(Z) %*% R_inv %*% X
  ZtRinvZ <- t(Z) %*% R_inv %*% Z
  
  # 分块矩阵，见分块矩阵公式
  ## G_inv 为随机效应方差矩阵
  LHS <- rbind(
    cbind(XtRinvX, XtRinvZ),
    cbind(ZtRinvX, ZtRinvZ + G_inv)
  )
  
  # MME右侧向量, 参考MME方程
  ## R_inv 为随机效应方程矩阵的逆矩阵
  RHS <- rbind(
    t(X) %*% R_inv %*% y,
    t(Z) %*% R_inv %*% y
  )
  
  # 求解MME（使用广义逆处理可能的奇异）
  ## 求解固定效应系数 b 和随机效应系数 u
  sol <- ginv(LHS) %*% RHS
  ## p为固定效应数量（含截距），p个系数
  beta_hat <- sol[1:p]
  ## q为随机效应数量，q个系数
  u_hat <- sol[(p+1):(p+q)]
  
  # ------------------------- EM更新方差组分 -------------------------
  # 残差
  e_hat <- y - X %*% beta_hat - Z %*% u_hat
  
  # 计算C矩阵（LHS的逆）中对应u的部分
  # 这里用广义逆的对应分块
  LHS_inv <- ginv(LHS)
  C_uu <- LHS_inv[(p+1):(p+q), (p+1):(p+q)]
  
  # 更新sigma2_e，公式参考理论部分
  sigma2_e_new <- (t(e_hat) %*% e_hat + sum(diag(Z %*% C_uu %*% t(Z)))) / n
  
  # 更新sigma2_u，公式参考理论部分
  sigma2_u_new <- (t(u_hat) %*% u_hat + sum(diag(C_uu))) / q
  
  # 记录
  history <- rbind(history, c(iter, sigma2_u_new, sigma2_e_new))
  
  # 检查收敛
  if(abs(sigma2_u_new - sigma2_u) < tol && abs(sigma2_e_new - sigma2_e) < tol){
    converged <- TRUE
  }
  
  # 更新
  sigma2_u <- as.numeric(sigma2_u_new)
  sigma2_e <- as.numeric(sigma2_e_new)
  
  cat(sprintf("Iter %d: sigma2_u = %.4f, sigma2_e = %.4f\n", iter, sigma2_u, sigma2_e))
}

# ------------------------- 输出结果 -------------------------
cat("\n========== 最终结果 ==========\n")
cat("收敛状态:", if(converged) "收敛" else "未收敛", "\n")
cat("估计的固定效应 (beta):\n")
print(beta_hat)
cat("\n估计的方差组分:\n")
cat("sigma2_u =", sigma2_u, "\n")
cat("sigma2_e =", sigma2_e, "\n")
cat("真实值: sigma2_u =", sigma2_u_true, ", sigma2_e =", sigma2_e_true, "\n")