多元线性回归列表上的 Wald 检验

Posted 2023-03-12

【中文标题】多元线性回归列表上的 Wald 检验

【英文标题】：Wald test on a list of multiple linear regressions

【发布时间】：2020-04-01 01:08:04

【问题描述】：

使用此处显示的使用 lm 的 69 个模型的新创建列表：Looping through many multiple regressions

我正在尝试运行 Wald 测试，但它似乎无法同时在 69 个模型上运行。只有当我指定对列表中的一个模型进行 Wald 测试时，它才有效。关于如何通过比较系数 1 和 2 获得整体 Wald 检验的任何想法？

install.packages("aod") 
library(aod) 
wald.test(b = coef(models[[1]]), Sigma = vcov(models[[1]]), Terms = 1:2) #this one works only 
wald.test(b = coef(models), Sigma = vcov(models,Terms = 1:2)) #this one does not work

【参考方案1】：

实际上，您可以使用other question 的解决方案一步完成 Wald 测试。只需将 $ 的 Wald 测试作为元素添加到模型输出中。 （注意我这里用的是reformulate()而不是as.formula()）。

n.ids <- 5
models <- lapply(1:n.ids, function(x) 
  fo <- reformulate(paste0("id", x, c(".PRE", ".POST")), response="growth_rate")
  fit <- lm(fo, output)
  fit$wald <- aod::wald.test(b=coef(fit), Sigma=vcov(fit), Terms=1:2)$result$chi2
  return(fit)
)

之后，您可以方便地提取模型拟合（此处为[[1]]^st 模型），

models[[1]]
# Call:
#   lm(formula = fo, data = output)
# 
# Coefficients:
#   (Intercept)      id1.PRE     id1.POST  
#     0.0006867    0.3621801   -0.6805282

采取summary()统计，

summary(models[[1]])
# Call:
#   lm(formula = fo, data = output)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -1.6332 -0.7522 -0.2040  0.3244  2.3479 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.0006867  0.4707891   0.001    0.999
# id1.PRE      0.3621801  0.7245488   0.500    0.632
# id1.POST    -0.6805282  0.4666070  -1.458    0.188
# 
# Residual standard error: 1.449 on 7 degrees of freedom
# Multiple R-squared:  0.2783,  Adjusted R-squared:  0.07208 
# F-statistic:  1.35 on 2 and 7 DF,  p-value: 0.3194

最后是 Wald 检验统计量：

models[[1]]$wald
#      chi2        df         P 
# 0.2602116 2.0000000 0.8780026

数据

output <- structure(list(growth_rate = c(-1.69916848732564, 0.580880770531803, 
-2.69289070966401, 1.35526046076076, -0.112984932782402, -0.462646659510964, 
0.719326790475801, -1.00142215628315, 2.47637095683882, 0.43656578200224
), id1.PRE = c(-0.658214571692787, 0.618757883859616, -1.00319408411405, 
-0.182500151494055, 0.900088109789603, 1.0196598975258, 0.412275530237581, 
0.392417622938614, -0.455987647372666, 0.102734423799357), id2.PRE = c(-0.170523375036079, 
0.686700233146119, 0.309985112368748, -1.90034614054823, -1.09227380863711, 
0.808287849081857, -0.499631929432237, 0.166532696960525, 0.302124067552245, 
0.387404901557692), id3.PRE = c(0.599764947203108, 0.89758476990722, 
-0.255477016875339, -0.356745234827032, 0.448024379737768, -0.256047154841515, 
0.585227234029825, -0.166333176678276, 0.823998488345663, -1.443644743169
), id4.PRE = c(-2.80604357500051, -1.09530031299481, 0.175000065646341, 
-0.392833839818346, -0.409796662979807, -0.195621477818208, -0.00577716957657597, 
-0.482953411765706, -2.02722251486943, 0.300223614906696), id5.PRE = c(-0.0622634352162204, 
0.040326767924506, 0.0497415705870156, 0.180785254746358, 1.33406548126792, 
0.0836709755262436, -0.73864218643137, -1.50512979963419, -0.515794856492169, 
-0.00132857851903661), id1.POST = c(1.31893320782006, -0.865022433359845, 
1.02422314119992, -1.5293010769327, 0.388204583082128, 0.026043363711214, 
1.91192868077842, -0.281576076956257, -0.430489296095949, -0.35410853516199
), id2.POST = c(-0.311676684306798, -0.368607205934452, -1.48144614489904, 
0.555325334696273, 0.570876384813286, -1.31976394074723, -0.612004459486653, 
1.22377107478453, 1.59353205696141, -0.690394114816965), id3.POST = c(-0.720115908724216, 
-0.628303814726474, 0.861931079614817, 0.462347223839036, -1.03796670332245, 
-1.1879031168125, 0.26141696844815, -1.42623852487712, 0.538672447311226, 
-2.1422820318952), id4.POST = c(-1.43103098130001, 0.70344192904054, 
-0.599438467172651, -2.6958955065305, -0.140567798927688, 0.393778135755608, 
0.931307986770656, 0.107100043420371, 0.17241808370283, -0.334264159331698
), id5.POST = c(-0.468527805538746, -0.231147053156926, 1.16752987470932, 
-0.025411116476825, 0.377661302479598, 1.60646402608108, -0.142579077289651, 
0.269743622121994, -0.562298443852192, 1.41564052399418)), row.names = c(NA, 
-10L), class = "data.frame")

【讨论】：

这太好了，非常感谢！但是你认为有可能从比较 PRE_id 系数和 POST_id 系数的 Wald 检验中获得总体 p 值吗？我正在尝试获得一个整体的单一 p 值，这可能吗？