Discussion :

[CupOfAppleTea]

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 :

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

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 :

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

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 :

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

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 :

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

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 :

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

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!!

[CupOfAppleTea]

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 :

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

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 :

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

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

Voici ce que j'obtiens:

Code :

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

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 :

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

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!