Comment effectuer une ANCOVA en R

17

Je souhaite effectuer une analyse ANCOVA des données concernant la densité des épiphytes végétaux. Dans un premier temps, je voudrais savoir s'il existe une différence de densité de plante entre deux pentes, une N et une S, mais j'ai d'autres données telles que l'altitude, l'ouverture du couvert et la hauteur de la plante hôte. Je sais que ma covariable devrait être les deux pentes (N et S). J'ai construit ce modèle qui fonctionne en R et bien que je ne sache pas s'il fonctionne bien. J'aimerais aussi savoir quelle est la différence si j'utilise le symbole +ou *.

model1 <- aov(density~slope+altitude+canopy+height)
summary(model1)
model1
r ancova Pauloc
la source
3
+ calculera uniquement les effets principaux, * estimera les interactions entre les facteurs liés à *. Les cadres ANCOVA n'évaluent généralement qu'un effet principal du facteur continu, mais les interactions entre tous les facteurs groupés.
russellpierce

Réponses:

14

L'outil de base pour cela est lm; notez que aovc'est un wrapper pour lm.

En particulier, si vous avez une variable de regroupement (facteur), g et une covariable continueX , le modèley ~ x + gcorrespondrait à un modèle ANCOVA à effets principaux, tandisy ~ x * gqu'il s'adapterait à un modèle qui inclut l'interaction avec la covariable. aovprendra les mêmes formules.

Portez une attention particulière Noteà l'aide dans aov.

Pour ce qui est de + vs *, russellpierce à peu près couvertures, mais je vous recommande de regarder ?lmet ?formulaet plus particulièrement l' article 11.1 du manuel Introduction à R qui vient avec R (ou vous pouvez le trouver en ligne si vous ne l' avez pas compris comment pour le trouver sur votre ordinateur, le plus facilement, cela implique de trouver le menu déroulant "Aide" dans R ou RStudio).

Glen_b -Reinstate Monica
la source
supposons que j'ai deux facteurs de groupe et deux covariables x 1 , x 2 , mon modèle étant y i j = μ + α i + η j + x i j 1 γ 1 + x i j 2 γ 2 + ϵ i j Est-ce que y ~ g_1 + g_2 + x_1 + x_2 fait le même tour? Est-ce que les valeurs F obtenues contre x_1 et x_2 testent γ 1 = 0 et γg1,g2X1,X2
yjej=μ+αje+ηj+Xjej1γ1+Xjej2γ2+ϵjej
γ1=0 respectivement? γ2=0
Sayan
Je ne sais pas comment j'ai raté ça. Oui. .... et si vous voulez tester les deux à la fois, ajustez-les avec et sans eux et passez les objets lm ajustés à anova(vous verrez bientôt si vous les donnez dans le mauvais ordre car certains SS seront négatifs si vous le faites) )
Glen_b -Reinstate Monica
10

Je recommande d'obtenir et de lire Découvrir les statistiques à l'aide de R par champ. Il a une belle section sur ANCOVA.

Pour exécuter ANCOVA dans R, chargez les packages suivants: 

car
compute.es
effects
ggplot2
multcomp
pastecs
WRS

Si vous utilisez lmou aov(j'utilise aov) assurez-vous de régler les contrastes à l'aide de la fonction "contrastes" avant de faire soit aovoulm . R utilise par défaut des contrastes non orthogonaux qui peuvent tout gâcher dans une ANCOVA. Si vous souhaitez définir des contrastes orthogonaux, utilisez:

contrasts(dataname$factorvariable)=contr.poly(# of levels, i.e. 3)

puis exécutez votre modèle en tant que 

model.1=aov(dv~covariate+factorvariable, data=dataname)

Pour afficher l'utilisation du modèle:

Anova(model.1, type="III")

Assurez-vous d'utiliser le "A" majuscule Anova ici et non anova. Cela donnera des résultats en utilisant le SS de type III.

summary.lm(model.1)donnera un autre résumé et comprend le R-sq. production. 

posth=glht(model.1, linfct=mcp(factorvariable="Tukey"))  ##gives the post-hoc Tukey analysis
summary(posth) ##shows the output in a nice format.

Si vous souhaitez tester l'homogénéité des pentes de régression, vous pouvez également inclure un terme d'interaction pour IV et covariable. Ce serait: 

model=aov(dv~covariate+IV+covariate:IV, data=dataname)

Si le terme d'interaction est significatif, vous n'avez pas d'homogénéité.

Dan Butorovich
la source
Pourquoi les contrastes non orthogonaux gâchent-ils tout?
Tintinthong
1
Pour répondre à la question ci-dessus sur "pourquoi les contrastes non orthogonaux gâchent tout". La réponse est que R par défaut est non orthogonal (c'est-à-dire la différence entre les moyennes), ce qui peut causer des problèmes si vous voulez voir la contribution de chaque IV séparément. Lorsque nous spécifions des contrastes orthogonaux, nous disons à R que nous voulons que les SS pour les IV soient entièrement partitionnés et ne se chevauchent pas. De cette façon, nous pouvons voir clairement et clairement la variation attribuée à chaque prédicteur. Si vous ne spécifiez pas, R par défaut à une approche plus libérale du contraste.
2
Pourquoi l'intérêt pour le SS de type III?
Frank Harrell
4

Voici une documentation complémentaire http://goo.gl/yxUZ1R de la procédure proposée par @Butorovich. De plus, mon observation est que lorsque la covariable est binaire, l'utilisation de résumé (lm.object) donnerait la même estimation IV que celle générée par Anova (lm.object, type = "III").

XX
la source
1
Il n'est pas clair que ce soit censé être une réponse. C'est ça? Si oui, veuillez modifier pour clarifier. S'il s'agit d'une question, posez-la en cliquant ASK QUESTIONsur le haut et posez-la là. Ensuite, nous pouvons vous aider correctement.
gung - Reinstate Monica
D'accord. Le message a été révisé en tant que réponse (complémentaire) à la précédente.
XX
3

Nous utilisons l'analyse de régression pour créer des modèles qui décrivent l'effet de la variation des variables prédictives sur la variable de réponse. Parfois, si nous avons une variable catégorielle avec des valeurs comme Oui / Non ou Masculin / Féminin, etc., l'analyse de régression simple donne plusieurs résultats pour chaque valeur de la variable catégorielle. Dans un tel scénario, nous pouvons étudier l'effet de la variable catégorielle en l'utilisant avec la variable prédictive et en comparant les lignes de régression pour chaque niveau de la variable catégorielle. Une telle analyse est appelée analyse de covariance également appelée ANCOVA.

Exemple
Considérez l' Rensemble de données intégré mtcars. Nous y observons que le champ amreprésente le type de transmission (automatique ou manuelle). Il s'agit d'une variable catégorielle avec les valeurs 0 et 1. La valeur en miles par gallon ( mpg) d'une voiture peut également en dépendre en plus de la valeur de la puissance en chevaux ( hp). Nous étudions l'effet de la valeur de amsur la régression entre mpget hp. Cela se fait en utilisant la aov()fonction suivie de la anova()fonction pour comparer les régressions multiples.

Les données d'entrée
Créer une trame de données contenant les champs mpg, hpet amde l'ensemble de données mtcars. Ici, nous prenons mpgcomme variable de réponse, hpcomme variable prédictive et amcomme variable catégorielle.

input <- mtcars[,c("am","mpg","hp")]
head(input)

Lorsque nous exécutons le code ci-dessus, cela produit le résultat suivant:

                  am  mpg  hp
Mazda RX4          1 21.0 110
Mazda RX4 Wag      1 21.0 110
Datsun 710         1 22.8  93
Hornet 4 Drive     0 21.4 110
Hornet Sportabout  0 18.7 175
Valiant            0 18.1 105

Analyse ANCOVA
Nous créons un modèle de régression en prenant hpcomme variable prédictive et mpgcomme variable de réponse en tenant compte de l'interaction entre amet hp.

Modèle avec interaction entre variable catégorielle et variable prédictive

Créer un modèle de régression1 

result1 <- aov(mpg~hp*am,data=mtcars)
summary(result1)

Lorsque nous exécutons le code ci-dessus, cela produit le résultat suivant:

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4  77.391 1.50e-09 ***
am           1  202.2   202.2  23.072 4.75e-05 ***
hp:am        1    0.0     0.0   0.001    0.981    
Residuals   28  245.4     8.8                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ce résultat montre que la puissance en chevaux et le type de transmission ont un effet significatif sur les miles par gallon car la valeur p dans les deux cas est inférieure à 0,05. Mais l'interaction entre ces deux variables n'est pas significative car la valeur de p est supérieure à 0,05.

Modèle sans interaction entre variable catégorielle et variable prédictive

Créer le modèle de régression2 

result2 <- aov(mpg~hp+am,data=mtcars)
summary(result2)

Lorsque nous exécutons le code ci-dessus, cela produit le résultat suivant:

            Df Sum Sq Mean Sq F value   Pr(>F)    
hp           1  678.4   678.4   80.15 7.63e-10 ***
am           1  202.2   202.2   23.89 3.46e-05 ***
Residuals   29  245.4     8.5                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ce résultat montre que la puissance en chevaux et le type de transmission ont un effet significatif sur les miles par gallon car la valeur p dans les deux cas est inférieure à 0,05.

Comparaison de deux modèles
Maintenant, nous pouvons comparer les deux modèles pour conclure si l'interaction des variables est vraiment insignifiante. Pour cela, nous utilisons la anova()fonction.

 anova(result1,result2)

 Model 1: mpg ~ hp * am
 Model 2: mpg ~ hp + am
   Res.Df    RSS Df  Sum of Sq     F Pr(>F)
 1     28 245.43                           
 2     29 245.44 -1 -0.0052515 6e-04 0.9806

Comme la valeur de p est supérieure à 0,05, nous concluons que l'interaction entre la puissance en chevaux et le type de transmission n'est pas significative. Ainsi, le kilométrage par gallon dépendra de manière similaire de la puissance en chevaux de la voiture en mode de transmission automatique et manuelle.

user101834
la source
Alors, qui est venu en premier, cette réponse ou ce post sur le point tutoriels? tutorialspoint.com/r/r_analysis_of_covariance.htm . Faut-il supprimer cette réponse comme plagiat? Ou les didacticiels pointaient-ils simplement de copier à partir d'ici?
John