# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
# Regresión (cont.) ----
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"
# Ahora hagamos modelos de regresión por separado para cada combinación de
# << LUZ >>, << SUELO >> y << HERBIVORIA >>. Son en total 8 combinaciones
# diferentes. Para cada una vamos a calcular su modelo de regresión y
# graficarlo junto con el intervalo de confianza
# Lo primero que debemos hacer es crear subconjuntos de datos para una de
# las combinaciones. Empecemos con la combinación << LUZ=='ML' >>, <<
# SUELO=='P' >> y << HERBIVORIA=='CH'
# MUCHA LUZ, SUELO POBRE, CON HERBIVORÍA
NO.FLORES <- NUM_FLORES[LUZ == "ML" & SUELO == "P" & HERBIVORIA == "CH"]
PEOS.FLORES <- FLORES[LUZ == "ML" & SUELO == "P" & HERBIVORIA == "CH"]
M <- lm(PEOS.FLORES ~ NO.FLORES)
summary(M)
##
## Call:
## lm(formula = PEOS.FLORES ~ NO.FLORES)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0540 -0.0347 -0.0213 0.0174 0.2374
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.03677 0.02744 1.34 0.191
## NO.FLORES 0.01862 0.00911 2.04 0.051 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0631 on 28 degrees of freedom
## Multiple R-squared: 0.13, Adjusted R-squared: 0.0987
## F-statistic: 4.17 on 1 and 28 DF, p-value: 0.0506
anova(M)
## Analysis of Variance Table
##
## Response: PEOS.FLORES
## Df Sum Sq Mean Sq F value Pr(>F)
## NO.FLORES 1 0.0166 0.01660 4.17 0.051 .
## Residuals 28 0.1113 0.00398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(PEOS.FLORES ~ NO.FLORES)
abline(M, col = "red")
SEQ <- seq(1, 6, length = 30)
ND <- data.frame(NO.FLORES = SEQ)
PRED <- predict(M, ND, se = TRUE)
T <- qt(0.025, 28)
lines(SEQ, PRED$fit - T * PRED$se.fit, lty = 2, col = "red")
lines(SEQ, PRED$fit + T * PRED$se.fit, lty = 2, col = "red")
lines(SEQ, PRED$fit - 1.96 * PRED$se.fit, lty = 2, col = "purple")
lines(SEQ, PRED$fit + 1.96 * PRED$se.fit, lty = 2, col = "purple")
abline(h = mean(PEOS.FLORES), col = "blue", lty = 3)
# Note como en este ejemplo al colocar una linea horizontal en el promedio
# de Y << PESO.FLORES >>, la hipótesis nula, esta cabe perfectamente
# dentro del intervalo de confianza lo cual es un indicio gráfico de la
# que relación es estadísticamente no significativa. Comparé esto con los
# siguientes ejemplos de las regresiones en las diferentes combinaciones
# de << LUZ >>, << SUELO >> Y << HERBIVORIA >> que se realizan mediante
# bucles anidados.
for (i in levels(LUZ)) {
for (j in levels(SUELO)) {
for (k in levels(HERBIVORIA)) {
NO.FLORES <- NUM_FLORES[LUZ == i & SUELO == j & HERBIVORIA == k]
PEOS.FLORES <- FLORES[LUZ == i & SUELO == j & HERBIVORIA == k]
M <- lm(PEOS.FLORES ~ NO.FLORES)
summary(M)
anova(M)
plot(PEOS.FLORES ~ NO.FLORES, main = paste(i, j, k))
abline(M, col = "red")
SEQ <- seq(min(NO.FLORES, na.rm = TRUE), max(NO.FLORES, na.rm = TRUE),
length = 30)
ND <- data.frame(NO.FLORES = SEQ)
PRED <- predict(M, ND, se = TRUE)
T <- qt(0.025, 28)
lines(SEQ, PRED$fit - T * PRED$se.fit, lty = 2, col = "red")
lines(SEQ, PRED$fit + T * PRED$se.fit, lty = 2, col = "red")
abline(h = mean(PEOS.FLORES, na.rm = TRUE), col = "blue", lty = 3)
}
}
}
## Regresión múltiple------
DATOS <- read.table("~/desktop/curso R 2012/TOLOACHE.csv", header = TRUE, sep = ",")
attach(DATOS)
## The following object(s) are masked from 'DATOS (position 3)':
##
## 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"
# Ahora tratemos de incluir varias variables explicativas en el modelo de
# regresión. Esto se vuelve una regresión múltiple. Consideremos el
# siguiente modelo con cinco variables explicativas
M <- lm(FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5878 -0.0444 -0.0003 0.0406 0.8259
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.039694 0.043096 0.92 0.3581
## TALLO 0.012660 0.004170 3.04 0.0027 **
## HOJAS -0.001311 0.006870 -0.19 0.8489
## RAIZ 0.042239 0.014315 2.95 0.0035 **
## AREA 0.001803 0.000328 5.50 1.1e-07 ***
## DUREZA -0.000877 0.000728 -1.21 0.2295
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.12 on 205 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.824, Adjusted R-squared: 0.819
## F-statistic: 192 on 5 and 205 DF, p-value: <2e-16
# Podemos observar que algunos de los coeficientes del predictor lineal
# son estadísticamente significativos, es decir son diferentes de cero. La
# varianza explicada es 0.8237. Con estas cinco variables explicativas
# tenemos además la posibilidad de incluir 10 interacciones en pares de
# estas variables. La notación para estos por ejemplo << TALLO:HOJAS >>
# donde los dos puntos es lo que determina que las variables que le
# flaquean son las que interactúan.
M <- lm(FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA + TALLO:HOJAS + TALLO:RAIZ +
TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ + HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA +
RAIZ:DUREZA + AREA:DUREZA)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA +
## TALLO:HOJAS + TALLO:RAIZ + TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ +
## HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA + RAIZ:DUREZA + AREA:DUREZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5757 -0.0443 -0.0057 0.0375 0.7928
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.16e-02 7.97e-02 0.77 0.44
## TALLO 2.46e-02 3.24e-02 0.76 0.45
## HOJAS -1.28e-02 4.41e-02 -0.29 0.77
## RAIZ 5.92e-02 1.04e-01 0.57 0.57
## AREA -1.33e-03 2.45e-03 -0.54 0.59
## DUREZA -4.27e-04 1.42e-03 -0.30 0.76
## TALLO:HOJAS 1.13e-03 2.01e-03 0.56 0.57
## TALLO:RAIZ -3.56e-03 3.52e-03 -1.01 0.31
## TALLO:AREA -1.47e-04 1.45e-04 -1.01 0.31
## TALLO:DUREZA 1.37e-04 3.65e-04 0.37 0.71
## HOJAS:RAIZ 2.07e-03 1.01e-02 0.21 0.84
## HOJAS:AREA 8.30e-05 2.44e-04 0.34 0.73
## HOJAS:DUREZA -3.65e-04 5.27e-04 -0.69 0.49
## RAIZ:AREA 7.10e-04 4.56e-04 1.56 0.12
## RAIZ:DUREZA -6.36e-04 1.21e-03 -0.52 0.60
## AREA:DUREZA 2.99e-05 3.59e-05 0.83 0.41
##
## Residual standard error: 0.12 on 195 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.832, Adjusted R-squared: 0.819
## F-statistic: 64.5 on 15 and 195 DF, p-value: <2e-16
# Note que hay una ligera mejora en el valor de varianza explicada, y la
# la probabilidad de los coeficientes ha cambiado.
# Además de las interacciones dobles también podemos incluir 10
# interacciones triples
M <- lm(FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA + TALLO:HOJAS + TALLO:RAIZ +
TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ + HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA +
RAIZ:DUREZA + AREA:DUREZA + TALLO:HOJAS:RAIZ + TALLO:HOJAS:AREA + TALLO:HOJAS:DUREZA +
TALLO:RAIZ:AREA + TALLO:RAIZ:DUREZA + TALLO:AREA:DUREZA + HOJAS:RAIZ:AREA +
HOJAS:RAIZ:DUREZA + HOJAS:AREA:DUREZA + RAIZ:AREA:DUREZA)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA +
## TALLO:HOJAS + TALLO:RAIZ + TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ +
## HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA + RAIZ:DUREZA + AREA:DUREZA +
## TALLO:HOJAS:RAIZ + TALLO:HOJAS:AREA + TALLO:HOJAS:DUREZA +
## TALLO:RAIZ:AREA + TALLO:RAIZ:DUREZA + TALLO:AREA:DUREZA +
## HOJAS:RAIZ:AREA + HOJAS:RAIZ:DUREZA + HOJAS:AREA:DUREZA +
## RAIZ:AREA:DUREZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5829 -0.0446 -0.0091 0.0439 0.7164
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.70e-02 1.70e-01 0.39 0.6948
## TALLO 2.05e-01 1.23e-01 1.67 0.0964 .
## HOJAS -1.98e-01 1.74e-01 -1.14 0.2564
## RAIZ -4.15e-01 3.78e-01 -1.10 0.2735
## AREA 3.93e-03 6.74e-03 0.58 0.5611
## DUREZA -1.01e-03 3.10e-03 -0.33 0.7447
## TALLO:HOJAS 2.96e-03 1.09e-02 0.27 0.7858
## TALLO:RAIZ -3.94e-02 1.66e-02 -2.38 0.0186 *
## TALLO:AREA 4.35e-04 1.14e-03 0.38 0.7025
## TALLO:DUREZA -4.88e-03 2.04e-03 -2.39 0.0181 *
## HOJAS:RAIZ 8.63e-02 4.49e-02 1.92 0.0563 .
## HOJAS:AREA -1.68e-03 2.28e-03 -0.74 0.4625
## HOJAS:DUREZA 3.49e-03 3.16e-03 1.11 0.2703
## RAIZ:AREA 5.31e-04 3.01e-03 0.18 0.8603
## RAIZ:DUREZA 1.36e-02 6.25e-03 2.18 0.0302 *
## AREA:DUREZA -6.23e-05 1.15e-04 -0.54 0.5898
## TALLO:HOJAS:RAIZ 1.33e-04 5.17e-04 0.26 0.7977
## TALLO:HOJAS:AREA -1.34e-04 6.46e-05 -2.08 0.0391 *
## TALLO:HOJAS:DUREZA 1.23e-04 1.62e-04 0.76 0.4505
## TALLO:RAIZ:AREA 3.36e-05 1.00e-04 0.34 0.7374
## TALLO:RAIZ:DUREZA 7.17e-04 2.86e-04 2.51 0.0130 *
## TALLO:AREA:DUREZA 8.76e-06 1.77e-05 0.49 0.6221
## HOJAS:RAIZ:AREA 3.90e-04 2.82e-04 1.38 0.1683
## HOJAS:RAIZ:DUREZA -2.29e-03 7.53e-04 -3.04 0.0027 **
## HOJAS:AREA:DUREZA 3.99e-05 3.36e-05 1.19 0.2365
## RAIZ:AREA:DUREZA -7.80e-05 4.06e-05 -1.92 0.0559 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.116 on 185 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.85, Adjusted R-squared: 0.83
## F-statistic: 42.1 on 25 and 185 DF, p-value: <2e-16
# Esto causa otro ligero aumento en la varianza explicada. Y todavía
# podemos agregar 5 interacciones de cuarto orden
M <- lm(FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA + TALLO:HOJAS + TALLO:RAIZ +
TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ + HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA +
RAIZ:DUREZA + AREA:DUREZA + TALLO:HOJAS:RAIZ + TALLO:HOJAS:AREA + TALLO:HOJAS:DUREZA +
TALLO:RAIZ:AREA + TALLO:RAIZ:DUREZA + TALLO:AREA:DUREZA + HOJAS:RAIZ:AREA +
HOJAS:RAIZ:DUREZA + HOJAS:AREA:DUREZA + RAIZ:AREA:DUREZA + TALLO:HOJAS:RAIZ:AREA +
TALLO:HOJAS:RAIZ:DUREZA + TALLO:HOJAS:AREA:DUREZA + TALLO:RAIZ:AREA:DUREZA +
HOJAS:RAIZ:AREA:DUREZA)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA +
## TALLO:HOJAS + TALLO:RAIZ + TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ +
## HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA + RAIZ:DUREZA + AREA:DUREZA +
## TALLO:HOJAS:RAIZ + TALLO:HOJAS:AREA + TALLO:HOJAS:DUREZA +
## TALLO:RAIZ:AREA + TALLO:RAIZ:DUREZA + TALLO:AREA:DUREZA +
## HOJAS:RAIZ:AREA + HOJAS:RAIZ:DUREZA + HOJAS:AREA:DUREZA +
## RAIZ:AREA:DUREZA + TALLO:HOJAS:RAIZ:AREA + TALLO:HOJAS:RAIZ:DUREZA +
## TALLO:HOJAS:AREA:DUREZA + TALLO:RAIZ:AREA:DUREZA + HOJAS:RAIZ:AREA:DUREZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5778 -0.0452 -0.0029 0.0405 0.7203
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.42e-01 2.07e-01 0.69 0.49
## TALLO 1.32e-01 2.56e-01 0.51 0.61
## HOJAS -1.99e-03 2.92e-01 -0.01 0.99
## RAIZ -4.23e-01 7.52e-01 -0.56 0.57
## AREA -1.83e-03 8.50e-03 -0.21 0.83
## DUREZA -2.18e-03 3.77e-03 -0.58 0.56
## TALLO:HOJAS -5.46e-04 2.47e-02 -0.02 0.98
## TALLO:RAIZ -1.68e-02 4.23e-02 -0.40 0.69
## TALLO:AREA 1.77e-03 2.89e-03 0.61 0.54
## TALLO:DUREZA -3.03e-03 5.15e-03 -0.59 0.56
## HOJAS:RAIZ 3.61e-03 1.03e-01 0.04 0.97
## HOJAS:AREA -3.33e-03 4.49e-03 -0.74 0.46
## HOJAS:DUREZA -1.73e-03 6.35e-03 -0.27 0.79
## RAIZ:AREA 4.29e-03 9.84e-03 0.44 0.66
## RAIZ:DUREZA 1.43e-02 1.40e-02 1.02 0.31
## AREA:DUREZA 4.18e-05 1.50e-04 0.28 0.78
## TALLO:HOJAS:RAIZ 2.26e-03 2.34e-03 0.97 0.34
## TALLO:HOJAS:AREA -1.51e-04 2.84e-04 -0.53 0.59
## TALLO:HOJAS:DUREZA 2.35e-04 4.56e-04 0.51 0.61
## TALLO:RAIZ:AREA -5.01e-04 4.90e-04 -1.02 0.31
## TALLO:RAIZ:DUREZA 1.12e-04 9.55e-04 0.12 0.91
## TALLO:AREA:DUREZA -2.71e-05 6.20e-05 -0.44 0.66
## HOJAS:RAIZ:AREA 1.04e-03 1.33e-03 0.78 0.44
## HOJAS:RAIZ:DUREZA -3.91e-04 2.14e-03 -0.18 0.86
## HOJAS:AREA:DUREZA 9.84e-05 9.69e-05 1.02 0.31
## RAIZ:AREA:DUREZA -1.57e-04 1.93e-04 -0.81 0.42
## TALLO:HOJAS:RAIZ:AREA 9.66e-07 1.17e-05 0.08 0.93
## TALLO:HOJAS:RAIZ:DUREZA -4.54e-05 4.38e-05 -1.04 0.30
## TALLO:HOJAS:AREA:DUREZA -1.34e-07 5.13e-06 -0.03 0.98
## TALLO:RAIZ:AREA:DUREZA 1.36e-05 1.27e-05 1.07 0.28
## HOJAS:RAIZ:AREA:DUREZA -1.81e-05 2.93e-05 -0.62 0.54
##
## Residual standard error: 0.117 on 180 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.853, Adjusted R-squared: 0.828
## F-statistic: 34.8 on 30 and 180 DF, p-value: <2e-16
# Y vuelve a aumentar la varianza explicada, y aún falta incluir la
# interacción de quinto orden!
M <- lm(FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA + TALLO:HOJAS + TALLO:RAIZ +
TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ + HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA +
RAIZ:DUREZA + AREA:DUREZA + TALLO:HOJAS:RAIZ + TALLO:HOJAS:AREA + TALLO:HOJAS:DUREZA +
TALLO:RAIZ:AREA + TALLO:RAIZ:DUREZA + TALLO:AREA:DUREZA + HOJAS:RAIZ:AREA +
HOJAS:RAIZ:DUREZA + HOJAS:AREA:DUREZA + RAIZ:AREA:DUREZA + TALLO:HOJAS:RAIZ:AREA +
TALLO:HOJAS:RAIZ:DUREZA + TALLO:HOJAS:AREA:DUREZA + TALLO:RAIZ:AREA:DUREZA +
HOJAS:RAIZ:AREA:DUREZA + TALLO:HOJAS:RAIZ:AREA:DUREZA)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO + HOJAS + RAIZ + AREA + DUREZA +
## TALLO:HOJAS + TALLO:RAIZ + TALLO:AREA + TALLO:DUREZA + HOJAS:RAIZ +
## HOJAS:AREA + HOJAS:DUREZA + RAIZ:AREA + RAIZ:DUREZA + AREA:DUREZA +
## TALLO:HOJAS:RAIZ + TALLO:HOJAS:AREA + TALLO:HOJAS:DUREZA +
## TALLO:RAIZ:AREA + TALLO:RAIZ:DUREZA + TALLO:AREA:DUREZA +
## HOJAS:RAIZ:AREA + HOJAS:RAIZ:DUREZA + HOJAS:AREA:DUREZA +
## RAIZ:AREA:DUREZA + TALLO:HOJAS:RAIZ:AREA + TALLO:HOJAS:RAIZ:DUREZA +
## TALLO:HOJAS:AREA:DUREZA + TALLO:RAIZ:AREA:DUREZA + HOJAS:RAIZ:AREA:DUREZA +
## TALLO:HOJAS:RAIZ:AREA:DUREZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5503 -0.0420 -0.0070 0.0386 0.7083
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.73e-01 2.39e-01 -0.72 0.471
## TALLO 2.45e-01 2.56e-01 0.96 0.341
## HOJAS 1.47e-01 2.94e-01 0.50 0.619
## RAIZ 9.37e-02 7.69e-01 0.12 0.903
## AREA 6.99e-03 9.09e-03 0.77 0.442
## DUREZA 4.61e-03 4.60e-03 1.00 0.318
## TALLO:HOJAS -1.78e-02 2.53e-02 -0.70 0.482
## TALLO:RAIZ -6.51e-02 4.60e-02 -1.42 0.159
## TALLO:AREA 1.70e-03 2.85e-03 0.59 0.553
## TALLO:DUREZA -4.81e-03 5.12e-03 -0.94 0.349
## HOJAS:RAIZ -7.61e-02 1.06e-01 -0.72 0.473
## HOJAS:AREA -7.62e-03 4.74e-03 -1.61 0.110
## HOJAS:DUREZA -4.67e-03 6.37e-03 -0.73 0.464
## RAIZ:AREA -1.60e-02 1.26e-02 -1.27 0.207
## RAIZ:DUREZA 7.95e-04 1.48e-02 0.05 0.957
## AREA:DUREZA -1.49e-04 1.66e-04 -0.90 0.371
## TALLO:HOJAS:RAIZ 8.49e-03 3.39e-03 2.50 0.013 *
## TALLO:HOJAS:AREA 1.64e-05 2.87e-04 0.06 0.955
## TALLO:HOJAS:DUREZA 5.21e-04 4.64e-04 1.12 0.263
## TALLO:RAIZ:AREA 4.60e-04 6.17e-04 0.75 0.457
## TALLO:RAIZ:DUREZA 1.13e-03 1.02e-03 1.10 0.274
## TALLO:AREA:DUREZA -3.16e-05 6.12e-05 -0.52 0.606
## HOJAS:RAIZ:AREA 3.90e-03 1.74e-03 2.25 0.026 *
## HOJAS:RAIZ:DUREZA 1.55e-03 2.25e-03 0.69 0.492
## HOJAS:AREA:DUREZA 1.85e-04 1.02e-04 1.82 0.071 .
## RAIZ:AREA:DUREZA 3.16e-04 2.68e-04 1.18 0.240
## TALLO:HOJAS:RAIZ:AREA -1.33e-04 5.47e-05 -2.43 0.016 *
## TALLO:HOJAS:RAIZ:DUREZA -1.74e-04 6.70e-05 -2.59 0.010 *
## TALLO:HOJAS:AREA:DUREZA -2.74e-06 5.17e-06 -0.53 0.597
## TALLO:RAIZ:AREA:DUREZA -7.42e-06 1.51e-05 -0.49 0.623
## HOJAS:RAIZ:AREA:DUREZA -8.28e-05 3.87e-05 -2.14 0.034 *
## TALLO:HOJAS:RAIZ:AREA:DUREZA 2.87e-06 1.15e-06 2.51 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.115 on 179 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.858, Adjusted R-squared: 0.833
## F-statistic: 34.9 on 31 and 179 DF, p-value: <2e-16
# Y esto genera otro ligero incremento en la varianza explicada. La forma
# simplificada de escribir este modelo saturado es usar asteriscos entre
# los nombres de las variables lo cual le indica a R que haga todas las
# posibles combinaciones, desde factores individuales hasta la interacción
# de quinto nivel.
M <- lm(FLORES ~ TALLO * HOJAS * RAIZ * AREA * DUREZA)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO * HOJAS * RAIZ * AREA * DUREZA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5503 -0.0420 -0.0070 0.0386 0.7083
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.73e-01 2.39e-01 -0.72 0.471
## TALLO 2.45e-01 2.56e-01 0.96 0.341
## HOJAS 1.47e-01 2.94e-01 0.50 0.619
## RAIZ 9.37e-02 7.69e-01 0.12 0.903
## AREA 6.99e-03 9.09e-03 0.77 0.442
## DUREZA 4.61e-03 4.60e-03 1.00 0.318
## TALLO:HOJAS -1.78e-02 2.53e-02 -0.70 0.482
## TALLO:RAIZ -6.51e-02 4.60e-02 -1.42 0.159
## HOJAS:RAIZ -7.61e-02 1.06e-01 -0.72 0.473
## TALLO:AREA 1.70e-03 2.85e-03 0.59 0.553
## HOJAS:AREA -7.62e-03 4.74e-03 -1.61 0.110
## RAIZ:AREA -1.60e-02 1.26e-02 -1.27 0.207
## TALLO:DUREZA -4.81e-03 5.12e-03 -0.94 0.349
## HOJAS:DUREZA -4.67e-03 6.37e-03 -0.73 0.464
## RAIZ:DUREZA 7.95e-04 1.48e-02 0.05 0.957
## AREA:DUREZA -1.49e-04 1.66e-04 -0.90 0.371
## TALLO:HOJAS:RAIZ 8.49e-03 3.39e-03 2.50 0.013 *
## TALLO:HOJAS:AREA 1.64e-05 2.87e-04 0.06 0.955
## TALLO:RAIZ:AREA 4.60e-04 6.17e-04 0.75 0.457
## HOJAS:RAIZ:AREA 3.90e-03 1.74e-03 2.25 0.026 *
## TALLO:HOJAS:DUREZA 5.21e-04 4.64e-04 1.12 0.263
## TALLO:RAIZ:DUREZA 1.13e-03 1.02e-03 1.10 0.274
## HOJAS:RAIZ:DUREZA 1.55e-03 2.25e-03 0.69 0.492
## TALLO:AREA:DUREZA -3.16e-05 6.12e-05 -0.52 0.606
## HOJAS:AREA:DUREZA 1.85e-04 1.02e-04 1.82 0.071 .
## RAIZ:AREA:DUREZA 3.16e-04 2.68e-04 1.18 0.240
## TALLO:HOJAS:RAIZ:AREA -1.33e-04 5.47e-05 -2.43 0.016 *
## TALLO:HOJAS:RAIZ:DUREZA -1.74e-04 6.70e-05 -2.59 0.010 *
## TALLO:HOJAS:AREA:DUREZA -2.74e-06 5.17e-06 -0.53 0.597
## TALLO:RAIZ:AREA:DUREZA -7.42e-06 1.51e-05 -0.49 0.623
## HOJAS:RAIZ:AREA:DUREZA -8.28e-05 3.87e-05 -2.14 0.034 *
## TALLO:HOJAS:RAIZ:AREA:DUREZA 2.87e-06 1.15e-06 2.51 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.115 on 179 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.858, Adjusted R-squared: 0.833
## F-statistic: 34.9 on 31 and 179 DF, p-value: <2e-16
# El gran problema de la regresión múltiple es definir cual será el modelo
# inicial, llegar al modelo mínimo. Son tantas combinaciones, que
# procedimiento de simplificación por pasos tomaría mucho tiempo
# Exploremos si podemos llegar a una simplificación con otro tipo de
# técnicas. Usaremos modelos de arboles por regresión. Este es un método
# divisivo que a partir de ir ordenando las variables explicativas busca
# dividir la base de datos en dos grupos analizando el comportamiento en
# la variable dependiente. La existencia de está variable dependiente es
# lo que diferencia estos modelos de árboles por regresión de los métodos
# de clasificación multivariados.
# Los modelos de árbol por regresión buscan maximizar la diferencias entre
# las medias de los grupos, y minimizar la devianza dentro de ellos para
# la variable de respuesta. Para cada división se realiza una búsqueda
# exhaustiva y en aquel punto que se maximice las diferencia entre las
# medias y se minimice la devianza, que por ahora podemos pensar
# simplemente como un análogo de la varianza, será el criterio y punto de
# corte para la división
# Para la siguiente sección debemos primero tener instalado los paquetes
# << tree >> y << gmodels >>. Una vez instalados estos paquetes
# activaremos el paquete << tree >> con la función << library() >>.
library(tree)
# Y ahora creamos un nuevo objeto que contenga solamente las variables que
# nos interesan. En este caso estas variables son << FLORES >>, << TALLO
# >>, << HOJAS >>, << RAIZ >>, << AREA >> y << DUREZA >> que son las
# mismas variables que estuvimos usando en la regresión múltiple. Para los
# modelos de árbol, la primer columna será siempre la variable de
# respuesta y el resto serán variables explicativas. El objeto que vamos a
# crear lo denominaremos << DATOS2 >>
DATOS2 <- data.frame(FLORES, TALLO, HOJAS, RAIZ, AREA, DUREZA)
# Y aplicamos la función << tree() >> a << DATOS2 >>
M <- tree(DATOS2)
# Esto genera el modelo de árbol, cuyo resultado numérico podemos ver
# llamando al modelo, en este caso << M >>
M
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 211 20.00 0.30
## 2) TALLO < 11.925 120 0.80 0.10
## 4) HOJAS < 6.29 108 0.40 0.08 *
## 5) HOJAS > 6.29 12 0.05 0.30 *
## 3) TALLO > 11.925 91 4.00 0.60
## 6) AREA < 85.985 65 2.00 0.50
## 12) TALLO < 16.74 31 0.30 0.40 *
## 13) TALLO > 16.74 34 1.00 0.60
## 26) AREA < 55.145 9 0.60 0.70 *
## 27) AREA > 55.145 25 0.60 0.50 *
## 7) AREA > 85.985 26 0.90 0.80
## 14) RAIZ < 4.705 7 0.10 0.60 *
## 15) RAIZ > 4.705 19 0.50 0.80 *
# En donde en cada fila de la tabla se muestra el número de nodo, la
# variable que se utilizó para la división y el criterio, el número de
# observaciones en el nodo, la devianza, y el valor promedio de la
# variable de respuesta. Los asteriscos denota a las ramas terminales del
# árbol,
# Con las funciones << plot() >> y << text() >> logramos una
# representación gráfica del modelo de árbol
plot(M)
text(M)
# Los modelos de regresión por árbol tiene la desventaja de que ajustan un
# único árbol a un conjunto de datos, y este arreglo puede ser singular
# por lo usualmente se critican los modelos por ser sobre estimados. Un
# procedimiento para disminuir este sobrestimación es podar el árbol. Esto
# se logra con la función << prune.tree() >> que tiene como argumentos el
# nombre del modelo y el número de nodos que se requiere representar en el
# modelo simplificado, el argumento << best= >>
M2 <- prune.tree(M, best = 4)
plot(M2)
text(M2)
# La gran desventaja es que nosotros debemos especificar a priori cuantas
# ramas queremos que queden en el árbol simplificado lo que resta
# objetividad al proceso de simplificación.
# Otra aproximación consiste en no generar un único árbol sino generar
# decenas o centenas de árboles. Esto se logra tomando para cada árbol
# solo una muestra aleatoria de los valores existentes, por ejemplo
# tomando el 80% de los datos totales se genera un árbol, tomando otra
# muestra de 80% de los datos al azar se genera otro árbol y así
# sucesivamente hasta completar unos 100 árboles o quizá 500 árboles. Con
# la información de cada árbol y teniendo la colección de árboles se
# reconstruye cuál sería el árbol más parsimonioso y esto genera
# automáticamente un proceso de poda en los árboles ya que aquellos nodos
# de divisiones muy inferiores son poco repetibles entre los diversos
# árboles se perderán en la reconstrucción. Esta técnica se conoce con el
# nombre bosque de árboles no sorprendentemente por el hecho de que
# generamos múltiples árboles lo cual da la idea de tener un bosque.
# Para usar esta técnica debemos primero cargarnos el archivo de funciones
# del curso que se llama << funciones.R >> colocando la ruta y el nombre
# del archivo dentro de la función << source() >> . La ruta debe
# actualizarse según el sistema que ustedes estén utilizando, y recuerden
# que siempre va entre comillas.
source("~/desktop/curso R 2012/funciones.r")
# Ya Con el archivo de funciones cargado podemos usar la función <<
# bosqeu() >> que tiene tres argumentos el nombre de la base de datos que
# en su primera columna contiene la variable de respuesta, el segundo
# argumento es el número de árboles a generar, y el tercer argumento es la
# proporción de datos a utilizar para generar cada árbol
bosque(DATOS2, 500, 0.8)
## [[1]]
## TALLO HOJAS RAIZ AREA DUREZA <leaf>
## 1 345 0 155 0 0 0
## 2 294 131 74 0 0 1
## 3 110 17 65 306 2 0
## 4 78 0 0 0 0 421
## 5 0 5 3 0 0 491
## 6 296 45 11 33 10 105
## 7 35 73 180 138 19 55
## 12 20 7 26 2 3 337
## 13 9 26 12 120 45 183
## 14 7 32 34 55 24 293
## 15 5 21 13 24 3 379
## 26 14 17 3 2 0 176
## 27 15 1 8 0 14 174
##
## [[2]]
## NODO VAR.DIV PUNTO.CORTE DEVIANZA N Y NA
## 1 1 TALLO 12.699 13.3820 168 0.29891 0.000794
## 2 2 TALLO 2.561 0.6975 97 0.09661 0.001012
## 3 3 AREA 105.017 3.4318 71 0.56945 0.002389
## 4 4 <leaf> NA 0.2166 82 0.07092 0.001193
## 5 5 <leaf> NA 0.1333 19 0.26923 0.012778
## 6 6 TALLO 17.122 1.4965 50 0.49018 0.006094
## 7 7 RAIZ 5.148 0.7265 21 0.76170 0.015529
## 12 12 <leaf> NA 0.2226 26 0.41172 0.010668
## 13 13 <leaf> NA 0.8417 26 0.58556 0.012331
## 14 14 <leaf> NA 0.2327 10 0.64167 0.022698
## 15 15 <leaf> NA 0.1524 10 0.87636 0.025793
##
## [[3]]
## _
## Mean_residual_deviance = 0.3560
## Mean explained deviance 0.9734
## [[1]]
## TALLO HOJAS RAIZ AREA DUREZA <leaf>
## 1 345 0 155 0 0 0
## 2 294 131 74 0 0 1
## 3 110 17 65 306 2 0
## 4 78 0 0 0 0 421
## 5 0 5 3 0 0 491
## 6 296 45 11 33 10 105
## 7 35 73 180 138 19 55
## 12 20 7 26 2 3 337
## 13 9 26 12 120 45 183
## 14 7 32 34 55 24 293
## 15 5 21 13 24 3 379
## 26 14 17 3 2 0 176
## 27 15 1 8 0 14 174
##
## [[2]]
## NODO VAR.DIV PUNTO.CORTE DEVIANZA N Y NA
## 1 1 TALLO 12.699 13.3820 168 0.29891 0.000794
## 2 2 TALLO 2.561 0.6975 97 0.09661 0.001012
## 3 3 AREA 105.017 3.4318 71 0.56945 0.002389
## 4 4 <leaf> NA 0.2166 82 0.07092 0.001193
## 5 5 <leaf> NA 0.1333 19 0.26923 0.012778
## 6 6 TALLO 17.122 1.4965 50 0.49018 0.006094
## 7 7 RAIZ 5.148 0.7265 21 0.76170 0.015529
## 12 12 <leaf> NA 0.2226 26 0.41172 0.010668
## 13 13 <leaf> NA 0.8417 26 0.58556 0.012331
## 14 14 <leaf> NA 0.2327 10 0.64167 0.022698
## 15 15 <leaf> NA 0.1524 10 0.87636 0.025793
##
## [[3]]
## _
## Mean_residual_deviance = 0.3560
## Mean explained deviance 0.9734
##
## [[4]]
## [1] 0.07092 0.26923 0.41172 0.58556 0.64167 0.87636
##
## [[5]]
## [1] 0.001193 0.012778 0.010668 0.012331 0.022698 0.025793
# Esto nos genera una salida similar a cuando generamos un solo árbol pero
# la parte gráfica contiene dos elementos, el árbol mismo y una
# representación gráfica de la variable de respuesta en barras con errores
# estándares para observar cómo difiere esta variable entre los grupos que
# se han formado.
# El modelo sugiere que son suficiente los valores del << TALLO >>, <<
# AREA >> y << RAIZ >> para describir la variable de respuesta que en este
# caso es el peso de las flores, << FLORES >>
# Con esta información podemos regresar al modelo de regresión múltiple
# que nos ocupaba e iniciar con un modelo mucho más sencillo en que sólo
# incluye tres variables << TALLO >>, << AREA >> y << RAIZ >>
M <- lm(FLORES ~ TALLO * AREA * RAIZ)
summary(M)
##
## Call:
## lm(formula = FLORES ~ TALLO * AREA * RAIZ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6004 -0.0447 -0.0077 0.0410 0.8318
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.32e-02 3.93e-02 0.84 0.399
## TALLO 2.32e-02 9.58e-03 2.42 0.017 *
## AREA 6.29e-04 1.37e-03 0.46 0.648
## RAIZ -1.94e-02 4.21e-02 -0.46 0.646
## TALLO:AREA -1.15e-04 1.24e-04 -0.93 0.353
## TALLO:RAIZ 6.93e-04 1.65e-03 0.42 0.676
## AREA:RAIZ 7.58e-04 5.44e-04 1.39 0.165
## TALLO:AREA:RAIZ -6.18e-06 2.08e-05 -0.30 0.766
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.119 on 203 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.827, Adjusted R-squared: 0.821
## F-statistic: 139 on 7 and 203 DF, p-value: <2e-16
# Tenemos dos rutas para terminar de simplificar este modo una es una ruta
# automatizada en R que la lleva acabo la función << step() >> en el que a
# partir del modelo inicial se realizará simplificación tomando el
# criterio de información de Akaike como referencia. El valor del criterio
# de información de Akaike podemos extraerlo directamente de cualquier
# modelo que permita una función de verosimilitud con la función en R
# llamada << AIC() >> o alternativamente con la función << extractAIC()
# >>. Note que estas funciones generarán valores diferentes de AIC para e
# mismo modelo. Por ejemplo para el modelo que tenemos ahora << AIC() >>
# genera un valor de -288.0979 mientras que << extractAIC() >> general un
# valor de -888.89. La función << step() >> utiliza esta segunda función,
# << extractAIC() >> para simplificar el modelo. La magnitud de esos
# valores de AIC es irrelevante, lo que es importantes que al comparar
# modelos usando una única función el valor más pequeño de AIC siempre
# denota un mejor modelo. Las diferencias entre las funciones que << AIC()
# >> y << extractAIC() >> radican en que usan diferentes constantes en el
# cálculo.
# Al Aplicar la función << step() >> podemos ver que concluimos con un
# modelo con tres interacciones principales y una interacción doble <<
# AREA:RAIZ >>. Al asignarla la función << step() >> a un objeto este
# objeto que se crea contiene el modelo simplificado y de este podemos
# directamente extraer el resumen del modelo como habíamos hecho con los
# otros modelos que directamente especificamos.
M2 <- step(M)
## Start: AIC=-888.9
## FLORES ~ TALLO * AREA * RAIZ
##
## Df Sum of Sq RSS AIC
## - TALLO:AREA:RAIZ 1 0.00126 2.9 -891
## <none> 2.9 -889
##
## Step: AIC=-890.8
## FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ + AREA:RAIZ
##
## Df Sum of Sq RSS AIC
## - TALLO:RAIZ 1 0.0013 2.90 -893
## - TALLO:AREA 1 0.0158 2.91 -892
## <none> 2.90 -891
## - AREA:RAIZ 1 0.0425 2.94 -890
##
## Step: AIC=-892.7
## FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + AREA:RAIZ
##
## Df Sum of Sq RSS AIC
## - TALLO:AREA 1 0.0145 2.91 -894
## <none> 2.90 -893
## - AREA:RAIZ 1 0.0451 2.94 -891
##
## Step: AIC=-893.6
## FLORES ~ TALLO + AREA + RAIZ + AREA:RAIZ
##
## Df Sum of Sq RSS AIC
## <none> 2.91 -894
## - AREA:RAIZ 1 0.0668 2.98 -891
## - TALLO 1 0.2362 3.15 -879
summary(M2)
##
## Call:
## lm(formula = FLORES ~ TALLO + AREA + RAIZ + AREA:RAIZ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.6027 -0.0434 -0.0078 0.0401 0.8366
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.040221 0.027192 1.48 0.141
## TALLO 0.015281 0.003739 4.09 6.3e-05 ***
## AREA 0.000203 0.000806 0.25 0.801
## RAIZ 0.020195 0.017558 1.15 0.251
## AREA:RAIZ 0.000312 0.000144 2.17 0.031 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.119 on 206 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.826, Adjusted R-squared: 0.823
## F-statistic: 245 on 4 and 206 DF, p-value: <2e-16
# La otra ruta de simplificación es la comparación entre modelos. Para
# realizar esto a partir del modelo inicial simplificamos quitando la
# interacción más compleja que contenga y una vez ajustado el segundo
# modelo se realiza una prueba de análisis de varianza para averiguar si
# la suma de cuadrados residual incremento significativamente. Por ejemplo
# podemos definir el modelo simplificado como el modelo ya existente y
# restándole la sintaxis << TALO:AREA:RAIZ >> que le indica a a la función
# <<lm() >> que elimine la interacción triple del modelo.
M2 <- lm(FLORES ~ TALLO * AREA * RAIZ - TALLO:AREA:RAIZ)
summary(M2)
##
## Call:
## lm(formula = FLORES ~ TALLO * AREA * RAIZ - TALLO:AREA:RAIZ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5996 -0.0433 -0.0082 0.0412 0.8324
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.026733 0.032670 0.82 0.414
## TALLO 0.023721 0.009378 2.53 0.012 *
## AREA 0.000816 0.001218 0.67 0.504
## RAIZ -0.012971 0.036128 -0.36 0.720
## TALLO:AREA -0.000125 0.000119 -1.06 0.293
## TALLO:RAIZ 0.000304 0.001013 0.30 0.764
## AREA:RAIZ 0.000639 0.000369 1.73 0.085 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.119 on 204 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.827, Adjusted R-squared: 0.822
## F-statistic: 163 on 6 and 204 DF, p-value: <2e-16
# La comparación entre modelos la llevamos acabo con la función << anova()
# >> donde los argumentos son los nombres de los dos modelos
anova(M2, M)
## Analysis of Variance Table
##
## Model 1: FLORES ~ TALLO * AREA * RAIZ - TALLO:AREA:RAIZ
## Model 2: FLORES ~ TALLO * AREA * RAIZ
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 204 2.9
## 2 203 2.9 1 0.00126 0.09 0.77
# Podemos observar que de esta simplificación de anova entre modelos nos
# indica que no hay un aumento significativo en la varianza residual ya
# que la probabilidad es mayor que el 0.05. En consecuencia podemos dejar
# esta interacción triple fuera del modelo y seguir simplificando ahora
# removiendo las interacciones dobles.
# Antes de proseguir simplificando veamos una sintaxis alternativa para la
# simplificación. La sintaxis preferida o quizá más fácil de utilizar para
# remover las interacciones es usando la función << update() >> que toma
# como argumentos el nombre del modelo, y <<< ~. >> como segundo
# argumento a partir del cual se restan las interacciones que se van a
# remover. A continuación repetimos la simplificación que ya realizamos
# pero ahora usando esta nueva función para remover la interacción triple.
M2 <- update(M, ~. - TALLO:AREA:RAIZ)
anova(M2, M)
## Analysis of Variance Table
##
## Model 1: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ + AREA:RAIZ
## Model 2: FLORES ~ TALLO * AREA * RAIZ
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 204 2.9
## 2 203 2.9 1 0.00126 0.09 0.77
summary(M2)
##
## Call:
## lm(formula = FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ +
## AREA:RAIZ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5996 -0.0433 -0.0082 0.0412 0.8324
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.026733 0.032670 0.82 0.414
## TALLO 0.023721 0.009378 2.53 0.012 *
## AREA 0.000816 0.001218 0.67 0.504
## RAIZ -0.012971 0.036128 -0.36 0.720
## TALLO:AREA -0.000125 0.000119 -1.06 0.293
## TALLO:RAIZ 0.000304 0.001013 0.30 0.764
## AREA:RAIZ 0.000639 0.000369 1.73 0.085 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.119 on 204 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.827, Adjusted R-squared: 0.822
## F-statistic: 163 on 6 and 204 DF, p-value: <2e-16
# Ahora si procedamos a remover las interacciones dobles. Para esto, del
# modelo sin interacción triple vamos a remover cada una de las
# interacciones dobles, una por una siempre del del mismo modelo. Al
# remover cada una de las interacciones dobles comparamos con la prueba de
# << anova() >> para decidir la interacción debe quedarse o puede ser
# removida del modelo.
# Por ejemplo removemos la interacción << TALLO:AREA >> y observamos que
# nos da un valor de probabilidad de 0.2927 con lo cual podemos concluir
# que esta interacción doble debe ser removida del modelo
M3 <- update(M2, ~. - TALLO:AREA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: FLORES ~ TALLO + AREA + RAIZ + TALLO:RAIZ + AREA:RAIZ
## Model 2: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ + AREA:RAIZ
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 205 2.91
## 2 204 2.90 1 0.0158 1.11 0.29
# Para la interacción entre << TALLO:RAIZ >> también notamos que el valor
# de probabilidad al comparar los modelos es de 0.7642 lo cual también nos
# permite remover esta interacción definitivamente del modelo
M3 <- update(M2, ~. - TALLO:RAIZ)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + AREA:RAIZ
## Model 2: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ + AREA:RAIZ
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 205 2.9
## 2 204 2.9 1 0.00128 0.09 0.76
# Para la otra interacción doble << RAIZ:AREA >> notamos que la
# probabilidad al comparar los modelos es de 0.085 lo cual es
# marginalmente no significativo y estrictamente debemos removerla del
# modelo
M3 <- update(M2, ~. - RAIZ:AREA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ
## Model 2: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ + AREA:RAIZ
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 205 2.94
## 2 204 2.90 1 0.0425 3 0.085 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ahora, ya habiendo decidido que las tres interacciones dobles se pueden
# eliminar al realizar pruebas separadas para cada una, ya podemos
# quitarlas en conjunto del modelo, y volvemos a comparar contra el modelo
# que si contiene las interacciones dobles. Al hacer esto observamos que
# eliminar las tres interacciones dobles no generó un aumento
# significativo en la varianza residual.
M3 <- update(M2, ~. - RAIZ:AREA - TALLO:RAIZ - TALLO:AREA)
anova(M3, M2)
## Analysis of Variance Table
##
## Model 1: FLORES ~ TALLO + AREA + RAIZ
## Model 2: FLORES ~ TALLO + AREA + RAIZ + TALLO:AREA + TALLO:RAIZ + AREA:RAIZ
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 207 2.98
## 2 204 2.90 3 0.0826 1.94 0.12
# Esto nos deja con un modelo bastante sencillo en el que sólo se
# encuentran los tres factores principales y el intercepto
summary(M3)
##
## Call:
## lm(formula = FLORES ~ TALLO + AREA + RAIZ)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5956 -0.0476 -0.0068 0.0435 0.8327
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.010429 0.014135 -0.74 0.46147
## TALLO 0.011577 0.003358 3.45 0.00068 ***
## AREA 0.001807 0.000325 5.55 8.6e-08 ***
## RAIZ 0.044855 0.013519 3.32 0.00107 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.12 on 207 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.822, Adjusted R-squared: 0.82
## F-statistic: 319 on 3 and 207 DF, p-value: <2e-16
# Ahora podemos proceder a graficar una de estas relaciones, por jemplo
# peso de las flores << FLORES >> y el peso de los tallos << TALLO >>.
par(mfrow = c(1, 1), mar = c(4, 4, 2, 1))
plot(FLORES ~ TALLO, xlab = "Peso del tall (g)", ylab = "Peso de las flores (g)",
col = "tomato3", pch = 19)
# La función << abline() >> en el caso de la regresión múltiple pierde su
# utilidad de pintar la linea de regresión mediante los valores del
# intercepto y la pendiente. Los estimados de pendiente que están en el
# predictor lineal para cada una de las variables corresponderían a la
# pendiente cuando el estimado de las otras dos variables es igual a cero,
# y al pintar la linea de regresión vemos que esta queda completamente
# fuera del área de los puntos en el gráfico.
abline(-0.0104286, 0.0115773)
# Resulta más útil dibujar la linea de regresión de una de estas variables
# cuando los estimados de las otras dos variables son iguales al promedio
# de las variables mismas. Esto únicamente lo podemos ahora hacer con la
# función << predict() >> para la cual como recordarán hay que generar una
# nueva serie de datos para los cuales se harán las predicciones. La
# secuencia de valores debemos incluirla en un << data.frame >> de la
# misma forma que lo hicimos antes sólo que ahora este << data.frame >>
# incluye otras dos variables cuyos valores estarán fijos en sus
# promedios.
SEQ <- seq(0, 30, length = 30)
ND <- data.frame(TALLO = SEQ, RAIZ = mean(RAIZ), AREA = mean(AREA, na.rm = T))
PRED <- predict(M3, ND)
lines(SEQ, PRED, col = "red")
# Recuerden que si deseamos graficar no sonó la línea de regresión sino
# además el intervalo de confianza al usar la función << predict() >>
# debemos solicitar que las observaciones vengan acompañadas del error
# estándar definiendo el argumento << se=TRUE >>. Y con esta único uso de
# la función podemos graficar la línea de regresión y los límites del
# intervalo de confianza.
PRED <- predict(M3, ND, se = TRUE)
T <- qt(0.025, 207)
EE <- PRED$se.fit
LR <- PRED$fit
plot(FLORES ~ TALLO, xlab = "Peso del tall (g)", ylab = "Peso de las flores (g)",
col = "tomato3", pch = 19)
lines(SEQ, LR, lty = 1, col = "blue")
lines(SEQ, LR + T * EE, lty = 2)
lines(SEQ, PRED$fit - T * PRED$se.fit, lty = 2)
