Discussion :

[CupOfAppleTea]

Bonjour à tous,

Je dispose d'un jeu de données (cf. pièce jointe) avec une variable dépendante (Res) et deux variables indépendantes (ModeF et T; avec T la variable décrivant les 2 temps auxquels les données ont été recueillies). Je souhaite procéder à la régression linéaire suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

m1=lm(Res~ModeF*T, a2)

Cependant, les résidus du modèle présentent une hétéroscédasticité qui croit en fonction de ModeF. Plus ModeF augmente, plus l'hétéroscédasticité de mes résidus est importante. Pour résoudre ce problème, je souhaite procéder à une régression par les moindres carrés pondérés via la fonction "gls" du package "nlme" de R. Cette fonction présente un argument "weights" qui permet d'allouer moins d'importance aux points pour lesquels la variabilité est forte et inversement. Cependant, je ne sais pas comment renseigner l'argument "weights"!

Dans l'aide de R, on peut trouver un certain nombre de fonction de variance (appelée "VarClasses") tel que varExp que j'ai tenté d'utiliser:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

m3=gls(Res~ModeF*T, a2, weights=varExp(a2$ModeF))

Mais j'obtiens le message d'erreur suivant:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
> m3=gls(Res~ModeF*T, a2, weights=varExp(a2$ModeF))
Error in Initialize.varExp(X[[1L]], ...) : 
  initial value for "varExp" should be of length 1

J'ai également tenté comme suit afin que le poids soit de moins en moins important au fur et à mesure que ModeF augmente:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

m3=gls(Res~ModeF*T, a2, weights=(1/a2$ModeF))

Mais j'obtiens cette fois-ci l'erreur suivante:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
> m3=gls(Res~ModeF*T, a2, weights=(1/a2$ModeF^2))
Error in varFunc(weights) : 
  can only construct "varFunc" object from another "varFunc" object, a formula, or a character string

Merci pour l'aide que vous pourrez m'apporter!

[CupOfAppleTea]

J'ai pu trouver des informations dans le livre de Zuur et al. Mixed Effects Models and Extensions in Ecology with R :
http://books.google.fr/books?id=vQUN...eights&f=false

J'ai donc procéder comme suit avec un premier modèle sans pondération (l'argument "weights" n'est pas renseigné), ce qui équivaut à une simple régression linéaire, et un second modèle avec la fonction de variance "varFixed" qui autorise une variabilité croissante des résidus en fonction de la variable ModeF:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
mLM=gls(Res~ModeF*T, a2)
mGLS=gls(Res~ModeF*T, a2, weights=varFixed(~ModeF))

Je procède ensuite à une ANOVA pour déterminer lequel de ces 2 modèles est le "meilleur":

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
anova(mLM,mGLS)
 
Model df      AIC      BIC    logLik
mLM      1  5 180.8594 191.9569 -85.42969
mGLS     2  5 179.7730 190.8706 -84.88650

Je peux donc conclure que le modèle mGLS est le meilleur des 2 (AIC plus faible) mais comment puis-je vérifier qu'il prend bien en compte l'hétéroscédasticité de mes données? Le test de Bartlett me renvoie la même valeurs pour ces 2 modèles:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
Bartlett test of homogeneity of variances
 
data:  residuals(mGLS) and interaction(a2$ModeF, a2$T, drop = TRUE)
Bartlett's K-squared = 35.2706, df = 11, p-value = 0.0002236
--------------------------------------------------------------
Bartlett test of homogeneity of variances
 
data:  residuals(mLM) and interaction(a2$ModeF, a2$T, drop = TRUE)
Bartlett's K-squared = 35.2706, df = 11, p-value = 0.0002236

Si je comprends bien, pondérer via l'argument "weights" permet de calculer l'erreur standard en considérant l'hétéroscédasticité et donc d'aboutir à des p-values interprétables, mais cela ne change rien aux résidus du modèles qui présentent toujours la même hétéroscédasticité. Il me manque une pièce du puzzle et j'aimerais comprendre ce que je fais!
Merci pour l'aide que vous pourrez m'apporter!