##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "reshape",
"doBy", "multcomp",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)
## lattice latticeExtra reshape doBy multcomp plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUE
source("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] splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.1 plyr_1.8.1 multcomp_1.3-3 TH.data_1.0-3
## [5] mvtnorm_0.9-99992 doBy_4.5-10 MASS_7.3-33 survival_2.37-7
## [9] reshape_0.8.5 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] knitr_1.5
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 grid_3.1.1 lme4_1.1-6
## [5] Matrix_1.1-4 methods_3.1.1 minqa_1.2.3 nlme_3.1-117
## [9] Rcpp_0.11.1 RcppEigen_0.3.2.1.2 sandwich_2.3-0 stringr_0.6.2
## [13] tools_3.1.1 zoo_1.7-11
## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)
Análise dos dados
##-----------------------------------------------------------------------------
## Experimento com doses de cama de frango como adubação em milho
## safrinha com 3 métodos de menejo.
da <-
read.table("http://www.leg.ufpr.br/~walmes/data/camafrango_prod.txt",
header=TRUE, sep="\t")
str(da)
## 'data.frame': 60 obs. of 13 variables:
## $ manejo: Factor w/ 3 levels "incorp","incorp120",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ cama : int 0 0 0 0 2 2 2 2 4 4 ...
## $ bloco : int 1 2 3 4 1 2 3 4 1 2 ...
## $ altes : num 38.4 69.8 52.8 39.2 60 54.6 55.2 50.8 62.6 54 ...
## $ altpl : num 109.2 131.4 103.8 86.8 126 ...
## $ maspl : num 53.6 72.6 83.6 55.6 97.6 60.6 72.6 81.6 88.6 82.6 ...
## $ stand : num 9375 31250 16406 12500 17188 ...
## $ mases : num 95.4 111.5 95.7 67.2 96.8 ...
## $ masgr : num 60.8 73.4 58.6 35 60.9 ...
## $ indice: num 29 28.5 24.6 22.2 23.9 ...
## $ prod : num 570 2293 961 438 1047 ...
## $ p100 : num 25 24.2 28.6 25.2 26.1 ...
## $ dens : num 0.718 0.741 0.703 0.716 0.689 0.73 0.682 0.73 0.739 0.712 ...
da$bloco <- factor(da$bloco)
## Tabelas de frequência.
with(da, ftable(cama~bloco+manejo))
## cama 0 2 4 8 16
## bloco manejo
## 1 incorp 1 1 1 1 1
## incorp120 1 1 1 1 1
## lanco120 1 1 1 1 1
## 2 incorp 1 1 1 1 1
## incorp120 1 1 1 1 1
## lanco120 1 1 1 1 1
## 3 incorp 1 1 1 1 1
## incorp120 1 1 1 1 1
## lanco120 1 1 1 1 1
## 4 incorp 1 1 1 1 1
## incorp120 1 1 1 1 1
## lanco120 1 1 1 1 1
##-----------------------------------------------------------------------------
## Ver.
xyplot(prod~cama, groups=manejo, data=da,
type=c("p","a"), auto.key=TRUE)
xyplot(prod~cama|manejo, groups=bloco, data=da,
type=c("p","a"), auto.key=TRUE)
xyplot(stand~cama, groups=manejo, data=da,
type=c("p","a"), auto.key=TRUE)
xyplot(prod/stand~cama, groups=manejo, data=da,
type=c("p","a"), auto.key=TRUE)
## Se fossem tratamentos aplicados após estabelecimento de stand...
##-----------------------------------------------------------------------------
## Ajuste do modelo.
## Modelo considerando cama como fator (modelo saturado).
m0 <- lm(prod~bloco+manejo*factor(cama), data=da)
par(mfrow=c(2,2))
plot(m0)
layout(1)
## Uma transformação Box-Cox confere melhoras?
MASS::boxcox(m0)
anova(m0)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 5927481 1975827 3.81 0.017 *
## manejo 2 1316030 658015 1.27 0.292
## factor(cama) 4 7034127 1758532 3.39 0.017 *
## manejo:factor(cama) 8 6755740 844468 1.63 0.146
## Residuals 42 21788948 518784
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Modelo considerando efeito linear para cama.
m1 <- lm(prod~bloco+manejo*cama, data=da)
par(mfrow=c(2,2))
plot(m1)
layout(1)
## Teste da falta de ajuste.
anova(m0, m1)
## Analysis of Variance Table
##
## Model 1: prod ~ bloco + manejo * factor(cama)
## Model 2: prod ~ bloco + manejo * cama
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 42 21788948
## 2 51 23291258 -9 -1502311 0.32 0.96
anova(m1)
## Analysis of Variance Table
##
## Response: prod
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 5927481 1975827 4.33 0.00859 **
## manejo 2 1316030 658015 1.44 0.24620
## cama 1 6335428 6335428 13.87 0.00049 ***
## manejo:cama 2 5952129 2976064 6.52 0.00302 **
## Residuals 51 23291258 456691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Partição das somas de quadrados dos efeitos lineares.
m2 <- aov(prod~bloco+manejo/cama, data=da)
coef(m2)
## (Intercept) bloco2 bloco3 bloco4
## 1341.71 29.69 -427.86 -726.09
## manejoincorp120 manejolanco120 manejoincorp:cama manejoincorp120:cama
## 916.74 -140.53 76.57 -18.27
## manejolanco120:cama
## 114.03
## Ocorrências.
oc <- paste0("manejo", levels(da$manejo), ":cama")
as <- m2$assign; as
## [1] 0 1 1 1 2 2 3 3 3
linear.ef <- lapply(oc, grep, x=names(coef(m2))[as==3])
names(linear.ef) <- abbreviate(levels(da$manejo))
linear.ef
## $incr
## [1] 1
##
## $i120
## [1] 2
##
## $l120
## [1] 3
## Somas de quadrados repartidas por manejo testando o efeito linear de
## cama de frango.
summary(m2, split=list("manejo:cama"=linear.ef))
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 5927481 1975827 4.33 0.0086 **
## manejo 2 1316030 658015 1.44 0.2462
## manejo:cama 3 12287556 4095852 8.97 7.1e-05 ***
## manejo:cama: incr 1 3752360 3752360 8.22 0.0060 **
## manejo:cama: i120 1 213624 213624 0.47 0.4971
## manejo:cama: l120 1 8321572 8321572 18.22 8.6e-05 ***
## Residuals 51 23291258 456691
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Coeficientes.
## Estão sob restrição.
summary(m1)
##
## Call:
## lm(formula = prod ~ bloco + manejo * cama, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1914.7 -420.4 36.7 393.2 1108.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1341.7 267.1 5.02 6.6e-06 ***
## bloco2 29.7 246.8 0.12 0.9047
## bloco3 -427.9 246.8 -1.73 0.0890 .
## bloco4 -726.1 246.8 -2.94 0.0049 **
## manejoincorp120 916.7 311.5 2.94 0.0049 **
## manejolanco120 -140.5 311.5 -0.45 0.6538
## cama 76.6 26.7 2.87 0.0060 **
## manejoincorp120:cama -94.8 37.8 -2.51 0.0153 *
## manejolanco120:cama 37.5 37.8 0.99 0.3261
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 676 on 51 degrees of freedom
## Multiple R-squared: 0.456, Adjusted R-squared: 0.371
## F-statistic: 5.35 on 8 and 51 DF, p-value: 6.79e-05
M <- model.matrix(m1)
attributes(M)
## $dim
## [1] 60 9
##
## $dimnames
## $dimnames[[1]]
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17"
## [18] "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30" "31" "32" "33" "34"
## [35] "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47" "48" "49" "50" "51"
## [52] "52" "53" "54" "55" "56" "57" "58" "59" "60"
##
## $dimnames[[2]]
## [1] "(Intercept)" "bloco2" "bloco3"
## [4] "bloco4" "manejoincorp120" "manejolanco120"
## [7] "cama" "manejoincorp120:cama" "manejolanco120:cama"
##
##
## $assign
## [1] 0 1 1 1 2 2 3 4 4
##
## $contrasts
## $contrasts$bloco
## [1] "contr.treatment"
##
## $contrasts$manejo
## [1] "contr.treatment"
## Estimativas com 0 e 2 ton de cama.
LSmeans(m1, effect="manejo", at=list(cama=0))
## estimate se df t.stat p.value lwr upr manejo cama
## 1 1060.6 220.3 51 4.815 1.356e-05 618.4 1503 incorp 0
## 2 1977.4 220.3 51 8.977 4.479e-12 1535.1 2420 incorp120 0
## 3 920.1 220.3 51 4.177 1.156e-04 477.9 1362 lanco120 0
LSmeans(m1, effect="manejo", at=list(cama=2))
## estimate se df t.stat p.value lwr upr manejo cama
## 1 1214 185.1 51 6.558 2.712e-08 842.2 1585 incorp 2
## 2 1941 185.1 51 10.487 2.487e-14 1569.3 2312 incorp120 2
## 3 1148 185.1 51 6.204 9.817e-08 776.6 1520 lanco120 2
## E as estimativas dos coeficientes de inclinação? Truque.
X1 <- LSmatrix(m1, effect="manejo", at=list(cama=1))
X0 <- LSmatrix(m1, effect="manejo", at=list(cama=0))
X <- X1-X0
rownames(X) <- levels(da$manejo)
g <- glht(m1, linfct=X)
summary(g, test=adjusted(type="fdr"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = prod ~ bloco + manejo * cama, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## incorp == 0 76.6 26.7 2.87 0.00903 **
## incorp120 == 0 -18.3 26.7 -0.68 0.49711
## lanco120 == 0 114.0 26.7 4.27 0.00026 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
confint(g, calpha=univariate_calpha())
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = prod ~ bloco + manejo * cama, data = da)
##
## Quantile = 2.008
## 95% confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## incorp == 0 76.571 22.942 130.199
## incorp120 == 0 -18.270 -71.898 35.359
## lanco120 == 0 114.028 60.400 167.657
##-----------------------------------------------------------------------------
## A mesma coisa de uma forma mais simples mudando a fórmula do modelo.
m2 <- lm(prod~0+manejo/cama+bloco, data=da,
contrasts=list(bloco=contr.sum))
summary(m2)
##
## Call:
## lm(formula = prod ~ 0 + manejo/cama + bloco, data = da, contrasts = list(bloco = contr.sum))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1914.7 -420.4 36.7 393.2 1108.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## manejoincorp 1060.6 220.3 4.81 1.4e-05 ***
## manejoincorp120 1977.4 220.3 8.98 4.5e-12 ***
## manejolanco120 920.1 220.3 4.18 0.00012 ***
## bloco1 281.1 151.1 1.86 0.06866 .
## bloco2 310.8 151.1 2.06 0.04487 *
## bloco3 -146.8 151.1 -0.97 0.33591
## manejoincorp:cama 76.6 26.7 2.87 0.00602 **
## manejoincorp120:cama -18.3 26.7 -0.68 0.49711
## manejolanco120:cama 114.0 26.7 4.27 8.6e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 676 on 51 degrees of freedom
## Multiple R-squared: 0.889, Adjusted R-squared: 0.869
## F-statistic: 45.2 on 9 and 51 DF, p-value: <2e-16
##-----------------------------------------------------------------------------
## Testando se as inclinações são iguais.
Xc <- apc(X)
h <- glht(m1, linfct=Xc)
summary(h, test=adjusted(type="fdr"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = prod ~ bloco + manejo * cama, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## incorp:1-incorp120:1 == 0 94.8 37.8 2.51 0.0229 *
## incorp:1-lanco120:1 == 0 -37.5 37.8 -0.99 0.3261
## incorp120:1-lanco120:1 == 0 -132.3 37.8 -3.50 0.0029 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)
##-----------------------------------------------------------------------------
## Gráfico dos observados, preditos e bandas.
## São os valores preditos sob **efeito do bloco 1!**
pred <- expand.grid(bloco="1",
manejo=levels(da$manejo),
cama=seq(0,16,by=2))
aux <- predict(m1, newdata=pred, interval="confidence")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame': 27 obs. of 6 variables:
## $ bloco : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## $ manejo: Factor w/ 3 levels "incorp","incorp120",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ cama : num 0 0 0 2 2 2 4 4 4 6 ...
## $ fit : num 1342 2258 1201 1495 2222 ...
## $ lwr : num 805 1722 665 1015 1742 ...
## $ upr : num 1878 2795 1737 1975 2702 ...
## Melhorar a legenda.
tx <- c("Incorporado", "Incorporado+NPK", "À lanço+NPK")
lgd <- simpleKey(columns=1, text=tx,
points=FALSE, lines=TRUE,
corner=c(0.03,0.97))
xyplot(fit~cama, groups=manejo, data=pred,
type="l", key=lgd,
xlab=expression("Cama de frango"~(ton~ha^{-1})),
ylab=expression("Produtividade"~(kg~ha^{-1})),
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.5,# fill=1,
prepanel=prepanel.cbH, panel.groups=panel.cbH, panel=panel.superpose)
##-----------------------------------------------------------------------------
## Como especificar um modelo com as restrições de que 1) b1_incorp120==0
## e 2) b1_incorp==b1_lanco120?
## Em 1).
da$cama.res <- ifelse(da$manejo=="incorp120", 0, 1)
## Em 2).
da$man.res <- da$manejo
levels(da$man.res)
## [1] "incorp" "incorp120" "lanco120"
levels(da$man.res) <- c("inclan", "inc120", "inclan")
levels(da$man.res)
## [1] "inclan" "inc120"
with(da, ftable(manejo~man.res))
## manejo incorp incorp120 lanco120
## man.res
## inclan 20 0 20
## inc120 0 20 0
m3 <- lm(prod~bloco+manejo+man.res:I(cama*cama.res), data=da)
coef(m3)
## (Intercept) bloco2
## 1229.34 29.69
## bloco3 bloco4
## -427.86 -726.09
## manejoincorp120 manejolanco120
## 919.49 84.22
## man.resinclan:I(cama * cama.res) man.resinc120:I(cama * cama.res)
## 95.30 NA
anova(m1, m3)
## Analysis of Variance Table
##
## Model 1: prod ~ bloco + manejo * cama
## Model 2: prod ~ bloco + manejo + man.res:I(cama * cama.res)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 51 23291258
## 2 53 23953867 -2 -662609 0.73 0.49
##-----------------------------------------------------------------------------
## Predição.
pred$cama.res <- ifelse(pred$manejo=="incorp120", 0, 1)
pred$man.res <- pred$manejo
levels(pred$man.res)
## [1] "incorp" "incorp120" "lanco120"
levels(pred$man.res) <- c("inclan", "inc120", "inclan")
aux <- predict(m3, newdata=pred, interval="confidence")
## Warning: prediction from a rank-deficient fit may be misleading
colnames(aux) <- paste0(colnames(aux), "r")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame': 27 obs. of 11 variables:
## $ bloco : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## $ manejo : Factor w/ 3 levels "incorp","incorp120",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ cama : num 0 0 0 2 2 2 4 4 4 6 ...
## $ fit : num 1342 2258 1201 1495 2222 ...
## $ lwr : num 805 1722 665 1015 1742 ...
## $ upr : num 1878 2795 1737 1975 2702 ...
## $ cama.res: num 1 0 1 1 0 1 1 0 1 1 ...
## $ man.res : Factor w/ 2 levels "inclan","inc120": 1 2 1 1 2 1 1 2 1 1 ...
## $ fitr : num 1229 2149 1314 1420 2149 ...
## $ lwrr : num 747 1722 831 968 1722 ...
## $ uprr : num 1712 2575 1796 1872 2575 ...
xyplot(fitr~cama, groups=manejo, data=pred,
type="l", key=lgd,
xlab=expression("Cama de frango"~(ton~ha^{-1})),
ylab=expression("Produtividade"~(kg~ha^{-1})),
ly=pred$lwrr, uy=pred$uprr, cty="bands", alpha=0.1, fill=1,
prepanel=prepanel.cbH, panel.groups=panel.cbH, panel=panel.superpose)
