# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Vamos a cargar la base de datos Toloache para hacer ahora análisis de
# varianza
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"
# También cargamos el archivo de funciones del curso para poder utilizar
# una de ellas que presenta gráficamente la suma de cuadrados A novas de 1
# y 2 vías y también para a ANOVAS anidadas
source("~/desktop/curso R 2012/funciones.r")
# Vamos también a crear un subconjunto de los datos usando la función <<
# subset() >> y escogiendo sólo aquellos datos Que tienen la
# característica sin herbivoría y sin micorriza y que están en el
# invernadero numero uno.La función utiliza como argumentos una base de
# datos << DATOS >> y una serie de atributos que van a definir al
# subconjunto.
DATOS2 <- subset(DATOS, HERBIVORIA == "SH" & MICORRIZA == "SIN" & INVERNADERO ==
1)
# Adjuntamos la nueva base de datos << DATOS2 >> y ahora tenemos los
# nombres de las columnas disponibles. Note que tantos << DATOS >> como <<
# DATOS2 >> tienen exactamente los mismos nombres en las columnas por lo
# cual al solicitar información de una de ellas R siempre mostrará la
# columna correspondiente en << DATOS >> y no en << DATOS2 >>.
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
# Para forzar a R a que la información que utilice sea la que está
# almacenada en << DATOS2 >> debemos primero remover la base de datos
# original llamada datos y también deshacer la función de adjuntar <<
# attach() >>. Para remover la base de datos usamos la función << rm() >>
# que es una abreviación de remover (remove) y para desvincular los
# nombres de la columnas usamos la función detach().
rm(DATOS)
detach(DATOS)
# Ahora usemos la función << sc() >> para hacer una representación gráfica
# de la suma de cuadrados en un ANOVA de dos días. Esta función tiene tres
# argumentos, la variable continua que se analizarían el análisis de
# varianza seguido de dos factores que representan los dos niveles. En
# este ejemplo usaremos la variable continua << FLORES >> y los factores
# << LUZ >> y << SUELO >>
sc(FLORES, LUZ, SUELO)
# Ahora si procedamos hacer los cálculos para un análisis de varianza. Qué
# es lo que representa el análisis de varianza, es la estimación de la
# distancia que existe entre las medias de los diferentes tratamientos o
# combinaciones de tratamientos a la media general. Bajo la hipótesis
# nula, la distancia de las medias de cualquier combinación de
# tratamientos a la media general debe ser cero ya que la hipótesis nula
# establece que todas las medias son iguales a las medias a la media
# general. Y precisamente el análisis de varianza mide este desvío de las
# medias de los tratamientos o combinación de tratamientos respecto a la
# media general. Una propiedad de los análisis de varianza es que la suma
# de cuadrados total, es decir los desvíos de los datos a la media general
# elevados al cuadrado y sumados sobre todo los datos, es aditiva, es
# decir, se puede repartir y partir la variase entre los diferentes
# tratamientos pero la suma de que todos los tratamientos o combinaciones
# de tratamiento siempre va a dar la suma de cuadrados total.
# Entonces para ejemplificar los cálculos de un análisis de varianza.
# Primero calculemos la media general de los datos de << FLORES >> que es
# el peso de las flores
MEDIA_G <- mean(FLORES)
# Ya con la media general podemos calcular fácilmente la suma de cuadrados
# total, es decir cada dato menos la media general, elevado al cuadrado y
# sumado sobre todo los datos. Esta suma de cuadrados está representado en
# el panel superior izquierdo de la figura.
TOTAL <- sum((FLORES - MEDIA_G)^2)
# Ahora procedamos a calcular la suma de cuadrados que correspondería al
# factor << LUZ >>. El factor << LUZ >> tiene dos niveles, mucha luz (ML)
# y poca luz (PL), y la suma de cuadrados del factor será la distancia
# entre el promedio de cada uno de los niveles del factor respecto a la
# media general elevada cuadrado y multiplicado por el número de datos que
# se usan en el cálculo de la media de cada uno de los dos niveles del
# factor, en este caso << LUZ >>, y una vez multiplicado se suman los
# valores resultante y eso nos da la suma de cuadrados del factor << LUZ
# >>. La función << tapply() >> La utilizamos para calcular el valor
# promedio de cada uno de los dos niveles del factor << LUZ >>, y a estas
# medias le restamos el objeto << MEDIA_G >> que es la media general, y a
# esas diferencias las elevamos al cuadrado, para multiplicarlas
# inmenditamente por el resultado de usar la función << table() >> con el
# factor << LUZ >> lo que nos da un conteo de cuantos datos hay dentro de
# cada nivel de luz, es decir cuantos datos se usan en el cálculo de la
# media de cada uno de los niveles. Finalmente la función << sum() >> que
# está por fuera de todas estos elementos ya utilizados nos da el
# resultado total de la suma de cuadrados del factor << LUZ >>. Esta suma
# de cuadrados está representada en el panel de medio del lado izquierdo.
# Note que en este caso las dos medias de los niveles de luz están
# prácticamente encima de la media general por lo cual gráficamente no se
# aprecia la distancia entre estas medias.
LZ <- sum(((tapply(FLORES, LUZ, mean) - MEDIA_G)^2) * table(LUZ))
# Exactamente con la misma lógica calculamos la suma de cuadrados para el
# factor << SUELO >>, usamos las mismas funciones y los mismos criterios.
# Esta suma de cuadrados está representada en el panel de medio del lado
# derecho. Este caso se puede apreciar claramente la distancia entre la
# media de cada nivel del factor suelo y la media general. Los números
# precedidos de una << x >> que aparecen junto a las medias representan el
# número total de datos que se usa en en el cálculo de la media, y es el
# valor por el cual se debe de multiplicar la suma de cuadrados.
SL <- sum(((tapply(FLORES, SUELO, mean) - MEDIA_G)^2) * table(SUELO))
# Para el caso de la interacción entre los factores << LUZ >> y << SUELO
# >>, lo cual quiere decir los efectos de las combinaciones de todos los
# niveles del factor << LUZ >> con todos los niveles del factor << SUELO
# >> lo cual en este ejemplo da cuatro combinaciones posibles. La única
# diferencia en estas funciones es el uso de la función << list() >>
# dentro de la función << tapply() >> lo cual nos permite calcular los
# promedios para cada una de las cuatro combinaciones entre los niveles de
# << LUZ >> y los niveles de << SUELO >>. Por otra parte una vez que ya
# hemos hecho estos cálculos y se hizo la suma sobre todas las diferencias
# de las medias elevada al cuadrado y ponderadas por el número de datos en
# cada combinación debemos restar la suma de cuadrados del factor << SUELO
# >> y el factor << LUZ >> para dejar únicamente la suma de cuadrados que
# correspondería a la combinación de los dos factores. Esta representación
# de suma de cuadrados está en la parte inferior del lado izquierdo.
INT <- sum(((tapply(FLORES, list(LUZ, SUELO), mean) - MEDIA_G)^2) * table(LUZ,
SUELO)) - SL - LZ
# Ya teniendo esta sumas de cuadrados solo nos resta calcular la suma de
# cuadrados residual lo cual podemos hacer fácilmente ya que será la suma
# de cuadrados total menos la suma de cuadrados del factor << LUZ >> menos
# la suma de cuadrados del factor << SUELO >> menos la suma de cuadrados
# de la interacción << LUZ:SUELO >>. La suma de cuadrados residual se
# encuentra representada en la parte superior del lado derecho.
RES <- TOTAL - LZ - SL - INT
# Ahora usemos la función << lm() >> para generar el análisis de varianza
# usando exactamente la misma sintaxis que habíamos usado para las
# regresiones.
M <- lm(FLORES ~ LUZ * SUELO)
# Podemos ver el resumen del modelo << summary() >> y ver que tiene
# exactamente la misma estructura que el de una regresión lineal.
summary(M)
##
## Call:
## lm(formula = FLORES ~ LUZ * SUELO)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2720 -0.0390 -0.0150 0.0086 0.5000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.04000 0.05615 0.71 0.48
## LUZPL -0.00857 0.08220 -0.10 0.92
## SUELOR 0.52200 0.09054 5.77 8.5e-06 ***
## LUZPL:SUELOR -0.07343 0.12652 -0.58 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.159 on 22 degrees of freedom
## Multiple R-squared: 0.728, Adjusted R-squared: 0.691
## F-statistic: 19.7 on 3 and 22 DF, p-value: 1.99e-06
# Ahora vemos una tabla de ANOVA con la función << anova() >> para
# comprobar que las sumas de cuadrados que calculamos son muy similares en
# este caso a las que calcula la función << lm >> . No son idénticas por
# el criterio de los redondeo que se utilizan en los diferentes cálculos.
anova(M)
## Analysis of Variance Table
##
## Response: FLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## LUZ 1 0.000 0.000 0.00 0.97
## SUELO 1 1.480 1.480 58.67 1.2e-07 ***
## LUZ:SUELO 1 0.008 0.008 0.34 0.57
## Residuals 22 0.555 0.025
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
rbind(LZ, SL, INT, RES)
## [,1]
## LZ 3.462e-05
## SL 1.470e+00
## INT 1.858e-02
## RES 5.550e-01
# Del predictor lineal en un ANOVA ---------------------
# El predicador lineal de un análisis de varianza tienen la misma
# estructura el de una regresión. A diferencia con la regresión el
# predicador lineal del ANOVA contiene la información necesaria para
# calcular los promedios de cualquier combinación entre los factores o
# niveles de factor que están incluidos en el análisis de varianza. El
# intercepto en el predictor lineal del ANOVA es el promedio de una de las
# combinaciones entre los factores y el resto de las betas en el
# predicador lineal del ANOVA son diferencias relativas a ese intercepto
# que nos servirán para calcular el promedio en cualquiera de las otras
# combinaciones. Estas diferencias sólo tienen sentido cuando se usan
# junto con el valor del intercepto.
# Especificando que deseo recuperar exclusivamente la tabla donde se
# encuentran los coeficientes del modelo poder después extraer la columna
# donde están los valores de las betas, los valores del predictor línea
summary(M)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.040000 0.05615 0.7123 4.837e-01
## LUZPL -0.008571 0.08220 -0.1043 9.179e-01
## SUELOR 0.522000 0.09054 5.7651 8.452e-06
## LUZPL:SUELOR -0.073429 0.12652 -0.5804 5.676e-01
# Almacenaremos las betas del predicador lineal en el objeto << PRE_L >>
# extrayendo la primera columna de la tabla recién vimos.
PRE_L <- summary(M)$coefficients[, 1]
# Es muy importante aprender a leer el predictor línea, esto es aprender a
# leer las etiquetas que tienen las diferentes betas dentro del predictor
# lineal
PRE_L
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040000 -0.008571 0.522000 -0.073429
# La primer beta como ya lo habíamos dicho es el intercepto y sabemos que
# es el promedio de una de las combinaciones de los factores, pero como
# todos los demás valores son relativos a este intercepto siempre vamos a
# considerar su valor en cualquier cálculo de medias.
# La forma más fácil de hacer los cálculos con el proyector lineal es
# generar un vector de unos y ceros que se multiplicará por el predictor
# lineal. Si la etiqueta de una de las becas no aparece en la combinación
# para la cual nosotros deseamos calcular el promedio entonces colocaremos
# un cero en este vector, y si la etiqueta de la beta si aparece en la
# combinación para la cual estamos calculando el promedio entonces
# colocaremos un uno.
# Construyamos este vector para calcular el promedio de la combinación de
# mucha luz en suelo pobre << ML_P >>. La primera beta que es el
# intercepto ya mencionamos que si la vamos a requerir siempre va ser
# requerida. La segunda beta está etiquetada para poca luz << LUZPL >> y
# nosotros queremos calcular la media para la combinación que incluye
# mucha luz << ML >> por lo tanto para esta segunda beta colocaremos un
# cero.La tercera beta está etiquetada << SUELOR >> suelo rico, también
# colocaremos cero ya nosotros queremos calcular la media para la
# combinación que incluye suelo pobre. Finalmente, la cuarta beta está
# etiquetada para una interacción <<LUZPL:SUELOR>>, poca luz en suelo
# rico. Esta combinación tampoco corresponde con la combinación para la
# cual nosotros estamos calculando la media, por lo tanto también
# colocaremos un cero. Así tenemos el vector para ponderarlas betas del
# predictor lineal para calcular la media en la combinación mucha luz y
# suelo pobre será 1,0,0,0.
ML_P <- c(1, 0, 0, 0)
# Hola solo tenemos que multiplicar el predictor lineal por este vector de
# ponderación y sumar los valores resultantes para obtener el promedio que
# deseamos. Lo cual nos da como resultado total sólo el valor del
# intercepto ya que las otras tres betas se multiplicaron por cero.
PRE_L * ML_P
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.04 0.00 0.00 0.00
sum(PRE_L * ML_P)
## [1] 0.04
# Podemos comprobar que este estimado corresponde en este caso exactamente
# al valor promedio estimado de los datos
tapply(FLORES, list(LUZ, SUELO), mean)
## P R
## ML 0.04000 0.562
## PL 0.03143 0.480
# Ahora calculemos este vector de ponderación para las otras
# combinaciones. En el caso de poca luz y suelo pobre el vector de
# ponderación será: uno para el intercepto, uno para la beta que está
# etiquetada como poca luz, cero para la beta que está etiquetada para
# suelo rico y cero para la beta que está etiquetada poca luz en suelo
# rico. No importa que en esta última beta parte de su etiqueta coincida
# con la combinación para la cual nosotros queremos calcular el promedio,
# debemos colocar cero porque la etiqueta completa no coincide con la
# combinación para la cual nosotros estamos calculando promedio.
PL_P <- c(1, 1, 0, 0)
# Multiplicamos este vector por el predictor lineal y sumamos para obtener
# el valor estimado del promedio para la combinación poca luz en suelo
# pobre
PRE_L * PL_P
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040000 -0.008571 0.000000 0.000000
sum(PRE_L * PL_P)
## [1] 0.03143
# Nuevamente comprobamos que este estimado se corresponde con el valor
# promedio estimado a partir de los datos
tapply(FLORES, list(LUZ, SUELO), mean)
## P R
## ML 0.04000 0.562
## PL 0.03143 0.480
# Aquí están los vectores de ponderación para las otras combinaciones de
# mucha luz en suelo rico y poca luz en suelo
ML_R <- c(1, 0, 1, 0)
PL_R <- c(1, 1, 1, 1)
# Multiplicamos estos vectores por el predictor lineal y sumamos para
# obtener el estimado del valor promedio para cada una de las dos
# combinaciones restantes
PRE_L * ML_R
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040 0.000 0.522 0.000
sum(PRE_L * ML_R)
## [1] 0.562
PRE_L * PL_R
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040000 -0.008571 0.522000 -0.073429
sum(PRE_L * PL_R)
## [1] 0.48
# Y estos estimados también se corresponden con los valores promedio
# estimados a partir de los datos.
tapply(FLORES, list(LUZ, SUELO), mean)
## P R
## ML 0.04000 0.562
## PL 0.03143 0.480
# Ahora bien ya teniendo los vectores de ponderación para cada una de las
# combinaciones podemos usarlos directamente para calcular el vector de
# ponderación que se requerirá por ejemplo para el nivel mucha luz
# independientemente del tipo de suelo. Esto es sencillamente el promedio
# de los dos vectores que incluyen a mucha luz, mucha luz en suelo rico y
# mucha luz en suelo pobre
ML <- (ML_P + ML_R)/2
# Y podemos observar que el vector de ponderación para el nivel mucha luz
# corresponde a uno en el intercepto, cero para la beta que está
# etiquetada poca luz, 1/2 para la beta etiqueta suelo rico y cero para la
# beta con la etiqueta de la combinación poca luz en suelo rico
ML
## [1] 1.0 0.0 0.5 0.0
# Multiplicamos el predictor lineal por el vector de ponderación de mucha
# luz y sumamos para obtener el estimado del valor promedio del peso de
# las flores de las plantas que crecieron en mucha luz independientemente
# del tipo de suelo.
PRE_L * ML
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040 0.000 0.261 0.000
sum(PRE_L * ML)
## [1] 0.301
# Comprobamos que este nuevo estimado para el nivel de mucha luz se
# corresponda con el valor calculado de los datos. Y notamos que realmente
# no hay una correspondencia, nosotros estimamos un valor de 0.301 cuando
# de los datos calculamos un valores 0.24. La razón principal de este
# desvío es que el modelo ajustado no cumple con los supuestos de
# normalidad y homocedasticidad, y en consecuencia los valores de las
# betas están sesgados.
tapply(FLORES, LUZ, mean)
## ML PL
## 0.2408 0.2385
# Con la misma lógica armamos los vectores para al nivel poca luz, y de
# los dos niveles de suelo: suelo pobre y suelo rico promediando los
# vectores iniciales de las combinaciones que correspondan.
PL <- (PL_P + PL_R)/2
SP <- (PL_P + ML_P)/2
SR <- (PL_R + ML_R)/2
# Multiplicamos estos vectores por el predictor lineal y sumamos para
# obtener el valor promedio estimado para cada uno de los niveles de estos
# factores
PRE_L * PL
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040000 -0.008571 0.261000 -0.036714
sum(PRE_L * PL)
## [1] 0.2557
# Estos valores también notaremos que no e corresponde con lo que podemos
# calcular directamente de los datos.
tapply(FLORES, LUZ, mean)
## ML PL
## 0.2408 0.2385
# A los niveles de suelo se observa también una incongruencia entre lo
# estimado a partir del predictor lineal y lo calculado a partir de los
# datos.
PRE_L * SP
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040000 -0.004286 0.000000 0.000000
sum(PRE_L * SP)
## [1] 0.03571
PRE_L * SR
## (Intercept) LUZPL SUELOR LUZPL:SUELOR
## 0.040000 -0.004286 0.522000 -0.036714
sum(PRE_L * SR)
## [1] 0.521
tapply(FLORES, SUELO, mean)
## P R
## 0.0360 0.5173
# Regresando a los supuestos del modelo, los residuales deben cumplir con
# el criterio de normalidad y homocedasticidad y esto lo debimos haber
# revisado antes y después haber realizado todos los cálculos que hicimos.
# Para revisar los supuestos del modelo en cuanto normalidad y
# homocedasticidad lo más fácil es dividir la pantalla gráfica en cuatro
# paneles dos filas y los columnas, y hacer un gráfico del modelo <<
# plot() >>
par(mfrow = c(2, 2))
plot(M)
# El panel superior izquierdo y el panel inferior izquierdo representan el
# criterio de la homocedasticidad. La distribución de los residuales en
# estos los gráficos debería de ser tal que los ámbitos de variación en
# esas pequeñas líneas verticales que se dibujan sean de la misma amplitud
# a todo lo largo de los valores ajustados. No es relevante que están
# separadas se encuentran esas líneas de puntos, la separación simplemente
# indicar diferencias entre las medias. Lo que es crítico es el ámbito de
# variación de cada una de estas líneas que debe ser similar al de las
# otras líneas.
# En este ejemplo claramente vemos que no se cumple con el principio de
# homocedasticidad. En el lado izquierdo del panel los residuales dibujan
# una pequeña línea vertical muy estrecha mientras los residuales en la
# parte derecha del gráfico dibujan una línea vertical mucho más amplia.
# En el panel superior derecho podemos revisar el supuesto de la
# normalidad. Aquí la distribución de puntos debe ser tal que se
# sobrepongan en la mayor parte sobre la línea dibujada que representa una
# distribución normal teórica. En los extremos los residuales siempre se
# saldrán de la línea eso no es de gran relevancia lo importante es que la
# parte central los puntos se encuentren sobre la línea o alrededor de
# ella de forma aleatoria.
# También podemos revisar la normalidad de los residuales utilizando ya
# sea la prueba Shapiro << shapiro.test() o la prueba de Kolmogrov-Smirnov
# << ks.test() >>. La prueba de Shapiro es un tanto cuanto exigente en
# cuanto a la normalida por lo cual se prefiere la exploración gráfica o
# la prueba Kolmogrov-Smirnov.
# Para la prueba de Shapiro sólo requerimos extraer los residuales a
# partir del modelo lo cual hacemos con el nombre del modelo $residuals <<
# M2$residuals >>
shapiro.test(M$residuals)
##
## Shapiro-Wilk normality test
##
## data: M$residuals
## W = 0.7463, p-value = 2.445e-05
# La prueba de Kolmogrov-Smirnov es más genérica, y sirve para comparar la
# distribución de dos grupos de datos o de un grupo de datos respecto a
# una distribución teórica. Por lo tanto debemos especificar en esta
# función los residuales y contra qué distribución queremos compararlos,
# en este caso contra una distribución probabilística normal.
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.398, p-value = 0.0005298
## alternative hypothesis: two-sided
# Como vemos es modelo ajustado no cumple con los supuestos de normalidad
# y homocedasticidad. Debemos transformar los datos para ajustar el modelo
# adecuadamente y entonces poder hacer predicciones más certeras que no
# estén sesgadas
# TRANSFORMANDO--------------
# La transformación que utilizaremos será el logaritmo. Como los datos
# están en gramos y son puras fracciones menores a uno primero
# multiplicaré los valores de los del peso de flores por 1000 para tenerlo
# en miligramos y posteriormente sacar el logaritmo.
LOGFLORES <- log(FLORES * 1000)
# Ahora ajustón modelo con nuestros datos transformados, revisamos los
# puestos del modelo y podemos verificar que se cumplen ambos. Ahora
# existe normalidad de los residuales y en los paneles de homocedasticidad
# se observan lineas muchos mas homogéneas en amplitud.
M2 <- lm(LOGFLORES ~ LUZ * SUELO)
par(mfrow = c(2, 2))
plot(M2)
shapiro.test(M2$residuals)
##
## Shapiro-Wilk normality test
##
## data: M2$residuals
## W = 0.9693, p-value = 0.6064
ks.test(M2$residuals, "pnorm")
## Warning: ties should not be present for the Kolmogorov-Smirnov test
##
## One-sample Kolmogorov-Smirnov test
##
## data: M2$residuals
## D = 0.1707, p-value = 0.4348
## alternative hypothesis: two-sided
# Ahora todos los estimados a partir del predictor lineal coincidirán
# mejor con los estimados directamente de los datos para los valores
# promedio. Primero extraemos en un vector el predictor lineal, es decir
# los valores de las betas.
PRE_L <- summary(M2)$coefficients[, 1]
sum(PRE_L * ML_P)
## [1] 3.346
sum(PRE_L * PL_P)
## [1] 3.284
sum(PRE_L * ML_R)
## [1] 6.261
sum(PRE_L * PL_R)
## [1] 6.088
tapply(LOGFLORES, list(LUZ, SUELO), mean)
## P R
## ML 3.346 6.261
## PL 3.284 6.088
sum(PRE_L * ML)
## [1] 4.804
sum(PRE_L * PL)
## [1] 4.686
tapply(LOGFLORES, LUZ, mean)
## ML PL
## 4.467 4.578
sum(PRE_L * SP)
## [1] 3.315
sum(PRE_L * SR)
## [1] 6.175
tapply(LOGFLORES, SUELO, mean)
## P R
## 3.317 6.167
tapply(LOGFLORES, list(LUZ, SUELO), mean, na.rm = TRUE)
## P R
## ML 3.346 6.261
## PL 3.284 6.088
# Graficando los resultados--------------- Para graficar los dos
# resultados tenemos dos o tres opciones gráficas de caja, gráficas de
# barras con errores o de líneas con errores
par(mfrow = c(1, 1))
boxplot(LOGFLORES ~ SUELO)
MEDIAS <- tapply(LOGFLORES, SUELO, mean)
EE <- tapply(LOGFLORES, SUELO, ee)
barplot.error(MEDIAS, EE)
lineplot.error(MEDIAS, EE, cex.pch = 0.5, lwd.x = 2, ylim = c(2, 7))
abline(h = 2, lwd = 2)
MEDIAS <- tapply(LOGFLORES, list(SUELO, LUZ), mean)
EE <- tapply(LOGFLORES, list(SUELO, LUZ), ee)
lineplot.error(MEDIAS, EE, cex.pch = 0.9, lwd.x = 2, ylim = c(2, 7), pch = c(19,
15), col.b = "grey")
legend(1, 7, c("ML", "PL"), pch = c(19, 15))
