Я - пользователь, более знакомый с R, и пытался оценить случайные уклоны (коэффициенты отбора) примерно для 35 особей в течение 5 лет для четырех переменных среды обитания. Переменная ответа - является ли место «использованным» (1) или «доступным» (0) местом обитания («использование» ниже).
Я использую Windows 64-битный компьютер.
В версии R 3.1.0 я использую данные и выражения ниже. PS, TH, RS и HW - фиксированные эффекты (стандартизированное, измеренное расстояние до типов среды обитания). Ime4 V 1,1-7.
str(dat)
'data.frame': 359756 obs. of 7 variables:
$ use : num 1 1 1 1 1 1 1 1 1 1 ...
$ Year : Factor w/ 5 levels "1","2","3","4",..: 4 4 4 4 4 4 4 4 3 4 ...
$ ID : num 306 306 306 306 306 306 306 306 162 306 ...
$ PS: num -0.32 -0.317 -0.317 -0.318 -0.317 ...
$ TH: num -0.211 -0.211 -0.211 -0.213 -0.22 ...
$ RS: num -0.337 -0.337 -0.337 -0.337 -0.337 ...
$ HW: num -0.0258 -0.19 -0.19 -0.19 -0.4561 ...
glmer(use ~ PS + TH + RS + HW +
(1 + PS + TH + RS + HW |ID/Year),
family = binomial, data = dat, control=glmerControl(optimizer="bobyqa"))
glmer дает мне оценки параметров для фиксированных эффектов, которые имеют смысл для меня, и случайные наклоны (которые я интерпретирую как коэффициенты отбора для каждого типа среды обитания) также имеют смысл, когда я качественно исследую данные. Логарифмическая вероятность для модели составляет -3050,8.
Тем не менее, большинство исследований в области экологии животных не используют R, потому что с данными о местонахождении животных пространственная автокорреляция может сделать стандартные ошибки склонными к ошибке I типа. В то время как R использует стандартные ошибки, основанные на модели, предпочтительными являются эмпирические (также по Хубер-Уайту или сэндвичу) стандартные ошибки.
Хотя R в настоящее время не предлагает эту опцию (насколько мне известно - ПОЖАЛУЙСТА, исправьте меня, если я ошибаюсь), SAS делает - хотя у меня нет доступа к SAS, коллега согласился позволить мне одолжить его компьютер, чтобы определить, есть ли стандартные ошибки значительно измениться при использовании эмпирического метода.
Во-первых, мы хотели убедиться, что при использовании стандартных ошибок на основе модели SAS будет давать оценки, аналогичные R - чтобы быть уверенным, что модель указана одинаково в обеих программах. Мне все равно, если они точно так же - просто похожи. Я попробовал (SAS V 9.2):
proc glimmix data=dat method=laplace;
class year id;
model use = PS TH RS HW / dist=bin solution ddfm=betwithin;
random intercept PS TH RS HW / subject = year(id) solution type=UN;
run;title;
Я также пробовал различные другие формы, такие как добавление строк
random intercept / subject = year(id) solution type=UN;
random intercept PS TH RS HW / subject = id solution type=UN;
Я пытался без указания
solution type = UN,
или комментируя
ddfm=betwithin;
Независимо от того, как мы определяем модель (и мы пробовали много способов), я не могу получить случайные наклоны в SAS, чтобы отдаленно напоминать те выходные данные из R - даже если фиксированные эффекты достаточно похожи. И когда я имею в виду другое, я имею в виду, что даже знаки не одинаковы. Логарифмическая правдоподобие -2 в SAS было 71344,94.
Я не могу загрузить свой полный набор данных; поэтому я сделал игрушечный набор данных только с записями трех человек. SAS дает мне вывод через несколько минут; в R это занимает более часа. Weird. С этим набором игрушек я теперь получаю разные оценки для фиксированных эффектов.
Мой вопрос: может ли кто-нибудь пролить свет на то, почему оценки случайных наклонов могут быть такими разными между R и SAS? Могу ли я что-нибудь сделать в R или SAS, чтобы изменить мой код так, чтобы вызовы давали похожие результаты? Я бы лучше изменил код в SAS, так как я «верю» моим оценкам R больше.
Я действительно обеспокоен этими различиями и хочу докопаться до сути этой проблемы!
Мой вывод из игрушечного набора данных, который использует только три из 35 человек в полном наборе данных для R и SAS, включен в виде JPEG.
РЕДАКТИРОВАТЬ И ОБНОВИТЬ:
Как помогло обнаружение @JakeWestfall, уклоны в SAS не включают в себя фиксированные эффекты. Когда я добавляю фиксированные эффекты, вот результат - сравнение уклонов R с уклонами SAS для одного фиксированного эффекта, «PS», между программами: (Коэффициент выбора = случайный уклон). Обратите внимание на увеличение вариации SAS.
ID
это не фактор в R; проверьте и посмотрите, изменит ли это что-нибудь.0
s и1
s,R
будет моделироваться вероятность ответа «1», в то время как SAS будет моделировать вероятность ответа «0». Чтобы сделать модель SAS вероятностью «1», вам нужно записать переменную ответа какuse(event='1')
. Конечно, даже не делая этого, я полагаю, что мы все же должны ожидать те же оценки дисперсий случайных эффектов, а также те же оценки фиксированных эффектов, хотя их знаки обратны.ranef()
функцию, а неcoef()
. Первый дает фактические случайные эффекты, а второй дает случайные эффекты плюс вектор с фиксированными эффектами. Это объясняет, почему цифры, показанные в вашем посте, отличаются друг от друга, но все еще остается существенное расхождение, которое я не могу полностью объяснить.Ответы:
Похоже, что я не должен был ожидать, что случайные наклоны будут одинаковыми между пакетами, в соответствии с Zhang et al 2011. В их статье « Подгонка обобщенных линейных моделей смешанных эффектов для двоичных ответов с использованием различных статистических пакетов» , они описывают:
Абстрактный:
Я надеюсь, что @BenBolker и команда рассмотрят мой вопрос как голосование за R, чтобы включить эмпирические стандартные ошибки и возможность квадратуры Гаусса-Эрмита для моделей с несколькими случайными коэффициентами наклона к блеску, поскольку я предпочитаю интерфейс R и хотел бы иметь возможность применять некоторые дальнейшие анализы в этой программе. К счастью, даже при том, что R и SAS не имеют сопоставимых значений для случайных наклонов, общие тенденции одинаковы. Спасибо всем за ваш вклад, я действительно ценю время и внимание, которое вы вложили в это!
источник
Смесь ответа и комментария / больше вопросов:
Я установил «игрушечный» набор данных с тремя различными вариантами оптимизатора. (* Примечание 1: вероятно, было бы более полезно для сравнительных целей сделать небольшой набор данных путем подвыборки из каждого года и идентификатора, а не путем подбора выборки переменных группировки. Как это, мы знаем, что GLMM не будет работать особенно хорошо с таким небольшим количеством уровней переменных группировки. Вы можете сделать это через что-то вроде:
Пакетный фитинг код:
Затем я прочитал результаты в новой сессии:
Прошедшее время и отклонение:
Эти отклонения значительно ниже отклонения, сообщенного OP от R (6101,7), и немного ниже отклонений, сообщенных OP от SAS (6078,9), хотя сравнение отклонений по пакетам не всегда целесообразно.
Я был действительно удивлен, что SAS сошлись всего в 100 оценках функций!
Времена варьируются от 17 минут (
nloptwrap
) до 80 минут (bobyqa
) на Macbook Pro, что соответствует опыту OP. Отклонение немного лучше дляnloptwrap
.Ответы выглядят совершенно иначе,
nloptwrap
хотя стандартные ошибки довольно велики ...(код здесь дает некоторые предупреждения о
year:id
том, что я должен выследить)Продолжение следует ... ?
источник