Cómo realizar un ANCOVA en R

Quiero realizar un análisis ANCOVA de datos relativos a la densidad de epífitas de plantas. Al principio, me gustaría saber si hay alguna diferencia en la densidad de la planta entre dos pendientes, una N y una S, pero tengo otros datos como la altitud, la apertura del dosel y la altura de la planta huésped. Sé que mi covariable tendría que ser las dos pendientes (N y S). Construí este modelo que se ejecuta en R y aunque no tengo idea si funciona bien. También me gustaría saber cuál es la diferencia si uso el símbolo +o *.

model1 <- aov(density~slope+altitude+canopy+height)
summary(model1)
model1

La herramienta básica para esto es lm; tenga en cuenta que aoves un contenedor para lm.

En particular, si tiene alguna variable de agrupación (factor), $g$ , y una covariable continua $x$ , el modelo y ~ x + gse ajustaría a un modelo ANCOVA de efectos principales, mientras que y ~ x * gse ajustaría a un modelo que incluye la interacción con la covariable. aovtomará las mismas fórmulas.

Presta especial atención a Notela ayuda en aov.

En cuanto a +frente *, russellpierce más o menos lo cubre, pero recomiendo nos fijamos en ?lmy ?formula, y muy especialmente la sección 11.1 del manual Introducción a R que viene con R (o se puede encontrar en línea si no se ha descubierto la manera para encontrarlo en su computadora; más fácilmente, esto implica encontrar el menú desplegable "Ayuda" en R o RStudio).

Recomiendo obtener y leer Discovering Statistics usando R by Field. Tiene una buena sección sobre ANCOVA.

Para ejecutar ANCOVA en R, cargue los siguientes paquetes:

car
compute.es
effects
ggplot2
multcomp
pastecs
WRS

Si está utilizando lmo aov(yo uso aov) asegúrese de establecer los contrastes utilizando la función "contrastes" antes de hacer cualquiera aovo lm. R usa contrastes no ortogonales por defecto que pueden estropear todo en un ANCOVA. Si desea establecer contrastes ortogonales, use:

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

luego ejecuta tu modelo como

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

Para ver el uso del modelo:

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

Asegúrese de usar mayúscula "A" Anovaaquí y no anova. Esto dará resultados con el tipo III SS.

summary.lm(model.1)dará otro resumen e incluye el R-sq. salida.

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

Si desea probar la homogeneidad de las pendientes de regresión, también puede incluir un término de interacción para IV y covariable. Eso sería:

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

Si el término de interacción es significativo, entonces no tiene homogeneidad.

Aquí hay una documentación complementaria http://goo.gl/yxUZ1R del procedimiento sugerido por @Butorovich. Además, mi observación es que cuando la covariable es binaria, el uso de resumen (lm.object) daría la misma estimación IV que la generada por Anova (lm.object, type = "III").

Utilizamos el análisis de regresión para crear modelos que describen el efecto de la variación en las variables predictoras sobre la variable de respuesta. A veces, si tenemos una variable categórica con valores como Sí / No o Masculino / Femenino, etc., el análisis de regresión simple da múltiples resultados para cada valor de la variable categórica. En tal escenario, podemos estudiar el efecto de la variable categórica usándola junto con la variable predictora y comparando las líneas de regresión para cada nivel de la variable categórica. Tal análisis se denomina Análisis de Covarianza también llamado ANCOVA.

Ejemplo
Considere el Rconjunto de datos incorporado mtcars. En él observamos que el campo amrepresenta el tipo de transmisión (automática o manual). Es una variable categórica con valores 0 y 1. El valor de millas por galón ( mpg) de un automóvil también puede depender de él además del valor de la potencia en caballos (hp ). Estudiamos el efecto del valor deam en la regresión entre mpgy hp. Se realiza utilizando la aov()función seguida de la anova()función para comparar las regresiones múltiples.

La entrada de datos
Crear una trama de datos que contiene los campos mpg, hpy amdel conjunto de datosmtcars . Aquí tomamosmpg como variable de respuesta, hpcomo variable predictiva y amcomo variable categórica.

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

Cuando ejecutamos el código anterior, produce el siguiente resultado:

                  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

Análisis ANCOVA
Creamos un modelo de regresión tomandohp como variable predictiva ympg como variable de respuesta teniendo en cuenta la interacción entre amy hp.

Modelo con interacción entre variable categórica y variable predictiva

Crear modelo de regresión1

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

Cuando ejecutamos el código anterior, produce el siguiente resultado:

            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

Este resultado muestra que tanto la potencia del caballo como el tipo de transmisión tienen un efecto significativo en millas por galón, ya que el valor p en ambos casos es menor a 0.05. Pero la interacción entre estas dos variables no es significativa ya que el valor p es más de 0.05.

Modelo sin interacción entre variable categórica y variable predictiva

Crea el modelo de regresión2

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

Cuando ejecutamos el código anterior, produce el siguiente resultado:

            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

Este resultado muestra que tanto la potencia del caballo como el tipo de transmisión tienen un efecto significativo en millas por galón, ya que el valor p en ambos casos es menor a 0.05.

Comparación de dos modelos
Ahora podemos comparar los dos modelos para concluir si la interacción de las variables es realmente insignificante. Para esto usamos la anova()función.

 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

Como el valor p es mayor que 0.05, concluimos que la interacción entre la potencia del caballo y el tipo de transmisión no es significativa. Por lo tanto, el kilometraje por galón dependerá de manera similar de la potencia del automóvil en el modo de transmisión automática y manual.