Modèle linéaire simple avec erreurs autocorrélées dans R [fermé]

19

Comment ajuster un modèle linéaire avec des erreurs autocorrélées dans R? En stata, j'utiliserais la praiscommande, mais je ne trouve pas d'équivalent R ...

Zach
la source

Réponses:

23

Jetez un oeil à gls(moindres carrés généralisés) du paquet nlme

Vous pouvez définir un profil de corrélation pour les erreurs dans la régression, par exemple ARMA, etc.:

 gls(Y ~ X, correlation=corARMA(p=1,q=1))

pour les erreurs ARMA (1,1).

Dr G
la source
1
Puis-je utiliser la fonction "prédire" sur un nouvel ensemble de données et conserver la même structure d'erreur? Comment la commande gls sait-elle dans quel ordre sont les observations?
Zach
27

En plus de la gls()fonction de nlme, vous pouvez également utiliser la arima()fonction dans le statspackage à l'aide de MLE. Voici un exemple avec les deux fonctions.

x <- 1:100
e <- 25*arima.sim(model=list(ar=0.3),n=100)
y <- 1 + 2*x + e

###Fit the model using gls()
require(nlme)
(fit1 <- gls(y~x, corr=corAR1(0.5,form=~1)))
Generalized least squares fit by REML
  Model: y ~ x 
  Data: NULL 
  Log-restricted-likelihood: -443.6371

Coefficients:
(Intercept)           x 
   4.379304    1.957357 

Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
      Phi 
0.3637263 
Degrees of freedom: 100 total; 98 residual
Residual standard error: 22.32908 

###Fit the model using arima()
(fit2 <- arima(y, xreg=x, order=c(1,0,0)))

Call:
arima(x = y, order = c(1, 0, 0), xreg = x)

Coefficients:
         ar1  intercept       x
      0.3352     4.5052  1.9548
s.e.  0.0960     6.1743  0.1060

sigma^2 estimated as 423.7:  log likelihood = -444.4,  aic = 896.81 

L'avantage de la fonction arima () est que vous pouvez adapter une plus grande variété de processus d'erreur ARMA. Si vous utilisez la fonction auto.arima () du package de prévisions, vous pouvez identifier automatiquement l'erreur ARMA:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)
Rob Hyndman
la source
1
À quoi sert le 0.5 dans "corr = corAR1 (0.5, form = ~ 1)?"
Zach
1
Il donne une valeur de départ pour l'optimisation. Cela ne fait presque aucune différence s'il est omis.
Rob Hyndman
1
+1 L' arimaoption semble plus différente de celle de Stata praisà première vue, mais elle est plus flexible et vous pouvez également l'utiliser tsdiagpour obtenir un bon aperçu de l'adéquation de votre hypothèse AR (1).
Wayne
1
Comment utiliser arima () si je régresse y uniquement sur une constante?
user43790
7

Utilisez la fonction gls du paquet nlme . Voici l'exemple.

##Generate data frame with regressor and AR(1) error. The error term is 
## \eps_t=0.3*\eps_{t-1}+v_t
df <- data.frame(x1=rnorm(100), err=filter(rnorm(100)/5,filter=0.3,method="recursive"))

##Create ther response
df$y <- 1 + 2*df$x + df$err

###Fit the model
gls(y~x, data=df, corr=corAR1(0.5,form=~1))

Generalized least squares fit by REML
  Model: y ~ x 
  Data: df 
  Log-restricted-likelihood: 9.986475

 Coefficients:
 (Intercept)           x 
   1.040129    2.001884 

 Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
     Phi 
 0.2686271 
Degrees of freedom: 100 total; 98 residual
Residual standard error: 0.2172698 

Étant donné que le modèle est ajusté en utilisant la probabilité maximale, vous devez fournir des valeurs de départ. La valeur de départ par défaut est 0, mais comme toujours, il est bon d'essayer plusieurs valeurs pour assurer la convergence.

Comme l'a souligné le Dr G, vous pouvez également utiliser d'autres structures de corrélation, à savoir ARMA.

Notez qu'en général, les estimations des moindres carrés sont cohérentes si la matrice de covariance des erreurs de régression n'est pas multiple de la matrice d'identité, donc si vous ajustez le modèle avec une structure de covariance spécifique, vous devez d'abord tester si elle est appropriée.

mpiktas
la source
À quoi sert le 0.5 dans "corr = corAR1 (0.5, form = ~ 1)?"
Zach
3

Vous pouvez utiliser la prévision sur la sortie gls. Voir? Predire.gls. Vous pouvez également spécifier l'ordre de l'observation par le terme "forme" dans la structure de corrélation. Par exemple:
corr=corAR1(form=~1)indique que l'ordre des données est celui qu'elles sont dans le tableau. corr=corAR1(form=~Year)indique que l'ordre est celui du facteur Année. Enfin la valeur "0,5" en corr=corAR1(0.5,form=~1)?est généralement fixée à la valeur du paramètre estimé pour représenter la structure de la variance (phi, en cas d'AR, thêta en cas de MA .. .). Il est facultatif de le configurer et de l'utiliser pour l'optimisation comme l'a mentionné Rob Hyndman.

Xochitl C.
la source
Les valeurs prédites ou ajustées seront bien différentes mais correctes (gls () versus Arima ())? Comme gls n'utilisera que des effets fixes tandis qu'Arima inclura les erreurs d'arima aux valeurs ajustées.
B_Miner