# 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 == "R" & NUM_FLORES > 0 & AREA > 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(NUM_FLORES ~ LUZ/(AREA * HOJAS), poisson("log"))
par(mfrow = c(2, 2), mar = c(2, 2, 1, 1))
plot(M)
ks.test(M$residuals, "pnorm")
##
## One-sample Kolmogorov-Smirnov test
##
## data: M$residuals
## D = 0.345, p-value = 2.799e-11
## alternative hypothesis: two-sided
shapiro.test(M$residuals)
##
## Shapiro-Wilk normality test
##
## data: M$residuals
## W = 0.9843, p-value = 0.2509
summary(M, test = "Chi")
##
## Call:
## glm(formula = NUM_FLORES ~ LUZ/(AREA * HOJAS), family = poisson("log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1640 -0.3680 -0.0724 0.4290 1.4586
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.844649 0.586063 3.15 0.0016 **
## LUZPL -0.454664 0.759482 -0.60 0.5494
## LUZML:AREA 0.002260 0.007929 0.28 0.7757
## LUZPL:AREA 0.004998 0.004850 1.03 0.3028
## LUZML:HOJAS 0.070639 0.058704 1.20 0.2289
## LUZPL:HOJAS 0.098306 0.057939 1.70 0.0898 .
## LUZML:AREA:HOJAS -0.000165 0.000753 -0.22 0.8271
## LUZPL:AREA:HOJAS -0.000379 0.000551 -0.69 0.4910
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 75.900 on 104 degrees of freedom
## Residual deviance: 43.094 on 97 degrees of freedom
## AIC: 509.2
##
## Number of Fisher Scoring iterations: 4
anova(M, test = "Chi", dispersion = 0.444268)
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: NUM_FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 104 75.9
## LUZ 1 11.42 103 64.5 4.0e-07 ***
## LUZ:AREA 2 9.93 101 54.6 1.4e-05 ***
## LUZ:HOJAS 2 10.94 99 43.6 4.5e-06 ***
## LUZ:AREA:HOJAS 2 0.53 97 43.1 0.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pchisq(11.4174/0.444268, 1, lower.tail = FALSE)
## [1] 3.99e-07
M2 <- update(M, ~. - LUZ:HOJAS:AREA)
anova(M2, M, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: NUM_FLORES ~ LUZ + LUZ:AREA + LUZ:HOJAS
## Model 2: NUM_FLORES ~ LUZ/(AREA * HOJAS)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 99 43.6
## 2 97 43.1 2 0.525 0.77
M3 <- update(M2, ~. - LUZ:HOJAS)
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: NUM_FLORES ~ LUZ + LUZ:AREA
## Model 2: NUM_FLORES ~ LUZ + LUZ:AREA + LUZ:HOJAS
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 101 54.6
## 2 99 43.6 2 10.9 0.0042 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M3 <- update(M2, ~. - LUZ:AREA)
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: NUM_FLORES ~ LUZ + LUZ:HOJAS
## Model 2: NUM_FLORES ~ LUZ + LUZ:AREA + LUZ:HOJAS
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 101 47.1
## 2 99 43.6 2 3.52 0.17
anova(M3, test = "Chi")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: NUM_FLORES
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 104 75.9
## LUZ 1 11.4 103 64.5 0.00073 ***
## LUZ:HOJAS 2 17.3 101 47.1 0.00017 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(M3, test = "Chi")
##
## Call:
## glm(formula = NUM_FLORES ~ LUZ + LUZ:HOJAS, family = poisson("log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9568 -0.3689 -0.0487 0.4456 1.6689
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9653 0.2325 8.45 <2e-16 ***
## LUZPL -0.2869 0.3187 -0.90 0.3682
## LUZML:HOJAS 0.0630 0.0241 2.61 0.0090 **
## LUZPL:HOJAS 0.0805 0.0247 3.26 0.0011 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 75.900 on 104 degrees of freedom
## Residual deviance: 47.139 on 101 degrees of freedom
## AIC: 505.2
##
## Number of Fisher Scoring iterations: 4
COLOR <- c("red", "green")
par(mfrow = c(1, 1), mar = c(4, 4, 1, 1))
plot(NUM_FLORES ~ HOJAS, type = "n")
k <- 0
for (i in levels(LUZ)) {
k <- k + 1
points(NUM_FLORES[LUZ == i] ~ HOJAS[LUZ == i], col = COLOR[k], pch = 15)
SEQ <- HOJAS[LUZ == i]
SEQ <- seq(min(SEQ), max(SEQ), length = 30)
ND <- data.frame(HOJAS = SEQ, LUZ = i)
PRED <- predict(M3, ND)
lines(SEQ, exp(PRED), col = COLOR[k], lwd = 3)
}
summary(M3)
##
## Call:
## glm(formula = NUM_FLORES ~ LUZ + LUZ:HOJAS, family = poisson("log"))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9568 -0.3689 -0.0487 0.4456 1.6689
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9653 0.2325 8.45 <2e-16 ***
## LUZPL -0.2869 0.3187 -0.90 0.3682
## LUZML:HOJAS 0.0630 0.0241 2.61 0.0090 **
## LUZPL:HOJAS 0.0805 0.0247 3.26 0.0011 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 75.900 on 104 degrees of freedom
## Residual deviance: 47.139 on 101 degrees of freedom
## AIC: 505.2
##
## Number of Fisher Scoring iterations: 4
P_ML <- c(0, 0, 1, 0)
P_PL <- c(0, 0, 0, 1)
DIF <- P_ML - P_PL
library(gmodels)
estimable(M3, rbind(P_ML, P_PL, DIF))
## Estimate Std. Error X^2 value DF Pr(>|X^2|)
## P_ML 0.06299 0.02412 6.8216 1 0.009006
## P_PL 0.08054 0.02471 10.6216 1 0.001118
## DIF -0.01755 0.03453 0.2583 1 0.611266
# Tablas de contingencia----------------------- Para poner a prueba si la
# frecuencia con que ocurren dos eventos son independientes se utilizan
# las tablas de contingencia. Estas tablas resumen los conteos realizados
# en donde cada conteo se puede asignar a la combinación de dos factores y
# cada factor tiene al menos dos niveles. En este ejemplo usaremos tres
# tipos de epífitas creciendo en tres tipos de forofitos.
# Primero generamos la matriz de observaciones
MAT <- matrix(c(5, 20, 25, 10, 8, 50, 12, 1, 2), 3, 3)
# Y para facilitar la identidad de las filas y columnas le agregamos
# nombres a las filas y columnas
rownames(MAT) <- c("Bromelias", "Helechos", "Orqid.")
colnames(MAT) <- c("Árboles", "Palmas", "Postes")
MAT
## Árboles Palmas Postes
## Bromelias 5 10 12
## Helechos 20 8 1
## Orqid. 25 50 2
mosaicplot(MAT, shade = T)
FREC <- as.vector(MAT)
EPI <- factor(rep(c("Bromelias", "Helechos", "Orqid."), 3))
FORO <- factor(rep(c("Árboles", "Palmas", "Postes"), each = 3))
data.frame(FREC, EPI, FORO)
## FREC EPI FORO
## 1 5 Bromelias Árboles
## 2 20 Helechos Árboles
## 3 25 Orqid. Árboles
## 4 10 Bromelias Palmas
## 5 8 Helechos Palmas
## 6 50 Orqid. Palmas
## 7 12 Bromelias Postes
## 8 1 Helechos Postes
## 9 2 Orqid. Postes
M <- glm(FREC ~ EPI * FORO, poisson)
summary(M)
##
## Call:
## glm(formula = FREC ~ EPI * FORO, family = poisson)
##
## Deviance Residuals:
## [1] 0 0 0 0 0 0 0 0 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.61e+00 4.47e-01 3.60 0.00032 ***
## EPIHelechos 1.39e+00 5.00e-01 2.77 0.00556 **
## EPIOrqid. 1.61e+00 4.90e-01 3.29 0.00102 **
## FOROPalmas 6.93e-01 5.48e-01 1.27 0.20569
## FOROPostes 8.75e-01 5.32e-01 1.64 0.10003
## EPIHelechos:FOROPalmas -1.61e+00 6.89e-01 -2.34 0.01953 *
## EPIOrqid.:FOROPalmas 6.98e-15 6.00e-01 0.00 1.00000
## EPIHelechos:FOROPostes -3.87e+00 1.15e+00 -3.35 0.00080 ***
## EPIOrqid.:FOROPostes -3.40e+00 9.07e-01 -3.75 0.00018 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1.1343e+02 on 8 degrees of freedom
## Residual deviance: 5.7732e-15 on 0 degrees of freedom
## AIC: 54.18
##
## Number of Fisher Scoring iterations: 3
anova(M, test = "Chi")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: FREC
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 8 113.4
## EPI 2 33.6 6 79.8 5.0e-08 ***
## FORO 2 37.7 4 42.1 6.5e-09 ***
## EPI:FORO 4 42.1 0 0.0 1.6e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M2 <- update(M, ~. - EPI:FORO)
anova(M2, M, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: FREC ~ EPI + FORO
## Model 2: FREC ~ EPI * FORO
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 4 42.1
## 2 0 0.0 4 42.1 1.6e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(M2)
##
## Call:
## glm(formula = FREC ~ EPI + FORO, family = poisson)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8 9
## -1.794 2.465 -0.751 -1.077 -1.945 1.626 3.873 -1.474 -2.738
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3175 0.2225 10.41 < 2e-16 ***
## EPIHelechos 0.0715 0.2674 0.27 0.789
## EPIOrqid. 1.0480 0.2237 4.69 2.8e-06 ***
## FOROPalmas 0.3075 0.1863 1.65 0.099 .
## FOROPostes -1.2040 0.2944 -4.09 4.3e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 113.432 on 8 degrees of freedom
## Residual deviance: 42.115 on 4 degrees of freedom
## AIC: 88.3
##
## Number of Fisher Scoring iterations: 6
B_A <- c(1, 0, 0, 0, 0)
H_A <- c(1, 1, 0, 0, 0)
O_A <- c(1, 0, 1, 0, 0)
B_PL <- c(1, 0, 0, 1, 0)
H_PL <- c(1, 1, 0, 1, 0)
O_PL <- c(1, 0, 1, 1, 0)
B_PS <- c(1, 0, 0, 0, 1)
H_PS <- c(1, 1, 0, 0, 1)
O_PS <- c(1, 0, 1, 0, 1)
CONT <- rbind(B_A, H_A, O_A, B_PL, H_PL, O_PL, B_PS, H_PS, O_PS)
library(gmodels)
estimable(M2, CONT, beta0 = log(FREC))
## beta0 Estimate Std. Error X^2 value DF Pr(>|X^2|)
## B_A 1.6094 2.318 0.2225 10.1249 1 1.463e-03
## H_A 2.9957 2.389 0.2167 7.8392 1 5.112e-03
## O_A 3.2189 3.365 0.1596 0.8439 1 3.583e-01
## B_PL 2.3026 2.625 0.2103 2.3505 1 1.252e-01
## H_PL 2.0794 2.696 0.2041 9.1362 1 2.506e-03
## O_PL 3.9120 3.673 0.1420 2.8328 1 9.236e-02
## B_PS 2.4849 1.114 0.3101 19.5525 1 9.787e-06
## H_PS 0.0000 1.185 0.3060 14.9974 1 1.077e-04
## O_PS 0.6931 2.162 0.2686 29.8896 1 4.574e-08
ESP <- exp(estimable(M2, CONT)[, 1])
ARBOLES <- rbind(FREC[1:3], ESP[1:3])
PALMAS <- rbind(FREC[4:6], ESP[4:6])
POSTES <- rbind(FREC[7:9], ESP[7:9])
par(mfrow = c(3, 1), mar = c(2, 4, 2, 1))
barplot(ARBOLES, beside = T, main = "Árboles", ylim = c(0, 35), ylab = "Frecuencia")
text(c(2, 5, 8), c(14, 24, 33), c("**", "**", "ns"), cex = c(2, 2, 1))
legend(1, 35, c("Observado", "Esperado"), fill = c("black", "grey"), bty = "n")
barplot(PALMAS, beside = T, main = "Palmas", ylim = c(0, 60), ylab = "Frecuencia")
text(c(2, 5, 8), c(20, 20, 55), c("ns", "**", "ns"), cex = c(1, 2, 1))
barplot(POSTES, beside = T, main = "Postes", ylim = c(0, 16), ylab = "Frecuencia")
text(c(2, 5, 8), c(14, 6, 11), c("***", "***", "***"), cex = c(2, 2, 2))
par(mfrow = c(1, 1))
mosaicplot(t(MAT), shade = T)
predict(M2, type = "response")
## 1 2 3 4 5 6 7 8 9
## 10.150 10.902 28.947 13.805 14.827 39.368 3.045 3.271 8.684
mosaicplot(chisq.test(MAT), shade = TRUE)
## Warning: Chi-squared approximation may be incorrect
## Error: (list) object cannot be coerced to type 'double'
chisq.test(MAT)$expected
## Warning: Chi-squared approximation may be incorrect
## Árboles Palmas Postes
## Bromelias 10.15 13.80 3.045
## Helechos 10.90 14.83 3.271
## Orqid. 28.95 39.37 8.684
chisq.test(MAT)$residuals
## Warning: Chi-squared approximation may be incorrect
## Árboles Palmas Postes
## Bromelias -1.6166 -1.024 5.132
## Helechos 2.7553 -1.773 -1.256
## Orqid. -0.7337 1.694 -2.268
# 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/lizards.csv", header = TRUE, sep = ",")
attach(DATOS)
names(DATOS)
## [1] "n" "sun" "height" "perch" "time" "species"
M <- glm(n ~ sun * height * perch * time * species, poisson)
summary(M)
##
## Call:
## glm(formula = n ~ sun * height * perch * time * species, family = poisson)
##
## Deviance Residuals:
## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [24] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [47] 0 0
##
## Coefficients:
## Estimate
## (Intercept) 1.39e+00
## sunSun 9.16e-01
## heightLow -2.57e+01
## perchNarrow -2.88e-01
## timeMid.day -1.39e+00
## timeMorning -6.93e-01
## speciesopalinus 6.20e-15
## sunSun:heightLow 2.45e+01
## sunSun:perchNarrow 6.45e-02
## heightLow:perchNarrow 2.46e+01
## sunSun:timeMid.day 2.08e+00
## sunSun:timeMorning 7.88e-01
## heightLow:timeMid.day 1.39e+00
## heightLow:timeMorning 6.93e-01
## perchNarrow:timeMid.day 2.88e-01
## perchNarrow:timeMorning 6.93e-01
## sunSun:speciesopalinus 5.88e-01
## heightLow:speciesopalinus 2.68e+01
## perchNarrow:speciesopalinus 5.11e-01
## timeMid.day:speciesopalinus 2.08e+00
## timeMorning:speciesopalinus 2.30e+00
## sunSun:heightLow:perchNarrow -2.41e+01
## sunSun:heightLow:timeMid.day -1.79e+00
## sunSun:heightLow:timeMorning -2.78e-01
## sunSun:perchNarrow:timeMid.day 4.05e-01
## sunSun:perchNarrow:timeMorning -1.60e-01
## heightLow:perchNarrow:timeMid.day -2.46e+01
## heightLow:perchNarrow:timeMorning -2.50e+01
## sunSun:heightLow:speciesopalinus -2.59e+01
## sunSun:perchNarrow:speciesopalinus -1.10e+00
## heightLow:perchNarrow:speciesopalinus -2.73e+01
## sunSun:timeMid.day:speciesopalinus -1.43e+00
## sunSun:timeMorning:speciesopalinus -1.76e+00
## heightLow:timeMid.day:speciesopalinus -2.48e+00
## heightLow:timeMorning:speciesopalinus -2.22e+00
## perchNarrow:timeMid.day:speciesopalinus -1.20e+00
## perchNarrow:timeMorning:speciesopalinus -1.83e+00
## sunSun:heightLow:perchNarrow:timeMid.day 2.38e+01
## sunSun:heightLow:perchNarrow:timeMorning 2.26e+01
## sunSun:heightLow:perchNarrow:speciesopalinus 2.64e+01
## sunSun:heightLow:timeMid.day:speciesopalinus 2.99e+00
## sunSun:heightLow:timeMorning:speciesopalinus 2.04e+00
## sunSun:perchNarrow:timeMid.day:speciesopalinus 1.18e+00
## sunSun:perchNarrow:timeMorning:speciesopalinus 1.42e+00
## heightLow:perchNarrow:timeMid.day:speciesopalinus 1.61e+00
## heightLow:perchNarrow:timeMorning:speciesopalinus 2.78e+01
## sunSun:heightLow:perchNarrow:timeMid.day:speciesopalinus -1.31e+00
## sunSun:heightLow:perchNarrow:timeMorning:speciesopalinus -2.53e+01
## Std. Error
## (Intercept) 5.00e-01
## sunSun 5.92e-01
## heightLow 1.15e+05
## perchNarrow 7.64e-01
## timeMid.day 1.12e+00
## timeMorning 8.66e-01
## speciesopalinus 7.07e-01
## sunSun:heightLow 1.15e+05
## sunSun:perchNarrow 8.99e-01
## heightLow:perchNarrow 1.15e+05
## sunSun:timeMid.day 1.18e+00
## sunSun:timeMorning 9.70e-01
## heightLow:timeMid.day 1.62e+05
## heightLow:timeMorning 1.62e+05
## perchNarrow:timeMid.day 1.61e+00
## perchNarrow:timeMorning 1.19e+00
## sunSun:speciesopalinus 8.10e-01
## heightLow:speciesopalinus 1.15e+05
## perchNarrow:speciesopalinus 1.02e+00
## timeMid.day:speciesopalinus 1.27e+00
## timeMorning:speciesopalinus 1.02e+00
## sunSun:heightLow:perchNarrow 1.15e+05
## sunSun:heightLow:timeMid.day 1.62e+05
## sunSun:heightLow:timeMorning 1.62e+05
## sunSun:perchNarrow:timeMid.day 1.70e+00
## sunSun:perchNarrow:timeMorning 1.34e+00
## heightLow:perchNarrow:timeMid.day 1.99e+05
## heightLow:perchNarrow:timeMorning 1.99e+05
## sunSun:heightLow:speciesopalinus 1.15e+05
## sunSun:perchNarrow:speciesopalinus 1.20e+00
## heightLow:perchNarrow:speciesopalinus 1.15e+05
## sunSun:timeMid.day:speciesopalinus 1.36e+00
## sunSun:timeMorning:speciesopalinus 1.15e+00
## heightLow:timeMid.day:speciesopalinus 1.62e+05
## heightLow:timeMorning:speciesopalinus 1.62e+05
## perchNarrow:timeMid.day:speciesopalinus 1.85e+00
## perchNarrow:timeMorning:speciesopalinus 1.43e+00
## sunSun:heightLow:perchNarrow:timeMid.day 1.99e+05
## sunSun:heightLow:perchNarrow:timeMorning 1.99e+05
## sunSun:heightLow:perchNarrow:speciesopalinus 1.15e+05
## sunSun:heightLow:timeMid.day:speciesopalinus 1.62e+05
## sunSun:heightLow:timeMorning:speciesopalinus 1.62e+05
## sunSun:perchNarrow:timeMid.day:speciesopalinus 1.98e+00
## sunSun:perchNarrow:timeMorning:speciesopalinus 1.64e+00
## heightLow:perchNarrow:timeMid.day:speciesopalinus 2.30e+05
## heightLow:perchNarrow:timeMorning:speciesopalinus 1.99e+05
## sunSun:heightLow:perchNarrow:timeMid.day:speciesopalinus 2.30e+05
## sunSun:heightLow:perchNarrow:timeMorning:speciesopalinus 1.99e+05
## z value Pr(>|z|)
## (Intercept) 2.77 0.0056
## sunSun 1.55 0.1214
## heightLow 0.00 0.9998
## perchNarrow -0.38 0.7064
## timeMid.day -1.24 0.2150
## timeMorning -0.80 0.4235
## speciesopalinus 0.00 1.0000
## sunSun:heightLow 0.00 0.9998
## sunSun:perchNarrow 0.07 0.9428
## heightLow:perchNarrow 0.00 0.9998
## sunSun:timeMid.day 1.76 0.0788
## sunSun:timeMorning 0.81 0.4163
## heightLow:timeMid.day 0.00 1.0000
## heightLow:timeMorning 0.00 1.0000
## perchNarrow:timeMid.day 0.18 0.8579
## perchNarrow:timeMorning 0.58 0.5603
## sunSun:speciesopalinus 0.73 0.4679
## heightLow:speciesopalinus 0.00 0.9998
## perchNarrow:speciesopalinus 0.50 0.6153
## timeMid.day:speciesopalinus 1.63 0.1028
## timeMorning:speciesopalinus 2.25 0.0246
## sunSun:heightLow:perchNarrow 0.00 0.9998
## sunSun:heightLow:timeMid.day 0.00 1.0000
## sunSun:heightLow:timeMorning 0.00 1.0000
## sunSun:perchNarrow:timeMid.day 0.24 0.8115
## sunSun:perchNarrow:timeMorning -0.12 0.9051
## heightLow:perchNarrow:timeMid.day 0.00 0.9999
## heightLow:perchNarrow:timeMorning 0.00 0.9999
## sunSun:heightLow:speciesopalinus 0.00 0.9998
## sunSun:perchNarrow:speciesopalinus -0.92 0.3597
## heightLow:perchNarrow:speciesopalinus 0.00 0.9998
## sunSun:timeMid.day:speciesopalinus -1.05 0.2928
## sunSun:timeMorning:speciesopalinus -1.53 0.1260
## heightLow:timeMid.day:speciesopalinus 0.00 1.0000
## heightLow:timeMorning:speciesopalinus 0.00 1.0000
## perchNarrow:timeMid.day:speciesopalinus -0.65 0.5143
## perchNarrow:timeMorning:speciesopalinus -1.28 0.1997
## sunSun:heightLow:perchNarrow:timeMid.day 0.00 0.9999
## sunSun:heightLow:perchNarrow:timeMorning 0.00 0.9999
## sunSun:heightLow:perchNarrow:speciesopalinus 0.00 0.9998
## sunSun:heightLow:timeMid.day:speciesopalinus 0.00 1.0000
## sunSun:heightLow:timeMorning:speciesopalinus 0.00 1.0000
## sunSun:perchNarrow:timeMid.day:speciesopalinus 0.60 0.5508
## sunSun:perchNarrow:timeMorning:speciesopalinus 0.86 0.3879
## heightLow:perchNarrow:timeMid.day:speciesopalinus 0.00 1.0000
## heightLow:perchNarrow:timeMorning:speciesopalinus 0.00 0.9999
## sunSun:heightLow:perchNarrow:timeMid.day:speciesopalinus 0.00 1.0000
## sunSun:heightLow:perchNarrow:timeMorning:speciesopalinus 0.00 0.9999
##
## (Intercept) **
## sunSun
## heightLow
## perchNarrow
## timeMid.day
## timeMorning
## speciesopalinus
## sunSun:heightLow
## sunSun:perchNarrow
## heightLow:perchNarrow
## sunSun:timeMid.day .
## sunSun:timeMorning
## heightLow:timeMid.day
## heightLow:timeMorning
## perchNarrow:timeMid.day
## perchNarrow:timeMorning
## sunSun:speciesopalinus
## heightLow:speciesopalinus
## perchNarrow:speciesopalinus
## timeMid.day:speciesopalinus
## timeMorning:speciesopalinus *
## sunSun:heightLow:perchNarrow
## sunSun:heightLow:timeMid.day
## sunSun:heightLow:timeMorning
## sunSun:perchNarrow:timeMid.day
## sunSun:perchNarrow:timeMorning
## heightLow:perchNarrow:timeMid.day
## heightLow:perchNarrow:timeMorning
## sunSun:heightLow:speciesopalinus
## sunSun:perchNarrow:speciesopalinus
## heightLow:perchNarrow:speciesopalinus
## sunSun:timeMid.day:speciesopalinus
## sunSun:timeMorning:speciesopalinus
## heightLow:timeMid.day:speciesopalinus
## heightLow:timeMorning:speciesopalinus
## perchNarrow:timeMid.day:speciesopalinus
## perchNarrow:timeMorning:speciesopalinus
## sunSun:heightLow:perchNarrow:timeMid.day
## sunSun:heightLow:perchNarrow:timeMorning
## sunSun:heightLow:perchNarrow:speciesopalinus
## sunSun:heightLow:timeMid.day:speciesopalinus
## sunSun:heightLow:timeMorning:speciesopalinus
## sunSun:perchNarrow:timeMid.day:speciesopalinus
## sunSun:perchNarrow:timeMorning:speciesopalinus
## heightLow:perchNarrow:timeMid.day:speciesopalinus
## heightLow:perchNarrow:timeMorning:speciesopalinus
## sunSun:heightLow:perchNarrow:timeMid.day:speciesopalinus
## sunSun:heightLow:perchNarrow:timeMorning:speciesopalinus
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 7.3756e+02 on 47 degrees of freedom
## Residual deviance: 3.3472e-10 on 0 degrees of freedom
## AIC: 259.2
##
## Number of Fisher Scoring iterations: 22
anova(M)
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 47 738
## sun 1 242.7 46 495
## height 1 49.6 45 445
## perch 1 28.4 44 417
## time 2 98.5 42 318
## species 1 165.7 41 153
## sun:height 1 0.9 40 152
## sun:perch 1 3.6 39 148
## height:perch 1 16.1 38 132
## sun:time 2 48.0 36 84
## height:time 2 1.7 34 82
## perch:time 2 1.6 32 81
## sun:species 1 6.5 31 74
## height:species 1 26.1 30 48
## perch:species 1 11.6 29 37
## time:species 2 11.5 27 25
## sun:height:perch 1 0.3 26 25
## sun:height:time 2 0.9 24 24
## sun:perch:time 2 2.7 22 21
## height:perch:time 2 1.8 20 19
## sun:height:species 1 2.3 19 17
## sun:perch:species 1 0.0 18 17
## height:perch:species 1 0.3 17 17
## sun:time:species 2 2.0 15 15
## height:time:species 2 1.6 13 13
## perch:time:species 2 0.0 11 13
## sun:height:perch:time 2 4.7 9 9
## sun:height:perch:species 1 2.8 8 6
## sun:height:time:species 2 0.9 6 5
## sun:perch:time:species 2 1.7 4 3
## height:perch:time:species 2 3.2 2 0
## sun:height:perch:time:species 2 0.0 0 0
M2 <- update(M, ~. - sun:height:perch:time:species)
## Warning: glm.fit: fitted rates numerically 0 occurred
anova(M2, M, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species +
## height:perch:time:species
## Model 2: n ~ sun * height * perch * time * species
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 2 2.18e-10
## 2 0 3.35e-10 2 -1.17e-10
anova(M2)
## Warning: glm.fit: fitted rates numerically 0 occurred
## Warning: glm.fit: fitted rates numerically 0 occurred
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 47 738
## sun 1 242.7 46 495
## height 1 49.6 45 445
## perch 1 28.4 44 417
## time 2 98.5 42 318
## species 1 165.7 41 153
## sun:height 1 0.9 40 152
## sun:perch 1 3.6 39 148
## height:perch 1 16.1 38 132
## sun:time 2 48.0 36 84
## height:time 2 1.7 34 82
## perch:time 2 1.6 32 81
## sun:species 1 6.5 31 74
## height:species 1 26.1 30 48
## perch:species 1 11.6 29 37
## time:species 2 11.5 27 25
## sun:height:perch 1 0.3 26 25
## sun:height:time 2 0.9 24 24
## sun:perch:time 2 2.7 22 21
## height:perch:time 2 1.8 20 19
## sun:height:species 1 2.3 19 17
## sun:perch:species 1 0.0 18 17
## height:perch:species 1 0.3 17 17
## sun:time:species 2 2.0 15 15
## height:time:species 2 1.6 13 13
## perch:time:species 2 0.0 11 13
## sun:height:perch:time 2 4.7 9 9
## sun:height:perch:species 1 2.8 8 6
## sun:height:time:species 2 0.9 6 5
## sun:perch:time:species 2 1.7 4 3
## height:perch:time:species 2 3.2 2 0
M3 <- update(M2, ~. - height:perch:time:species)
## Warning: glm.fit: fitted rates numerically 0 occurred
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species +
## height:perch:time:species
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 4 3.22
## 2 2 0.00 2 3.22 0.2
M3 <- update(M2, ~. - sun:perch:time:species)
## Warning: glm.fit: fitted rates numerically 0 occurred
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + height:perch:time:species
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species +
## height:perch:time:species
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 4 0.81
## 2 2 0.00 2 0.81 0.67
M3 <- update(M2, ~. - sun:height:time:species)
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:perch:time:species + height:perch:time:species
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species +
## height:perch:time:species
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 4 0.442
## 2 2 0.000 2 0.442 0.8
M3 <- update(M2, ~. - sun:height:perch:species)
## Warning: glm.fit: fitted rates numerically 0 occurred
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:time:species + sun:perch:time:species + height:perch:time:species
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species +
## height:perch:time:species
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 3 2.71
## 2 2 0.00 1 2.71 0.1 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M3 <- update(M2, ~. - height:perch:time:species - sun:perch:time:species - sun:height:time:species -
sun:height:perch:species)
anova(M3, M2, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time +
## sun:height:perch:species + sun:height:time:species + sun:perch:time:species +
## height:perch:time:species
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 9 8.57
## 2 2 0.00 7 8.57 0.28
anova(M3)
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 47 738
## sun 1 242.7 46 495
## height 1 49.6 45 445
## perch 1 28.4 44 417
## time 2 98.5 42 318
## species 1 165.7 41 153
## sun:height 1 0.9 40 152
## sun:perch 1 3.6 39 148
## height:perch 1 16.1 38 132
## sun:time 2 48.0 36 84
## height:time 2 1.7 34 82
## perch:time 2 1.6 32 81
## sun:species 1 6.5 31 74
## height:species 1 26.1 30 48
## perch:species 1 11.6 29 37
## time:species 2 11.5 27 25
## sun:height:perch 1 0.3 26 25
## sun:height:time 2 0.9 24 24
## sun:perch:time 2 2.7 22 21
## height:perch:time 2 1.8 20 19
## sun:height:species 1 2.3 19 17
## sun:perch:species 1 0.0 18 17
## height:perch:species 1 0.3 17 17
## sun:time:species 2 2.0 15 15
## height:time:species 2 1.6 13 13
## perch:time:species 2 0.0 11 13
## sun:height:perch:time 2 4.7 9 9
M4 <- update(M3, ~. - perch:time:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 8.59
## 2 9 8.57 2 0.0182 0.99
M4 <- update(M3, ~. - height:time:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## perch:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 9.97
## 2 9 8.57 2 1.4 0.5
M4 <- update(M3, ~. - sun:time:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + height:time:species +
## perch:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 11 10.90
## 2 9 8.57 2 2.33 0.31
M4 <- update(M3, ~. - height:perch:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + sun:time:species + height:time:species +
## perch:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 10 8.88
## 2 9 8.57 1 0.302 0.58
M4 <- update(M3, ~. - sun:perch:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## height:perch:species + sun:time:species + height:time:species +
## perch:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 10 8.59
## 2 9 8.57 1 0.0143 0.9
M4 <- update(M3, ~. - sun:height:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:perch:species +
## height:perch:species + sun:time:species + height:time:species +
## perch:time:species + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 10 11.87
## 2 9 8.57 1 3.29 0.07 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M4 <- update(M3, ~. - perch:time:species - height:time:species - sun:time:species -
height:perch:species - sun:perch:species - sun:height:species)
anova(M4, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 18 14.20
## 2 9 8.57 9 5.63 0.78
anova(M4)
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 47 738
## sun 1 242.7 46 495
## height 1 49.6 45 445
## perch 1 28.4 44 417
## time 2 98.5 42 318
## species 1 165.7 41 153
## sun:height 1 0.9 40 152
## sun:perch 1 3.6 39 148
## height:perch 1 16.1 38 132
## sun:time 2 48.0 36 84
## height:time 2 1.7 34 82
## perch:time 2 1.6 32 81
## sun:species 1 6.5 31 74
## height:species 1 26.1 30 48
## perch:species 1 11.6 29 37
## time:species 2 11.5 27 25
## sun:height:perch 1 0.3 26 25
## sun:height:time 2 0.9 24 24
## sun:perch:time 2 2.7 22 21
## height:perch:time 2 1.8 20 19
## sun:height:perch:time 2 5.3 18 14
M5 <- update(M4, ~. - time:species)
anova(M5, M4, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + sun:height:perch + sun:height:time +
## sun:perch:time + height:perch:time + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 20 25.8
## 2 18 14.2 2 11.6 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M5 <- update(M4, ~. - perch:species)
anova(M5, M4, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + time:species + sun:height:perch + sun:height:time +
## sun:perch:time + height:perch:time + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 19 27.3
## 2 18 14.2 1 13.1 0.00029 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M5 <- update(M4, ~. - height:species)
anova(M5, M3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## perch:species + time:species + sun:height:perch + sun:height:time +
## sun:perch:time + height:perch:time + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:species +
## sun:perch:species + height:perch:species + sun:time:species +
## height:time:species + perch:time:species + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 19 36.3
## 2 9 8.6 10 27.7 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
M5 <- update(M4, ~. - sun:species)
anova(M5, M4, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + height:species +
## perch:species + time:species + sun:height:perch + sun:height:time +
## sun:perch:time + height:perch:time + sun:height:perch:time
## Model 2: n ~ sun + height + perch + time + species + sun:height + sun:perch +
## height:perch + sun:time + height:time + perch:time + sun:species +
## height:species + perch:species + time:species + sun:height:perch +
## sun:height:time + sun:perch:time + height:perch:time + sun:height:perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 19 21.9
## 2 18 14.2 1 7.69 0.0056 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(M4)
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev
## NULL 47 738
## sun 1 242.7 46 495
## height 1 49.6 45 445
## perch 1 28.4 44 417
## time 2 98.5 42 318
## species 1 165.7 41 153
## sun:height 1 0.9 40 152
## sun:perch 1 3.6 39 148
## height:perch 1 16.1 38 132
## sun:time 2 48.0 36 84
## height:time 2 1.7 34 82
## perch:time 2 1.6 32 81
## sun:species 1 6.5 31 74
## height:species 1 26.1 30 48
## perch:species 1 11.6 29 37
## time:species 2 11.5 27 25
## sun:height:perch 1 0.3 26 25
## sun:height:time 2 0.9 24 24
## sun:perch:time 2 2.7 22 21
## height:perch:time 2 1.8 20 19
## sun:height:perch:time 2 5.3 18 14
summary(M4)
##
## Call:
## glm(formula = n ~ sun + height + perch + time + species + sun:height +
## sun:perch + height:perch + sun:time + height:time + perch:time +
## sun:species + height:species + perch:species + time:species +
## sun:height:perch + sun:height:time + sun:perch:time + height:perch:time +
## sun:height:perch:time, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4536 -0.3348 -0.0002 0.2343 1.3751
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 0.610 0.446 1.37
## sunSun 1.832 0.467 3.93
## heightLow -0.555 0.511 -1.09
## perchNarrow 0.529 0.523 1.01
## timeMid.day -0.693 0.542 -1.28
## timeMorning 0.403 0.481 0.84
## speciesopalinus 1.208 0.354 3.42
## sunSun:heightLow -0.809 0.557 -1.45
## sunSun:perchNarrow -0.711 0.594 -1.20
## heightLow:perchNarrow -1.654 0.914 -1.81
## sunSun:timeMid.day 1.178 0.537 2.20
## sunSun:timeMorning -0.424 0.482 -0.88
## heightLow:timeMid.day -0.621 0.668 -0.93
## heightLow:timeMorning -1.012 0.577 -1.76
## perchNarrow:timeMid.day -0.465 0.751 -0.62
## perchNarrow:timeMorning -0.593 0.624 -0.95
## sunSun:speciesopalinus -0.847 0.322 -2.63
## heightLow:speciesopalinus 1.130 0.257 4.40
## perchNarrow:speciesopalinus -0.763 0.211 -3.61
## timeMid.day:speciesopalinus 0.964 0.282 3.42
## timeMorning:speciesopalinus 0.737 0.299 2.46
## sunSun:heightLow:perchNarrow 1.715 1.058 1.62
## sunSun:heightLow:timeMid.day 0.597 0.757 0.79
## sunSun:heightLow:timeMorning 1.209 0.693 1.74
## sunSun:perchNarrow:timeMid.day 1.227 0.826 1.49
## sunSun:perchNarrow:timeMorning 0.945 0.734 1.29
## heightLow:perchNarrow:timeMid.day -16.852 3111.564 -0.01
## heightLow:perchNarrow:timeMorning 1.659 1.102 1.51
## sunSun:heightLow:perchNarrow:timeMid.day 16.068 3111.564 0.01
## sunSun:heightLow:perchNarrow:timeMorning -2.251 1.287 -1.75
## Pr(>|z|)
## (Intercept) 0.17153
## sunSun 8.6e-05 ***
## heightLow 0.27711
## perchNarrow 0.31185
## timeMid.day 0.20144
## timeMorning 0.40254
## speciesopalinus 0.00063 ***
## sunSun:heightLow 0.14696
## sunSun:perchNarrow 0.23123
## heightLow:perchNarrow 0.07046 .
## sunSun:timeMid.day 0.02811 *
## sunSun:timeMorning 0.37911
## heightLow:timeMid.day 0.35273
## heightLow:timeMorning 0.07915 .
## perchNarrow:timeMid.day 0.53603
## perchNarrow:timeMorning 0.34158
## sunSun:speciesopalinus 0.00858 **
## heightLow:speciesopalinus 1.1e-05 ***
## perchNarrow:speciesopalinus 0.00031 ***
## timeMid.day:speciesopalinus 0.00062 ***
## timeMorning:speciesopalinus 0.01373 *
## sunSun:heightLow:perchNarrow 0.10513
## sunSun:heightLow:timeMid.day 0.43024
## sunSun:heightLow:timeMorning 0.08113 .
## sunSun:perchNarrow:timeMid.day 0.13726
## sunSun:perchNarrow:timeMorning 0.19760
## heightLow:perchNarrow:timeMid.day 0.99568
## heightLow:perchNarrow:timeMorning 0.13227
## sunSun:heightLow:perchNarrow:timeMid.day 0.99588
## sunSun:heightLow:perchNarrow:timeMorning 0.08039 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 737.555 on 47 degrees of freedom
## Residual deviance: 14.205 on 18 degrees of freedom
## AIC: 237.5
##
## Number of Fisher Scoring iterations: 16
OE <- rbind(n, predict(M5, type = "response"))
barplot(OE, beside = T)
MM <- M <- glm(n ~ sun + height + perch + time + species, poisson)
ETIQ <- character()
for (i in 1:48) {
for (j in 2:6) {
ETIQ[i] <- paste(ETIQ[i], substr(DATOS[i, j], 1, 2), sep = "")
}
}
write.table(ETIQ, "~/desktop/curso R 2012/ETIQ.csv", sep = ",")
CONT <- read.table("~/desktop/curso R 2012/CONTRLAGARTAS.csv", sep = ",", header = T,
row.names = 1)
colnames(CONT) <- rownames(summary(MM)$coefficients)
V <- as.vector(tapply(n, list(species, sun), sum))
estimable(MM, CONT[49:52, ], beta0 = log(V[c(4, 2, 3, 1)]))
## beta0 Estimate Std. Error X^2 value DF Pr(>|X^2|)
## Suop 5.835 3.2133 0.05791 2049.6 1 0
## Shop 4.489 1.7264 0.10390 706.8 1 0
## Sugr 4.771 2.0375 0.09248 873.5 1 0
## Shgr 2.708 0.5507 0.12646 291.0 1 0
V <- as.vector(tapply(n, list(species, height), sum))
estimable(MM, CONT[53:56, ], beta0 = log(V[c(4, 2, 3, 1)]))
## beta0 Estimate Std. Error X^2 value DF Pr(>|X^2|)
## Loop 5.170 2.1665 0.08491 1251.5 1 0
## Hiop 5.541 2.7731 0.07018 1555.8 1 0
## Logr 3.135 0.9908 0.11140 370.7 1 0
## Higr 4.700 1.5974 0.10062 951.2 1 0
V <- as.vector(tapply(n, list(species, perch), sum))
estimable(MM, CONT[57:60, ], beta0 = log(V[c(2, 4, 1, 3)]))
## beta0 Estimate Std. Error X^2 value DF Pr(>|X^2|)
## Brop 5.652 2.697 0.07183 1693.1 1 0
## Naop 4.984 2.243 0.08262 1100.6 1 0
## Brgr 4.094 1.521 0.10177 639.2 1 0
## Nrgr 4.290 1.067 0.10966 864.2 1 0
V <- as.vector(tapply(n, list(species, time), sum))
estimable(MM, CONT[61:66, ], beta0 = log(V[c(2, 4, 6, 1, 3, 5)]))
## beta0 Estimate Std. Error X^2 value DF Pr(>|X^2|)
## Afop 4.174 1.9116 0.11039 420.2 1 0
## Miop 5.416 2.9896 0.07386 1079.4 1 0
## Moop 4.949 2.5084 0.08718 783.5 1 0
## Afgr 3.497 0.7358 0.13185 438.4 1 0
## Migr 4.143 1.8138 0.10322 509.3 1 0
## Mogr 3.611 1.3326 0.11313 405.5 1 0