Альтернативы одностороннему ANOVA для гетероскедастических данных

36

У меня есть данные от 3 групп биомассы водорослей ( , , ), которые содержат неравные размеры выборки ( , , ), и я хотел бы сравнить, если эти группы принадлежат к одной популяции.B C n A = 15 n B = 13 n C = 12ABCnA=15nB=13nC=12

Односторонний ANOVA определенно был бы подходящим вариантом, однако при проведении тестов нормальности на моих данных основной проблемой является гетероскедотность. Мои необработанные данные без какого-либо преобразования дали соотношение отклонений ( ), которое намного выше критического значения ( F _ {\ rm crit} = 4,16 ), и поэтому я не могу выполнить односторонний ANOVA ,Fmax=19.1Fcrit=4.16

Я также попытался преобразование, чтобы нормализовать мои данные. Даже после испытаний различных преобразований (log, square root, square) самое низкое значение Fmax полученное после преобразования с преобразованием log10 составляло 7.16 , что было все еще выше по сравнению с Fcrit .

Может кто-нибудь здесь посоветовать мне, куда идти отсюда? Я не могу думать о других методах преобразования, чтобы нормализовать данные. Есть ли альтернативы односторонней ANOVA?

PS: мои необработанные данные ниже:

A: 0.178 0.195 0.225 0.294 0.315 0.341 0.36  0.363 0.371 0.398 0.407 0.409 0.432 
   0.494 0.719
B: 0.11  0.111 0.204 0.416 0.417 0.441 0.492 0.965 1.113 1.19  1.233 1.505 1.897
C: 0.106 0.114 0.143 0.435 0.448 0.51  0.576 0.588 0.608 0.64  0.658 0.788 0.958
Рик Л.
источник
2
F-тест ясно показывает, что группы не принадлежат к одной и той же популяции. Следующим шагом будет интерпретация этого результата, начиная с четкой визуализации и количественного описания данных с разбивкой по группам.
whuber
Есть также тесты перестановки в пакете монеты. Для ANOVA это будет функция oneway_test. Пример Quick-R: statmethods.net/stats/resampling.html
user16263

Ответы:

73

Есть несколько вариантов, доступных при работе с гетероскедастическими данными. К сожалению, ни один из них не гарантированно всегда работает. Вот несколько вариантов, с которыми я знаком:

  • преобразования
  • Уэлч АНОВА
  • взвешенные наименьшие квадраты
  • устойчивая регрессия
  • гетероскедастичность соответствует стандартным ошибкам
  • начальная загрузка
  • Тест Крускала-Уоллиса
  • порядковая логистическая регрессия

Обновление: Вот демонстрация R некоторых способов подгонки линейной модели (например, ANOVA или регрессии), когда у вас есть гетероскедастичность / гетерогенность дисперсии.

Давайте начнем с просмотра ваших данных. Для удобства я загрузил их в два названных фрейма данных my.data(который структурирован, как указано выше, с одним столбцом на группу) и stacked.data(который имеет два столбца: valuesс числами и indс индикатором группы).

Мы можем официально проверить на гетероскедастичность с помощью теста Левена:

library(car)
leveneTest(values~ind, stacked.data)
# Levene's Test for Homogeneity of Variance (center = median)
#       Df F value   Pr(>F)   
# group  2  8.1269 0.001153 **
#       38                    

Конечно же, у вас есть гетероскедастичность. Мы проверим, каковы дисперсии групп. Эмпирическое правило заключается в том, что линейные модели достаточно устойчивы к неоднородности дисперсии, если максимальная дисперсия не более чем в Раза превышает минимальную дисперсию, поэтому мы также найдем это соотношение: 4×

apply(my.data, 2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01734578 0.33182844 0.06673060 
var(my.data$B, na.rm=T) / var(my.data$A, na.rm=T)
# [1] 19.13021

Ваши дисперсии существенно отличаются, с самым большим, B, будучи самый маленький . Это проблемный уровень гетероскедсатичности. 19×A

Вы думали использовать преобразования, такие как лог или квадратный корень, чтобы стабилизировать дисперсию. В некоторых случаях это будет работать, но преобразования типа Бокса-Кокса стабилизируют дисперсию, сжимая данные асимметрично, либо сжимая их вниз с максимальными сжимаемыми данными, либо сжимая их вверх с наименьшими сжимаемыми данными наиболее. Таким образом, вам нужно, чтобы дисперсия ваших данных изменялась со средним значением для оптимальной работы. Ваши данные имеют огромную разницу в дисперсии, но относительно небольшую разницу между средними и средними значениями, т. Е. Распределения в основном перекрываются. В качестве учебного упражнения мы можем создать некоторые parallel.universe.data, добавив ко всем значениям и к.72.7B.7CЧтобы показать, как это будет работать:

parallel.universe.data = with(my.data, data.frame(A=A, B=B+2.7, C=C+.7))
apply(parallel.universe.data,       2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01734578 0.33182844 0.06673060 
apply(log(parallel.universe.data),  2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.12750634 0.02631383 0.05240742 
apply(sqrt(parallel.universe.data), 2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01120956 0.02325107 0.01461479 
var(sqrt(parallel.universe.data$B), na.rm=T) / 
var(sqrt(parallel.universe.data$A), na.rm=T)
# [1] 2.074217

Использование преобразования квадратного корня довольно хорошо стабилизирует эти данные. Вы можете увидеть улучшение для данных параллельной вселенной здесь:

введите описание изображения здесь

Вместо того, чтобы просто пытаться использовать разные преобразования, более систематический подход заключается в оптимизации параметра Бокса-Кокса (хотя обычно рекомендуется округлить его до ближайшего интерпретируемого преобразования). В вашем случае допустимы либо квадратный корень , либо лог , хотя на самом деле ни один из них не работает. Для данных параллельной вселенной лучше использовать квадратный корень: λ = .5 λ = 0λλ=.5λ=0

boxcox(values~ind, data=stacked.data,    na.action=na.omit)
boxcox(values~ind, data=stacked.pu.data, na.action=na.omit) 

введите описание изображения здесь

Поскольку этот случай представляет собой ANOVA (т. Е. Не имеет непрерывных переменных), один из способов справиться с неоднородностью - это использовать поправку Уэлча к знаменателям степеней свободы в тесте (nb , вместо дробного значения ): Fdf = 19.445df = 38

oneway.test(values~ind, data=stacked.data, na.action=na.omit, var.equal=FALSE)
#  One-way analysis of means (not assuming equal variances)
# 
# data:  values and ind
# F = 4.1769, num df = 2.000, denom df = 19.445, p-value = 0.03097

Более общий подход заключается в использовании взвешенных наименьших квадратов . Поскольку некоторые группы ( B) разбросаны больше, данные в этих группах предоставляют меньше информации о расположении среднего значения, чем данные в других группах. Мы можем позволить модели включить это, предоставляя вес для каждой точки данных. Обычной системой является использование обратной величины групповой дисперсии в качестве веса:

wl = 1 / apply(my.data, 2, function(x){ var(x, na.rm=T) })
stacked.data$w = with(stacked.data, ifelse(ind=="A", wl[1], 
                                           ifelse(ind=="B", wl[2], wl[3])))
w.mod = lm(values~ind, stacked.data, na.action=na.omit, weights=w)
anova(w.mod)
# Response: values
#           Df Sum Sq Mean Sq F value  Pr(>F)  
# ind        2   8.64  4.3201  4.3201 0.02039 *
# Residuals 38  38.00  1.0000                  

Это дает немного отличающиеся значения и от невзвешенного ANOVA ( , ), но оно хорошо учитывает неоднородность: pFp4.50890.01749

введите описание изображения здесь

Однако взвешенные наименьшие квадраты не являются панацеей. Один неприятный факт заключается в том, что весовые коэффициенты верны, и это означает, что они известны априори. Он также не учитывает ненормальность (например, перекос) или выбросы. Использование весов, оцененных по вашим данным, часто будет работать нормально, особенно если у вас достаточно данных для оценки отклонения с разумной точностью (это аналогично идее использования таблицы вместо таблицы, когда у вас есть илиt 50 100 Nzt50100степеней свободы), ваши данные достаточно нормальны, и у вас, похоже, нет никаких выбросов. К сожалению, у вас относительно мало данных (13 или 15 на группу), некоторые искажены и, возможно, некоторые выбросы. Я не уверен, что они достаточно плохи, чтобы из них что-то сделать, но вы могли бы смешать взвешенные наименьшие квадраты с помощью надежных методов. Вместо того, чтобы использовать дисперсию в качестве меры разброса (которая чувствительна к выбросам, особенно с низким ), вы можете использовать обратную величину интерквартильного диапазона (на которую не влияют до 50% выбросов в каждой группе). Эти веса можно затем объединить с устойчивой регрессией, используя другую функцию потерь, такую ​​как биквадрат Тьюки: N

1 / apply(my.data,       2, function(x){ var(x, na.rm=T) })
#         A         B         C 
# 57.650907  3.013606 14.985628 
1 / apply(my.data,       2, function(x){ IQR(x, na.rm=T) })
#        A        B        C 
# 9.661836 1.291990 4.878049 

rw = 1 / apply(my.data, 2, function(x){ IQR(x, na.rm=T) })
stacked.data$rw = with(stacked.data, ifelse(ind=="A", rw[1], 
                                            ifelse(ind=="B", rw[2], rw[3])))
library(robustbase)
w.r.mod = lmrob(values~ind, stacked.data, na.action=na.omit, weights=rw)
anova(w.r.mod, lmrob(values~1,stacked.data,na.action=na.omit,weights=rw), test="Wald")
# Robust Wald Test Table
# 
# Model 1: values ~ ind
# Model 2: values ~ 1
# Largest model fitted by lmrob(), i.e. SM
# 
#   pseudoDf Test.Stat Df Pr(>chisq)  
# 1       38                          
# 2       40    6.6016  2    0.03685 *

Вес здесь не такой экстремальный. Прогнозируемые средства группы незначительно отличаются ( A: WLS 0.36673, надежным 0.35722; B: ВМНК 0.77646, надежным 0.70433; C: WLS 0.50554, надежный 0.51845), с помощью Bи Cтого меньше тянут экстремальных значения.

В эконометрике очень популярна стандартная ошибка Хубера-Уайта («сэндвич») . Как и поправка Уэлча, это не требует, чтобы вы знали априорные отклонения, и не требует, чтобы вы оценивали веса по вашим данным и / или зависели от модели, которая может быть неверной. С другой стороны, я не знаю, как включить это в ANOVA, а это означает, что вы получаете их только для тестов отдельных фиктивных кодов, что в данном случае кажется мне менее полезным, но я все равно продемонстрирую их:

library(sandwich)
mod = lm(values~ind, stacked.data, na.action=na.omit)
sqrt(diag(vcovHC(mod)))
# (Intercept)        indB        indC 
#  0.03519921  0.16997457  0.08246131 
2*(1-pt(coef(mod) / sqrt(diag(vcovHC(mod))), df=38))
#  (Intercept)         indB         indC 
# 1.078249e-12 2.087484e-02 1.005212e-01 

Функция vcovHCрассчитывает гетероскедастичную согласованную дисперсионно-ковариационную матрицу для ваших бета-версий (фиктивных кодов), то есть то, что обозначают буквы в функции. Чтобы получить стандартные ошибки, вы извлекаете основную диагональ и берете квадратные корни. Чтобы получить тесты для ваших бета-версий, вы делите оценки коэффициентов на SE и сравниваете результаты с соответствующим распределением (а именно, распределением с вашими остаточными степенями свободы). т тttt

RСпециально для пользователей @TomWenseleers отмечает в комментариях ниже, что функция ? Anova в carпакете может принимать white.adjustаргумент для получения значения для фактора, использующего ошибки, согласованные с гетероскедастичностью. p

Anova(mod, white.adjust=TRUE)
# Analysis of Deviance Table (Type II tests)
# 
# Response: values
#           Df      F  Pr(>F)  
# ind        2 3.9946 0.02663 *
# Residuals 38                 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Вы можете попытаться получить эмпирическую оценку того, как выглядит фактическое распределение выборки вашей тестовой статистики с помощью начальной загрузки . Во-первых, вы создаете истинный ноль, делая все группы равными. Затем вы производите повторную выборку с заменой и рассчитываете свою статистику теста ( ) для каждой начальной загрузки, чтобы получить эмпирическую оценку распределения выборки под нулевым значением с вашими данными независимо от их статуса в отношении нормальности или однородности. Доля того распределения выборки, которое является таким же экстремальным или более экстремальным, чем наблюдаемая вами статистика теста, является значением: F pFFp

mod    = lm(values~ind, stacked.data, na.action=na.omit)
F.stat = anova(mod)[1,4]
  # create null version of the data
nullA  = my.data$A - mean(my.data$A)
nullB  = my.data$B - mean(my.data$B, na.rm=T)
nullC  = my.data$C - mean(my.data$C, na.rm=T)
set.seed(1)
F.vect = vector(length=10000)
for(i in 1:10000){
  A = sample(na.omit(nullA), 15, replace=T)
  B = sample(na.omit(nullB), 13, replace=T)
  C = sample(na.omit(nullC), 13, replace=T)
  boot.dat  = stack(list(A=A, B=B, C=C))
  boot.mod  = lm(values~ind, boot.dat)
  F.vect[i] = anova(boot.mod)[1,4]
}
1-mean(F.stat>F.vect)
# [1] 0.0485

В некотором смысле, самозагрузка является окончательным подходом с уменьшенным допущением для проведения анализа параметров (например, средних), но она предполагает, что ваши данные являются хорошим представлением совокупности, то есть у вас есть разумный размер выборки. Так как ваши маленькие, это может быть менее надежным. Вероятно, окончательная защита от ненормальности и неоднородности состоит в использовании непараметрического теста. Основной непараметрической версией ANOVA является тест Крускала-Уоллиса : n

kruskal.test(values~ind, stacked.data, na.action=na.omit)
#  Kruskal-Wallis rank sum test
# 
# data:  values by ind
# Kruskal-Wallis chi-squared = 5.7705, df = 2, p-value = 0.05584

Хотя тест Крускала-Уоллиса, безусловно, является лучшей защитой от ошибок типа I, его можно использовать только с одной категориальной переменной (т. Е. Без непрерывных предикторов или факторных планов), и он имеет наименьшую мощность из всех обсуждаемых стратегий. Другой непараметрический подход заключается в использовании порядковой логистической регрессии . Многим это кажется странным, но вам нужно только предположить, что ваши ответные данные содержат достоверную порядковую информацию, что они, безусловно, делают, иначе любая другая стратегия, указанная выше, также недействительна:

library(rms)
olr.mod = orm(values~ind, stacked.data)
olr.mod
# Model Likelihood          Discrimination          Rank Discrim.    
# Ratio Test                 Indexes                Indexes       
# Obs            41    LR chi2      6.63    R2                  0.149    rho     0.365
# Unique Y       41    d.f.            2    g                   0.829
# Median Y    0.432    Pr(> chi2) 0.0363    gr                  2.292
# max |deriv| 2e-04    Score chi2   6.48    |Pr(Y>=median)-0.5| 0.179
#                      Pr(> chi2) 0.0391

Это может быть не ясно из результатов, но проверка модели в целом, которая в данном случае является проверкой ваших групп, является chi2ниже Discrimination Indexes. Перечислены две версии: тест отношения правдоподобия и тест оценки. Тест отношения правдоподобия обычно считается лучшим. Это дает значение . p0.0363

Gung - Восстановить Монику
источник
1
Любите ваши ответы! Теперь есть проблеск дома, чтобы иметь дело с моими данными! :) Из всех этих подходов, которые вы бы порекомендовали для моих данных, особенно с точки зрения предоставления достаточной статистической мощности? Я не слишком заинтересован в использовании непараметрического подхода, то есть теста Крускала-Уоллиса, так как это может означать, что я могу допустить ошибку типа I? Поправьте меня если я ошибаюсь.
Рик Л.
2
Здесь я поработаю над некоторыми материалами, когда у меня будет такая возможность, но нет, тест Крускала-Уоллиса, вероятно, обеспечит вам наилучшую защиту от ошибок I типа. ОТО, это может быть не самый мощный вариант у вас есть.
gung - Восстановить Монику
4
Обратите внимание, что в библиотеке carтакже есть опция для white.adjust=Tработы с гетероскелетом с использованием скорректированных по белому стандарту ошибок, скорректированных на гетероскедастичность
Том Венселерс
3
Да , это только упоминает его lm«с, но это также , кажется, работает для aov» s (варианты white.adjustявляются white.adjust=c(FALSE, TRUE, "hc3", "hc0", "hc1", "hc2", "hc4")- для получения дополнительной информации см ?hccm)
Том Wenseleers
1
@Nakx, нет, не то, что я знаю. Это эмпирическое правило; вряд ли это будет в статье в статистической литературе. Это возможно в некоторых вводных учебниках.
gung - Восстановить Монику