REPEATED statement de SAS et {nlme}

[Florian Heinze]

Bonsoir à tous,

Je dois transcrire un code SAS en un code R. Je travaille sur des données longitudinales et je veux spécifier la structure de covariance des résidus (en l'occurrence un processus auto-régressif d'ordre 1).

Le code SAS est le suivant :

Code : Tout sélectionner

Proc mixed data=df;
class id td;
model Y=t /solution;
repeated td/ subject=id type=AR(1) LOCAL;

l'option LOCAL permettant d'obtenir une erreur indépendante en plus de l'erreur autocorrélée.

Le package nlme est censé me permettre de spécifier ce type de structure (ce que ne permet pas lme4).
Je pensais réussir avec la fonction lme() mais je n'y suis pas parvenu. Après d'âpres recherches sur la toile, il se trouve que c'est gls() qu'il faut utiliser "Be sure to switch to gls when you remove the random effects." lu sur stackoverflow http://goo.gl/5FzCd

Code : Tout sélectionner

gls(Y~t, data=df, na.action = (na.omit), method = "REML",
    correlation=corAR1(form=~1|id))

Avec cette syntaxe j'obtiens la même estimation (que la syntaxe SAS sans l'option LOCAL) pour le phi (AR1), mais j'ai une différence au niveau de l'erreur.

Quelqu'un peut-il m'éclairer... J'ai l'impression de ne plus comprendre ce que je fais, et pire, je n'y arrive pas



Cordialement

[Renaud Lancelot]

Est-ce que SAS renvoie des estimateurs ML ou REML ? Il nous faudrait une jeu de données (et la sortie SAS correspondante) pour reproduire les résultats. Voir les post-it pour savoir comment faire.

[Stéphane Laurent]

SAS fait du REML par défaut.

Je pense que Florian devrait trouver les mêmes résultats, mais que les termes d'erreur sont présentés de façon différente. Je ne peux rien dire sans voir la sortie car je n'ai pas beaucoup l'habitude de gls(), mais je ne serais pas étonné par exemple qu'il retourne un écart-type alors que SAS retourne une variance.

[Florian Heinze]

Bonjour,

Malheureusement, je ne suis pas autorisé à mettre en ligne le jeu de données.

Sinon, oui SAS renvoie des REML.

Donc pour le code SAS suivant

Code : Tout sélectionner

Proc mixed data=df;
class id td;
model Y=t /s;
repeated td/ subject=id type=AR(1);

dont la sortie est

Code : Tout sélectionner

AR(1) 0.7246
Residual 0.1426

Voici ma syntaxe R

Code : Tout sélectionner

require(nlme)
gls(Y~t, data=df, na.action = (na.omit), method = "REML", correlation=corAR1(form=~as.factor(td)|id))

et la sortie correspondante

Code : Tout sélectionner

Generalized least squares fit by REML
  Model: Y ~ t
  Data: df
  Log-restricted-likelihood: -137

Coefficients:
(Intercept)           t
        4.5         0.1

Correlation Structure: AR(1)
 Formula: ~as.factor(td) | id
 Parameter estimate(s):
Phi
0.7
Degrees of freedom: 950 total; 948 residual
Residual standard error: 0.4

alors qu'avec le code SAS suivant avec option LOCAL :

Code : Tout sélectionner

Proc mixed data=df;
class id td;
model Y=t /s;
repeated td/ subject=id type=AR(1) LOCAL;

on obtient la sortie SAS :

Code : Tout sélectionner

Variance 0.1003
AR(1) 0.9755
Residual 0.03952

J'ignore si ça peut aider sans jeu de données, mais je ne peux renseigner que ça.

Merci.

[Stéphane Laurent]

Pourquoi tout est arrondi à 0.1 dans la sortie de gls() ?

Comme je le disais je ne serais pas étonné que SAS donne la variance mais gls() l'écart-type :

Code : Tout sélectionner

> sqrt(0.1426)
[1] 0.3776242

N'est-ce pas le 0.4 de gls() ?

[Florian Heinze]

Bonjour Stéphane,

Merci pour tes commentaires.

Effectivement, j'ai copié-collé une sortie avec option sur les digits. Sinon ça donne :

Code : Tout sélectionner

> gls(Y~t, data=df, na.action = (na.omit), method = "REML", correlation=corAR1(form=~as.factor(td)|id))
Generalized least squares fit by REML
  Model: Y ~ t
  Data: df
  Log-restricted-likelihood: -137

Coefficients:
(Intercept)           t
      4.519       0.124

Correlation Structure: AR(1)
 Formula: ~as.factor(td) | id
 Parameter estimate(s):
  Phi
0.712
Degrees of freedom: 950 total; 948 residual
Residual standard error: 0.375

Tu penses que c'est donc bon ? (même si je ne sais toujours pas comment faire avec LOCAL, savoir faire sans est déjà important !)

[Stéphane Laurent]

C'est fort proche. Pour te convaincre que SAS et R font bien la même chose tu devrais essayer avec d'autres données. Il devrait y avoir des cas où les résultats sont exactement les mêmes. Enfin je n'en sais rien, j'ai jamais essayé.
C'est quoi ce LOCAL ?

[Florian Heinze]

Ok ! Merci !

LOCAL je l'ai expliqué dans mon post de départ et ça apparaît sous l'intitulé "variance" dans la sortie correspondante.

[Florian Heinze]

Plus précisément selon l'aide SAS :

LOCAL
requests that a diagonal matrix be added to . With just the LOCAL option, this diagonal matrix equals , and becomes an additional variance parameter that PROC MIXED profiles out of the likelihood provided that you do not specify the NOPROFILE option in the PROC MIXED statement. The LOCAL option is useful if you want to add an observational error to a time series structure (Jones and Boadi-Boateng 1991) or a nugget effect to a spatial structure (Cressie 1993).

Comment implémenter ça sous R ? Je ne le sais pas.

[Jean-baptiste Woillard]

Dans SAS PROC MIXED, il y a des effets aléatoires ce que tu n'a plus lorsque tu utilise gls, or dans tes données longitudinales, tu veux peut être garder un effet aléatoire sur le temps. Je pense qu'il vaut mieux que tu utilise lme avec correlation=corAR1.
jb

[Florian Heinze]

Bonjour,

En fait je traduis un code SAS en un code R. Il s'agit d'un cours où tous les cas de figure sont envisagés. Mais bien sûr, au final on montre que l'on garde le repeated statement en plus du random statement avec lme().

Merci en tout cas,

Bonne fin de journée.

[Stéphane Laurent]

Dans SAS PROC MIXED, il y a des effets aléatoires

Non pas ici, il n'y a pas de ligne RANDOM.

[Jean-baptiste Woillard]

C'est quand même étrange que dans une procédure effets mixtes on puisse ne mettre que des effets fixes!!!!!!
Je suis bien content de ne pas travailler dans SAS!!!

[Florian Heinze]

oui, SAS, c'est vraiment tout un poème :-)

[Stéphane Laurent]

Au fait aujourd'hui j'ai vu un truc sur le ouebbe : asreml http://www.r-project.org/conferences/useR-2007/program/presentations/butler.pdf

Apparemment c'est un package R qui permettrait d'ajuster des modèles mixtes avec plus de possibilités que lme() et lmer(). Mais c'est payant !