##-----------------------------------------------------------------------------
## 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
##-----------------------------------------------------------------------------
## Dados de experimento com soja sob 3 níveis de umidade do solo (massa
## de água/massa de solo) e níveis de potássio fornecidos na adubação. A
## hipótese sob estudo é que presença de potássio dá suporte para a
## cultura superar a ocorrência de períodos sem chuva. Experimento
## conduzido em casa de vegetação com 5 blocos, 2 plantas por
## u.e. (vaso).
## Para acessar o artigo (pdf online).
## browseURL("http://www.scielo.br/pdf/rca/v43n2/a03v43n2.pdf")
da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
header=TRUE, sep="\t", dec=",")
str(da)
## 'data.frame': 75 obs. of 10 variables:
## $ potassio: int 0 30 60 120 180 0 30 60 120 180 ...
## $ agua : num 37.5 37.5 37.5 37.5 37.5 50 50 50 50 50 ...
## $ bloco : Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ rengrao : num 14.6 21.5 24.6 21.9 28.1 ...
## $ pesograo: num 10.7 13.5 15.8 12.8 14.8 ...
## $ kgrao : num 15.1 17.1 19.1 18.1 19.1 ...
## $ pgrao : num 1.18 0.99 0.82 0.85 0.88 1.05 1.08 0.74 1.01 1.01 ...
## $ ts : int 136 159 156 171 190 140 193 200 208 237 ...
## $ nvi : int 22 2 0 2 0 20 6 6 7 10 ...
## $ nv : int 56 62 66 68 82 63 86 94 86 97 ...
xyplot(rengrao~potassio, groups=agua, data=da,
type=c("p","a"))
xyplot(rengrao~potassio|agua, data=da,
type=c("p","a"))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
da$A <- factor(da$agua)
levels(da$A)
## [1] "37.5" "50" "62.5"
da$K <- factor(da$potassio)
levels(da$K)
## [1] "0" "30" "60" "120" "180"
## Remove obs 74 (uma planta que definou).
da <- da[-74,]
m0 <- lm(rengrao~bloco+A*K, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: rengrao
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 4 100 25 4.89 0.0019 **
## A 2 447 223 43.55 4.6e-12 ***
## K 4 2592 648 126.39 < 2e-16 ***
## A:K 8 286 36 6.97 2.6e-06 ***
## Residuals 55 282 5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Polinômio de segundo grau para o efeito de potássio.
m1 <- lm(rengrao~bloco+A*(potassio+I(potassio^2)), data=da)
par(mfrow=c(2,2)); plot(m1); layout(1)
## Testa o abandono dos termos.
anova(m0, m1)
## Analysis of Variance Table
##
## Model 1: rengrao ~ bloco + A * K
## Model 2: rengrao ~ bloco + A * (potassio + I(potassio^2))
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 282
## 2 61 352 -6 -70.5 2.29 0.048 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1)
## Analysis of Variance Table
##
## Response: rengrao
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 4 100 25 4.34 0.0037 **
## A 2 447 223 38.64 1.4e-11 ***
## potassio 1 2000 2000 346.12 < 2e-16 ***
## I(potassio^2) 1 555 555 96.06 3.8e-14 ***
## A:potassio 2 245 123 21.22 1.0e-07 ***
## A:I(potassio^2) 2 8 4 0.65 0.5248
## Residuals 61 352 6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Com b2 comum aos níveis de água.
## m1 <- lm(rengrao~bloco+A*potassio+I(potassio^2), data=da)
##-----------------------------------------------------------------------------
## Predição com bandas.
pred <- expand.grid(bloco="I",
A=levels(da$A),
potassio=seq(0,180,by=3))
aux <- predict(m1, newdata=pred, interval="confidence")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame': 183 obs. of 6 variables:
## $ bloco : Factor w/ 1 level "I": 1 1 1 1 1 1 1 1 1 1 ...
## $ A : Factor w/ 3 levels "37.5","50","62.5": 1 2 3 1 2 3 1 2 3 1 ...
## $ potassio: num 0 0 0 3 3 3 6 6 6 9 ...
## $ fit : num 14.9 16.7 15.8 15.5 17.5 ...
## $ lwr : num 12.7 14.5 13.6 13.4 15.4 ...
## $ upr : num 17.1 19 18.1 17.6 19.6 ...
## Melhorar a legenda.
tx <- paste("Umidade de ", levels(da$A), "%", sep="")
lgd <- simpleKey(columns=1, text=tx,
points=FALSE, lines=TRUE,
corner=c(0.03,0.97))
xyplot(fit~potassio, groups=A, data=pred,
type="l", key=lgd,
xlab=expression("Potássio no solo"~(mg~dm^{-3})),
ylab=expression("Rendimento de grãos"~(g~vaso^{-1})),
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.1, fill=1,
prepanel=prepanel.cbH, panel.groups=panel.cbH, panel=panel.superpose)
## Como calcular os pontos de máximo? São dados por x_max = -0.5*b1/b2.
##-----------------------------------------------------------------------------
## Testar se o rendimento sem potássio é o mesmo para as águas.
X <- LSmatrix(m1, effect="A", at=list("potassio"=0, "I(potassio^2)"=0))
rownames(X) <- levels(da$A)
g <- glht(m1, linfct=X)
confint(g)
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = rengrao ~ bloco + A * (potassio + I(potassio^2)),
## data = da)
##
## Quantile = 2.452
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 37.5 == 0 14.356 12.016 16.696
## 50 == 0 16.194 13.854 18.534
## 62.5 == 0 15.276 12.919 17.632
## Contrastes par a par.
Xc <- apc(X)
summary(glht(m1, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = rengrao ~ bloco + A * (potassio + I(potassio^2)),
## data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 37.5:0:0-50:0:0 == 0 -1.838 1.349 -1.36 0.37
## 37.5:0:0-62.5:0:0 == 0 -0.920 1.354 -0.68 0.78
## 50:0:0-62.5:0:0 == 0 0.918 1.354 0.68 0.78
## (Adjusted p values reported -- single-step method)
## O teste simultâneo.
formula(m1)
## rengrao ~ bloco + A * (potassio + I(potassio^2))
m2 <- update(m1, formula=.~.-A)
coef(m2)
## (Intercept) blocoII blocoIII blocoIV
## 15.8304341 1.0660000 -0.8606667 -1.4433333
## blocoV potassio I(potassio^2) potassio:A50
## -1.5385025 0.1701733 -0.0006452 0.0922962
## potassio:A62.5 I(potassio^2):A50 I(potassio^2):A62.5
## 0.1074427 -0.0004149 -0.0002182
anova(m2, m1)
## Analysis of Variance Table
##
## Model 1: rengrao ~ bloco + potassio + I(potassio^2) + A:potassio + A:I(potassio^2)
## Model 2: rengrao ~ bloco + A * (potassio + I(potassio^2))
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 63 363
## 2 61 352 2 10.7 0.93 0.4
anova(m2)
## Analysis of Variance Table
##
## Response: rengrao
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 4 100 25 4.35 0.0036 **
## potassio 1 1979 1979 343.26 < 2e-16 ***
## I(potassio^2) 1 536 536 92.96 5.1e-14 ***
## potassio:A 2 693 346 60.07 2.5e-15 ***
## I(potassio^2):A 2 36 18 3.14 0.0502 .
## Residuals 63 363 6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Predição com bandas.
pred <- expand.grid(bloco="I",
A=levels(da$A),
potassio=seq(0,180,by=3))
aux <- predict(m2, newdata=pred, interval="confidence")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame': 183 obs. of 6 variables:
## $ bloco : Factor w/ 1 level "I": 1 1 1 1 1 1 1 1 1 1 ...
## $ A : Factor w/ 3 levels "37.5","50","62.5": 1 2 3 1 2 3 1 2 3 1 ...
## $ potassio: num 0 0 0 3 3 3 6 6 6 9 ...
## $ fit : num 15.8 15.8 15.8 16.3 16.6 ...
## $ lwr : num 14.3 14.3 14.3 14.8 15.1 ...
## $ upr : num 17.4 17.4 17.4 17.9 18.1 ...
xyplot(fit~potassio, groups=A, data=pred,
type="l", key=lgd,
xlab=expression("Potássio no solo"~(mg~dm^{-3})),
ylab=expression("Rendimento de grãos"~(g~vaso^{-1})),
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.1, fill=1,
prepanel=prepanel.cbH, panel.groups=panel.cbH, panel=panel.superpose)
##-----------------------------------------------------------------------------
## Coeficientes das equações. Mudar a formula do modelo para ter uma
## interpretação de parâmetros por estrato.
formula(m2)
## rengrao ~ bloco + potassio + I(potassio^2) + A:potassio + A:I(potassio^2)
m3 <- lm(rengrao~A:(potassio+I(potassio^2))+bloco, data=da,
contrasts=list(bloco=contr.sum))
summary(m3)
##
## Call:
## lm(formula = rengrao ~ A:(potassio + I(potassio^2)) + bloco,
## data = da, contrasts = list(bloco = contr.sum))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.080 -1.471 0.076 1.602 4.798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.275134 0.551591 27.69 < 2e-16 ***
## bloco1 0.555301 0.555573 1.00 0.3214
## bloco2 1.621301 0.555573 2.92 0.0049 **
## bloco3 -0.305366 0.555573 -0.55 0.5845
## bloco4 -0.888033 0.555573 -1.60 0.1150
## A37.5:potassio 0.170173 0.022111 7.70 1.2e-10 ***
## A50:potassio 0.262470 0.022111 11.87 < 2e-16 ***
## A62.5:potassio 0.277616 0.022817 12.17 < 2e-16 ***
## A37.5:I(potassio^2) -0.000645 0.000128 -5.02 4.5e-06 ***
## A50:I(potassio^2) -0.001060 0.000128 -8.25 1.3e-11 ***
## A62.5:I(potassio^2) -0.000863 0.000131 -6.58 1.1e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.4 on 63 degrees of freedom
## Multiple R-squared: 0.902, Adjusted R-squared: 0.886
## F-statistic: 58 on 10 and 63 DF, p-value: <2e-16
## O R² é o quadrado da correlação entre observado e predito.
cor(da$rengrao, fitted(m2))^2
## [1] 0.902
## R² separado por água.
ddply(data.frame(o=da$rengrao, f=fitted(m2)), .(da$A),
summarise, R2=cor(o, f)^2)
## da$A R2
## 1 37.5 0.7854
## 2 50 0.8530
## 3 62.5 0.9458
##-----------------------------------------------------------------------------
## Os pontos de máximo.
## Buscas essas ocorrências de texto.
oc <- paste0("^A", levels(da$A), ":")
## Posições de ocorrência.
c3 <- coef(m3)
po <- sapply(oc, grep, x=names(c3))
apply(po, 2,
function(i){
-0.5*c3[i[1]]/c3[i[2]]
})
## ^A37.5: ^A50: ^A62.5:
## 131.9 123.8 160.8
##-----------------------------------------------------------------------------
## E se ajustar um modelo linear-platô? E se for um quadrático-platô?