Как проверить разницу в средствах между 5 группами? Разница между группами не одинакова. Если бы они были равны, я мог бы использовать ANOVA [migrated]

Я хотел бы рекомендовать обобщенные наименьшие квадраты в качестве опции. Точно так же вы можете создать модель для среднего значения данных с учетом ваших предикторов, вы также можете создать модель для отклонения. Функция gls() в пакете nlme позволяет нам это сделать.

Я продемонстрирую, как показано ниже:

set.seed(1984)
n <- 25
# Create heteroskedastic data, balanced design
dat <- data.frame(
  y = c(
    rnorm(n, 1, 1), rnorm(n, 0, 1), rnorm(n, .5, 1), 
    rnorm(n, 0, 3), rnorm(n, -.5, 4)),
  groups = sort(rep(LETTERS[1:5], n))
)
boxplot(y ~ groups, dat)

# Assuming all cases have identical variances
summary(aov(y ~ groups, dat))

             Df Sum Sq Mean Sq F value Pr(>F)  
groups        4   69.0  17.255   3.475 0.0101 *
Residuals   120  595.9   4.966                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Modeling the error variances
library(nlme)
# This is your standard linear model
coef(summary(fit.gls.baseline <- gls(y ~ groups, dat))))

                 Value Std.Error   t-value      p-value
(Intercept)  1.1784414 0.4457005  2.644021 0.0092875812
groupsB     -0.9748277 0.6303157 -1.546571 0.1246000669
groupsC     -0.8980932 0.6303157 -1.424831 0.1568017302
groupsD     -1.2903790 0.6303157 -2.047195 0.0428223361
groupsE     -2.3076197 0.6303157 -3.661054 0.0003747902

sigma(fit.gls.baseline) # get residual standard deviation
[1] 2.228502

# Next we fit a heteroskedastic model, hidden some output
# varIdent says the variances are identical for cases that have the same value of group
summary(fit.gls <- gls(
  y ~ groups, dat, weights = varIdent(form = ~ 1 | groups)))

Generalized least squares fit by REML
  Model: y ~ groups 
  Data: dat 
       AIC      BIC    logLik
  479.5354 507.4103 -229.7677

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | groups 
 Parameter estimates:
        A         B         C         D         E 
1.0000000 0.9741804 0.8698579 3.0062743 3.7867899 

Coefficients:
                 Value Std.Error   t-value p-value
(Intercept)  1.1784414 0.1951410  6.038923  0.0000
groupsB     -0.9748277 0.2724316 -3.578248  0.0005
groupsC     -0.8980932 0.2586375 -3.472401  0.0007
groupsD     -1.2903790 0.6182517 -2.087142  0.0390
groupsE     -2.3076197 0.7642898 -3.019299  0.0031

Residual standard error: 0.975705 

Несколько заметок об этой модели. Вы найдете, что остаточное стандартное отклонение (ошибка) намного меньше для гетероскедастической модели по сравнению с той же оценкой для более ранней модели. Теперь взгляните на раздел функции Variance, и вы увидите, что группа A имеет стандартное отклонение 1. Эти значения являются относительными стандартными отклонениями. Таким образом, $ .976 $ для модели - это остаточная SD группы A. $ .974 \ times .976 $ - значение для группы B. Группы D $ (3 \ times .976) $ и E $ (3.79 \ раз .976) $ имеют гораздо большие остаточные стандартные отклонения.

Также давайте получим гетероскедастическую версию F -test:

# marginal for type III sum of squares, only makes a difference for the
# test of the grouping variable if we have more than one grouping variable
anova(fit.gls, type = "marginal")
Denom. DF: 120 
            numDF  F-value p-value
(Intercept)     1 36.46859  <.0001
groups          4  5.51477   4e-04

Достаточно что не все группы имеют одно и то же среднее значение, $ F (4, 120) = 5,51, p & lt; .001 $.

. Теперь, чтобы сравнить гетероскедастическую модель с стандартной моделью, обратите внимание, что коэффициенты примерно одинаковы, однако стандартные ошибки намного меньше для гетероскедастической модели, за исключением групп D и E . Значит, наши p -значения также меньше. Глобальный подход к сопоставлению моделей был бы следующим:

anova(fit.gls, fit.gls.baseline) # Model fit improved
                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit.gls              1 10 479.5354 507.4103 -229.7677                        
fit.gls.baseline     2  6 560.9588 577.6837 -274.4794 1 vs 2 89.42343  <.0001

. Эти результаты свидетельствуют о том, что гетероскедастическая модель была лучше, ниже AIC, BIC и статистически значимым критерием отношения правдоподобия.

Теперь мы делаем не имеют доступа к процессу генерации данных в реальности, но мы всегда можем построить boxplt выше. Предполагая, что мы хотим выбрать более простую модель, мы можем предположить, что группы A, B и C имеют одинаковые дисперсии, а группы D и E - разные. Затем мы можем попробовать:

# Create a variable that classifies the groups according to variance
dat$var.groups <- (dat$groups %in% c("D", "E")) + 0
# then run:
summary(fit.gls.parsim <- gls(
  y ~ groups, dat, weights = varIdent(form = ~ 1 | var.groups)))

Generalized least squares fit by REML
  Model: y ~ groups 
  Data: dat 
       AIC      BIC    logLik
  475.3162 494.8287 -230.6581

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | var.groups 
 Parameter estimates:
       0        1 
1.000000 3.600023 

Coefficients:
                 Value Std.Error   t-value p-value
(Intercept)  1.1784414 0.1853218  6.358893  0.0000
groupsB     -0.9748277 0.2620846 -3.719516  0.0003
groupsC     -0.8980932 0.2620846 -3.426730  0.0008
groupsD     -1.2903790 0.6924235 -1.863569  0.0648
groupsE     -2.3076197 0.6924235 -3.332671  0.0011

# A model comparison between our previous best and this new model
anova(fit.gls, fit.gls.parsim)
               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit.gls            1 10 479.5354 507.4103 -229.7677                        
fit.gls.parsim     2  7 475.3162 494.8287 -230.6581 1 vs 2 1.780859  0.6191

Тест отношения правдоподобия не может отличить эту более простую модель, которая имеет еще три параметра и нашу раннюю гетероскедастическую модель. Также AIC и BIC склоняются к этой модели. Поэтому мы можем выбрать, что эта модель идет вперед, поскольку мы не знаем правды. Теория может также предложить некоторые вещи о различиях в групповых отклонениях. Можно видеть, как конкретная проблема гетероскедастичности может иметь существенный интерес сама по себе.

F -тест снова будет:

anova(fit.gls.parsim, type = "marginal")
Denom. DF: 120 
            numDF  F-value p-value
(Intercept)     1 40.43552  <.0001
groups          4  6.04201   2e-04

Теперь, следующая вещь может быть контрастом. В моей работе мне все равно, потому что я делаю здесь самые простые вещи. Есть, вероятно, лучшие способы борьбы с контрастами, см. Документацию по пакетам emmeans для получения дополнительных параметров. Я использую настройки по умолчанию.

library(emmeans)
pairs(emmeans(fit.gls.parsim, "groups"))
 contrast   estimate        SE  df t.ratio p.value
 A - B     0.9748277 0.2620846 120   3.720  0.0028
 A - C     0.8980932 0.2620846 120   3.427  0.0073
 A - D     1.2903790 0.6924235 120   1.864  0.3427
 A - E     2.3076197 0.6924235 120   3.333  0.0099
 B - C    -0.0767345 0.2620846 120  -0.293  0.9984
 B - D     0.3155513 0.6924235 120   0.456  0.9910
 B - E     1.3327920 0.6924235 120   1.925  0.3099
 C - D     0.3922858 0.6924235 120   0.567  0.9796
 C - E     1.4095265 0.6924235 120   2.036  0.2555
 D - E     1.0172407 0.9435106 120   1.078  0.8174

P value adjustment: tukey method for comparing a family of 5 estimates 

Эти контрасты используют информацию из модели гетероскедастики, поскольку результаты будут разными, если вы использовали модель fit.gls.baseline.

Надеюсь, что демонстрация выше показывает, как вы можете учитывать гетероскедастичность, используя обобщенные наименьшие квадраты при анализе дисперсии. По сравнению со многими другими подходами этот подход не просто рассматривает гетероскедастичность как неприятность. Можно также включить контрольные переменные, как в ANCOVA.

Как и в стороне, если вы также подозреваете, что в ANCOVA существует также гетеросквадичность по управляющим переменным, то функция glmmTMB() в пакете glmmTMB может справиться с этой ситуацией. Вы просто указываете уравнение регрессии для дисперсии с использованием аргумента dispformula =.