Les modèles résiduels autocorrélés restent-ils même dans les modèles avec des structures de corrélation appropriées, et comment sélectionner les meilleurs modèles?

Le contexte

Cette question utilise R, mais concerne des problèmes statistiques généraux.

J'analyse les effets des facteurs de mortalité (% de mortalité due aux maladies et au parasitisme) sur le taux de croissance de la population de papillons au fil du temps, où les populations de larves ont été échantillonnées à partir de 12 sites une fois par an pendant 8 ans. Les données sur le taux de croissance démographique affichent une tendance cyclique claire mais irrégulière au fil du temps.

Les résidus d'un modèle linéaire généralisé simple (taux de croissance ~% maladie +% parasitisme + année) ont montré une tendance cyclique tout aussi claire mais irrégulière dans le temps. Par conséquent, des modèles de moindres carrés généralisés de la même forme ont également été ajustés aux données avec des structures de corrélation appropriées pour gérer l'autocorrélation temporelle, par exemple la symétrie composée, l'ordre de processus autorégressif 1 et les structures de corrélation moyenne mobile autorégressives.

Les modèles contenaient tous les mêmes effets fixes, ont été comparés en utilisant AIC et ont été ajustés par REML (pour permettre la comparaison de différentes structures de corrélation par AIC). J'utilise le package R nlme et la fonction gls.

question 1

Les valeurs résiduelles des modèles GLS affichent toujours des modèles cycliques presque identiques lorsqu'ils sont tracés en fonction du temps. De tels modèles resteront-ils toujours, même dans les modèles qui tiennent compte avec précision de la structure d'autocorrélation?

J'ai simulé des données simplifiées mais similaires dans R ci-dessous ma deuxième question, qui montre le problème basé sur ma compréhension actuelle des méthodes nécessaires pour évaluer les modèles autocorrélés temporellement dans les résidus du modèle , que je connais maintenant faux (voir réponse).

question 2

J'ai ajusté les modèles GLS avec toutes les structures de corrélation plausibles possibles à mes données, mais aucun n'est en réalité bien mieux adapté que le GLM sans aucune structure de corrélation: un seul modèle GLS est légèrement meilleur (score AIC = 1,8 inférieur), tandis que tous les autres ont des valeurs AIC plus élevées. Cependant, ce n'est le cas que lorsque tous les modèles sont équipés de REML, pas ML où les modèles GLS sont clairement beaucoup mieux, mais je comprends que dans les livres de statistiques, vous ne devez utiliser REML que pour comparer des modèles avec différentes structures de corrélation et les mêmes effets fixes pour des raisons Je ne détaillerai pas ici.

Étant donné la nature clairement autocorrélée temporellement des données, si aucun modèle n'est même modérément meilleur que le GLM simple, quelle est la façon la plus appropriée de décider quel modèle utiliser pour l'inférence, en supposant que j'utilise une méthode appropriée (je veux éventuellement utiliser AIC pour comparer différentes combinaisons de variables)?

«Simulation» du premier trimestre explorant les schémas résiduels dans les modèles avec et sans structures de corrélation appropriées

Générez une variable de réponse simulée avec un effet cyclique de «temps» et un effet linéaire positif de «x»:

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)

y devrait afficher une tendance cyclique dans le temps avec une variation aléatoire:

plot(time,y)

Et une relation linéaire positive avec 'x' avec une variation aléatoire:

plot(x,y)

Créez un modèle additif linéaire simple de "y ~ temps + x":

require(nlme)
m1 <- gls(y ~ time + x, method="REML")

Le modèle affiche des modèles cycliques clairs dans les résidus lorsqu'il est tracé en fonction du «temps», comme on pourrait s'y attendre:

plot(time, m1$residuals)

Et ce qui devrait être un manque clair et clair de tout motif ou tendance dans les résidus lorsqu'il est tracé contre `` x '':

plot(x, m1$residuals)

Un modèle simple de "y ~ temps + x" qui comprend une structure de corrélation autorégressive d'ordre 1 devrait mieux correspondre aux données que le modèle précédent en raison de la structure d'autocorrélation, lorsqu'il est évalué à l'aide de l'AIC:

m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")
AIC(m1,m2)

Cependant, le modèle devrait toujours afficher des résidus autocorrélés presque identiquement «temporellement»:

plot(time, m2$residuals)

Merci beaucoup pour tout conseil.

Q1

Vous faites deux mauvaises choses ici. Le premier est généralement une mauvaise chose; ne fouillez pas en général dans les objets du modèle et n'extrayez pas les composants. Apprenez à utiliser les fonctions d'extraction, dans ce cas resid(). Dans ce cas, vous obtenez quelque chose d'utile mais si vous aviez un autre type d'objet de modèle, tel qu'un GLM glm(), mod$residualscontiendrait alors des résidus de travail de la dernière itération IRLS et vous ne le souhaitez généralement pas !

La deuxième chose que vous faites mal, c'est quelque chose qui m'a aussi rattrapé. Les résidus que vous avez extraits (et que vous auriez également extraits si vous les aviez utilisés resid()) sont les résidus bruts ou de réponse. Essentiellement , cela est la différence entre les valeurs ajustées et les valeurs observées de la réponse, en tenant compte des effets fixes termes seulement . Ces valeurs contiendront la même autocorrélation résiduelle que celle de m1car les effets fixes (ou si vous préférez, le prédicteur linéaire) sont les mêmes dans les deux modèles ( ~ time + x).

Pour obtenir des résidus qui incluent le terme de corrélation que vous avez spécifié, vous avez besoin des résidus normalisés . Vous les obtenez en faisant:

resid(m1, type = "normalized")

Ceci (et d'autres types de résidus disponibles) est décrit dans ?residuals.gls:

type: an optional character string specifying the type of residuals
      to be used. If ‘"response"’, the "raw" residuals (observed -
      fitted) are used; else, if ‘"pearson"’, the standardized
      residuals (raw residuals divided by the corresponding
      standard errors) are used; else, if ‘"normalized"’, the
      normalized residuals (standardized residuals pre-multiplied
      by the inverse square-root factor of the estimated error
      correlation matrix) are used. Partial matching of arguments
      is used, so only the first character needs to be provided.
      Defaults to ‘"response"’.

A titre de comparaison, voici les ACF du brut (réponse) et des résidus normalisés

layout(matrix(1:2))
acf(resid(m2))
acf(resid(m2, type = "normalized"))
layout(1)

Pour voir pourquoi cela se produit et où les résidus bruts n'incluent pas le terme de corrélation, considérez le modèle que vous avez ajusté

où

et $\Lambda$$\hat{\rho}$$\rho^{|d|}$$d$

Les résidus bruts, les valeurs renvoyées par défaut resid(m2)proviennent uniquement de la partie du prédicteur linéaire, donc de ce bit

$\Lambda$

Q2

Il semble que vous essayez d'ajuster une tendance non linéaire avec une fonction linéaire timeet de tenir compte du manque d'ajustement à la "tendance" avec un AR (1) (ou d'autres structures). Si vos données sont quelque chose comme les données d'exemple que vous donnez ici, j'adapterais un GAM pour permettre une fonction lisse des covariables. Ce modèle serait

$\Lambda = \mathbf{I}$

library("mgcv")
m3 <- gam(y ~ s(time) + s(x), select = TRUE, method = "REML")

où select = TRUEapplique un retrait supplémentaire pour permettre au modèle de supprimer l'un ou l'autre des termes du modèle.

Ce modèle donne

> summary(m3)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(time) + s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  23.1532     0.7104   32.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
          edf Ref.df      F  p-value    
s(time) 8.041      9 26.364  < 2e-16 ***
s(x)    1.922      9  9.749 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

et a des termes fluides qui ressemblent à ceci:

Les résidus de ce modèle se comportent également mieux (résidus bruts)

acf(resid(m3))

Maintenant un mot d'avertissement; il y a un problème avec le lissage des séries temporelles dans la mesure où les méthodes qui déterminent le niveau de fluidité ou de ondulation des fonctions supposent que les données sont indépendantes. Concrètement, cela signifie que la fonction lisse de time ( s(time)) pourrait correspondre à des informations qui sont vraiment des erreurs autocorrélées aléatoires et pas seulement la tendance sous-jacente. Par conséquent, vous devez être très prudent lors de l'ajustement des lisseurs aux données de séries chronologiques.

Il existe plusieurs façons de contourner ce problème, mais l'une d'elles consiste à passer à l'ajustement du modèle via gamm()lequel les appels lme()internes et qui vous permettent d'utiliser l' correlationargument que vous avez utilisé pour le gls()modèle. Voici un exemple

mm1 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML")
mm2 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML", correlation = corAR1(form = ~ time))

Notez que je dois corriger les degrés de liberté s(time)car il y a un problème d'identification avec ces données. Le modèle pourrait être un ondulés(time)$\rho = 0$s(time)$\rho >> .5$

Le modèle avec l'AR (1) ne représente pas une amélioration significative par rapport au modèle sans l'AR (1):

> anova(mm1$lme, mm2$lme)
        Model df      AIC      BIC    logLik   Test   L.Ratio p-value
mm1$lme     1  9 301.5986 317.4494 -141.7993                         
mm2$lme     2 10 303.4168 321.0288 -141.7084 1 vs 2 0.1817652  0.6699

Si nous regardons l'estimation de $ \ hat {\ rho}}, nous voyons

> intervals(mm2$lme)
....

 Correlation structure:
         lower      est.     upper
Phi -0.2696671 0.0756494 0.4037265
attr(,"label")
[1] "Correlation structure:"

où Phiest ce que j'ai appelé$\rho$$\rho$