Estatística Aplicada à Ciência do Solo github.com/walmes/EACS

Propriedades Físicas do Solo para Manejos da Cultura do Eucalipto

Descrição e Análise Exploratória

Os dados no data frame rp_eucal são resultados de um experimento em delineamento de blocos casualizados onde estudou-se o efeito de 4 tipos de manejo sobre as propriedades físicas do solo sob o cultivo de eucalipto: densidade do solo, resistência à penetração, umidade em 0 e 6 kPa. Em cada parcela, amostras de solo indeformadas foram extraídas de 9 camadas de 5 cm de espessura a partir da superfície no mesmo ponto (medidas repetidas), sendo assim, explorados 45 cm de profundidade.

# Nomes curtos são mais fáceis de manuzear.
rp <- rp_eucal
rp$dumid <- rp$umid0 - rp$umid6
str(rp)

## 'data.frame':    216 obs. of  8 variables:
##  $ manejo: Factor w/ 4 levels "BS","BL","BSL",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ cam   : num  5 5 5 5 5 5 10 10 10 10 ...
##  $ bloc  : Factor w/ 6 levels "I","II","III",..: 1 2 3 4 5 6 1 2 3 4 ...
##  $ umid0 : num  0.36 0.37 0.36 0.33 0.42 0.36 0.33 0.34 0.35 0.32 ...
##  $ umid6 : num  0.26 0.23 0.19 0.19 0.24 0.21 0.2 0.18 0.16 0.25 ...
##  $ rp    : num  1.72 1.74 1.81 1.02 0.77 2.06 1.12 1.28 0.86 1.5 ...
##  $ dens  : num  1.68 1.68 1.62 1.58 1.5 1.63 1.65 1.67 1.48 1.64 ...
##  $ dumid : num  0.1 0.14 0.17 0.14 0.18 0.15 0.13 0.16 0.19 0.07 ...

# Visualização do desenho experimental.
xtabs(~manejo + cam, data = rp)

##       cam
## manejo 5 10 15 20 25 30 35 40 45
##    BS  6  6  6  6  6  6  6  6  6
##    BL  6  6  6  6  6  6  6  6  6
##    BSL 6  6  6  6  6  6  6  6  6
##    S   6  6  6  6  6  6  6  6  6

Tabela 1: Nomenclatura usada para os níveis de manejo de solo, descrição, sistema (Sist), uso de braquiária nas entrelinhas (Braq), porcentagem de adubação no sulco (Sulc) e a lanço (Lanç).

Nomenclatura	Descrição	Sist	Braq	Sulc	Lanç
BS	Manejo conservacionista do solo com cultivo de braquiária nas entrelinhas e adubação total (100%) no sulco de plantio.	Conserv.	Sim	100	0
BL	Manejo conservacionista do solo com cultivo de braquiária nas entrelinhas e adubação total (100%) a lanço.	Conserv.	Sim	0	100
BSL	Manejo conservacionista do solo com cultivo de braquiária nas entrelinhas e adubação metade (50%) no sulco de plantio e metade (50%) a lanço (BSL = BS + BL).	Conserv.	Sim	50	50
S	Manejo convencional do solo sem o cultivo de braquiária na entrelinhas e adubação total (100%) no sulco de plantio (testemunha).	Convenc.	Não	100	0

TODO modelo: considerar erros correlacionados entre camadas de uma parcela. Blocos podem ser de efeito aleatório.

Figura 1: Diagrama de dispersão da umidade (g g\(^{-1}\)) em 0 e 6 kPa em função da profundidade no perfil do solo (cm) para cada manejo do solo. A linhas suaves foram adicionadas para marcar a tendência sobre a média dos dados.

Pela Figura 1, a umidade do solo em 0 kPa parece não depender da profundidade do solo e ser a mesma entre os manejos. Visualmente, o valor médio seria de 0.37 g g\(^{-1}\). A umidade em 6 kPa tem o mesmo comportamento que a em 0 kPa. A diferença de umidade mantém-se praticamente constante em um valor de 0.15 ao longo da profundidade do solo e é a mesma entre os manejos.

Para o manejo BSL, na profundidade 45 cm, tem-se uma observação de umidade em 0 kPa e uma em 6 kPa atípicas, com valores acima de 0.45 de umidade. Essa observação pertence ao bloco II.

# Cela dos possíveis outliers.
subset(rp, manejo == "BSL" & cam == 45)

##     manejo cam bloc umid0 umid6   rp dens dumid
## 157    BSL  45    I  0.39  0.24 0.68 1.53  0.15
## 158    BSL  45   II  0.55  0.46 1.02 1.52  0.09
## 159    BSL  45  III  0.26  0.09 0.35 1.44  0.17
## 160    BSL  45   IV  0.43  0.20 0.56 1.46  0.23
## 161    BSL  45    V  0.34  0.23 1.23 1.63  0.11
## 162    BSL  45   VI  0.36  0.26 0.85 1.56  0.10

Figura 2: Diagrama de dispersão de cada uma das respostas em função da profundidade no perfil do solo (cm) para cada manejo do solo. A linhas suaves foram adicionadas para marcar a tendência sobre a média dos dados. Pontos com as mesmas cores pertencem ao mesmo bloco.

A Figura 2 exibe cada uma das respostas, incluindo a diferença de umidade entre 0 e 6 kPa (dumid), em função da profundidade para cada manejo. Visualmente, o maior efeito da profundudade e diferença entre manejos está na densidade do solo e na resistência à penetração, sendo o manejo BSL o mais distinto dos demais.

Método de Análise dos Dados

Conforme o delineamento considerado para realização do experimento e aquisição dos dados, o modelo proposto é

\[ \mu_{ijk} = \mu + \text{BLC}_i + \text{MAN}_j \times s(\text{CAM}_k) + u_{ij} + e_{ijk} \]

em que

• \(\mu_{ijk}\) é o valor esperado para a resposta em uma parcela no bloco \(i\) com o manejo \(j\) tomada na camada \(k\);
• \(\mu\) é uma constante pertencente à todas as observações que, se e somente se na ausência de efeito dos demais termos do modelo, corresponde ao valor esparado da resposta;
• \(\text{BLC}_i\) representa o efeito do bloco \(i\);
• \(\text{MAN}_j \times s(\text{CAM}_k)\) representa o efeito do manejo \(j\), o efeito da camada do solo \(k\) a interação entre estes termos. Por camada ser um fator métrico, para representar seu efeito é usado uma função contínua \(s()\) que pode ser um polinômio ou um base splines;
• \(u_{ij}\) é o erro experimental em nível de parcela que se supõe ser independente e ter distribuição de média 0 e variância constante.
• \(e_{ijk}\) é o erro experimental em nível de amostra de solo que se supõe ser independente e ter distribuição de média 0 e variância constante. A suposição de independencia pode não ser verificada haja visto que amostras de solo próximas tenham correlação.

O efeito de bloco (\(BLC\)) foi considerado como aleatório. Para o efeito de camada considerou-se base splines cúbico de 3 graus de liberdade.

Para o erro em nível de amostra, considerou-se inicialmente uma estrutura de correlação exponencial \[ \rho(d) = \exp\left\{\frac{-d}{\lambda}\right\} \] em que \(d\) é a distância entre duas camadas e \(\lambda\) é o parâmetro relacionado ao alcance. Na prática, o alcançe é a distância na qual a correlação é 5% que corresponde à \(d = 3\rho \approx 0.05\).

Para acomodar o efeito de parcela, foi criada a variável ue concatenando os rótulos de bloco e manejo. Para fazer análise de modelos mistos com a função lme() é recomendado criar um data frame de classe groupedData.

Para cada resposta, três variações do modelo acima descrito foram ajustados:

• m0: modelo em que \(e_{ijk}\) é considerado independente dentro da parcela \(ij\);
• m1: modelo em que \(e_{ijk}\) é considerado autocorrelacionado dentro da parcela \(ij\) descrito pela função exponencial. O modelo m0 está encaixado no modelo m1;
• m2: modelo igual ao m1 porém com efeito de camada descrito pelo polinômio de grau 1 (reta). O modelo m2 está encaixado no modelo m2.

O modelo foi ajustado aos dados pelo método da máxima verossimilhança (ML). Testes entre modelos encaixados foram feitos pelo teste da razão de verossimilhanças considerando um nível de significância de 5%. O resultado de ajuste é apresentado pelo gráfico com as curvas ajustadas acompanhadas das bandas de confiança.

# Ordena e cria UE para representar parcelas (unidades experimentais).
rp <- rp[with(rp, order(bloc, manejo, cam)), ]
rp$ue <- with(rp, interaction(bloc, manejo, drop = TRUE))
rownames(rp) <- NULL
# nlevels(rp$ue)

# Talvez um outlier na umidade. Obs 63.
# subset(rp, umid0 > 0.5 & umid6 > 0.4)

# Cria o `groupedData` para ter `augPred` do ajuste.
rpg <- groupedData(rp ~ cam | ue,
                   data = rp,
                   order.groups = FALSE)

Umidade em 0 kPa

# Modelo que assume independência entre camadas.
m0 <- lme(umid0 ~ manejo * ns(cam),
          data = rpg[-63, ],
          random = ~1 | bloc/ue,
          method = "ML")

# Modelo que expecifica correlação exponencial.
m1 <- update(m0,
             correlation = corExp(value = 5, form = ~cam))
anova(m0, m1)

##    Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## m0     1 11 -885.7386 -848.6615 453.8693                        
## m1     2 12 -896.7342 -856.2865 460.3671 1 vs 2 12.99565   3e-04

# Comparação visual dos ajustes.
# plot(comparePred(m0, m1), layout = c(4, NA))

# Modelo reduzido em termos fixos.
m2 <- update(m1, fixed = . ~ manejo * cam)
anova(m1, m2)

##    Model df       AIC       BIC   logLik
## m1     1 12 -896.7342 -856.2865 460.3671
## m2     2 12 -896.7342 -856.2865 460.3671

Para a resposta umid0, foi verificado que o modelo que assume correlação nula entre camadas (m0) foi rejeitado em favor do modelo que especifica uma correlação do tipo exponencial (m1). O modelo com efeito linear para a camada (m2) mostrou-se tão adequado quanto o modelo com base splines (m1). Neste modelo final (m2) não foi verificado, pela tabela de testes de Wald¹, efeito de manejos e camada do solo.

# Resultados pelo modelo final.
anova(m2)

##             numDF denDF   F-value p-value
## (Intercept)     1   187 11669.435  <.0001
## manejo          3    15     0.266  0.8486
## cam             1   187     0.272  0.6028
## manejo:cam      3   187     0.217  0.8847

summary(m2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: rpg[-63, ] 
##         AIC       BIC   logLik
##   -896.7342 -856.2865 460.3671
## 
## Random effects:
##  Formula: ~1 | bloc
##         (Intercept)
## StdDev: 0.004654178
## 
##  Formula: ~1 | ue %in% bloc
##          (Intercept)   Residual
## StdDev: 1.475662e-06 0.02956587
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~cam | bloc/ue 
##  Parameter estimate(s):
##    range 
## 4.267685 
## Fixed effects: umid0 ~ manejo + cam + manejo:cam 
##                    Value   Std.Error  DF  t-value p-value
## (Intercept)    0.3635329 0.011090569 187 32.77856  0.0000
## manejoBL      -0.0148198 0.015443507  15 -0.95961  0.3525
## manejoBSL     -0.0036543 0.015496270  15 -0.23582  0.8168
## manejoS       -0.0100951 0.015443507  15 -0.65368  0.5232
## cam           -0.0000289 0.000379967 187 -0.07600  0.9395
## manejoBL:cam   0.0003646 0.000537354 187  0.67853  0.4983
## manejoBSL:cam -0.0000181 0.000543281 187 -0.03336  0.9734
## manejoS:cam    0.0001611 0.000537354 187  0.29972  0.7647
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.2632625 -0.6875996 -0.0630175  0.6154105  3.4004376 
## 
## Number of Observations: 215
## Number of Groups: 
##         bloc ue %in% bloc 
##            6           24

# Estimativa de alcance.
r <- coef(m2$modelStruct$corStruct, unconstrained = FALSE)

# Um alcance de 95% (correlação 5%) está na distância de 3 * r.
3 * r

##    range 
## 12.80305

Figura 3: Correlação entre camadas em função da distância (cm) entre elas no perfil do solo.

A estimativa de alcance² foi de 12.8030545 cm (Figura 3). O gráfico mostra o decaimento da correlação entre camadas em função da distância entre elas segundo a função exponencial considerada no ajuste.

Figura 4: Valores observados de umidade do solo (g g\(^{-1}\)) em 0 kPa em função da profundidade do solo (cm) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%).

Na Figura 4 estão os valores observados e ajustados pelo modelo. As bandas de confiança podem conter uma linha horizontal o que uma confirmalção visual de que o modelo nulo (reta horizontal ou apenas intercepto) não é um modelo rejeitado para descrever o comportamento do a umidade do solo em relação à profundidade.

Umidade em 6 kPa

# Modelo que assume independência entre camadas.
m0 <- lme(umid6 ~ manejo * bs(cam),
          data = rpg[-63, ],
          random = ~1 | bloc/ue,
          method = "ML")

# Modelo que expecifica correlação exponencial.
m1 <- update(m0,
             correlation = corExp(value = 2, form = ~cam))
anova(m0, m1)

##    Model df       AIC       BIC   logLik   Test     L.Ratio p-value
## m0     1 19 -820.0369 -755.9948 429.0184                           
## m1     2 20 -818.0465 -750.6338 429.0233 1 vs 2 0.009660833  0.9217

# Comparação visual dos ajustes.
# plot(comparePred(m0, m1), layout = c(4, NA))

# Modelo reduzido em termos fixos.
m2 <- update(m1, fixed = . ~ manejo * cam)
anova(m1, m2)

##    Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## m1     1 20 -818.0465 -750.6338 429.0233                        
## m2     2 12 -826.2502 -785.8026 425.1251 1 vs 2 7.796312  0.4536

Para a resposta umid6, assim como ocorre para umid0, foi verificado que o modelo que assume correlação nula entre camadas (m0) foi rejeitado em favor do modelo que especifica uma correlação do tipo exponencial (m1). O modelo com efeito linear para a camada (m2) mostrou-se tão adequado quanto o modelo com base splines (m1). Neste modelo final (m2) não foi verificado à 5% (mas foi à 10%), pela tabela de testes de Wald, efeito de manejos e camada do solo.

# Resultados pelo modelo final.
anova(m2)

##             numDF denDF   F-value p-value
## (Intercept)     1   187 1321.1439  <.0001
## manejo          3    15    1.5170  0.2508
## cam             1   187    1.9907  0.1599
## manejo:cam      3   187    1.6733  0.1742

summary(m2)

## Linear mixed-effects model fit by maximum likelihood
##  Data: rpg[-63, ] 
##         AIC       BIC   logLik
##   -826.2502 -785.8026 425.1251
## 
## Random effects:
##  Formula: ~1 | bloc
##         (Intercept)
## StdDev: 0.009712083
## 
##  Formula: ~1 | ue %in% bloc
##         (Intercept)   Residual
## StdDev:  0.01427053 0.03138025
## 
## Correlation Structure: Exponential spatial correlation
##  Formula: ~cam | bloc/ue 
##  Parameter estimate(s):
##    range 
## 1.638731 
## Fixed effects: umid6 ~ manejo + cam + manejo:cam 
##                     Value   Std.Error  DF   t-value p-value
## (Intercept)    0.21550037 0.012154388 187 17.730253  0.0000
## manejoBL      -0.03090931 0.016211157  15 -1.906669  0.0759
## manejoBSL     -0.03905405 0.016262737  15 -2.401444  0.0297
## manejoS       -0.01645470 0.016211157  15 -1.015023  0.3262
## cam           -0.00040723 0.000347687 187 -1.171260  0.2430
## manejoBL:cam   0.00075401 0.000491703 187  1.533470  0.1269
## manejoBSL:cam  0.00101789 0.000497856 187  2.044548  0.0423
## manejoS:cam    0.00086151 0.000491703 187  1.752100  0.0814
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.534675713 -0.471287440 -0.002513793  0.578471935  3.018267539 
## 
## Number of Observations: 215
## Number of Groups: 
##         bloc ue %in% bloc 
##            6           24

# Estimativa de alcance.
r <- coef(m2$modelStruct$corStruct, unconstrained = FALSE)

# Um alcance de 95% (correlação 5%) está na distância de 3 * r.
3 * r

##    range 
## 4.916194

Figura 5: Correlação entre camadas em função da distância (cm) entre elas no perfil do solo.

A estimativa de alcance foi de 4.916194 cm (Figura 5). O gráfico mostra o decaimento da correlação entre camadas em função da distância entre elas segundo a função exponencial considerada no ajuste.

Figura 6: Valores observados de umidade do solo (g g\(^{-1}\)) em 6 kPa em função da profundidade do solo (cm) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%).

Na Figura 6 estão os valores observados e ajustados pelo modelo. Exceto para o manejo BSL, as bandas de confiança podem conter uma linha horizontal. Para o BSL verifica-se um efeito significativo de camada.

Diferença de Umidade entre 0 e 6 kPa

# Modelo que assume independência entre camadas.
m0 <- lme(dumid ~ manejo * bs(cam),
          data = rpg,
          random = ~1 | bloc/ue,
          method = "ML")

# Modelo que expecifica correlação exponencial.
m1 <- update(m0,
             correlation = corExp(value = 5, form = ~cam))
anova(m0, m1)

##    Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## m0     1 19 -788.1014 -723.9711 413.0507                        
## m1     2 20 -787.7677 -720.2621 413.8838 1 vs 2 1.666299  0.1968

# Comparação visual dos ajustes.
# plot(comparePred(m0, m1), layout = c(4, NA))

# Modelo reduzido em termos fixos.
m2 <- update(m0, fixed = . ~ manejo * cam)
anova(m1, m2)

##    Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## m1     1 20 -787.7677 -720.2621 413.8838                        
## m2     2 11 -784.3167 -747.1886 403.1584 1 vs 2 21.45096  0.0108

Para a resposta dumid, ao contrário de umid0 e umid6, foi verificado que o modelo que assume correlação nula entre camadas (m0) não foi rejeitado em favor do modelo que especifica uma correlação do tipo exponencial (m1). O modelo com efeito linear para a camada (m2) não se mostrou superior ao modelo com base splines (m1). Neste modelo final (m0), foi detectado interação entre manejo e camada do solo.

# Resultados pelo modelo final.
anova(m0)

##                numDF denDF  F-value p-value
## (Intercept)        1   180 924.7991  <.0001
## manejo             3    15   0.7134  0.5591
## bs(cam)            3   180   0.9726  0.4069
## manejo:bs(cam)     9   180   2.4708  0.0111

summary(m0)

## Linear mixed-effects model fit by maximum likelihood
##  Data: rpg 
##         AIC       BIC   logLik
##   -788.1014 -723.9711 413.0507
## 
## Random effects:
##  Formula: ~1 | bloc
##          (Intercept)
## StdDev: 1.854727e-06
## 
##  Formula: ~1 | ue %in% bloc
##         (Intercept)   Residual
## StdDev:  0.02200506 0.03266135
## 
## Fixed effects: dumid ~ manejo * bs(cam) 
##                          Value  Std.Error  DF   t-value p-value
## (Intercept)         0.14777778 0.01587525 180  9.308689  0.0000
## manejoBL            0.00319865 0.02245100  15  0.142473  0.8886
## manejoBSL           0.06249158 0.02245100  15  2.783466  0.0139
## manejoS             0.00203704 0.02245100  15  0.090733  0.9289
## bs(cam)1           -0.01901395 0.03951766 180 -0.481151  0.6310
## bs(cam)2            0.06290524 0.03013678 180  2.087325  0.0383
## bs(cam)3            0.00010101 0.01889130 180  0.005347  0.9957
## manejoBL:bs(cam)1   0.09389557 0.05588641 180  1.680115  0.0947
## manejoBSL:bs(cam)1 -0.12171931 0.05588641 180 -2.177977  0.0307
## manejoS:bs(cam)1    0.01008124 0.05588641 180  0.180388  0.8571
## manejoBL:bs(cam)2  -0.09785901 0.04261984 180 -2.296090  0.0228
## manejoBSL:bs(cam)2 -0.04816525 0.04261984 180 -1.130113  0.2599
## manejoS:bs(cam)2   -0.04884506 0.04261984 180 -1.146064  0.2533
## manejoBL:bs(cam)3   0.01804714 0.02671633 180  0.675510  0.5002
## manejoBSL:bs(cam)3 -0.06306397 0.02671633 180 -2.360503  0.0193
## manejoS:bs(cam)3   -0.01659933 0.02671633 180 -0.621318  0.5352
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.70541000 -0.63681371  0.01178798  0.59659032  3.38556100 
## 
## Number of Observations: 216
## Number of Groups: 
##         bloc ue %in% bloc 
##            6           24

L <- rbind(c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
           c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0),
           c(0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0),
           c(0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1))
rownames(L) <- levels(rp$manejo)

# Teste para o efeito de profundidade dentro de cada camada.
summary(glht(m0, linfct = L))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lme.formula(fixed = dumid ~ manejo * bs(cam), data = rpg, random = ~1 | 
##     bloc/ue, method = "ML")
## 
## Linear Hypotheses:
##          Estimate Std. Error z value Pr(>|z|)    
## BS == 0   0.04399    0.04871   0.903 0.838885    
## BL == 0   0.05808    0.04871   1.192 0.654184    
## BSL == 0 -0.18896    0.04871  -3.879 0.000419 ***
## S == 0   -0.01137    0.04871  -0.233 0.998839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Figura 7: Valores observados para a diferença entre umidade do solo em 0 e 6 kPa (g g\(^{-1}\)) em função da profundidade do solo (cm) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%).

Na Figura 7 estão os valores observados e ajustados pelo modelo. Exceto para o manejo BSL, as bandas de confiança podem conter uma linha horizontal. Para o BSL verifica-se um efeito significativo de camada descrito pelo base splines.

Densidade do Solo

# Modelo que assume independência entre camadas.
m0 <- lme(dens ~ manejo * bs(cam),
          data = rpg,
          random = ~1 | bloc/ue,
          method = "ML")

# Modelo que expecifica correlação exponencial.
m1 <- update(m0,
             correlation = corExp(value = 5, form = ~cam))
anova(m0, m1)

##    Model df       AIC       BIC   logLik   Test     L.Ratio p-value
## m0     1 19 -553.9514 -489.8211 295.9757                           
## m1     2 20 -551.9514 -484.4458 295.9757 1 vs 2 1.46656e-10       1

# Comparação visual dos ajustes.
# plot(comparePred(m0, m1), layout = c(4, NA))

# Modelo reduzido em termos fixos.
m2 <- update(m0, fixed = . ~ manejo * cam)
anova(m1, m2)

##    Model df       AIC       BIC   logLik   Test  L.Ratio p-value
## m1     1 20 -551.9514 -484.4458 295.9757                        
## m2     2 11 -531.1674 -494.0394 276.5837 1 vs 2 38.78397  <.0001

Para a resposta dens, assim como ocorre com dumid, foi verificado que o modelo que assume correlação nula entre camadas (m0) não foi rejeitado em favor do modelo que especifica uma correlação do tipo exponencial (m1). O modelo com efeito linear para a camada (m2) não se mostrou superior ao modelo com base splines (m1). Neste modelo final (m0), foi detectado interação entre manejo e camada do solo.

# Resultados pelo modelo final.
anova(m0)

##                numDF denDF  F-value p-value
## (Intercept)        1   180 63786.54  <.0001
## manejo             3    15     0.59  0.6331
## bs(cam)            3   180    12.44  <.0001
## manejo:bs(cam)     9   180     2.47  0.0111

summary(m0)

## Linear mixed-effects model fit by maximum likelihood
##  Data: rpg 
##         AIC       BIC   logLik
##   -553.9514 -489.8211 295.9757
## 
## Random effects:
##  Formula: ~1 | bloc
##          (Intercept)
## StdDev: 1.356123e-06
## 
##  Formula: ~1 | ue %in% bloc
##         (Intercept)   Residual
## StdDev:  0.02182441 0.05877596
## 
## Fixed effects: dens ~ manejo * bs(cam) 
##                         Value  Std.Error  DF  t-value p-value
## (Intercept)         1.5973569 0.02489236 180 64.17056  0.0000
## manejoBL           -0.0070202 0.03520312  15 -0.19942  0.8446
## manejoBSL          -0.1152525 0.03520312  15 -3.27393  0.0051
## manejoS            -0.0447475 0.03520312  15 -1.27112  0.2230
## bs(cam)1            0.0636000 0.07111427 180  0.89434  0.3723
## bs(cam)2           -0.1221763 0.05423284 180 -2.25281  0.0255
## bs(cam)3           -0.0602694 0.03399596 180 -1.77284  0.0779
## manejoBL:bs(cam)1  -0.0785265 0.10057077 180 -0.78081  0.4359
## manejoBSL:bs(cam)1  0.2478307 0.10057077 180  2.46424  0.0147
## manejoS:bs(cam)1    0.1240131 0.10057077 180  1.23309  0.2191
## manejoBL:bs(cam)2   0.1119224 0.07669682 180  1.45928  0.1462
## manejoBSL:bs(cam)2  0.1313324 0.07669682 180  1.71236  0.0886
## manejoS:bs(cam)2    0.0622062 0.07669682 180  0.81107  0.4184
## manejoBL:bs(cam)3   0.0003030 0.04807755 180  0.00630  0.9950
## manejoBSL:bs(cam)3  0.1030303 0.04807755 180  2.14300  0.0335
## manejoS:bs(cam)3    0.0774747 0.04807755 180  1.61145  0.1088
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.63678001 -0.59719794  0.01121062  0.63810960  2.68200552 
## 
## Number of Observations: 216
## Number of Groups: 
##         bloc ue %in% bloc 
##            6           24

L <- rbind(c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
           c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0),
           c(0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0),
           c(0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1))
rownames(L) <- levels(rp$manejo)

# Teste para o efeito de profundidade dentro de cada camada.
summary(glht(m0, linfct = L))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lme.formula(fixed = dens ~ manejo * bs(cam), data = rpg, random = ~1 | 
##     bloc/ue, method = "ML")
## 
## Linear Hypotheses:
##          Estimate Std. Error z value Pr(>|z|)    
## BS == 0  -0.11885    0.08766  -1.356 0.537102    
## BL == 0  -0.08515    0.08766  -0.971 0.800122    
## BSL == 0  0.36335    0.08766   4.145 0.000136 ***
## S == 0    0.14485    0.08766   1.652 0.339338    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Figura 8: Valores observados para a densidade do solo (Mg m\(^{-3}\)) em função da profundidade do solo (cm) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%).

Na Figura 8 estão os valores observados e ajustados pelo modelo. Exceto para o manejo BSL, as bandas de confiança podem conter uma linha horizontal. Para o BSL verifica-se um efeito significativo de camada descrito pelo base splines.

Resistência à Penetração

# Modelo que assume independência entre camadas.
m0 <- lme(rp ~ manejo * bs(cam),
          data = rpg,
          random = ~1 | bloc/ue,
          method = "ML")

# Modelo que expecifica correlação exponencial.
m1 <- update(m0,
             correlation = corExp(value = 5, form = ~cam))
anova(m0, m1)

##    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m0     1 19 274.3968 338.5271 -118.1984                        
## m1     2 20 263.0088 330.5144 -111.5044 1 vs 2 13.38793   3e-04

# Comparação visual dos ajustes.
# plot(comparePred(m0, m1), layout = c(4, NA))

# Modelo reduzido em termos fixos.
m2 <- update(m0, fixed = . ~ manejo * cam)
anova(m1, m2)

##    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m1     1 20 263.0088 330.5144 -111.5044                        
## m2     2 11 279.8954 317.0235 -128.9477 1 vs 2 34.88657   1e-04

Para a resposta rp, foi verificado que o modelo que assume correlação nula entre camadas (m0) foi rejeitado em favor do modelo que especifica uma correlação do tipo exponencial (m1). O modelo com efeito linear para a camada (m2) não se mostrou superior ao modelo com base splines (m1). Neste modelo final (m0), foi detectado interação entre manejo e camada do solo.

# Resultados pelo modelo final.
anova(m0)

##                numDF denDF  F-value p-value
## (Intercept)        1   180 396.9692  <.0001
## manejo             3    15   0.3794  0.7692
## bs(cam)            3   180   7.7807  0.0001
## manejo:bs(cam)     9   180   2.5507  0.0088

summary(m0)

## Linear mixed-effects model fit by maximum likelihood
##  Data: rpg 
##        AIC      BIC    logLik
##   274.3968 338.5271 -118.1984
## 
## Random effects:
##  Formula: ~1 | bloc
##          (Intercept)
## StdDev: 2.341164e-05
## 
##  Formula: ~1 | ue %in% bloc
##         (Intercept)  Residual
## StdDev:   0.1768159 0.3949478
## 
## Fixed effects: rp ~ manejo * bs(cam) 
##                         Value Std.Error  DF   t-value p-value
## (Intercept)         1.4101852 0.1724357 180  8.178034  0.0000
## manejoBL            0.1726263 0.2438609  15  0.707888  0.4899
## manejoBSL          -0.5529798 0.2438609  15 -2.267603  0.0386
## manejoS            -0.2069865 0.2438609  15 -0.848789  0.4093
## bs(cam)1           -0.5906729 0.4778557 180 -1.236090  0.2180
## bs(cam)2           -0.8280092 0.3644202 180 -2.272128  0.0243
## bs(cam)3           -0.6394613 0.2284374 180 -2.799284  0.0057
## manejoBL:bs(cam)1  -0.6600000 0.6757900 180 -0.976635  0.3301
## manejoBSL:bs(cam)1  0.4412714 0.6757900 180  0.652971  0.5146
## manejoS:bs(cam)1   -0.0713799 0.6757900 180 -0.105624  0.9160
## manejoBL:bs(cam)2   0.2352862 0.5153679 180  0.456540  0.6486
## manejoBSL:bs(cam)2  1.5523264 0.5153679 180  3.012074  0.0030
## manejoS:bs(cam)2    0.3705767 0.5153679 180  0.719053  0.4730
## manejoBL:bs(cam)3  -0.0292929 0.3230593 180 -0.090674  0.9279
## manejoBSL:bs(cam)3  0.5512121 0.3230593 180  1.706226  0.0897
## manejoS:bs(cam)3    0.2800337 0.3230593 180  0.866818  0.3872
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.2963317 -0.6023723 -0.1081951  0.4749592  4.1769385 
## 
## Number of Observations: 216
## Number of Groups: 
##         bloc ue %in% bloc 
##            6           24

L <- rbind(c(0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
           c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0),
           c(0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0),
           c(0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1))
rownames(L) <- levels(rp$manejo)

# Teste para o efeito de profundidade dentro de cada camada.
summary(glht(m0, linfct = L))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lme.formula(fixed = rp ~ manejo * bs(cam), data = rpg, random = ~1 | 
##     bloc/ue, method = "ML")
## 
## Linear Hypotheses:
##          Estimate Std. Error z value Pr(>|z|)    
## BS == 0   -2.0581     0.5890  -3.494   0.0019 ** 
## BL == 0   -2.5122     0.5890  -4.265 7.99e-05 ***
## BSL == 0   0.4867     0.5890   0.826   0.8777    
## S == 0    -1.4789     0.5890  -2.511   0.0473 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Figura 9: Valores observados para a resistência a penetração do solo (MPa) em função da profundidade do solo (cm) para cada manejo com a curva ajustada acompanhada das bandas de confiança (95%).

Na Figura 9 estão os valores observados e ajustados pelo modelo. Apenas para o manejo BSL, as bandas de confiança pode conter uma linha horizontal. Para os demais manejos verifica-se um efeito significativo de camada descrito pelo base splines.

Informações da Sessão

## Atualizado em 27 de dezembro de 2016.
## 
## R version 3.3.2 (2016-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets 
## [7] base     
## 
## other attached packages:
##  [1] multcomp_1.4-6      TH.data_1.0-7       MASS_7.3-45        
##  [4] survival_2.39-5     mvtnorm_1.0-5       nlme_3.1-128       
##  [7] plyr_1.8.4          reshape2_1.4.1      captioner_2.2.3    
## [10] EACS_0.0-2          wzRfun_0.65         latticeExtra_0.6-28
## [13] RColorBrewer_1.1-2  lattice_0.20-34     rmarkdown_1.0      
## [16] knitr_1.13         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.5      magrittr_1.5     roxygen2_5.0.1  
##  [4] devtools_1.12.0  highr_0.6        stringr_1.0.0   
##  [7] tools_3.3.2      rpanel_1.1-3     grid_3.3.2      
## [10] withr_1.0.2      htmltools_0.3.5  yaml_2.1.13     
## [13] digest_0.6.9     Matrix_1.2-7.1   formatR_1.4     
## [16] codetools_0.2-15 memoise_1.0.0    evaluate_0.9    
## [19] sandwich_2.3-4   doBy_4.5-15      stringi_1.1.1   
## [22] methods_3.3.2    zoo_1.7-13

---

1. Em modelos mistos, o teste para os efeitos fixos é o teste F de Wald e não o teste F de um quadro de ANOVA. A interpretação permanece a mesma.↩

2. Alcance é definida como sendo a distância para uma correlação de 5%.↩

---

Estatística Aplicada à Ciência do Solo

github.com/walmes/EACS