контекст
В этом вопросе используется R, но речь идет об общих статистических вопросах.
Я анализирую влияние факторов смертности (% смертности от болезней и паразитов) на скорость роста популяции моли с течением времени, когда популяция личинок отбиралась из 12 мест один раз в год в течение 8 лет. Данные о темпах роста населения показывают четкую, но нерегулярную циклическую тенденцию во времени.
Остатки от простой обобщенной линейной модели (скорость роста ~% заболевания +% паразитизм + год) демонстрировали аналогичную четкую, но нерегулярную циклическую тенденцию во времени. Следовательно, обобщенные модели наименьших квадратов той же формы также были адаптированы к данным с соответствующими корреляционными структурами, чтобы иметь дело с временной автокорреляцией, например, составной симметрией, авторегрессионным порядком обработки 1 и авторегрессионными корреляционными структурами скользящего среднего.
Все модели содержали одинаковые фиксированные эффекты, сравнивались с использованием AIC и устанавливались REML (чтобы позволить сравнение различных структур корреляции с помощью AIC). Я использую пакет R NLME и функцию GLS.
Вопрос 1
Остатки моделей GLS по-прежнему отображают практически идентичные циклические закономерности при построении графика в зависимости от времени. Будут ли такие шаблоны всегда оставаться, даже в моделях, которые точно учитывают автокорреляционную структуру?
Я смоделировал некоторые упрощенные, но похожие данные в R ниже моего второго вопроса, который показывает проблему, основанную на моем текущем понимании методов, необходимых для оценки временных автокоррелированных моделей в остатках модели , которые, как я теперь знаю, неверны (см. Ответ).
вопрос 2
Я приспособил модели GLS со всеми возможными вероятными структурами корреляции к моим данным, но на самом деле ни одна из них не является существенно более подходящей, чем GLM, без какой-либо структуры корреляции: только одна модель GLS незначительно лучше (оценка AIC = 1,8 ниже), в то время как все остальные имеют более высокие значения AIC. Тем не менее, это только тот случай, когда все модели соответствуют REML, а не ML, где модели GLS явно намного лучше, но я понимаю, что из книг статистики вы должны использовать REML только для сравнения моделей с различными структурами корреляции и одинаковыми фиксированными эффектами по причинам. Я не буду подробно здесь.
Учитывая явно временную автокорреляционную природу данных, если ни одна модель даже не в несколько раз лучше, чем простая GLM, какой самый подходящий способ решить, какую модель использовать для вывода, предполагая, что я использую соответствующий метод (я в конечном итоге хочу использовать AIC для сравнения различных переменных комбинаций)?
Q1 «симуляция», исследующая остаточные модели в моделях с соответствующими корреляционными структурами и без них
Создайте смоделированную переменную отклика с циклическим эффектом «время» и положительным линейным эффектом «х»:
time <- 1:50
x <- sample(rep(1:25,each=2),50)
y <- rnorm(50,5,5) + (5 + 15*sin(2*pi*time/25)) + (x/1)
у должен отображаться циклический тренд в течение «времени» со случайным изменением:
plot(time,y)
И положительные линейные отношения с «х» со случайным изменением:
plot(x,y)
Создайте простую линейную аддитивную модель "y ~ time + x":
require(nlme)
m1 <- gls(y ~ time + x, method="REML")
Модель отображает четкие циклические закономерности в остатках при построении в зависимости от «времени», как и следовало ожидать:
plot(time, m1$residuals)
И что должно быть приятным, ясным отсутствием какого-либо паттерна или тренда в остатках при построении против 'x':
plot(x, m1$residuals)
Простая модель "y ~ time + x", которая включает авторегрессивную корреляционную структуру порядка 1, должна соответствовать данным намного лучше, чем предыдущая модель, из-за автокорреляционной структуры, при оценке с использованием AIC:
m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")
AIC(m1,m2)
Однако модель должна по-прежнему отображать практически идентичные «временные» автокоррелированные остатки:
plot(time, m2$residuals)
Большое спасибо за любой совет.
Ответы:
Q1
Вы делаете две вещи неправильно здесь. Первое - вообще плохая вещь; вообще не копайтесь в модельных объектах и не вырывайте компоненты. Научитесь использовать функции экстрактора, в этом случае
resid()
. В этом случае вы получаете что-то полезное, но если у вас был другой тип объекта модели, такой как GLMglm()
, то онmod$residuals
будет содержать рабочие остатки из последней итерации IRLS и вам вообще не нужен!Второе, что вы делаете неправильно, это то, что меня тоже застало врасплох. Остатки, которые вы извлекли (а также извлекли бы, если бы использовали
resid()
), являются необработанными или ответными остатками. По сути, это разница между подобранными значениями и наблюдаемыми значениями отклика, принимая во внимание только условия с фиксированными эффектами . Эти значения будут содержать ту же остаточную автокорреляцию, что и из-m1
за того, что фиксированные эффекты (или, если хотите, линейный предиктор) одинаковы в двух моделях (~ time + x
).Чтобы получить остатки, которые включают указанный вами коэффициент корреляции, вам нужны нормализованные остатки. Вы получаете это, делая:
Это (и другие доступные типы остатков) описано в
?residuals.gls
:Для сравнения, вот ACF необработанных (ответ) и нормализованных остатков
Чтобы понять, почему это происходит, и где необработанные остатки не включают термин корреляции, рассмотрите модель, которую вы выбрали
где
Необработанные остатки, по умолчанию возвращаемые
resid(m2)
из только части линейного предиктора, следовательно, из этого битаQ2
Похоже, что вы пытаетесь согласовать нелинейный тренд с линейной функцией
time
и учитывать отсутствие соответствия «тренду» с AR (1) (или другими структурами). Если ваши данные похожи на пример данных, которые вы приводите здесь, я бы использовал GAM, чтобы обеспечить гладкую функцию ковариат. Эта модель будетгде
select = TRUE
применяется некоторая дополнительная усадка, чтобы позволить модели удалить любой из терминов из модели.Эта модель дает
и имеет гладкие условия, которые выглядят так:
Остатки от этой модели также лучше себя ведут (необработанные остатки)
Теперь слово предостережения; Существует проблема сглаживания временных рядов, заключающаяся в том, что методы, определяющие, насколько гладкими или волнистыми являются функции, предполагают независимость данных. На практике это означает, что сглаженная функция time (
s(time)
) может соответствовать информации, которая представляет собой действительно случайную автокоррелированную ошибку, а не только основной тренд. Следовательно, вы должны быть очень осторожны при подборе сглаживателей к данным временных рядов.Есть несколько способов обойти это, но один из них состоит в том, чтобы перейти к подгонке модели, с помощью
gamm()
которой осуществляетсяlme()
внутренний вызов и которая позволяет вам использоватьcorrelation
аргумент, который вы использовали дляgls()
модели. Вот примерs(time)
s(time)
s(time)
Модель с AR (1) не представляет значительного улучшения по сравнению с моделью без AR (1):
Если мы посмотрим на оценку для $ \ hat {\ rho}}, мы увидим
Phi
источник