# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Nuevamente cargamos la base de datos << TOLOACHES.CSV >>
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"
# Hacemos un subconjunto para los datos en suelo pobre y quitamos los
# valores ausentes de la variable << FLORES >> que usaremos como variable
# de respuesta
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
# Desvinculamos y removemos la base de daos original
detach(DATOS)
rm(DATOS)
# Transformamos la variable de respuesta en logaritmo después de haberla
# multiplicado por mil, ahora está expresada en ln de miligramos
LOGFLORES <- log(FLORES * 1000)
# Con esto compararemos dos sintaxis diferentes dentro de los modelos. La
# primera es la sintaxis ya conocida, usando el asterisco entre los
# factores que explican. Esta notación como vimos indica que se ajuste el
# modelo en función de los factores << LUZ >>, << MICORRIZA >> y su
# interacción << LUZ:MICORRIZA >> y en esta interacción existen cuatro
# posibles combinaciones, << ML:CON >>, << ML:SIN >>, << PL:CON >> y <<
# PL:SIN >>. Este es un diseño factorial y su uso implica que cualquier
# comparación entre las cuatro posibles combinaciones de la interacción
# son válidas y sobre todo que tienen sentido biológico. Además implica
# que cada uno de los factores individualmente tambien tienen sentido
# biológico para explicar nuestros resultados.
M1 <- lm(LOGFLORES ~ LUZ * MICORRIZA)
# Note como en la tabla de ANOVA se muestran, los factores individuales y
# la interacción con sus respectivas sumas de cuadrados, grados de
# libertad, cuadrados medios, estimados de F, y sus probabilidades
anova(M1)
## Analysis of Variance Table
##
## Response: LOGFLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## LUZ 1 10.8 10.76 18.78 3.4e-05 ***
## MICORRIZA 1 2.8 2.82 4.93 0.029 *
## LUZ:MICORRIZA 1 5.3 5.28 9.22 0.003 **
## Residuals 102 58.4 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# La segunda notación es muy similar, simplemente se sustituye el
# asterisco por una diagonal. Esto se define como un modelo anidado, en
# donde el factor << MICORRIZA >> está anidado dentro del factor << LUZ
# >>. Esto quiere decir que, el factor << MICORRIZA >> por si solo no
# tiene sentido biológico, y que los efectos de este factor deben ser
# siempre considerados dentro de cada uno de los niveles del factor << LUZ
# >>, << ML >> y << PL >>. Esto también implica que en la interacción ya
# nos son válidas o carecen de sentido biológico las comparaciones entre
# de los diferentes niveles del factor << LUZ >>. En este modelo, dentro
# de la interacción las únicas comparaciones que tienen sentido son << CON
# >> y << SIN >> (micorriza) dentro de cada nivel de luz, es decir solo se
# debe comparar << ML:CON >> vs << ML:SIN >>, y << PL:CON >> vs << PL:SIN
# >>.
M2 <- lm(LOGFLORES ~ LUZ/MICORRIZA)
# Note como en la tabla de ANOVA ya no aparece el factor << MICORRIZA >>
# como variable principal, solo aparece en la interacción. También note el
# cambio en los grados de libertad en la interacción de este nuevo modelo
# comparado con la interacción del modelo previo. Ahora tenemos dos grados
# de libertad, ya que el grado de libertad que existía para el factor <<
# MICORRIZA >> en el modelo factorial se suma al de la interacción en el
# modelo anidado. Lo mismo ocurre con la suma de cuadrados, en el modelo
# anidado la suma de cuadrados de la interacción es mayor que en el modelo
# factorial por que incluye lo de la interacción mas lo del factor <<
# MICORRIZA >> como variable principal.
anova(M2)
## Analysis of Variance Table
##
## Response: LOGFLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## LUZ 1 10.8 10.76 18.78 3.4e-05 ***
## LUZ:MICORRIZA 2 8.1 4.05 7.07 0.0013 **
## Residuals 102 58.4 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Como ahora las comparaciones entre << CON >> y << SIN >> micorriza solo
# tienen sentido dentro de cada nivel de << LUZ >> en el resumen del
# modelo, en específico, dentro del predictor lineal, ahora existe una
# beta etiquetada para < ML:SIN >> y otro para << PL:SIN >>, y no existe
# una beta solo etiquetada para << SIN >> como ocurren con el modelo
# factorial
summary(M2)
##
## Call:
## lm(formula = LOGFLORES ~ LUZ/MICORRIZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6255 -0.5474 0.0476 0.5638 1.4607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.653 0.141 33.10 < 2e-16 ***
## LUZPL -1.103 0.214 -5.15 1.3e-06 ***
## LUZML:MICORRIZASIN -0.725 0.197 -3.68 0.00038 ***
## LUZPL:MICORRIZASIN 0.175 0.221 0.79 0.43117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.757 on 102 degrees of freedom
## Multiple R-squared: 0.244, Adjusted R-squared: 0.222
## F-statistic: 11 on 3 and 102 DF, p-value: 2.63e-06
summary(M1)
##
## Call:
## lm(formula = LOGFLORES ~ LUZ * MICORRIZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6255 -0.5474 0.0476 0.5638 1.4607
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.653 0.141 33.10 < 2e-16 ***
## LUZPL -1.103 0.214 -5.15 1.3e-06 ***
## MICORRIZASIN -0.725 0.197 -3.68 0.00038 ***
## LUZPL:MICORRIZASIN 0.900 0.296 3.04 0.00304 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.757 on 102 degrees of freedom
## Multiple R-squared: 0.244, Adjusted R-squared: 0.222
## F-statistic: 11 on 3 and 102 DF, p-value: 2.63e-06
# Note que entre estos dos predictores lineales el número total de betas
# no cambia, pero si sus etiquetas. Esto nos lleva a que los vectores de
# ponderación para las combinaciones << LUZ >> con << MICORRIZA >> cambien
# según sea el modelo utilizados.
# Estos son los vectores de ponderación para modelo anidado
ML_CM <- c(1, 0, 0, 0)
ML_SM <- c(1, 0, 1, 0)
PL_CM <- c(1, 1, 0, 0)
PL_SM <- c(1, 1, 0, 1)
# Una vez hecho estos vectores, se procede de la misma forma que ya lo
# hemos hecho para obtener los vectores de ponderación para los niveles en
# los factores principales, en este caso solo en el factor << LUZ >> ya
# que como mencionamos el factor << MICORRIZA >> carece de sentido como
# factor principal
ML <- (ML_CM + ML_SM)/2
PL <- (PL_CM + PL_SM)/2
# Las diferencias entre los niveles de << MICORRIZA >> dentro de cada
# nivel de luz también e hacen de la misma forma que ya lo hemos hecho
DIF_ML_CM_ML_S <- ML_CM - ML_SM
DIF_PL_CM_PL_S <- PL_CM - PL_SM
# Lo mismo entre los dos niveles de << LUZ >>
DIF_PL_ML <- PL - ML
# Y los contrastes los obtenemos con la función << estimable() >> que se
# encuentra en la biblioteca << gmodels >>
library(gmodels)
estimable(M2, ML_CM)
## Estimate Std. Error t value DF Pr(>|t|)
## (1 0 0 0) 4.653 0.1406 33.1 102 0
rbind(ML_CM, ML_SM, PL_CM, PL_SM)
## [,1] [,2] [,3] [,4]
## ML_CM 1 0 0 0
## ML_SM 1 0 1 0
## PL_CM 1 1 0 0
## PL_SM 1 1 0 1
estimable(M2, rbind(ML_CM, ML_SM, PL_CM, PL_SM))
## Estimate Std. Error t value DF Pr(>|t|)
## ML_CM 4.653 0.1406 33.10 102 0
## ML_SM 3.928 0.1382 28.42 102 0
## PL_CM 3.550 0.1614 22.00 102 0
## PL_SM 3.725 0.1514 24.60 102 0
estimable(M2, rbind(DIF_ML_CM_ML_S, DIF_PL_CM_PL_S, DIF_PL_ML))
## Estimate Std. Error t value DF Pr(>|t|)
## DIF_ML_CM_ML_S 0.7248 0.1971 3.6770 102 3.786e-04
## DIF_PL_CM_PL_S -0.1749 0.2213 -0.7903 102 4.312e-01
## DIF_PL_ML -0.6531 0.1482 -4.4077 102 2.589e-05
# Podemos extraer las medias y errores estimados del predictor lineal
MEDIAS <- estimable(M2, rbind(ML_CM, ML_SM, PL_CM, PL_SM))[, 1]
ERROR <- estimable(M2, rbind(ML_CM, ML_SM, PL_CM, PL_SM))[, 2]
# y graficarlos
source("~/desktop/curso R 2012/funciones.r")
barplot.error(MEDIAS, ERROR, col = 1:4)
legend(5, 4.5, c("ML:Con", "ML:Sin", "PL:COM", "PL:SIN"), fill = 1:4)
# Dada la naturaleza del analisis, << MICORRIZA >> anidado en << LUZ >> es
# mejor representar un gráfico con dos pares de barras, en donde cada
# pareja de barras son los niveles del factor anidado y los dos grupos de
# barras son os niveles del factor que anida, en este caso << LUZ >>. Para
# esto necesitamos reacomodar las medias y errores en un matriz de la
# forma siguiente
MEDIAS <- matrix(MEDIAS, 2, 2)
ERROR <- matrix(ERROR, 2, 2)
# En este arreglo, cada fila de la matriz será un grupo de barras. Por el
# orden en que incluimos los vectores de ponderación en la función <<
# estimable() >> los grupos estarían definido para cada nivel del factor
# micorriza.
MEDIAS
## [,1] [,2]
## [1,] 4.653 3.550
## [2,] 3.928 3.725
ERROR
## [,1] [,2]
## [1,] 0.1406 0.1614
## [2,] 0.1382 0.1514
# Sin embargo como deseamos hacer grupos por los niveles del factor << LUZ
# >> al momento de hacer el gráfico transpondremos la matriz con la
# función << t() >> , es decir las filas se vuelven columnas y las
# columnas se vuelven filas. Además usamos el argumento << func='los' >>
# para representar el gráfico en la escala original
barplot.error(t(MEDIAS), t(ERROR), col = c("grey40", "grey80"), func = "log",
names.arg = c("ML", "PL"))
legend(5, 120, c("Con", "Sin"), fill = c("grey40", "grey80"))
