Примечание: этот вопрос является репостом, так как мой предыдущий вопрос пришлось удалить по юридическим причинам.
Сравнивая PROC MIXED из SAS с функцией lme
из nlme
пакета в R, я наткнулся на некоторые довольно запутанные различия. Более конкретно, степени свободы в разных тестах различаются между PROC MIXED
и lme
, и я задавался вопросом, почему.
Начните со следующего набора данных (код R приведен ниже):
- ind: коэффициент, обозначающий человека, где проводится измерение
- fac: орган, где производится измерение
- trt: фактор, указывающий на лечение
- у: некоторая непрерывная переменная ответа
Идея состоит в том, чтобы построить следующие простые модели:
y ~ trt + (ind)
: ind
как случайный фактор
y ~ trt + (fac(ind))
: fac
вложенный ind
как случайный фактор
Обратите внимание, что последняя модель должна вызывать особенности, так как y
для каждой комбинации ind
и есть только 1 значение fac
.
Первая модель
В SAS я строю следующую модель:
PROC MIXED data=Data;
CLASS ind fac trt;
MODEL y = trt /s;
RANDOM ind /s;
run;
Согласно учебным пособиям, та же модель с использованием R nlme
должна быть:
> require(nlme)
> options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))
> m2<-lme(y~trt,random=~1|ind,data=Data)
Обе модели дают одинаковые оценки для коэффициентов и их SE, но при проведении F-теста на эффект trt
они используют различное количество степеней свободы:
SAS :
Type 3 Tests of Fixed Effects
Effect Num DF Den DF F Value Pr > F
trt 1 8 0.89 0.3724
R :
> anova(m2)
numDF denDF F-value p-value
(Intercept) 1 8 70.96836 <.0001
trt 1 6 0.89272 0.3812
Вопрос 1: В чем разница между обоими тестами? Оба приспособлены, используя REML, и используют те же самые контрасты.
ПРИМЕЧАНИЕ: я пробовал разные значения для опции DDFM = (включая BETWITHIN, которая теоретически должна давать те же результаты, что и lme)
Вторая модель
В САС:
PROC MIXED data=Data;
CLASS ind fac trt;
MODEL y = trt /s;
RANDOM fac(ind) /s;
run;
Эквивалентная модель в R должна быть:
> m4<-lme(y~trt,random=~1|ind/fac,data=Data)
В этом случае есть некоторые очень странные различия:
- R подходит без жалоб, в то время как SAS отмечает, что последний гессиан не является положительно определенным (что меня немного не удивляет, см. Выше)
- SE на коэффициенты отличаются (меньше в SAS)
- Опять же, F-тест использовал другое количество DF (фактически, в SAS это количество = 0)
Выход SAS:
Effect trt Estimate Std Error DF t Value Pr > |t|
Intercept 0.8863 0.1192 14 7.43 <.0001
trt Cont -0.1788 0.1686 0 -1.06 .
R выход:
> summary(m4)
...
Fixed effects: y ~ trt
Value Std.Error DF t-value p-value
(Intercept) 0.88625 0.1337743 8 6.624963 0.0002
trtCont -0.17875 0.1891855 6 -0.944840 0.3812
...
(Обратите внимание, что в этом случае тесты F и T эквивалентны и используют один и тот же DF.)
Интересно, что при использовании lme4
в R модель даже не подходит:
> require(lme4)
> m4r <- lmer(y~trt+(1|ind/fac),data=Data)
Error in function (fr, FL, start, REML, verbose) :
Number of levels of a grouping factor for the random effects
must be less than the number of observations
Вопрос 2 : В чем разница между этими моделями с вложенными факторами? Правильно ли они указаны, и если да, то почему результаты так отличаются?
Имитация данных в R:
Data <- structure(list(y = c(1.05, 0.86, 1.02, 1.14, 0.68, 1.05, 0.22,
1.07, 0.46, 0.65, 0.41, 0.82, 0.6, 0.49, 0.68, 1.55), ind = structure(c(1L,
2L, 3L, 1L, 3L, 4L, 4L, 2L, 5L, 6L, 7L, 8L, 6L, 5L, 7L, 8L), .Label = c("1",
"2", "3", "4", "5", "6", "7", "8"), class = "factor"), fac = structure(c(1L,
1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L), .Label = c("l",
"r"), class = "factor"), trt = structure(c(2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Cont",
"Treat"), class = "factor")), .Names = c("y", "ind", "fac", "trt"
), row.names = c(NA, -16L), class = "data.frame")
Имитация данных:
y ind fac trt
1.05 1 l Treat
0.86 2 l Treat
1.02 3 l Treat
1.14 1 r Treat
0.68 3 r Treat
1.05 4 l Treat
0.22 4 r Treat
1.07 2 r Treat
0.46 5 r Cont
0.65 6 l Cont
0.41 7 l Cont
0.82 8 l Cont
0.60 6 r Cont
0.49 5 l Cont
0.68 7 r Cont
1.55 8 r Cont
источник
Ответы:
Для первого вопроса, метод по умолчанию в SAS, чтобы найти df, не очень умен; он ищет термины в случайном эффекте, которые синтаксически включают фиксированный эффект, и использует их. В этом случае, так
trt
как не найден вind
, он не делает правильные вещи. Я никогда не пробовалBETWITHIN
и не знаю деталей, но либо опция Satterthwaite (satterth
), либо использованиеind*trt
в качестве случайного эффекта дают правильные результаты.Что касается второго вопроса, ваш код SAS не совсем соответствует вашему R-коду; он имеет только термин для
fac*ind
, в то время как код R имеет термин для обоихind
иfac*ind
. (См. Выходные данные Variance Components, чтобы увидеть это.) Добавление этого дает одинаковый SE дляtrt
всех моделей в Q1 и Q2 (0.1892).Как вы заметили, это странная модель для подбора, так как
fac*ind
термин имеет одно наблюдение для каждого уровня, поэтому эквивалентен термину ошибки. Это отражено в результатах SAS, гдеfac*ind
термин имеет нулевую дисперсию. Это также то, что сообщение об ошибке от lme4 говорит вам; причина ошибки заключается в том, что вы, скорее всего, что-то неправильно определили, поскольку вы включаете термин ошибки в модель двумя различными способами. Интересно, что есть одна небольшая разница в модели NLME; он каким-то образом находит член отклонения дляfac*ind
члена в дополнение к члену ошибки, но вы заметите, что сумма этих двух отклонений равна члену ошибки как из SAS, так и из nlme безfac*ind
члена. Тем не менее, SE дляtrt
остается таким же (0,1892), как вложеноtrt
вind
таким образом, эти более низкие условия отклонения не влияют на это.Наконец, общее замечание о степенях свободы в этих моделях: они рассчитываются после подбора модели, и поэтому различия в степенях свободы между различными программами или опциями программы не обязательно означают, что модель подгоняется по-разному. Для этого нужно взглянуть на оценки параметров, как параметров фиксированного эффекта, так и параметров ковариации.
Кроме того, с помощью t- и F приближения с заданным числом степеней свободы является довольно спорным. Мало того, что есть несколько способов приблизить df, некоторые полагают, что практика в любом случае это не очень хорошая идея. Несколько слов совета:
Если все сбалансировано, сравните результаты с традиционным методом наименьших квадратов, как они должны согласиться. Если он близок к сбалансированному, вычислите его самостоятельно (при условии баланса), чтобы вы могли убедиться, что те, которые вы используете, находятся в правильном поле.
Если у вас большой размер выборки, степени свободы не имеют большого значения, так как распределения приближаются к нормальному и хи-квадрат.
Проверьте методы Дуга Бейтса для вывода. Его старый метод основан на моделировании MCMC; его новый метод основан на профилировании вероятности.
источник