如何在nlme非线性混合效应模型中计算杠杆值与Cook距离？

问题背景 #

使用nlme包构建非线性混合效应模型后，计算杠杆值和Cook距离时出现以下报错：

杠杆值计算报错 #

执行lev<-hat(model.matrix(model))时触发：

Error in x$terms %||% attr(x, "terms") %||% stop("no terms component nor attribute") : no terms component nor attribute

Cook距离计算报错 #

执行cd<-cooks.distance(model)时触发：

Error in cooks.distance.lme(model) : not implemented for "nlme" objects"

模型代码如下：

formula = log_Lobs ~ log(150*((1 + ((150/Lt_1)^(1/exp(p))-1)*exp(-exp(k)*td/365))^(-exp(p))))
model <- do.call(nlme,
                    list(formula,
                         fixed = c(p ~ 1, k ~ 1 + season2),
                         random = k ~ 1 | id,
                         data = data_select,
                         start = list(fixed = c(p, k)),
                         na.action = na.exclude,
                         control=list(maxIter=1e6, msMaxIter = 1e6, msVerbose = TRUE)
                    ))

解决方案 #

一、计算杠杆值 #

nlme模型的model.matrix()返回对象不包含terms属性，无法直接用hat()处理。可通过两种方式计算：

方法1：基于固定效应设计矩阵

提取固定效应的设计矩阵，结合协方差矩阵构建投影矩阵计算杠杆值：

# 提取固定效应设计矩阵
X <- model.matrix(model, type = "fixed")
# 计算杠杆值
lev <- hat(t(X) %*% solve(vcov(model)) %*% X) / nrow(X)

方法2：基于模型预测矩阵

利用nlme模型内置的预测矩阵（predMatrix）计算，更适配非线性混合效应模型场景：

# 获取模型预测矩阵
P <- model$predMatrix
# 计算杠杆值
lev <- diag(P %*% t(P))

二、计算Cook距离 #

nlme未内置Cook距离计算函数，可通过两种方式实现：

方法1：基于残差与杠杆值的公式计算

Cook距离核心公式为：$D_i = \frac{r_i^2}{p} \times \frac{h_{ii}}{1 - h_{ii}}$，其中$r_i$为标准化残差，$h_{ii}$为杠杆值，$p$为固定效应参数个数。代码实现：

# 获取皮尔逊标准化残差
res <- residuals(model, type = "pearson")
# 计算杠杆值（用方法1或2均可）
X <- model.matrix(model, type = "fixed")
lev <- hat(t(X) %*% solve(vcov(model)) %*% X) / nrow(X)
# 获取固定效应参数数量
p <- length(fixef(model))
# 计算Cook距离
cd <- (res^2 / p) * (lev / (1 - lev))

方法2：逐观测删除重拟合计算

通过循环删除每个观测后重新拟合模型，计算参数估计的变化程度（适合小数据集，计算量较大）：

# 初始化Cook距离向量
cd <- numeric(nrow(data_select))
# 循环处理每个观测
for (i in 1:nrow(data_select)) {
  data_temp <- data_select[-i, ]
  # 尝试拟合删除当前观测后的模型，用原模型参数作为初始值加速拟合
  model_temp <- try(do.call(nlme,
                    list(formula,
                         fixed = c(p ~ 1, k ~ 1 + season2),
                         random = k ~ 1 | id,
                         data = data_temp,
                         start = list(fixed = fixef(model)),
                         na.action = na.exclude,
                         control=list(maxIter=1e6, msMaxIter = 1e6, msVerbose = FALSE)
                    )), silent = TRUE)
  # 若拟合成功则计算Cook距离，失败标记为NA
  if (!inherits(model_temp, "try-error")) {
    cd[i] <- t(fixef(model) - fixef(model_temp)) %*% solve(vcov(model)) %*% (fixef(model) - fixef(model_temp)) / p
  } else {
    cd[i] <- NA
  }
}

内容的提问来源于stack exchange，提问作者TKH_9