# ROGER GUEVARA (roger.guevara@inecol.edu.mx)
# INSTITUTO DE ECOLOGIA A.C RED DE BIOLOGÍA EVOLUTIVA
# ESTADISTICA COMPUTACIONAL EN R
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, SUELO == "P" & FLORES > 0)
attach(DATOS2)
## The following object(s) are masked from 'DATOS':
##
## ABORTOS, AFE, ALT_TOT, ARBUSCULOS, AREA, CAPSULAS,
## DIAS_A_LA_FLOR, DUREZA, ESPORAS, FLORES, FRUT_INMADUROS,
## FRUT_MADUROS, HERBIVORIA, HIFAS, HOJAS, HOJAS_TOT, INMADUROS,
## INVERNADERO, LONG_HIFAS, LUZ, MICORRIZA, NOTAS, NUM,
## NUM_FLORES, PESO_HOJAS, PESO_TOTAL, RAIZ, SEMILLAS, SUELO,
## TALLO, TOT_FRUTOS, TOTAL, VESICULAS
detach(DATOS)
rm(DATOS)
M <- glm(FLORES ~ 1)
summary(M)
##
## Call:
## glm(formula = FLORES ~ 1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.0647 -0.0447 -0.0147 0.0253 0.2553
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.07472 0.00561 13.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.003338)
##
## Null deviance: 0.35044 on 105 degrees of freedom
## Residual deviance: 0.35044 on 105 degrees of freedom
## AIC: -300.7
##
## Number of Fisher Scoring iterations: 2
var(FLORES) * (length(FLORES) - 1)
## [1] 0.3504
# Poisson
M <- glm(NUM_FLORES ~ 1, poisson)
summary(M)
##
## Call:
## glm(formula = NUM_FLORES ~ 1, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.062 -0.307 -0.307 0.328 1.897
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9029 0.0621 14.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 59.252 on 104 degrees of freedom
## Residual deviance: 59.252 on 104 degrees of freedom
## (1 observation deleted due to missingness)
## AIC: 345.1
##
## Number of Fisher Scoring iterations: 4
MU <- mean(NUM_FLORES, na.rm = T)
2 * sum(NUM_FLORES * log(NUM_FLORES/MU) - (NUM_FLORES - MU), na.rm = T)
## [1] 59.25
2 * sum(NUM_FLORES * log(NUM_FLORES/MU), na.rm = T)
## [1] 59.25
# Binomial
Y <- cbind(FRUT_MADUROS, FRUT_INMADUROS)
M <- glm(Y ~ 1, binomial)
summary(M)
##
## Call:
## glm(formula = Y ~ 1, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.45 -1.22 0.00 1.13 2.26
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.108 0.134 0.8 0.42
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 146.12 on 93 degrees of freedom
## AIC: 213.7
##
## Number of Fisher Scoring iterations: 3
M <- glm(Y ~ LUZ * MICORRIZA * TALLO, binomial("logit"))
summary(M)
##
## Call:
## glm(formula = Y ~ LUZ * MICORRIZA * TALLO, family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.286 -0.939 0.000 0.853 2.647
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.255 0.820 1.53 0.1256
## LUZPL -2.199 1.067 -2.06 0.0394 *
## MICORRIZASIN -2.272 1.025 -2.22 0.0266 *
## TALLO -0.921 0.484 -1.90 0.0570 .
## LUZPL:MICORRIZASIN 3.184 1.502 2.12 0.0340 *
## LUZPL:TALLO 2.205 1.173 1.88 0.0602 .
## MICORRIZASIN:TALLO 2.016 0.687 2.93 0.0034 **
## LUZPL:MICORRIZASIN:TALLO -1.980 1.896 -1.04 0.2964
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 126.83 on 86 degrees of freedom
## AIC: 208.4
##
## Number of Fisher Scoring iterations: 4
anova(M, test = "Chi", dispersion = 1.4747)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 93 146
## LUZ 1 0.98 92 145 0.416
## MICORRIZA 1 6.07 91 139 0.043 *
## TALLO 1 1.01 90 138 0.407
## LUZ:MICORRIZA 1 0.30 89 138 0.652
## LUZ:TALLO 1 1.80 88 136 0.269
## MICORRIZA:TALLO 1 8.05 87 128 0.019 *
## LUZ:MICORRIZA:TALLO 1 1.08 86 127 0.392
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(M, dispersion = 1.4747)
##
## Call:
## glm(formula = Y ~ LUZ * MICORRIZA * TALLO, family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.286 -0.939 0.000 0.853 2.647
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.255 0.995 1.26 0.207
## LUZPL -2.199 1.296 -1.70 0.090 .
## MICORRIZASIN -2.272 1.244 -1.83 0.068 .
## TALLO -0.921 0.587 -1.57 0.117
## LUZPL:MICORRIZASIN 3.184 1.824 1.75 0.081 .
## LUZPL:TALLO 2.205 1.425 1.55 0.122
## MICORRIZASIN:TALLO 2.016 0.835 2.42 0.016 *
## LUZPL:MICORRIZASIN:TALLO -1.980 2.303 -0.86 0.390
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.475)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 126.83 on 86 degrees of freedom
## AIC: 208.4
##
## Number of Fisher Scoring iterations: 4
M2 <- update(M, ~. - LUZ:MICORRIZA:TALLO)
anova(M2, M, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO + MICORRIZA:TALLO
## Model 2: Y ~ LUZ * MICORRIZA * TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 87 128
## 2 86 127 1 1.08 0.3
summary(M2)
##
## Call:
## glm(formula = Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO +
## MICORRIZA:TALLO, family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.084 -0.986 0.000 0.847 2.635
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.052 0.789 1.33 0.1825
## LUZPL -1.673 0.928 -1.80 0.0713 .
## MICORRIZASIN -1.922 0.960 -2.00 0.0452 *
## TALLO -0.795 0.463 -1.72 0.0860 .
## LUZPL:MICORRIZASIN 1.882 0.819 2.30 0.0216 *
## LUZPL:TALLO 1.484 0.921 1.61 0.1073
## MICORRIZASIN:TALLO 1.764 0.635 2.78 0.0055 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 127.91 on 87 degrees of freedom
## AIC: 207.5
##
## Number of Fisher Scoring iterations: 4
anova(M2, test = "Chi", dispersion = 127.91/87)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 93 146
## LUZ 1 0.98 92 145 0.415
## MICORRIZA 1 6.07 91 139 0.042 *
## TALLO 1 1.01 90 138 0.407
## LUZ:MICORRIZA 1 0.30 89 138 0.652
## LUZ:TALLO 1 1.80 88 136 0.268
## MICORRIZA:TALLO 1 8.05 87 128 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M3 <- update(M2, ~. - LUZ:MICORRIZA)
anova(M3, M2, test = "Chi", dispersion = 127.91/87)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA + TALLO + LUZ:TALLO + MICORRIZA:TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO + MICORRIZA:TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 88 133
## 2 87 128 1 5.38 0.056 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M3 <- update(M2, ~. - LUZ:TALLO)
anova(M3, M2, test = "Chi", dispersion = 127.91/87)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + MICORRIZA:TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO + MICORRIZA:TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 88 131
## 2 87 128 1 2.66 0.18
M3 <- update(M2, ~. - TALLO:MICORRIZA)
anova(M3, M2, test = "Chi", dispersion = 127.91/87)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO + MICORRIZA:TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 88 136
## 2 87 128 1 8.05 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M3 <- update(M2, ~. - LUZ:TALLO - LUZ:MICORRIZA)
anova(M3, M2, test = "Chi", dispersion = 127.91/87)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA + TALLO + MICORRIZA:TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO + LUZ:MICORRIZA + LUZ:TALLO + MICORRIZA:TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 89 135
## 2 87 128 2 7.33 0.083 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(M3, test = "Chi", dispersion = 135.24/89)
##
## Call:
## glm(formula = Y ~ LUZ + MICORRIZA + TALLO + MICORRIZA:TALLO,
## family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.215 -1.011 0.000 0.939 2.703
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.41118 0.68374 -0.60 0.55
## LUZPL 0.49771 0.46859 1.06 0.29
## MICORRIZASIN -0.09607 0.71802 -0.13 0.89
## TALLO 0.00629 0.42048 0.01 0.99
## MICORRIZASIN:TALLO 0.78622 0.58288 1.35 0.18
##
## (Dispersion parameter for binomial family taken to be 1.52)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 135.24 on 89 degrees of freedom
## AIC: 210.8
##
## Number of Fisher Scoring iterations: 4
anova(M3, test = "Chi", dispersion = 135.24/89)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 93 146
## LUZ 1 0.98 92 145 0.423
## MICORRIZA 1 6.07 91 139 0.046 *
## TALLO 1 1.01 90 138 0.414
## MICORRIZA:TALLO 1 2.83 89 135 0.172
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M4 <- update(M3, ~. - TALLO:MICORRIZA)
anova(M4, M3, test = "Chi", dispersion = 135.24/89)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA + TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO + MICORRIZA:TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 90 138
## 2 89 135 1 2.83 0.17
summary(M4, test = "Chi", dispersion = 138.07/90)
##
## Call:
## glm(formula = Y ~ LUZ + MICORRIZA + TALLO, family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.439 -0.981 0.000 0.986 2.593
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.784 0.630 -1.24 0.213
## LUZPL 0.490 0.473 1.04 0.300
## MICORRIZASIN 0.756 0.357 2.12 0.034 *
## TALLO 0.294 0.362 0.81 0.417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.534)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 138.07 on 90 degrees of freedom
## AIC: 211.7
##
## Number of Fisher Scoring iterations: 3
anova(M4, test = "Chi", dispersion = 138.07/90)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 93 146
## LUZ 1 0.98 92 145 0.425
## MICORRIZA 1 6.07 91 139 0.047 *
## TALLO 1 1.01 90 138 0.417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M5 <- update(M4, ~. - MICORRIZA)
anova(M5, M4, test = "Chi", dispersion = 138.07/90)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 91 145
## 2 90 138 1 7.02 0.032 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M5 <- update(M4, ~. - TALLO)
anova(M5, M4, test = "Chi", dispersion = 138.07/90)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + MICORRIZA
## Model 2: Y ~ LUZ + MICORRIZA + TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 91 139
## 2 90 138 1 1.01 0.42
M5 <- update(M4, ~. - LUZ)
anova(M5, M4, test = "Chi", dispersion = 138.07/90)
## Analysis of Deviance Table
##
## Model 1: Y ~ MICORRIZA + TALLO
## Model 2: Y ~ LUZ + MICORRIZA + TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 91 140
## 2 90 138 1 1.66 0.3
M5 <- update(M4, ~. - TALLO - LUZ)
anova(M5, M4, test = "Chi", dispersion = 138.07/90)
## Analysis of Deviance Table
##
## Model 1: Y ~ MICORRIZA
## Model 2: Y ~ LUZ + MICORRIZA + TALLO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 92 140
## 2 90 138 2 1.71 0.57
summary(M5, test = "Chi", dispersion = 139.78/92)
##
## Call:
## glm(formula = Y ~ MICORRIZA, family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.73 -1.07 0.00 1.00 2.57
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.249 0.242 -1.03 0.305
## MICORRIZASIN 0.683 0.336 2.03 0.042 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.519)
##
## Null deviance: 146.12 on 93 degrees of freedom
## Residual deviance: 139.78 on 92 degrees of freedom
## AIC: 209.4
##
## Number of Fisher Scoring iterations: 3
anova(M5, test = "Chi", dispersion = 139.78/92)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 93 146
## MICORRIZA 1 6.35 92 140 0.041 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ND <- data.frame(MICORRIZA = c("CON", "SIN"))
library(boot)
inv.logit(predict(M5, ND))
## 1 2
## 0.4381 0.6068
predict(M5, ND, se = T)
## $fit
## 1 2
## -0.2489 0.4340
##
## $se.fit
## 1 2
## 0.1967 0.1893
##
## $residual.scale
## [1] 1
inv.logit(-0.248896 - 1.985802 * 0.1966934)
## [1] 0.3454
inv.logit(-0.248896 + 1.985802 * 0.1966934)
## [1] 0.5354
inv.logit(0.4340385 - 1.985802 * 0.1892701)
## [1] 0.5145
inv.logit(0.4340385 + 1.985802 * 0.1892701)
## [1] 0.6921
inv.logit(0.4340385) - 0.5145423
## [1] 0.0923
inv.logit(0.4340385) - 0.6920864
## [1] -0.08525
MEDIAS <- predict(M5, ND, se = T)$fit
ERROR <- predict(M5, ND, se = T)$se.fit
source("~/desktop/curso R 2012/funciones.r")
barplot.error(MEDIAS, ERROR)
barplot.error(MEDIAS, ERROR, func = "logit")
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 5)':
##
## 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)
Y <- cbind(FRUT_MADUROS, FRUT_INMADUROS)
M <- glm(Y ~ LUZ/(MICORRIZA * TALLO), binomial("logit"))
summary(M)
##
## Call:
## glm(formula = Y ~ LUZ/(MICORRIZA * TALLO), family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.040 -0.868 -0.113 0.929 3.332
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.6524 0.6137 -2.69 0.0071 **
## LUZPL 2.0901 0.8418 2.48 0.0130 *
## LUZML:MICORRIZASIN 4.1729 0.9033 4.62 3.8e-06 ***
## LUZPL:MICORRIZASIN -1.1311 0.7455 -1.52 0.1292
## LUZML:TALLO 0.0708 0.0297 2.38 0.0171 *
## LUZPL:TALLO -0.0158 0.0377 -0.42 0.6757
## LUZML:MICORRIZASIN:TALLO -0.2112 0.0457 -4.62 3.8e-06 ***
## LUZPL:MICORRIZASIN:TALLO 0.0606 0.0464 1.31 0.1919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 202.33 on 108 degrees of freedom
## Residual deviance: 168.38 on 101 degrees of freedom
## AIC: 472.2
##
## Number of Fisher Scoring iterations: 4
par(mfrow = c(2, 2))
plot(M)
ks.test(M$residuals, "pnorm")
##
## One-sample Kolmogorov-Smirnov test
##
## data: M$residuals
## D = 0.089, p-value = 0.3538
## alternative hypothesis: two-sided
# Binomial
Y <- cbind(FRUT_MADUROS, FRUT_INMADUROS)
M <- glm(Y ~ LUZ/(MICORRIZA * HOJAS), binomial("logit"))
summary(M)
##
## Call:
## glm(formula = Y ~ LUZ/(MICORRIZA * HOJAS), family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8758 -1.0215 -0.0558 0.8576 2.8865
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.5244 0.5492 -0.95 0.33969
## LUZPL 2.4908 0.8790 2.83 0.00460 **
## LUZML:MICORRIZASIN 3.9114 1.0557 3.70 0.00021 ***
## LUZPL:MICORRIZASIN -3.2254 0.9484 -3.40 0.00067 ***
## LUZML:HOJAS 0.0328 0.0573 0.57 0.56724
## LUZPL:HOJAS -0.2106 0.0803 -2.62 0.00875 **
## LUZML:MICORRIZASIN:HOJAS -0.4117 0.1124 -3.66 0.00025 ***
## LUZPL:MICORRIZASIN:HOJAS 0.3587 0.1080 3.32 0.00089 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 202.33 on 108 degrees of freedom
## Residual deviance: 166.40 on 101 degrees of freedom
## AIC: 470.2
##
## Number of Fisher Scoring iterations: 4
par(mfrow = c(2, 2))
plot(M)
ks.test(M$residuals, "pnorm")
##
## One-sample Kolmogorov-Smirnov test
##
## data: M$residuals
## D = 0.0884, p-value = 0.3618
## alternative hypothesis: two-sided
par(mfrow = c(1, 1))
anova(M, test = "Chi", dispersion = 166.4/101)
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 108 202
## LUZ 1 6.86 107 196 0.04127 *
## LUZ:MICORRIZA 2 1.03 105 194 0.73217
## LUZ:HOJAS 2 2.75 103 192 0.43442
## LUZ:MICORRIZA:HOJAS 2 25.29 101 166 0.00046 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M2 <- update(M, ~. - LUZ:MICORRIZA:HOJAS)
anova(M2, M, test = "Chi", dispersion = 166.4/101)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + LUZ:MICORRIZA + LUZ:HOJAS
## Model 2: Y ~ LUZ/(MICORRIZA * HOJAS)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 103 192
## 2 101 166 2 25.3 0.00046 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
k <- 0
par(mfrow = c(2, 2))
for (i in levels(LUZ)) {
for (j in levels(MICORRIZA)) {
k <- k + 1
plot(c(min(HOJAS), max(HOJAS)), c(0, 1), type = "n", main = paste(i,
j))
HOJA2 <- HOJAS[LUZ == i & MICORRIZA == j]
P <- Y[LUZ == i & MICORRIZA == j, 1]/rowSums(Y)[LUZ == i & MICORRIZA ==
j]
points(HOJA2, P)
text(12, 1, paste("b =", round(PEND[k, 1], 2)))
HOJA2 <- seq(min(HOJA2), max(HOJA2), length = 30)
ND <- data.frame(HOJAS = HOJA2, LUZ = i, MICORRIZA = j)
PRED <- predict(M, ND)
lines(HOJA2, inv.logit(PRED))
}
}
## Error: object 'PEND' not found
MLCON <- c(0, 0, 0, 0, 1, 0, 0, 0)
MLSIN <- c(0, 0, 0, 0, 1, 0, 1, 0)
PLCON <- c(0, 0, 0, 0, 0, 1, 0, 0)
PLSIN <- c(0, 0, 0, 0, 0, 1, 0, 1)
libary(gmodels)
## Error: could not find function "libary"
PEND <- estimable(M, rbind(MLCON, MLSIN, PLCON, PLSIN))
## Error: could not find function "estimable"
DMLCONMLSIN <- MLCON - MLSIN
D_PLCONPLSIN <- PLCON - PLSIN
estimable(M, rbind(DMLCONMLSIN, D_PLCONPLSIN))
## Error: could not find function "estimable"
summary(M)$coefficients
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.52438 0.54921 -0.9548 0.3396907
## LUZPL 2.49076 0.87902 2.8336 0.0046031
## LUZML:MICORRIZASIN 3.91144 1.05572 3.7050 0.0002114
## LUZPL:MICORRIZASIN -3.22540 0.94838 -3.4010 0.0006715
## LUZML:HOJAS 0.03277 0.05727 0.5721 0.5672435
## LUZPL:HOJAS -0.21064 0.08034 -2.6217 0.0087495
## LUZML:MICORRIZASIN:HOJAS -0.41175 0.11239 -3.6634 0.0002488
## LUZPL:MICORRIZASIN:HOJAS 0.35867 0.10797 3.3220 0.0008938
# BINOMIAL (BINARIO)
DATOS <- read.table("~/desktop/curso R 2012/TOLOACHE.csv", header = TRUE, sep = ",")
attach(DATOS)
## The following object(s) are masked from 'DATOS2 (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
## The following object(s) are masked from 'DATOS2 (position 5)':
##
## 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 == "P" & FLORES > 0)
attach(DATOS2)
## The following object(s) are masked from 'DATOS':
##
## ABORTOS, AFE, ALT_TOT, ARBUSCULOS, AREA, CAPSULAS,
## DIAS_A_LA_FLOR, DUREZA, ESPORAS, FLORES, FRUT_INMADUROS,
## FRUT_MADUROS, HERBIVORIA, HIFAS, HOJAS, HOJAS_TOT, INMADUROS,
## INVERNADERO, LONG_HIFAS, LUZ, MICORRIZA, NOTAS, NUM,
## NUM_FLORES, PESO_HOJAS, PESO_TOTAL, RAIZ, SEMILLAS, SUELO,
## TALLO, TOT_FRUTOS, TOTAL, VESICULAS
## 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
## The following object(s) are masked from 'DATOS2 (position 6)':
##
## 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)
names(DATOS2)
## [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"
Y <- (FRUT_MADUROS > 0) * 1
M <- glm(Y ~ LUZ/(MICORRIZA * AREA), binomial("logit"))
summary(M)
##
## Call:
## glm(formula = Y ~ LUZ/(MICORRIZA * AREA), family = binomial("logit"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.587 -1.127 0.596 0.948 1.515
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6539 1.4900 1.78 0.075 .
## LUZPL -3.4032 1.9327 -1.76 0.078 .
## LUZML:MICORRIZASIN -4.6065 2.3094 -1.99 0.046 *
## LUZPL:MICORRIZASIN -2.0013 2.1341 -0.94 0.348
## LUZML:AREA -0.0747 0.0436 -1.71 0.086 .
## LUZPL:AREA 0.0332 0.0416 0.80 0.425
## LUZML:MICORRIZASIN:AREA 0.1853 0.0803 2.31 0.021 *
## LUZPL:MICORRIZASIN:AREA 0.1141 0.0891 1.28 0.200
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 140.50 on 105 degrees of freedom
## Residual deviance: 125.25 on 98 degrees of freedom
## AIC: 141.2
##
## Number of Fisher Scoring iterations: 4
par(mfrow = c(2, 2))
plot(M)
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.4774, p-value < 2.2e-16
## alternative hypothesis: two-sided
par(mfrow = c(1, 1))
anova(M, test = "Chi")
## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Y
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 105 140
## LUZ 1 0.26 104 140 0.610
## LUZ:MICORRIZA 2 2.57 102 138 0.277
## LUZ:AREA 2 3.92 100 134 0.141
## LUZ:MICORRIZA:AREA 2 8.50 98 125 0.014 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M2 <- update(M, ~. - LUZ:MICORRIZA:AREA)
anova(M2, M, test = "Chi", dispersion = 125.25/98)
## Analysis of Deviance Table
##
## Model 1: Y ~ LUZ + LUZ:MICORRIZA + LUZ:AREA
## Model 2: Y ~ LUZ/(MICORRIZA * AREA)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 100 134
## 2 98 125 2 8.5 0.036 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MLCON <- c(0, 0, 0, 0, 1, 0, 0, 0)
MLSIN <- c(0, 0, 0, 0, 1, 0, 1, 0)
PLCON <- c(0, 0, 0, 0, 0, 1, 0, 0)
PLSIN <- c(0, 0, 0, 0, 0, 1, 0, 1)
library(gmodels)
PEND <- estimable(M, rbind(MLCON, MLSIN, PLCON, PLSIN))
DMLCONMLSIN <- MLCON - MLSIN
D_PLCONPLSIN <- PLCON - PLSIN
estimable(M, rbind(DMLCONMLSIN, D_PLCONPLSIN))
## Estimate Std. Error X^2 value DF Pr(>|X^2|)
## DMLCONMLSIN -0.1853 0.08031 5.321 1 0.02108
## D_PLCONPLSIN -0.1141 0.08908 1.640 1 0.20036
k <- 0
par(mfrow = c(2, 2))
for (i in levels(LUZ)) {
for (j in levels(MICORRIZA)) {
k <- k + 1
plot(c(min(AREA), max(AREA)), c(0, 1), type = "n", main = paste(i, j))
AREA2 <- AREA[LUZ == i & MICORRIZA == j]
P <- Y[LUZ == i & MICORRIZA == j]
points(AREA2, P)
text(50, 0.8, paste("b =", round(PEND[k, 1], 2)))
AREA2 <- seq(min(AREA2), max(AREA2), length = 30)
ND <- data.frame(AREA = AREA2, LUZ = i, MICORRIZA = j)
PRED <- predict(M, ND, se = T)
T <- qt(0.025, length(P))
lines(AREA2, inv.logit(PRED$fit))
lines(AREA2, inv.logit(PRED$fit - T * PRED$se.fit), lty = 2)
lines(AREA2, inv.logit(PRED$fit + T * PRED$se.fit), lty = 2)
}
}
