Как выполнить t-критерий Велча на наклонах из двух регрессий с R?

Я не собираюсь точно отвечать на ваш вопрос, но вместо этого я отвечу на более общий вопрос о том, как в 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 - взвешенный метод наименьших квадратов

Я не уверен, что это то же самое, что и 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 - надежная регрессия, сравнение через бутстрап

Второй вариант - использовать 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. Бутстрап более надежен и делает ненужными многие из старых названных статистических тестов для работы с конкретными ситуациями.