Análise exploratória
library(MRDCr)
# help(soja)
ls("package:MRDCr")
## [1] "apc" "calc_mean_cmp"
## [3] "calc_mean_gcnt" "calc_var_cmp"
## [5] "cambras" "capdesfo"
## [7] "capmosca" "cholV_eta"
## [9] "cmp" "conftemp"
## [11] "confterm" "convergencez"
## [13] "dcmp" "dgcnt"
## [15] "dpgnz" "gcnt"
## [17] "led" "llcmp"
## [19] "llgcnt" "llpgnz"
## [21] "nematoide" "ninfas"
## [23] "panel.beeswarm" "panel.cbH"
## [25] "panel.groups.segplot" "peixe"
## [27] "pgnz" "postura"
## [29] "predict.mle2" "prepanel.cbH"
## [31] "seguros" "soja"
library(lattice)
str(soja)
## 'data.frame': 75 obs. of 5 variables:
## $ K : int 0 30 60 120 180 0 30 60 120 180 ...
## $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 1 1 1 1 2 2 2 2 2 ...
## $ bloc: Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ngra: int 136 159 156 171 190 140 193 200 208 237 ...
## $ nvag: int 56 62 66 68 82 63 86 94 86 97 ...
xtabs(~umid + K, data = soja[-75, ])
## K
## umid 0 30 60 120 180
## 37,5 5 5 5 5 5
## 50 5 5 5 5 5
## 62,5 5 5 5 5 4
xyplot(nvag + ngra ~ K, groups = umid, outer = TRUE,
data = soja[-74, ],
type = c("p", "a"), scales = "free",
ylab = NULL,
xlab = expression("Nível de potássio aplicado"~(mg~dm^{-3})),
auto.key = list(title = "Umidade do solo (%)",
cex.title = 1,
columns = 3),
strip = strip.custom(
factor.levels = c("Número de vagens",
"Número de grãos")))
soja <- soja[-74, ]
GLM Poisson para número de vagens viáveis por vaso
# Considerar K categórico.
soja <- transform(soja, K = factor(K))
m0 <- glm(nvag ~ bloc + umid * K, data = soja, family = poisson)
par(mfrow = c(2, 2))
plot(m0); layout(1)
#-----------------------------------------------------------------------
deviance(m0)
## [1] 65.27834
df.residual(m0)
## [1] 55
summary(m0)
##
## Call:
## glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9681 -0.6393 -0.0195 0.4854 2.3306
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.95369 0.06892 57.363 < 2e-16 ***
## blocII -0.02928 0.04091 -0.716 0.47414
## blocIII -0.07265 0.04136 -1.756 0.07902 .
## blocIV -0.12544 0.04194 -2.991 0.00278 **
## blocV -0.10795 0.04297 -2.512 0.01199 *
## umid50 0.13404 0.08765 1.529 0.12619
## umid62,5 0.21656 0.08602 2.518 0.01181 *
## K30 0.27427 0.08493 3.229 0.00124 **
## K60 0.30797 0.08432 3.652 0.00026 ***
## K120 0.32883 0.08395 3.917 8.97e-05 ***
## K180 0.25540 0.08528 2.995 0.00275 **
## umid50:K30 0.06322 0.11557 0.547 0.58433
## umid62,5:K30 -0.10747 0.11536 -0.932 0.35152
## umid50:K60 0.16561 0.11370 1.457 0.14521
## umid62,5:K60 0.10735 0.11220 0.957 0.33869
## umid50:K120 0.14920 0.11338 1.316 0.18818
## umid62,5:K120 0.11839 0.11404 1.038 0.29922
## umid50:K180 0.30370 0.11361 2.673 0.00751 **
## umid62,5:K180 0.19838 0.11256 1.762 0.07800 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 322.684 on 73 degrees of freedom
## Residual deviance: 65.278 on 55 degrees of freedom
## AIC: 557.23
##
## Number of Fisher Scoring iterations: 4
anova(m0, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: nvag
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 73 322.68
## bloc 4 14.284 69 308.40 0.006442 **
## umid 2 92.911 67 215.49 < 2.2e-16 ***
## K 4 136.060 63 79.43 < 2.2e-16 ***
## umid:K 8 14.150 55 65.28 0.077926 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Predição com intervalos de confiança.
library(doBy)
library(multcomp)
X <- LSmatrix(m0, effect = c("umid", "K"))
pred <- attr(X, "grid")
pred <- transform(pred,
K = as.integer(as.character(K)),
umid = factor(umid))
# Quantil normal para fazer um IC de 95%.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
# Preditos pela Poisson.
# aux <- predict(m0, newdata = pred$pois, se.fit = TRUE)
# aux <- exp(aux$fit + outer(aux$se.fit, qn, FUN = "*"))
# pred$pois <- cbind(pred$pois, aux)
aux <- confint(glht(m0, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred <- cbind(pred, exp(aux))
str(pred)
## 'data.frame': 15 obs. of 5 variables:
## $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 2 3 1 2 3 1 2 3 1 ...
## $ K : int 0 0 0 30 30 30 60 60 60 120 ...
## $ fit : num 48.7 55.7 60.5 64.1 78.1 ...
## $ lwr : num 43 49.6 54.1 57.5 70.7 ...
## $ upr : num 55.3 62.7 67.7 71.5 86.3 ...
xyplot(fit ~ K | umid, data = pred,
layout = c(NA, 1), as.table = TRUE,
ylab = "Número de vagens por vaso",
xlab = expression("Nível de potássio aplicado"~(mg~dm^{-3})),
xlim = extendrange(range(pred$K), f = 0.1),
ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
prepanel = prepanel.cbH,
panel = panel.cbH)
#-----------------------------------------------------------------------
# Comparações múltiplas.
L <- by(X, INDICES = pred$umid, FUN = as.matrix)
names(L) <- levels(soja$umid)
L <- lapply(L, "rownames<-", levels(soja$K))
K <- lapply(L, apc)
apc(L[[1]])
## (Intercept) blocII blocIII blocIV blocV umid50 umid62,5 K30
## 0-30 0 0 0 0 0 0 0 -1
## 0-60 0 0 0 0 0 0 0 0
## 0-120 0 0 0 0 0 0 0 0
## 0-180 0 0 0 0 0 0 0 0
## 30-60 0 0 0 0 0 0 0 1
## 30-120 0 0 0 0 0 0 0 1
## 30-180 0 0 0 0 0 0 0 1
## 60-120 0 0 0 0 0 0 0 0
## 60-180 0 0 0 0 0 0 0 0
## 120-180 0 0 0 0 0 0 0 0
## K60 K120 K180 umid50:K30 umid62,5:K30 umid50:K60
## 0-30 0 0 0 0 0 0
## 0-60 -1 0 0 0 0 0
## 0-120 0 -1 0 0 0 0
## 0-180 0 0 -1 0 0 0
## 30-60 -1 0 0 0 0 0
## 30-120 0 -1 0 0 0 0
## 30-180 0 0 -1 0 0 0
## 60-120 1 -1 0 0 0 0
## 60-180 1 0 -1 0 0 0
## 120-180 0 1 -1 0 0 0
## umid62,5:K60 umid50:K120 umid62,5:K120 umid50:K180
## 0-30 0 0 0 0
## 0-60 0 0 0 0
## 0-120 0 0 0 0
## 0-180 0 0 0 0
## 30-60 0 0 0 0
## 30-120 0 0 0 0
## 30-180 0 0 0 0
## 60-120 0 0 0 0
## 60-180 0 0 0 0
## 120-180 0 0 0 0
## umid62,5:K180
## 0-30 0
## 0-60 0
## 0-120 0
## 0-180 0
## 30-60 0
## 30-120 0
## 30-180 0
## 60-120 0
## 60-180 0
## 120-180 0
lapply(K,
FUN = function(k) {
summary(glht(model = m0, linfct = k),
test = adjusted(type = "fdr"))
})
## $`37,5`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 0-30 == 0 -0.27427 0.08493 -3.229 0.004137 **
## 0-60 == 0 -0.30797 0.08432 -3.652 0.001300 **
## 0-120 == 0 -0.32883 0.08395 -3.917 0.000897 ***
## 0-180 == 0 -0.25540 0.08528 -2.995 0.006865 **
## 30-60 == 0 -0.03369 0.07828 -0.430 0.811949
## 30-120 == 0 -0.05456 0.07788 -0.701 0.719952
## 30-180 == 0 0.01887 0.07931 0.238 0.811949
## 60-120 == 0 -0.02087 0.07721 -0.270 0.811949
## 60-180 == 0 0.05256 0.07866 0.668 0.719952
## 120-180 == 0 0.07343 0.07826 0.938 0.696218
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
##
##
## $`50`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 0-30 == 0 -0.337496 0.078369 -4.306 4.15e-05 ***
## 0-60 == 0 -0.473581 0.076265 -6.210 1.77e-09 ***
## 0-120 == 0 -0.478036 0.076200 -6.273 1.77e-09 ***
## 0-180 == 0 -0.559104 0.075056 -7.449 9.39e-13 ***
## 30-60 == 0 -0.136086 0.069208 -1.966 0.07037 .
## 30-120 == 0 -0.140540 0.069136 -2.033 0.07012 .
## 30-180 == 0 -0.221608 0.067873 -3.265 0.00219 **
## 60-120 == 0 -0.004454 0.066741 -0.067 0.94679
## 60-180 == 0 -0.085522 0.065432 -1.307 0.23870
## 120-180 == 0 -0.081068 0.065356 -1.240 0.23870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
##
##
## $`62,5`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = nvag ~ bloc + umid * K, family = poisson, data = soja)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 0-30 == 0 -0.166800 0.078062 -2.137 0.046595 *
## 0-60 == 0 -0.415317 0.074020 -5.611 6.71e-08 ***
## 0-120 == 0 -0.447218 0.077181 -5.794 3.43e-08 ***
## 0-180 == 0 -0.453784 0.073463 -6.177 6.53e-09 ***
## 30-60 == 0 -0.248517 0.070512 -3.524 0.000707 ***
## 30-120 == 0 -0.280417 0.073823 -3.799 0.000291 ***
## 30-180 == 0 -0.286984 0.069927 -4.104 0.000101 ***
## 60-120 == 0 -0.031900 0.069536 -0.459 0.718229
## 60-180 == 0 -0.038466 0.065384 -0.588 0.695404
## 120-180 == 0 -0.006566 0.068942 -0.095 0.924124
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
GLM Poisson para número de grãos por vaso
m1 <- glm(ngra ~ bloc + umid * K, data = soja, family = poisson)
par(mfrow = c(2, 2))
plot(m1); layout(1)
#-----------------------------------------------------------------------
deviance(m1)
## [1] 124.7194
df.residual(m1)
## [1] 55
summary(m1)
##
## Call:
## glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6297 -1.0707 -0.0873 0.7428 3.2810
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.80196 0.04477 107.254 < 2e-16 ***
## blocII -0.01939 0.02655 -0.730 0.465260
## blocIII -0.03663 0.02667 -1.373 0.169673
## blocIV -0.10559 0.02715 -3.889 0.000101 ***
## blocV -0.09313 0.02788 -3.340 0.000837 ***
## umid50 0.13245 0.05692 2.327 0.019968 *
## umid62,5 0.18548 0.05623 3.299 0.000972 ***
## K30 0.29799 0.05486 5.432 5.56e-08 ***
## K60 0.34434 0.05432 6.339 2.32e-10 ***
## K120 0.36493 0.05409 6.746 1.52e-11 ***
## K180 0.29542 0.05489 5.383 7.35e-08 ***
## umid50:K30 0.04344 0.07482 0.581 0.561495
## umid62,5:K30 -0.13669 0.07527 -1.816 0.069349 .
## umid50:K60 0.11559 0.07361 1.570 0.116354
## umid62,5:K60 0.09174 0.07289 1.259 0.208190
## umid50:K120 0.11860 0.07329 1.618 0.105637
## umid62,5:K120 0.16266 0.07367 2.208 0.027242 *
## umid50:K180 0.28832 0.07327 3.935 8.33e-05 ***
## umid62,5:K180 0.21569 0.07285 2.961 0.003071 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 761.43 on 73 degrees of freedom
## Residual deviance: 124.72 on 55 degrees of freedom
## AIC: 681.34
##
## Number of Fisher Scoring iterations: 4
anova(m1, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: ngra
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 73 761.43
## bloc 4 28.34 69 733.09 1.063e-05 ***
## umid 2 185.06 67 548.03 < 2.2e-16 ***
## K 4 380.32 63 167.71 < 2.2e-16 ***
## umid:K 8 42.99 55 124.72 8.830e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#-----------------------------------------------------------------------
# Predição com intervalos de confiança.
X <- LSmatrix(m1, effect = c("umid", "K"))
pred <- attr(X, "grid")
pred <- transform(pred,
K = as.integer(as.character(K)),
umid = factor(umid))
# Quantil normal para fazer um IC de 95%.
qn <- qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1)
aux <- confint(glht(m1, linfct = X),
calpha = univariate_calpha())$confint
colnames(aux)[1] <- "fit"
pred <- cbind(pred, exp(aux))
str(pred)
## 'data.frame': 15 obs. of 5 variables:
## $ umid: Factor w/ 3 levels "37,5","50","62,5": 1 2 3 1 2 3 1 2 3 1 ...
## $ K : int 0 0 0 30 30 30 60 60 60 120 ...
## $ fit : num 116 132 139 156 186 ...
## $ lwr : num 107 122 129 145 174 ...
## $ upr : num 126 143 150 167 198 ...
xyplot(fit ~ K | umid, data = pred,
layout = c(NA, 1), as.table = TRUE,
ylab = "Número de grãos por vaso",
xlab = expression("Nível de potássio aplicado"~(mg~dm^{-3})),
xlim = extendrange(range(pred$K), f = 0.1),
ly = pred$lwr, uy = pred$upr, cty = "bars", length = 0,
prepanel = prepanel.cbH,
panel = panel.cbH)
#-----------------------------------------------------------------------
# Comparações múltiplas.
L <- by(X, INDICES = pred$umid, FUN = as.matrix)
names(L) <- levels(soja$umid)
L <- lapply(L, "rownames<-", levels(soja$K))
K <- lapply(L, apc)
apc(L[[1]])
## (Intercept) blocII blocIII blocIV blocV umid50 umid62,5 K30
## 0-30 0 0 0 0 0 0 0 -1
## 0-60 0 0 0 0 0 0 0 0
## 0-120 0 0 0 0 0 0 0 0
## 0-180 0 0 0 0 0 0 0 0
## 30-60 0 0 0 0 0 0 0 1
## 30-120 0 0 0 0 0 0 0 1
## 30-180 0 0 0 0 0 0 0 1
## 60-120 0 0 0 0 0 0 0 0
## 60-180 0 0 0 0 0 0 0 0
## 120-180 0 0 0 0 0 0 0 0
## K60 K120 K180 umid50:K30 umid62,5:K30 umid50:K60
## 0-30 0 0 0 0 0 0
## 0-60 -1 0 0 0 0 0
## 0-120 0 -1 0 0 0 0
## 0-180 0 0 -1 0 0 0
## 30-60 -1 0 0 0 0 0
## 30-120 0 -1 0 0 0 0
## 30-180 0 0 -1 0 0 0
## 60-120 1 -1 0 0 0 0
## 60-180 1 0 -1 0 0 0
## 120-180 0 1 -1 0 0 0
## umid62,5:K60 umid50:K120 umid62,5:K120 umid50:K180
## 0-30 0 0 0 0
## 0-60 0 0 0 0
## 0-120 0 0 0 0
## 0-180 0 0 0 0
## 30-60 0 0 0 0
## 30-120 0 0 0 0
## 30-180 0 0 0 0
## 60-120 0 0 0 0
## 60-180 0 0 0 0
## 120-180 0 0 0 0
## umid62,5:K180
## 0-30 0
## 0-60 0
## 0-120 0
## 0-180 0
## 30-60 0
## 30-120 0
## 30-180 0
## 60-120 0
## 60-180 0
## 120-180 0
lapply(K,
FUN = function(k) {
summary(glht(model = m1, linfct = k),
test = adjusted(type = "fdr"))
})
## $`37,5`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 0-30 == 0 -0.297991 0.054856 -5.432 1.84e-07 ***
## 0-60 == 0 -0.344337 0.054324 -6.339 1.16e-09 ***
## 0-120 == 0 -0.364931 0.054094 -6.746 1.52e-10 ***
## 0-180 == 0 -0.295424 0.054886 -5.383 1.84e-07 ***
## 30-60 == 0 -0.046345 0.050060 -0.926 0.443
## 30-120 == 0 -0.066939 0.049811 -1.344 0.298
## 30-180 == 0 0.002567 0.050670 0.051 0.960
## 60-120 == 0 -0.020594 0.049224 -0.418 0.751
## 60-180 == 0 0.048913 0.050093 0.976 0.443
## 120-180 == 0 0.069507 0.049844 1.394 0.298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
##
##
## $`50`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 0-30 == 0 -0.34143 0.05087 -6.711 4.82e-11 ***
## 0-60 == 0 -0.45993 0.04968 -9.258 < 2e-16 ***
## 0-120 == 0 -0.48353 0.04945 -9.777 < 2e-16 ***
## 0-180 == 0 -0.58374 0.04855 -12.024 < 2e-16 ***
## 30-60 == 0 -0.11850 0.04506 -2.630 0.01068 *
## 30-120 == 0 -0.14210 0.04481 -3.171 0.00253 **
## 30-180 == 0 -0.24231 0.04381 -5.531 6.36e-08 ***
## 60-120 == 0 -0.02360 0.04345 -0.543 0.58707
## 60-180 == 0 -0.12381 0.04241 -2.919 0.00501 **
## 120-180 == 0 -0.10022 0.04215 -2.378 0.01936 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
##
##
## $`62,5`
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = ngra ~ bloc + umid * K, family = poisson, data = soja)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 0-30 == 0 -0.16130 0.05153 -3.130 0.0025 **
## 0-60 == 0 -0.43608 0.04860 -8.972 < 2e-16 ***
## 0-120 == 0 -0.52759 0.05001 -10.550 < 2e-16 ***
## 0-180 == 0 -0.51111 0.04791 -10.668 < 2e-16 ***
## 30-60 == 0 -0.27478 0.04635 -5.928 5.11e-09 ***
## 30-120 == 0 -0.36629 0.04782 -7.659 3.73e-14 ***
## 30-180 == 0 -0.34981 0.04562 -7.667 3.73e-14 ***
## 60-120 == 0 -0.09152 0.04465 -2.050 0.0505 .
## 60-180 == 0 -0.07504 0.04229 -1.774 0.0844 .
## 120-180 == 0 0.01648 0.04389 0.375 0.7073
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)