Experimento em arranjo fatorial com resposta binomial
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## lattice latticeExtra doBy multcomp reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUEsource("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.3 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-7
## [5] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 grid_3.1.1 htmltools_0.2.6
## [6] Matrix_1.1-4 Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1
## [11] yaml_2.1.13 zoo_1.7-10## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)Importação dos dados
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "caiaue.txt"
## url <- "http://www.leg.ufpr.br/~walmes/data/caiaue.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
names(da) <- substr(names(da), 1, 4)
## Convertendo para fator.
da$peri <- factor(da$peri)
## Cria uma variável umidade numérica com os pontos médios.
levels(da$umid)## [1] "18-19" "19-20" "20-21" "21-22" "22-23"da$umi <- as.numeric(as.character(
factor(da$umid, labels=c(18:22+0.5))))
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 45 obs. of 10 variables:
## $ peri: Factor w/ 3 levels "55","75","100": 1 1 1 1 1 1 1 1 1 1 ...
## $ umid: Factor w/ 5 levels "18-19","19-20",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ trat: int 1 1 1 2 2 2 3 3 3 4 ...
## $ rep : int 1 2 3 1 2 3 1 2 3 1 ...
## $ germ: num 34.52 36 30.79 6.03 5.26 ...
## $ pcon: num 11.61 6.86 3.16 23 18.97 ...
## $ ivg : num 5.81 6.04 4.61 14.24 10.61 ...
## $ sg : int 107 126 117 241 194 221 272 344 331 188 ...
## $ sng : int 203 224 263 159 175 167 183 96 90 205 ...
## $ umi : num 18.5 18.5 18.5 19.5 19.5 19.5 20.5 20.5 20.5 21.5 ...## Os dados são de um experimento para verificar o efeito a adubação com
## potássio nos componentes de produção de soja sob diferentes níveis de
## umidade.
##-----------------------------------------------------------------------------
## Análise exploratória.
## Proporção amostral.
da$pa <- with(da, sg/(sg+sng))
xyplot(pa~umi, groups=peri, data=da,
type=c("p","a"), auto.key=TRUE, jitter.x=TRUE)
Especificação e estimação do modelo
\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Binomial}(\pi_{ij}, n_{ijk}) \newline g(\pi_{ij}) &= \mu+\alpha_i+\tau_j \end{aligned} \]
em que \(\alpha_i\) é o efeito período \(i\), \(\tau_j\) o efeito da umidade \(j\) sobre uma função \(g\) do parâmetro \(\pi\) que representa a probabilidade de germinação. A variável resposta \(Y\) representa o número de sementes germinadas e não germinadas cujo total é \(n_{ijk}\). \(\mu\) é a estimativa na ausência do efeito dos termos mencioandos. \(\pi_{ij}\) é a estimativa para as observações na combinação \(ij\). A função \(g\) é a chamada função de ligação. Para distribuição binomial é comum usar a função de ligação logit, a probit e a complementar log-log.
##-----------------------------------------------------------------------------
## Espeficicação do modelo, considerando distribuição normal para a
## proporção amostral.
m0 <- lm(pa~peri*umid, data=da)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Apresenta relação mádia variância condizente com a suposição de
## distribuição binomial onde V(Y) ~= p*(1-p).
##-----------------------------------------------------------------------------
## Análise considerando distribuição binomial (ou quasibinomial).
## Quasibinomial estima parâmetro de dispersão. Link default é logit.
m1 <- glm(cbind(sg, sng)~peri*umid, data=da,
family=quasibinomial(link="logit"))
summary(m1)##
## Call:
## glm(formula = cbind(sg, sng) ~ peri * umid, family = quasibinomial(link = "logit"),
## data = da)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.3424 -2.0733 -0.3986 2.3595 6.8903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6788 0.2476 -2.741 0.010219 *
## peri75 -0.4342 0.3600 -1.206 0.237236
## peri100 -4.0652 1.2871 -3.158 0.003605 **
## umid19-20 0.9483 0.3338 2.841 0.008013 **
## umid20-21 1.6213 0.3390 4.782 4.31e-05 ***
## umid21-22 0.9524 0.3365 2.830 0.008219 **
## umid22-23 0.3132 0.3319 0.944 0.352824
## peri75:umid19-20 -0.8575 0.4926 -1.741 0.091931 .
## peri100:umid19-20 1.5202 1.3591 1.119 0.272197
## peri75:umid20-21 0.6679 0.4995 1.337 0.191280
## peri100:umid20-21 3.4288 1.3267 2.585 0.014862 *
## peri75:umid21-22 1.0923 0.4933 2.214 0.034549 *
## peri100:umid21-22 4.7038 1.3290 3.539 0.001330 **
## peri75:umid22-23 1.7087 0.4841 3.529 0.001366 **
## peri100:umid22-23 5.7034 1.3337 4.276 0.000178 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 14.23902)
##
## Null deviance: 5162.40 on 44 degrees of freedom
## Residual deviance: 433.37 on 30 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5summary(m1)$dispersion## [1] 14.23902## Diagnóstico.
par(mfrow=c(2,2)); plot(m1); layout(1)
## Relação média variância incorporada no modelo.
## Quadro de análise de deviance.
anova(m1, test="F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(sg, sng)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 44 5162.4
## peri 2 143.34 42 5019.1 5.0332 0.01304 *
## umid 4 3008.32 38 2010.7 52.8182 3.711e-13 ***
## peri:umid 8 1577.38 30 433.4 13.8473 3.552e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Comparações múltiplas.
LSmeans(m1, effect=c("peri","umid"))## estimate se z.stat p.value lwr upr peri umid
## 1 -0.6787584 0.2476271 -2.741050 6.124311e-03 -1.1640987 -0.1934182 55 18-19
## 2 -1.1129495 0.2613337 -4.258729 2.055928e-05 -1.6251541 -0.6007448 75 18-19
## 3 -4.7439654 1.2630939 -3.755830 1.727681e-04 -7.2195839 -2.2683470 100 18-19
## 4 0.2695547 0.2238907 1.203957 2.286063e-01 -0.1692630 0.7083724 55 19-20
## 5 -1.0221706 0.2507289 -4.076796 4.566046e-05 -1.5135902 -0.5307510 75 19-20
## 6 -2.2754341 0.3744295 -6.077069 1.223991e-09 -3.0093025 -1.5415657 100 19-20
## 7 0.9425024 0.2315687 4.070077 4.699766e-05 0.4886361 1.3963688 55 20-21
## 8 1.1761882 0.2574933 4.567840 4.927748e-06 0.6715107 1.6808658 75 20-21
## 9 0.3061227 0.2230354 1.372529 1.698987e-01 -0.1310187 0.7432641 100 20-21
## 10 0.2736083 0.2278246 1.200961 2.297665e-01 -0.1729196 0.7201363 55 21-22
## 11 0.9317470 0.2486554 3.747141 1.788616e-04 0.4443913 1.4191027 75 21-22
## 12 0.9122480 0.2398229 3.803841 1.424698e-04 0.4422038 1.3822922 100 21-22
## 13 -0.3655424 0.2209583 -1.654350 9.805638e-02 -0.7986127 0.0675279 55 22-23
## 14 0.9089812 0.2365510 3.842644 1.217159e-04 0.4453498 1.3726125 75 22-23
## 15 1.2726443 0.2705582 4.703773 2.553972e-06 0.7423600 1.8029287 100 22-23## Os erros padrões são diferentes mesmo o experimento sendo balanceado
## porque na binomial o existe relação média variãncia.
## Desdobrar umidade dentro de período.
X <- LSmatrix(m1, effect=c("peri","umid"))
grid <- equallevels(attr(X, "grid"), da)
## Divide a matriz de acordo com os níveis de peri.
L <- by(X, INDICES=grid$peri, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(da$umid))
grid <- lapply(L, apmc, model=m0, focus="umid", test="fdr")
grid## $`55`
## umid estimate lwr upr cld
## 1 18-19 0.3376853 0.2373469 0.4380238 c
## 2 19-20 0.5659443 0.4656058 0.6662828 ab
## 3 20-21 0.7219479 0.6216094 0.8222863 a
## 4 21-22 0.5702004 0.4698620 0.6705389 ab
## 5 22-23 0.4250422 0.3247037 0.5253807 bc
##
## $`75`
## umid estimate lwr upr cld
## 1 18-19 0.2448451 0.1445067 0.3451836 b
## 2 19-20 0.2590378 0.1586993 0.3593762 b
## 3 20-21 0.7643027 0.6639643 0.8646412 a
## 4 21-22 0.7207804 0.6204419 0.8211189 a
## 5 22-23 0.7138861 0.6135477 0.8142246 a
##
## $`100`
## umid estimate lwr upr cld
## 1 18-19 0.009118541 -0.091219914 0.1094570 c
## 2 19-20 0.092481969 -0.007856486 0.1928204 c
## 3 20-21 0.578140872 0.477802417 0.6784793 a
## 4 21-22 0.713216054 0.612877599 0.8135545 ab
## 5 22-23 0.778372671 0.678034216 0.8787111 bgrid <- ldply(grid); names(grid)[1] <- "peri"
grid <- equallevels(grid, da)
segplot(umid~lwr+upr, data=grid, centers=estimate,
groups=grid$peri, pch=grid$peri, horizontal=FALSE,
f=0.1, panel=panel.segplotBy,
ylab="Log da razão de chances",
key=list(text=list(levels(da$peri)),
points=list(pch=1:3)))
## Passando para escala da probabilidade de germinação.
ci <- sapply(grid[,c("estimate","lwr","upr")], m1$family$linkinv)
colnames(ci) <- paste0("p", colnames(ci))
grid <- cbind(grid, ci)
segplot(umid~plwr+pupr, data=grid, centers=pestimate,
groups=grid$peri, pch=grid$peri, horizontal=FALSE,
f=0.1, panel=panel.segplotBy,
ylab="Probabilidade de germinação",
xlab="Faixa de umidade no armazenamento (%)",
key=list(text=list(levels(da$peri)),
points=list(pch=1:3)))
## Análise do IVG##-----------------------------------------------------------------------------
## Ver o gráfico.
xyplot(ivg~umid, groups=peri, data=da,
type=c("p","a"), auto.key=TRUE)
## Soma 0.01 para ter valores de y>0.
m0 <- lm(ivg+0.01~peri*umid, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Relação média variância.
MASS::boxcox(m0)
abline(v=0.5, col=2)
## Transformando com raíz quadrada.
m0 <- lm(sqrt(ivg+0.01)~peri*umid, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de anova.
anova(m0)## Analysis of Variance Table
##
## Response: sqrt(ivg + 0.01)
## Df Sum Sq Mean Sq F value Pr(>F)
## peri 2 7.356 3.6778 26.453 2.388e-07 ***
## umid 4 68.573 17.1432 123.304 < 2.2e-16 ***
## peri:umid 8 16.830 2.1038 15.132 1.298e-08 ***
## Residuals 30 4.171 0.1390
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Desdobrar umidade dentro de período.
X <- LSmatrix(m1, effect=c("peri","umid"))
grid <- equallevels(attr(X, "grid"), da)
## Divide a matriz de acordo com os níveis de peri.
L <- by(X, INDICES=grid$peri, FUN=as.matrix)
L <- lapply(L, "rownames<-", levels(da$umid))
grid <- lapply(L, apmc, model=m0, focus="umid", test="fdr")
grid## $`55`
## umid estimate lwr upr cld
## 1 18-19 2.340520 1.900866 2.780174 c
## 2 19-20 3.516003 3.076349 3.955657 b
## 3 20-21 4.823600 4.383946 5.263254 a
## 4 21-22 3.858843 3.419189 4.298497 b
## 5 22-23 3.345226 2.905573 3.784880 b
##
## $`75`
## umid estimate lwr upr cld
## 1 18-19 1.878450 1.438796 2.318104 b
## 2 19-20 2.060367 1.620713 2.500021 b
## 3 20-21 5.056444 4.616790 5.496098 a
## 4 21-22 4.459289 4.019635 4.898943 a
## 5 22-23 5.024509 4.584855 5.464163 a
##
## $`100`
## umid estimate lwr upr cld
## 1 18-19 0.2552285 -0.1844254 0.6948824 d
## 2 19-20 1.2093229 0.7696689 1.6489768 c
## 3 20-21 3.6264047 3.1867508 4.0660586 a
## 4 21-22 4.1382976 3.6986437 4.5779515 ab
## 5 22-23 4.6951948 4.2555409 5.1348488 bgrid <- ldply(grid); names(grid)[1] <- "peri"
grid <- equallevels(grid, da)
segplot(umid~lwr+upr, data=grid, centers=estimate,
groups=grid$peri, pch=grid$peri, horizontal=FALSE,
f=0.1, panel=panel.segplotBy,
ylab=expression(IVG~(y[t]==(y+0.01)^{1/2})),
key=list(text=list(levels(da$peri)),
points=list(pch=1:3)))
##-----------------------------------------------------------------------------
## Análise com pacote ExpDes.
require(ExpDes)## Loading required package: ExpDes
##
## Attaching package: 'ExpDes'
##
## The following object is masked from 'package:MASS':
##
## ginvwith(da, fat2.crd(factor1=peri, factor2=umid,
resp=sqrt(ivg+0.01), quali=c(TRUE,TRUE),
mcomp="tukey"))## ------------------------------------------------------------------------
## Legend:
## FACTOR 1: F1
## FACTOR 2: F2
## ------------------------------------------------------------------------
##
##
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## F1 2 7.356 3.6778 26.453 2.3882e-07
## F2 4 68.573 17.1432 123.304 0.0000e+00
## F1*F2 8 16.830 2.1038 15.132 1.2985e-08
## Residuals 30 4.171 0.1390
## Total 44 96.929
## ------------------------------------------------------------------------
## CV = 11.12 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.7373667
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
##
##
## Significant interaction: analyzing the interaction
## ------------------------------------------------------------------------
##
## Analyzing F1 inside of each level of F2
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr.Fc
## F2 4 68.57272 17.14318 123.3036 0
## F2:F1 18-19 2 7.19680 3.59840 25.8818 0
## F2:F1 19-20 2 8.16392 4.08196 29.3598 0
## F2:F1 20-21 2 3.53251 1.76625 12.7039 1e-04
## F2:F1 21-22 2 0.54167 0.27083 1.948 0.1602
## F2:F1 22-23 2 4.75085 2.37543 17.0854 0
## Residuals 30 4.17097 0.13903
## Total 44 96.92944 2.20294
## ------------------------------------------------------------------------
##
##
##
## F1 inside of the level 18-19 of F2
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 1 2.34052
## a 2 1.87845
## b 3 0.2552285
## ------------------------------------------------------------------------
##
##
## F1 inside of the level 19-20 of F2
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 1 3.516003
## b 2 2.060367
## c 3 1.209323
## ------------------------------------------------------------------------
##
##
## F1 inside of the level 20-21 of F2
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 5.056444
## a 1 4.8236
## b 3 3.626405
## ------------------------------------------------------------------------
##
##
## F1 inside of the level 21-22 of F2
##
## According to the F test, the means of this factor are statistical equal.
## Levels Means
## 1 1 3.858843
## 2 2 4.459289
## 3 3 4.138298
## ------------------------------------------------------------------------
##
##
## F1 inside of the level 22-23 of F2
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 2 5.024509
## a 3 4.695195
## b 1 3.345226
## ------------------------------------------------------------------------
##
##
##
## Analyzing F2 inside of each level of F1
## ------------------------------------------------------------------------
## ------------------------------------------------------------------------
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr.Fc
## F1 2 7.35558 3.67779 26.4528 0
## F1:F2 55 4 9.65931 2.41483 17.3688 0
## F1:F2 75 4 30.53141 7.63285 54.8998 0
## F1:F2 100 4 45.21217 11.30304 81.2979 0
## Residuals 30 4.17097 0.13903
## Total 44 96.92944 2.20294
## ------------------------------------------------------------------------
##
##
##
## F2 inside of the level 55 of F1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 4.8236
## b 4 3.858843
## b 2 3.516003
## b 5 3.345226
## c 1 2.34052
## ------------------------------------------------------------------------
##
##
## F2 inside of the level 75 of F1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 3 5.056444
## a 5 5.024509
## a 4 4.459289
## b 2 2.060367
## b 1 1.87845
## ------------------------------------------------------------------------
##
##
## F2 inside of the level 100 of F1
## ------------------------------------------------------------------------
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 5 4.695195
## ab 4 4.138298
## b 3 3.626405
## c 2 1.209323
## d 1 0.2552285
## ------------------------------------------------------------------------