# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Primero cargamos la base de datos Toloaches
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 de los datos filtrando sólo aquellos datos del
# suelo pobre y eliminamos los Valores ausentes,NA, de la variable Flores
# con un argumento flores mayor hacer
DATOS2 <- subset(DATOS, SUELO == "P" & FLORES > 0)
# Adjuntamos esta nueva base de datos en su conjunto que hemos llamado a
# tus dos
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 la base de datos original con la función << detach() >>
detach(DATOS)
# Y removemos la base de datos original << rm() >> Así desvinculando y
# removiendo la base de datos original evitamos que se puedan utilizar los
# datos en este objeto. Recuerda que ambos objetos contiene exactamente
# los mismos nombres en sus columnas y podríamos fácilmente equivocaras. R
# por precedencia siempre utilizará primero las columnas en el objeto <<
# DATOS >> ya que fue el primero que se creó.
rm(DATOS)
# Ahora sí ya con el subconjunto creado y habiendo eliminado y desvincula
# el objeto << DATOS >> vamos a ajustar un modelo lineal con los objetos
# << FLORES >> como variable dependiente en función de los factores << LUZ
# >> , << HERBIVORIA >> y << MICORRIZA >>
M <- lm(FLORES ~ LUZ * HERBIVORIA * MICORRIZA)
# Ya con el modelo ajustado lo primero que debemos hacer es cerciorarnos
# que se cumplan los supuestos, es decir que haya homocedasticidad y
# normalidad en los residuales. Para esto utilizaremos la función << par()
# >> con el argumento << mfrow = c(2,2) >> para dividir la pantalla
# gráfica y poder desplegar los cuatro paneles concernientes con los
# residuales del modelo con la función << plot() >>
par(mfrow = c(2, 2))
plot(M)
# En el panel superior izquierdo podemos observar una tendencia en forma
# de embudo. Los residuales al lado izquierdo de este panel presentan poca
# variación y esta se va incrementando gradualmente hacia el lado derecho
# un patrón similar se observan el panel inferior izquierda. La normalidad
# parece ser adecuada a juzgar por el gráfico en el panel superior
# derecho, pero utilizaremos la prueba de Kolmogorov-Smirnov para
# verificarlo
ks.test(M$residuals, "pnorm")
## Warning: ties should not be present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: M$residuals
## D = 0.4732, p-value < 2.2e-16
## alternative hypothesis: two-sided
# Aunque los residuales aparentaba ser normales en el gráfico la prueba de
# Kolmogorov-Smirnov detecta un desvío significativo de la normalidad, un
# desvío fuerte ya que el valor de P es muy pequeño: 2.2 a la -16.
# Transformaremos nuevamente los datos usando el logaritmo de << FLORES >>
# multiplicado por mil. De esta forma, << FLORES >> que estaba en
# fracciones de gramo el multiplicarlo por 1000 lo colocaremos en
# miligramos y tomaremos el logaritmo de estos miligramos como variable de
# respuesta.
LOGFLORES <- log(FLORES * 1000)
# Ajustamos nuevamente el modelo lineal usando el logaritmo de << FLORES
# >> como variable dependiente y verificamos inmediatamente lo residuales
# de la misma forma que procedimos anteriormente
M <- lm(LOGFLORES ~ LUZ * HERBIVORIA * MICORRIZA)
par(mfrow = c(2, 2))
plot(M)
# Podemos observar que la distribución de los residuales en los paneles
# del lado izquierdo, superior e inferior una mejora ya que la amplitud de
# las lineas de puntos que se define es aproximadamente la misma a lo
# largo de todo el eje horizontal, el eje X, que son los valores ajustados
# del modelo.
# La prueba de Kolmogorov-Smirnov muestra ahora que no existe un desvío
# significativo de la normalidad para los residuales. Sin embargo si
# usásemos la prueba de Shapiro observaríamos un desvío significativo de
# la normalidad, aunque cabe mencionar que este desvió es apenas
# marginalmente significativo.
ks.test(M$residuals, "pnorm")
## Warning: ties should not be present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: M$residuals
## D = 0.108, p-value = 0.1687
## alternative hypothesis: two-sided
shapiro.test(M$residuals)
##
## Shapiro-Wilk normality test
##
## data: M$residuals
## W = 0.9747, p-value = 0.0402
# Una vez que sea cumplido con los supuestos del modelo procedemos
# explorar la tabla de ANOVA
anova(M)
## Analysis of Variance Table
##
## Response: LOGFLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## LUZ 1 10.8 10.76 18.24 4.5e-05 ***
## HERBIVORIA 1 0.0 0.01 0.01 0.9205
## MICORRIZA 1 2.8 2.82 4.78 0.0312 *
## LUZ:HERBIVORIA 1 0.5 0.50 0.86 0.3573
## LUZ:MICORRIZA 1 5.2 5.24 8.89 0.0036 **
## HERBIVORIA:MICORRIZA 1 0.1 0.14 0.24 0.6284
## LUZ:HERBIVORIA:MICORRIZA 1 0.0 0.00 0.00 0.9795
## Residuals 98 57.8 0.59
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Podemos posteriormente utilizar la función << step() >> para simplificar
# el modelo y eliminar de este los factores que no son significativos.
M2_STEP <- step(M)
## Start: AIC=-48.22
## LOGFLORES ~ LUZ * HERBIVORIA * MICORRIZA
##
## Df Sum of Sq RSS AIC
## - LUZ:HERBIVORIA:MICORRIZA 1 0.00039 57.8 -50.2
## <none> 57.8 -48.2
##
## Step: AIC=-50.22
## LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA +
## HERBIVORIA:MICORRIZA
##
## Df Sum of Sq RSS AIC
## - HERBIVORIA:MICORRIZA 1 0.14 58.0 -52.0
## - LUZ:HERBIVORIA 1 0.46 58.3 -51.4
## <none> 57.8 -50.2
## - LUZ:MICORRIZA 1 5.27 63.1 -43.0
##
## Step: AIC=-51.97
## LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA
##
## Df Sum of Sq RSS AIC
## - LUZ:HERBIVORIA 1 0.47 58.4 -53.1
## <none> 58.0 -52.0
## - LUZ:MICORRIZA 1 5.24 63.2 -44.8
##
## Step: AIC=-53.12
## LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:MICORRIZA
##
## Df Sum of Sq RSS AIC
## - HERBIVORIA 1 0.00 58.4 -55.1
## <none> 58.4 -53.1
## - LUZ:MICORRIZA 1 5.28 63.7 -45.9
##
## Step: AIC=-55.11
## LOGFLORES ~ LUZ + MICORRIZA + LUZ:MICORRIZA
##
## Df Sum of Sq RSS AIC
## <none> 58.4 -55.1
## - LUZ:MICORRIZA 1 5.28 63.7 -47.9
# Explorando el resultado de la simplificación vía la función << strep()
# >> tenemos que quedaron como factores principales << LUZ >> y <<
# MICORRIZA >> y también fue significativa la interacción << LUZ:MICORRIZA
# >>
anova(M2_STEP)
## 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
# Ahora haremos una simplificación a mano eliminando término por término
# del modelo empezando siempre por las interacciones más complejas en este
# caso la interacción triple y de ser posible eliminar esta interacción
# continuaremos con las interacciones dobles.
# Utilizamos la función << update() >> para crear la simplificación
M2 <- update(M, ~. - LUZ:HERBIVORIA:MICORRIZA)
# Y comparamos entre el modelo más complejo que contiene la interacción
# triple y el modelo simplificado que ya no contiene la interacción triple
# con la función << anova() >>. Este análisis de varianza entre modelos
# nos indica que el aumento en la suma de cuadrados residual una vez que
# removimos interacción triple no fue significativo, es decir que ese
# cambio en la suma de cuadrados es muy probable que ser observado
# simplemente azar. En este caso la probabilidad del cambio observado es
# de 0.9041, es decir 90% de las veces sólo por azar encontraremos un
# cambio de tal magnitud o menor en la suma de cuadrados.
anova(M2, M)
## Analysis of Variance Table
##
## Model 1: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA +
## HERBIVORIA:MICORRIZA
## Model 2: LOGFLORES ~ LUZ * HERBIVORIA * MICORRIZA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 99 57.8
## 2 98 57.8 1 0.00039 0 0.98
# Una vez eliminada la interacción triple ya tenemos un modelo más simple
# que hemos llamado << M2 >> y ahora pondremos a prueba si podemos
# eliminar cada una de las tres interacciones dobles que existen. Cada una
# de esas interacciones dobles se prueban contra el mismo modelo, << M2 >>
# Vamos a podemos eliminar la interacción << LUZ:HERBIVORIA >>
M3 <- update(M2, ~. - LUZ:HERBIVORIA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:MICORRIZA + HERBIVORIA:MICORRIZA
## Model 2: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA +
## HERBIVORIA:MICORRIZA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 100 58.3
## 2 99 57.8 1 0.456 0.78 0.38
# También vamos a poder eliminar la interacción << MICORRIZA:HERBIVORIA >>
M3 <- update(M2, ~. - MICORRIZA:HERBIVORIA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA
## Model 2: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA +
## HERBIVORIA:MICORRIZA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 100 58.0
## 2 99 57.8 1 0.139 0.24 0.63
# Pero no vamos a poder eliminar la interacción << LUZ:MICORRIZA << ya que
# la probabilidad de cambio es de 0.003375, es decir una probabilidad
# baja. El cambio observado en la suma de cuadrados residual al eliminar
# esta interacción no es un producto del azar.
M3 <- update(M2, ~. - LUZ:MICORRIZA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + HERBIVORIA:MICORRIZA
## Model 2: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA +
## HERBIVORIA:MICORRIZA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 100 63.1
## 2 99 57.8 1 5.27 9.02 0.0034 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Una vez realizadas las pruebas para cada una de las tres interacciones
# dobles procedemos a remover del modelo las dos interacciones dobles que
# se indicó podrían ser removidas del modelo,<< MICORRIZA:HERBIVORIA >> y
# << LUZ:HERBIVORIA >> Y verificamos mediante la prueba << anova() >> que
# efectivamente remover estas dos interacciones no genera un cambio
# significativo en la suma de cuadrados residual.
M3 <- update(M2, ~. - MICORRIZA:HERBIVORIA - LUZ:HERBIVORIA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:MICORRIZA
## Model 2: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:HERBIVORIA + LUZ:MICORRIZA +
## HERBIVORIA:MICORRIZA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 101 58.4
## 2 99 57.8 2 0.608 0.52 0.6
# Podemos observar la tabla de ANOVA de este modelo que ya sólo contiene
# una interacción doble y vemos que de los factores listados, <<
# HERBIVORIA >> parece no tener efecto significativo. Dado que este
# factor no forma parte de la interacción significativa podemos intentar
# removerlo del modelo. Los factores << LUZ >> y << MICORRIZA >> aunque no
# fuese significativos en la tabla no podríamos intentar removerlos del
# modelo ya que estos forman parte de la interacción doble que tiene
# efectos significativos dentro del modelo.
anova(M3)
## Analysis of Variance Table
##
## Response: LOGFLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## LUZ 1 10.8 10.76 18.60 3.8e-05 ***
## HERBIVORIA 1 0.0 0.01 0.01 0.9197
## MICORRIZA 1 2.8 2.82 4.87 0.0295 *
## LUZ:MICORRIZA 1 5.3 5.28 9.13 0.0032 **
## Residuals 101 58.4 0.58
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Removemos el factor principal << HERBIVORIA >> con la misma función <<
# update() >> y comparamos entre modelos
M4 <- update(M3, ~. - HERBIVORIA)
anova(M4, M3)
## Analysis of Variance Table
##
## Model 1: LOGFLORES ~ LUZ + MICORRIZA + LUZ:MICORRIZA
## Model 2: LOGFLORES ~ LUZ + HERBIVORIA + MICORRIZA + LUZ:MICORRIZA
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 102 58.4
## 2 101 58.4 1 0.0013 0 0.96
# Observamos que el cambio en la suma de cuadrados residual no es
# significativa por lo tanto dejamos el factor << HERBIVORIA >> fuera del
# modelo.
# Observamos la tabla de ANOVA de este modelo, el modelo final, ya que no
# vamos a poder remover ninguno de los otros factores.
anova(M4)
## 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
# Si comparamos este modelo simplificado a mano con el modelo simplificado
# con la función << step() >> vemos que coinciden perfectamente. No
# siempre es el caso de estas coincidencias, muchas veces la función <<
# step() >> es demasiado conservadora en remover algunos términos del
# modelo por lo que se sugiere hacerlo la simplificación a mano o a partir
# de dónde queda la función esté verificar que ya no es posible
# simplificar mas.
anova(M2_STEP)
## 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
# Ya con el modelo simplificado podemos explorar el predictor lineal, los
# coeficientes de cada una de las betas del modelo
summary(M4)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6529 0.1406 33.102 2.276e-56
## LUZPL -1.1029 0.2140 -5.154 1.254e-06
## MICORRIZASIN -0.7248 0.1971 -3.677 3.786e-04
## LUZPL:MICORRIZASIN 0.8997 0.2963 3.036 3.043e-03
# Ahora procederemos a hacer los vectores de ponderación para utilizar el
# predictor lineal colocando simplemente unos y ceros en las betas que
# corresponden según lo que queremos calcular. Primero vamos a extraer en
# un objeto llamado << BETAS >> los valores con sus nombres de este
# predictor lineal, que están almacenadas en la columna 1 de la tabla de
# coeficients, por eso usamos la notación con << [, 1] >>
BETAS <- summary(M4)$coefficients[, 1]
BETAS
## (Intercept) LUZPL MICORRIZASIN
## 4.6529 -1.1029 -0.7248
## LUZPL:MICORRIZASIN
## 0.8997
# Ahora si procedamos a crear los vectores de ponderación para la primera
# combinación, << ML:CON >> tenemos que colocar uno en el intercepto, cero
# en la beta que está etiquetada para << LUZPL >>, cero en la beta este
# etiquetada << MICORRIZASIN >>, y cero en la beta que está etiquetada <<
# LUZPL:MICORRIZASIN >>
ML_CM <- c(1, 0, 0, 0)
# Para la siguiente combinación <<ML:SM>> procedemos de manera similar,
# uno en el intercepto, cero en la beta etiquetada << LUZPL >>, uno en la
# beta etiquetada << MICORRIZASIN >>, y cero para la combinación <<
# LUZPL:MICORRIZASIN >>
ML_SM <- c(1, 0, 1, 0)
# Con la misma lógica hacemos los vectores de ponderación para las
# siguientes dos combinaciones, ahora las que están en poca luz con
# micorriza y sin micorriza.
PL_CM <- c(1, 1, 0, 0)
PL_SM <- c(1, 1, 1, 1)
# A partir de esas combinaciones podemos calcular ahora los vectores de
# ponderación para los grupos más grandes. Por ejemplo si promediamos las
# dos combinaciones que incluyen a mucha luz obtendremos el vector de
# ponderación para el nivel mucha luz, independientemente de si tenían o
# no micorrizas las plantas. De igual forma si promediamos las dos
# combinaciones en poca luz obtendremos vector de ponderación para el
# nivel de poca luz, si promediamos las dos combinaciones con micorriza
# obtendremos el vector el nivel con micorriza, y si promediamos las dos
# combinaciones sin micorriza tendremos el vector de ponderación para el
# nivel sin micorriza
ML <- (ML_CM + ML_SM)/2
PL <- (PL_CM + PL_SM)/2
CM <- (ML_CM + PL_CM)/2
SM <- (ML_SM + PL_SM)/2
# Cada uno de estos vectores de ponderación se utiliza para calcular los
# valores promedio correspondiente al multiplicar los valores de las betas
# por le vector de ponderación. A continuación ilustramos el ejemplo para
# el nivel de mucha luz. Hemos colocado en una tabla las << BETAS >> el
# vector de ponderación para mucha luz, y el producto entre las betas y
# este vector de ponderación
TABLA <- cbind(BETAS, ML, BETAS * ML)
TABLA
## BETAS ML
## (Intercept) 4.6529 1.0 4.6529
## LUZPL -1.1029 0.0 0.0000
## MICORRIZASIN -0.7248 0.5 -0.3624
## LUZPL:MICORRIZASIN 0.8997 0.0 0.0000
# Y sumando la columna tres de esta tabla, que es el producto de las betas
# por el vector de ponderación, obtendremos el valor promedio para el
# nivel mucha luz
sum(TABLA[, 3])
## [1] 4.29
# Comprobamos calculando el promedio de los datos y observamos que para el
# nivel de mucha luz, el promedio de los datos crudos es muy similar el
# valor al calculado con el predictor lineal.
tapply(LOGFLORES, LUZ, mean)
## ML PL
## 4.284 3.643
# Un ejemplo más, en este caso para el nivel sin micorriza. Se sigue
# exactamente la misma lógica, unir en una tabla las betas, el vector de
# ponderación y su producto. Sumamos la tercera columna para obtener el
# valor promedio en el nivel sin micorriza.
TABLA <- cbind(BETAS, SM, BETAS * SM)
sum(TABLA[, 3])
## [1] 3.826
# También comprobamos este cálculo y vemos que se parece mucho al valor
# calculado con los datos crudos
tapply(LOGFLORES, MICORRIZA, mean)
## CON SIN
## 4.177 3.836
### CONTRASTES, DIFERENCIAS ENTRE NIVELES-------------------------
# Así como podemos obtener los vectores de ponderación de los niveles de
# los factores principales promediando las combinaciones correspondientes,
# también podemos obtener los vectores de ponderación para estimar la
# diferencias entres cualquier par de medias. Por ejemplo para obtener el
# vector de ponderación para la diferencia entre la media de la
# combinación << ML:CM >> y la combinación << ML:SM >> simplemente debemos
# restar los dos vectores correspondientes. El orden que se haga la resta
# es irrelevante ya que simplemente cambiará el sigo de la diferencias
# estimada una vez que el vector de ponderación obtenido se multiplique
# por los valores de las betas y se sumen estas betas ponderadas. Así
# calculamos a diferencias << ML_CM - ML_SM >>
DIF_ML_CM_ML_S <- ML_CM - ML_SM
TABLA <- cbind(BETAS, DIF_ML_CM_ML_S, BETAS * DIF_ML_CM_ML_S)
TABLA
## BETAS DIF_ML_CM_ML_S
## (Intercept) 4.6529 0 0.0000
## LUZPL -1.1029 0 0.0000
## MICORRIZASIN -0.7248 -1 0.7248
## LUZPL:MICORRIZASIN 0.8997 0 0.0000
sum(TABLA[, 3])
## [1] 0.7248
# Y la comparamos con la resta hecha en el orden inverso y observamos que
# el valor es el mismo solo que con el signo opuesto, ML_SM - ML_CM
DIF_ML_CM_ML_S <- ML_SM - ML_CM
TABLA <- cbind(BETAS, DIF_ML_CM_ML_S, BETAS * DIF_ML_CM_ML_S)
TABLA
## BETAS DIF_ML_CM_ML_S
## (Intercept) 4.6529 0 0.0000
## LUZPL -1.1029 0 0.0000
## MICORRIZASIN -0.7248 1 -0.7248
## LUZPL:MICORRIZASIN 0.8997 0 0.0000
sum(TABLA[, 3])
## [1] -0.7248
# Podemos calcular otras diferencias
DIF_PL_CM_PL_S <- PL_CM - PL_SM
DIF_ML_CM_ML_S <- ML_CM - ML_SM
DIF_PL_ML <- PL - ML
# Estimando diferencias y sus errores estándares----------
# La parte mas útil del predictor lineal es utilizarlo para calcular las
# diferencias y sus errores estándares asociados para poder establecer si
# la estas son estadísticamente distintas de cero, es decir si la
# diferencia es estadísticamente significativa. Para lograr esto de forma
# sencilla vamos a utilizar la función << estimable() >> en la biblioteca
# << gmodels >>
library(gmodels)
# La función << estimables() >> utiliza dos argumentos, el modelo << M2 >>
# en este caso, y el vector de ponderación, el cual puede ser el vector de
# ponderación para una media o deferencia entre las medias, o como veremos
# mas adelante el vector de ponderación de una pendiente o diferencias
# entre pendientes. Así podemos calcular la media para << ML_CM >>, y la
# función << estimable() >> nos arrojará el estimado junto con su error
# estándar. La función estimable hace el cálculo de multiplicar el vector
# de ponderación por las betas del modelo y sumarlo, además calcula el
# error estándar.
estimable(M4, ML_CM)
## Estimate Std. Error t value DF Pr(>|t|)
## (1 0 0 0) 4.653 0.1406 33.1 102 0
# Ahora con la función << rbind() >> podemos unir varios vectores de
# ponderación en diferentes filas
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 1 1
# Y podemos incluir esta función dentro de la función << estimable() >>
# para calcular los estimados y errores estándares de mas de un vector a
# la vez
estimable(M4, 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
# Ahora hagamos las diferencias entre algunos vectores
estimable(M4, 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
# Ahora podemos graficar los valores estimados del modelo junto con sus
# errores estándares. De la función << estimable() >> podemos extraer la
# primer columna que contiene las medias
MEDIAS <- estimable(M4, rbind(ML_CM, ML_SM, PL_CM, PL_SM))[, 1]
# Y colocamos las cuatro medias en formato de matriz, dos columnas y dos
# filas
MEDIAS <- matrix(MEDIAS, 2, 2)
# Y en otro objeto podemos extraer los errores estándares que están en la
# columna 2
ERROR <- estimable(M4, rbind(ML_CM, ML_SM, PL_CM, PL_SM))[, 2]
# Y los errores también los colocamos en el formato de matriz
ERROR <- matrix(ERROR, 2, 2)
# Cargamos las funciones del curso
source("~/desktop/curso R 2012/funciones.r")
# Y con la función << barplot.error() >> representamos las cuatro medias
# con sus errores estándares. Este gráfico está en escala logarítmica que
# es la transformación que utilizamos de los datos
par(mfrow = c(1, 1))
barplot.error(t(MEDIAS), t(ERROR), names.arg = c("Mucha luz", "Poca luz"), col = c("grey40",
"grey80"), xlab = "", ylab = "ln Peso de las flores (mg)", ylim = c(3, 5))
# Añadimos una legenda para identificar las barras
legend(5, 4.8, c("Con micorriza", "Sin micorriza"), fill = c("grey40", "grey80"))
# Usndo el argumento << func = 'log' >> dentro de la función del gráfico
# podemos representar el resultado en la escala original
barplot.error(t(MEDIAS), t(ERROR), col = c("grey40", "grey80"), func = "log",
names.arg = c("Mucha luz", "Poca luz"), xlab = "", ylab = "Peso de las flores (mg)",
ylim = c(20, 140))
legend(5, 120, c("Con micorriza", "Sin micorriza"), fill = c("grey40", "grey80"))
