# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# ANOVAS DE PARCELAS DIVIDIDAS, split-plot
# Vamos a utilizar las bases de datos << splityield,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"
DATOS2 <- subset(DATOS, LUZ == "ML" & 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
detach(DATOS)
rm(DATOS)
# DEVIANZA CON ERROR GAUSIANO ------- Los cálculos de la devianza en los
# GLM con error normal y función de liga de identidad son idénticos que
# los cálculos de la suma de cuadrados
# Asi en el modelo numlo la devianza residual es igual a la suma de
# cuadrados de la media
par(mar = c(4, 4, 3, 1))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "", ylab = "", main = "Gausiano",
cex.main = 4)
text(0.5, 0.5, expression(paste(Sigma, "(", Y, " - ", mu, ")")^2), cex = 4)
M <- glm(FLORES ~ 1)
anova(M)
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 9.78
sum((FLORES - mean(FLORES))^2)
## [1] 9.776
# Y de la misma forma se calcula la devianza para en factor, en este caso
# el factor SUELO
M <- glm(FLORES ~ SUELO)
anova(M)
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 9.78
## SUELO 1 6.84 111 2.94
sum((FLORES[SUELO == "P"] - mean(FLORES[SUELO == "P"]))^2) + sum((FLORES[SUELO ==
"R"] - mean(FLORES[SUELO == "R"]))^2)
## [1] 2.938
# DEVIANZA CON ERROR GAMMA ------ Aquí los cálculos son diferentes se
# utiliza como mediad de dispersión
par(mar = c(4, 4, 3, 1))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "", ylab = "", main = "Gamma",
cex.main = 4)
text(0.5, 0.5, expression(paste(2, Sigma, frac(Y, mu), " - ", "ln", frac(Y,
mu), " - 1")), cex = 4)
M <- glm(FLORES ~ 1, Gamma)
anova(M)
## Analysis of Deviance Table
##
## Model: Gamma, link: inverse
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 127
2 * sum((FLORES/mean(FLORES)) - log(FLORES/mean(FLORES)) - 1)
## [1] 126.6
# Y de la misma forma se calcula la devianza para el factor, en este caso
# el factor SUELO
M <- glm(FLORES ~ SUELO, Gamma)
anova(M)
## Analysis of Deviance Table
##
## Model: Gamma, link: inverse
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 126.6
## SUELO 1 85.1 111 41.5
2 * sum((FLORES[SUELO == "P"]/mean(FLORES[SUELO == "P"])) - log(FLORES[SUELO ==
"P"]/mean(FLORES[SUELO == "P"])) - 1) + 2 * sum((FLORES[SUELO == "R"]/mean(FLORES[SUELO ==
"R"])) - log(FLORES[SUELO == "R"]/mean(FLORES[SUELO == "R"])) - 1)
## [1] 41.52
# DEVIANZA CON ERROR GAUSIANO INVERSO ----- Aquí los cálculos son
# diferentes se utiliza como medida de dispersión
par(mar = c(4, 4, 3, 1))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "", ylab = "", main = "Gausiano inverso",
cex.main = 3)
text(0.5, 0.5, expression(paste(Sigma, frac(paste("(", Y, " - ", mu, ")")^2,
paste(mu^2, Y)))), cex = 4)
M <- glm(FLORES ~ 1, inverse.gaussian)
anova(M)
## Analysis of Deviance Table
##
## Model: inverse.gaussian, link: 1/mu^2
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 973
sum((FLORES - mean(FLORES))^2/(mean(FLORES)^2 * FLORES))
## [1] 973.3
# Y de la misma forma se calcula la devianza para el factor, en este caso
# el factor SUELO
M <- glm(FLORES ~ SUELO, inverse.gaussian)
anova(M)
## Analysis of Deviance Table
##
## Model: inverse.gaussian, link: 1/mu^2
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 973
## SUELO 1 372 111 601
sum((FLORES[SUELO == "P"] - mean(FLORES[SUELO == "P"]))^2/(mean(FLORES[SUELO ==
"P"])^2 * FLORES[SUELO == "P"])) + sum((FLORES[SUELO == "R"] - mean(FLORES[SUELO ==
"R"]))^2/(mean(FLORES[SUELO == "R"])^2 * FLORES[SUELO == "R"]))
## [1] 601
# DEVIANZA CON ERROR POISSON ------ Aquí los cálculos son diferentes se
# utiliza como medida de dispersión
par(mar = c(4, 4, 3, 1))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "", ylab = "", main = "Poisson",
cex.main = 3)
text(0.5, 0.5, expression(paste(2, Sigma, "y ln", frac(y, mu))), cex = 4)
M <- glm(NUM_FLORES ~ 1, poisson)
anova(M)
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: NUM_FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 464
2 * sum(NUM_FLORES * log(NUM_FLORES/mean(NUM_FLORES)))
## [1] 464.1
# Y de la misma forma se calcula la devianza para el factor, en este caso
# el factor SUELO
M <- glm(NUM_FLORES ~ SUELO, poisson)
anova(M)
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: NUM_FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 112 464
## SUELO 1 396 111 69
2 * sum(NUM_FLORES[SUELO == "P"] * log(NUM_FLORES[SUELO == "P"]/mean(NUM_FLORES[SUELO ==
"P"]))) + 2 * sum(NUM_FLORES[SUELO == "R"] * log(NUM_FLORES[SUELO == "R"]/mean(NUM_FLORES[SUELO ==
"R"])))
## [1] 68.54
# DEVIANZA CON ERROR BINOMIAL ------ Aquí los cálculos son diferentes se
# utiliza como medida de dispersión
par(mar = c(4, 4, 3, 1))
plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "", ylab = "", main = "Binomial",
cex.main = 3)
text(0.5, 0.5, expression(paste(2, Sigma, "y ln", frac(y, mu), " + n-y ln",
frac(n - y, n - mu))), cex = 3.8)
Y <- cbind(FRUT_MADUROS, FRUT_INMADUROS)
ESTOS <- which(Y[, 1] > 0 & Y[, 2] > 0)
Y <- Y[ESTOS, ]
N <- rowSums(Y)
P <- sum(Y[, 1])/sum(N)
MU <- N * P
YFM <- Y[, 1]
YFIM <- Y[, 2]
M <- glm(Y ~ 1, binomial)
anova(M)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 77 92
2 * sum(YFM * log(YFM/MU) + YFIM * log(YFIM/(N - MU)))
## [1] 92.01
# Y de la misma forma se calcula la devianza para el factor, en este caso
# el factor SUELO
SL <- SUELO[ESTOS]
M <- glm(Y ~ SL, binomial)
ESTOS2 <- which(SL == "P")
YSP <- Y[ESTOS2, ]
NSP <- rowSums(Y[ESTOS2, ])
PSP <- sum(YSP[, 1])/sum(NSP)
MUSP <- NSP * PSP
YFMSP <- Y[ESTOS2, 1]
YFIMSP <- Y[ESTOS2, 2]
ESTOS3 <- which(SL == "R")
YSR <- Y[ESTOS3, ]
NSR <- rowSums(Y[ESTOS3, ])
PSR <- sum(YSR[, 1])/sum(NSR)
MUSR <- NSR * PSR
YFMSR <- Y[ESTOS3, 1]
YFIMSR <- Y[ESTOS3, 2]
anova(M)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 77 92
## SL 1 3.02 76 89
2 * sum(YFMSP * log(YFMSP/MUSP) + YFIMSP * log(YFIMSP/(NSP - MUSP))) + 2 * sum(YFMSR *
log(YFMSR/MUSR) + YFIMSR * log(YFIMSR/(NSR - MUSR)))
## [1] 88.99
#### GLM ------
DATOS <- read.table("~/desktop/curso R 2012/TOLOACHE.csv", header = TRUE, sep = ",")
attach(DATOS)
## The following object(s) are masked from 'DATOS2':
##
## 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
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 == "R" & 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
## The following object(s) are masked from 'DATOS2 (position 4)':
##
## 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
detach(DATOS)
rm(DATOS)
# Primero ajustamos el modelo GLM con error gausiano
M <- glm(FLORES ~ LUZ/(HERBIVORIA * MICORRIZA))
# Note la falta de homocedasticidad
par(mfrow = c(2, 2))
plot(M)
# Ahora el GLM con error Gamma para corregir la falta de homocedasticidad
M2 <- glm(FLORES ~ LUZ/(HERBIVORIA * MICORRIZA), Gamma(link = "inverse"))
# Note la mejora en la homocedasticidad
par(mfrow = c(1, 1))
par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(M2)
# Verifiamos la normalidad de los residuales
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.1238, p-value = 0.07074
## alternative hypothesis: two-sided
shapiro.test(M2$residuals)
##
## Shapiro-Wilk normality test
##
## data: M2$residuals
## W = 0.9461, p-value = 0.000242
# Simplificamos el modelo via << step() >>
MSTEP <- step(M)
## Start: AIC=-19.99
## FLORES ~ LUZ/(HERBIVORIA * MICORRIZA)
##
## Df Deviance AIC
## - LUZ:HERBIVORIA:MICORRIZA 2 4.58 -22.1
## <none> 4.50 -20.0
##
## Step: AIC=-22.11
## FLORES ~ LUZ + LUZ:HERBIVORIA + LUZ:MICORRIZA
##
## Df Deviance AIC
## <none> 4.58 -22.1
## - LUZ:HERBIVORIA 2 4.89 -18.9
## - LUZ:MICORRIZA 2 5.10 -14.5
anova(MSTEP)
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 108 5.85
## LUZ 1 0.450 107 5.40
## LUZ:HERBIVORIA 2 0.303 105 5.10
## LUZ:MICORRIZA 2 0.516 103 4.58
# Simplificamos el modelo por ANOVA entre modelos
M3 <- update(M2, ~. - LUZ:HERBIVORIA:MICORRIZA)
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: FLORES ~ LUZ + LUZ:HERBIVORIA + LUZ:MICORRIZA
## Model 2: FLORES ~ LUZ/(HERBIVORIA * MICORRIZA)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 103 19.8
## 2 101 19.7 2 0.0995 0.74
M4 <- update(M3, ~. - LUZ:HERBIVORIA)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: FLORES ~ LUZ + LUZ:MICORRIZA
## Model 2: FLORES ~ LUZ + LUZ:HERBIVORIA + LUZ:MICORRIZA
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 105 21.2
## 2 103 19.8 2 1.45 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M4 <- update(M3, ~. - LUZ:MICORRIZA)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: FLORES ~ LUZ + LUZ:HERBIVORIA
## Model 2: FLORES ~ LUZ + LUZ:HERBIVORIA + LUZ:MICORRIZA
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 105 21.6
## 2 103 19.8 2 1.83 0.0036 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(M3)$coefficients
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.54814 0.1506 10.2772 1.788e-17
## LUZPL 1.30539 0.2801 4.6607 9.466e-06
## LUZML:HERBIVORIASH -0.07707 0.1854 -0.4158 6.784e-01
## LUZPL:HERBIVORIASH -0.70079 0.2406 -2.9131 4.389e-03
## LUZML:MICORRIZASIN 0.46833 0.1940 2.4142 1.754e-02
## LUZPL:MICORRIZASIN -0.53693 0.2366 -2.2690 2.535e-02
## Los vectores de ponderación, parecen mas complicados pero siguen la
## misma lógica que vimos anteriormente. Note los valores de 0.5, ya que
## esta interacción << HERBIVORÍA:MICORRIZA >> está anidada en << LUZ >>
ML_CONH <- c(1, 0, 0, 0, 0.5, 0)
ML_SINH <- c(1, 0, 1, 0, 0.5, 0)
DIF1 <- ML_CONH - ML_SINH
PL_CONH <- c(1, 1, 0, 0, 0, 0.5)
PL_SINH <- c(1, 1, 0, 1, 0, 0.5)
DIF2 <- PL_CONH - PL_SINH
ML_CONM <- c(1, 0, 0.5, 0, 0, 0)
ML_SINM <- c(1, 0, 0.5, 0, 1, 0)
DIF3 <- ML_CONM - ML_SINM
PL_CONM <- c(1, 1, 0, 0.5, 0, 0)
PL_SINM <- c(1, 1, 0, 0.5, 0, 1)
DIF4 <- PL_CONM - PL_SINM
library(gmodels)
estimable(M3, rbind(DIF1, DIF2, DIF3, DIF4))
## Estimate Std. Error t value DF Pr(>|t|)
## DIF1 0.07707 0.1854 0.4158 103 0.678423
## DIF2 0.70079 0.2406 2.9131 103 0.004389
## DIF3 -0.46833 0.1940 -2.4142 103 0.017536
## DIF4 0.53693 0.2366 2.2690 103 0.025352
# Para graficar extraemos las medias HERBIVORIA
MEDIAS <- estimable(M3, rbind(ML_CONH, ML_SINH, PL_CONH, PL_SINH))[, 1]
# Las acomodamos en una matriz ene l que cada fila será un grupo de barras
# en el gráfico
MEDIAS <- matrix(MEDIAS, 2, 2, byrow = T)
# Hacemos los mismo con los errores
ERROR <- estimable(M3, rbind(ML_CONH, ML_SINH, PL_CONH, PL_SINH))[, 2]
ERROR <- matrix(ERROR, 2, 2, byrow = T)
# Y graficamos
source("~/desktop/curso R 2012/funciones.r")
par(mfrow = c(1, 1), mar = c(4, 4, 3, 1))
barplot.error(MEDIAS, ERROR, func = "inverse", ylim = c(0, 0.8), col = c("grey30",
"grey80"), names.arg = c("Mucha luz", "Poca Luz"), ylab = "Flores (g)",
xlab = "", main = "Herbivoría")
legend(2, 0.8, c("Con herbivoría", "Sin herbivoría"), fill = c("grey30", "grey80"))
# Para el gráfico de las micorrizas se sigue exactamente los mismos pasos
# MICORRIZA
MEDIAS <- estimable(M3, rbind(ML_CONM, ML_SINM, PL_CONM, PL_SINM))[, 1]
MEDIAS <- matrix(MEDIAS, 2, 2, byrow = T)
ERROR <- estimable(M3, rbind(ML_CONM, ML_SINM, PL_CONM, PL_SINM))[, 2]
ERROR <- matrix(ERROR, 2, 2, byrow = T)
par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
barplot.error(MEDIAS, ERROR, func = "inverse", ylim = c(0, 0.8), col = c("grey30",
"grey80"), names.arg = c("Mucha luz", "Poca Luz"), ylab = "Flores (g)",
xlab = "", main = "Micorriza")
legend(7, 0.8, c("Con micorriza", "Sin micorriza"), fill = c("grey30", "grey80"))
