|
Estatística Aplicada à Ciência do Solo
|
Modelagem da capacidade de suporte de carga de latossolos com diferente mineralogia
Walmes Marques Zeviani & Carla Eloize Carducci
16 de janeiro de 2020
Preconsolidation curve
Exploratory data analysis
#-----------------------------------------------------------------------
# Load packages.
rm(list = objects())
library(EACS)
library(lmerTest)
library(lme4)
library(emmeans)
library(car)
library(tidyverse)
# To compare curves with Wald F test statistic.
compare_curves <- function(index,
levels = NULL,
L = L,
ords = list(ord0, ord1)) {
J <- L[index, , drop = FALSE]
K <- wzRfun::apc(J, lev = levels)
f_test <- apply(K[, , drop = FALSE],
MARGIN = 1,
FUN = function(Ki) {
U <- lapply(ords,
function(ord) {
Ki %*% diag(ord)
})
U <- do.call(rbind, U)
f_test <-
linearHypothesis(model = m1,
hypothesis.matrix = U,
test = "F")
tail(f_test[, c("F", "Pr(>F)")], n = 1)
})
f_test <- do.call(what = rbind, args = f_test)
cbind(contrast = rownames(f_test), f_test)
}
# Text to be used as axis label.
axis_text <- list(
posi = "Position",
solo = "Soil",
prof = "Depth (cm)",
tensao = expression(Log[10] ~ (Psi * "m, kPa") - "Matric tension"),
thetai = expression("Water content" ~ (cm^{3} ~ cm^{-3})),
pcc = expression(Log[10] ~ (sigma * "p, kPa") - "Preconsolidation pressure"),
pr = expression("Penetration resistance" ~ (MPa)),
ds = expression("Bulk density" ~ (g ~ cm^{-3})),
macro = expression("Macroporosity" ~ (cm^3 ~ cm^{-3})),
micro = expression("Microporosity" ~ (cm^3 ~ cm^{-3})),
total = expression("Total porosity" ~ (cm^3 ~ cm^{-3})))# Dataset.
da <- as_tibble(cafe_pedotrans)
str(da)## Classes 'tbl_df', 'tbl' and 'data.frame': 180 obs. of 11 variables:
## $ solo : Factor w/ 2 levels "LVAd","LVd": 2 2 2 2 2 2 2 2 2 2 ...
## $ tensao: int 6 6 6 6 6 6 6 6 6 6 ...
## $ posi : Factor w/ 2 levels "Entrelinha","Linha": 1 1 1 1 1 1 1 1 1 2 ...
## $ prof : Factor w/ 3 levels "0.00-0.05 m",..: 1 1 1 2 2 2 3 3 3 1 ...
## $ rep : int 1 2 3 1 2 3 1 2 3 1 ...
## $ dsi : num 0.84 0.82 0.99 0.84 0.95 0.99 0.83 0.92 0.91 0.86 ...
## $ dsp : num 0.9 0.82 0.99 0.99 0.99 1.04 0.92 0.95 0.97 0.87 ...
## $ thetai: num 0.506 0.481 0.465 0.451 0.421 ...
## $ ppc : num 125.2 86.4 207.5 256.7 120.3 ...
## $ trt : Factor w/ 12 levels "LVAd_Entrelinha_0.00-0.05 m",..: 2 2 2 6 6 6 10 10 10 4 ...
## $ ue : Factor w/ 36 levels "LVAd_Entrelinha_0.00-0.05 m_1",..: 2 14 26 6 18 30 10 22 34 4 ...# Creates variables to represent different experimental units sizes.
da <- da %>%
mutate(ue = NULL,
ue1 = interaction(solo, rep),
ue2 = interaction(solo, rep, posi),
ue3 = interaction(solo, rep, posi, prof))
# Recode factor levels.
levels(da$posi) <- c("Interrow", "Row")
s <- c(KH = "Kaolinitic Haplustox", GA = "Gibbsitic Acrustox")
levels(da$solo) <- names(s)
levels(da$prof)## [1] "0.00-0.05 m" "0.10-0.15 m" "0.20-0.25 m"# Recode levels.
da$trt <- with(da, interaction(solo, posi, as.integer(prof), sep = "·"))
# Numeric summary.
summary(da)## solo tensao posi prof
## KH:90 Min. : 6.0 Interrow:90 0.00-0.05 m:60
## GA:90 1st Qu.: 33.0 Row :90 0.10-0.15 m:60
## Median : 100.0 0.20-0.25 m:60
## Mean : 927.8
## 3rd Qu.:1500.0
## Max. :3000.0
##
## rep dsi dsp thetai
## Min. :1 Min. :0.3300 Min. :0.810 Min. :0.0800
## 1st Qu.:1 1st Qu.:0.8700 1st Qu.:0.930 1st Qu.:0.2951
## Median :2 Median :0.9400 Median :0.990 Median :0.3498
## Mean :2 Mean :0.9418 Mean :1.008 Mean :0.3219
## 3rd Qu.:3 3rd Qu.:1.0000 3rd Qu.:1.080 3rd Qu.:0.3902
## Max. :3 Max. :1.1800 Max. :1.240 Max. :0.5064
##
## ppc trt ue1
## Min. : 51.55 KH·Interrow·1:15 LVAd.1:30
## 1st Qu.: 99.76 GA·Interrow·1:15 LVd.1 :30
## Median :168.53 KH·Row·1 :15 LVAd.2:30
## Mean :214.15 GA·Row·1 :15 LVd.2 :30
## 3rd Qu.:319.83 KH·Interrow·2:15 LVAd.3:30
## Max. :535.24 GA·Interrow·2:15 LVd.3 :30
## (Other) :90
## ue2 ue3
## LVAd.1.Entrelinha:15 LVAd.1.Entrelinha.0.00-0.05 m: 5
## LVd.1.Entrelinha :15 LVd.1.Entrelinha.0.00-0.05 m : 5
## LVAd.2.Entrelinha:15 LVAd.2.Entrelinha.0.00-0.05 m: 5
## LVd.2.Entrelinha :15 LVd.2.Entrelinha.0.00-0.05 m : 5
## LVAd.3.Entrelinha:15 LVAd.3.Entrelinha.0.00-0.05 m: 5
## LVd.3.Entrelinha :15 LVd.3.Entrelinha.0.00-0.05 m : 5
## (Other) :90 (Other) :150# Count of experimental points.
ftable(tensao ~ solo + posi + prof, data = da)## tensao 6 33 100 1500 3000
## solo posi prof
## KH Interrow 0.00-0.05 m 3 3 3 3 3
## 0.10-0.15 m 3 3 3 3 3
## 0.20-0.25 m 3 3 3 3 3
## Row 0.00-0.05 m 3 3 3 3 3
## 0.10-0.15 m 3 3 3 3 3
## 0.20-0.25 m 3 3 3 3 3
## GA Interrow 0.00-0.05 m 3 3 3 3 3
## 0.10-0.15 m 3 3 3 3 3
## 0.20-0.25 m 3 3 3 3 3
## Row 0.00-0.05 m 3 3 3 3 3
## 0.10-0.15 m 3 3 3 3 3
## 0.20-0.25 m 3 3 3 3 3# Creates variables in the log scale.
da$logtensao <- log10(da$tensao)
da$logppc <- log10(da$ppc)
# Experimental units profile in each experimental point.
ggplot(data = da,
mapping = aes(x = tensao, y = ppc, color = rep, group = rep)) +
facet_wrap(facets = ~trt) +
geom_point(show.legend = FALSE) +
geom_line(show.legend = FALSE) +
scale_x_log10() +
scale_y_log10() +
labs(x = axis_text$tensao,
y = axis_text$pcc)
# Mean profile in each experimental point. The curve is the second order
# polynomial model fitted to the sample data only to check if it
# captures properly the trend in data.
ggplot(data = da,
mapping = aes(x = tensao, y = ppc)) +
facet_wrap(facets = ~trt) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
geom_smooth(method = "lm",
formula = y ~ poly(x, degree = 2),
se = FALSE) +
scale_x_log10() +
scale_y_log10() +
labs(x = axis_text$tensao,
y = axis_text$pcc)
# Differnce in position mean profile conditional to other factors.
ggplot(data = da,
mapping = aes(x = tensao, y = ppc, color = posi)) +
facet_grid(facets = prof ~ solo) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
scale_x_log10() +
scale_y_log10() +
labs(x = axis_text$tensao,
y = axis_text$pcc,
color = axis_text$posi)
# Differnce in soil mean profile conditional to other factors.
ggplot(data = da,
mapping = aes(x = tensao, y = ppc, color = solo)) +
facet_grid(facets = prof ~ posi) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
scale_x_log10() +
scale_y_log10() +
labs(x = axis_text$tensao,
y = axis_text$pcc,
color = axis_text$solo)
Model specification and fitting
Based on exploratory data analysis, the effect of log of matric tension on log of preconsolidation pressure (response variable) can be properly described by a second order polynomial in each experimental condition (combination of soil, position and depth levels). So, log of preconsolidation pressure was modelled as a function of soil type, sampling position, soil depth and log of matric tension. All interactions (until those of 4th way) were included in the model. The effect of different sizes experimental units was accounted by random effets terms at intercept level. Likelihood ratio tests were performed to check the relevance random effects terms and Wald F tests were performed check the relevance of fixed effect terms, so those not relevant terms were abandoned getting a simpler model. The complete model is
\[ \begin{align*} y_{ijklm} &= \mu + S_i + P_j + D_k + (\beta_1 T_l + \beta_2 T^2_l) + \\ &\quad (SP)_{ij} + \cdots + (SPD)_{ijk} + \cdots + (\beta_{1ijk} T_l + \beta_{2ijk} T_l^2) + \\ &\quad e_{im} + e_{ijm} + e_{ijkm} + e_{ijklm}\\ e_{im} &\sim \text{Normal}(0, \sigma^2_a) \qquad [\text{6 whole plots}]\\ e_{ijm} &\sim \text{Normal}(0, \sigma^2_b) \qquad [\text{12 sub plots}]\\ e_{ijkm} &\sim \text{Normal}(0, \sigma^2_c) \qquad [\text{36 sub sub plots}]\\ e_{ijklm} &\sim \text{Normal}(0, \sigma^2) \qquad [\text{180 data points}]\\ \end{align*} \] where
- \(y_{ijklm}\) is the response at the \(ijklm\) experimental condition.
- \(\mu\) is a constant that belongs to all experimental condition.
- \(S_i\) accounts for the soil effect.
- \(P_j\) accounts for the sampling position effect.
- \(D_k\) accounts for the soil depth effect.
- \(\beta_1\) and \(\beta_2\) accounts for the linear em quadratic effect of matric tension.
- All two and three way interaction of the above terms are included.
- \(e_{im}\) accounts for the experimental error at whole plots, \(e_{ijm}\) at sub plots, \(e_{ijlm}\) at sub sub plots and \(e_{ijklm}\) at data point level.
This model were fitted by restricted maximum likelihood (REML). Non relevant termos were abandoned to get a simpler model to inspect the research hypothesis.
# 4th degree interaction fixed effects and random effects at intercept.
m0 <- lmer(logppc ~ solo * posi * prof * (logtensao + I(logtensao^2)) +
(1 | ue1) + (1 | ue2) + (1 | ue3),
data = da)## boundary (singular) fit: see ?isSingular# Variance of the random effects.
VarCorr(m0)## Groups Name Std.Dev.
## ue3 (Intercept) 0.026517
## ue2 (Intercept) 0.000000
## ue1 (Intercept) 0.033548
## Residual 0.104107# Test for the fixed effect terms of the model.
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## solo 0.00437 0.00437 1 97.733 0.4031
## posi 0.00854 0.00854 1 128.959 0.7879
## prof 0.10519 0.05259 2 128.959 4.8526
## logtensao 0.14217 0.14217 1 120.000 13.1175
## I(logtensao^2) 0.77136 0.77136 1 120.000 71.1703
## solo:posi 0.04223 0.04223 1 128.959 3.8964
## solo:prof 0.01163 0.00582 2 128.959 0.5367
## posi:prof 0.00919 0.00460 2 128.959 0.4241
## solo:logtensao 0.01946 0.01946 1 120.000 1.7958
## solo:I(logtensao^2) 0.03373 0.03373 1 120.000 3.1123
## posi:logtensao 0.00083 0.00083 1 120.000 0.0768
## posi:I(logtensao^2) 0.00607 0.00607 1 120.000 0.5602
## prof:logtensao 0.12341 0.06170 2 120.000 5.6931
## prof:I(logtensao^2) 0.12517 0.06258 2 120.000 5.7743
## solo:posi:prof 0.04402 0.02201 2 128.959 2.0309
## solo:posi:logtensao 0.01736 0.01736 1 120.000 1.6016
## solo:posi:I(logtensao^2) 0.01139 0.01139 1 120.000 1.0510
## solo:prof:logtensao 0.03073 0.01537 2 120.000 1.4178
## solo:prof:I(logtensao^2) 0.03214 0.01607 2 120.000 1.4829
## posi:prof:logtensao 0.01316 0.00658 2 120.000 0.6071
## posi:prof:I(logtensao^2) 0.01654 0.00827 2 120.000 0.7630
## solo:posi:prof:logtensao 0.05301 0.02651 2 120.000 2.4455
## solo:posi:prof:I(logtensao^2) 0.05128 0.02564 2 120.000 2.3658
## Pr(>F)
## solo 0.5269878
## posi 0.3763748
## prof 0.0092916 **
## logtensao 0.0004299 ***
## I(logtensao^2) 8.627e-14 ***
## solo:posi 0.0505281 .
## solo:prof 0.5859657
## posi:prof 0.6552362
## solo:logtensao 0.1827518
## solo:I(logtensao^2) 0.0802473 .
## posi:logtensao 0.7821859
## posi:I(logtensao^2) 0.4556471
## prof:logtensao 0.0043439 **
## prof:I(logtensao^2) 0.0040335 **
## solo:posi:prof 0.1353884
## solo:posi:logtensao 0.2081328
## solo:posi:I(logtensao^2) 0.3073318
## solo:prof:logtensao 0.2462730
## solo:prof:I(logtensao^2) 0.2310982
## posi:prof:logtensao 0.5466039
## posi:prof:I(logtensao^2) 0.4685294
## solo:posi:prof:logtensao 0.0909915 .
## solo:posi:prof:I(logtensao^2) 0.0982377 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model summary.
# summary(m0)
# No significant 4 degree interaction.
# No significant 3 degree interaction.
# No significant 2 degree with soil.
# No significant 2 degree with posi.
# Fitting a smaller model with relevant terms.
m1 <- update(m0,
. ~ (solo + posi + prof) *
(logtensao + I(logtensao^2)) +
solo:posi +
(1 | ue1) + (1 | ue2) + (1 | ue3))## boundary (singular) fit: see ?isSingularanova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: da
## Models:
## m1: logppc ~ solo + posi + prof + logtensao + I(logtensao^2) + (1 |
## m1: ue1) + (1 | ue2) + (1 | ue3) + solo:logtensao + solo:I(logtensao^2) +
## m1: posi:logtensao + posi:I(logtensao^2) + prof:logtensao + prof:I(logtensao^2) +
## m1: solo:posi
## m0: logppc ~ solo * posi * prof * (logtensao + I(logtensao^2)) +
## m0: (1 | ue1) + (1 | ue2) + (1 | ue3)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 20 -259.85 -195.99 149.93 -299.85
## m0 40 -247.10 -119.38 163.55 -327.10 27.25 20 0.1284anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.00437 0.00437 1 104.976 0.4008 0.5280597
## posi 0.00851 0.00851 1 144.954 0.7802 0.3785342
## prof 0.10484 0.05242 2 144.954 4.8050 0.0095380
## logtensao 0.14217 0.14217 1 134.007 13.0319 0.0004315
## I(logtensao^2) 0.77136 0.77136 1 134.007 70.7059 5.404e-14
## solo:logtensao 0.01946 0.01946 1 134.007 1.7841 0.1839121
## solo:I(logtensao^2) 0.03373 0.03373 1 134.007 3.0920 0.0809610
## posi:logtensao 0.00083 0.00083 1 134.007 0.0763 0.7828295
## posi:I(logtensao^2) 0.00607 0.00607 1 134.007 0.5565 0.4569652
## prof:logtensao 0.12341 0.06170 2 134.007 5.6560 0.0043836
## prof:I(logtensao^2) 0.12517 0.06258 2 134.007 5.7367 0.0040694
## solo:posi 0.05139 0.05139 1 25.998 4.7109 0.0392918
##
## solo
## posi
## prof **
## logtensao ***
## I(logtensao^2) ***
## solo:logtensao
## solo:I(logtensao^2) .
## posi:logtensao
## posi:I(logtensao^2)
## prof:logtensao **
## prof:I(logtensao^2) **
## solo:posi *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From this model, by Wald F tests, there weren’t 4th degree signficant interactions neither 3rd degree. All two way interactions were kept except those involving depth and soil or depth and position. So, depth and soil had only additive effect, the same for depth and position.
Curve comparisons
# Estimated marginal means with covariate set at 1.
emm <- emmeans(m1,
specs = ~solo + posi + prof,
at = list(logtensao = 1))
# Get experimental points and matrix to calculate marginal means.
grid <- attr(emm, "grid")
L <- attr(emm, "linfct")
# Fixed effects estimates.
b <- fixef(m1)
# Zero order terms (intercept).
ord0 <- names(b) %>%
grepl(pattern = "logtensao") %>% `!`()
# First order terms (slope).
ord1 <- names(b) %>%
grepl(pattern = "logtensao$")
# Second order terms (curvatures).
ord2 <- names(b) %>%
grepl(pattern = "I\\(logtensao\\^2\\)$")
# Get the equation coefficients: b0 + b1 * x + b2 * x^2.
grid$b0 <- c(L %*% (b * ord0))
grid$b1 <- c(L %*% (b * ord1))
grid$b2 <- c(L %*% (b * ord2))
grid <- grid %>%
mutate(label = sprintf("'%0.2f' %s '%0.3f' * x %s '%0.4f' * x^2",
b0,
ifelse(b1 >= 0, "+", "-"),
abs(b1),
ifelse(b2 >= 0, "+", "-"),
abs(b2)))
# Table of equation coefficients for each experimental point.
grid %>%
select(solo, posi, prof, label)## solo posi prof label
## 1 KH Interrow 0.00-0.05 m '2.13' - '0.339' * x + '0.1335' * x^2
## 2 GA Interrow 0.00-0.05 m '2.23' - '0.212' * x + '0.0959' * x^2
## 3 KH Row 0.00-0.05 m '2.09' - '0.365' * x + '0.1495' * x^2
## 4 GA Row 0.00-0.05 m '2.11' - '0.238' * x + '0.1119' * x^2
## 5 KH Interrow 0.10-0.15 m '2.24' - '0.330' * x + '0.1186' * x^2
## 6 GA Interrow 0.10-0.15 m '2.34' - '0.203' * x + '0.0810' * x^2
## 7 KH Row 0.10-0.15 m '2.20' - '0.356' * x + '0.1346' * x^2
## 8 GA Row 0.10-0.15 m '2.22' - '0.229' * x + '0.0970' * x^2
## 9 KH Interrow 0.20-0.25 m '1.90' + '0.004' * x + '0.0503' * x^2
## 10 GA Interrow 0.20-0.25 m '2.00' + '0.131' * x + '0.0126' * x^2
## 11 KH Row 0.20-0.25 m '1.87' - '0.022' * x + '0.0662' * x^2
## 12 GA Row 0.20-0.25 m '1.88' + '0.105' * x + '0.0286' * x^2The above table has the estimated curve coefficients for each experimental point resulted of combining soil, position and depth levels (12 equations at all). These fitted curves has equation in the form of \[ y = \beta_0 + \beta_1 x + \beta_2 x^2 \] where \(y\) is the (conditional mean of) log 10 of preconsolidation pressure, \(x\) is the log 10 of matric tension and \(\beta_0\), \(\beta_1\) and \(\beta_2\) are intercept, slope and curvature parameters of the second order polynomial.
# Conditions to get estimated values.
pred <- with(da,
expand.grid(solo = levels(solo),
posi = levels(posi),
prof = levels(prof),
logtensao = seq(0.5, 3.5, length = 31),
KEEP.OUT.ATTRS = FALSE))
# Model matrix for prediction.
X <- model.matrix(formula(terms(m1))[-2], data = pred)
# Checks.
stopifnot(ncol(X) == length(fixef(m1)))
stopifnot(all(colnames(X) == names(fixef(m1))))
# Predicted value.
pred$fit <- X %*% fixef(m1)
# Standard error of predicted values.
pred$se <- sqrt(diag(X %*% tcrossprod(vcov(m1), X)))
# Quantile of t distribution.
qtl <- qt(p = 0.975, df = df.residual(m1))
# Confidence limits for predicted value.
pred$lwr <- pred$fit - qtl * pred$se
pred$upr <- pred$fit + qtl * pred$se
# Adds the equation to annotate on graphics.
grid <- grid %>%
mutate(xpos = 2.5,
ypos = ifelse(solo == "GA", 500, 400) - 30,
hjust = c(0),
vjust = c(0))
# Fitted curves with confidence bands and observed data.
# Soil in evidence.
ggplot(data = da,
mapping = aes(x = tensao, y = ppc, color = solo)) +
facet_grid(facets = prof ~ posi) +
geom_point() +
scale_x_log10() +
scale_y_log10() +
geom_line(data = pred,
mapping = aes(x = 10^logtensao,
y = 10^fit,
color = solo)) +
geom_smooth(data = pred,
mapping = aes(x = 10^logtensao,
y = 10^fit,
ymin = 10^lwr,
ymax = 10^upr),
stat = "identity",
show.legend = FALSE) +
geom_text(data = grid,
mapping = aes(x = xpos,
y = ypos,
label = label,
hjust = hjust,
vjust = vjust),
show.legend = FALSE,
parse = TRUE) +
labs(x = axis_text$tensao,
y = axis_text$pcc,
color = axis_text$solo)
grid <- grid %>%
mutate(xpos = 2.5,
ypos = ifelse(posi == "Interrow", 500, 400) - 30)
# Fitted curves with confidence bands and observed data.
# Position in evidence.
ggplot(data = da,
mapping = aes(x = tensao, y = ppc, color = posi)) +
facet_grid(facets = prof ~ solo) +
geom_point() +
scale_x_log10() +
scale_y_log10() +
geom_line(data = pred,
mapping = aes(x = 10^logtensao,
y = 10^fit,
color = posi)) +
geom_smooth(data = pred,
mapping = aes(x = 10^logtensao,
y = 10^fit,
ymin = 10^lwr,
ymax = 10^upr),
stat = "identity",
show.legend = FALSE) +
geom_text(data = grid,
mapping = aes(x = xpos,
y = ypos,
label = label,
hjust = hjust,
vjust = vjust),
show.legend = FALSE,
parse = TRUE) +
labs(x = axis_text$tensao,
y = axis_text$pcc,
color = axis_text$posi)
The above graphics displays the estimated quadratic curve fitted (solid line) in each combination of soil, position and depth to describe log of preconsolidation pressure as a function of log matric tension. The shaded region wrapping the estimated curve is the confidence band for the fitted conditional mean. The equations inside each cell plot are in log-log.
By visualizing the curves, it can be noticed that curves are more dissimilar between soils at interrow than in row position. Also, curves are more dissimilar between sampling positions for GA soil than for KH.
# Tests of soil curves are equal in each posi:prof combination.
grid %>%
mutate(row = 1:n()) %>%
group_by(posi, prof) %>%
do({
compare_curves(.$row, .$solo, L, list(ord0, ord1, ord2))
}) %>%
as.data.frame()## posi prof contrast F Pr(>F)
## 1 Interrow 0.00-0.05 m KH-GA 8.532725 0.0003851223
## 2 Interrow 0.10-0.15 m KH-GA 8.532725 0.0003851223
## 3 Interrow 0.20-0.25 m KH-GA 8.532725 0.0003851223
## 4 Row 0.00-0.05 m KH-GA 4.252112 0.0139741381
## 5 Row 0.10-0.15 m KH-GA 4.252112 0.0139741381
## 6 Row 0.20-0.25 m KH-GA 4.252112 0.0139741381# Tests of posi curves are equal in each soil:prof combination.
grid %>%
mutate(row = 1:n()) %>%
group_by(solo, prof) %>%
do({
compare_curves(.$row, .$posi, L, list(ord0, ord1, ord2))
}) %>%
as.data.frame()## solo prof contrast F Pr(>F)
## 1 KH 0.00-0.05 m Interrow-Row 2.670002 0.080247817
## 2 KH 0.10-0.15 m Interrow-Row 2.670002 0.080247817
## 3 KH 0.20-0.25 m Interrow-Row 2.670002 0.080247817
## 4 GA 0.00-0.05 m Interrow-Row 5.804931 0.006322117
## 5 GA 0.10-0.15 m Interrow-Row 5.804931 0.006322117
## 6 GA 0.20-0.25 m Interrow-Row 5.804931 0.006322117Wald F test were applied to test with two curves are equal. The null hypothesis states that all three curve parameters (\(\beta_0\), \(\beta_1\), \(\beta_2\)) are simultaneous equal for a pair of curves \(a\) and \(b\) \[ H_0: \begin{bmatrix} \beta_{0a} \\ \beta_{1a} \\ \beta_{2a} \end{bmatrix} = \begin{bmatrix} \beta_{0b} \\ \beta_{1b} \\ \beta_{2b} \end{bmatrix}. \]
The tables above shows the F statistics and associate p-value to test the equality of two curves using the Wald test. Curves from soils were compared for each position and depth combination (6 conditions). The F statistic is the same for all depth at the same position to compare soils because soil and depth had only additive effect. The greater difference, as showed by the graphics, were between soil curves at interrow position (p < 0.01) than row position (p < 0.05).
The same test were applied to test equality of curves between position for each soil and depth combination (6 conditions). The F statistic is the same for all depth at the same soil to compare positions because position and depth had only additive effect. The greater difference, as indicated by the graphics, were between position curves at GA soil (p < 0.01). As suggested by overlapping of confidence bands between positions in the KH soil, according two the F test, these curves are not different at 5% but at 10% (p < 0.10).
Resistance and Water Content
Resistance to penetration
# Data structure.
str(da)## 'data.frame': 60 obs. of 10 variables:
## $ solo : Factor w/ 2 levels "KH","GA": 2 2 2 2 2 2 2 2 2 2 ...
## $ posi : Factor w/ 2 levels "Interrow","Row": 1 1 1 1 1 1 1 1 1 1 ...
## $ pr : num 0.183 0.43 0.633 1.446 0.524 ...
## $ tensao : int 6 6 6 33 33 33 100 100 100 1500 ...
## $ thetai : num 0.413 0.402 0.414 0.4 0.402 ...
## $ trt : Factor w/ 4 levels "KH·Interrow",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ logtensao: num 0.778 0.778 0.778 1.519 1.519 ...
## $ rept : int 1 2 3 1 2 3 1 2 3 1 ...
## $ ue1 : Factor w/ 6 levels "KH.1","GA.1",..: 2 4 6 2 4 6 2 4 6 2 ...
## $ ue2 : Factor w/ 12 levels "KH.1.Interrow",..: 2 4 6 2 4 6 2 4 6 2 ...# Cell frequency.
ftable(xtabs(~solo + posi + tensao, data = da))## tensao 6 33 100 1500 3000
## solo posi
## KH Interrow 3 3 3 3 3
## Row 3 3 3 3 3
## GA Interrow 3 3 3 3 3
## Row 3 3 3 3 3# Differnce in position mean profile conditional to other factors.
ggplot(data = da,
mapping = aes(x = tensao, y = pr, color = posi)) +
facet_grid(facets = ~solo) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
scale_x_log10() +
labs(x = axis_text$tensao,
y = axis_text$pr,
color = axis_text$posi)
# Differnce in position mean profile conditional to other factors.
ggplot(data = filter(da, tensao < 3000),
mapping = aes(x = tensao, y = pr, color = posi)) +
facet_grid(facets = ~solo) +
geom_point() +
geom_smooth(method = "lm",
formula = y ~ poly(x, degree = 2),
se = FALSE) +
scale_x_log10() +
labs(x = axis_text$tensao,
y = axis_text$pr,
color = axis_text$posi)
# Differnce in soil mean profile conditional to other factors.
ggplot(data = da,
mapping = aes(x = tensao, y = pr, color = solo)) +
facet_grid(facets = ~posi) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
scale_x_log10() +
labs(x = axis_text$tensao,
y = axis_text$pr,
color = axis_text$solo)
Linear mixed effects model were used to model soil resistance to penetration as function of experimental factors: soil, position and matric tension. As suggested by exploratory data analysis, the effect of matric tension (at log 10 scale) on penetration resistance can be properly described by a second order degree polynomial in the range of tension from 6 to 1500 kPa.
The statistical model is \[ \begin{align*} y_{ijkl} &= \mu + S_i + P_j + \beta_1 T_k + \beta_2 T_k^2 + \\ &\qquad (SP)_{ij} + (\beta_{1i} T_k + \beta_{2i} T_k^2) + (\beta_{1j} T_k + \beta_{2j} T_k^2) + \\ &\qquad (\beta_{1ij} T_k + \beta_{2ij} T_k^2) + \\ &\qquad e_{il} + e_{ijl} + e_{ijkl} \\ e_{il} &\sim \text{Normal}(0, \sigma^2_a) \qquad [\text{6 whole plots}]\\ e_{ijl} &\sim \text{Normal}(0, \sigma^2_b) \qquad [\text{12 sub plots}]\\ e_{ijkl} &\sim \text{Normal}(0, \sigma^2) \qquad [\text{48 data points}]\\ \end{align*} \]
- \(y_{ijkl}\) is response observed at the experimental condition \(ijkl\).
- \(\mu\) is a constant that belongs to all observations.
- \(S_i\) is the effect of soil \(i = 1, 2\).
- \(P_j\) is the effect of position \(j = 1, 2\).
- \(\beta_1\) and \(\beta_2\) are the first and second order effect of log 10 matric tension.
- The following terms are two and three way interactions among soil, position and matric tension.
- \(e_{il}\) is the experimental error at experimental unit \(il\), \(e_{ijl}\) is the experimental error at experimental unit \(ijl\) and \(e_{ijkl}\) is the experimental error at sample unit \(ijkl\).
This model were estimated by restricted maximum likelihood (REML).
# Random intercept, slope and curvature for random effects.
# Interaction until the 4th degree.
m0 <- lmer(pr ~ solo * posi * (logtensao + I(logtensao^2)) +
(1 | ue1) + (1 | ue2),
data = subset(da, tensao < 3000))
# Variance of the random effects.
VarCorr(m0)## Groups Name Std.Dev.
## ue2 (Intercept) 0.046621
## ue1 (Intercept) 0.027758
## Residual 0.223862# Test for the fixed effects terms.
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## solo 0.15192 0.15192 1 30.144 3.0316
## posi 0.07303 0.07303 1 29.975 1.4573
## logtensao 0.00788 0.00788 1 27.999 0.1572
## I(logtensao^2) 0.24008 0.24008 1 27.999 4.7906
## solo:posi 0.02487 0.02487 1 29.975 0.4963
## solo:logtensao 0.34764 0.34764 1 27.999 6.9370
## solo:I(logtensao^2) 0.56398 0.56398 1 27.999 11.2538
## posi:logtensao 0.19194 0.19194 1 27.999 3.8300
## posi:I(logtensao^2) 0.21903 0.21903 1 27.999 4.3707
## solo:posi:logtensao 0.01804 0.01804 1 27.999 0.3599
## solo:posi:I(logtensao^2) 0.03522 0.03522 1 27.999 0.7029
## Pr(>F)
## solo 0.091858 .
## posi 0.236793
## logtensao 0.694734
## I(logtensao^2) 0.037117 *
## solo:posi 0.486574
## solo:logtensao 0.013596 *
## solo:I(logtensao^2) 0.002295 **
## posi:logtensao 0.060380 .
## posi:I(logtensao^2) 0.045761 *
## solo:posi:logtensao 0.553372
## solo:posi:I(logtensao^2) 0.408924
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model summary.
# summary(m0)
# Fitting a smaller model with relevant terms.
m1 <- update(m0,
. ~ (posi + solo) *
(logtensao + I(logtensao^2)) +
# posi:solo +
(1 | ue1) + (1 | ue2))## boundary (singular) fit: see ?isSingularanova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: subset(da, tensao < 3000)
## Models:
## m1: pr ~ posi + solo + logtensao + I(logtensao^2) + (1 | ue1) + (1 |
## m1: ue2) + posi:logtensao + posi:I(logtensao^2) + solo:logtensao +
## m1: solo:I(logtensao^2)
## m0: pr ~ solo * posi * (logtensao + I(logtensao^2)) + (1 | ue1) +
## m0: (1 | ue2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 12 12.972 35.426 5.5141 -11.028
## m0 15 11.137 39.205 9.4313 -18.863 7.8344 3 0.04956 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Test for the fixed effects terms.
anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## posi 0.07182 0.07182 1 33.217 1.4263 0.240823
## solo 0.15002 0.15002 1 33.217 2.9792 0.093633
## logtensao 0.00788 0.00788 1 30.000 0.1565 0.695229
## I(logtensao^2) 0.24008 0.24008 1 30.000 4.7677 0.036955
## posi:logtensao 0.19194 0.19194 1 30.000 3.8117 0.060284
## posi:I(logtensao^2) 0.21903 0.21903 1 30.000 4.3497 0.045621
## solo:logtensao 0.34764 0.34764 1 30.000 6.9037 0.013423
## solo:I(logtensao^2) 0.56398 0.56398 1 30.000 11.1999 0.002213
##
## posi
## solo .
## logtensao
## I(logtensao^2) *
## posi:logtensao .
## posi:I(logtensao^2) *
## solo:logtensao *
## solo:I(logtensao^2) **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1All three way interactions were dropped beacause were not relevant to the model. Also, the two way interaction between soil and position were not relevant and, then, dropped. In summary, soil and position had only additive effect, but they had interaction with first and second order polynomial terms for tension effect.
# Estimated marginal means with covariate set at 1. The corresponding
# matrix will be used to get que equation of each experimental
# condition.
emm <- emmeans(m1,
specs = ~solo + posi,
at = list(logtensao = 1))
# Get experimental points and matrix to calculate marginal means.
grid <- attr(emm, "grid")
L <- attr(emm, "linfct")
# Fixed effects estimates.
b <- fixef(m1)
# Zero order terms (intercept).
ord0 <- names(b) %>%
grepl(pattern = "logtensao") %>% `!`()
# First order terms (slope).
ord1 <- names(b) %>%
grepl(pattern = "logtensao$")
# Second order terms (curvatures).
ord2 <- names(b) %>%
grepl(pattern = "I\\(logtensao\\^2\\)$")
# Checking: all row sum must be 1.
# rowSums(cbind(ord0, ord1, ord2))
# Get the equation coefficients: b0 + b1 * x + b2 * x^2.
grid$b0 <- c(L %*% (b * ord0))
grid$b1 <- c(L %*% (b * ord1))
grid$b2 <- c(L %*% (b * ord2))
# Equations.
grid <- grid %>%
mutate(label = sprintf("'%0.2f' %s '%0.3f' * x %s '%0.4f' * x^2",
b0,
ifelse(b1 >= 0, "+", "-"),
abs(b1),
ifelse(b2 >= 0, "+", "-"),
abs(b2)))
# Table of equation coefficients for each experimental point.
grid %>%
select(solo, posi, label)## solo posi label
## 1 KH Interrow '0.44' - '0.214' * x + '0.1674' * x^2
## 2 GA Interrow '-0.18' + '0.835' * x - '0.1578' * x^2
## 3 KH Row '0.87' - '0.992' * x + '0.3700' * x^2
## 4 GA Row '0.25' + '0.056' * x + '0.0448' * x^2# Checking: are the equation curves correct? Compare with `predict()`.
grid_curves <- grid %>%
group_by(posi, solo) %>%
do({
logtensao <- seq(0.5, 3.5, length = 31)
fit <- .$b0 + .$b1 * logtensao + .$b2 * logtensao^2
data.frame(logtensao, fit)
})
# grid_curves# Conditions to get estimated values.
pred <- with(da,
expand.grid(solo = levels(solo),
posi = levels(posi),
logtensao = seq(0.5, 3.5, length = 31),
KEEP.OUT.ATTRS = FALSE))
# Model matrix for prediction.
X <- model.matrix(formula(terms(m1))[-2], data = pred)
# Checking: length and names.
stopifnot(ncol(X) == length(fixef(m1)))
stopifnot(all(colnames(X) == names(fixef(m1))))
# Predicted value.
pred$fit <- X %*% fixef(m1)
# Standard error of predicted values.
pred$se <- sqrt(diag(X %*% tcrossprod(vcov(m1), X)))
# Quantile of t distribution.
qtl <- qt(p = 0.975, df = df.residual(m1))
# Confidence limits for predicted value.
pred$lwr <- pred$fit - qtl * pred$se
pred$upr <- pred$fit + qtl * pred$se
# Adds the equation to annotate on graphics.
grid <- grid %>%
mutate(xpos = 2.5,
ypos = ifelse(solo == "GA", 2.10, 2.02) - 0.1,
hjust = c(0),
vjust = c(0))
# Fitted curves with confidence bands and observed data.
# Soil in evidence.
ggplot(data = filter(da, tensao < 3000),
mapping = aes(x = tensao, y = pr, color = solo)) +
facet_grid(facets = ~posi) +
geom_point() +
scale_x_log10() +
# geom_line(data = pred,
# mapping = aes(x = 10^logtensao,
# y = fit,
# color = solo)) +
geom_smooth(data = pred,
mapping = aes(x = 10^logtensao,
y = fit,
ymin = lwr,
ymax = upr),
stat = "identity",
show.legend = FALSE) +
geom_text(data = grid,
mapping = aes(x = xpos,
y = ypos,
label = label,
hjust = hjust,
vjust = vjust),
show.legend = FALSE,
parse = TRUE) +
# geom_line(data = grid_curves,
# mapping = aes(x = 10^logtensao,
# y = fit,
# group = solo),
# linetype = 2,
# color = "black") +
labs(x = axis_text$tensao,
y = axis_text$pr,
color = axis_text$solo)
grid <- grid %>%
mutate(xpos = 2.5,
ypos = ifelse(posi == "Interrow", 2.10, 1.88))
# Fitted curves with confidence bands and observed data.
# Position in evidence.
gg1 <-
ggplot(data = filter(da, tensao < 3000),
mapping = aes(x = tensao, y = pr, color = posi)) +
facet_grid(facets = ~solo) +
geom_point() +
scale_x_log10() +
geom_line(data = pred,
mapping = aes(x = 10^logtensao,
y = fit,
color = posi)) +
geom_smooth(data = pred,
mapping = aes(x = 10^logtensao,
y = fit,
ymin = lwr,
ymax = upr),
stat = "identity",
show.legend = FALSE) +
geom_text(data = grid,
mapping = aes(x = xpos,
y = ypos,
label = label,
hjust = hjust,
vjust = vjust),
show.legend = FALSE,
parse = TRUE) +
# geom_line(data = grid_curves,
# mapping = aes(x = 10^logtensao,
# y = fit,
# group = posi),
# linetype = 2,
# color = "black") +
labs(x = axis_text$tensao,
y = axis_text$pr,
color = axis_text$posi)
gg1
# Tests of soil curves are equal in each posi:prof combination.
grid %>%
mutate(row = 1:n()) %>%
group_by(posi) %>%
do({
compare_curves(.$row, .$solo, L, list(ord0, ord1, ord2))
}) %>%
as.data.frame()## posi contrast F Pr(>F)
## 1 Interrow KH-GA 7.415602 0.003546246
## 2 Row KH-GA 7.415602 0.003546246# Tests of posi curves are equal in each soil:prof combination.
grid %>%
mutate(row = 1:n()) %>%
group_by(solo) %>%
do({
compare_curves(.$row, .$posi, L, list(ord0, ord1, ord2))
}) %>%
as.data.frame()## solo contrast F Pr(>F)
## 1 KH Interrow-Row 2.669739 0.08250887
## 2 GA Interrow-Row 2.669739 0.08250887According to Wald F test to compare curve parameters, the null hypothesis of curve equality between soils in each soil were rejected (p < 0.01). The result of the hypothesis test were the same for both position because soil and position did not interact. So, the difference between curves due to soil don’t depend of sampling position.
According to Wald F test to compare curve parameters, there isn’t enought evidence at 5% to reject the null hypothesis of curve equality between sampling position (p < 0.10). The result of the hypothesis test were the same for both soils because soil and position did not interact. So, the difference between curves due to sampling position don’t depend of soil.
Water content
# Differnce in position mean profile conditional to other factors.
ggplot(data = da,
mapping = aes(x = tensao, y = thetai, color = posi)) +
facet_grid(facets = ~solo) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
scale_x_log10() +
labs(x = axis_text$tensao,
y = axis_text$thetai,
color = axis_text$posi)
# Differnce in position mean profile conditional to other factors.
ggplot(data = filter(da, tensao < 3000),
mapping = aes(x = tensao, y = thetai, color = posi)) +
facet_grid(facets = ~solo) +
geom_point() +
geom_smooth(method = "lm",
formula = y ~ poly(x, degree = 2),
se = FALSE) +
scale_x_log10() +
labs(x = axis_text$tensao,
y = axis_text$thetai,
color = axis_text$posi)
# Differnce in soil mean profile conditional to other factors.
ggplot(data = da,
mapping = aes(x = tensao, y = thetai, color = solo)) +
facet_grid(facets = ~posi) +
geom_point() +
stat_summary(geom = "line", fun.y = "mean") +
scale_x_log10() +
labs(x = axis_text$tensao,
y = axis_text$thetai,
color = axis_text$solo)
Linear mixed effects model were used to model water content as function of experimental factors: soil, position and matric tension. As suggested by exploratory data analysis, the effect of matric tension (at log 10 scale) on soil water content can be properly described by a second order degree polynomial in the range of tension from 6 to 1500 kPa.
The statistical model is the same applied to resistance to penetration, \[ \begin{align*} y_{ijkl} &= \mu + S_i + P_j + \beta_1 T_k + \beta_2 T_k^2 + \\ &\qquad (SP)_{ij} + (\beta_{1i} T_k + \beta_{2i} T_k^2) + (\beta_{1j} T_k + \beta_{2j} T_k^2) + \\ &\qquad (\beta_{1ij} T_k + \beta_{2ij} T_k^2) + \\ &\qquad e_{il} + e_{ijl} + e_{ijkl} \\ e_{il} &\sim \text{Normal}(0, \sigma^2_a) \qquad [\text{6 whole plots}]\\ e_{ijl} &\sim \text{Normal}(0, \sigma^2_b) \qquad [\text{12 sub plots}]\\ e_{ijkl} &\sim \text{Normal}(0, \sigma^2) \qquad [\text{48 data points}]\\ \end{align*} \]
- \(y_{ijkl}\) is response observed at the experimental condition \(ijkl\).
- \(\mu\) is a constant that belongs to all observations.
- \(S_i\) is the effect of soil \(i = 1, 2\).
- \(P_j\) is the effect of position \(j = 1, 2\).
- \(\beta_1\) and \(\beta_2\) are the first and second order effect of log 10 matric tension.
- The following terms are two and three way interactions among soil, position and matric tension.
- \(e_{il}\) is the experimental error at experimental unit \(il\), \(e_{ijl}\) is the experimental error at experimental unit \(ijl\) and \(e_{ijkl}\) is the experimental error at sample unit \(ijkl\).
This model were estimated by restricted maximum likelihood (REML).
# Random intercept, slope and curvature for random effects.
# Interaction until the 4th degree.
m0 <- lmer(thetai ~ solo * posi * (logtensao + I(logtensao^2)) +
(1 | ue1) + (1 | ue2),
data = subset(da, tensao < 3000))
# Variance of the random effects.
VarCorr(m0)## Groups Name Std.Dev.
## ue2 (Intercept) 0.0079023
## ue1 (Intercept) 0.0055402
## Residual 0.0166652# Test for the fixed effects terms.
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## solo 0.0035795 0.0035795 1 31.603 12.8884
## posi 0.0002328 0.0002328 1 30.870 0.8382
## logtensao 0.0062800 0.0062800 1 28.000 22.6121
## I(logtensao^2) 0.0005569 0.0005569 1 28.000 2.0053
## solo:posi 0.0002694 0.0002694 1 30.870 0.9700
## solo:logtensao 0.0010077 0.0010077 1 28.000 3.6283
## solo:I(logtensao^2) 0.0002967 0.0002967 1 28.000 1.0685
## posi:logtensao 0.0000280 0.0000280 1 28.000 0.1009
## posi:I(logtensao^2) 0.0000818 0.0000818 1 28.000 0.2945
## solo:posi:logtensao 0.0001672 0.0001672 1 28.000 0.6021
## solo:posi:I(logtensao^2) 0.0000566 0.0000566 1 28.000 0.2037
## Pr(>F)
## solo 0.001104 **
## posi 0.367014
## logtensao 5.415e-05 ***
## I(logtensao^2) 0.167784
## solo:posi 0.332345
## solo:logtensao 0.067125 .
## solo:I(logtensao^2) 0.310138
## posi:logtensao 0.753121
## posi:I(logtensao^2) 0.591671
## solo:posi:logtensao 0.444287
## solo:posi:I(logtensao^2) 0.655226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Model summary.
# summary(m0)
# Fitting a smaller model with relevant terms.
m1 <- update(m0, . ~ posi + solo * logtensao +
(1 | ue1) + (1 | ue2))
anova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: subset(da, tensao < 3000)
## Models:
## m1: thetai ~ posi + solo + logtensao + (1 | ue1) + (1 | ue2) + solo:logtensao
## m0: thetai ~ solo * posi * (logtensao + I(logtensao^2)) + (1 | ue1) +
## m0: (1 | ue2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 8 -235.36 -220.39 125.68 -251.36
## m0 15 -230.77 -202.70 130.38 -260.77 9.4064 7 0.2248# Test for the fixed effects terms.
anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## posi 0.005157 0.005157 1 4.999 17.496 0.008634 **
## solo 0.008781 0.008781 1 22.249 29.787 1.687e-05 ***
## logtensao 0.090524 0.090524 1 34.009 307.089 < 2.2e-16 ***
## solo:logtensao 0.006326 0.006326 1 34.009 21.459 5.119e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1According F test, all interactions were not relevant except soil with the linear effect of tension. So, the final model has all main effect of soil, position and (linear) tension, and also the soil and (linear) tension interaction.
# Estimated marginal means with covariate set at 1. The corresponding
# matrix will be used to get que equation of each experimental
# condition.
emm <- emmeans(m1,
specs = ~solo + posi,
at = list(logtensao = 1))
# Get experimental points and matrix to calculate marginal means.
grid <- attr(emm, "grid")
L <- attr(emm, "linfct")
# Fixed effects estimates.
b <- fixef(m1)
# Zero order terms (intercept).
ord0 <- names(b) %>%
# grep(pattern = "logtensao", value = TRUE, invert = TRUE)
grepl(pattern = "logtensao") %>% `!`()
# First order terms (slope).
ord1 <- names(b) %>%
# grep(pattern = "logtensao$", value = TRUE, invert = FALSE)
grepl(pattern = "logtensao$")
# Checking: all row sum must be 1.
# rowSums(cbind(ord0, ord1))
# Get the equation coefficients: b0 + b1 * x.
grid$b0 <- c(L %*% (b * ord0))
grid$b1 <- c(L %*% (b * ord1))
grid <- grid %>%
mutate(label = sprintf("'%0.2f' %s '%0.3f' * x",
b0,
ifelse(b1 >= 0, "+", "-"),
abs(b1)))
# Table of equation coefficients for each experimental point.
grid %>%
select(solo, posi, label)## solo posi label
## 1 KH Interrow '0.52' - '0.063' * x
## 2 GA Interrow '0.44' - '0.037' * x
## 3 KH Row '0.49' - '0.063' * x
## 4 GA Row '0.42' - '0.037' * x# Checking: are the equation curves correct? Compare with `predict()`.
grid_curves <- grid %>%
group_by(posi, solo) %>%
do({
logtensao <- seq(0.5, 3.5, length = 31)
fit <- .$b0 + .$b1 * logtensao
data.frame(logtensao, fit)
})
# grid_curves# Conditions to get estimated values.
pred <- with(da,
expand.grid(solo = levels(solo),
posi = levels(posi),
logtensao = seq(0.5, 3.5, length = 31),
KEEP.OUT.ATTRS = FALSE))
# Model matrix for prediction.
X <- model.matrix(formula(terms(m1))[-2], data = pred)
# Checking: length and names.
stopifnot(ncol(X) == length(fixef(m1)))
stopifnot(all(colnames(X) == names(fixef(m1))))
# Predicted value.
pred$fit <- X %*% fixef(m1)
# Standard error of predicted values.
pred$se <- sqrt(diag(X %*% tcrossprod(vcov(m1), X)))
# Quantile of t distribution.
qtl <- qt(p = 0.975, df = df.residual(m1))
# Confidence limits for predicted value.
pred$lwr <- pred$fit - qtl * pred$se
pred$upr <- pred$fit + qtl * pred$se
# Adds the equation to annotate on graphics.
grid <- grid %>%
mutate(xpos = 2.5,
ypos = ifelse(solo == "GA", 0.29, 0.30) + 0,
hjust = c(0),
vjust = c(0))
# Fitted curves with confidence bands and observed data.
# Soil in evidence.
ggplot(data = filter(da, tensao < 3000),
mapping = aes(x = tensao, y = thetai, color = solo)) +
facet_grid(facets = ~posi) +
geom_point() +
scale_x_log10() +
# geom_line(data = pred,
# mapping = aes(x = 10^logtensao,
# y = fit,
# color = solo)) +
geom_smooth(data = pred,
mapping = aes(x = 10^logtensao,
y = fit,
ymin = lwr,
ymax = upr),
stat = "identity",
show.legend = FALSE) +
geom_text(data = grid,
mapping = aes(x = xpos,
y = ypos,
label = label,
hjust = hjust,
vjust = vjust),
show.legend = FALSE,
parse = TRUE) +
# geom_line(data = grid_curves,
# mapping = aes(x = 10^logtensao,
# y = fit,
# group = solo),
# linetype = 2,
# color = "black") +
labs(x = axis_text$tensao,
y = axis_text$thetai,
color = axis_text$solo)
grid <- grid %>%
mutate(xpos = 2.5,
ypos = ifelse(posi == "Interrow", 0.28, 0.30))
# Fitted curves with confidence bands and observed data.
# Position in evidence.
gg2 <-
ggplot(data = filter(da, tensao < 3000),
mapping = aes(x = tensao, y = thetai, color = posi)) +
facet_grid(facets = ~solo) +
geom_point() +
scale_x_log10() +
geom_line(data = pred,
mapping = aes(x = 10^logtensao,
y = fit,
color = posi)) +
geom_smooth(data = pred,
mapping = aes(x = 10^logtensao,
y = fit,
ymin = lwr,
ymax = upr),
stat = "identity",
show.legend = FALSE) +
geom_text(data = grid,
mapping = aes(x = xpos,
y = ypos,
label = label,
hjust = hjust,
vjust = vjust),
show.legend = FALSE,
parse = TRUE) +
# geom_line(data = grid_curves,
# mapping = aes(x = 10^logtensao,
# y = fit,
# group = posi),
# linetype = 2,
# color = "black") +
labs(x = axis_text$tensao,
y = axis_text$thetai,
color = axis_text$posi)
gg2
# Tests of soil curves are equal in each posi:prof combination.
grid %>%
mutate(row = 1:n()) %>%
group_by(posi) %>%
do({
compare_curves(.$row, .$solo, L, list(ord0, ord1))
}) %>%
as.data.frame()## posi contrast F Pr(>F)
## 1 Interrow KH-GA 13.96336 0.001828974
## 2 Row KH-GA 13.96336 0.001828974# Tests of posi curves are equal in each soil:prof combination.
grid %>%
mutate(row = 1:n()) %>%
group_by(solo) %>%
do({
compare_curves(.$row, .$posi, L, list(ord0, ord1))
}) %>%
as.data.frame()## solo contrast F Pr(>F)
## 1 KH Interrow-Row 17.49569 0.008631706
## 2 GA Interrow-Row 17.49569 0.008631706By Wald F test, the null hypothesis of equality between soil curves at each position were rejected (p < 0.01). The test statistic is the same for each position because soil and position do not interact. The null hypothesis of equality between position curves at each soil were also rejected (p < 0.01). The test statistic is the same for each soil because soil and position do not interact.
Combining the results
gridExtra::grid.arrange(gg1, gg2, ncol = 1)
Porosity and bulk density
# Data structure.
str(da)## Classes 'tbl_df', 'tbl' and 'data.frame': 36 obs. of 11 variables:
## $ solo : Factor w/ 2 levels "KH","GA": 1 1 1 1 1 1 1 1 1 1 ...
## $ posi : Factor w/ 2 levels "Interrow","Row": 1 1 1 1 1 1 1 1 1 2 ...
## $ prof : Factor w/ 3 levels "0-5","10-15",..: 1 1 1 2 2 2 3 3 3 1 ...
## $ rept : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
## $ total: num 0.665 0.65 0.615 0.653 0.605 ...
## $ macro: num 0.242 0.26 0.16 0.287 0.205 ...
## $ micro: num 0.42 0.39 0.46 0.37 0.4 0.47 0.39 0.42 0.46 0.37 ...
## $ ds : num 0.83 0.953 0.912 0.868 0.925 ...
## $ trt : Factor w/ 12 levels "KH·Interrow·0-5",..: 1 1 1 5 5 5 9 9 9 3 ...
## $ ue1 : Factor w/ 6 levels "KH.1","GA.1",..: 1 3 5 1 3 5 1 3 5 1 ...
## $ ue2 : Factor w/ 12 levels "KH.1.Interrow",..: 1 3 5 1 3 5 1 3 5 7 ...# Frequency of cell combination.
ftable(xtabs(~solo + posi + prof, data = da))## prof 0-5 10-15 20-25
## solo posi
## KH Interrow 3 3 3
## Row 3 3 3
## GA Interrow 3 3 3
## Row 3 3 3db <- da %>%
gather(key = "variable", value = "value", total, macro, micro, ds)
ggplot(data = db,
mapping = aes(x = prof, y = value, color = posi, group = posi)) +
facet_grid(facets = variable ~ solo, scale = "free_y") +
geom_jitter(width = 0.1, height = 0) +
stat_summary(geom = "line", fun.y = "mean")
# To keep results.
tb_results <- list()Density and porosity
Soil density
Linear mixed effects model were used to model the responses:
- Soil density,
- Total porosity,
- Macroporosity and
- Microporosity,
as function of experimental factors: soil, position and depth.
The statistical model assumes a split split plot experimental design, \[ \begin{align*} y_{ijkl} &= \mu + S_i + P_j + D_k + \\ &\qquad (SP)_{ij} + (SD)_{ik} + (PD)_{jk} + (SPD)_{ijk} + \\ &\qquad e_{il} + e_{ijl} + e_{ijkl} \\ e_{il} &\sim \text{Normal}(0, \sigma^2_a) \qquad [\text{6 whole plots}]\\ e_{ijl} &\sim \text{Normal}(0, \sigma^2_b) \qquad [\text{12 sub plots}]\\ e_{ijkl} &\sim \text{Normal}(0, \sigma^2) \qquad [\text{36 data points}]\\ \end{align*} \]
- \(y_{ijkl}\) is response observed at the experimental condition \(ijkl\).
- \(\mu\) is a constant that belongs to all observations.
- \(S_i\) is the effect of soil \(i = 1, 2\).
- \(P_j\) is the effect of position \(j = 1, 2\).
- \(D_k\) is the effect of depth \(k = 1, 2, 3\).
- The following terms are two and three way interactions among soil, position and depth.
- \(e_{il}\) is the experimental error at experimental unit \(il\), \(e_{ijl}\) is the experimental error at experimental unit \(ijl\) and \(e_{ijkl}\) is the experimental error at sample unit \(ijkl\).
This model were estimated by restricted maximum likelihood (REML).
m0 <- lmer(ds ~ solo * posi * prof + (1 | ue1) + (1 | ue2),
data = da)## boundary (singular) fit: see ?isSingularVarCorr(m0)## Groups Name Std.Dev.
## ue2 (Intercept) 7.5996e-07
## ue1 (Intercept) 1.4535e-02
## Residual 4.0264e-02# summary(m0)
anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.033802 0.033802 1 4 20.8497 0.01029 *
## posi 0.003959 0.003959 1 20 2.4417 0.13383
## prof 0.007211 0.003605 2 20 2.2239 0.13425
## solo:posi 0.000127 0.000127 1 20 0.0781 0.78280
## solo:prof 0.017423 0.008712 2 20 5.3735 0.01356 *
## posi:prof 0.009902 0.004951 2 20 3.0540 0.06959 .
## solo:posi:prof 0.002376 0.001188 2 20 0.7328 0.49303
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- update(m0,
. ~ (solo + posi + prof)^2 - solo:posi +
(1 | ue1) + (1 | ue2),
data = da)## boundary (singular) fit: see ?isSingularanova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: da
## Models:
## m1: ds ~ solo + posi + prof + (1 | ue1) + (1 | ue2) + solo:prof +
## m1: posi:prof
## m0: ds ~ solo * posi * prof + (1 | ue1) + (1 | ue2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 12 -114.02 -95.020 69.011 -138.02
## m0 15 -110.25 -86.499 70.126 -140.25 2.2305 3 0.526anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.031662 0.031662 1 4 20.8499 0.010291 *
## posi 0.003959 0.003959 1 23 2.6067 0.120047
## prof 0.007211 0.003605 2 23 2.3742 0.115513
## solo:prof 0.017423 0.008712 2 23 5.7368 0.009523 **
## posi:prof 0.009902 0.004951 2 23 3.2605 0.056673 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm <- emmeans(m1, ~prof | solo)
as.data.frame(emm)## prof solo emmean SE df lower.CL upper.CL
## 1 0-5 KH 0.8975000 0.01814423 13.51316 0.8584525 0.9365475
## 2 10-15 KH 0.8691667 0.01814423 13.51316 0.8301192 0.9082141
## 3 20-25 KH 0.9004167 0.01814423 13.51316 0.8613692 0.9394641
## 4 0-5 GA 0.9233333 0.01814423 13.51316 0.8842859 0.9623808
## 5 10-15 GA 1.0025000 0.01814423 13.51316 0.9634525 1.0415475
## 6 20-25 GA 0.9866667 0.01814423 13.51316 0.9476192 1.0257141results <- multcomp::cld(emm)
results## solo = KH:
## prof emmean SE df lower.CL upper.CL .group
## 10-15 0.869 0.0181 13.5 0.830 0.908 1
## 0-5 0.897 0.0181 13.5 0.858 0.937 1
## 20-25 0.900 0.0181 13.5 0.861 0.939 1
##
## solo = GA:
## prof emmean SE df lower.CL upper.CL .group
## 0-5 0.923 0.0181 13.5 0.884 0.962 1
## 20-25 0.987 0.0181 13.5 0.948 1.026 2
## 10-15 1.002 0.0181 13.5 0.963 1.042 2
##
## Results are averaged over the levels of: posi
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 3 estimates
## significance level used: alpha = 0.05results <- results %>%
as.data.frame() %>%
mutate(cld = c("a", "a", "a", "b", "a", "a"))
gg1 <-
ggplot(data = results,
mapping = aes(x = prof,
y = emmean,
ymin = lower.CL,
ymax = upper.CL)) +
facet_wrap(facets = ~solo) +
geom_point() +
geom_errorbar(width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.2f %s",
emmean,
cld)),
angle = 90, vjust = -0.5) +
labs(x = axis_text$prof,
y = axis_text$ds)
tb_results$ds <- resultsTotal porosity
m0 <- lmer(total ~ solo * posi * prof + (1 | ue1) + (1 | ue2),
data = da)
VarCorr(m0)## Groups Name Std.Dev.
## ue2 (Intercept) 0.0065953
## ue1 (Intercept) 0.0061089
## Residual 0.0164226anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.0034358 0.0034358 1 4.0001 12.7392 0.02339 *
## posi 0.0052585 0.0052585 1 3.9997 19.4974 0.01155 *
## prof 0.0003962 0.0001981 2 16.0002 0.7345 0.49525
## solo:posi 0.0009903 0.0009903 1 3.9997 3.6718 0.12784
## solo:prof 0.0007858 0.0003929 2 16.0002 1.4567 0.26230
## posi:prof 0.0005108 0.0002554 2 16.0002 0.9469 0.40864
## solo:posi:prof 0.0002545 0.0001273 2 16.0002 0.4718 0.63226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- update(m0,
. ~ solo + posi +
(1 | ue1) + (1 | ue2),
data = da)
anova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: da
## Models:
## m1: total ~ solo + posi + (1 | ue1) + (1 | ue2)
## m0: total ~ solo * posi * prof + (1 | ue1) + (1 | ue2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 6 -175.98 -166.47 93.988 -187.98
## m0 15 -170.89 -147.13 100.443 -200.89 12.91 9 0.1667anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.0033245 0.0033245 1 3.9846 12.741 0.02354 *
## posi 0.0033158 0.0033158 1 4.9834 12.707 0.01623 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm <- emmeans(m1, ~posi | solo)
as.data.frame(emm)## posi solo emmean SE df lower.CL upper.CL
## 1 Interrow KH 0.6333333 0.007192112 7.459632 0.6165368 0.6501298
## 2 Row KH 0.6627778 0.007192112 7.459632 0.6459813 0.6795743
## 3 Interrow GA 0.6036111 0.007192112 7.459632 0.5868146 0.6204076
## 4 Row GA 0.6330556 0.007192112 7.459632 0.6162590 0.6498521results <- multcomp::cld(emm)
results## solo = KH:
## posi emmean SE df lower.CL upper.CL .group
## Interrow 0.633 0.00719 7.46 0.617 0.65 1
## Row 0.663 0.00719 7.46 0.646 0.68 2
##
## solo = GA:
## posi emmean SE df lower.CL upper.CL .group
## Interrow 0.604 0.00719 7.46 0.587 0.62 1
## Row 0.633 0.00719 7.46 0.616 0.65 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05results <- results %>%
as.data.frame() %>%
mutate(cld = c("b", "a", "b", "a"))
gg2 <-
ggplot(data = results,
mapping = aes(x = posi,
y = emmean,
ymin = lower.CL,
ymax = upper.CL)) +
facet_wrap(facets = ~solo) +
geom_point() +
geom_errorbar(width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.2f %s",
emmean,
cld)),
angle = 90, vjust = -0.5) +
labs(x = axis_text$posi,
y = axis_text$total)
tb_results$total <- resultsMacro porosity
m0 <- lmer(macro ~ solo * posi * prof + (1 | ue1) + (1 | ue2),
data = da)
VarCorr(m0)## Groups Name Std.Dev.
## ue2 (Intercept) 0.018951
## ue1 (Intercept) 0.017649
## Residual 0.030146anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.0009331 0.0009331 1 3.9996 1.0268 0.36823
## posi 0.0150311 0.0150311 1 3.9993 16.5399 0.01527 *
## prof 0.0012931 0.0006465 2 16.0011 0.7114 0.50583
## solo:posi 0.0018112 0.0018112 1 3.9993 1.9930 0.23089
## solo:prof 0.0048389 0.0024194 2 16.0011 2.6623 0.10044
## posi:prof 0.0001500 0.0000750 2 16.0011 0.0825 0.92117
## solo:posi:prof 0.0002931 0.0001465 2 16.0011 0.1612 0.85246
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- update(m0,
. ~ posi +
(1 | ue1) + (1 | ue2),
data = da)
anova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: da
## Models:
## m1: macro ~ posi + (1 | ue1) + (1 | ue2)
## m0: macro ~ solo * posi * prof + (1 | ue1) + (1 | ue2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 5 -128.44 -120.526 69.222 -138.44
## m0 15 -121.19 -97.441 75.597 -151.19 12.75 10 0.238anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## posi 0.012143 0.012143 1 5 13.801 0.01378 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm <- emmeans(m1, ~posi)
as.data.frame(emm)## posi emmean SE df lower.CL upper.CL
## 1 Interrow 0.2018056 0.01318235 9.45958 0.1722044 0.2314067
## 2 Row 0.2622222 0.01318235 9.45958 0.2326211 0.2918234results <- multcomp::cld(emm)
results## posi emmean SE df lower.CL upper.CL .group
## Interrow 0.202 0.0132 9.46 0.172 0.231 1
## Row 0.262 0.0132 9.46 0.233 0.292 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05results <- results %>%
as.data.frame() %>%
mutate(cld = c("b", "a"))
gg3 <-
ggplot(data = results,
mapping = aes(x = posi,
y = emmean,
ymin = lower.CL,
ymax = upper.CL)) +
geom_point() +
geom_errorbar(width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.2f %s",
emmean,
cld)),
angle = 90, vjust = -0.5) +
labs(x = axis_text$posi,
y = axis_text$macro)
tb_results$macro <- resultsMicro porosity
m0 <- lmer(micro ~ solo * posi * prof + (1 | ue1) + (1 | ue2),
data = da)
VarCorr(m0)## Groups Name Std.Dev.
## ue2 (Intercept) 0.0179609
## ue1 (Intercept) 0.0090285
## Residual 0.0209170anova(m0)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## solo 0.00031051 0.00031051 1 4.0000 0.7097 0.44697
## posi 0.00252182 0.00252182 1 4.0007 5.7639 0.07428
## prof 0.00057222 0.00028611 2 15.9994 0.6539 0.53335
## solo:posi 0.00005535 0.00005535 1 4.0007 0.1265 0.74005
## solo:prof 0.00160556 0.00080278 2 15.9994 1.8348 0.19169
## posi:prof 0.00015000 0.00007500 2 15.9994 0.1714 0.84399
## solo:posi:prof 0.00100556 0.00050278 2 15.9994 1.1491 0.34172
##
## solo
## posi .
## prof
## solo:posi
## solo:prof
## posi:prof
## solo:posi:prof
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- update(m0,
. ~ posi +
(1 | ue1) + (1 | ue2),
data = da)
anova(m0, m1)## refitting model(s) with ML (instead of REML)## Data: da
## Models:
## m1: micro ~ posi + (1 | ue1) + (1 | ue2)
## m0: micro ~ solo * posi * prof + (1 | ue1) + (1 | ue2)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 5 -154.56 -146.65 82.281 -164.56
## m0 15 -145.08 -121.32 87.538 -175.08 10.514 10 0.3966anova(m1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## posi 0.0030066 0.0030066 1 5.0005 6.983 0.04583 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emm <- emmeans(m1, ~posi)
as.data.frame(emm)## posi emmean SE df lower.CL upper.CL
## 1 Interrow 0.4161111 0.009043868 9.569438 0.3958365 0.4363857
## 2 Row 0.3861111 0.009043868 9.569438 0.3658365 0.4063857results <- multcomp::cld(emm)
results## posi emmean SE df lower.CL upper.CL .group
## Row 0.386 0.00904 9.57 0.366 0.406 1
## Interrow 0.416 0.00904 9.57 0.396 0.436 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05results <- results %>%
as.data.frame() %>%
mutate(cld = c("b", "a"))
gg4 <-
ggplot(data = results,
mapping = aes(x = posi,
y = emmean,
ymin = lower.CL,
ymax = upper.CL)) +
geom_point() +
geom_errorbar(width = 0.1) +
geom_text(mapping = aes(label = sprintf("%0.2f %s",
emmean,
cld)),
angle = 90, vjust = -0.5) +
labs(x = axis_text$posi,
y = axis_text$micro)
tb_results$micro <- resultsCombining results
tb_results %>%
map(select, -.group, -df)## $ds
## prof solo emmean SE lower.CL upper.CL cld
## 1 10-15 KH 0.8691667 0.01814423 0.8301192 0.9082141 a
## 2 0-5 KH 0.8975000 0.01814423 0.8584525 0.9365475 a
## 3 20-25 KH 0.9004167 0.01814423 0.8613692 0.9394641 a
## 4 0-5 GA 0.9233333 0.01814423 0.8842859 0.9623808 b
## 5 20-25 GA 0.9866667 0.01814423 0.9476192 1.0257141 a
## 6 10-15 GA 1.0025000 0.01814423 0.9634525 1.0415475 a
##
## $total
## posi solo emmean SE lower.CL upper.CL cld
## 1 Interrow KH 0.6333333 0.007192112 0.6165368 0.6501298 b
## 2 Row KH 0.6627778 0.007192112 0.6459813 0.6795743 a
## 3 Interrow GA 0.6036111 0.007192112 0.5868146 0.6204076 b
## 4 Row GA 0.6330556 0.007192112 0.6162590 0.6498521 a
##
## $macro
## posi emmean SE lower.CL upper.CL cld
## 1 Interrow 0.2018056 0.01318235 0.1722044 0.2314067 b
## 2 Row 0.2622222 0.01318235 0.2326211 0.2918234 a
##
## $micro
## posi emmean SE lower.CL upper.CL cld
## 1 Row 0.3861111 0.009043868 0.3658365 0.4063857 b
## 2 Interrow 0.4161111 0.009043868 0.3958365 0.4363857 a# tb <- tb_results %>%
# bind_rows()
# tb
gridExtra::grid.arrange(grobs = list(gg1, gg2, gg3, gg4),
layout_matrix = rbind(c(1, 1),
c(2, 2),
c(3, 4)))
Session information
## Atualizado em 16 de janeiro de 2020.
##
## R version 3.6.2 (2019-12-12)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## 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] stats graphics grDevices datasets utils methods
## [7] base
##
## other attached packages:
## [1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3
## [4] purrr_0.3.2 readr_1.3.1 tidyr_0.8.3
## [7] tibble_2.1.3 ggplot2_3.2.0 tidyverse_1.2.1
## [10] car_3.0-3 carData_3.0-2 emmeans_1.3.5.1
## [13] lmerTest_3.1-0 lme4_1.1-21 Matrix_1.2-18
## [16] EACS_0.0-8 wzRfun_0.91 captioner_2.2.3
## [19] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-38
## [22] rmarkdown_1.14 knitr_1.23
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 minqa_1.2.4 colorspace_1.4-1
## [4] rio_0.5.16 rprojroot_1.3-2 estimability_1.3
## [7] fs_1.3.1 rstudioapi_0.10 remotes_2.1.0
## [10] mvtnorm_1.0-11 lubridate_1.7.4 xml2_1.2.0
## [13] codetools_0.2-16 splines_3.6.2 pkgload_1.0.2
## [16] zeallot_0.1.0 jsonlite_1.6 nloptr_1.2.1
## [19] pbkrtest_0.4-7 broom_0.5.2 compiler_3.6.2
## [22] httr_1.4.0 backports_1.1.4 assertthat_0.2.1
## [25] lazyeval_0.2.2 cli_1.1.0 htmltools_0.4.0
## [28] prettyunits_1.0.2 tools_3.6.2 coda_0.19-3
## [31] gtable_0.3.0 glue_1.3.1 reshape2_1.4.3
## [34] Rcpp_1.0.3 cellranger_1.1.0 vctrs_0.2.0
## [37] nlme_3.1-140 xfun_0.8 ps_1.3.0
## [40] openxlsx_4.1.0.1 testthat_2.2.0 rvest_0.3.4
## [43] devtools_2.1.0 MASS_7.3-51.5 zoo_1.8-6
## [46] scales_1.0.0 hms_0.5.0 parallel_3.6.2
## [49] sandwich_2.5-1 yaml_2.2.0 curl_4.0
## [52] memoise_1.1.0 gridExtra_2.3 stringi_1.4.3
## [55] desc_1.2.0 boot_1.3-23 pkgbuild_1.0.3
## [58] zip_2.0.3 rlang_0.4.0 pkgconfig_2.0.2
## [61] evaluate_0.14 labeling_0.3 processx_3.4.1
## [64] tidyselect_0.2.5 plyr_1.8.4 magrittr_1.5
## [67] R6_2.4.0 multcompView_0.1-7 generics_0.0.2
## [70] multcomp_1.4-10 pillar_1.4.2 haven_2.1.1
## [73] foreign_0.8-71 withr_2.1.2 survival_2.44-1.1
## [76] abind_1.4-5 modelr_0.1.4 crayon_1.3.4
## [79] usethis_1.5.1 grid_3.6.2 readxl_1.3.1
## [82] data.table_1.12.2 callr_3.3.1 digest_0.6.21
## [85] xtable_1.8-4 numDeriv_2016.8-1.1 munsell_0.5.0
## [88] sessioninfo_1.1.1
Estatística Aplicada à Ciência do Solo |
github.com/walmes/EACS |