Dav*_*rit 5 regression r linear-regression t-test
我对具有相同自变量的两组进行回归。然后，我想测试两个回归的斜率是否显着不同。
我读到，当两组之间的样本量和方差不相等时，建议进行韦尔奇 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 中编写斜率的韦尔奇 t 检验代码？
我不太会准确地回答你的问题，而是更普遍的问题是，在 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)
然而，这个假设是不可接受的，所以我们需要一个替代方案。我认为关于交叉验证的讨论很有用。
我不确定这是否与韦尔奇 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
第二种选择是使用 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即使你使用例如它也可以使用。引导程序更加强大，并且使许多用于处理特定情况的旧命名统计测试变得多余。