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=",")
names(da) <- substr(names(da), 1, 4)
str(da)
## 'data.frame': 75 obs. of 10 variables:
## $ pota: 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 ...
## $ bloc: Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ reng: num 14.6 21.5 24.6 21.9 28.1 ...
## $ peso: num 10.7 13.5 15.8 12.8 14.8 ...
## $ kgra: num 15.1 17.1 19.1 18.1 19.1 ...
## $ pgra: 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(reng~pota, groups=agua, data=da,
type=c("p","a"))
xyplot(reng~pota|agua, data=da,
type=c("p","a"))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
da$A <- factor(da$agua)
da$K <- factor(da$pota)
## Remove obs 74 (uma planta que definou).
da <- da[-74,]
m0 <- lm(reng~bloc+A*K, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: reng
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 100.39 25.10 4.8948 0.001896 **
## A 2 446.53 223.26 43.5456 4.617e-12 ***
## K 4 2592.13 648.03 126.3925 < 2.2e-16 ***
## A:K 8 286.00 35.75 6.9726 2.629e-06 ***
## Residuals 55 281.99 5.13
## ---
## 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(reng~bloc+A*(pota+I(pota^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: reng ~ bloc + A * K
## Model 2: reng ~ bloc + A * (pota + I(pota^2))
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 281.99
## 2 61 352.46 -6 -70.47 2.2907 0.04804 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1)
## Analysis of Variance Table
##
## Response: reng
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 100.39 25.10 4.3434 0.003718 **
## A 2 446.53 223.26 38.6399 1.443e-11 ***
## pota 1 1999.89 1999.89 346.1175 < 2.2e-16 ***
## I(pota^2) 1 555.07 555.07 96.0640 3.821e-14 ***
## A:pota 2 245.17 122.58 21.2152 1.013e-07 ***
## A:I(pota^2) 2 7.53 3.77 0.6516 0.524783
## Residuals 61 352.46 5.78
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Com b2 comum aos níveis de água.
m2 <- lm(reng~bloc+A*pota+I(pota^2), data=da)
anova(m2, m0)
## Analysis of Variance Table
##
## Model 1: reng ~ bloc + A * pota + I(pota^2)
## Model 2: reng ~ bloc + A * K
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 63 359.99
## 2 55 281.99 8 78 1.9017 0.07823 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Predição com bandas.
pred <- expand.grid(bloc="I",
A=levels(da$A),
pota=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:
## $ bloc: 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 ...
## $ pota: num 0 0 0 3 3 3 6 6 6 9 ...
## $ fit : num 14.4 17.2 15.9 15.1 17.9 ...
## $ lwr : num 12.4 15.2 13.9 13.1 15.9 ...
## $ upr : num 16.4 19.2 17.8 17 19.8 ...
## 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~pota, 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(m2, effect="A",
at=list("pota"=0, "I(pota^2)"=0))
rownames(X) <- levels(da$A)
g <- glht(m2, linfct=X)
confint(g)
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = reng ~ bloc + A * pota + I(pota^2), data = da)
##
## Quantile = 2.446
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 37.5 == 0 13.8704 11.8586 15.8823
## 50 == 0 16.6523 14.6404 18.6641
## 62.5 == 0 15.3039 13.2943 17.3136
## Contrastes par a par.
Xc <- apc(X)
summary(glht(m2, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = reng ~ bloc + A * pota + I(pota^2), data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 37.5:0:0-50:0:0 == 0 -2.782 1.060 -2.625 0.0289 *
## 37.5:0:0-62.5:0:0 == 0 -1.433 1.060 -1.352 0.3720
## 50:0:0-62.5:0:0 == 0 1.348 1.060 1.272 0.4162
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## O teste simultâneo, intercepto comum.
formula(m1)
## reng ~ bloc + A * (pota + I(pota^2))
m3 <- update(m2, formula=.~.-A)
coef(m3)
## (Intercept) blocII blocIII blocIV blocV pota
## 15.8313206978 1.0660000000 -0.8606666667 -1.4433333333 -1.5406099373 0.2022515736
## I(pota^2) pota:A50 pota:A62.5
## -0.0008561129 0.0291812865 0.0742331466
anova(m3, m2)
## Analysis of Variance Table
##
## Model 1: reng ~ bloc + pota + I(pota^2) + A:pota
## Model 2: reng ~ bloc + A * pota + I(pota^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 65 399.38
## 2 63 359.99 2 39.384 3.4462 0.03799 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m2)
## Analysis of Variance Table
##
## Response: reng
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 100.39 25.10 4.392 0.003397 **
## A 2 446.53 223.26 39.072 9.223e-12 ***
## pota 1 1999.89 1999.89 349.988 < 2.2e-16 ***
## I(pota^2) 1 555.07 555.07 97.138 2.200e-14 ***
## A:pota 2 245.17 122.58 21.453 7.841e-08 ***
## Residuals 63 359.99 5.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Coeficientes das equações. Mudar a formula do modelo para ter uma
## interpretação de parâmetros por estrato.
## Usar contraste tipo soma para blocos.
formula(m2)
## reng ~ bloc + A * pota + I(pota^2)
m3 <- lm(reng~0+A/pota+I(pota^2)+bloc, data=da,
contrasts=list(bloc=contr.sum))
summary(m3)
##
## Call:
## lm(formula = reng ~ 0 + A/pota + I(pota^2) + bloc, data = da,
## contrasts = list(bloc = contr.sum))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7702 -1.4154 -0.0617 1.5246 5.1096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## A37.5 1.387e+01 8.225e-01 16.864 < 2e-16 ***
## A50 1.665e+01 8.225e-01 20.246 < 2e-16 ***
## A62.5 1.530e+01 8.216e-01 18.627 < 2e-16 ***
## I(pota^2) -8.561e-04 8.602e-05 -9.952 1.51e-14 ***
## bloc1 5.557e-01 5.531e-01 1.005 0.31890
## bloc2 1.622e+00 5.531e-01 2.932 0.00469 **
## bloc3 -3.050e-01 5.531e-01 -0.551 0.58331
## bloc4 -8.876e-01 5.531e-01 -1.605 0.11353
## A37.5:pota 2.129e-01 1.732e-02 12.293 < 2e-16 ***
## A50:pota 2.210e-01 1.732e-02 12.757 < 2e-16 ***
## A62.5:pota 2.763e-01 1.748e-02 15.810 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.39 on 63 degrees of freedom
## Multiple R-squared: 0.9927, Adjusted R-squared: 0.9914
## F-statistic: 778.4 on 11 and 63 DF, p-value: < 2.2e-16
## Prova de que são o mesmo modelo.
deviance(m3); deviance(m2)
## [1] 359.9933
## [1] 359.9933
## O R² é o quadrado da correlação entre observado e predito.
cor(da$reng, fitted(m2))^2
## [1] 0.9028891
## R² separado por água.
ddply(data.frame(o=da$reng, f=fitted(m2)), .(da$A),
summarise, R2=cor(o, f)^2)
## da$A R2
## 1 37.5 0.7851939
## 2 50 0.8522193
## 3 62.5 0.9457677
##-----------------------------------------------------------------------------
## Os pontos de máximo.
## Buscas essas ocorrências de texto.
oc <- paste0("^A", levels(da$A), ":")
## Estimativas e suas posições de ordem de efeito.
c3 <- coef(m3)
a <- m3$assign
cbind(c3, a)
## c3 a
## A37.5 13.8704427563 1
## A50 16.6522875839 1
## A62.5 15.3039401494 1
## I(pota^2) -0.0008561256 2
## bloc1 0.5556840320 3
## bloc2 1.6216840320 3
## bloc3 -0.3049826347 3
## bloc4 -0.8876493013 3
## A37.5:pota 0.2129359777 4
## A50:pota 0.2209687363 4
## A62.5:pota 0.2762725227 4
## Estimativas do ponto de máximo.
ptmax <- -0.5*c3[a==4]/c3[a==2]
xyplot(fit~pota, 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)+
layer(
panel.abline(v=ptmax, lty=2))
##-----------------------------------------------------------------------------
## E se ajustar um modelo linear-platô? E se for um quadrático-platô?