Análise de experimento em parcela subsubdividida na profundidade
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "doBy", "multcomp",
"reshape", "plyr", "nlme", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra doBy multcomp reshape
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr nlme wzRfun
## TRUE TRUE TRUE
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 grid stats graphics grDevices utils datasets
## [9] base
##
## other attached packages:
## [1] wzRfun_0.3 nlme_3.1-117 plyr_1.8.1 reshape_0.8.5
## [5] multcomp_1.3-7 TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11
## [9] MASS_7.3-34 survival_2.37-7 gridExtra_0.9.1 latticeExtra_0.6-26
## [13] RColorBrewer_1.0-5 lattice_0.20-29 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 htmltools_0.2.6 Matrix_1.1-4
## [6] Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1 yaml_2.1.13
## [11] zoo_1.7-10
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/sistema_gesso_solo.txt",
header=TRUE, sep="\t")
## str(da)
##-----------------------------------------------------------------------------
## Editar.
da <- transform(da,
bloco=as.factor(bloco),
gesso=as.factor(gesso),
Prof=as.factor(prof),
parc=interaction(bloco, sistema),
subp=interaction(bloco, sistema, gesso))
str(da)
## 'data.frame': 80 obs. of 21 variables:
## $ sistema: Factor w/ 2 levels "PC","PD": 2 2 2 2 2 2 2 2 2 2 ...
## $ gesso : Factor w/ 2 levels "0","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ prof : int 5 5 5 5 10 10 10 10 15 15 ...
## $ bloco : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
## $ Al : num 0 0 0 0 0 0 0 0 0 0 ...
## $ MO : num 60.4 42.4 51.7 52.4 31.7 39.9 32.3 40.7 28.7 31 ...
## $ pHKCl : num 5.2 5.6 6.4 4.2 6.7 6.2 6.6 3.9 6.2 5.3 ...
## $ pHCaCl : num 5.5 5.7 6.1 4.3 6.2 6.3 6.5 4 6.2 5.8 ...
## $ pHH2O : num 6 6.4 7 5 7.2 7.1 7.4 4.6 7.2 6.4 ...
## $ Ca : num 69.9 77.3 91.4 79.4 89.2 86.9 95.5 90.5 76.6 57.7 ...
## $ Mg : num 32.3 30.3 28.6 30.4 36 32.6 26 31 28.3 19.6 ...
## $ Hal : int 45 29 18 30 16 19 16 17 20 34 ...
## $ S : num 8.2 2.8 3.8 4.9 7.2 2.7 7.2 2.7 6.6 7.8 ...
## $ P : int 19 17 22 16 9 4 8 8 4 2 ...
## $ K : num 4.8 4.7 8 6.9 1.7 2.3 3.3 2 0.9 1.2 ...
## $ SB : num 107 112 128 117 127 ...
## $ CTC : num 152 141 146 147 143 ...
## $ SAT : int 70 79 88 80 89 87 89 88 84 70 ...
## $ Prof : Factor w/ 5 levels "5","10","15",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ parc : Factor w/ 8 levels "1.PC","2.PC",..: 5 6 7 8 5 6 7 8 5 6 ...
## $ subp : Factor w/ 16 levels "1.PC.0","2.PC.0",..: 5 6 7 8 5 6 7 8 5 6 ...
##-----------------------------------------------------------------------------
## Ver.
da$y <- da$P
xyplot(y~prof|sistema, groups=gesso, data=da, type=c("p","a"),
auto.key=TRUE, jitter.x=TRUE)
##-----------------------------------------------------------------------------
## Ajuste do modelo de parcela subdivididas (assume independência entre
## observações de uma mesma parcela na profundidade).
m0 <- lme(y~bloco+sistema*gesso*Prof,
random=~1|parc/subp, data=da,
method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])
bs <- unlist(ranef(m0)[2])
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r), qqmath(bp), qqmath(bs))
## Embora esteja considerando só a parte de efeito fixo serve de indicação.
MASS::boxcox(lm(formula(m0), data=da))
## Embora esteja considerando só a parte de efeito aleatório um nível
## anterior ao de resíduo, serve de indicação.
MASS::boxcox(lm(y~subp, data=da))
##-----------------------------------------------------------------------------
## Modelo com a variável transformada. Refazer a análise dos resíduos.
m0 <- update(m0, log(y)~.)
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])
bs <- unlist(ranef(m0)[2])
grid.arrange(ncol=2,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r), qqmath(bp), qqmath(bs))
## Teste de Wald para os termos de efeito fixo.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 48 1103.3533 <.0001
## bloco 3 3 4.2924 0.1312
## sistema 1 3 32.3413 0.0108
## gesso 1 6 0.8181 0.4006
## Prof 4 48 197.1857 <.0001
## sistema:gesso 1 6 0.0078 0.9324
## sistema:Prof 4 48 6.5600 0.0003
## gesso:Prof 4 48 4.3878 0.0042
## sistema:gesso:Prof 4 48 1.4311 0.2381
##-----------------------------------------------------------------------------
## Mas e a correlação serial ao longo do perfil?
plot(ACF(m0))
S <- split(data.frame(r=r, prof=da$prof), da$subp)
S <- lapply(S,
function(i){
i$r[order(i$prof)]
})
S <- do.call(rbind, S)
splom(S, type=c("p","smooth"))
round(cor(S), 2)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.00 0.13 -0.38 -0.42 0.69
## [2,] 0.13 1.00 0.30 -0.10 0.11
## [3,] -0.38 0.30 1.00 -0.05 -0.33
## [4,] -0.42 -0.10 -0.05 1.00 0.16
## [5,] 0.69 0.11 -0.33 0.16 1.00
## Como se comporta dentro de cada sistema de cultivo?
f <- factor(do.call(rbind, strsplit(rownames(S), "\\."))[,2])
by(S, f, cor)
## INDICES: PC
## V1 V2 V3 V4 V5
## V1 1.00000000 0.03150293 -0.6331372 -0.68370638 0.68729735
## V2 0.03150293 1.00000000 0.5394903 -0.24830273 0.17403248
## V3 -0.63313724 0.53949026 1.0000000 0.47557690 -0.12569855
## V4 -0.68370638 -0.24830273 0.4755769 1.00000000 -0.03007535
## V5 0.68729735 0.17403248 -0.1256985 -0.03007535 1.00000000
## -------------------------------------------------------------------
## INDICES: PD
## V1 V2 V3 V4 V5
## V1 1.0000000 0.23038025 -0.1248384 -0.11471696 0.70587232
## V2 0.2303802 1.00000000 0.1788678 -0.02528559 0.08041909
## V3 -0.1248384 0.17886779 1.0000000 -0.43709352 -0.60037585
## V4 -0.1147170 -0.02528559 -0.4370935 1.00000000 0.41162070
## V5 0.7058723 0.08041909 -0.6003758 0.41162070 1.00000000
## Estrutura de correlação não estruturada (tem k*(k-1)/2 parâmetros).
m1 <- update(m0,
correlation=corSymm(form=~as.integer(Prof)|parc/subp))
summary(m1)
## Linear mixed-effects model fit by maximum likelihood
## Data: da
## AIC BIC logLik
## 64.52099 150.2739 3.739507
##
## Random effects:
## Formula: ~1 | parc
## (Intercept)
## StdDev: 0.02511193
##
## Formula: ~1 | subp %in% parc
## (Intercept) Residual
## StdDev: 0.07776823 0.3540443
##
## Correlation Structure: General
## Formula: ~as.integer(Prof) | parc/subp
## Parameter estimate(s):
## Correlation:
## 1 2 3 4
## 2 0.634
## 3 0.053 0.407
## 4 -0.378 -0.198 -0.036
## 5 0.940 0.659 0.129 -0.079
## Fixed effects: log(y) ~ bloco + sistema + gesso + Prof + sistema:gesso + sistema:Prof + gesso:Prof + sistema:gesso:Prof
## Value Std.Error DF t-value p-value
## (Intercept) 3.396698 0.2245999 48 15.12333 0.0000
## bloco2 -0.389777 0.1048206 3 -3.71852 0.0338
## bloco3 -0.152128 0.1048206 3 -1.45132 0.2426
## bloco4 -0.374328 0.1048206 3 -3.57113 0.0375
## sistemaPD -0.257319 0.3043842 3 -0.84538 0.4600
## gesso2 -0.119140 0.3036564 6 -0.39235 0.7084
## Prof10 -0.115323 0.1793306 48 -0.64308 0.5232
## Prof15 -1.552273 0.2886479 48 -5.37774 0.0000
## Prof20 -2.821067 0.3481911 48 -8.10206 0.0000
## Prof25 -2.821067 0.0729021 48 -38.69665 0.0000
## sistemaPD:gesso2 0.619822 0.4294350 6 1.44334 0.1990
## sistemaPD:Prof10 -0.859397 0.2536118 48 -3.38863 0.0014
## sistemaPD:Prof15 -0.491614 0.4082097 48 -1.20432 0.2344
## sistemaPD:Prof20 0.257319 0.4924166 48 0.52256 0.6037
## sistemaPD:Prof25 -0.089254 0.1030991 48 -0.86571 0.3910
## gesso2:Prof10 -0.493850 0.2536118 48 -1.94727 0.0574
## gesso2:Prof15 -0.142202 0.4082097 48 -0.34835 0.7291
## gesso2:Prof20 0.694787 0.4924166 48 1.41097 0.1647
## gesso2:Prof25 0.119140 0.1030991 48 1.15559 0.2536
## sistemaPD:gesso2:Prof10 -0.432368 0.3586613 48 -1.20551 0.2339
## sistemaPD:gesso2:Prof15 -0.878340 0.5772958 48 -1.52147 0.1347
## sistemaPD:gesso2:Prof20 -1.094102 0.6963822 48 -1.57112 0.1227
## sistemaPD:gesso2:Prof25 -0.619822 0.1458042 48 -4.25106 0.0001
## Correlation:
## (Intr) bloco2 bloco3 bloco4 sstmPD gesso2 Prof10 Prof15 Prof20
## bloco2 -0.233
## bloco3 -0.233 0.500
## bloco4 -0.233 0.500 0.500
## sistemaPD -0.678 0.000 0.000 0.000
## gesso2 -0.676 0.000 0.000 0.000 0.499
## Prof10 -0.399 0.000 0.000 0.000 0.295 0.295
## Prof15 -0.643 0.000 0.000 0.000 0.474 0.475 0.611
## Prof20 -0.775 0.000 0.000 0.000 0.572 0.573 0.385 0.564
## Prof25 -0.162 0.000 0.000 0.000 0.120 0.120 0.285 0.285 0.624
## sistemaPD:gesso2 0.478 0.000 0.000 0.000 -0.705 -0.707 -0.209 -0.336 -0.405
## sistemaPD:Prof10 0.282 0.000 0.000 0.000 -0.417 -0.209 -0.707 -0.432 -0.272
## sistemaPD:Prof15 0.454 0.000 0.000 0.000 -0.671 -0.336 -0.432 -0.707 -0.399
## sistemaPD:Prof20 0.548 0.000 0.000 0.000 -0.809 -0.405 -0.272 -0.399 -0.707
## sistemaPD:Prof25 0.115 0.000 0.000 0.000 -0.169 -0.085 -0.201 -0.202 -0.441
## gesso2:Prof10 0.282 0.000 0.000 0.000 -0.208 -0.418 -0.707 -0.432 -0.272
## gesso2:Prof15 0.454 0.000 0.000 0.000 -0.335 -0.672 -0.432 -0.707 -0.399
## gesso2:Prof20 0.548 0.000 0.000 0.000 -0.404 -0.811 -0.272 -0.399 -0.707
## gesso2:Prof25 0.115 0.000 0.000 0.000 -0.085 -0.170 -0.201 -0.202 -0.441
## sistemaPD:gesso2:Prof10 -0.200 0.000 0.000 0.000 0.295 0.295 0.500 0.306 0.192
## sistemaPD:gesso2:Prof15 -0.321 0.000 0.000 0.000 0.474 0.475 0.306 0.500 0.282
## sistemaPD:gesso2:Prof20 -0.388 0.000 0.000 0.000 0.572 0.573 0.192 0.282 0.500
## sistemaPD:gesso2:Prof25 -0.081 0.000 0.000 0.000 0.120 0.120 0.142 0.143 0.312
## Prof25 ssPD:2 sPD:P10 sPD:P15 sPD:P20 sPD:P25 g2:P10 g2:P15
## bloco2
## bloco3
## bloco4
## sistemaPD
## gesso2
## Prof10
## Prof15
## Prof20
## Prof25
## sistemaPD:gesso2 -0.085
## sistemaPD:Prof10 -0.201 0.295
## sistemaPD:Prof15 -0.202 0.475 0.611
## sistemaPD:Prof20 -0.441 0.573 0.385 0.564
## sistemaPD:Prof25 -0.707 0.120 0.285 0.285 0.624
## gesso2:Prof10 -0.201 0.295 0.500 0.306 0.192 0.142
## gesso2:Prof15 -0.202 0.475 0.306 0.500 0.282 0.143 0.611
## gesso2:Prof20 -0.441 0.573 0.192 0.282 0.500 0.312 0.385 0.564
## gesso2:Prof25 -0.707 0.120 0.142 0.143 0.312 0.500 0.285 0.285
## sistemaPD:gesso2:Prof10 0.142 -0.418 -0.707 -0.432 -0.272 -0.201 -0.707 -0.432
## sistemaPD:gesso2:Prof15 0.143 -0.672 -0.432 -0.707 -0.399 -0.202 -0.432 -0.707
## sistemaPD:gesso2:Prof20 0.312 -0.811 -0.272 -0.399 -0.707 -0.441 -0.272 -0.399
## sistemaPD:gesso2:Prof25 0.500 -0.170 -0.201 -0.202 -0.441 -0.707 -0.201 -0.202
## g2:P20 g2:P25 sPD:2:P10 sPD:2:P15 sPD:2:P20
## bloco2
## bloco3
## bloco4
## sistemaPD
## gesso2
## Prof10
## Prof15
## Prof20
## Prof25
## sistemaPD:gesso2
## sistemaPD:Prof10
## sistemaPD:Prof15
## sistemaPD:Prof20
## sistemaPD:Prof25
## gesso2:Prof10
## gesso2:Prof15
## gesso2:Prof20
## gesso2:Prof25 0.624
## sistemaPD:gesso2:Prof10 -0.272 -0.201
## sistemaPD:gesso2:Prof15 -0.399 -0.202 0.611
## sistemaPD:gesso2:Prof20 -0.707 -0.441 0.385 0.564
## sistemaPD:gesso2:Prof25 -0.441 -0.707 0.285 0.285 0.624
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.05309288 -0.62628125 0.01908487 0.49551168 2.21988742
##
## Number of Observations: 80
## Number of Groups:
## parc subp %in% parc
## 8 16
anova(m0, m1)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 26 87.31893 149.2516 -17.659463
## m1 2 36 64.52099 150.2739 3.739507 1 vs 2 42.79794 <.0001
## Estrutura de continous AR(1).
m2 <- update(m0, correlation=corCAR1(form=~prof|parc/subp))
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 26 87.31893 149.2516 -17.65946
## m2 2 27 88.12685 152.4416 -17.06342 1 vs 2 1.192079 0.2749
AIC(m1)
## [1] 64.52099
AIC(m2)
## [1] 88.12685
##-----------------------------------------------------------------------------
## Veja como os resultados mudam completamente.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 48 1103.3533 <.0001
## bloco 3 3 4.2924 0.1312
## sistema 1 3 32.3413 0.0108
## gesso 1 6 0.8181 0.4006
## Prof 4 48 197.1857 <.0001
## sistema:gesso 1 6 0.0078 0.9324
## sistema:Prof 4 48 6.5600 0.0003
## gesso:Prof 4 48 4.3878 0.0042
## sistema:gesso:Prof 4 48 1.4311 0.2381
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 48 23831.520 <.0001
## bloco 3 3 6.390 0.0810
## sistema 1 3 14.009 0.0333
## gesso 1 6 38.130 0.0008
## Prof 4 48 2276.191 <.0001
## sistema:gesso 1 6 26.938 0.0020
## sistema:Prof 4 48 18.882 <.0001
## gesso:Prof 4 48 8.227 <.0001
## sistema:gesso:Prof 4 48 5.384 0.0012
##-----------------------------------------------------------------------------
## Qual foi a correlação estimada?
## Embora seja aceitável a correlação de 63% entre 1-2 e de 40% entre
## 2-3, as correlações negativas com a camada 4 as correlações
## superiores à 60% de entre 1-5 e 2-5 parecem bem estranhas. Será que o
## fato dos sistemas serem diferentes quando ao revolvimento do solo
## (PD=não mexe, PC=mexe) não justifica uma especificação de covariância
## ao longo do perfil que seja separada por sistema?
summary(m1)$modelStruct$corStruct
## Correlation structure of class corSymm representing
## Correlation:
## 1 2 3 4
## 2 0.634
## 3 0.053 0.407
## 4 -0.378 -0.198 -0.036
## 5 0.940 0.659 0.129 -0.079
##-----------------------------------------------------------------------------
## O interesse é saber como o gesso altera o nível dos nutrientes em
## cada sistema de manejo. Então obter valores preditos para
## profundidade em cada nível de gesso e sistema.
X <- LSmatrix(lm(formula(m1), data=da),
effect=c("sistema","gesso","Prof"))
grid <- attr(X, "grid")
## Médias com IC.
ci <- confint(glht(m1, linfct=X), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
grid <- equallevels(grid, da)
grid
## sistema gesso Prof Estimate lwr upr
## 1 PC 0 5 3.167640e+00 2.811560262 3.5237203
## 2 PD 0 5 2.910321e+00 2.554240864 3.2664009
## 3 PC 2 5 3.048500e+00 2.692420005 3.4045800
## 4 PD 2 5 3.411003e+00 3.054922549 3.7670826
## 5 PC 0 10 3.052317e+00 2.696236873 3.4083969
## 6 PD 0 10 1.935601e+00 1.579520494 2.2916805
## 7 PC 2 10 2.439326e+00 2.083246249 2.7954063
## 8 PD 2 10 1.510064e+00 1.153983667 1.8661437
## 9 PC 0 15 1.615367e+00 1.259287033 1.9714471
## 10 PD 0 15 8.664340e-01 0.510353965 1.2225140
## 11 PC 2 15 1.354025e+00 0.997945089 1.7101051
## 12 PD 2 15 3.465736e-01 -0.009506421 0.7026536
## 13 PC 0 20 3.465736e-01 -0.009506421 0.7026536
## 14 PD 0 20 3.465736e-01 -0.009506421 0.7026536
## 15 PC 2 20 9.222199e-01 0.566139852 1.2782999
## 16 PD 2 20 4.479399e-01 0.091859856 0.8040199
## 17 PC 0 25 3.465736e-01 -0.009506421 0.7026536
## 18 PD 0 25 2.720046e-15 -0.356080011 0.3560800
## 19 PC 2 25 3.465736e-01 -0.009506421 0.7026536
## 20 PD 2 25 2.664535e-15 -0.356080011 0.3560800
##-----------------------------------------------------------------------------
## Representação.
## l <- sapply(levels(grid$gesso),
## function(l){
## eval(substitute(expression(gesso~(ton~ha^{-1})),
## list(gesso=l)))
## })
l <- levels(grid$gesso)
key <- list(type="o", divide=1, columns=length(l),
title=expression("Gesso agrícola"~(ton~ha^{-1})),
cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("plot.symbol")$col),
text=list(l))
## Para que os pontos sejam corretamente representados, deve-se
## reordenar os dados.
grid <- grid[with(grid, order(sistema, Prof, gesso)),]
segplot(Prof~lwr+upr|sistema, horizontal=FALSE,
centers=Estimate, groups=grid$gesso, f=0.05, layout=c(NA,1),
data=grid, draw=FALSE, panel=panel.segplotBy, key=key,
pch=as.integer(grid$gesso),
ylab="log Fósforo",
xlab="Camada do solo (cm)")
##-----------------------------------------------------------------------------
## Como ficariam os resultados se fosse assumido independência?
grie <- attr(X, "grid")
ci <- confint(glht(m0, linfct=X), calpha=univariate_calpha())
grie <- cbind(grie, ci$confint)
grie <- equallevels(grie, da)
grie <- grie[with(grie, order(sistema, Prof, gesso)),]
grid$modelo <- "correlação"
grie$modelo <- "independente"
gri <- rbind(grid, grie)
gri$modelo <- factor(gri$modelo)
levels(gri$modelo)
## [1] "correlação" "independente"
l <- levels(gri$modelo)
key <- list(type="o", divide=1, columns=length(l),
title="Modelo",
cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("plot.symbol")$col),
text=list(l))
gri <- gri[with(gri, order(sistema, gesso, Prof, modelo)),]
segplot(Prof~lwr+upr|sistema*gesso, horizontal=FALSE,
centers=Estimate, groups=gri$modelo, f=0.05,
data=gri, draw=FALSE, panel=panel.segplotBy, key=key,
pch=as.integer(gri$modelo),
ylab="log Fósforo",
xlab="Camada do solo (cm)")
##-----------------------------------------------------------------------------
## Ligando as médias para dar inteção de continuidade.
d <- (as.integer(gri$modelo)-median(seq_along(l)))/5
key <- list(type="o", divide=1, columns=length(l),
title="Modelo",
cex.title=1.1,
lines=list(pch=1:length(l),
col=trellis.par.get("superpose.symbol")$col[1:length(l)]),
text=list(l))
xyplot(Estimate~Prof|sistema*gesso, type="o",
data=gri,
groups=modelo,
ly=gri$lwr, uy=gri$upr,
cty="bars", desloc=d, length=0.02,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose,
key=key, pch=as.integer(gri$modelo),
ylab="Fósforo",
xlab="Camada do solo (cm)")
