机器算法验证 - nlme gls计算的beta系数的SEs计算算法终于固定了吗？ - 吾爱随笔录

我可以在这里阅读： https ://www4.stat.ncsu.edu/~davidian/st732/examples/dental_pa.R 和这里： https ://math.unm.edu/~luyan/stat57918/week14.pdf

那：

警告：gls() 中有一个错误，它不能正确计算 betahat 的基于模型的采样协方差矩阵！因此， gls() 报告的基于模型的标准错误与正确值有一点“偏差”。此外， gls() 不计算稳健的三明治协方差矩阵，因此不提供获得稳健（或经验）标准误差的选项。下面的函数 robust.cov 将计算正确的基于模型的标准误差以及稳健的三明治标准误差。

有谁知道，如果这是固定的？这可能吗，像 nlme 一样古老并被数千名用户使用的软件包返回的值仍然不正确并误导用户？或者这是固定的？

在您看来，应该是链接材料中提出的功能

#  Function to compute correct and robust standard errors from a gls()
#  model object

library(Matrix)
library(magic)

#  u <- model object from fit, m <- total # individuals

robust.cov <- function(u,m){
  form <-  formula(u)
  mf <- model.frame(form,getData(u))
  Xmat <- model.matrix(form,mf)
  Vlist <-  as.list(1:m)
  for (i in 1:m){
    Vlist[[i]] <- getVarCov(u,individual=i,type="marginal")
  }
  V <- Reduce(adiag,Vlist)
  Vinv <- solve(V)
  Sig.model <- solve(t(Xmat)%*%Vinv%*%Xmat)
  resid <- diag(residuals(u,type="response"))
  ones.list <- lapply(Vlist,FUN=function(u){matrix(1,nrow(u),ncol(u))})
  Ones <- Reduce(adiag,ones.list)
  meat <- t(Xmat)%*%Vinv%*%resid%*%Ones%*%resid%*%Vinv%*%Xmat
  Sig.robust <- Sig.model%*%meat%*%Sig.model
  se.robust <- sqrt(diag(Sig.robust))
  se.model <- sqrt(diag(Sig.model))
  return(list(Sig.model=Sig.model,se.model=se.model,Sig.robust=Sig.robust,se.robust=se.robust))
}

或此包中列出的方法：https ://cran.r-project.org/web/packages/clubSandwich/clubSandwich.pdf在处理 GLS 时仍然使用？

编辑：我刚刚检查过。似乎 gls 仍在计算不符合 SAS 的 SE。但与 lme、lme4 和 glmmTMB 是一致的。

使用具有鲁棒估计类型“CR0”的 clubSandwich R 包给出的结果与使用上面介绍的 robust.cov() 函数获得的结果相同。

在我的情况下，lme4 和 gls() 的 Kenward-Roger SE 仅略有不同，并且与稳健估计的结果相去甚远。

                 GLS_original GLS_corrected GLS_KR_Lava2 lmer_original lmer_Sandwich
(Intercept)         0.8215916     0.7537342    0.8441051     0.8215916     0.7537342
Program2            1.0905540     1.0028205    1.1204378     1.0905540     1.0028205
Program3            1.1022808     0.9526228    1.1324858     1.1022808     0.9526228
Timef4              0.3867933     0.2142392    0.3973924     0.3867933     0.2142392
Timef6              0.3867933     0.3408973    0.3973924     0.3867933     0.3408973
Timef8              0.3867933     0.4023248    0.3973924     0.3867933     0.4023248
Timef10             0.3867933     0.5517639    0.3973924     0.3867933     0.5517639
Timef12             0.3867933     0.6185211    0.3973924     0.3867933     0.6185211
Timef14             0.3867933     0.6371859    0.3973924     0.3867933     0.6371859
Program2:Timef4     0.5134169     0.2743518    0.5274858     0.5134169     0.2743518
Program3:Timef4     0.5189377     0.2896523    0.5331578     0.5189377     0.2896523
Program2:Timef6     0.5134169     0.4411397    0.5274858     0.5134169     0.4411397
Program3:Timef6     0.5189377     0.4069225    0.5331578     0.5189377     0.4069225
Program2:Timef8     0.5134169     0.5008398    0.5274858     0.5134169     0.5008398
Program3:Timef8     0.5189377     0.5407081    0.5331578     0.5189377     0.5407081
Program2:Timef10    0.5134169     0.6413605    0.5274858     0.5134169     0.6413605
Program3:Timef10    0.5189377     0.6832411    0.5331578     0.5189377     0.6832411
Program2:Timef12    0.5134169     0.7148875    0.5274858     0.5134169     0.7148875
Program3:Timef12    0.5189377     0.7409071    0.5331578     0.5189377     0.7409071
Program2:Timef14    0.5134169     0.7529247    0.5274858     0.5134169     0.7529247
Program3:Timef14    0.5189377     0.7729042    0.5331578     0.5189377     0.7729042

正如我们所看到的，Sandwich 与 KR 完全相反（关于未校正的）

现在，我想知道该选择哪一个：原始 SE、坚固的三明治 SE 或 Kenward-Roger SE...

编辑：好的，我想我得到了答案。这是模型和经验 SE 之间的区别。默认情况下，R 中的所有线性模型都返回基于模型的模型，以及 GEE（取决于库）——三明治（经验）SE。

这是与软件无关的话题。人们分析 SAS 和 Stata 内部/之间的类似问题。 http://www.biostat.umn.edu/~wguan/class/PUBH7402/notes/lecture9_GEEvsRE_linear.pdf https://www.stata.com/support/faqs/statistics/xtgee-versus-proc-genmod/

建议是对较少数量的集群使用基于模型的 SE。也可以使用适当校正的三明治。他们有很多，但莫雷尔等人。似乎是最常用的校正。

我仍然不确定要报告哪一个，所以我将两者都报告，并根据更保守的方法进行推断，只要类型 1 错误更严重，这似乎更安全。