1 Exploratory data analysis

#-----------------------------------------------------------------------
# Load packages.

library(lattice)
library(latticeExtra)
library(gridExtra)
library(tidyverse)

#-----------------------------------------------------------------------
# Load pre-processed data.

rm(list = ls())
load(file = "high_soybean.RData")
load(file = "swrc.RData")
load(file = "phis-chem.RData")
phch <- dc
swrc <- params
rm(dc, params)
ls()

## [1] "high_soybean" "phch"         "swrc"

#-----------------------------------------------------------------------
# Merge tables.

# apropos("join")
da <- full_join(swrc,
                phch %>% select(-prod),
                by = c("id", "cam")) %>%
    rename(Sindex = "S.x",
           S = "S.y",
           Iindex = "I")

# Values of thr outside the interpretable interval.
da$thr[da$thr < -0.5] <- NA
da$paw[da$paw > 100] <- NA

names(da)

##  [1] "id"       "prod"     "idc"      "cam"      "ds"       "poros"    "thr"     
##  [8] "ths"      "tha"      "thn"      "cad"      "U"        "Iindex"   "Sindex"  
## [15] "fwc"      "macrop"   "paw"      "glic"     "sulf"     "phos"     "FLF"     
## [22] "OLF"      "MO2"      "HF"       "rp"       "dens"     "poro"     "clay"    
## [29] "sand"     "pH_H2O"   "pH_CaCl2" "P"        "K"        "Ca_Mg"    "Ca"      
## [36] "Mg"       "Al"       "MO"       "base"     "CTC"      "V"        "CaMg"    
## [43] "CaK"      "MgK"      "Ca_"      "Mg_"      "K_"       "H_"       "m_"      
## [50] "S"        "dp"

db <- da %>%
    select(-cad) %>%
    gather(key = "variable",
           value = "value",
           -(1:4))

ggplot(data = db,
       mapping = aes(x = value, y = prod, color = cam)) +
    geom_point() +
    geom_smooth(span = 1) +
    facet_wrap(facet = ~variable,
               scales = "free_x",
               ncol = 3)

db_stats <- db %>%
    filter(!is.na(value)) %>%
    group_by(variable, cam) %>%
    summarise_at(.vars = "value",
                 .funs = funs(n = length,
                              mean = mean,
                              sd = sd,
                              min = min,
                              q25 = quantile(., 0.25),
                              median = median,
                              q75 = quantile(., 0.75),
                              max = max))

# db_stats %>% print(n = Inf)
kable(db_stats)

variable	cam	n	mean	sd	min	q25	median	q75	max
Al	L1	62	0.0000000	0.0000000	0.0000000	0.0000000	0.0000000	0.0000000	0.0000000
Al	L2	62	0.2854839	0.2031286	0.0000000	0.0000000	0.3000000	0.4000000	0.7000000
Al	L3	63	0.4587302	0.1421436	0.0000000	0.4000000	0.5000000	0.5500000	0.7000000
base	L1	62	4.8751613	0.7284044	3.1600000	4.3250000	4.8300000	5.4650000	6.4500000
base	L2	62	2.2408065	0.7452479	0.9100000	1.6950000	2.1900000	2.6375000	4.2400000
base	L3	63	1.3326984	0.4812926	0.6900000	1.0550000	1.2200000	1.4500000	3.7800000
Ca	L1	62	3.3774194	0.5074589	2.2000000	3.0000000	3.3000000	3.7750000	4.5000000
Ca	L2	62	1.5177419	0.5314762	0.6000000	1.1000000	1.5000000	1.8000000	2.9000000
Ca	L3	63	0.8793651	0.3336788	0.4000000	0.7000000	0.8000000	1.0000000	2.6000000
Ca_	L1	62	36.6629032	4.0584352	28.2000000	33.8250000	36.6500000	39.3750000	45.7000000
Ca_	L2	62	21.8403226	6.9718746	8.9000000	17.0000000	21.4000000	27.7500000	42.6000000
Ca_	L3	63	15.7412698	4.9946951	7.0000000	12.5500000	14.8000000	18.2500000	35.4000000
CaK	L1	62	14.2062903	4.5505900	7.0700000	10.7150000	13.3750000	16.1200000	28.7500000
CaK	L2	62	13.4198387	6.0681675	4.6700000	8.2900000	12.3700000	16.9200000	28.6200000
CaK	L3	63	11.2011111	5.7432100	3.9500000	7.1250000	9.4500000	14.8050000	29.1000000
CaMg	L1	62	2.7288710	0.0820656	2.5000000	2.6700000	2.7300000	2.7975000	2.9200000
CaMg	L2	62	2.5435484	0.1646566	2.0000000	2.5000000	2.5700000	2.6000000	3.0000000
CaMg	L3	63	2.4685714	0.2585621	2.0000000	2.3300000	2.4000000	2.6000000	3.0000000
Ca_Mg	L1	62	4.6161290	0.6914450	3.0000000	4.1000000	4.5000000	5.1000000	6.1000000
Ca_Mg	L2	62	2.1112903	0.7213614	0.8000000	1.6000000	2.1000000	2.5000000	4.0000000
Ca_Mg	L3	63	1.2396825	0.4675053	0.6000000	1.0000000	1.1000000	1.4000000	3.6000000
clay	L1	62	590.0806452	123.2816711	151.0000000	562.2500000	610.0000000	666.2500000	752.0000000
clay	L2	62	649.5645161	123.3395448	195.0000000	624.7500000	682.0000000	727.5000000	770.0000000
clay	L3	63	680.8888889	121.4088150	215.0000000	651.5000000	710.0000000	762.0000000	798.0000000
CTC	L1	62	9.2106452	0.9173570	6.8200000	8.6950000	9.2550000	9.9600000	12.0500000
CTC	L2	62	6.9456452	0.9309181	4.4000000	6.3500000	6.9450000	7.4575000	9.3400000
CTC	L3	63	5.6009524	0.8906245	3.3900000	5.0650000	5.6000000	6.0050000	8.5800000
dens	L1	65	1.2487692	0.1438424	0.9400000	1.1450000	1.2450000	1.3550000	1.6800000
dens	L2	65	1.3988462	0.1245861	1.1950000	1.3150000	1.3550000	1.4450000	1.7000000
dens	L3	65	1.3433846	0.1201858	1.1400000	1.2650000	1.3150000	1.3900000	1.7100000
dp	L1	65	2.6670769	0.0736911	2.5200000	2.6100000	2.6700000	2.7100000	2.8300000
dp	L2	65	2.6756923	0.0666607	2.5500000	2.6200000	2.6700000	2.7300000	2.8500000
dp	L3	65	2.6840000	0.0999187	2.4500000	2.6300000	2.6700000	2.7300000	3.1200000
ds	L1	61	1.2448361	0.1444298	0.9400000	1.1450000	1.2350000	1.3450000	1.6800000
ds	L2	60	1.3997500	0.1229042	1.1950000	1.3187500	1.3550000	1.4450000	1.7000000
ds	L3	62	1.3371774	0.1192484	1.1400000	1.2612500	1.3075000	1.3887500	1.7100000
FLF	L1	62	15.2062903	3.8947414	2.2400000	12.7325000	15.3050000	17.7000000	24.5200000
fwc	L1	61	0.4032167	0.0487485	0.1836859	0.3893453	0.4084403	0.4261199	0.5136955
fwc	L2	60	0.3714857	0.0514525	0.1833500	0.3487409	0.3797535	0.4027146	0.4986296
fwc	L3	62	0.3867923	0.0482252	0.2198367	0.3806986	0.3989490	0.4136445	0.4791805
glic	L1	62	219.8145161	60.2411189	101.7000000	176.5250000	218.3500000	256.3750000	366.0000000
H_	L1	62	47.0693548	5.9241542	32.1000000	43.2500000	47.3500000	51.3500000	58.8000000
H_	L2	62	63.5451613	7.0268963	39.4000000	59.2500000	63.3000000	68.0500000	76.3000000
H_	L3	63	67.7857143	5.2128519	47.2000000	65.0000000	68.6000000	70.5500000	77.4000000
HF	L1	62	14.8433871	5.7456852	2.1000000	9.4675000	16.2850000	18.5300000	30.1000000
Iindex	L1	61	0.8649539	0.4427030	0.1465554	0.5647225	0.7729931	1.0244809	2.7781969
Iindex	L2	60	1.1149772	0.4110527	0.0612063	0.9395932	1.0906270	1.2637338	3.2110558
Iindex	L3	62	0.9065923	0.2190766	0.3191542	0.7745131	0.9300537	1.0674940	1.3020688
K	L1	62	100.9951613	31.1912526	41.5000000	78.2250000	102.4500000	121.3500000	182.1000000
K	L2	62	50.2790323	23.1682594	22.5000000	32.6250000	42.6000000	60.2000000	119.1000000
K	L3	63	36.2777778	16.0416338	14.2000000	22.6500000	33.6000000	47.0500000	76.6000000
K_	L1	62	2.8177419	0.8690534	1.3000000	2.2000000	2.8000000	3.4750000	5.1000000
K_	L2	62	1.8758065	0.8403133	0.8000000	1.2000000	1.6000000	2.4000000	4.1000000
K_	L3	63	1.6920635	0.7775689	0.6000000	1.1000000	1.5000000	2.2500000	3.3000000
m_	L1	62	0.0000000	0.0000000	0.0000000	0.0000000	0.0000000	0.0000000	0.0000000
m_	L2	62	13.2822581	10.7206179	0.0000000	0.0000000	13.2500000	20.2000000	43.5000000
m_	L3	63	27.2428571	10.3650734	0.0000000	21.7000000	26.2000000	35.3000000	50.2000000
macrop	L1	61	0.1402806	0.0635959	0.0184741	0.1093753	0.1352847	0.1625038	0.3674942
macrop	L2	60	0.1266810	0.0755259	0.0143186	0.0943649	0.1068904	0.1295038	0.3972585
macrop	L3	62	0.1313529	0.0753326	0.0308195	0.0950858	0.1110790	0.1372371	0.3784978
Mg	L1	62	1.2387097	0.1867363	0.8000000	1.1000000	1.2000000	1.4000000	1.6000000
Mg	L2	62	0.5935484	0.1915176	0.2000000	0.5000000	0.6000000	0.7000000	1.1000000
Mg	L3	63	0.3603175	0.1362395	0.2000000	0.3000000	0.3000000	0.4000000	1.0000000
Mg_	L1	62	13.4419355	1.5004424	10.4000000	12.5000000	13.4000000	14.6500000	17.3000000
Mg_	L2	62	8.5467742	2.4832368	3.0000000	6.8250000	8.2500000	10.0750000	15.8000000
Mg_	L3	63	6.4412698	2.1039621	3.2000000	5.2000000	6.0000000	7.5500000	14.2000000
MgK	L1	62	5.2170968	1.7057371	2.5700000	3.9550000	4.8250000	5.9675000	10.3200000
MgK	L2	62	5.2774194	2.3368516	1.8400000	3.6175000	4.9950000	6.8150000	11.2300000
MgK	L3	63	4.6406349	2.5000709	1.3200000	2.7100000	4.0500000	6.3800000	11.6400000
MO	L1	62	40.2337097	4.3202399	28.7000000	37.8000000	39.9000000	43.3000000	48.4000000
MO	L2	62	28.3516129	4.7591411	12.3000000	25.6000000	28.7000000	31.2000000	41.0000000
MO	L3	63	20.8746032	4.9048682	8.7000000	17.7500000	20.6000000	23.4000000	36.8000000
MO2	L1	62	35.9933871	5.2523028	20.9400000	32.7075000	36.8950000	38.6300000	45.8600000
OLF	L1	62	5.9437097	2.9852260	1.6500000	4.1775000	5.2150000	6.7450000	14.8900000
P	L1	62	27.0112903	13.2252628	4.6000000	16.4000000	27.7000000	35.9000000	62.3000000
P	L2	62	5.5661290	6.9195362	1.2000000	2.4000000	3.6000000	5.8250000	48.4000000
P	L3	63	1.5888889	1.6119624	0.3000000	0.6000000	1.2000000	1.7000000	8.5000000
paw	L1	60	19.7375446	8.0094138	6.9307587	14.8262876	17.6268704	21.5781597	49.2945718
paw	L2	59	14.6904188	6.1304679	7.4843779	11.1247969	12.7748898	15.4862305	38.3221924
paw	L3	62	28.7130720	7.7263692	11.4345888	24.7453375	27.5985496	30.5619583	62.0684911
pH_CaCl2	L1	62	5.3032258	0.2032262	4.9000000	5.2000000	5.3000000	5.4750000	5.7000000
pH_CaCl2	L2	62	4.6419355	0.2808350	4.2000000	4.5000000	4.5500000	4.8000000	5.6000000
pH_CaCl2	L3	63	4.4111111	0.1884667	4.1000000	4.3000000	4.4000000	4.5000000	5.3000000
pH_H2O	L1	62	6.0612903	0.2075263	5.6000000	5.9000000	6.1000000	6.2000000	6.5000000
pH_H2O	L2	62	5.4129032	0.2802130	4.9000000	5.2000000	5.3000000	5.6000000	6.3000000
pH_H2O	L3	63	5.2126984	0.2011997	4.8000000	5.1000000	5.2000000	5.3000000	6.1000000
phos	L1	62	695.8225806	154.0763643	324.0000000	594.5250000	709.5500000	773.1500000	1021.7000000
poro	L1	65	0.5316154	0.0530559	0.3750000	0.5100000	0.5350000	0.5600000	0.6700000
poro	L2	65	0.4765385	0.0469368	0.3700000	0.4500000	0.4900000	0.5050000	0.5650000
poro	L3	65	0.4990769	0.0438471	0.3900000	0.4750000	0.5100000	0.5300000	0.5700000
poros	L1	61	0.5434973	0.0771235	0.3750000	0.5100000	0.5350000	0.5650000	0.8350000
poros	L2	60	0.4981667	0.0858949	0.3700000	0.4587500	0.4900000	0.5100000	0.7650000
poros	L3	62	0.5181452	0.0726615	0.3950000	0.4862500	0.5150000	0.5337500	0.7650000
rp	L1	65	1.5200000	0.7619768	0.3100000	0.9400000	1.4150000	1.8550000	3.7050000
rp	L2	65	2.0975385	0.5313104	1.1600000	1.7450000	2.0550000	2.4000000	3.9000000
rp	L3	65	1.3952308	0.3947551	0.8550000	1.1100000	1.3300000	1.5200000	2.6900000
S	L1	62	13.1870968	6.2749171	7.0000000	9.0000000	10.7000000	15.4000000	32.9000000
S	L2	62	21.8225806	10.9268870	8.8000000	12.4500000	18.6500000	28.9000000	50.0000000
S	L3	63	43.7650794	18.5706164	13.0000000	30.6000000	46.6000000	54.6000000	81.4000000
sand	L1	62	294.7903226	149.7967078	143.0000000	192.0000000	246.5000000	330.2500000	826.0000000
sand	L2	62	257.5000000	143.6219708	121.0000000	171.0000000	215.5000000	277.2500000	786.0000000
sand	L3	63	237.1587302	137.3266892	111.0000000	151.0000000	192.0000000	262.0000000	763.0000000
Sindex	L1	61	0.1104692	0.0492827	0.0312534	0.0750110	0.1077605	0.1302144	0.3362753
Sindex	L2	60	0.0702067	0.0217119	0.0394904	0.0546909	0.0668596	0.0814770	0.1474765
Sindex	L3	62	0.0832693	0.0221066	0.0444308	0.0685574	0.0815847	0.0978324	0.1517703
sulf	L1	62	89.7741935	35.2941083	34.6000000	65.5000000	78.5000000	113.0000000	199.5000000
tha	L1	61	-0.4278274	0.3488027	-1.1103833	-0.6242082	-0.4393956	-0.2526655	0.8591424
tha	L2	60	-0.6707908	0.3944306	-1.5703649	-0.8884000	-0.7078085	-0.4538719	0.6071773
tha	L3	62	-0.5983923	0.3145982	-1.0553702	-0.8085614	-0.6375017	-0.5081995	0.6542074
thn	L1	61	2.1813124	0.6024322	1.0676862	1.7792934	2.2258327	2.5098557	4.0378579
thn	L2	60	2.1516299	0.6011177	1.0288813	1.7765692	2.1127272	2.4374379	4.2054640
thn	L3	62	2.2797442	0.4657286	1.3617824	1.9436296	2.3589031	2.5636814	3.6494193
thr	L1	60	0.1686951	0.0735586	-0.0993956	0.1297807	0.1922512	0.2184501	0.2739950
thr	L2	59	0.1652974	0.0583768	-0.0262822	0.1301350	0.1859772	0.2056553	0.2623583
thr	L3	62	0.1847529	0.0442600	0.0623564	0.1670140	0.1944361	0.2122446	0.2732352
ths	L1	61	0.4446750	0.0964380	0.1406799	0.3854681	0.4381136	0.4899636	0.7253569
ths	L2	60	0.3452483	0.0615550	0.2113470	0.3079910	0.3515361	0.3784418	0.5336300
ths	L3	62	0.3734625	0.0640256	0.1945221	0.3417181	0.3863289	0.4141717	0.4931953
U	L1	61	0.3311645	0.0685891	0.1093369	0.2922323	0.3295222	0.3699335	0.5243368
U	L2	60	0.2692889	0.0531174	0.1399618	0.2390070	0.2807919	0.2985244	0.3957378
U	L3	62	0.2931210	0.0521210	0.1525955	0.2770946	0.3074810	0.3249296	0.4003916
V	L1	62	52.9293548	5.9247148	41.1600000	48.6875000	52.6200000	56.7700000	67.8700000
V	L2	62	32.2612903	9.6908863	13.5500000	25.5900000	31.7400000	39.6075000	60.6000000
V	L3	63	23.8847619	7.2394561	12.1700000	19.2950000	22.9100000	27.1100000	52.7700000

2 Quantile regression

2.1 Define functions

# Para fazer regressão quantílica.
library(quantreg)
# ls("package:quantreg")

# O que vai aparecer escrito no eixo horizontal.
lay <- c(L1 = "[0, 10] cm",
         L2 = "[10, 20] cm",
         L3 = "[20, 40] cm")

#' Função que ajusta a regressão quantílica e retorna vários artefatos
#' do ajuste.
#'
#' A regressão é feita usando a combinação linear dos valores das
#' camadas dado pela análise de componentes principais.  Os pessos são
#' os carregamentos do primeiro fator.
#'
#' A regressão quantílica é para o quantil 95%.
#'
#' O preditor linear maximal é o polinômio de segundo grau.
fit_model <- function(variable = "CTC",
                      formula = prod ~ x + I(x^2),
                      layer = NULL,
                      xlab = variable) {
    # Prepare the dataset. ---------------------------------------------
    dd <- da %>%
        mutate(x = get(variable)) %>%
        select(id, cam, prod, x) %>%
        spread(key = "cam", value = "x") %>%
        na.omit(.)
    # dd

    # Principal component analysis. ------------------------------------
    # pc <- with(dd, {
    #     princomp(x = cbind(L1, L2, L3))
    # })
    # # summary(pc)    # Cummulated variances.
    # # pc$loadings[,] # Variable loadings (weights).
    # # biplot(pc)     # Biplot of loadings and scores.
    # # Synthetic variable: the 1st score from PCA.
    # # Use positive loadings.
    # s <- ifelse(sum(sign(pc$loadings[, 1])) < -1,
    #             yes = -1,
    #             no = 1)
    # # dd$x <- s * pc$scores[, 1]
    # eqn <- s * pc$loadings[, 1]
    # dd$x <- c(with(dd, cbind(L1, L2, L3)) %*% cbind(eqn))

    # Choose a layer to be used. ---------------------------------------
    if (!is.null(layer)) {
        dd[, "x"] <- dd[, layer]
    } else {
        dd[, "x"] <- dd[, "L2"]
    }

    # Fit the 90th quantile regression. -------------------------------------
    m0 <- rq(formula,
             tau = 0.9,
             data = dd)
    # summary(m0)

    # Null model and LRT. ----------------------------------------------
    m00 <- update(m0, . ~ 1)
    # anova(m0, m00)

    # Grid of values for prediction. -----------------------------------
    minmax <- extendrange(dd$x, f = 0.05)
    grid <- with(dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
    grid$prod <- predict(m0, newdata = grid)

    # Plot the results. ------------------------------------------------
    gg <- ggplot(data = dd,
                 mapping = aes(x = x, y = prod)) +
        geom_point() +
        geom_line(data = grid,
                  mapping = aes(x = x, y = prod),
                  col = "tomato") +
        ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
        xlab(label = xlab)

    # List of objects. -------------------------------------------------
    return(list(dd = dd,
                # pc = pc,
                m0 = m0,
                # m00 = m00,
                anova = anova(m0, m00),
                grid = grid,
                gg = gg))
}

#' Função que faz o gráfico conforme o tipo de função.  Se for uma
#' função quadrática côncava, um par de valores é determinado.  Se for
#' uma função linear, apenas um limite é determinado.
plot_curve <- function(fit_obj, tol = 0.1) {
    # Coeficients.
    limits <- list(y_max = NA,
                   y_tol = NA,
                   x_tol_inf = NA,
                   x_tol_sup = NA)
    gg <- fit_obj$gg
    coefs <- coef(fit_obj$m0)
    if (length(coefs) == 3 && coefs[3] < 0) { #-------------------------
        coefs <- structure(coefs,
                           .Names = letters[3:1])
        x_crit <- with(as.list(coefs), {
            -b/(2 * a)
        })
        y_crit <- x_crit^(0:2) %*% coefs
        roots <- with(as.list(coefs), {
            delta <- sqrt(b^2 - 4 * a * (c - y_crit + tol))
            suppressWarnings({
                (-b + c(-1, 1) * delta)/(2 * a)
            })
        })
        limits <- c(y_max = y_crit,
                    y_tol = y_crit - tol,
                    x_tol_inf = min(roots),
                    x_tol_sup = max(roots))
        gg <- fit_obj$gg +
            geom_segment(data = data.frame(x = x_crit,
                                           y = y_crit),
                         mapping = aes(x = x,
                                       y = -Inf,
                                       xend = x,
                                       yend = y),
                         col = "green",
                         linetype = 3) +
            geom_segment(data = data.frame(x = x_crit,
                                           y = y_crit),
                         mapping = aes(x = x,
                                       y = y,
                                       xend = -Inf,
                                       yend = y),
                         col = "green",
                         linetype = 3) +
            geom_hline(yintercept = y_crit - tol,
                       col = "blue",
                       linetype = 3) +
            geom_segment(data = data.frame(x = roots,
                                           y = y_crit - tol),
                         mapping = aes(x = x,
                                       y = -Inf,
                                       xend = x,
                                       yend = y),
                         col = "blue",
                         linetype = 3)
    } else if (length(coefs) == 2) { #----------------------------------
        i <- which.max(fit_obj$grid$prod)
        beta <- coef(fit_obj$m0)["x"]
        y_max <- fit_obj$grid$prod[i]
        x_max <- fit_obj$grid$x[i]
        x_tol <- x_max - tol/beta
        names(x_tol) <- NULL
        limits <- c(y_max = y_max,
                    y_tol = y_max - tol,
                    x_tol_inf = ifelse(x_tol < x_max, x_tol, -Inf),
                    x_tol_sup = ifelse(x_tol < x_max, Inf, x_tol))
        col <- c("green", "blue")
        gg <- fit_obj$gg +
            geom_vline(xintercept = c(x_max, x_tol),
                       col = col,
                       linetype = 3) +
            geom_hline(yintercept = c(y_max, y_max - tol),
                       col = col,
                       linetype = 3)
    }
    return(list(limits, gg))
}

2.2 Residual water content (thr)

# Legendas.
variable <- "thr"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Residual water content in" ~ lay ~ (g ~ g^{-1}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.4581317   0.6379005 
## 
## Degrees of freedom: 52 total; 50 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    5.111472   -2.785691 
## 
## Degrees of freedom: 52 total; 50 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.4446509   0.6333507 
## 
## Degrees of freedom: 52 total; 50 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       50  0.9732 0.3286
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       50  1.3604  0.249
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       50  0.0602 0.8071

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1           L2        L3
## y_max     4.6345243  5.224888890 4.6243826
## y_tol     4.5345243  5.124888890 4.5243826
## x_tol_inf 0.1197562         -Inf 0.1258888
## x_tol_sup       Inf -0.004816508       Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.3 Water content at saturation (ths)

# Legendas.
variable <- "ths"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Water content at saturation in" ~ lay ~ (g ~ g^{-1}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    4.423477    0.295334 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.3853115   0.5377847 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.7521541  -0.4862805 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0497 0.8245
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0605 0.8067
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0614 0.8052

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1        L2        L3
## y_max     4.6248344 4.6304979 4.6461934
## y_tol     4.5248344 4.5304979 4.5461934
## x_tol_inf 0.3431947 0.2699711      -Inf
## x_tol_sup       Inf       Inf 0.4235428

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.4 Soil S-index (slope at inflexion point)

# Legendas.
variable <- "Sindex"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Soil S-index in" ~ lay,
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(2, 2, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##    3.761762   11.595362  -28.926111 
## 
## Degrees of freedom: 54 total; 51 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##    3.030303   38.087871 -230.139148 
## 
## Degrees of freedom: 54 total; 51 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    4.034551    5.930447 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value    Pr(>F)    
## 1  2       51  16.198 3.579e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value   Pr(>F)    
## 1  2       51  28.929 4.01e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  1.4082 0.2408

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1         L2        L3
## y_max     4.9237954 4.60618182 4.9664459
## y_tol     4.8237954 4.50618182 4.8664459
## x_tol_inf 0.1416337 0.06190452 0.1402752
## x_tol_sup 0.2592277 0.10359479       Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.5 Soil I-index (tension at inflexion point)

# Legendas.
variable <- "Iindex"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Soil I-index in" ~ lay,
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 2, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.6834845  -0.1550112 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##   3.9067991   0.8769461  -0.2508731 
## 
## Degrees of freedom: 54 total; 51 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.7032942  -0.1745639 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.6084 0.4389
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value   Pr(>F)   
## 1  2       51  5.5333 0.006686 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.1461 0.7038

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1       L2       L3
## y_max     4.6542855 4.673157 4.636356
## y_tol     4.5542855 4.573157 4.536356
## x_tol_inf      -Inf 1.116434     -Inf
## x_tol_sup 0.8334819 2.379142 0.956314

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.6 Macroporosity

expression(“Macroporosity” ~ (m^3 ~ m^{-3}))

# Legendas.
variable <- "macrop"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Macroporosity in" ~ lay ~ (m^3 ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    4.301982    1.900857 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    4.252295    2.214228 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.5765138  -0.0450462 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  1.1464 0.2893
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  1.0424  0.312
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52   4e-04 0.9851

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1       L2       L3
## y_max     4.6890028 5.174312 4.575677
## y_tol     4.5890028 5.074312 4.475677
## x_tol_inf 0.1509955 0.371243     -Inf
## x_tol_sup       Inf      Inf 2.238525

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.7 Field water capacity (microporosity)

# Legendas.
variable <- "fwc"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Field water capacity (microporosity) in" ~ lay ~ (m^3 ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 2)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.3855800   0.4035359 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.4003447   0.3837563 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##  -0.5027545  29.0108513 -41.0261828 
## 
## Degrees of freedom: 54 total; 51 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0289 0.8657
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0213 0.8845
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  2       51  4.1452 0.02148 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1        L2        L3
## y_max     4.5995331 4.5728024 4.6258575
## y_tol     4.4995331 4.4728024 4.5258575
## x_tol_inf 0.2823866 0.1888117 0.3041944
## x_tol_sup       Inf       Inf 0.4029358

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.8 Plant available water

# Legendas.
variable <- "paw"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Plant available water in" ~ lay ~ (m^3 ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 2, 2)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
##  4.687158260 -0.005587313 
## 
## Degrees of freedom: 52 total; 50 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##  3.05008122  0.15477010 -0.00312627 
## 
## Degrees of freedom: 52 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x       I(x^2) 
##  3.651998889  0.069766506 -0.001316386 
## 
## Degrees of freedom: 52 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       50  0.1215 0.7288
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value   Pr(>F)   
## 1  2       49  7.9454 0.001026 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value   Pr(>F)   
## 1  2       49  7.8172 0.001131 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1        L2        L3
## y_max      4.660269  4.965605  4.576379
## y_tol      4.560269  4.865605  4.476379
## x_tol_inf      -Inf 19.097448 17.783435
## x_tol_sup 22.710259 30.408858 35.215080

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.9 Soil resistence penetration

# Legendas.
variable <- "rp"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Soil resistence penetration in" ~ lay ~ (MPa),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.7612571  -0.1496104 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.9428660  -0.1938144 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.7178419  -0.1246512 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52   0.631 0.4306
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.5237 0.4725
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.1688 0.6829

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1       L2       L3
## y_max     4.7370577 4.744594 4.622702
## y_tol     4.6370577 4.644594 4.522702
## x_tol_inf      -Inf     -Inf     -Inf
## x_tol_sup 0.8301528 1.538957 1.565489

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.10 Soil bulk density

# Legendas.
variable <- "dens"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Soil bulk density in" ~ lay ~ (Mg ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##      6.0336     -1.1800 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.8721846  -0.2215385 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    2.418875    1.625000 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  2.6008 0.1129
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0463 0.8304
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  1       52  3.4939 0.06723 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                 L1       L2       L3
## y_max     4.937085 4.613040 5.243937
## y_tol     4.837085 4.513040 5.143937
## x_tol_inf     -Inf     -Inf 1.676962
## x_tol_sup 1.013996 1.621139      Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.11 Particle density

# Legendas.
variable <- "dp"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Particle density in" ~ lay ~ (Mg ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 2)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##        9.92       -2.00 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   -2.822000    2.733333 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##  -23.873408   20.694955   -3.763237 
## 
## Degrees of freedom: 54 total; 51 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  1.0413 0.3122
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value    Pr(>F)    
## 1  1       52  12.195 0.0009865 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df    F value    Pr(>F)    
## 1  2       51 4355418852 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##              L1       L2       L3
## y_max     4.888 5.006267 4.578241
## y_tol     4.788 4.906267 4.478241
## x_tol_inf  -Inf 2.827415 2.586610
## x_tol_sup 2.566      Inf 2.912634

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.12 Total porosity

# Legendas.
variable <- "poro"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Total porosity in" ~ lay ~ (m ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    3.948000    1.066667 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##    3.602769    1.938462 
## 
## Degrees of freedom: 54 total; 52 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##       4.719      -0.300 
## 
## Degrees of freedom: 54 total; 52 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.1913 0.6636
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.4738 0.4943
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       52  0.0137 0.9074

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                L1        L2        L3
## y_max     4.62800 4.6965462 4.6046250
## y_tol     4.52800 4.5965462 4.5046250
## x_tol_inf 0.54375 0.5126627      -Inf
## x_tol_sup     Inf       Inf 0.7145833

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.13 Clay content

# Legendas.
variable <- "clay"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Clay content in" ~ lay ~ (m^3 ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##   (Intercept)             x 
##  4.6358571429 -0.0001071429 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.4680000000 0.0001212121 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.4624847251 0.0001221996 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0058 0.9398
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0053  0.942
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0058 0.9396

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                    L1         L2       L3
## y_max        4.622898   4.564818 4.563562
## y_tol        4.522898   4.464818 4.463562
## x_tol_inf        -Inf -26.250000 8.816667
## x_tol_sup 1054.283333        Inf      Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.14 Sand content

# Legendas.
variable <- "sand"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Sand content in" ~ lay ~ (m^3 ~ m^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.546647e+00 9.022556e-05 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.5470608696 0.0001043478 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.5478909091 0.0001090909 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0058 0.9394
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0058 0.9397
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0057 0.9401

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                    L1          L2          L3
## y_max        4.624254    4.632532    4.634684
## y_tol        4.524254    4.532532    4.534684
## x_tol_inf -248.183333 -139.233333 -121.066667
## x_tol_sup         Inf         Inf         Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.15 pH on H₂O

# Legendas.
variable <- "pH_H2O"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute(pH ~ H[2] * O ~ "in" ~ lay,
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   6.4426667  -0.3133333 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##      2.5845      0.3750 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##     -0.8410      1.0225 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.2755  0.602
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.9074 0.3455
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  1       49   3.129 0.08314 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                 L1       L2       L3
## y_max     4.702100 4.973250 5.462712
## y_tol     4.602100 4.873250 5.362712
## x_tol_inf     -Inf 6.103333 6.067200
## x_tol_sup 5.874149      Inf      Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.16 Phosphorus (P)

# Legendas.
variable <- "P"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Phosphorus in" ~ lay ~ (P * ", " ~ mg ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
##  4.583773585 -0.001509434 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x 
##  4.580000000 -0.005714286 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.62644737 -0.03605263 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0295 0.8643
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0262 0.8721
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.4641 0.4989

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1        L2       L3
## y_max      4.581185  4.586629 4.630413
## y_tol      4.481185  4.486629 4.530413
## x_tol_inf      -Inf      -Inf     -Inf
## x_tol_sup 67.965000 16.340000 2.663723

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.17 Potassium (K)

# Legendas.
variable <- "K"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Potassium in" ~ lay ~ (K * ", " * mg ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.3693441   0.0020954 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.06457143  0.01142857 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
## 4.332044444 0.006074074 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.4172 0.5214
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  1       49  6.4705 0.01417 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49    1.06 0.3083

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                   L1        L2       L3
## y_max       4.765647  5.305714  4.81627
## y_tol       4.665647  5.205714  4.71627
## x_tol_inf 141.406423 99.850000 63.25659
## x_tol_sup        Inf       Inf      Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.18 Calcium (Ca)

# Legendas.
variable <- "Ca"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Calcium in" ~ lay ~ (Ca * ", " * cmolc ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 2, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##       4.356       0.068 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##   3.4820000   1.2178022  -0.2989011 
## 
## Degrees of freedom: 51 total; 48 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.3748571   0.2314286 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0783 0.7807
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value    Pr(>F)    
## 1  2       48  9.8237 0.0002653 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.4457 0.5075

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                 L1       L2       L3
## y_max     4.648400 4.722412 4.734729
## y_tol     4.548400 4.622412 4.634729
## x_tol_inf 2.829412 1.458722 1.122901
## x_tol_sup      Inf 2.615543      Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.19 Magnesium (Mg)

# Legendas.
variable <- "Mg"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Magnesium in" ~ lay ~ (Mg * ", " * cmolc ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(2, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x      I(x^2) 
##  -0.2171429   7.8392857  -3.1785714 
## 
## Degrees of freedom: 51 total; 48 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##       3.861       1.230 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##       4.278       0.740 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value   Pr(>F)   
## 1  2       48  6.5901 0.002959 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  1       49   3.798 0.05705 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.6378 0.4284

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                 L1       L2        L3
## y_max     4.616349 5.269350 4.7368000
## y_tol     4.516349 5.169350 4.6368000
## x_tol_inf 1.055774 1.063699 0.4848649
## x_tol_sup 1.410518      Inf       Inf

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.20 Aluminium (Al)

# Legendas.
variable <- "Al"
xlab <- lapply(lay[2:3],
               FUN = function(l) {
                   substitute("Aluminium in" ~ lay ~ (cmolc ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1)],
               layer = names(lay)[2:3],
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay)[2:3])

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##       4.722      -0.375 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##        4.85       -0.64 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.3641  0.549
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.5555 0.4596

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L2      L3
## y_max     4.7351250 4.87240
## y_tol     4.6351250 4.77240
## x_tol_inf      -Inf    -Inf
## x_tol_sup 0.2316667 0.12125

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.21 Sulfur (S)

# Legendas.
variable <- "S"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("Sulfur in" ~ lay ~ (mg ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.69871739 -0.01456522 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##   4.8705283  -0.0154717 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.74627129 -0.00384858 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.9494 0.3347
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  1.0907 0.3014
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.3196 0.5744

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1        L2        L3
## y_max      4.615623  4.763774  4.709152
## y_tol      4.515623  4.663774  4.609152
## x_tol_inf      -Inf      -Inf      -Inf
## x_tol_sup 12.570672 13.363415 35.628607

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.22 CTC

# Legendas.
variable <- "CTC"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("CTC pH 7 in" ~ lay ~ (cmolc ~ dm^{-3}),
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(1, 1, 1)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.21056604  0.03773585 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.24195364  0.03907285 
## 
## Degrees of freedom: 51 total; 49 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
## (Intercept)           x 
##  4.64226804 -0.01237113 
## 
## Degrees of freedom: 51 total; 49 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.0672 0.7966
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49  0.1096 0.7421
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       49   0.042 0.8385

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                 L1       L2        L3
## y_max     4.611358 4.616545  4.603367
## y_tol     4.511358 4.516545  4.503367
## x_tol_inf 7.971000 7.027678      -Inf
## x_tol_sup      Inf      Inf 11.227833

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.23 V (%)

# Legendas.
variable <- "V"
xlab <- lapply(lay,
               FUN = function(l) {
                   substitute("V (%) in" ~ lay,
                              list(lay = l))
               })

# Ajustes.
fits <- mapply(FUN = "fit_model",
               formula = c(prod ~ x, prod ~ x + I(x^2))[c(2, 2, 2)],
               layer = names(lay),
               xlab = xlab,
               MoreArgs = list(variable = variable),
               SIMPLIFY = FALSE)
fits <- setNames(fits, names(lay))

# Resumo dos ajustes em cada camada.
lapply(fits, FUN = "[[", "m0")

## $L1
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x       I(x^2) 
##  0.929495705  0.139878108 -0.001324919 
## 
## Degrees of freedom: 51 total; 48 residual
## 
## $L2
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##   (Intercept)             x        I(x^2) 
##  3.2505750775  0.0722525403 -0.0008855925 
## 
## Degrees of freedom: 51 total; 48 residual
## 
## $L3
## Call:
## rq(formula = formula, tau = 0.9, data = dd)
## 
## Coefficients:
##  (Intercept)            x       I(x^2) 
##  3.194872195  0.088272402 -0.001268733 
## 
## Degrees of freedom: 51 total; 48 residual

# Teste conta o modelo nulo.
lapply(fits, FUN = "[[", "anova")

## $L1
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  2       48  2.8937 0.06508 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L2
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value    Pr(>F)    
## 1  2       48  11.308 9.465e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $L3
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x + I(x^2)
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  2       48  2.8056 0.07041 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Limites de suficiência.
plt <- lapply(fits, FUN = plot_curve)
sapply(plt, "[[", 1)

##                  L1        L2        L3
## y_max      4.621398  4.724286  4.730265
## y_tol      4.521398  4.624286  4.630265
## x_tol_inf 44.099711 30.167010 25.909621
## x_tol_sup 61.475133 51.419657 43.665609

# Gráficos.
do.call(function(...) grid.arrange(..., nrow = 1),
        lapply(plt, "[[", 2))

2.24 glic

fit <- list()

xlab <- "Glycosidase"
fit$dd <- da %>%
    select(id, cam, prod, glic) %>%
    na.omit(.) %>%
    rename(x = glic)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.4509508197 0.0004371585 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       57  0.0407 0.8408

plot_curve(fit)

## [[1]]
##      y_max      y_tol  x_tol_inf  x_tol_sup 
##   4.616728   4.516728 150.465000        Inf 
## 
## [[2]]

2.25 Phosphatase

fit <- list()

xlab <- "Phosphatase"
fit$dd <- da %>%
    select(id, cam, prod, phos) %>%
    na.omit(.) %>%
    rename(x = phos)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
##  (Intercept)            x 
## 4.4610264635 0.0001202887 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       57   0.108 0.7436

plot_curve(fit)

## [[1]]
##      y_max      y_tol  x_tol_inf  x_tol_sup 
##   4.588122   4.488122 225.251667        Inf 
## 
## [[2]]

2.26 Sulfatase

fit <- list()

xlab <- "Sulfatase"
fit$dd <- da %>%
    select(id, cam, prod, sulf) %>%
    na.omit(.) %>%
    rename(x = sulf)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
## (Intercept)           x 
## 4.105056478 0.004983389 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value  Pr(>F)  
## 1  1       57  2.9143 0.09324 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_curve(fit)

## [[1]]
##      y_max      y_tol  x_tol_inf  x_tol_sup 
##   5.140331   5.040331 187.678333        Inf 
## 
## [[2]]

2.27 Organic matter

fit <- list()

xlab <- expression("Organic matter" ~ (g ~ kg^{-1}))
fit$dd <- da %>%
    select(id, cam, prod, MO2) %>%
    na.omit(.) %>%
    rename(x = MO2)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
## (Intercept)           x 
## 4.432560386 0.003220612 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       57  0.0254 0.8739

plot_curve(fit)

## [[1]]
##     y_max     y_tol x_tol_inf x_tol_sup 
##  4.584271  4.484271 16.056000       Inf 
## 
## [[2]]

2.28 Free light fractions of SOC

fit <- list()

xlab <- expression("Free light fractions of SOC" ~ (g ~ kg^{-1}))
fit$dd <- da %>%
    select(id, cam, prod, FLF) %>%
    na.omit(.) %>%
    rename(x = FLF)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
## (Intercept)           x 
## 4.431929777 0.005671843 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       57  0.0547 0.8159

plot_curve(fit)

## [[1]]
##     y_max     y_tol x_tol_inf x_tol_sup 
##  4.577322  4.477322  8.003048       Inf 
## 
## [[2]]

2.29 Ocluded light fractions of SOC

fit <- list()

xlab <- expression("Ocluded light fractions of SOC" ~ (g ~ kg^{-1}))
fit$dd <- da %>%
    select(id, cam, prod, OLF) %>%
    na.omit(.) %>%
    rename(x = OLF)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
## (Intercept)           x 
## 4.492522659 0.004531722 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       57  0.0149 0.9034

plot_curve(fit)

## [[1]]
##     y_max     y_tol x_tol_inf x_tol_sup 
##  4.563000  4.463000 -6.514667       Inf 
## 
## [[2]]

2.30 Heavy fractions of SOC

fit <- list()

xlab <- expression("Heavy fractions of SOC" ~ (g ~ kg^{-1}))
fit$dd <- da %>%
    select(id, cam, prod, HF) %>%
    na.omit(.) %>%
    rename(x = HF)

fit$m0 <- rq(c(prod ~ x, prod ~ x + I(x^2))[[1]],
             tau = 0.9,
             data = fit$dd)
fit$m00 <- update(fit$m0, . ~ 1)
fit$anova <- anova(fit$m0, fit$m00)

minmax <- extendrange(fit$dd$x, f = 0.05)
fit$grid <- with(fit$dd,
                 data.frame(x = seq(minmax[1],
                                    minmax[2],
                                    length.out = 150)))
fit$grid$prod <- predict(fit$m0, newdata = fit$grid)

fit$gg <- ggplot(data = fit$dd,
                 mapping = aes(x = x, y = prod)) +
    geom_point() +
    geom_line(data = fit$grid,
              mapping = aes(x = x, y = prod),
              col = "tomato") +
    ylab(label = expression("Soybean productivity" ~ (t ~ ha^{-1}))) +
    xlab(label = xlab)

fit[c("m0", "anova")]

## $m0
## Call:
## rq(formula = c(prod ~ x, prod ~ x + I(x^2))[[1]], tau = 0.9, 
##     data = fit$dd)
## 
## Coefficients:
##  (Intercept)            x 
##  4.570333919 -0.004920914 
## 
## Degrees of freedom: 59 total; 57 residual
## 
## $anova
## Quantile Regression Analysis of Deviance Table
## 
## Model 1: prod ~ x
## Model 2: prod ~ 1
##   Df Resid Df F value Pr(>F)
## 1  1       57  0.2176 0.6427

plot_curve(fit)

## [[1]]
##     y_max     y_tol x_tol_inf x_tol_sup 
##  4.566889  4.466889      -Inf 21.021429 
## 
## [[2]]

Quantile regression to determine suffiency frontier to soybean crop production under no-till

Milson E. Serafim & Walmes M. Zeviani

Last modified at 2018-07-08 15:32:41

1 Exploratory data analysis

2 Quantile regression

2.1 Define functions

2.2 Residual water content (thr)

2.3 Water content at saturation (ths)

2.4 Soil S-index (slope at inflexion point)

2.5 Soil I-index (tension at inflexion point)

2.6 Macroporosity

2.7 Field water capacity (microporosity)

2.8 Plant available water

2.9 Soil resistence penetration

2.10 Soil bulk density

2.11 Particle density

2.12 Total porosity

2.13 Clay content

2.14 Sand content

2.15 pH on H₂O

2.16 Phosphorus (P)

2.17 Potassium (K)

2.18 Calcium (Ca)

2.19 Magnesium (Mg)

2.20 Aluminium (Al)

2.21 Sulfur (S)

2.22 CTC

2.23 V (%)

2.24 glic

2.25 Phosphatase

2.26 Sulfatase

2.27 Organic matter

2.28 Free light fractions of SOC

2.29 Ocluded light fractions of SOC

2.30 Heavy fractions of SOC

Quantile regression to determine suffiency frontier to soybean crop production under no-till

Milson E. Serafim & Walmes M. Zeviani

Last modified at 2018-07-08 15:32:41

1 Exploratory data analysis

2 Quantile regression

2.1 Define functions

2.2 Residual water content (thr)

2.3 Water content at saturation (ths)

2.4 Soil S-index (slope at inflexion point)

2.5 Soil I-index (tension at inflexion point)

2.6 Macroporosity

2.7 Field water capacity (microporosity)

2.8 Plant available water

2.9 Soil resistence penetration

2.10 Soil bulk density

2.11 Particle density

2.12 Total porosity

2.13 Clay content

2.14 Sand content

2.15 pH on H2O

2.16 Phosphorus (P)

2.17 Potassium (K)

2.18 Calcium (Ca)

2.19 Magnesium (Mg)

2.20 Aluminium (Al)

2.21 Sulfur (S)

2.22 CTC

2.23 V (%)

2.24 glic

2.25 Phosphatase

2.26 Sulfatase

2.27 Organic matter

2.28 Free light fractions of SOC

2.29 Ocluded light fractions of SOC

2.30 Heavy fractions of SOC

2.15 pH on H₂O