r - 如何使用 R 对两个回归的斜率进行韦尔奇 t 检验？

Question

我对具有相同自变量的两组进行回归。然后，我想测试两个回归的斜率是否显着不同。

我读到，当两组之间的样本量和方差不相等时，建议执行韦尔奇 t 检验。我找到了 t.test() 函数，但是我没有在斜坡上应用它。

Data <- data.frame(
  gender = sample (c("men", "women"), 2000, replace = TRUE),
  var1 = sample (c("value1", "value2"), 2000, replace = TRUE),
  var2 = sample (c("valueA", "valueB"), 2000, replace = TRUE),
  var3 = sample (c("valueY", "valueZ"), 2000, replace = TRUE),
  y = sample(0:10, 2000, replace = TRUE)
)

我的两个回归:

lm.male <- lm(y ~ var1 + var2 + var3, data = subset(Data, gender == "men"))
summary(lm.male)

lm.women <- lm(y ~ var1 + var2 + var3, data = subset(Data, gender == "women"))
summary(lm.women)

使用 Stata，我将使用 suet 和 test 函数来执行测试。

有人知道如何在 R 中编写斜率的 Welch t 检验吗？

Answer 1

我不太会准确回答你的问题，而是更普遍的问题是，在 R 中，我如何测试响应变量中疑似方差不等的两组之间斜率差异的假设。

概述

有几个选项，我将介绍其中两个。所有好的选择都涉及将两个数据集组合成一个建模策略，并面对一个“完整”模型，其中包括性别和斜率的交互效应，以及“无交互”模型，该模型具有附加的性别效应，但相同其他变量的斜率。

如果我们准备假设两个性别组的方差相同，我们只需使用普通最小二乘法将我们的两个模型拟合到组合数据中，并使用经典的 F 检验:

Data <- data.frame(
  gender = sample (c("men", "women"), 2000, replace = TRUE),
  var1 = sample (c("value1", "value2"), 2000, replace = TRUE),
  var2 = sample (c("valueA", "valueB"), 2000, replace = TRUE),
  var3 = sample (c("valueY", "valueZ"), 2000, replace = TRUE),
  y = sample(0:10, 2000, replace = TRUE)
)


lm_full <- lm(y ~ (var1 + var2 + var3) * gender, data = Data)
lm_nointeraction <- lm(y ~ var1 + var2 + var3 + gender, data = Data)

# If the variance were equal we could just do an F-test:
anova(lm_full, lm_nointeraction)

但是，这个假设是 Not Acceptable ，因此我们需要一个替代方案。我认为this discussion on Cross-Validated很有用。

选项 1 - 加权最小二乘法

我不确定这是否与韦尔奇的 t 检验相同；我怀疑这是它的更高层次的概括。这是解决该问题的一种非常简单的参数方法。基本上，我们只是同时对响应的方差及其平均值进行建模。然后在拟合过程(变成迭代)中，我们对预计具有较高方差(即更多随机性)的点给予较小的权重。 nlme 包中的 gls 函数 - 广义最小二乘法 - 为我们做到了这一点。

# Option 1 - modelling variance, and making weights inversely proportional to it
library(nlme)
gls_full <- gls(y ~ (var1 + var2 + var3) * gender, data = Data, weights = varPower())
gls_nointeraction <- gls(y ~ var1 + var2 + var3 + gender, data = Data, weights = varPower())

# test the two models against eachother (preferred):
AIC(gls_full, gls_nointeraction) # lower value wins

# or test individual interactions:
summary(gls_full)$tTable

选项 2 - 稳健回归，通过 bootstrap 进行比较

第二个选项是使用 M 估计，它的设计目的是对数据中组内的不等方差具有鲁棒性。比较两个模型的稳健回归的良好做法是选择某种验证统计数据，并使用 Bootstrap 来查看哪个模型在该统计数据上平均表现更好。

这有点复杂，但这是一个使用模拟数据的有效示例:

# Option 2 - use robust regression and the bootstrap
library(MASS)
library(boot)
rlm_full <- rlm(y ~ (var1 + var2 + var3) * gender, data = Data)
rlm_nointeraction <- rlm(y ~ var1 + var2 + var3 + gender, data = Data)

# Could just test to see which one fits best (lower value wins)
AIC(rlm_full, rlm_nointeraction)

# or - preferred - use the bootstrap to validate each model and pick the best one.
# First we need a function to give us a performance statistic on how good
# a model is at predicting values compared to actuality.  Let's use root 
#  mean squared error:
RMSE <- function(predicted, actual){
  sqrt(mean((actual - predicted) ^ 2))
}

# This function takes a dataset full_data, "scrambled" by the resampling vector i.
# It fits the model to the resampled/scrambled version of the data, and uses this
# to predict the values of y in the full original unscrambled dataset.  This is
# described as the "simple bootstrap" in Harrell *Regression Modeling Strategies*,
# buiolding on Efron and Tibshirani.
simple_bootstrap <- function(full_data, i){
  sampled_data <- full_data[i, ]

  rlm_full <- rlm(y ~ (var1 + var2 + var3) * gender, data = sampled_data)
  rlm_nointeraction <- rlm(y ~ var1 + var2 + var3 + gender, data = sampled_data)

  pred_full <- predict(rlm_full, newdata = full_data)
  pred_nointeraction <- predict(rlm_nointeraction, newdata = full_data)

  rmse_full <- RMSE(pred_full, full_data$y)
  rmse_nointeraction <- RMSE(pred_nointeraction, full_data$y)
  return(rmse_full - rmse_nointeraction)
}

rlm_boot <- boot(Data, statistic = simple_bootstrap, R = 500, strata = Data$gender)

# Confidence interval for the improvement from the full model, compared to the one with no interaction:
boot.ci(rlm_boot, type = "perc")

结论

以上一项或其中一项是合适的。当我怀疑方差中的方差时，我通常会认为 Bootstrap 是推理的一个重要方面。例如，即使您使用 nlme::gls 也可以使用它。 Bootstrap 更加强大，并且使许多用于处理特定情况的旧命名统计测试变得多余。