4

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

我读过，当两组之间的样本量和方差不相等时，建议进行韦尔奇的 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，我会使用suetandtest函数来执行测试。

有谁知道如何为 R 中的斜率编写 Welch 的 t 检验？

4

1 回答 1

4

我不会完全准确地回答您的问题，而是更一般的问题，即在 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)

然而，这个假设是不可接受的，所以我们需要一个替代方案。我认为这个关于交叉验证的讨论很有用。

选项 1 - 加权最小二乘

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

# 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 估计，它被设计为对数据内组中的不等方差具有鲁棒性。比较两个模型的稳健回归的良好做法是选择某种验证统计量并使用引导程序来查看平均而言哪个模型在该统计量上表现更好。

这有点复杂，但这里有一个模拟数据的工作示例：

# 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")

结论

上述一种或一种将是合适的。当我怀疑方差中的方差时，我通常会认为引导程序是推理的一个重要方面。即使您使用它也可以使用它nlme::gls。bootstrap 更加健壮，并且使许多用于处理特定情况的旧命名统计测试变得多余。

于 2017-01-14T06:22:27.897 回答