LM / Unbalanced data / Interactions

Bonjour,
Je dispose d'un jeu de données de température de feuille, avec pour chaque feuille l'indication East/West selon qu'elle se trouve du coté est ou ouest de l'arbre. En même temps a été mesurée la température de l'air, données dont je dispose également. J'ai construit le modèle linéaire suivant:

Code:

---

1 2	a=read.table("StatsDescript.txt", header=TRUE) mod=lm(Tleaf~Tair*orientation, na.action="na.exclude", data=a)

---

Cependant, comme mon jeu de données n'est pas équilibré (unbalanced data), l'ordre dans lequel j'ajoute mes variables au modèles change le résultat.

Voici pour le premier modèle:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13

> a=read.table("StatsDescript.txt", header=TRUE) > mod1=lm(Tleaf~Tair*orientation, na.action="na.exclude", data=a) > anova(mod1) Analysis of Variance Table   Response: Tleaf * * * * * * * * * * * * * * * * Df Sum Sq Mean Sq * F value * * * Pr(>F) * * Tair * * * * * * * * * * * * * *1 * 5419.9 *5419.9 * *2371.4658 < 2e-16 *** orientation * * * * * * * * *1 * *6.3 * * * *6.3 * * * *2.7769 * * * 0.09653 . * Tair:orientation * * * * * 1 * *11.1 * * *11.1 * * *4.8493 * * * *0.02831 * * Residuals * * * * * * * * * * * * *349 * * *797.6 * * *2.3 * * * * * * * * * * * --- Signif. codes: *0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

Voici pour le second modèle:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12

> mod2=lm(Tleaf~orientation*Tair, na.action="na.exclude", data=a) > anova(mod2) Analysis of Variance Table   Response: Tleaf * * * * * * * * * * * * * * * * Df Sum Sq Mean Sq * F value * * Pr(>F) * * orientation * * * * * * * * *1 * 12.0 * *12.0 * * * *5.2509 * * *0.02253 * * Tair * * * * * * * * * * * * * *1 5414.3 *5414.3 * 2368.9918 *< 2e-16 *** orientation:Tair * * * * * 1 * 11.1 * *11.1 * * *4.8493 * * * *0.02831 * * Residuals * * * * * * * * * * * * 349 * * 797.6 * * 2.3 * * * * * * * * * * * --- Signif. codes: *0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

Comme vous pouvez le voir, la variable "orientation" (East Vs. West) n'est pas significative dans le 1er modèle mais elle l'est dans le second.
D'après ce que j'ai lu (notamment sur ce site http://goanna.cs.rmit.edu.au/~fscholer/anova.php) il me faut calculer la somme des carrés (SS) par la méthode de type III du fait du jeu de données déséquilibré et de l'interaction entre mes deux variables explicatives.
J'ai donc récupéré le script suivant:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

> a=read.table("StatsDescript.txt", header=TRUE) > options(contrasts = c("contr.sum","contr.poly")) > mod=lm(Tleaf~Tair*orientation, na.action="na.exclude", data=a) > drop1(mod, .~., test="F") Single term deletions   Model: Tleaf ~ Tair * orientation * * * * * * * * * * * * * * * *Df Sum of Sq * *RSS * * AIC * * * F value * * * Pr(>F) * * <none> * * * * * * * * * * * * * * * * * * * * *797.6 * 295.76 * * * * * * * * * * * Tair * * * * * * * * * * 1 * *5158.6 * *5956.3 1003.48 2257.1256 *< 2e-16 *** orientation * * * * * * * *1 * * * 14.1 * * * 811.7 *299.94 * *6.1672 * * 0.01348 * * Tair:orientation * * * * *1 * * *11.1 * * * 808.7 * 298.63 * *4.8493 * * 0.02831 * * --- Signif. codes: *0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

La fonction drop() permet de déterminer le caractère significatif ou non de chacune des variables de manière indépendante à l'ordre dans lequel elles sont spécifiées dans le modèle.
J'aimerais maintenant pouvoir récupérer les estimateurs de chacun de mes paramètres. Est-ce une erreur si je procède ensuite à un summary:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

> summary(mod)   Call: lm(formula = Tleaf ~ Tair * orientation, * *na.action = "na.exclude")   Residuals: * *Min * * *1Q *Median * * *3Q * * Max -2.9928 -1.1355 -0.0953 *0.9155 *5.2725   Coefficients: * * * * * * * * * * * * * * * * Estimate Std. Error * * * *t value * * Pr(>|t|) * * (Intercept) * * * * * * * * * * * 4.77483 * *0.42897 *11.131 * <2e-16 *** Tair * * * * * * * * * * * * * * * * 0.89850 * *0.01891 *47.509 * <2e-16 *** orientation1 * * * * * * * * * * 1.06529 * *0.42897 * 2.483 * 0.0135 * * Tair:orientation1 * * * * * * -0.04165 * *0.01891 *-2.202 * 0.0283 * * --- Signif. codes: *0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 1.512 on 349 degrees of freedom *(2 observations deleted due to missingness) Multiple R-squared: *0.8721, Adjusted R-squared: *0.871 F-statistic: * 793 on 3 and 349 DF, *p-value: < 2.2e-16

---

summary() me renvoie la valeur de l'estimateur pour chacun des paramètres. Cependant, cette fonction traite les facteurs par niveau (cf. orientation dont on a l'effet que d'un niveau, East ou West je ne sais pas pour ce cas). Et est-ce correct de faire un summary après une anova sur le modèle.
Bref, j'ai peur de mélanger plusieurs choses et j'aimerais avoir votre avis. Comment obtenir la valeur des estimateurs autrement? Et sans considérer les variables factorielles par niveau?
Merci!!

Un autre moyen d'exposer mon problème:

J'ai besoin d'analyser un jeu de données non-équilibré via une régression linéaire:

Code:

---

modJuin=lm(TleafMax~TairMax*orientation, na.action="na.exclude", data=aJuin)

---

"TairMax" est une variable continue numérique et "orientation" est une variable factorielle à 2 niveaux. J'ai donc utilisé drop1() pour tester le modèle et gérer le caractère non-équilibré de mon jeu de données:

Code:

---

drop1(modJuin, .~., test="F")

---

Voici ce que j'obtiens:

Code:

---

1 2 3 4 5 6 7 8 9 10 11

Single term deletions   Model: TleafMax ~ TairMax * orientation Df Sum of Sq RSS AIC F value Pr(>F) <none> 797.6 295.76 TairMax 1 5158.6 5956.3 1003.48 2257.1256 < 2e-16 *** orientation 1 14.1 811.7 299.94 6.1672 0.01348 * TairMax:orientation 1 11.1 808.7 298.63 4.8493 0.02831 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

---

Maintenant, je sais quelles variables mènent à des modèles qui diffèrent significativement lorsqu'on les retire. Apparemment chacunes d'elles (TairMax, orientation) ainsi que leur interaction (TairMax*orientation). J'aimerais ensuite obtenir les estimateurs des paramètres du modèle. J'ai donc utilisé summary() comme ceci:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

> summary.lm(modJuin)   Call: lm(formula = TleafMax ~ TairMax * orientation, data = aJuin, na.action = "na.exclude")   Residuals: Min 1Q Median 3Q Max -2.9928 -1.1355 -0.0953 0.9155 5.2725   Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.77483 0.42897 11.131 <2e-16 *** TairMax 0.89850 0.01891 47.509 <2e-16 *** orientation1 1.06529 0.42897 2.483 0.0135 * TairMax:orientation1 -0.04165 0.01891 -2.202 0.0283 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   Residual standard error: 1.512 on 349 degrees of freedom (2 observations deleted due to missingness) Multiple R-squared: 0.8721, Adjusted R-squared: 0.871 F-statistic: 793 on 3 and 349 DF, p-value: < 2.2e-16

---

Voici ma question: Est-il correct d'obtenir les estimateurs des paramètres de mon modèle linéaire via summary() après avoir identifié les variables qui modifient significativement mon modèle via drop1()?
Egalement, je n'utilise pas summary() pour déterminer la significativité de mes variables car un collègue m'a fait une critique l'autre jour: "comme tu peux le voir dans l'output de summary(), tes facteurs sont traités par niveaux. En effet, l'un des niveaux de ton facteur n'apparait pas dans l'output et la significativité n'est donnée que pour les autres niveaux". A-t'il raison? Dans mon cas, cela s'illustre sur la variable orientation dont seul le niveaux "1" apparait dans le summary() (voir ci-dessus). S'il a tord, j'aurais tendance à utiliser directement summary() à la place de drop1() puisque summary() est indépendant de l'ordre de spécification des variables dans le modèle linéaire.

Je précise que je suis au courant des problèmes liés au fait d'identifier les effets individuels des variables lorsqu'un interaction existe. Merci pour votre aide!