Atributos químicos do solo sob cultivo de abacaxi em função da aplicação de gesso e cobertura do solo
Walmes Marques Zeviani,
Milson Evaldo Serafim,
Definições da sessão
##-----------------------------------------------------------------------------
## Definições da sessão.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape")
sapply(pkg, require, character.only=TRUE)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
## Loading required package: doBy
## Loading required package: survival
## Loading required package: splines
## Loading required package: MASS
## Loading required package: multcomp
## Loading required package: mvtnorm
## Loading required package: reshape
## Loading required package: plyr
##
## Attaching package: 'reshape'
##
## The following objects are masked from 'package:plyr':
##
## rename, round_any
## lattice latticeExtra doBy multcomp reshape
## TRUE TRUE TRUE TRUE TRUE
##-----------------------------------------------------------------------------
## Funções extras.
u <- Sys.info()["user"]; u
## user
## "walmes"
if(u=="walmes"){
source("/home/walmes/Dropbox/Public/apc.R")
source("/home/walmes/Dropbox/Public/equallevels.R")
source("/home/walmes/Dropbox/Public/segplotby.R")
} else {
source("http://dl.dropboxusercontent.com/u/48140237/apc.R")
source("http://dl.dropboxusercontent.com/u/48140237/equallevels.R")
source("http://dl.dropboxusercontent.com/u/48140237/segplotby.R")
}
##-----------------------------------------------------------------------------
## Gráficos monocromáticos.
## names(lattice.options())
trellis.device(color=FALSE)
Leitura e análise preliminar dos dados
##-----------------------------------------------------------------------------
## Ler.
list.files(pattern="\\.txt$")
## [1] "abacaxisolo.txt" "resultados.txt"
url <- "http://www.leg.ufpr.br/~walmes/data/abacaxisolo.txt"
da <- read.table("abacaxisolo.txt", header=TRUE, sep="\t")
str(da)
## 'data.frame': 64 obs. of 12 variables:
## $ cober: Factor w/ 2 levels "com cobertura",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ ges : Factor w/ 2 levels "com gesso","sem gesso": 2 2 2 2 2 2 2 2 2 2 ...
## $ vari : Factor w/ 4 levels "Havai","IAC Fantástico",..: 4 4 4 4 1 1 1 1 2 2 ...
## $ blo : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ph : num 4.8 5.6 5 5 4.5 4.8 4.8 5.1 4.7 4.4 ...
## $ p : num 3.16 3.74 14.81 33.97 3.53 ...
## $ k : num 0.57 0.54 0.31 0.68 0.65 0.68 0.39 0.72 0.84 0.52 ...
## $ ca : num 1.1 1.5 1.5 1.3 1.2 1.1 1.2 1.1 1.1 1.1 ...
## $ mg : num 1.3 1.1 1 0.8 1 1.2 1.4 1.5 1.3 1.2 ...
## $ mo : num 13.8 14.8 16 14.1 20.7 ...
## $ ctc : num 2.97 3.24 2.91 2.98 2.95 3.18 3.09 3.42 3.34 2.92 ...
## $ v : num 56.2 57.5 57.5 57.2 57.9 ...
head(da)
## cober ges vari blo ph p k ca mg mo ctc v
## 1 sem cobertura sem gesso Pérola 1 4.8 3.16 0.57 1.1 1.3 13.82 2.97 56.18
## 2 sem cobertura sem gesso Pérola 2 5.6 3.74 0.54 1.5 1.1 14.76 3.24 57.48
## 3 sem cobertura sem gesso Pérola 3 5.0 14.81 0.31 1.5 1.0 16.01 2.91 57.52
## 4 sem cobertura sem gesso Pérola 4 5.0 33.97 0.68 1.3 0.8 14.13 2.98 57.25
## 5 sem cobertura sem gesso Havai 1 4.5 3.53 0.65 1.2 1.0 20.72 2.95 57.86
## 6 sem cobertura sem gesso Havai 2 4.8 33.97 0.68 1.1 1.2 14.44 3.18 58.94
##-----------------------------------------------------------------------------
## Editar.
da$blo <- factor(da$blo)
##-----------------------------------------------------------------------------
## Resumo descritivo.
with(da, ftable(vari~cober+ges))
## vari Havai IAC Fantástico Imperial Pérola
## cober ges
## com cobertura com gesso 4 4 4 4
## sem gesso 4 4 4 4
## sem cobertura com gesso 4 4 4 4
## sem gesso 4 4 4 4
##-----------------------------------------------------------------------------
## Presença de missings.
colSums(apply(da[,-c(1:4)], 2, is.na))
## ph p k ca mg mo ctc v
## 0 0 0 0 5 0 0 0
## lapply(da[,-c(1:4)],
## function(i){
## with(da, xtabs(!is.na(i)~vari+interaction(ges, cober)))
## })
pH
##-----------------------------------------------------------------------------
## Preliminar.
xyplot(ph~cober|vari, groups=ges, data=da,
jitter.x=TRUE, type=c("p","a"), layout=c(NA,1))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(ph~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: ph
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 0.039 0.013 0.19 0.9032
## vari 3 0.439 0.146 2.12 0.1110
## ges 1 0.214 0.214 3.10 0.0852 .
## cober 1 0.544 0.544 7.87 0.0074 **
## vari:ges 3 0.297 0.099 1.43 0.2460
## vari:cober 3 0.159 0.053 0.77 0.5178
## ges:cober 1 0.056 0.056 0.82 0.3710
## vari:ges:cober 3 0.137 0.046 0.66 0.5811
## Residuals 45 3.108 0.069
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Modelo com interações abandonadas.
m1 <- update(m0, .~blo+vari+ges+cober)
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: ph ~ blo + vari + ges + cober
## Model 2: ph ~ blo + vari * ges * cober
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 3.76
## 2 45 3.11 10 0.649 0.94 0.51
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## LSmeans(m0, effect="cober")
lsmc <- LSmeans(m1, effect="cober")
lsmc <- cbind(lsmc$grid, lsmc$coef)
lsmc <- equallevels(lsmc, da)
lsmc
## cober estimate se df t.stat p.value lwr upr
## 1 com cobertura 4.922 0.0462 55 106.5 2.106e-65 4.829 5.014
## 2 sem cobertura 4.738 0.0462 55 102.5 1.702e-64 4.645 4.830
## LSmeans(m0, effect="ges")
lsmg <- LSmeans(m1, effect="ges")
lsmg <- cbind(lsmg$grid, lsmg$coef)
lsmg <- equallevels(lsmg, da)
lsmg
## ges estimate se df t.stat p.value lwr upr
## 1 com gesso 4.772 0.0462 55 103.3 1.146e-64 4.679 4.864
## 2 sem gesso 4.888 0.0462 55 105.8 3.091e-65 4.795 4.980
## LSmeans(m0, effect="vari")
lsmv <- LSmeans(m1, effect="vari")
lsmv <- cbind(lsmv$grid, lsmv$coef)
lsmv <- equallevels(lsmv, da)
lsmv
## vari estimate se df t.stat p.value lwr upr
## 1 Havai 4.850 0.06534 55 74.22 7.823e-57 4.719 4.981
## 2 IAC Fantástico 4.713 0.06534 55 72.12 3.745e-56 4.582 4.843
## 3 Imperial 4.813 0.06534 55 73.65 1.194e-56 4.682 4.943
## 4 Pérola 4.944 0.06534 55 75.66 2.757e-57 4.813 5.075
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression("pH")
p1 <-
segplot(cober~lwr+upr, centers=estimate,
data=lsmc, draw=FALSE,
xlab=xlab, ylab="Cobertura do solo")
p2 <-
segplot(ges~lwr+upr, centers=estimate,
data=lsmg, draw=FALSE,
xlab=xlab, ylab="Aplicação de gesso")
p3 <-
segplot(vari~lwr+upr, centers=estimate,
data=lsmv, draw=FALSE,
xlab=xlab, ylab="Variedade de abacaxi")
plot(p1, split=c(1,1,1,3), more=TRUE)
plot(p2, split=c(1,2,1,3), more=TRUE)
plot(p3, split=c(1,3,1,3), more=FALSE)
Fósforo (P)
##-----------------------------------------------------------------------------
## Preliminar.
xyplot(log(p)~cober|vari, groups=ges, data=da,
jitter.x=TRUE, type=c("p","a"), layout=c(NA,1))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(p~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
MASS::boxcox(m0)
abline(v=1/3, col=2)
m0 <- lm(p^(1/3)~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: p^(1/3)
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 3.0 1.0 1.81 0.16
## vari 3 1.6 0.5 0.98 0.41
## ges 1 10.8 10.8 19.91 5.4e-05 ***
## cober 1 59.6 59.6 109.69 1.2e-13 ***
## vari:ges 3 1.4 0.5 0.85 0.48
## vari:cober 3 0.6 0.2 0.36 0.78
## ges:cober 1 0.5 0.5 0.99 0.33
## vari:ges:cober 3 1.7 0.6 1.04 0.38
## Residuals 45 24.4 0.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Modelo com interações abandonadas.
m1 <- update(m0, .~blo+vari+ges+cober)
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: p^(1/3) ~ blo + vari + ges + cober
## Model 2: p^(1/3) ~ blo + vari * ges * cober
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 28.6
## 2 45 24.4 10 4.2 0.77 0.65
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## LSmeans(m0, effect="cober")
lsmc <- LSmeans(m1, effect="cober")
lsmc <- cbind(lsmc$grid, lsmc$coef)
lsmc <- equallevels(lsmc, da)
lsmc[,c("estimate","lwr","upr")] <- lsmc[,c("estimate","lwr","upr")]^3
lsmc
## cober estimate se df t.stat p.value lwr upr
## 1 com cobertura 74.51 0.1276 55 32.99 6.356e-38 61.737 88.93
## 2 sem cobertura 11.83 0.1276 55 17.86 1.480e-24 8.275 16.27
## LSmeans(m0, effect="ges")
lsmg <- LSmeans(m1, effect="ges")
lsmg <- cbind(lsmg$grid, lsmg$coef)
lsmg <- equallevels(lsmg, da)
lsmg[,c("estimate","lwr","upr")] <- lsmg[,c("estimate","lwr","upr")]^3
lsmg
## ges estimate se df t.stat p.value lwr upr
## 1 com gesso 48.79 0.1276 55 28.64 9.799e-35 39.25 59.77
## 2 sem gesso 22.72 0.1276 55 22.20 3.964e-29 17.10 29.44
## LSmeans(m0, effect="vari")
lsmv <- LSmeans(m1, effect="vari")
lsmv <- cbind(lsmv$grid, lsmv$coef)
lsmv <- equallevels(lsmv, da)
lsmv[,c("estimate","lwr","upr")] <- lsmv[,c("estimate","lwr","upr")]^3
lsmv
## vari estimate se df t.stat p.value lwr upr
## 1 Havai 28.42 0.1804 55 16.92 1.863e-23 19.47 39.77
## 2 IAC Fantástico 40.10 0.1804 55 18.97 8.416e-26 28.68 54.19
## 3 Imperial 38.39 0.1804 55 18.70 1.684e-25 27.32 52.10
## 4 Pérola 30.51 0.1804 55 17.32 6.230e-24 21.10 42.37
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression("Fósforo"~(mg~dm^{-3}))
p1 <-
segplot(cober~lwr+upr, centers=estimate,
data=lsmc, draw=FALSE,
xlab=xlab, ylab="Cobertura do solo")
p2 <-
segplot(ges~lwr+upr, centers=estimate,
data=lsmg, draw=FALSE,
xlab=xlab, ylab="Aplicação de gesso")
p3 <-
segplot(vari~lwr+upr, centers=estimate,
data=lsmv, draw=FALSE,
xlab=xlab, ylab="Variedade de abacaxi")
plot(p1, split=c(1,1,1,3), more=TRUE)
plot(p2, split=c(1,2,1,3), more=TRUE)
plot(p3, split=c(1,3,1,3), more=FALSE)
Potássio (K)
##-----------------------------------------------------------------------------
## Preliminar.
xyplot(k~cober|vari, groups=ges, data=da,
jitter.x=TRUE, type=c("p","a"), layout=c(NA,1))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(k~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: k
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 0.096 0.0321 1.63 0.196
## vari 3 0.060 0.0199 1.01 0.397
## ges 1 0.050 0.0501 2.54 0.118
## cober 1 0.024 0.0244 1.24 0.271
## vari:ges 3 0.008 0.0028 0.14 0.935
## vari:cober 3 0.152 0.0506 2.57 0.066 .
## ges:cober 1 0.031 0.0311 1.58 0.216
## vari:ges:cober 3 0.004 0.0013 0.06 0.979
## Residuals 45 0.886 0.0197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Modelo com interações abandonadas.
m1 <- update(m0, .~blo+vari+ges+cober)
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: k ~ blo + vari + ges + cober
## Model 2: k ~ blo + vari * ges * cober
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 1.081
## 2 45 0.886 10 0.195 0.99 0.47
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## LSmeans(m0, effect="cober")
lsmc <- LSmeans(m1, effect="cober")
lsmc <- cbind(lsmc$grid, lsmc$coef)
lsmc <- equallevels(lsmc, da)
lsmc
## cober estimate se df t.stat p.value lwr upr
## 1 com cobertura 0.5594 0.02479 55 22.57 1.758e-29 0.5097 0.6090
## 2 sem cobertura 0.5984 0.02479 55 24.14 5.998e-31 0.5488 0.6481
## LSmeans(m0, effect="ges")
lsmg <- LSmeans(m1, effect="ges")
lsmg <- cbind(lsmg$grid, lsmg$coef)
lsmg <- equallevels(lsmg, da)
lsmg
## ges estimate se df t.stat p.value lwr upr
## 1 com gesso 0.5509 0.02479 55 22.23 3.740e-29 0.5013 0.6006
## 2 sem gesso 0.6069 0.02479 55 24.48 2.962e-31 0.5572 0.6565
## LSmeans(m0, effect="vari")
lsmv <- LSmeans(m1, effect="vari")
lsmv <- cbind(lsmv$grid, lsmv$coef)
lsmv <- equallevels(lsmv, da)
lsmv
## vari estimate se df t.stat p.value lwr upr
## 1 Havai 0.5356 0.03505 55 15.28 1.909e-21 0.4654 0.6059
## 2 IAC Fantástico 0.6212 0.03505 55 17.72 2.125e-24 0.5510 0.6915
## 3 Imperial 0.5850 0.03505 55 16.69 3.484e-23 0.5147 0.6553
## 4 Pérola 0.5737 0.03505 55 16.37 8.502e-23 0.5035 0.6440
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression("Potássio"~(cmol[c]~dm^{-3}))
p1 <-
segplot(cober~lwr+upr, centers=estimate,
data=lsmc, draw=FALSE,
xlab=xlab, ylab="Cobertura do solo")
p2 <-
segplot(ges~lwr+upr, centers=estimate,
data=lsmg, draw=FALSE,
xlab=xlab, ylab="Aplicação de gesso")
p3 <-
segplot(vari~lwr+upr, centers=estimate,
data=lsmv, draw=FALSE,
xlab=xlab, ylab="Variedade de abacaxi")
plot(p1, split=c(1,1,1,3), more=TRUE)
plot(p2, split=c(1,2,1,3), more=TRUE)
plot(p3, split=c(1,3,1,3), more=FALSE)
Cálcio (Ca)
##-----------------------------------------------------------------------------
## Preliminar.
xyplot(ca~cober|vari, groups=ges, data=da,
jitter.x=TRUE, type=c("p","a"), layout=c(NA,1))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(ca~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
MASS::boxcox(m0)
abline(v=1/3, col=2)
m0 <- lm(ca^(1/3)~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: ca^(1/3)
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 0.010 0.003 0.56 0.643
## vari 3 0.040 0.013 2.32 0.088 .
## ges 1 0.004 0.004 0.66 0.422
## cober 1 0.948 0.948 164.94 <2e-16 ***
## vari:ges 3 0.034 0.011 1.98 0.130
## vari:cober 3 0.028 0.009 1.64 0.193
## ges:cober 1 0.035 0.035 6.08 0.018 *
## vari:ges:cober 3 0.027 0.009 1.58 0.207
## Residuals 45 0.259 0.006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Quadro de anova particionado.
m0 <- aov(ca^(1/3)~blo+vari*ges+ges/cober+vari:ges:cober, data=da)
c0 <- coef(m0)
a <- m0$assign
c0 <- c0[a==5]
s <- paste0("ges", levels(da$ges), ":cober")
s <- sapply(s, grep, x=names(c0), simplify=FALSE)
names(s) <- levels(da$ges)
summary(m0, split=list("ges:cober"=s))
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 0.010 0.003 0.56 0.643
## vari 3 0.040 0.013 2.32 0.088 .
## ges 1 0.004 0.004 0.66 0.422
## vari:ges 3 0.034 0.011 1.98 0.130
## ges:cober 2 0.983 0.491 85.51 4.7e-16 ***
## ges:cober: com gesso 1 0.673 0.673 117.17 4.1e-14 ***
## ges:cober: sem gesso 1 0.310 0.310 53.85 3.2e-09 ***
## vari:ges:cober 6 0.056 0.009 1.61 0.166
## Residuals 45 0.259 0.006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##-----------------------------------------------------------------------------
## Modelo com interações abandonadas.
m1 <- update(m0, .~blo+vari+ges*cober)
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: ca^(1/3) ~ blo + vari + ges + cober + ges:cober
## Model 2: ca^(1/3) ~ blo + vari * ges + ges/cober + vari:ges:cober
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 54 0.348
## 2 45 0.259 9 0.0898 1.74 0.11
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## LSmeans(m0, effect="cober")
lsmc <- LSmeans(m1, effect=c("cober","ges"))
lsmc <- cbind(lsmc$grid, lsmc$coef)
lsmc <- equallevels(lsmc, da)
lsmc[,c("estimate","lwr","upr")] <- lsmc[,c("estimate","lwr","upr")]^3
lsmc
## cober ges estimate se df t.stat p.value lwr upr
## 1 com cobertura com gesso 2.181 0.02008 54 64.58 8.253e-53 1.9838 2.390
## 2 sem cobertura com gesso 1.020 0.02008 54 50.13 5.749e-47 0.9025 1.147
## 3 com cobertura sem gesso 1.882 0.02008 54 61.48 1.127e-51 1.7040 2.072
## 4 sem cobertura sem gesso 1.118 0.02008 54 51.69 1.134e-47 0.9932 1.254
## LSmeans(m0, effect="vari")
lsmv <- LSmeans(m1, effect="vari")
lsmv <- cbind(lsmv$grid, lsmv$coef)
lsmv <- equallevels(lsmv, da)
lsmv[,c("estimate","lwr","upr")] <- lsmv[,c("estimate","lwr","upr")]^3
lsmv
## vari estimate se df t.stat p.value lwr upr
## 1 Havai 1.660 0.02008 54 58.97 1.043e-50 1.497 1.836
## 2 IAC Fantástico 1.388 0.02008 54 55.55 2.490e-49 1.243 1.544
## 3 Imperial 1.465 0.02008 54 56.56 9.560e-50 1.315 1.627
## 4 Pérola 1.484 0.02008 54 56.80 7.623e-50 1.332 1.647
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression("Cálcio"~(cmol[c]~dm^{-3}))
p1 <-
segplot(cober~lwr+upr|ges, centers=estimate,
data=lsmc, draw=FALSE,
xlab=xlab, ylab="Cobertura do solo")
l <- levels(lsmc$cober)
key <- list(corner=c(0.97,0.97),
type="o", divide=1,
lines=list(lty=1, pch=1:length(l)),
text=list(l))
p2 <-
segplot(ges~lwr+upr, centers=estimate,
data=lsmc, draw=FALSE, key=key,
by=lsmc$cober, f=0.1, pch=as.integer(lsmc$cober),
xlab=xlab, ylab="Aplicação de gesso",
panel=segplot.by)
p3 <-
segplot(vari~lwr+upr, centers=estimate,
data=lsmv, draw=FALSE,
xlab=xlab, ylab="Variedade de abacaxi")
plot(p1, split=c(1,1,1,3), more=TRUE)
plot(p2, split=c(1,2,1,3), more=TRUE)
plot(p3, split=c(1,3,1,3), more=FALSE)
Magnésio (Mg)
##-----------------------------------------------------------------------------
## Preliminar.
xyplot(mg~cober|vari, groups=ges, data=da[!is.na(da$mg),],
jitter.x=TRUE, type=c("p","a"), layout=c(NA,1))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(mg~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)
layout(1)
## Quadro de anova.
anova(m0)
## Analysis of Variance Table
##
## Response: mg
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 0.18 0.059 0.59 0.62433
## vari 3 0.14 0.045 0.45 0.71693
## ges 1 0.32 0.322 3.21 0.08067 .
## cober 1 1.36 1.365 13.64 0.00066 ***
## vari:ges 3 0.30 0.101 1.01 0.39663
## vari:cober 3 0.89 0.297 2.96 0.04351 *
## ges:cober 1 0.11 0.114 1.14 0.29232
## vari:ges:cober 3 0.65 0.216 2.16 0.10753
## Residuals 40 4.00 0.100
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Quadro de anova particionado.
m0 <- aov(mg~blo+vari*ges+vari/cober+vari:ges:cober, data=da)
c0 <- coef(m0)
a <- m0$assign
c0 <- c0[a==5]
s <- paste0("vari", levels(da$vari), ":cober")
s <- sapply(s, grep, x=names(c0), simplify=FALSE)
names(s) <- levels(da$vari)
summary(m0, split=list("vari:cober"=s))
## Df Sum Sq Mean Sq F value Pr(>F)
## blo 3 0.18 0.059 0.59 0.62433
## vari 3 0.14 0.045 0.45 0.71693
## ges 1 0.32 0.322 3.21 0.08067 .
## vari:ges 3 0.30 0.099 0.99 0.40584
## vari:cober 4 2.26 0.565 5.65 0.00106 **
## vari:cober: Havai 1 1.51 1.511 15.09 0.00038 ***
## vari:cober: IAC Fantástico 1 0.68 0.681 6.80 0.01276 *
## vari:cober: Imperial 1 0.00 0.002 0.02 0.88948
## vari:cober: Pérola 1 0.07 0.067 0.67 0.41686
## vari:ges:cober 4 0.76 0.191 1.91 0.12819
## Residuals 40 4.00 0.100
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 5 observations deleted due to missingness
##-----------------------------------------------------------------------------
## Modelo com interações abandonadas.
m1 <- update(m0, .~blo+vari*cober+ges)
anova(m1, m0)
## Analysis of Variance Table
##
## Model 1: mg ~ blo + vari + cober + ges + vari:cober
## Model 2: mg ~ blo + vari * ges + vari/cober + vari:ges:cober
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 47 5.2
## 2 40 4.0 7 1.2 1.71 0.13
##-----------------------------------------------------------------------------
## Comparações múltiplas.
## LSmeans(m0, effect="cober")
lsmc <- LSmeans(m1, effect=c("cober","vari"))
lsmc <- cbind(lsmc$grid, lsmc$coef)
lsmc <- equallevels(lsmc, da)
lsmc
## cober vari estimate se df t.stat p.value lwr upr
## 1 com cobertura Havai 0.6675 0.1375 47 4.854 1.376e-05 0.3909 0.9441
## 2 sem cobertura Havai 1.3000 0.1176 47 11.051 1.153e-14 1.0634 1.5366
## 3 com cobertura IAC Fantástico 0.9250 0.1176 47 7.864 4.074e-10 0.6884 1.1616
## 4 sem cobertura IAC Fantástico 1.3375 0.1176 47 11.370 4.326e-15 1.1009 1.5741
## 5 com cobertura Imperial 1.0092 0.1264 47 7.984 2.692e-10 0.7549 1.2635
## 6 sem cobertura Imperial 1.0375 0.1176 47 8.820 1.567e-11 0.8009 1.2741
## 7 com cobertura Pérola 1.0175 0.1375 47 7.399 2.034e-09 0.7409 1.2941
## 8 sem cobertura Pérola 1.1875 0.1176 47 10.095 2.358e-13 0.9509 1.4241
## LSmeans(m0, effect="ges")
lsmg <- LSmeans(m1, effect="ges")
lsmg <- cbind(lsmg$grid, lsmg$coef)
lsmg <- equallevels(lsmg, da)
lsmg
## ges estimate se df t.stat p.value lwr upr
## 1 com gesso 1.118 0.06443 47 17.35 4.528e-22 0.9885 1.248
## 2 sem gesso 1.002 0.05994 47 16.72 2.046e-21 0.8817 1.123
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression("Magnésio"~(cmol[c]~dm^{-3}))
p1 <-
segplot(cober~lwr+upr|vari, centers=estimate,
data=lsmc, draw=FALSE,
xlab=xlab, ylab="Cobertura do solo")
l <- levels(lsmc$cober)
key <- list(corner=c(0.97,0.97),
type="o", divide=1,
lines=list(lty=1, pch=1:length(l)),
text=list(l))
p2 <-
segplot(vari~lwr+upr, centers=estimate,
data=lsmc, draw=FALSE, key=key,
by=lsmc$cober, f=0.1, pch=as.integer(lsmc$cober),
xlab=xlab, ylab="Variedade de abacaxi",
panel=segplot.by)
p3 <-
segplot(ges~lwr+upr, centers=estimate,
data=lsmg, draw=FALSE,
xlab=xlab, ylab="Aplicação de gesso")
plot(p1, split=c(1,1,1,3), more=TRUE)
plot(p2, split=c(1,2,1,3), more=TRUE)
plot(p3, split=c(1,3,1,3), more=FALSE)
Matéria orgânica (MO)
##-----------------------------------------------------------------------------
## Preliminar.
xyplot(mo~cober|vari, groups=ges, data=da[!is.na(da$mg),],
jitter.x=TRUE, type=c("p","a"), layout=c(NA,1))
##-----------------------------------------------------------------------------
## Modelo.
m0 <- lm(mo~blo+vari*ges*cober, data=da)
## Diagnóstico dos resíduos.
par(mfrow=c(2,2))
plot(m0)