# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Nuevamente vamos a cargar la base de datos << TOLOACHE.CSV >> para esta
# sesión y hacemos un subconjunto para los datos de suelo pobre y vamos a
# quitar los valores ausentes de la variable << FLORES >>
DATOS <- read.table("~/desktop/curso R 2012/TOLOACHE.csv", header = TRUE, sep = ",")
attach(DATOS)
names(DATOS)
## [1] "NOTAS" "NUM" "INVERNADERO" "LUZ"
## [5] "MICORRIZA" "SUELO" "HERBIVORIA" "RAIZ"
## [9] "TALLO" "HOJAS" "FLORES" "CAPSULAS"
## [13] "SEMILLAS" "INMADUROS" "PESO_TOTAL" "NUM_FLORES"
## [17] "DIAS_A_LA_FLOR" "FRUT_INMADUROS" "FRUT_MADUROS" "TOT_FRUTOS"
## [21] "ABORTOS" "ALT_TOT" "HOJAS_TOT" "DUREZA"
## [25] "AREA" "PESO_HOJAS" "AFE" "LONG_HIFAS"
## [29] "HIFAS" "ESPORAS" "VESICULAS" "ARBUSCULOS"
## [33] "TOTAL"
DATOS2 <- subset(DATOS, SUELO == "P" & FLORES > 0)
attach(DATOS2)
## The following object(s) are masked from 'DATOS':
##
## ABORTOS, AFE, ALT_TOT, ARBUSCULOS, AREA, CAPSULAS,
## DIAS_A_LA_FLOR, DUREZA, ESPORAS, FLORES, FRUT_INMADUROS,
## FRUT_MADUROS, HERBIVORIA, HIFAS, HOJAS, HOJAS_TOT, INMADUROS,
## INVERNADERO, LONG_HIFAS, LUZ, MICORRIZA, NOTAS, NUM,
## NUM_FLORES, PESO_HOJAS, PESO_TOTAL, RAIZ, SEMILLAS, SUELO,
## TALLO, TOT_FRUTOS, TOTAL, VESICULAS
# Recuerden que debemos desvincular y remover la base de datos original <<
# DATOS >>
detach(DATOS)
rm(DATOS)
LOGFLORES <- log(FLORES * 1000)
# Ahora vamos a ajustar el siguiente modelo línea. El logaritmo de la
# variable << FLORES >> que ya habíamos utilizado antes con el factor luz
# y anidado en este factor la interacción entre micorriza y el peso de los
# tallos. Este modelo al incluir una combinación de predictores
# categóricos (<< LUZ >> y << MICORRIZA >>) y continuos ( << PESO_HOJA >>:
# peso de las hojas) representa un análisis de covarianza en donde además
# de los ya analizados intercepto y sus diferencias, es decir entre las
# medias que ya hemos calculado en ejercicios anteriores s podemos ahora
# también calcular la diferencia entre pendientes de diferentes rectas
M3 <- lm(LOGFLORES ~ LUZ/(MICORRIZA * PESO_HOJAS))
# Este modelo tiene exactamente los mismos supuestos que los modelos
# anteriores, que los residuales se distribuyen normalmente y que sean
# homocedásticos. Estos supuestos los verificamos con la función del
# modelo correspondiente, y con la función de Kolmogorov-Smirnov para la
# normalidad
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(M3)
ks.test(M3$residuals, "pnorm")
## Warning: ties should not be present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: M3$residuals
## D = 0.1156, p-value = 0.1174
## alternative hypothesis: two-sided
# Como podemos apreciar apreciar tanto en los gráficos como la prueba de
# Kolmogorov-Smirnov se cumplen los supuestos del modelo, la
# homocedasticidad y la normalidad que lo residuales.
anova(M3)
## Analysis of Variance Table
##
## Response: LOGFLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## LUZ 1 10.8 10.76 22.89 6.1e-06 ***
## LUZ:MICORRIZA 2 8.1 4.05 8.62 0.00036 ***
## LUZ:PESO_HOJAS 2 8.4 4.22 8.98 0.00026 ***
## LUZ:MICORRIZA:PESO_HOJAS 2 3.9 1.97 4.18 0.01810 *
## Residuals 98 46.1 0.47
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ahora exploramos la tabla de coeficientes donde encontramos en los
# estimados los valores de las betas del predictor lineal. En las
# etiquetas correspondientes a estas betas podemos rápidamente diferenciar
# dos tipos algunas de ellas se caracterizan por estar siempre asociadas a
# la variable continua, en este caso << PESO_HOJAS >>. Todas estas betas
# que están asociadas a la variable continua son betas relativas a las
# pendientes de las lineas de rregresion. Por otra parte todas las otras
# betas que NO están asociadas a la variable continua son betas relativas
# al intercepto, es decir a las medias.
summary(M3)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4116 0.2764 15.9617 5.209e-29
## LUZPL -1.5675 0.3694 -4.2432 5.005e-05
## LUZML:MICORRIZASIN -1.4494 0.4084 -3.5489 5.958e-04
## LUZPL:MICORRIZASIN 0.8798 0.4334 2.0298 4.509e-02
## LUZML:PESO_HOJAS 3.5346 3.5930 0.9838 3.277e-01
## LUZPL:PESO_HOJAS 13.0499 3.6373 3.5878 5.225e-04
## LUZML:MICORRIZASIN:PESO_HOJAS 18.9283 7.3026 2.5920 1.100e-02
## LUZPL:MICORRIZASIN:PESO_HOJAS -13.0212 10.1594 -1.2817 2.030e-01
# Una vez que hemos diferenciado entre estos dos tipos de betas es fácil
# construir el vector de ponderación para calcular los interceptos o las
# pendientes de cada una de las combinaciones como se muestra a
# continuación.
# Vamos a construir los contrastes dentro de la matriz. Será más fácil
# seguir las instrucciones para generar la matriz una vez que visualicen
# el resultado. El orden en que están ingresados los datos en la matriz es
# fila por fila, y esto lo hacemos definiendo en la función << matrix() >>
# el argumento << byrow= T >>. Esta Matriz tiene ocho filas
# correspondientes a cada una de las combinaciones que vamos a calcular,
# cuatro combinaciones para los interceptos y las mismas cuatro
# combinaciones para las pendientes. Esta matriz también tiene ocho
# columnas ya que existe ocho betas en el predictor lineal
CONTRASTES <- matrix(c(1, rep(0, 7), 1, 0, 1, rep(0, 5), 1, 1, rep(0, 6), 1,
1, 0, 1, rep(0, 4), rep(0, 4), 1, 0, 0, 0, rep(0, 4), 1, 0, 1, 0, rep(0,
5), 1, 0, 0, rep(0, 5), 1, 0, 1), nrow = 8, ncol = 8, byrow = T)
# Con las funciones << rownames() >> y << colnames() >> asignamos los
# nombres de las filas y las columnas. En las filas tenemos cada una de
# las combinaciones para las cuales estamos creando vector de ponderación,
# y en las columnas tenemos las etiquetas del predictor lineal
rownames(CONTRASTES) <- c("I_ML&CM", "I_ML&SM", "I_PL&CM", "I_PL&SM", "P_ML&CM",
"P_ML&SM", "P_PL&CM", "P_PL&SM")
# Las etiquetas del predictor lineal deben coincidir con los nombres que
# tiene asignados cada una de las betas, y esto lo podemos hacer
# fácilmente extrayendo los nombres de las betas desde la tabla de
# coeficientes usando la misma función << rownames() >> pero esta vez para
# asignarlos como nombres de las columnas
colnames(CONTRASTES) <- rownames(summary(M3)$coefficients)
CONTRASTES
## (Intercept) LUZPL LUZML:MICORRIZASIN LUZPL:MICORRIZASIN
## I_ML&CM 1 0 0 0
## I_ML&SM 1 0 1 0
## I_PL&CM 1 1 0 0
## I_PL&SM 1 1 0 1
## P_ML&CM 0 0 0 0
## P_ML&SM 0 0 0 0
## P_PL&CM 0 0 0 0
## P_PL&SM 0 0 0 0
## LUZML:PESO_HOJAS LUZPL:PESO_HOJAS LUZML:MICORRIZASIN:PESO_HOJAS
## I_ML&CM 0 0 0
## I_ML&SM 0 0 0
## I_PL&CM 0 0 0
## I_PL&SM 0 0 0
## P_ML&CM 1 0 0
## P_ML&SM 1 0 1
## P_PL&CM 0 1 0
## P_PL&SM 0 1 0
## LUZPL:MICORRIZASIN:PESO_HOJAS
## I_ML&CM 0
## I_ML&SM 0
## I_PL&CM 0
## I_PL&SM 0
## P_ML&CM 0
## P_ML&SM 0
## P_PL&CM 0
## P_PL&SM 1
# Habilitamos la biblioteca << gmodels >> dos para poder usar la función
# << estimable() >>
library(gmodels)
estimable(M3, CONTRASTES)
## Estimate Std. Error t value DF Pr(>|t|)
## I_ML&CM 4.41155 0.2764 15.961749 98 0.000e+00
## I_ML&SM 2.96217 0.3007 9.851800 98 2.220e-16
## I_PL&CM 2.84407 0.2451 11.603605 98 0.000e+00
## I_PL&SM 3.72383 0.3575 10.417474 98 0.000e+00
## P_ML&CM 3.53462 3.5930 0.983764 98 3.277e-01
## P_ML&SM 22.46288 6.3576 3.533224 98 6.280e-04
## P_PL&CM 13.04985 3.6373 3.587816 98 5.225e-04
## P_PL&SM 0.02866 9.4860 0.003021 98 9.976e-01
# y obtenemis los estimados para los cuatro interceptos y las cuatro
# pendientes correspondientes
# Alternativamente podríamos construir esta matriz de vectores de
# ponderación en una hoja de cálculo como Excel con los nombres de cada
# una de las combinaciones en la primer columna y los nombres de las betas
# como encabezados en el resto de las columnas. Una vez que tenemos
# preparada esta tabla de Excel en un formato que puedan leer R, ya sea <<
# CSV >> o << TXT >> importamos la tabla de la misma forma que importamos
# la base de datos original, y añadimos un argumento más <<row.names=1 >>
# que indica que indica que la primer columna tienen las etiquetas de las
# filas.
CONTRASTES <- read.table("~/desktop/curso R 2012/CONTRASTES.csv", header = TRUE,
sep = ",", row.names = 1)
# Y sustituimos los nombres de las columnas que hayamos usado para hacer
# la tabla en excel por los nombres exactos de las betas del predictor
# lineal, y esto lo hacemos de la misma forma que lo hicimos
# anteriormente.
colnames(CONTRASTES) <- rownames(summary(M3)$coefficients)
# Como observarán la tabla preparada en Excel contiene además de los
# cuatro estimados interceptos y los cuatro estimados de pendientes, las
# diferencias entre pendientes interceptas y los vectores para los niveles
# superiores, por ejemplo de mucha luz, poca luz, y la diferencia entre
# los niveles de luz
CONTRASTES
## (Intercept) LUZPL LUZML:MICORRIZASIN LUZPL:MICORRIZASIN
## I_ML&CM 1 0 0.0 0.0
## I_ML&SM 1 0 1.0 0.0
## I_PL&CM 1 1 0.0 0.0
## I_PL&SM 1 1 0.0 1.0
## P_ML&CM 0 0 0.0 0.0
## P_ML&SM 0 0 0.0 0.0
## P_PL&CM 0 0 0.0 0.0
## P_PL&SM 0 0 0.0 0.0
## I_D_1_2 0 0 -1.0 0.0
## I_D_3_4 0 0 0.0 -1.0
## P_D_5_6 0 0 0.0 0.0
## P_D_7_8 0 0 0.0 0.0
## ML 1 0 0.5 0.0
## PL 1 1 0.0 0.5
## D_ML_PL 0 -1 0.5 -0.5
## LUZML:PESO_HOJAS LUZPL:PESO_HOJAS LUZML:MICORRIZASIN:PESO_HOJAS
## I_ML&CM 0 0 0
## I_ML&SM 0 0 0
## I_PL&CM 0 0 0
## I_PL&SM 0 0 0
## P_ML&CM 1 0 0
## P_ML&SM 1 0 1
## P_PL&CM 0 1 0
## P_PL&SM 0 1 0
## I_D_1_2 0 0 0
## I_D_3_4 0 0 0
## P_D_5_6 0 0 -1
## P_D_7_8 0 0 0
## ML 0 0 0
## PL 0 0 0
## D_ML_PL 0 0 0
## LUZPL:MICORRIZASIN:PESO_HOJAS
## I_ML&CM 0
## I_ML&SM 0
## I_PL&CM 0
## I_PL&SM 0
## P_ML&CM 0
## P_ML&SM 0
## P_PL&CM 0
## P_PL&SM 1
## I_D_1_2 0
## I_D_3_4 0
## P_D_5_6 0
## P_D_7_8 -1
## ML 0
## PL 0
## D_ML_PL 0
estimable(M3, CONTRASTES)
## Estimate Std. Error t value DF Pr(>|t|)
## I_ML&CM 4.41155 0.2764 15.961749 98 0.000e+00
## I_ML&SM 2.96217 0.3007 9.851800 98 2.220e-16
## I_PL&CM 2.84407 0.2451 11.603605 98 0.000e+00
## I_PL&SM 3.72383 0.3575 10.417474 98 0.000e+00
## P_ML&CM 3.53462 3.5930 0.983764 98 3.277e-01
## P_ML&SM 22.46288 6.3576 3.533224 98 6.280e-04
## P_PL&CM 13.04985 3.6373 3.587816 98 5.225e-04
## P_PL&SM 0.02866 9.4860 0.003021 98 9.976e-01
## I_D_1_2 1.44938 0.4084 3.548902 98 5.958e-04
## I_D_3_4 -0.87976 0.4334 -2.029807 98 4.509e-02
## P_D_5_6 -18.92827 7.3026 -2.591975 98 1.100e-02
## P_D_7_8 13.02119 10.1594 1.281687 98 2.030e-01
## ML 3.68686 0.2042 18.055084 98 0.000e+00
## PL 3.28395 0.2167 15.153668 98 0.000e+00
## D_ML_PL 0.40292 0.2978 1.353160 98 1.791e-01
# Como podemos observar la tabla de contrastes los cuatro interceptos de
# las combinaciones son distintos a cero, mientras que en las pendientes
# tenemos que dos son estadísticamente distintas a cero las de las
# combinaciones mucha luz sin micorriza y la de poca luz con micorriza
# mientras que las otras dos combinaciones las pendientes no son distintas
# a cero.
# Para las diferencias entre los interceptos dentro de cada nivel de luz
# podemos observar que con y sin micorriza son estadísticamente diferentes
# en ambos niveles de luz, mientras que para las pendientes sólo hay
# diferencias dentro de la combinación de mucha luz. También podemos
# observar una diferencia general entre los niveles de luz la cual está
# representada en la última fila esta tarde contrastes
# Para realizar el gráfico correspondiente especialmente para ilustrar las
# pendientes utilizaremos nuevamente los repetidos “loops” que hemos
# estado usando.
par(mar = c(3, 3, 1, 1), mfrow = c(1, 1))
plot(LOGFLORES ~ PESO_HOJAS, type = "n")
COL <- c("cyan", "blue", "green", "red")
k <- 0
for (i in levels(LUZ)) {
for (j in levels(MICORRIZA)) {
k <- k + 1
points(LOGFLORES[LUZ == i & MICORRIZA == j] ~ PESO_HOJAS[LUZ == i &
MICORRIZA == j], pch = 19, col = COL[k])
PH <- PESO_HOJAS[LUZ == i & MICORRIZA == j]
PH <- seq(min(PH), max(PH), length = 30)
ND <- data.frame(LUZ = i, MICORRIZA = j, PESO_HOJAS = PH)
PRED <- predict(M3, ND, se = TRUE)
lines(PH, PRED$fit, col = COL[k])
}
}
legend(0.1, 3.5, c("MLC", "MLS", "PLC", "PLS"), col = c("cyan", "blue", "green",
"red"), pch = 19, ncol = 2, y.intersp = 0.8, x.intersp = 0.3, bty = "n")
