Autor: Álvaro Alonso Fernández
Departamento de Ciencias de la Vida
Universidad de Alcalá (España)


Un buen café siempre puede ayudar a manejar R…..


1 Introducción

Vamos a realizar un análisis ANCOVA, que es un análisis en el cual la variable explicativa es categórica (por ejemplo varios tratamientos en un bioensayo, varias concentraciones diferentes de un fertilizante, o diferentes medicamentos, etc.) y además tenemos una o varias covariables, es decir una variable continua que puede influir en el resultado de la variable dependiente (la eficacia de un tratamiento puede depender de la edad del paciente, por ejemplo).

2 Nuestros datos

Veamos como son nuestros datos:

#Fijamos el directorio de trabajo
setwd(dir = "F:/R/MARKDOWN/ANCOVAyPERMUTATION")

#leemos el fichero de datos
dataper<-read.csv("lmperm.csv", sep=";")
str(dataper)
## 'data.frame':    50 obs. of  4 variables:
##  $ animal   : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ treatment: int  1 1 1 1 1 1 1 1 1 1 ...
##  $ behaviour: num  12.3 10.6 5.9 15.6 12.3 16.9 21.2 20 21.3 20 ...
##  $ sizemm   : num  4.5 4.9 5.6 3.9 3.9 6.1 5.8 4.2 4.6 4.1 ...
head(dataper)
##   animal treatment behaviour sizemm
## 1      1         1      12.3    4.5
## 2      2         1      10.6    4.9
## 3      3         1       5.9    5.6
## 4      4         1      15.6    3.9
## 5      5         1      12.3    3.9
## 6      6         1      16.9    6.1

Tenemos una variable dependiente cuantitativa (behaviour), unas concentraciones de tóxico (treatment) y el tamaño de los animales (sizemm).

3 Representación gráfica

Representamos gráficamente nuestros datos:

par(mfrow=c(1,1))
boxplot(behaviour ~ treatment,
        data = dataper, xlab="Treatment", ylab="behaviour")

Parece que los dos últimos tratamientos (4 y 5) presentan un comportamiento diferente al del control (1), pero vamos a analizarlo por medio de ANCOVA.

4 ANCOVA

Vamos a hacer un ANCOVA (Análisis de la Covarianza, fusión de ANOVA y regresión lineal múltiple). Es un modelo lineal general con una variable cuantitativa dependiente y uno o más factores:

library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
modancova <- aov(behaviour ~ sizemm + treatment, data=dataper)
summary(modancova)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## sizemm       1    382     382   1.432   0.237    
## treatment    1   9929    9929  37.205 1.9e-07 ***
## Residuals   47  12543     267                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Claramente la concentración de tóxico afecta a la variable dependiente comportamiento (behaviour) pero el tamaño de los animales parece no influir en este respuesta. Se pueden añadir más covariables después de la primera, pero el orden importa. Se debe poner la que sea más importante al principio.

5 ANCOVA y normalidad de los residuos

Vamos a ver si los residuos de nuestro modelo son normales, ya que es un requisito del ANCOVA:

qqnorm(resid(modancova), main="Normal Q-Q Plot")
qqline(resid(modancova), col="red")

Los datos parece que no son normales (los puntos se encuentran muy alejados de la línea). Vamos a optar por un modelo no paramétrico equivalente a ANCOVA pero que no requiere de normalidad en los residuos.


que raro que unos datos no salgan normales en ecología….


6 PERMUTATION test

Optamos por un test de permutaciones con aovp:

library(lmPerm)
modpermut <- aovp(behaviour ~ sizemm + treatment, data=dataper)
## [1] "Settings:  unique SS : numeric variables centered"
"Settings:  unique SS : numeric variables centered"
## [1] "Settings:  unique SS : numeric variables centered"
summary(modpermut)
## Component 1 :
##             Df R Sum Sq R Mean Sq Iter Pr(Prob)    
## sizemm       1    557.4     557.4   83   0.5542    
## treatment    1   9928.7    9928.7 5000   <2e-16 ***
## Residuals   47  12542.7     266.9                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La interpretación de los resultados es la misma que para el ANCOVA.

Espero que sea de utilidad.


R produce muchos dolores de cabeza….


7 CRÉDITOS

Álvaro Alonso Fernández

Departamento de Ciencias de la Vida

Universidad de Alcalá (España)