Soil water rentention curve
of soybean crop under no-till
Milson E. Serafim & Walmes M. Zeviani
Last modified at 2018-08-05 12:19:31
Abstract
These analysis are part of the Milson’s posdoc. The objective is to find the frontier of nutrient sufficiency to high soybean productions. This document does preliminary analysis on the soil water retention curve and its relation with the soybean productivity. Simple and multiple linear regression is applied to explore those relations.
- 1 Exploratory data analysis
- 2 Fitting SWRC model
- 3 Regression for the soybean yield on each SWRC parameter
- 3.1 Calculate values based on soil density and porosity
- 3.2 Exploratory data analysis
- 3.3 S index
- 3.4 Inflexion point (I index)
- 3.5 Macroporosity
- 3.6 Field water capacity (microporosity)
- 3.7 Plant available water
- 3.8 Water content at saturation (ths)
- 3.9 Residual water content (thr)
- 3.10 \(\alpha\) empirical shape parameter (tha)
- 3.11 \(n\) empirical shape parameter (thn)
- 4 Regression for the soybean yield on all SWRC parameters
- 5 Regression for the soybean yield on all physical and hydrical variables
- 6 Additional relations between soil variables
1 Exploratory data analysis
#-----------------------------------------------------------------------
# Packages.
library(lattice)
library(latticeExtra)
library(gridExtra)
library(nlme)
library(tidyverse)#-----------------------------------------------------------------------
# Load pre-processed data.
rm(list = ls())
load(file = "high_soybean.RData")
str(high_soybean)## List of 4
## $ atrib :'data.frame': 394 obs. of 39 variables:
## ..$ faz : chr [1:394] "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" ...
## ..$ tal : chr [1:394] "AG2" "AG2" "AG2" "AG2" ...
## ..$ cam : chr [1:394] "0-0.1" "0-0.1" "0.1-0.2" "0.1-0.2" ...
## ..$ rpt : int [1:394] 1 2 1 2 1 2 1 2 1 2 ...
## ..$ prod : num [1:394] 3.89 NA NA NA NA ...
## ..$ glic : num [1:394] 275 NA NA NA NA ...
## ..$ sulf : num [1:394] 200 NA NA NA NA ...
## ..$ phos : num [1:394] 922 NA NA NA NA ...
## ..$ FLF : num [1:394] 14.4 NA NA NA NA ...
## ..$ OLF : num [1:394] 6.21 NA NA NA NA NA 3.1 NA NA NA ...
## ..$ MO2 : num [1:394] 45.5 NA NA NA NA ...
## ..$ HF : num [1:394] 24.8 NA NA NA NA NA 30.1 NA NA NA ...
## ..$ rp : num [1:394] 2.09 1.75 2.39 2.9 2.35 1.14 0.54 1.26 0.94 2.65 ...
## ..$ dens : num [1:394] 1.23 1.19 1.28 1.3 1.23 1.2 1.01 1.08 1.15 1.24 ...
## ..$ poro_e : num [1:394] 0.53 0.54 0.53 0.52 0.5 0.51 0.61 0.59 0.57 0.53 ...
## ..$ clay : int [1:394] 653 NA 725 NA 775 NA 644 NA 724 NA ...
## ..$ sand : int [1:394] 162 NA 139 NA 120 NA 175 NA 138 NA ...
## ..$ pH_H2O : num [1:394] 6.2 NA 5.3 NA 5.1 NA 6.1 NA 5.2 NA ...
## ..$ pH_CaCl2: num [1:394] 5.4 NA 4.5 NA 4.4 NA 5.2 NA 4.5 NA ...
## ..$ P : num [1:394] 16.4 NA 2.1 NA 0.9 NA 14.4 NA 1.7 NA ...
## ..$ K : num [1:394] 97.3 NA 76.9 NA 67.7 NA 78.6 NA 25.6 NA ...
## ..$ Ca_Mg : num [1:394] 5.5 NA 2.1 NA 1.1 NA 5.1 NA 1.8 NA ...
## ..$ Ca : num [1:394] 4 NA 1.5 NA 0.8 NA 3.7 NA 1.3 NA ...
## ..$ Mg : num [1:394] 1.5 NA 0.6 NA 0.3 NA 1.4 NA 0.5 NA ...
## ..$ Al : num [1:394] 0 NA 0.4 NA 0.5 NA 0 NA 0.4 NA ...
## ..$ MO : num [1:394] 43.3 NA 31.2 NA 21.3 NA 44.5 NA 29.5 NA ...
## ..$ base : num [1:394] 5.75 NA 2.3 NA 1.27 NA 5.3 NA 1.87 NA ...
## ..$ CTC : num [1:394] 10.15 NA 7.6 NA 5.77 ...
## ..$ V : num [1:394] 56.6 NA 30.2 NA 22.1 ...
## ..$ CaMg : num [1:394] 2.67 NA 2.5 NA 2.67 NA 2.64 NA 2.6 NA ...
## ..$ CaK : num [1:394] 16.03 NA 7.61 NA 4.61 ...
## ..$ MgK : num [1:394] 6.01 NA 3.04 NA 1.73 NA 6.95 NA 7.62 NA ...
## ..$ Ca_ : num [1:394] 39.4 NA 19.7 NA 13.9 NA 36.6 NA 17.9 NA ...
## ..$ Mg_ : num [1:394] 14.8 NA 7.9 NA 5.2 NA 13.9 NA 6.9 NA ...
## ..$ K_ : num [1:394] 2.5 NA 2.6 NA 3 NA 2 NA 0.9 NA ...
## ..$ H_ : num [1:394] 43.4 NA 64.5 NA 69.3 NA 47.5 NA 68.8 NA ...
## ..$ m_ : num [1:394] 0 NA 14.8 NA 28.2 NA 0 NA 17.7 NA ...
## ..$ S : num [1:394] 21.1 NA 33.9 NA 56 NA 10.8 NA 32.7 NA ...
## ..$ dp : num [1:394] 2.59 2.59 2.73 2.73 2.45 2.45 2.63 2.63 2.65 2.65 ...
## $ swrc :'data.frame': 390 obs. of 15 variables:
## ..$ faz : chr [1:390] "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" ...
## ..$ tal : chr [1:390] "AG2" "AG2" "AG2" "AG2" ...
## ..$ cam : chr [1:390] "0-0.1" "0-0.1" "0.1-0.2" "0.1-0.2" ...
## ..$ rpt : int [1:390] 1 2 1 2 1 2 1 2 1 2 ...
## ..$ S0 : num [1:390] 0.441 0.445 0.39 0.381 0.421 0.437 0.579 0.518 0.477 0.402 ...
## ..$ S1 : num [1:390] 0.425 0.418 0.35 0.367 0.396 0.41 0.527 0.501 0.431 0.398 ...
## ..$ S2 : num [1:390] 0.402 0.404 0.347 0.358 0.387 0.408 0.457 0.459 0.422 0.389 ...
## ..$ S4 : num [1:390] 0.376 0.384 0.34 0.349 0.379 0.403 0.407 0.422 0.395 0.378 ...
## ..$ S6 : num [1:390] 0.366 0.368 0.335 0.34 0.369 0.394 0.389 0.406 0.384 0.366 ...
## ..$ S8 : num [1:390] 0.357 0.361 0.328 0.331 0.359 0.381 0.376 0.4 0.376 0.359 ...
## ..$ S10 : num [1:390] 0.354 0.356 0.321 0.325 0.352 0.37 0.355 0.395 0.371 0.352 ...
## ..$ S33 : num [1:390] 0.314 0.319 0.286 0.283 0.294 0.29 0.314 0.359 0.331 0.302 ...
## ..$ S66 : num [1:390] 0.296 0.314 0.283 0.276 0.287 0.282 0.291 0.304 0.319 0.284 ...
## ..$ S300 : num [1:390] 0.27 0.287 0.253 0.253 0.261 0.259 0.272 0.297 0.284 0.269 ...
## ..$ S1500: num [1:390] 0.223 0.23 0.225 0.222 0.237 0.242 0.261 0.245 0.284 0.264 ...
## $ fract :'data.frame': 390 obs. of 10 variables:
## ..$ faz : chr [1:390] "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" ...
## ..$ tal : chr [1:390] "AG2" "AG2" "AG2" "AG2" ...
## ..$ cam : chr [1:390] "0-0.1" "0-0.1" "0.1-0.2" "0.1-0.2" ...
## ..$ rpt : int [1:390] 1 2 1 2 1 2 1 2 1 2 ...
## ..$ lightfree : num [1:390] 14.4 NA NA NA NA ...
## ..$ filter1 : num [1:390] 0.13 NA NA NA NA ...
## ..$ filterC1 : num [1:390] 0.274 NA NA NA NA ...
## ..$ lightocluded: num [1:390] 6.21 NA NA NA NA NA 3.1 NA NA NA ...
## ..$ filter2 : num [1:390] 0.128 NA NA NA NA ...
## ..$ filterC2 : num [1:390] 0.19 NA NA NA NA ...
## $ legend:List of 4
## ..$ abbrev : chr [1:40] "prod" "glic" "phos" "sulf" ...
## ..$ short_label : chr [1:40] "Soybean ~ productivity" "beta-glucosidase" "Acid ~ phosphatase" "Arylsulfatase" ...
## ..$ measure_unit: chr [1:40] "(t~ha^{-1})" "(mu*g~p-nitrophenol~g^{-1}~soil~h^{-1})" "(mu*g~p-nitrophenol~g^{-1}~soil~h^{-1})" "(mu*g~p-nitrophenol~g^{-1}~soil~h^{-1})" ...
## ..$ long_label : chr [1:40] "Soybean ~ productivity" "beta-glucosidase" "Acid ~ phosphatase" "Arylsulfatase" ...#-----------------------------------------------------------------------
# Tidy data to the exploratory data analysis.
# str(high_soybean)
swrc0 <- as_tibble(high_soybean$swrc)
str(swrc0)## Classes 'tbl_df', 'tbl' and 'data.frame': 390 obs. of 15 variables:
## $ faz : chr "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" ...
## $ tal : chr "AG2" "AG2" "AG2" "AG2" ...
## $ cam : chr "0-0.1" "0-0.1" "0.1-0.2" "0.1-0.2" ...
## $ rpt : int 1 2 1 2 1 2 1 2 1 2 ...
## $ S0 : num 0.441 0.445 0.39 0.381 0.421 0.437 0.579 0.518 0.477 0.402 ...
## $ S1 : num 0.425 0.418 0.35 0.367 0.396 0.41 0.527 0.501 0.431 0.398 ...
## $ S2 : num 0.402 0.404 0.347 0.358 0.387 0.408 0.457 0.459 0.422 0.389 ...
## $ S4 : num 0.376 0.384 0.34 0.349 0.379 0.403 0.407 0.422 0.395 0.378 ...
## $ S6 : num 0.366 0.368 0.335 0.34 0.369 0.394 0.389 0.406 0.384 0.366 ...
## $ S8 : num 0.357 0.361 0.328 0.331 0.359 0.381 0.376 0.4 0.376 0.359 ...
## $ S10 : num 0.354 0.356 0.321 0.325 0.352 0.37 0.355 0.395 0.371 0.352 ...
## $ S33 : num 0.314 0.319 0.286 0.283 0.294 0.29 0.314 0.359 0.331 0.302 ...
## $ S66 : num 0.296 0.314 0.283 0.276 0.287 0.282 0.291 0.304 0.319 0.284 ...
## $ S300 : num 0.27 0.287 0.253 0.253 0.261 0.259 0.272 0.297 0.284 0.269 ...
## $ S1500: num 0.223 0.23 0.225 0.222 0.237 0.242 0.261 0.245 0.284 0.264 ...# Soil bulk density (`dens`) and total porosity (`poro_e`) are needed to
# calculate variables.
poro <- as_tibble(high_soybean$atrib) %>%
select(faz, tal, cam, rpt, dens, poro_e) %>%
rename(ds = dens, poros = poro_e)
str(poro)## Classes 'tbl_df', 'tbl' and 'data.frame': 394 obs. of 6 variables:
## $ faz : chr "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" ...
## $ tal : chr "AG2" "AG2" "AG2" "AG2" ...
## $ cam : chr "0-0.1" "0-0.1" "0.1-0.2" "0.1-0.2" ...
## $ rpt : int 1 2 1 2 1 2 1 2 1 2 ...
## $ ds : num 1.23 1.19 1.28 1.3 1.23 1.2 1.01 1.08 1.15 1.24 ...
## $ poros: num 0.53 0.54 0.53 0.52 0.5 0.51 0.61 0.59 0.57 0.53 ...swrc0 <- inner_join(poro, swrc0)
# Creates an ID.
swrc0 <- swrc0 %>%
mutate(id = gsub(x = paste(faz, tal, sep = "_"),
pattern = "[^[:alnum:]_]",
replacement = ""),
cam = fct_recode(cam,
"I" = "0-0.1",
"II" = "0.1-0.2",
"III" = "0.2-0.4"),
idc = paste(id, cam))
# From wide to long format.
swrc1 <- swrc0 %>%
gather(key = "tens",
value = "umid",
grep("^S", names(swrc0), value = TRUE))
# Creates a numeric variable for tension and recode levels of layer.
swrc1 <- swrc1 %>%
mutate(tens = as.integer(sub("^S", "", tens)))
str(swrc1)## Classes 'tbl_df', 'tbl' and 'data.frame': 4290 obs. of 10 variables:
## $ faz : chr "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" "FAZ. ARGEMIRA" ...
## $ tal : chr "AG2" "AG2" "AG2" "AG2" ...
## $ cam : Factor w/ 3 levels "I","II","III": 1 1 2 2 3 3 1 1 2 2 ...
## $ rpt : int 1 2 1 2 1 2 1 2 1 2 ...
## $ ds : num 1.23 1.19 1.28 1.3 1.23 1.2 1.01 1.08 1.15 1.24 ...
## $ poros: num 0.53 0.54 0.53 0.52 0.5 0.51 0.61 0.59 0.57 0.53 ...
## $ id : chr "FAZARGEMIRA_AG2" "FAZARGEMIRA_AG2" "FAZARGEMIRA_AG2" "FAZARGEMIRA_AG2" ...
## $ idc : chr "FAZARGEMIRA_AG2 I" "FAZARGEMIRA_AG2 I" "FAZARGEMIRA_AG2 II" "FAZARGEMIRA_AG2 II" ...
## $ tens : int 0 0 0 0 0 0 0 0 0 0 ...
## $ umid : num 0.441 0.445 0.39 0.381 0.421 0.437 0.579 0.518 0.477 0.402 ...# Remove cases where umid is a non valid value.
swrc1 <- swrc1 %>%
filter(!is.na(umid),
between(umid, left = 0, right = 0.9))
# Replace 0 tension by a small value 0.1.
swrc1$tens[swrc1$tens == 0] <- 0.1
# ID with few tension.
id_rem <- swrc1 %>%
group_by(idc) %>%
summarize(u = length(unique(tens))) %>%
filter(u < 9) %>%
.$idc
unique(id_rem)## character(0)# Remove id with few values of tens x umid pairs.
swrc1 <- swrc1 %>%
filter(!(idc %in% id_rem))
#-----------------------------------------------------------------------
# Aggregate data using mean of replicates on each soil layer.
swrc2 <- swrc1 %>%
group_by(id, cam, tens) %>%
summarize(umid = mean(umid, na.rm = TRUE)) %>%
mutate(idc = paste(id, cam))
swrc2## # A tibble: 2,079 x 5
## # Groups: id, cam [189]
## id cam tens umid idc
## <chr> <fct> <dbl> <dbl> <chr>
## 1 FAZARGEMIRA_AG10 I 0.100 0.368 FAZARGEMIRA_AG10 I
## 2 FAZARGEMIRA_AG10 I 1.00 0.358 FAZARGEMIRA_AG10 I
## 3 FAZARGEMIRA_AG10 I 2.00 0.353 FAZARGEMIRA_AG10 I
## 4 FAZARGEMIRA_AG10 I 4.00 0.342 FAZARGEMIRA_AG10 I
## 5 FAZARGEMIRA_AG10 I 6.00 0.332 FAZARGEMIRA_AG10 I
## 6 FAZARGEMIRA_AG10 I 8.00 0.326 FAZARGEMIRA_AG10 I
## 7 FAZARGEMIRA_AG10 I 10.0 0.323 FAZARGEMIRA_AG10 I
## 8 FAZARGEMIRA_AG10 I 33.0 0.260 FAZARGEMIRA_AG10 I
## 9 FAZARGEMIRA_AG10 I 66.0 0.254 FAZARGEMIRA_AG10 I
## 10 FAZARGEMIRA_AG10 I 300. 0.238 FAZARGEMIRA_AG10 I
## # ... with 2,069 more rows#-----------------------------------------------------------------------
# Perform exploratory data analysis.
swrc0 %>% count(faz)## # A tibble: 14 x 2
## faz n
## <chr> <int>
## 1 FAZ. ARGEMIRA 48
## 2 FAZ. CALIFORNIA 6
## 3 FAZ. FLORIDA 6
## 4 FAZ. GMC 12
## 5 FAZ. JOTA BASSO 102
## 6 FAZ. LIBERDADE 6
## 7 FAZ. MONICA 30
## 8 FAZ. PORTO SEGURO 30
## 9 FAZ. SANGA DO ENGENHO 18
## 10 FAZ. SANTA MARIA 36
## 11 FAZ. SANTA TEREZINHA 18
## 12 FAZ. SÃO MIGUEL 24
## 13 FAZ. SUCURI 42
## 14 FAZ. ZANCANARO 12print(swrc0 %>% count(faz, tal), n = Inf)## # A tibble: 65 x 3
## faz tal n
## <chr> <chr> <int>
## 1 FAZ. ARGEMIRA AG10 6
## 2 FAZ. ARGEMIRA AG2 6
## 3 FAZ. ARGEMIRA AG3 6
## 4 FAZ. ARGEMIRA AG5 6
## 5 FAZ. ARGEMIRA AG6 6
## 6 FAZ. ARGEMIRA AG9 6
## 7 FAZ. ARGEMIRA STI5 6
## 8 FAZ. ARGEMIRA STI6 6
## 9 FAZ. CALIFORNIA CA4 6
## 10 FAZ. FLORIDA F1 6
## 11 FAZ. GMC MC1 6
## 12 FAZ. GMC MC2 6
## 13 FAZ. JOTA BASSO V10 6
## 14 FAZ. JOTA BASSO V16 6
## 15 FAZ. JOTA BASSO V17 6
## 16 FAZ. JOTA BASSO V18 6
## 17 FAZ. JOTA BASSO V19 6
## 18 FAZ. JOTA BASSO V2 6
## 19 FAZ. JOTA BASSO V20 6
## 20 FAZ. JOTA BASSO V21 7
## 21 FAZ. JOTA BASSO V22 5
## 22 FAZ. JOTA BASSO V23 6
## 23 FAZ. JOTA BASSO V24 6
## 24 FAZ. JOTA BASSO V3 6
## 25 FAZ. JOTA BASSO V33 6
## 26 FAZ. JOTA BASSO V4 6
## 27 FAZ. JOTA BASSO V5 6
## 28 FAZ. JOTA BASSO V6 6
## 29 FAZ. JOTA BASSO V7 6
## 30 FAZ. LIBERDADE L33A 6
## 31 FAZ. MONICA M14 6
## 32 FAZ. MONICA M15 6
## 33 FAZ. MONICA M25 6
## 34 FAZ. MONICA M28 6
## 35 FAZ. MONICA M30 6
## 36 FAZ. PORTO SEGURO PS1 6
## 37 FAZ. PORTO SEGURO PS3 6
## 38 FAZ. PORTO SEGURO PSP 6
## 39 FAZ. PORTO SEGURO PSP4 6
## 40 FAZ. PORTO SEGURO Q1 6
## 41 FAZ. SANGA DO ENGENHO SE7 6
## 42 FAZ. SANGA DO ENGENHO SE8 6
## 43 FAZ. SANGA DO ENGENHO SE9 6
## 44 FAZ. SANTA MARIA SM10 6
## 45 FAZ. SANTA MARIA SM31 6
## 46 FAZ. SANTA MARIA SM33 6
## 47 FAZ. SANTA MARIA SM35 6
## 48 FAZ. SANTA MARIA SM43 6
## 49 FAZ. SANTA MARIA SM49 6
## 50 FAZ. SANTA TEREZINHA ST310 6
## 51 FAZ. SANTA TEREZINHA STC1 6
## 52 FAZ. SANTA TEREZINHA STD1 6
## 53 FAZ. SÃO MIGUEL SMG1 6
## 54 FAZ. SÃO MIGUEL SMG13 6
## 55 FAZ. SÃO MIGUEL SMG4 6
## 56 FAZ. SÃO MIGUEL SMG5 6
## 57 FAZ. SUCURI SU1 6
## 58 FAZ. SUCURI SU12 6
## 59 FAZ. SUCURI SU19 6
## 60 FAZ. SUCURI SU20 6
## 61 FAZ. SUCURI SU6 6
## 62 FAZ. SUCURI SU8 6
## 63 FAZ. SUCURI SU9 6
## 64 FAZ. ZANCANARO ZN01 6
## 65 FAZ. ZANCANARO ZN10 6ggplot(data = swrc1,
mapping = aes(x = ds, y = poros)) +
geom_point(pch = 1) +
facet_wrap(facets = ~cam)
ggplot(data = swrc2,
mapping = aes(x = tens, y = umid, color = cam)) +
geom_point(pch = 1) +
geom_line() +
# geom_jitter() +
# geom_smooth(se = FALSE) +
# scale_x_log10() +
scale_x_log10(
breaks = scales::trans_breaks("log10",
function(x) 10^x),
labels = scales::trans_format("log10",
scales::math_format(10^.x))
) +
# annotation_logticks() +
facet_wrap(facets = ~id, ncol = 5) +
theme(strip.text.x = element_text(size = 5),
legend.position = "top")
2 Fitting SWRC model
#-----------------------------------------------------------------------
# Fitting one curve to all data to get good start values for van
# Genuchten SWRC model.
# Creates the log10 version of tens.
swrc2 <- swrc2 %>% mutate(ltens = log10(tens))
# The van Genuchten most used SWRC model.
vangen <- umid ~ thr + (ths - thr)/(1 + exp(tha + ltens)^thn)^(1 - 1/thn)
# Optimizing the model parameters.
n0 <- nls(formula = vangen,
data = swrc2,
start = list(thr = 0.2,
ths = 0.7,
tha = 1,
thn = 1.6))
coef(n0)## thr ths tha thn
## 0.1847449 0.3866983 -0.5447355 2.0363967plot(umid ~ ltens, data = swrc2)
with(as.list(coef(n0)), {
curve(thr + (ths - thr)/(1 + exp(tha + ltens)^thn)^(1 - 1/thn),
xname = "ltens",
add = TRUE)
})
ATTENTION
- The SWRC can be fitted for each site x layer combination.
- However, the estimated parameters should be aggregated to the productivity resolution level, that is, a value for each site.
- There are, at least, two ways of aggregating:
- Get an average of estimated parameters by layer in each site, that is, fit separed curves for each layer in a site then get the mean (or median) for each SWRC parameter estimate.
- Fit the curve for all data of a site, regardless of layer, and use this.
- The difference should be small, but as a mean of a nonlinear function is different of a nonlinear function of the means, they aren’t be the same.
- Which one is the most appropriate, I can’t say because I think it is a domain specific question not a statistical one.
- So, I chose fit curves with all data available in a site, regardless of layers.
2.1 Fitting for each site \(\times\) layer
#-----------------------------------------------------------------------
# Fitting for all idc.
vangen_g <-
umid ~ thr + (ths - thr)/(1 + exp(tha + ltens)^thn)^(1 - 1/thn) | idc
# Start values.
start <- coef(n0)
n1 <- nlsList(model = vangen_g,
start = start,
data = swrc2)
p <- plot(intervals(n1),
layout = c(4, NA))
update(p,
xlab = "Estimates",
ylab = "Experimental unit")
# Cases with non convergence.
i <- rownames(coef(n1))[!complete.cases(coef(n1))]
if (length(i)) {
# These did not converged.
ggplot(data = swrc2 %>% filter(idc %in% i),
mapping = aes(x = tens, y = umid)) +
geom_point(pch = 1) +
geom_line() +
scale_x_log10() +
facet_wrap(facets = ~idc)
}
#-----------------------------------------------------------------------
# Calculate S and I indexes.
estimates <- coef(n1)
params <- cbind(idc = rownames(estimates),
as.data.frame(estimates))
rownames(params) <- NULL
params <- as_tibble(params)
# Keep only complete cases.
params <- params %>% filter(complete.cases(.)) %>% droplevels()
# S and I indexes.
params <- within(params, {
S <- -(ths - thr) * thn * (1 + 1/(1 - 1/thn))^(-(1 - 1/thn) - 1)
I <- -tha - log(1 - 1/thn)/thn
U <- thr + (ths - thr)/(1 + exp(tha + I)^thn)^(1 - 1/thn)
cad <- U - thr
})
print(params, n = Inf)## # A tibble: 182 x 9
## idc thr ths tha thn cad U I S
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FAZARGEMIRA_AG10 I 0.226 0.363 -1.05 3.19 0.0740 0.300 1.17 -0.0960
## 2 FAZARGEMIRA_AG10 II 0.211 0.352 -0.878 2.23 0.0796 0.291 1.14 -0.0632
## 3 FAZARGEMIRA_AG10 III 0.215 0.413 -0.589 2.11 0.113 0.328 0.894 -0.0820
## 4 FAZARGEMIRA_AG2 I 0.107 0.464 -0.171 1.35 0.237 0.344 1.18 -0.0651
## 5 FAZARGEMIRA_AG2 II 0.168 0.390 -0.569 1.52 0.139 0.307 1.28 -0.0537
## 6 FAZARGEMIRA_AG2 III 0.240 0.419 -1.02 2.77 0.0982 0.338 1.18 -0.106
## 7 FAZARGEMIRA_AG3 I 0.248 0.565 -0.255 2.09 0.181 0.429 0.565 -0.130
## 8 FAZARGEMIRA_AG3 II 0.262 0.437 -0.693 2.27 0.0987 0.361 0.950 -0.0803
## 9 FAZARGEMIRA_AG3 III 0.273 0.492 -0.517 1.94 0.127 0.400 0.889 -0.0808
## 10 FAZARGEMIRA_AG5 I 0.199 0.476 0.162 1.66 0.168 0.367 0.394 -0.0795
## 11 FAZARGEMIRA_AG5 II 0.204 0.342 -0.485 1.66 0.0837 0.288 1.04 -0.0396
## 12 FAZARGEMIRA_AG5 III 0.213 0.380 -0.964 2.07 0.0960 0.309 1.28 -0.0679
## 13 FAZARGEMIRA_AG6 I 0.238 0.392 -0.940 3.11 0.0834 0.322 1.06 -0.105
## 14 FAZARGEMIRA_AG6 II 0.222 0.359 -0.918 2.44 0.0769 0.298 1.13 -0.0694
## 15 FAZARGEMIRA_AG6 III 0.222 0.421 -0.469 2.08 0.114 0.336 0.783 -0.0811
## 16 FAZARGEMIRA_AG9 I 0.208 0.429 -0.283 1.84 0.130 0.338 0.707 -0.0752
## 17 FAZARGEMIRA_AG9 II 0.224 0.355 -0.713 2.41 0.0731 0.297 0.936 -0.0651
## 18 FAZARGEMIRA_AG9 III 0.209 0.426 -0.251 2.02 0.125 0.334 0.589 -0.0849
## 19 FAZARGEMIRA_STI5 I 0.274 0.725 -0.253 2.52 0.250 0.524 0.454 -0.237
## 20 FAZARGEMIRA_STI5 III 0.254 0.440 -0.557 2.48 0.103 0.357 0.764 -0.0959
## 21 FAZARGEMIRA_STI6 I 0.115 0.656 0.193 1.48 0.343 0.457 0.562 -0.125
## 22 FAZARGEMIRA_STI6 III 0.230 0.466 -0.724 2.48 0.131 0.361 0.931 -0.122
## 23 FAZCALIFORNIA_CA4 I 0.128 0.335 -0.739 2.17 0.118 0.246 1.02 -0.0892
## 24 FAZCALIFORNIA_CA4 II 0.123 0.304 -0.374 1.80 0.107 0.231 0.827 -0.0591
## 25 FAZCALIFORNIA_CA4 III 0.160 0.328 -0.820 2.53 0.0932 0.253 1.02 -0.0887
## 26 FAZFLORIDA_F1 I 0.157 0.322 -0.886 2.57 0.0911 0.248 1.08 -0.0890
## 27 FAZFLORIDA_F1 II 0.114 0.238 -1.18 2.29 0.0702 0.184 1.43 -0.0581
## 28 FAZFLORIDA_F1 III 0.127 0.231 -0.841 2.44 0.0577 0.185 1.06 -0.0523
## 29 FAZGMC_MC1 I 0.221 0.592 -0.315 2.32 0.208 0.429 0.559 -0.175
## 30 FAZGMC_MC1 II 0.172 0.325 -1.57 4.21 0.0811 0.253 1.63 -0.147
## 31 FAZGMC_MC1 III 0.192 0.342 -0.736 1.87 0.0878 0.280 1.14 -0.0523
## 32 FAZGMC_MC2 I 0.214 0.490 -0.286 2.24 0.156 0.370 0.550 -0.125
## 33 FAZGMC_MC2 II 0.211 0.335 -0.886 3.06 0.0668 0.278 1.01 -0.0825
## 34 FAZGMC_MC2 III 0.207 0.408 -0.710 2.64 0.111 0.318 0.891 -0.112
## 35 FAZJOTABASSO_V10 I 0.238 0.488 -0.524 2.09 0.143 0.382 0.836 -0.102
## 36 FAZJOTABASSO_V10 II 0.205 0.375 -0.449 1.89 0.1000 0.305 0.849 -0.0603
## 37 FAZJOTABASSO_V10 III 0.193 0.390 -0.821 1.94 0.114 0.307 1.19 -0.0722
## 38 FAZJOTABASSO_V16 I 0.0882 0.420 -0.861 1.40 0.216 0.304 1.75 -0.0674
## 39 FAZJOTABASSO_V16 II 0.212 0.363 -0.542 2.06 0.0866 0.298 0.865 -0.0605
## 40 FAZJOTABASSO_V16 III 0.165 0.413 -0.662 1.81 0.146 0.312 1.10 -0.0824
## 41 FAZJOTABASSO_V17 I 0.236 0.509 -0.239 2.44 0.152 0.388 0.455 -0.138
## 42 FAZJOTABASSO_V17 II 0.219 0.405 -0.743 2.49 0.103 0.322 0.950 -0.0963
## 43 FAZJOTABASSO_V17 III 0.234 0.399 -0.529 2.16 0.0936 0.328 0.819 -0.0704
## 44 FAZJOTABASSO_V18 I -0.0316 0.476 0.0411 1.22 0.361 0.330 1.35 -0.0682
## 45 FAZJOTABASSO_V18 II 0.206 0.385 -0.439 1.99 0.104 0.310 0.789 -0.0686
## 46 FAZJOTABASSO_V18 III 0.211 0.380 -0.473 2.12 0.0962 0.308 0.774 -0.0704
## 47 FAZJOTABASSO_V19 I 0.229 0.384 -0.179 2.26 0.0876 0.316 0.439 -0.0707
## 48 FAZJOTABASSO_V19 II 0.198 0.399 -0.635 2.02 0.116 0.314 0.975 -0.0783
## 49 FAZJOTABASSO_V19 III 0.227 0.433 -0.570 2.85 0.113 0.339 0.722 -0.126
## 50 FAZJOTABASSO_V20 I 0.187 0.325 -0.888 2.23 0.0777 0.265 1.16 -0.0614
## 51 FAZJOTABASSO_V20 II 0.162 0.292 -0.843 1.78 0.0773 0.239 1.31 -0.0418
## 52 FAZJOTABASSO_V20 III 0.112 0.385 -0.433 1.51 0.172 0.284 1.15 -0.0655
## 53 FAZJOTABASSO_V21 I 0.0769 0.392 0.294 1.29 0.215 0.292 0.865 -0.0508
## 54 FAZJOTABASSO_V21 II 0.0790 0.343 -0.451 1.35 0.176 0.255 1.45 -0.0486
## 55 FAZJOTABASSO_V21 III 0.163 0.397 -0.166 1.65 0.142 0.306 0.732 -0.0662
## 56 FAZJOTABASSO_V22 I 0.188 0.471 -0.341 2.01 0.163 0.351 0.682 -0.110
## 57 FAZJOTABASSO_V22 II 0.180 0.305 -0.994 2.21 0.0705 0.251 1.27 -0.0553
## 58 FAZJOTABASSO_V22 III 0.119 0.380 -0.118 1.40 0.170 0.288 1.01 -0.0532
## 59 FAZJOTABASSO_V23 I 0.172 0.421 -0.529 1.71 0.149 0.322 1.04 -0.0750
## 60 FAZJOTABASSO_V23 II 0.209 0.364 -0.655 2.17 0.0876 0.297 0.941 -0.0664
## 61 FAZJOTABASSO_V23 III 0.193 0.403 -0.860 2.05 0.121 0.314 1.19 -0.0835
## 62 FAZJOTABASSO_V24 I 0.147 0.438 -0.439 1.66 0.177 0.323 0.991 -0.0839
## 63 FAZJOTABASSO_V24 II -0.0263 0.331 -0.258 1.18 0.263 0.236 1.85 -0.0409
## 64 FAZJOTABASSO_V24 III 0.201 0.381 -0.799 2.64 0.0992 0.300 0.979 -0.100
## 65 FAZJOTABASSO_V2 II 0.210 0.371 -0.395 1.82 0.0950 0.305 0.830 -0.0539
## 66 FAZJOTABASSO_V2 III 0.192 0.434 -0.415 1.71 0.145 0.337 0.929 -0.0730
## 67 FAZJOTABASSO_V33 I 0.156 0.474 -0.0622 1.91 0.185 0.342 0.451 -0.114
## 68 FAZJOTABASSO_V33 II 0.126 0.306 -0.792 1.95 0.105 0.231 1.16 -0.0665
## 69 FAZJOTABASSO_V33 III 0.148 0.369 -0.176 1.89 0.129 0.277 0.573 -0.0784
## 70 FAZJOTABASSO_V3 I 0.142 0.540 0.859 1.35 0.264 0.406 0.147 -0.0729
## 71 FAZJOTABASSO_V3 III 0.130 0.493 0.654 1.36 0.240 0.370 0.319 -0.0685
## 72 FAZJOTABASSO_V4 I -0.0994 0.453 -0.606 1.22 0.394 0.295 2.02 -0.0727
## 73 FAZJOTABASSO_V4 II 0.126 0.430 0.379 1.33 0.203 0.329 0.661 -0.0542
## 74 FAZJOTABASSO_V4 III 0.208 0.424 -0.786 2.02 0.125 0.332 1.12 -0.0845
## 75 FAZJOTABASSO_V5 I 0.119 0.425 -0.118 1.38 0.201 0.319 1.06 -0.0595
## 76 FAZJOTABASSO_V5 II 0.203 0.367 -0.736 2.40 0.0920 0.295 0.961 -0.0811
## 77 FAZJOTABASSO_V5 III 0.204 0.414 -0.806 2.40 0.117 0.321 1.03 -0.104
## 78 FAZJOTABASSO_V6 I 0.234 0.495 -0.501 2.54 0.145 0.379 0.697 -0.139
## 79 FAZJOTABASSO_V6 II 0.200 0.434 -0.238 1.81 0.139 0.339 0.682 -0.0776
## 80 FAZJOTABASSO_V6 III 0.216 0.419 -0.642 2.37 0.113 0.329 0.874 -0.0982
## 81 FAZJOTABASSO_V7 I 0.200 0.440 -0.178 2.16 0.137 0.336 0.466 -0.103
## 82 FAZJOTABASSO_V7 II 0.197 0.388 -0.582 2.01 0.110 0.307 0.925 -0.0741
## 83 FAZJOTABASSO_V7 III 0.205 0.398 -0.717 1.96 0.112 0.317 1.08 -0.0724
## 84 FAZLIBERDADE_L33A I 0.213 0.488 -0.706 2.87 0.150 0.363 0.855 -0.170
## 85 FAZMONICA_M14 I 0.217 0.577 -0.309 2.35 0.201 0.419 0.544 -0.173
## 86 FAZMONICA_M14 II 0.203 0.371 -0.951 2.79 0.0922 0.295 1.11 -0.101
## 87 FAZMONICA_M14 III 0.203 0.412 -0.730 2.39 0.116 0.320 0.957 -0.102
## 88 FAZMONICA_M15 II 0.169 0.534 0.607 1.55 0.227 0.396 0.0612 -0.0920
## 89 FAZMONICA_M15 III 0.122 0.234 -0.522 2.58 0.0616 0.184 0.712 -0.0602
## 90 FAZMONICA_M25 I 0.191 0.342 -1.07 2.13 0.0860 0.277 1.36 -0.0636
## 91 FAZMONICA_M25 II 0.206 0.299 -1.05 2.44 0.0520 0.258 1.26 -0.0471
## 92 FAZMONICA_M25 III 0.190 0.374 -0.623 1.80 0.109 0.299 1.07 -0.0605
## 93 FAZMONICA_M28 I 0.222 0.499 -0.341 2.33 0.155 0.377 0.582 -0.131
## 94 FAZMONICA_M28 III 0.118 0.195 -0.877 3.13 0.0414 0.159 1.000 -0.0526
## 95 FAZMONICA_M30 I 0.195 0.518 -0.357 1.82 0.191 0.386 0.797 -0.108
## 96 FAZMONICA_M30 II 0.198 0.420 -0.652 1.84 0.131 0.329 1.08 -0.0758
## 97 FAZMONICA_M30 III 0.216 0.395 -0.598 2.41 0.0998 0.316 0.821 -0.0886
## 98 FAZPORTOSEGURO_PS1 I 0.193 0.382 -0.387 1.88 0.110 0.304 0.791 -0.0661
## 99 FAZPORTOSEGURO_PS1 II 0.137 0.290 -0.455 1.69 0.0925 0.229 0.985 -0.0453
## 100 FAZPORTOSEGURO_PS1 III 0.189 0.329 -0.784 2.49 0.0780 0.267 0.991 -0.0727
## 101 FAZPORTOSEGURO_PS3 I 0.245 0.452 -0.715 2.93 0.113 0.358 0.857 -0.132
## 102 FAZPORTOSEGURO_PS3 II 0.216 0.363 -0.766 2.49 0.0817 0.297 0.972 -0.0763
## 103 FAZPORTOSEGURO_PS3 III 0.207 0.367 -0.763 2.48 0.0891 0.296 0.971 -0.0827
## 104 FAZPORTOSEGURO_PSP4 I 0.167 0.396 -0.168 1.66 0.139 0.306 0.727 -0.0654
## 105 FAZPORTOSEGURO_PSP4 II 0.127 0.332 -0.702 1.88 0.120 0.247 1.11 -0.0718
## 106 FAZPORTOSEGURO_PSP4 III 0.174 0.331 -0.930 1.94 0.0911 0.266 1.30 -0.0579
## 107 FAZPORTOSEGURO_PSP I 0.0681 0.141 -0.483 2.17 0.0413 0.109 0.769 -0.0313
## 108 FAZPORTOSEGURO_PSP II 0.128 0.233 -0.739 2.25 0.0593 0.187 1.00 -0.0476
## 109 FAZPORTOSEGURO_PSP III 0.134 0.234 -0.585 2.21 0.0568 0.191 0.857 -0.0444
## 110 FAZPORTOSEGURO_Q1 I 0.167 0.367 -0.624 2.12 0.114 0.281 0.926 -0.0833
## 111 FAZPORTOSEGURO_Q1 II 0.133 0.278 -0.866 1.61 0.0888 0.222 1.47 -0.0395
## 112 FAZPORTOSEGURO_Q1 III 0.185 0.307 -1.06 2.81 0.0666 0.251 1.21 -0.0735
## 113 FAZSANGADOENGENHO_SE7 I 0.212 0.478 -0.770 2.44 0.148 0.360 0.985 -0.134
## 114 FAZSANGADOENGENHO_SE7 III 0.206 0.425 -0.937 2.40 0.122 0.328 1.16 -0.108
## 115 FAZSANGADOENGENHO_SE8 I 0.161 0.385 -0.00223 2.10 0.128 0.289 0.312 -0.0923
## 116 FAZSANGADOENGENHO_SE8 II 0.186 0.376 -0.660 1.78 0.113 0.299 1.13 -0.0610
## 117 FAZSANGADOENGENHO_SE8 III 0.229 0.406 -0.809 2.96 0.0962 0.326 0.948 -0.114
## 118 FAZSANGADOENGENHO_SE9 I -0.646 0.428 -0.195 1.07 0.899 0.253 2.78 -0.0572
## 119 FAZSANGADOENGENHO_SE9 II 0.149 0.360 -1.14 2.41 0.118 0.267 1.37 -0.105
## 120 FAZSANGADOENGENHO_SE9 III 0.223 0.392 -0.569 2.49 0.0942 0.317 0.775 -0.0879
## 121 FAZSANTAMARIA_SM10 I 0.222 0.607 -0.371 2.41 0.215 0.437 0.595 -0.191
## 122 FAZSANTAMARIA_SM10 II 0.226 0.420 -0.701 2.22 0.109 0.336 0.969 -0.0863
## 123 FAZSANTAMARIA_SM10 III 0.225 0.387 -0.763 2.65 0.0892 0.315 0.942 -0.0908
## 124 FAZSANTAMARIA_SM31 II 0.223 0.346 -1.10 2.67 0.0677 0.291 1.27 -0.0695
## 125 FAZSANTAMARIA_SM31 III 0.221 0.340 -0.968 3.07 0.0644 0.285 1.10 -0.0798
## 126 FAZSANTAMARIA_SM33 I 0.209 0.426 -0.662 2.56 0.120 0.329 0.856 -0.117
## 127 FAZSANTAMARIA_SM33 III 0.201 0.366 -0.731 2.65 0.0906 0.292 0.910 -0.0920
## 128 FAZSANTAMARIA_SM35 I 0.231 0.484 -0.563 2.47 0.141 0.372 0.773 -0.130
## 129 FAZSANTAMARIA_SM35 II 0.171 0.333 -0.874 2.52 0.0902 0.261 1.07 -0.0854
## 130 FAZSANTAMARIA_SM43 I 0.227 0.436 -0.724 2.93 0.114 0.341 0.866 -0.133
## 131 FAZSANTAMARIA_SM43 II 0.195 0.356 -0.962 2.49 0.0893 0.285 1.17 -0.0831
## 132 FAZSANTAMARIA_SM43 III 0.223 0.414 -0.880 2.58 0.106 0.329 1.07 -0.103
## 133 FAZSANTAMARIA_SM49 I 0.259 0.625 -0.522 4.04 0.194 0.453 0.592 -0.336
## 134 FAZSANTAMARIA_SM49 II -0.0118 0.399 -0.330 1.22 0.293 0.281 1.73 -0.0549
## 135 FAZSANTAMARIA_SM49 III 0.194 0.374 -0.733 2.41 0.101 0.295 0.955 -0.0895
## 136 FAZSANTATEREZINHA_ST310 I 0.179 0.503 -0.316 1.86 0.191 0.369 0.731 -0.112
## 137 FAZSANTATEREZINHA_ST310 II 0.206 0.343 -1.14 3.16 0.0739 0.280 1.26 -0.0949
## 138 FAZSANTATEREZINHA_ST310 III 0.201 0.370 -0.834 2.69 0.0931 0.294 1.01 -0.0968
## 139 FAZSANTATEREZINHA_STC1 I -0.0299 0.478 -0.581 1.29 0.347 0.317 1.73 -0.0827
## 140 FAZSANTATEREZINHA_STC1 II 0.196 0.355 -0.734 1.73 0.0955 0.291 1.23 -0.0491
## 141 FAZSANTATEREZINHA_STC1 III 0.204 0.383 -0.971 2.14 0.102 0.306 1.26 -0.0761
## 142 FAZSANTATEREZINHA_STD1 I 0.218 0.421 -0.593 2.55 0.113 0.330 0.788 -0.109
## 143 FAZSANTATEREZINHA_STD1 II 0.206 0.340 -0.865 2.32 0.0753 0.281 1.11 -0.0631
## 144 FAZSANTATEREZINHA_STD1 III 0.210 0.365 -0.633 2.33 0.0870 0.297 0.875 -0.0735
## 145 FAZSÃOMIGUEL_SMG13 I 0.211 0.443 -0.541 2.52 0.129 0.340 0.742 -0.122
## 146 FAZSÃOMIGUEL_SMG13 II 0.203 0.358 -1.14 3.62 0.0829 0.286 1.23 -0.126
## 147 FAZSÃOMIGUEL_SMG13 III 0.199 0.415 -0.569 2.51 0.120 0.319 0.772 -0.113
## 148 FAZSÃOMIGUEL_SMG1 I 0.188 0.432 -0.656 2.34 0.137 0.324 0.894 -0.117
## 149 FAZSÃOMIGUEL_SMG1 II 0.183 0.351 -0.794 2.41 0.0939 0.277 1.02 -0.0838
## 150 FAZSÃOMIGUEL_SMG1 III 0.188 0.342 -0.779 2.96 0.0838 0.271 0.919 -0.0986
## 151 FAZSÃOMIGUEL_SMG4 I 0.196 0.402 -0.559 2.88 0.112 0.308 0.707 -0.128
## 152 FAZSÃOMIGUEL_SMG4 II 0.0537 0.321 -0.808 1.43 0.172 0.226 1.65 -0.0568
## 153 FAZSÃOMIGUEL_SMG4 III 0.175 0.320 -0.505 1.95 0.0841 0.259 0.872 -0.0539
## 154 FAZSÃOMIGUEL_SMG5 I 0.194 0.560 -0.325 2.49 0.203 0.397 0.531 -0.189
## 155 FAZSÃOMIGUEL_SMG5 II 0.132 0.392 -0.650 1.91 0.152 0.284 1.04 -0.0934
## 156 FAZSÃOMIGUEL_SMG5 III 0.195 0.429 -0.591 2.99 0.127 0.322 0.727 -0.152
## 157 FAZSUCURI_SU12 I 0.202 0.349 -0.920 3.42 0.0790 0.281 1.02 -0.112
## 158 FAZSUCURI_SU12 II 0.188 0.316 -1.01 2.52 0.0712 0.259 1.22 -0.0674
## 159 FAZSUCURI_SU12 III 0.173 0.411 -0.224 1.73 0.142 0.315 0.720 -0.0735
## 160 FAZSUCURI_SU19 I 0.0978 0.336 -0.287 3.23 0.128 0.226 0.401 -0.170
## 161 FAZSUCURI_SU19 II 0.0589 0.211 -0.580 2.59 0.0842 0.143 0.768 -0.0830
## 162 FAZSUCURI_SU19 III 0.0624 0.232 -0.569 3.65 0.0902 0.153 0.657 -0.139
## 163 FAZSUCURI_SU1 I 0.162 0.337 -1.11 1.78 0.104 0.266 1.57 -0.0564
## 164 FAZSUCURI_SU1 II 0.144 0.344 -0.359 1.49 0.127 0.270 1.11 -0.0464
## 165 FAZSUCURI_SU1 III 0.179 0.385 -0.345 1.82 0.121 0.300 0.782 -0.0687
## 166 FAZSUCURI_SU20 I 0.105 0.360 -0.335 2.34 0.143 0.248 0.573 -0.122
## 167 FAZSUCURI_SU20 II 0.0724 0.215 -0.302 1.89 0.0832 0.156 0.701 -0.0503
## 168 FAZSUCURI_SU20 III 0.0699 0.237 0.124 1.67 0.101 0.171 0.427 -0.0481
## 169 FAZSUCURI_SU6 II 0.202 0.307 -0.897 2.39 0.0583 0.261 1.12 -0.0513
## 170 FAZSUCURI_SU6 III 0.173 0.418 0.203 1.61 0.150 0.323 0.399 -0.0665
## 171 FAZSUCURI_SU8 I 0.130 0.355 -0.251 2.51 0.124 0.255 0.453 -0.117
## 172 FAZSUCURI_SU8 II 0.105 0.256 -0.265 1.95 0.0878 0.193 0.634 -0.0560
## 173 FAZSUCURI_SU8 III 0.0939 0.317 -0.136 1.68 0.135 0.229 0.674 -0.0654
## 174 FAZSUCURI_SU9 I 0.0876 0.259 -0.535 2.48 0.0953 0.183 0.743 -0.0883
## 175 FAZSUCURI_SU9 II 0.0799 0.213 -0.624 2.36 0.0747 0.155 0.858 -0.0644
## 176 FAZSUCURI_SU9 III 0.0874 0.269 -0.549 2.35 0.102 0.189 0.784 -0.0875
## 177 FAZZANCANARO_ZN01 I 0.0432 0.493 -0.519 1.39 0.294 0.337 1.43 -0.0897
## 178 FAZZANCANARO_ZN01 II 0.184 0.444 0.154 1.69 0.157 0.341 0.374 -0.0772
## 179 FAZZANCANARO_ZN01 III 0.180 0.353 -0.232 2.39 0.0967 0.277 0.459 -0.0847
## 180 FAZZANCANARO_ZN10 I 0.122 0.404 -0.508 2.27 0.159 0.281 0.763 -0.130
## 181 FAZZANCANARO_ZN10 II 0.113 0.238 -1.16 2.80 0.0682 0.182 1.32 -0.0747
## 182 FAZZANCANARO_ZN10 III 0.140 0.314 -0.845 2.77 0.0953 0.236 1.01 -0.103xyplot(umid ~ ltens | as.factor(idc),
data = swrc2 %>% filter(idc %in% levels(params$idc)),
as.table = TRUE,
strip = FALSE,
xlab = expression(Log[e] ~ "tension (kPa)"),
ylab = expression("Soil moisture" ~ (dm^{3} ~ dm^{-3}))) +
latticeExtra::layer({
p <- as.list(params[packet.number(), ])
with(p, {
# Curva ajustada.
panel.curve(
thr + (ths - thr)/(1 + exp(tha + x)^thn)^(1 - 1/thn))
# Tensão, umidade e inclinação no ponto de inflexão.
panel.points(x = I, y = U, pch = 4)
panel.segments(I, 0, I, U, lty = 2)
panel.segments(-5, U, I, U, lty = 2)
panel.abline(a = U - S * I, b = S, lty = 3)
})
})
# Keep in an object with descriptive name.
params_idc <- paramssplom(params[, -1])
2.2 Fitting for each site
#-----------------------------------------------------------------------
# Fitting for all id.
vangen_g <-
umid ~ thr + (ths - thr)/(1 + exp(tha + ltens)^thn)^(1 - 1/thn) | id
# Start values.
start <- coef(n0)
n1 <- nlsList(model = vangen_g,
start = start,
data = swrc2)
p <- plot(intervals(n1),
layout = c(4, NA))
update(p,
xlab = "Estimativa",
ylab = "Unidade experimental")
# Cases with non convergence.
i <- rownames(coef(n1))[!complete.cases(coef(n1))]
if (length(i)) {
# These did not converged.
ggplot(data = swrc2 %>% filter(id %in% i),
mapping = aes(x = tens, y = umid)) +
geom_point(pch = 1) +
stat_summary(mapping = aes(group = 1),
fun.y = mean,
colour = "tomato",
geom = "line",
group = 1) +
scale_x_log10() +
facet_wrap(facets = ~id)
}#-----------------------------------------------------------------------
# Calculate S and I indexes.
estimates <- coef(n1)
# Transforma de matriz para tabela.
params <- cbind(id = rownames(estimates),
as.data.frame(estimates))
rownames(params) <- NULL
params <- as_tibble(params)
# Keep only complete cases.
params <- params %>% filter(complete.cases(.)) %>% droplevels()
# S and I indexes.
params <- within(params, {
S <- -(ths - thr) * thn * (1 + 1/(1 - 1/thn))^(-(1 - 1/thn) - 1)
I <- -tha - log(1 - 1/thn)/thn
U <- thr + (ths - thr)/(1 + exp(tha + I)^thn)^(1 - 1/thn)
cad <- U - thr
})
print(params, n = Inf)## # A tibble: 65 x 9
## id thr ths tha thn cad U I S
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FAZARGEMIRA_AG10 0.217 0.376 -0.832 2.36 0.0889 0.306 1.07 -0.0765
## 2 FAZARGEMIRA_AG2 0.199 0.423 -0.678 1.74 0.134 0.332 1.17 -0.0698
## 3 FAZARGEMIRA_AG3 0.259 0.499 -0.414 2.02 0.139 0.397 0.752 -0.0938
## 4 FAZARGEMIRA_AG5 0.201 0.396 -0.375 1.70 0.118 0.318 0.895 -0.0585
## 5 FAZARGEMIRA_AG6 0.227 0.391 -0.763 2.41 0.0912 0.319 0.985 -0.0814
## 6 FAZARGEMIRA_AG9 0.213 0.404 -0.364 1.99 0.111 0.323 0.717 -0.0729
## 7 FAZARGEMIRA_STI5 0.241 0.552 -0.257 1.99 0.180 0.421 0.606 -0.119
## 8 FAZARGEMIRA_STI6 0.143 0.496 -0.289 1.63 0.216 0.358 0.875 -0.0975
## 9 FAZCALIFORNIA_CA4 0.140 0.320 -0.717 2.20 0.102 0.242 0.993 -0.0792
## 10 FAZFLORIDA_F1 0.133 0.264 -0.950 2.41 0.0730 0.206 1.17 -0.0650
## 11 FAZGMC_MC1 0.189 0.426 -0.545 2.04 0.137 0.325 0.876 -0.0940
## 12 FAZGMC_MC2 0.210 0.408 -0.596 2.49 0.110 0.320 0.803 -0.103
## 13 FAZJOTABASSO_V10 0.210 0.418 -0.574 1.94 0.121 0.331 0.949 -0.0763
## 14 FAZJOTABASSO_V16 0.176 0.398 -0.667 1.72 0.133 0.310 1.17 -0.0678
## 15 FAZJOTABASSO_V17 0.228 0.438 -0.437 2.27 0.118 0.346 0.691 -0.0963
## 16 FAZJOTABASSO_V18 0.192 0.412 -0.364 1.74 0.131 0.323 0.857 -0.0682
## 17 FAZJOTABASSO_V19 0.220 0.404 -0.492 2.36 0.103 0.323 0.726 -0.0888
## 18 FAZJOTABASSO_V2 0.152 0.405 -0.159 1.45 0.162 0.314 0.965 -0.0556
## 19 FAZJOTABASSO_V20 0.159 0.334 -0.704 1.76 0.104 0.263 1.18 -0.0551
## 20 FAZJOTABASSO_V21 0.122 0.377 -0.103 1.43 0.164 0.286 0.948 -0.0541
## 21 FAZJOTABASSO_V22 0.168 0.385 -0.424 1.76 0.129 0.297 0.901 -0.0685
## 22 FAZJOTABASSO_V23 0.194 0.396 -0.683 1.94 0.117 0.311 1.06 -0.0742
## 23 FAZJOTABASSO_V24 0.172 0.382 -0.602 1.85 0.124 0.295 1.02 -0.0718
## 24 FAZJOTABASSO_V3 0.120 0.468 0.612 1.32 0.234 0.355 0.474 -0.0595
## 25 FAZJOTABASSO_V33 0.139 0.383 -0.231 1.83 0.144 0.283 0.665 -0.0821
## 26 FAZJOTABASSO_V4 0.138 0.432 -0.434 1.47 0.187 0.325 1.22 -0.0662
## 27 FAZJOTABASSO_V5 0.196 0.399 -0.658 2.02 0.117 0.313 0.995 -0.0794
## 28 FAZJOTABASSO_V6 0.220 0.448 -0.492 2.26 0.129 0.349 0.752 -0.104
## 29 FAZJOTABASSO_V7 0.196 0.410 -0.397 1.93 0.124 0.321 0.776 -0.0780
## 30 FAZLIBERDADE_L33A 0.213 0.488 -0.706 2.87 0.150 0.363 0.855 -0.170
## 31 FAZMONICA_M14 0.202 0.456 -0.518 2.23 0.144 0.346 0.784 -0.114
## 32 FAZMONICA_M15 0.151 0.370 0.00697 1.80 0.130 0.280 0.445 -0.0717
## 33 FAZMONICA_M25 0.197 0.337 -0.916 2.06 0.0807 0.277 1.24 -0.0564
## 34 FAZMONICA_M28 0.169 0.345 -0.487 2.41 0.0985 0.267 0.708 -0.0878
## 35 FAZMONICA_M30 0.206 0.444 -0.512 1.96 0.138 0.344 0.876 -0.0892
## 36 FAZPORTOSEGURO_PS1 0.175 0.332 -0.578 1.99 0.0909 0.266 0.928 -0.0601
## 37 FAZPORTOSEGURO_PS3 0.223 0.394 -0.741 2.67 0.0941 0.317 0.918 -0.0965
## 38 FAZPORTOSEGURO_PSP 0.110 0.202 -0.627 2.23 0.0521 0.162 0.895 -0.0412
## 39 FAZPORTOSEGURO_PSP4 0.149 0.353 -0.542 1.72 0.122 0.272 1.05 -0.0621
## 40 FAZPORTOSEGURO_Q1 0.167 0.317 -0.821 2.12 0.0853 0.253 1.12 -0.0625
## 41 FAZSANGADOENGENHO_SE7 0.172 0.456 -0.832 1.92 0.166 0.337 1.22 -0.103
## 42 FAZSANGADOENGENHO_SE8 0.201 0.380 -0.596 2.41 0.100 0.301 0.819 -0.0888
## 43 FAZSANGADOENGENHO_SE9 0.185 0.392 -0.747 1.91 0.121 0.306 1.13 -0.0748
## 44 FAZSANTAMARIA_SM10 0.223 0.472 -0.518 2.32 0.140 0.363 0.763 -0.117
## 45 FAZSANTAMARIA_SM31 0.222 0.343 -1.02 2.84 0.0659 0.288 1.18 -0.0737
## 46 FAZSANTAMARIA_SM33 0.188 0.368 -0.853 2.61 0.0994 0.287 1.04 -0.0992
## 47 FAZSANTAMARIA_SM35 0.200 0.409 -0.672 2.44 0.117 0.316 0.888 -0.105
## 48 FAZSANTAMARIA_SM43 0.215 0.403 -0.825 2.59 0.104 0.318 1.01 -0.103
## 49 FAZSANTAMARIA_SM49 0.221 0.468 -0.513 2.55 0.137 0.357 0.709 -0.132
## 50 FAZSANTATEREZINHA_ST310 0.196 0.401 -0.754 2.30 0.115 0.311 1.00 -0.0955
## 51 FAZSANTATEREZINHA_STC1 0.166 0.404 -0.798 1.66 0.144 0.311 1.36 -0.0678
## 52 FAZSANTATEREZINHA_STD1 0.210 0.376 -0.661 2.36 0.0924 0.303 0.893 -0.0799
## 53 FAZSÃOMIGUEL_SMG1 0.186 0.375 -0.731 2.50 0.105 0.291 0.936 -0.0982
## 54 FAZSÃOMIGUEL_SMG13 0.201 0.406 -0.705 2.53 0.114 0.314 0.904 -0.108
## 55 FAZSÃOMIGUEL_SMG4 0.171 0.347 -0.594 2.13 0.100 0.271 0.892 -0.0741
## 56 FAZSÃOMIGUEL_SMG5 0.178 0.461 -0.465 2.38 0.158 0.336 0.694 -0.138
## 57 FAZSUCURI_SU1 0.161 0.354 -0.564 1.65 0.117 0.279 1.13 -0.0546
## 58 FAZSUCURI_SU12 0.190 0.358 -0.717 2.25 0.0949 0.285 0.979 -0.0760
## 59 FAZSUCURI_SU19 0.0744 0.259 -0.448 3.14 0.0998 0.174 0.570 -0.127
## 60 FAZSUCURI_SU20 0.0834 0.270 -0.227 1.99 0.108 0.191 0.576 -0.0714
## 61 FAZSUCURI_SU6 0.172 0.347 -0.0650 1.56 0.108 0.281 0.726 -0.0444
## 62 FAZSUCURI_SU8 0.114 0.310 -0.206 2.03 0.113 0.227 0.540 -0.0770
## 63 FAZSUCURI_SU9 0.0853 0.247 -0.566 2.41 0.0903 0.176 0.789 -0.0803
## 64 FAZZANCANARO_ZN01 0.185 0.423 -0.314 1.89 0.139 0.324 0.715 -0.0838
## 65 FAZZANCANARO_ZN10 0.122 0.319 -0.746 2.37 0.110 0.232 0.977 -0.0955xyplot(umid ~ ltens | as.factor(id),
data = swrc2 %>% filter(id %in% levels(params$id)),
as.table = TRUE,
strip = FALSE,
xlab = expression(Log[e] ~ "tension (kPa)"),
ylab = expression("Soil moisture" ~ (dm^{3} ~ dm^{-3}))) +
latticeExtra::layer({
p <- as.list(params[packet.number(), ])
with(p, {
# Curva ajustada.
panel.curve(
thr + (ths - thr)/(1 + exp(tha + x)^thn)^(1 - 1/thn))
# Tensão, umidade e inclinação no ponto de inflexão.
panel.points(x = I, y = U, pch = 4)
panel.segments(I, 0, I, U, lty = 2)
panel.segments(-5, U, I, U, lty = 2)
panel.abline(a = U - S * I, b = S, lty = 3)
})
})
params_id <- paramssplom(params[, -1])
3 Regression for the soybean yield on each SWRC parameter
3.1 Calculate values based on soil density and porosity
# params_idc
params <-
swrc0 %>%
select(idc, id, cam, ds, poros) %>%
group_by(idc) %>%
summarize(id = id[1],
cam = cam[1],
ds = mean(ds, na.rm = TRUE),
poros = mean(poros, na.rm = TRUE)) %>%
inner_join(params_idc)
str(params)## Classes 'tbl_df', 'tbl' and 'data.frame': 182 obs. of 13 variables:
## $ idc : chr "FAZARGEMIRA_AG10 I" "FAZARGEMIRA_AG10 II" "FAZARGEMIRA_AG10 III" "FAZARGEMIRA_AG2 I" ...
## $ id : chr "FAZARGEMIRA_AG10" "FAZARGEMIRA_AG10" "FAZARGEMIRA_AG10" "FAZARGEMIRA_AG2" ...
## $ cam : Factor w/ 3 levels "I","II","III": 1 2 3 1 2 3 1 2 3 1 ...
## $ ds : num 1.36 1.35 1.22 1.21 1.29 ...
## $ poros: num 0.51 0.495 0.525 0.535 0.525 0.505 0.6 0.55 0.56 0.565 ...
## $ thr : num 0.226 0.211 0.215 0.107 0.168 ...
## $ ths : num 0.363 0.352 0.413 0.464 0.39 ...
## $ tha : num -1.048 -0.878 -0.589 -0.171 -0.569 ...
## $ thn : num 3.19 2.23 2.11 1.35 1.52 ...
## $ cad : num 0.074 0.0796 0.113 0.2369 0.1393 ...
## $ U : num 0.3 0.291 0.328 0.344 0.307 ...
## $ I : num 1.166 1.144 0.894 1.181 1.279 ...
## $ S : num -0.096 -0.0632 -0.082 -0.0651 -0.0537 ...params <- params %>%
mutate(fwc = U * ds,
macrop = poros - fwc,
paw = (U - thr) * ds * 10 * ifelse(cam == "III", 20, 10),
S = -1 * S)
# Table with productivity.
atrib <- as_tibble(high_soybean$atrib)
names(atrib)## [1] "faz" "tal" "cam" "rpt" "prod" "glic" "sulf"
## [8] "phos" "FLF" "OLF" "MO2" "HF" "rp" "dens"
## [15] "poro_e" "clay" "sand" "pH_H2O" "pH_CaCl2" "P" "K"
## [22] "Ca_Mg" "Ca" "Mg" "Al" "MO" "base" "CTC"
## [29] "V" "CaMg" "CaK" "MgK" "Ca_" "Mg_" "K_"
## [36] "H_" "m_" "S" "dp"# Select and filter.
atrib <- atrib %>%
select(faz, tal, cam, rpt, prod) %>%
mutate(id = gsub(x = paste(faz, tal, sep = "_"),
pattern = "[^[:alnum:]_]",
replacement = ""),
idc = paste(id, cam)) %>%
select(id, prod) %>%
filter(!is.na(prod))
str(atrib)## Classes 'tbl_df', 'tbl' and 'data.frame': 69 obs. of 2 variables:
## $ id : chr "FAZARGEMIRA_AG2" "FAZARGEMIRA_AG3" "FAZARGEMIRA_AG5" "FAZARGEMIRA_AG6" ...
## $ prod: num 3.89 4.12 4.25 4.12 4.38 ...# Add the productivity.
params <- inner_join(atrib, params)
# More descriptive labels.
params <- params %>%
mutate(cam = fct_recode(cam,
"L1" = "I",
"L2" = "II",
"L3" = "III"))
# Remove outliers.
i <- which(params$S > 0.3)
params$S[i] <- NA
i <- which(params$paw > 140)
params$paw[i] <- NA
ylab <- expression("Soybean productivity" ~ (t ~ ha^{-1}))
mp <- atrib %>% pull(prod) %>% mean(na.rm = TRUE)
mp## [1] 4.1293043.2 Exploratory data analysis
# Brings variables in another table.
# str(high_soybean$atrib)
phy <- high_soybean$atrib %>%
as_tibble() %>%
filter(faz != "") %>%
select(faz, tal, cam, prod,
rp, dp, poro_e, clay, sand, MO2, FLF, OLF, HF) %>%
group_by(faz, tal, cam) %>%
summarise_if(is.numeric, "mean", na.rm = TRUE) %>%
mutate(id = gsub(x = paste(faz, tal, sep = "_"),
pattern = "[^[:alnum:]_]",
replacement = ""),
cam = fct_recode(cam,
"L1" = "0-0.1",
"L2" = "0.1-0.2",
"L3" = "0.2-0.4"))
phy## # A tibble: 195 x 14
## # Groups: faz, tal [65]
## faz tal cam prod rp dp poro_e clay sand MO2 FLF OLF HF id
## <chr> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 FAZ.… AG10 L1 4.01 2.00 2.78 0.510 661. 211. 42.3 16.0 2.83 23.4 FAZA…
## 2 FAZ.… AG10 L2 NaN 2.40 2.67 0.495 722. 174. NaN NaN NaN NaN FAZA…
## 3 FAZ.… AG10 L3 NaN 0.915 2.57 0.525 749. 156. NaN NaN NaN NaN FAZA…
## 4 FAZ.… AG2 L1 3.89 1.92 2.59 0.535 653. 162. 45.5 14.4 6.21 24.8 FAZA…
## 5 FAZ.… AG2 L2 NaN 2.64 2.73 0.525 725. 139. NaN NaN NaN NaN FAZA…
## 6 FAZ.… AG2 L3 NaN 1.74 2.45 0.505 775. 120. NaN NaN NaN NaN FAZA…
## 7 FAZ.… AG3 L1 4.12 0.900 2.63 0.600 644. 175. 45.3 12.1 3.10 30.1 FAZA…
## 8 FAZ.… AG3 L2 NaN 1.80 2.65 0.550 724. 138. NaN NaN NaN NaN FAZA…
## 9 FAZ.… AG3 L3 NaN 1.21 2.59 0.560 773. 127. NaN NaN NaN NaN FAZA…
## 10 FAZ.… AG5 L1 4.25 1.66 2.68 0.565 713. 179. 38.7 16.0 6.86 15.8 FAZA…
## # ... with 185 more rows# Legend.
leg <- list("prod" = '"Soybean productivity"',
"rp" = '"Penetration resistence"',
"ds" = '"Bulk density"',
"dp" = '"Particle density"',
"poros" = '"Total porosity"',
"macrop" = '"Macroporosity"',
"fwc" = '"Field water capacity"',
"ths" = '"WC at saturation"',
"thr" = '"Residual WC"',
"I" = '"Inflexion point"',
"paw" = '"Plant available water"',
"thn" = '"SWRC" ~ n ~ "parameter"',
"tha" = '"SWRC" ~ alpha ~ "parameter"',
"S" = '"S index"',
# "cad" = '"Available WC"',
# "U" = '"WC at inflexion"',
"clay" = '"Clay"',
"sand" = '"Sand"',
"MO2" = '"SOM"',
"FLF" = '"FLF of SOM"',
"OLF" = '"OLF of SOM"',
"HF" = '"HF of SOM"')
db <- left_join(params[, -2], phy)
# str(db)
v <- names(leg)
dc <- db %>%
gather(key = "variable", value = "value", v) %>%
filter(is.finite(value)) %>%
mutate(variable = factor(variable, levels = v))
dc## # A tibble: 3,033 x 10
## id idc cam cad U faz tal poro_e variable value
## <chr> <chr> <fct> <dbl> <dbl> <chr> <chr> <dbl> <fct> <dbl>
## 1 FAZARGEMIRA_AG2 FAZARGEMIRA… L1 0.237 0.344 FAZ. ARG… AG2 0.535 prod 3.89
## 2 FAZARGEMIRA_AG3 FAZARGEMIRA… L1 0.181 0.429 FAZ. ARG… AG3 0.600 prod 4.12
## 3 FAZARGEMIRA_AG5 FAZARGEMIRA… L1 0.168 0.367 FAZ. ARG… AG5 0.565 prod 4.25
## 4 FAZARGEMIRA_AG6 FAZARGEMIRA… L1 0.0834 0.322 FAZ. ARG… AG6 0.540 prod 4.12
## 5 FAZARGEMIRA_AG9 FAZARGEMIRA… L1 0.130 0.338 FAZ. ARG… AG9 0.550 prod 4.38
## 6 FAZARGEMIRA_AG10 FAZARGEMIRA… L1 0.0740 0.300 FAZ. ARG… AG10 0.510 prod 4.01
## 7 FAZARGEMIRA_STI5 FAZARGEMIRA… L1 0.250 0.524 FAZ. ARG… STI5 0.835 prod 3.88
## 8 FAZARGEMIRA_STI6 FAZARGEMIRA… L1 0.343 0.457 FAZ. ARG… STI6 0.625 prod 4.14
## 9 FAZJOTABASSO_V3 FAZJOTABASS… L1 0.264 0.406 FAZ. JOT… V3 0.565 prod 4.23
## 10 FAZJOTABASSO_V4 FAZJOTABASS… L1 0.394 0.295 FAZ. JOT… V4 0.535 prod 4.04
## # ... with 3,023 more rows# A labeller function.
label_legend <- function(labels) {
legend <- attr(label_legend, "leg")
labels <- lapply(labels, as.character)
x <- lapply(unname(labels),
lapply,
function(values) {
values <- legend[[match(x = values,
table = names(legend))]]
# print(values)
c(parse(text = as.character(values)))
})
return(x)
}
class(label_legend) <- c("function", "labeller") # IMPORTANT.
attr(label_legend, "legend") <- leg
ggplot(data = dc,
mapping = aes(x = cam,
y = value)) +
geom_boxplot(fill = "gray90") +
stat_summary(mapping = aes(group = 1),
geom = "point",
fun.y = mean,
shape = 4) +
stat_summary(mapping = aes(group = 1),
geom = "line",
fun.y = mean) +
xlab(label = "Soil layer") +
ylab(label = "Values of each response variable") +
facet_wrap(facets = ~variable,
scales = "free",
ncol = 4,
labeller = labeller(variable = label_legend)) +
theme_bw()
dc_stats <- dc %>%
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))
kable(dc_stats)| variable | cam | n | mean | sd | min | q25 | median | q75 | max |
|---|---|---|---|---|---|---|---|---|---|
| prod | L1 | 61 | 4.1168033 | 0.3373218 | 3.3690000 | 3.9300000 | 4.1220000 | 4.3260000 | 4.9740000 |
| rp | L1 | 61 | 1.5119672 | 0.7156817 | 0.3100000 | 0.9750000 | 1.4150000 | 1.8550000 | 3.2750000 |
| rp | L2 | 58 | 2.0636207 | 0.5605165 | 0.9850000 | 1.7100000 | 1.9475000 | 2.3987500 | 3.9000000 |
| rp | L3 | 63 | 1.3938889 | 0.4038139 | 0.8550000 | 1.0875000 | 1.3300000 | 1.5150000 | 2.6900000 |
| ds | L1 | 61 | 1.2448361 | 0.1444298 | 0.9400000 | 1.1450000 | 1.2350000 | 1.3450000 | 1.6800000 |
| ds | L2 | 58 | 1.4005172 | 0.1243287 | 1.1950000 | 1.3212500 | 1.3550000 | 1.4437500 | 1.7000000 |
| ds | L3 | 63 | 1.3384127 | 0.1186885 | 1.1400000 | 1.2625000 | 1.3150000 | 1.3900000 | 1.7100000 |
| dp | L1 | 61 | 2.6686885 | 0.0717281 | 2.5300000 | 2.6100000 | 2.6700000 | 2.7100000 | 2.8300000 |
| dp | L2 | 58 | 2.6796552 | 0.0639618 | 2.5700000 | 2.6400000 | 2.6700000 | 2.7300000 | 2.8500000 |
| dp | L3 | 63 | 2.6853968 | 0.1011997 | 2.4500000 | 2.6300000 | 2.6700000 | 2.7300000 | 3.1200000 |
| poros | L1 | 61 | 0.5434973 | 0.0771235 | 0.3750000 | 0.5100000 | 0.5350000 | 0.5650000 | 0.8350000 |
| poros | L2 | 58 | 0.4990517 | 0.0870393 | 0.3700000 | 0.4600000 | 0.4900000 | 0.5100000 | 0.7650000 |
| poros | L3 | 63 | 0.5173810 | 0.0723280 | 0.3950000 | 0.4825000 | 0.5150000 | 0.5325000 | 0.7650000 |
| macrop | L1 | 61 | 0.1402805 | 0.0635960 | 0.0184709 | 0.1093753 | 0.1352847 | 0.1625038 | 0.3674941 |
| macrop | L2 | 58 | 0.1238411 | 0.0722021 | 0.0143186 | 0.0935764 | 0.1068904 | 0.1285664 | 0.3972585 |
| macrop | L3 | 63 | 0.1301759 | 0.0753043 | 0.0308192 | 0.0942957 | 0.1110611 | 0.1352419 | 0.3784978 |
| fwc | L1 | 61 | 0.4032167 | 0.0487485 | 0.1836859 | 0.3893453 | 0.4084403 | 0.4261199 | 0.5136955 |
| fwc | L2 | 58 | 0.3752106 | 0.0457329 | 0.2411629 | 0.3609451 | 0.3803481 | 0.4032940 | 0.4986296 |
| fwc | L3 | 63 | 0.3872050 | 0.0479468 | 0.2198367 | 0.3809953 | 0.3991494 | 0.4133674 | 0.4791808 |
| ths | L1 | 61 | 0.4446750 | 0.0964380 | 0.1406798 | 0.3854681 | 0.4381136 | 0.4899636 | 0.7253569 |
| ths | L2 | 58 | 0.3450356 | 0.0620951 | 0.2113470 | 0.3090173 | 0.3515361 | 0.3759856 | 0.5336300 |
| ths | L3 | 63 | 0.3733374 | 0.0635150 | 0.1945221 | 0.3417250 | 0.3853031 | 0.4140684 | 0.4931956 |
| thr | L1 | 61 | 0.1553391 | 0.1272870 | -0.6460182 | 0.1281111 | 0.1911035 | 0.2175876 | 0.2739949 |
| thr | L2 | 58 | 0.1654414 | 0.0588757 | -0.0262791 | 0.1290156 | 0.1868742 | 0.2057042 | 0.2623583 |
| thr | L3 | 63 | 0.1850131 | 0.0439502 | 0.0623564 | 0.1689792 | 0.1947094 | 0.2119418 | 0.2732351 |
| I | L1 | 61 | 0.8649532 | 0.4427001 | 0.1465555 | 0.5647225 | 0.7729930 | 1.0244808 | 2.7781596 |
| I | L2 | 58 | 1.0731385 | 0.3067811 | 0.0612062 | 0.9370119 | 1.0754989 | 1.2533480 | 1.8544357 |
| I | L3 | 63 | 0.9066542 | 0.2173034 | 0.3191510 | 0.7746363 | 0.9293324 | 1.0638919 | 1.3020688 |
| paw | L1 | 60 | 19.7375538 | 8.0094275 | 6.9307563 | 14.8262872 | 17.6268721 | 21.5781533 | 49.2942457 |
| paw | L2 | 58 | 14.7423353 | 6.1708437 | 7.4843783 | 10.9985724 | 12.9510854 | 15.4865102 | 38.3214888 |
| paw | L3 | 63 | 28.6642381 | 7.6736460 | 11.4345892 | 24.8153413 | 27.5614759 | 30.5117860 | 62.0686709 |
| thn | L1 | 61 | 2.1813122 | 0.6024322 | 1.0676900 | 1.7792934 | 2.2258327 | 2.5098557 | 4.0378574 |
| thn | L2 | 58 | 2.1454776 | 0.5598745 | 1.1794053 | 1.7815309 | 2.1127252 | 2.4305181 | 4.2054472 |
| thn | L3 | 63 | 2.2855702 | 0.4642661 | 1.3617810 | 1.9437218 | 2.3655854 | 2.5763630 | 3.6494193 |
| tha | L1 | 61 | -0.4278274 | 0.3488026 | -1.1103832 | -0.6242082 | -0.4393957 | -0.2526656 | 0.8591422 |
| tha | L2 | 58 | -0.6750400 | 0.3708154 | -1.5703628 | -0.8835285 | -0.7078087 | -0.4623732 | 0.6071776 |
| tha | L3 | 63 | -0.6005003 | 0.3124999 | -1.0553702 | -0.8077760 | -0.6417096 | -0.5112255 | 0.6542138 |
| S | L1 | 60 | 0.1067058 | 0.0398915 | 0.0312534 | 0.0744891 | 0.1063172 | 0.1298299 | 0.2369213 |
| S | L2 | 58 | 0.0693430 | 0.0209180 | 0.0394904 | 0.0543412 | 0.0664534 | 0.0809262 | 0.1474759 |
| S | L3 | 63 | 0.0834072 | 0.0219549 | 0.0444308 | 0.0686133 | 0.0820496 | 0.0974872 | 0.1517702 |
| clay | L1 | 61 | 591.0491803 | 124.8885405 | 151.0000000 | 562.0000000 | 614.0000000 | 669.0000000 | 752.0000000 |
| clay | L2 | 58 | 639.4310345 | 128.4912601 | 195.0000000 | 618.5000000 | 675.5000000 | 725.7500000 | 770.0000000 |
| clay | L3 | 63 | 679.9047619 | 123.7085172 | 215.0000000 | 659.0000000 | 710.0000000 | 762.0000000 | 798.0000000 |
| sand | L1 | 61 | 294.1803279 | 151.5610447 | 143.0000000 | 189.0000000 | 248.0000000 | 328.0000000 | 826.0000000 |
| sand | L2 | 58 | 269.6206897 | 148.4886229 | 124.0000000 | 173.2500000 | 234.0000000 | 280.7500000 | 786.0000000 |
| sand | L3 | 63 | 236.8412698 | 139.2885386 | 111.0000000 | 151.0000000 | 190.0000000 | 253.5000000 | 763.0000000 |
| MO2 | L1 | 61 | 35.8181967 | 5.2122665 | 20.9400000 | 32.6300000 | 36.7100000 | 38.6600000 | 45.8600000 |
| FLF | L1 | 61 | 15.2509836 | 3.9019566 | 2.2400000 | 12.7400000 | 15.3200000 | 17.7600000 | 24.5200000 |
| OLF | L1 | 61 | 5.9586885 | 2.9978717 | 1.6500000 | 4.2000000 | 5.2000000 | 6.4000000 | 14.8900000 |
| HF | L1 | 61 | 14.6085246 | 5.6791528 | 2.1000000 | 9.4200000 | 15.3900000 | 18.5000000 | 30.1000000 |
# Save the tidy dataset.
params## # A tibble: 182 x 17
## id prod idc cam ds poros thr ths tha thn cad U I S
## <chr> <dbl> <chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FAZA… 3.89 FAZA… L1 1.21 0.535 0.107 0.464 -0.171 1.35 0.237 0.344 1.18 0.0651
## 2 FAZA… 3.89 FAZA… L2 1.29 0.525 0.168 0.390 -0.569 1.52 0.139 0.307 1.28 0.0537
## 3 FAZA… 3.89 FAZA… L3 1.21 0.505 0.240 0.419 -1.02 2.77 0.0982 0.338 1.18 0.106
## 4 FAZA… 4.12 FAZA… L1 1.04 0.600 0.248 0.565 -0.255 2.09 0.181 0.429 0.565 0.130
## 5 FAZA… 4.12 FAZA… L2 1.19 0.550 0.262 0.437 -0.693 2.27 0.0987 0.361 0.950 0.0803
## 6 FAZA… 4.12 FAZA… L3 1.14 0.560 0.273 0.492 -0.517 1.94 0.127 0.400 0.889 0.0808
## 7 FAZA… 4.25 FAZA… L1 1.17 0.565 0.199 0.476 0.162 1.66 0.168 0.367 0.394 0.0795
## 8 FAZA… 4.25 FAZA… L2 1.36 0.490 0.204 0.342 -0.485 1.66 0.0837 0.288 1.04 0.0396
## 9 FAZA… 4.25 FAZA… L3 1.30 0.520 0.213 0.380 -0.964 2.07 0.0960 0.309 1.28 0.0679
## 10 FAZA… 4.12 FAZA… L1 1.28 0.540 0.238 0.392 -0.940 3.11 0.0834 0.322 1.06 0.105
## # ... with 172 more rows, and 3 more variables: fwc <dbl>, macrop <dbl>, paw <dbl>save(params, file = "swrc.RData")3.3 S index
xlab <- expression("S index")
params$x <- params$S
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 4.9633 -108.16
## poly(L1, degree = 2) 2 0.46854 5.4318 -107.47 2.1240 0.13138
## poly(L2, degree = 2) 2 0.58830 5.5516 -106.33 2.6669 0.08043 .
## poly(L3, degree = 2) 2 0.36157 5.3248 -108.50 1.6391 0.20554
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.66989 -0.17191 -0.02583 0.24382 0.58270
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11263 0.04605 89.299 <2e-16 ***
## poly(L1, degree = 2)1 0.72242 0.44892 1.609 0.1146
## poly(L1, degree = 2)2 -0.35454 0.39192 -0.905 0.3705
## poly(L2, degree = 2)1 -0.26955 0.44002 -0.613 0.5432
## poly(L2, degree = 2)2 -0.80020 0.38480 -2.080 0.0433 *
## poly(L3, degree = 2)1 -0.13940 0.40301 -0.346 0.7310
## poly(L3, degree = 2)2 0.70128 0.39784 1.763 0.0847 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3321 on 45 degrees of freedom
## Multiple R-squared: 0.2393, Adjusted R-squared: 0.1379
## F-statistic: 2.36 on 6 and 45 DF, p-value: 0.04569m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63551 -0.22381 -0.01233 0.16144 0.72027
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8934 0.2130 18.277 <2e-16 ***
## L1 4.0732 1.6729 2.435 0.0187 *
## L2 -3.8441 2.8686 -1.340 0.1865
## L3 0.8668 2.3685 0.366 0.7160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3454 on 48 degrees of freedom
## Multiple R-squared: 0.1222, Adjusted R-squared: 0.06729
## F-statistic: 2.226 on 3 and 48 DF, p-value: 0.09714m2 <- update(m0, . ~ L1, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L1, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.68595 -0.19640 -0.01904 0.19943 0.77336
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8156 0.1425 26.769 <2e-16 ***
## L1 2.9230 1.3212 2.212 0.0315 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3448 on 50 degrees of freedom
## Multiple R-squared: 0.08916, Adjusted R-squared: 0.07095
## F-statistic: 4.895 on 1 and 50 DF, p-value: 0.03154# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L1, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.4 Inflexion point (I index)
xlab <- expression("Inflexion point (kPa)")
params$x <- params$I
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5.6394 -104.75
## poly(L1, degree = 2) 2 0.63522 6.2746 -103.09 2.5907 0.08587 .
## poly(L2, degree = 2) 2 0.10701 5.7464 -107.75 0.4364 0.64898
## poly(L3, degree = 2) 2 0.07912 5.7185 -108.01 0.3227 0.72584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73389 -0.20509 -0.02254 0.16358 0.88670
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11787 0.04809 85.620 <2e-16 ***
## poly(L1, degree = 2)1 -0.80259 0.35699 -2.248 0.0294 *
## poly(L1, degree = 2)2 -0.09171 0.36064 -0.254 0.8004
## poly(L2, degree = 2)1 -0.28236 0.39033 -0.723 0.4731
## poly(L2, degree = 2)2 -0.17606 0.38059 -0.463 0.6458
## poly(L3, degree = 2)1 -0.27645 0.40228 -0.687 0.4954
## poly(L3, degree = 2)2 -0.11513 0.37906 -0.304 0.7627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3501 on 46 degrees of freedom
## Multiple R-squared: 0.1456, Adjusted R-squared: 0.03416
## F-statistic: 1.307 on 6 and 46 DF, p-value: 0.2734m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70840 -0.21195 -0.00454 0.17694 0.86175
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.6248 0.2589 17.861 <2e-16 ***
## L1 -0.2528 0.1047 -2.414 0.0195 *
## L2 -0.1310 0.1752 -0.748 0.4581
## L3 -0.1485 0.2483 -0.598 0.5525
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3409 on 49 degrees of freedom
## Multiple R-squared: 0.1372, Adjusted R-squared: 0.08438
## F-statistic: 2.597 on 3 and 49 DF, p-value: 0.06283m2 <- update(m0, . ~ L1, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L1, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70579 -0.17354 -0.01795 0.19745 0.79421
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3540 0.1034 42.101 <2e-16 ***
## L1 -0.2632 0.1029 -2.556 0.0136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3387 on 51 degrees of freedom
## Multiple R-squared: 0.1136, Adjusted R-squared: 0.09621
## F-statistic: 6.536 on 1 and 51 DF, p-value: 0.01359# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L1, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.5 Macroporosity
xlab <- expression("Macroporosity" ~ (m^3 ~ m^{-3}))
params$x <- params$macrop
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5.8066 -103.20
## poly(L1, degree = 2) 2 0.24834 6.0550 -104.98 0.9837 0.3817
## poly(L2, degree = 2) 2 0.29714 6.1038 -104.55 1.1770 0.3173
## poly(L3, degree = 2) 2 0.21619 6.0228 -105.26 0.8563 0.4314summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75477 -0.23942 0.02531 0.26291 0.73124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11787 0.04880 84.378 <2e-16 ***
## poly(L1, degree = 2)1 0.49949 0.37325 1.338 0.187
## poly(L1, degree = 2)2 -0.12470 0.36511 -0.342 0.734
## poly(L2, degree = 2)1 0.02624 0.38334 0.068 0.946
## poly(L2, degree = 2)2 -0.60010 0.39129 -1.534 0.132
## poly(L3, degree = 2)1 -0.33046 0.37300 -0.886 0.380
## poly(L3, degree = 2)2 0.37040 0.37952 0.976 0.334
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3553 on 46 degrees of freedom
## Multiple R-squared: 0.1203, Adjusted R-squared: 0.005514
## F-statistic: 1.048 on 6 and 46 DF, p-value: 0.4073m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73785 -0.28456 -0.00999 0.23129 0.79557
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.89560 0.19106 20.389 <2e-16 ***
## L1 2.32107 1.32192 1.756 0.0854 .
## L2 -0.03841 0.72692 -0.053 0.9581
## L3 -0.56119 0.69595 -0.806 0.4239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.354 on 49 degrees of freedom
## Multiple R-squared: 0.06962, Adjusted R-squared: 0.01266
## F-statistic: 1.222 on 3 and 49 DF, p-value: 0.3116m2 <- update(m0, . ~ L1, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L1, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72318 -0.27253 0.00019 0.24055 0.83112
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8322 0.1717 22.316 <2e-16 ***
## L1 2.2153 1.2786 1.733 0.0892 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3496 on 51 degrees of freedom
## Multiple R-squared: 0.05559, Adjusted R-squared: 0.03707
## F-statistic: 3.002 on 1 and 51 DF, p-value: 0.08921# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L1, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.6 Field water capacity (microporosity)
xlab <- expression("Field water capacity (microporosity)" ~ (m^3 ~ m^{-3}))
params$x <- params$fwc
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
dd## # A tibble: 53 x 5
## id prod L1 L2 L3
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 FAZARGEMIRA_AG10 4.01 0.408 0.391 0.400
## 2 FAZARGEMIRA_AG2 3.89 0.417 0.396 0.411
## 3 FAZARGEMIRA_AG3 4.12 0.449 0.431 0.456
## 4 FAZARGEMIRA_AG5 4.25 0.430 0.391 0.401
## 5 FAZARGEMIRA_AG6 4.12 0.412 0.397 0.430
## 6 FAZARGEMIRA_AG9 4.38 0.413 0.404 0.403
## 7 FAZCALIFORNIA_CA4 4.03 0.334 0.333 0.364
## 8 FAZFLORIDA_F1 3.83 0.343 0.298 0.317
## 9 FAZGMC_MC1 3.64 0.434 0.363 0.380
## 10 FAZGMC_MC2 3.75 0.437 0.398 0.399
## # ... with 43 more rowsm0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 6.4536 -97.60
## poly(L1, degree = 2) 2 0.035513 6.4891 -101.31 0.1266 0.8814
## poly(L2, degree = 2) 2 0.043406 6.4970 -101.24 0.1547 0.8571
## poly(L3, degree = 2) 2 0.044755 6.4984 -101.23 0.1595 0.8530summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71172 -0.26510 0.01399 0.18495 0.88770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11787 0.05145 80.036 <2e-16 ***
## poly(L1, degree = 2)1 0.24583 0.52892 0.465 0.644
## poly(L1, degree = 2)2 0.06136 0.41034 0.150 0.882
## poly(L2, degree = 2)1 -0.09763 1.20183 -0.081 0.936
## poly(L2, degree = 2)2 0.35732 0.82706 0.432 0.668
## poly(L3, degree = 2)1 -0.22859 1.36398 -0.168 0.868
## poly(L3, degree = 2)2 -0.20081 0.70600 -0.284 0.777
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3746 on 46 degrees of freedom
## Multiple R-squared: 0.02224, Adjusted R-squared: -0.1053
## F-statistic: 0.1744 on 6 and 46 DF, p-value: 0.9824m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71528 -0.25871 -0.00075 0.19474 0.89644
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4210 0.5286 8.364 5.34e-11 ***
## L1 0.8091 1.4601 0.554 0.582
## L2 0.8812 2.4175 0.364 0.717
## L3 -2.4575 3.0062 -0.817 0.418
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3637 on 49 degrees of freedom
## Multiple R-squared: 0.01779, Adjusted R-squared: -0.04235
## F-statistic: 0.2958 on 3 and 49 DF, p-value: 0.8282m2 <- update(m0, . ~ L3, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73325 -0.20804 0.00464 0.20936 0.86330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.4491 0.5085 8.750 9.97e-12 ***
## L3 -0.8548 1.3060 -0.655 0.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3582 on 51 degrees of freedom
## Multiple R-squared: 0.00833, Adjusted R-squared: -0.01111
## F-statistic: 0.4284 on 1 and 51 DF, p-value: 0.5157# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L3, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.7 Plant available water
xlab <- expression("Plant available water" ~ (mm))
params$x <- params$paw
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
scale_x_continuous(limits = c(NA, 75)) + # <-- outlier.
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
filter(x < 100) %>% # <-- outlier.
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5.7245 -100.74
## poly(L1, degree = 2) 2 0.042221 5.7668 -104.36 0.1659 0.8476
## poly(L2, degree = 2) 2 0.124504 5.8491 -103.62 0.4894 0.6162
## poly(L3, degree = 2) 2 0.148504 5.8731 -103.41 0.5837 0.5620summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74278 -0.19611 0.00482 0.22336 0.73964
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.13227 0.04946 83.546 <2e-16 ***
## poly(L1, degree = 2)1 -0.22084 0.40041 -0.552 0.584
## poly(L1, degree = 2)2 0.07715 0.36512 0.211 0.834
## poly(L2, degree = 2)1 0.20970 0.38773 0.541 0.591
## poly(L2, degree = 2)2 -0.32796 0.37219 -0.881 0.383
## poly(L3, degree = 2)1 -0.20137 0.35937 -0.560 0.578
## poly(L3, degree = 2)2 -0.33842 0.37003 -0.915 0.365
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3567 on 45 degrees of freedom
## Multiple R-squared: 0.05047, Adjusted R-squared: -0.07613
## F-statistic: 0.3987 on 6 and 45 DF, p-value: 0.876m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73744 -0.23211 0.00482 0.19905 0.83345
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.228327 0.286094 14.780 <2e-16 ***
## L1 -0.001071 0.006421 -0.167 0.868
## L2 0.004479 0.008611 0.520 0.605
## L3 -0.004941 0.008455 -0.584 0.562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.352 on 48 degrees of freedom
## Multiple R-squared: 0.01347, Adjusted R-squared: -0.04818
## F-statistic: 0.2185 on 3 and 48 DF, p-value: 0.8831m2 <- update(m0, . ~ L3, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75425 -0.22785 -0.00921 0.22403 0.85574
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.281802 0.241636 17.720 <2e-16 ***
## L3 -0.005228 0.008281 -0.631 0.531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3459 on 50 degrees of freedom
## Multiple R-squared: 0.00791, Adjusted R-squared: -0.01193
## F-statistic: 0.3987 on 1 and 50 DF, p-value: 0.5307# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L3, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.8 Water content at saturation (ths)
xlab <- expression("Water content at saturation" ~ (g ~ g^{-1}))
params$x <- params$ths
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 6.3056 -98.829
## poly(L1, degree = 2) 2 0.060762 6.3664 -102.321 0.2216 0.8021
## poly(L2, degree = 2) 2 0.062404 6.3680 -102.307 0.2276 0.7973
## poly(L3, degree = 2) 2 0.146976 6.4526 -101.608 0.5361 0.5886summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72715 -0.31769 0.02935 0.16558 0.84744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11787 0.05086 80.970 <2e-16 ***
## poly(L1, degree = 2)1 0.34916 0.52449 0.666 0.509
## poly(L1, degree = 2)2 -0.04444 0.41676 -0.107 0.916
## poly(L2, degree = 2)1 0.43014 0.86265 0.499 0.620
## poly(L2, degree = 2)2 0.08117 0.58679 0.138 0.891
## poly(L3, degree = 2)1 -0.73825 0.87461 -0.844 0.403
## poly(L3, degree = 2)2 -0.03891 0.53126 -0.073 0.942
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3702 on 46 degrees of freedom
## Multiple R-squared: 0.04467, Adjusted R-squared: -0.07994
## F-statistic: 0.3584 on 6 and 46 DF, p-value: 0.9013m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72787 -0.32079 0.03086 0.16740 0.84112
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.2286 0.3436 12.308 <2e-16 ***
## L1 0.5447 0.7758 0.702 0.486
## L2 1.1259 1.6738 0.673 0.504
## L3 -1.9637 1.5311 -1.282 0.206
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3588 on 49 degrees of freedom
## Multiple R-squared: 0.04403, Adjusted R-squared: -0.0145
## F-statistic: 0.7523 on 3 and 49 DF, p-value: 0.5263m2 <- update(m0, . ~ L3, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73809 -0.21217 0.01494 0.21539 0.88604
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3133 0.3352 12.868 <2e-16 ***
## L3 -0.5254 0.8915 -0.589 0.558
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3585 on 51 degrees of freedom
## Multiple R-squared: 0.006764, Adjusted R-squared: -0.01271
## F-statistic: 0.3473 on 1 and 51 DF, p-value: 0.5582# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L3, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.9 Residual water content (thr)
xlab <- expression("Residual water content" ~ (g ~ g^{-1}))
params$x <- params$thr
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
scale_x_continuous(limits = c(-0.25, NA)) + # <-- outlier.
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
filter(x > -0.25) %>% # <-- outlier.
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5.8384 -99.713
## poly(L1, degree = 2) 2 0.062361 5.9007 -103.161 0.2403 0.7874
## poly(L2, degree = 2) 2 0.086471 5.9248 -102.949 0.3332 0.7183
## poly(L3, degree = 2) 2 0.088375 5.9267 -102.932 0.3406 0.7132summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.66195 -0.22680 -0.02659 0.16745 0.79239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.13227 0.04995 82.728 <2e-16 ***
## poly(L1, degree = 2)1 0.21626 0.41028 0.527 0.601
## poly(L1, degree = 2)2 0.24401 0.43611 0.560 0.579
## poly(L2, degree = 2)1 -0.38644 0.48861 -0.791 0.433
## poly(L2, degree = 2)2 -0.25312 0.48384 -0.523 0.603
## poly(L3, degree = 2)1 0.08841 0.57493 0.154 0.878
## poly(L3, degree = 2)2 0.32353 0.39353 0.822 0.415
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3602 on 45 degrees of freedom
## Multiple R-squared: 0.0316, Adjusted R-squared: -0.09753
## F-statistic: 0.2447 on 6 and 45 DF, p-value: 0.959m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73335 -0.18101 0.00245 0.20478 0.81634
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.1551 0.2187 19.003 <2e-16 ***
## L1 0.2278 0.6916 0.329 0.743
## L2 -0.5516 0.9843 -0.560 0.578
## L3 0.1617 1.3936 0.116 0.908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3529 on 48 degrees of freedom
## Multiple R-squared: 0.008368, Adjusted R-squared: -0.05361
## F-statistic: 0.135 on 3 and 48 DF, p-value: 0.9387m2 <- update(m0, . ~ L2, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L2, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74585 -0.20532 0.01083 0.20822 0.82909
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.1997 0.1380 30.440 <2e-16 ***
## L2 -0.4137 0.7938 -0.521 0.605
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3463 on 50 degrees of freedom
## Multiple R-squared: 0.005402, Adjusted R-squared: -0.01449
## F-statistic: 0.2716 on 1 and 50 DF, p-value: 0.6046# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L2, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.10 \(\alpha\) empirical shape parameter (tha)
xlab <- expression("SWRC" ~ alpha ~ "parameter")
params$x <- params$tha
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5.6115 -105.01
## poly(L1, degree = 2) 2 0.31123 5.9227 -106.15 1.2756 0.2889
## poly(L2, degree = 2) 2 0.58338 6.1949 -103.77 2.3911 0.1028
## poly(L3, degree = 2) 2 0.23040 5.8419 -106.88 0.9444 0.3963summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76098 -0.15400 0.00765 0.17102 0.69786
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11787 0.04798 85.832 <2e-16 ***
## poly(L1, degree = 2)1 0.06537 0.36222 0.180 0.8576
## poly(L1, degree = 2)2 -0.57036 0.36050 -1.582 0.1205
## poly(L2, degree = 2)1 0.22628 0.35861 0.631 0.5312
## poly(L2, degree = 2)2 -0.74944 0.35667 -2.101 0.0411 *
## poly(L3, degree = 2)1 0.21216 0.36098 0.588 0.5596
## poly(L3, degree = 2)2 -0.44049 0.35896 -1.227 0.2260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3493 on 46 degrees of freedom
## Multiple R-squared: 0.1498, Adjusted R-squared: 0.03894
## F-statistic: 1.351 on 6 and 46 DF, p-value: 0.2545m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74723 -0.22501 0.02776 0.21500 0.85079
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.22516 0.16880 25.030 <2e-16 ***
## L1 0.04166 0.16602 0.251 0.803
## L2 0.07526 0.15399 0.489 0.627
## L3 0.05904 0.20008 0.295 0.769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3652 on 49 degrees of freedom
## Multiple R-squared: 0.009811, Adjusted R-squared: -0.05081
## F-statistic: 0.1618 on 3 and 49 DF, p-value: 0.9215m2 <- update(m0, . ~ L3, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75288 -0.21762 -0.00384 0.21808 0.85378
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.16528 0.13063 31.887 <2e-16 ***
## L3 0.07625 0.19450 0.392 0.697
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3592 on 51 degrees of freedom
## Multiple R-squared: 0.003005, Adjusted R-squared: -0.01654
## F-statistic: 0.1537 on 1 and 51 DF, p-value: 0.6967# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L3, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
3.11 \(n\) empirical shape parameter (thn)
xlab <- expression("SWRC" ~ n ~ "parameter")
params$x <- params$thn
ggplot(data = params,
mapping = aes(x = x, y = prod)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
dd <- params %>%
select(id, prod, cam, x) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# dd
m0 <- lm(prod ~
poly(L1, degree = 2) +
poly(L2, degree = 2) +
poly(L3, degree = 2),
data = dd)
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) + poly(L3,
## degree = 2)
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 5.7309 -103.89
## poly(L1, degree = 2) 2 0.58865 6.3196 -102.71 2.3624 0.1055
## poly(L2, degree = 2) 2 0.29838 6.0293 -105.20 1.1975 0.3112
## poly(L3, degree = 2) 2 0.04803 5.7789 -107.45 0.1928 0.8253summary(m0)##
## Call:
## lm(formula = prod ~ poly(L1, degree = 2) + poly(L2, degree = 2) +
## poly(L3, degree = 2), data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.66114 -0.18346 -0.03934 0.17005 0.69694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.11787 0.04848 84.933 <2e-16 ***
## poly(L1, degree = 2)1 0.76680 0.37331 2.054 0.0457 *
## poly(L1, degree = 2)2 -0.35082 0.38372 -0.914 0.3654
## poly(L2, degree = 2)1 -0.49921 0.38795 -1.287 0.2046
## poly(L2, degree = 2)2 -0.26645 0.36416 -0.732 0.4681
## poly(L3, degree = 2)1 0.19766 0.36026 0.549 0.5859
## poly(L3, degree = 2)2 -0.09940 0.35829 -0.277 0.7827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.353 on 46 degrees of freedom
## Multiple R-squared: 0.1317, Adjusted R-squared: 0.01848
## F-statistic: 1.163 on 6 and 46 DF, p-value: 0.3423m1 <- update(m0, . ~ L1 + L2 + L3)
summary(m1)##
## Call:
## lm(formula = prod ~ L1 + L2 + L3, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63706 -0.26232 -0.01741 0.20970 0.73247
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.79842 0.31342 12.119 2.35e-16 ***
## L1 0.16759 0.08145 2.058 0.045 *
## L2 -0.09305 0.08732 -1.066 0.292
## L3 0.06814 0.11172 0.610 0.545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3485 on 49 degrees of freedom
## Multiple R-squared: 0.09813, Adjusted R-squared: 0.04291
## F-statistic: 1.777 on 3 and 49 DF, p-value: 0.1638m2 <- update(m0, . ~ L1, data = dd)
summary(m2)##
## Call:
## lm(formula = prod ~ L1, data = dd)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.62444 -0.26996 -0.02548 0.20765 0.80644
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.77981 0.17609 21.465 <2e-16 ***
## L1 0.15575 0.07811 1.994 0.0515 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3465 on 51 degrees of freedom
## Multiple R-squared: 0.07232, Adjusted R-squared: 0.05413
## F-statistic: 3.976 on 1 and 51 DF, p-value: 0.05152# par(mfrow = c(2, 2))
# plot(m2)
# layout(1)
gg <- ggplot(data = dd,
mapping = aes(x = L1, y = prod)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = xlab) +
ylab(label = ylab)gg
4 Regression for the soybean yield on all SWRC parameters
# str(params)
dd <- params %>%
select(id, prod, cam, ths, thr, S, I, fwc, paw, macrop) %>%
gather(key = "var", value = "value", ths:macrop) %>%
mutate(varcam = paste0(var, cam)) %>%
select(-cam, -var) %>%
spread(key = varcam, value = value)
str(dd)## Classes 'tbl_df', 'tbl' and 'data.frame': 65 obs. of 23 variables:
## $ id : chr "FAZARGEMIRA_AG10" "FAZARGEMIRA_AG2" "FAZARGEMIRA_AG3" "FAZARGEMIRA_AG5" ...
## $ prod : num 4.01 3.89 4.12 4.25 4.12 ...
## $ fwcL1 : num 0.408 0.417 0.449 0.43 0.412 ...
## $ fwcL2 : num 0.391 0.396 0.431 0.391 0.397 ...
## $ fwcL3 : num 0.4 0.411 0.456 0.401 0.43 ...
## $ IL1 : num 1.166 1.181 0.565 0.394 1.064 ...
## $ IL2 : num 1.14 1.28 0.95 1.04 1.13 ...
## $ IL3 : num 0.894 1.181 0.889 1.281 0.783 ...
## $ macropL1: num 0.102 0.118 0.151 0.135 0.128 ...
## $ macropL2: num 0.1041 0.1291 0.1186 0.0986 0.0932 ...
## $ macropL3: num 0.1251 0.0943 0.1036 0.1189 0.0849 ...
## $ pawL1 : num 10.1 28.7 19 19.7 10.7 ...
## $ pawL2 : num 10.7 18 11.8 11.4 10.2 ...
## $ pawL3 : num 27.6 23.9 29 25 29.1 ...
## $ SL1 : num 0.096 0.0651 0.1303 0.0795 0.1049 ...
## $ SL2 : num 0.0632 0.0537 0.0803 0.0396 0.0694 ...
## $ SL3 : num 0.082 0.106 0.0808 0.0679 0.0811 ...
## $ thrL1 : num 0.226 0.107 0.248 0.199 0.238 ...
## $ thrL2 : num 0.211 0.168 0.262 0.204 0.222 ...
## $ thrL3 : num 0.215 0.24 0.273 0.213 0.222 ...
## $ thsL1 : num 0.363 0.464 0.565 0.476 0.392 ...
## $ thsL2 : num 0.352 0.39 0.437 0.342 0.359 ...
## $ thsL3 : num 0.413 0.419 0.492 0.38 0.421 ...dd %>% complete.cases(.) %>% table()## .
## FALSE TRUE
## 14 51# m0 <- lm(prod ~ . -id, data = dd, na.action = na.omit)
m0 <- lm(prod ~ . -id, data = na.omit(dd))
# Type III hypothesis test.
drop1(m0, test = "F")## Single term deletions
##
## Model:
## prod ~ (id + fwcL1 + fwcL2 + fwcL3 + IL1 + IL2 + IL3 + macropL1 +
## macropL2 + macropL3 + pawL1 + pawL2 + pawL3 + SL1 + SL2 +
## SL3 + thrL1 + thrL2 + thrL3 + thsL1 + thsL2 + thsL3) - id
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 4.2622 -82.583
## fwcL1 1 0.083592 4.3458 -83.593 0.5688 0.4568
## fwcL2 1 0.157758 4.4200 -82.730 1.0734 0.3087
## fwcL3 1 0.261170 4.5234 -81.550 1.7770 0.1929
## IL1 1 0.213802 4.4760 -82.087 1.4547 0.2375
## IL2 1 0.000405 4.2627 -84.579 0.0028 0.9585
## IL3 1 0.037268 4.2995 -84.140 0.2536 0.6184
## macropL1 1 0.012500 4.2747 -84.434 0.0850 0.7726
## macropL2 1 0.194894 4.4571 -82.303 1.3260 0.2589
## macropL3 1 0.002255 4.2645 -84.557 0.0153 0.9023
## pawL1 1 0.079629 4.3419 -83.639 0.5418 0.4676
## pawL2 1 0.006579 4.2688 -84.505 0.0448 0.8339
## pawL3 1 0.157233 4.4195 -82.736 1.0698 0.3095
## SL1 1 0.112720 4.3750 -83.252 0.7669 0.3884
## SL2 1 0.115215 4.3775 -83.223 0.7839 0.3832
## SL3 1 0.047081 4.3093 -84.023 0.3203 0.5758
## thrL1 1 0.102241 4.3645 -83.375 0.6956 0.4111
## thrL2 1 0.000541 4.2628 -84.577 0.0037 0.9521
## thrL3 1 0.145544 4.4078 -82.871 0.9903 0.3279
## thsL1 1 0.016291 4.2785 -84.389 0.1108 0.7416
## thsL2 1 0.123593 4.3858 -83.126 0.8409 0.3667
## thsL3 1 0.089597 4.3518 -83.523 0.6096 0.4413# The same p-values in this case.
summary(m0)##
## Call:
## lm(formula = prod ~ . - id, data = na.omit(dd))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63186 -0.19456 -0.00795 0.20424 0.55444
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.73292 2.32740 2.034 0.0512 .
## fwcL1 9.28432 12.31087 0.754 0.4568
## fwcL2 -10.56898 10.20134 -1.036 0.3087
## fwcL3 15.49251 11.62199 1.333 0.1929
## IL1 -0.32758 0.27160 -1.206 0.2375
## IL2 -0.02300 0.43782 -0.053 0.9585
## IL3 -0.18492 0.36724 -0.504 0.6184
## macropL1 1.43271 4.91277 0.292 0.7726
## macropL2 -1.10623 0.96066 -1.152 0.2589
## macropL3 -0.13760 1.11080 -0.124 0.9023
## pawL1 -0.11254 0.15289 -0.736 0.4676
## pawL2 -0.03985 0.18836 -0.212 0.8339
## pawL3 -0.13358 0.12915 -1.034 0.3095
## SL1 3.42066 3.90597 0.876 0.3884
## SL2 -4.96332 5.60581 -0.885 0.3832
## SL3 -3.32867 5.88122 -0.566 0.5758
## thrL1 -14.50732 17.39383 -0.834 0.4111
## thrL2 -1.38649 22.86065 -0.061 0.9521
## thrL3 -28.96710 29.10900 -0.995 0.3279
## thsL1 2.17015 6.51839 0.333 0.7416
## thsL2 10.14001 11.05763 0.917 0.3667
## thsL3 12.08657 15.48021 0.781 0.4413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3834 on 29 degrees of freedom
## Multiple R-squared: 0.285, Adjusted R-squared: -0.2328
## F-statistic: 0.5504 on 21 and 29 DF, p-value: 0.9199# Stepwise model selection.
m1 <- step(m0,
direction = "both",
trace = FALSE,
k = c(AIC = 2, BIC = log(nrow(dd)))[1])
# Final model by AIC minimization.
summary(m1)##
## Call:
## lm(formula = prod ~ fwcL2 + thsL2, data = na.omit(dd))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73378 -0.25560 -0.00164 0.22372 0.73033
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7955 0.4646 10.322 8.87e-14 ***
## fwcL2 -5.3873 2.7509 -1.958 0.0560 .
## thsL2 3.9288 2.0767 1.892 0.0646 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.339 on 48 degrees of freedom
## Multiple R-squared: 0.07487, Adjusted R-squared: 0.03633
## F-statistic: 1.942 on 2 and 48 DF, p-value: 0.15455 Regression for the soybean yield on all physical and hydrical variables
dd <- db %>%
select(id, cam, names(leg))
# str(dd)
dd <- dd %>%
gather(key = "var", value = "val", rp:sand) %>%
mutate(varcam = paste(var, cam, sep = ".")) %>%
select(id, varcam, val) %>%
spread(key = "varcam", value = "val")
# names(dd)
dd <- inner_join(dd,
db %>% select(id, prod, MO2:HF) %>% filter(!is.na(prod)),
by = "id") %>%
na.omit()
dd## # A tibble: 51 x 51
## id clay.L1 clay.L2 clay.L3 dp.L1 dp.L2 dp.L3 ds.L1 ds.L2 ds.L3 fwc.L1 fwc.L2 fwc.L3
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 FAZA… 661. 722. 749. 2.78 2.67 2.57 1.36 1.35 1.22 0.408 0.391 0.400
## 2 FAZA… 653. 725. 775. 2.59 2.73 2.45 1.21 1.29 1.21 0.417 0.396 0.411
## 3 FAZA… 644. 724. 773. 2.63 2.65 2.59 1.04 1.19 1.14 0.449 0.431 0.456
## 4 FAZA… 713. 759. 788. 2.68 2.66 2.73 1.17 1.36 1.30 0.430 0.391 0.401
## 5 FAZA… 676. 752. 719. 2.81 2.59 2.61 1.28 1.33 1.28 0.412 0.397 0.430
## 6 FAZA… 669. 734. 755. 2.73 2.64 2.69 1.22 1.36 1.20 0.413 0.404 0.403
## 7 FAZC… 489. 533. 563. 2.59 2.67 2.62 1.36 1.44 1.44 0.334 0.333 0.364
## 8 FAZF… 414. 448. 469. 2.54 2.64 2.81 1.38 1.62 1.71 0.343 0.298 0.317
## 9 FAZG… 677. 682. 690. 2.65 2.74 2.78 1.01 1.44 1.36 0.434 0.363 0.380
## 10 FAZG… 682. 689. 701. 2.63 2.60 2.79 1.18 1.43 1.25 0.437 0.398 0.399
## # ... with 41 more rows, and 38 more variables: I.L1 <dbl>, I.L2 <dbl>, I.L3 <dbl>,
## # macrop.L1 <dbl>, macrop.L2 <dbl>, macrop.L3 <dbl>, paw.L1 <dbl>, paw.L2 <dbl>,
## # paw.L3 <dbl>, poros.L1 <dbl>, poros.L2 <dbl>, poros.L3 <dbl>, rp.L1 <dbl>,
## # rp.L2 <dbl>, rp.L3 <dbl>, sand.L1 <dbl>, sand.L2 <dbl>, sand.L3 <dbl>, S.L1 <dbl>,
## # S.L2 <dbl>, S.L3 <dbl>, tha.L1 <dbl>, tha.L2 <dbl>, tha.L3 <dbl>, thn.L1 <dbl>,
## # thn.L2 <dbl>, thn.L3 <dbl>, thr.L1 <dbl>, thr.L2 <dbl>, thr.L3 <dbl>, ths.L1 <dbl>,
## # ths.L2 <dbl>, ths.L3 <dbl>, prod <dbl>, MO2 <dbl>, FLF <dbl>, OLF <dbl>, HF <dbl># Model with all variables.
m0 <- lm(prod ~ ., data = dd[, -1])
summary(m0)##
## Call:
## lm(formula = prod ~ ., data = dd[, -1])
##
## Residuals:
## 1 2 3 4 5 6 7 8
## 0.0180171 -0.0945767 0.0020493 0.1615630 0.0142825 -0.0395378 -0.0748722 0.0050342
## 9 10 11 12 13 14 15 16
## -0.0354766 0.0312685 -0.2393481 0.0393703 -0.1211268 0.0542560 0.0262594 -0.0206635
## 17 18 19 20 21 22 23 24
## -0.0149006 0.0430683 -0.0108894 -0.0013548 -0.0496344 0.0131474 -0.0441891 0.1340352
## 25 26 27 28 29 30 31 32
## -0.0511132 0.0581303 -0.0064238 0.0186329 0.0259232 -0.1771054 -0.0250224 0.1077580
## 33 34 35 36 37 38 39 40
## -0.0229518 -0.0563769 0.0642484 0.1066889 0.0293836 -0.0459304 0.0303431 0.1105967
## 41 42 43 44 45 46 47 48
## 0.0152758 0.0233562 -0.0425460 0.0515014 0.0023443 -0.0245073 -0.0001946 -0.0189823
## 49 50 51
## 0.0694665 -0.0114304 -0.0268457
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.482e+01 3.717e+01 2.013 0.1003
## clay.L1 1.443e-04 4.349e-03 0.033 0.9748
## clay.L2 6.110e-03 3.757e-03 1.626 0.1649
## clay.L3 -6.903e-03 4.214e-03 -1.638 0.1624
## dp.L1 1.340e+01 9.078e+00 1.476 0.2001
## dp.L2 1.238e+00 1.099e+00 1.127 0.3109
## dp.L3 -1.860e+00 1.669e+00 -1.114 0.3158
## ds.L1 -3.080e+01 1.846e+01 -1.669 0.1560
## ds.L2 -1.437e+01 6.965e+00 -2.063 0.0941 .
## ds.L3 5.350e+00 7.250e+00 0.738 0.4937
## fwc.L1 -8.429e+01 6.129e+01 -1.375 0.2275
## fwc.L2 3.110e+01 3.334e+01 0.933 0.3938
## fwc.L3 -3.971e+01 3.952e+01 -1.005 0.3611
## I.L1 -6.581e+00 1.804e+00 -3.648 0.0148 *
## I.L2 -3.755e-01 1.359e+00 -0.276 0.7934
## I.L3 -6.244e-01 4.215e+00 -0.148 0.8880
## macrop.L1 -7.678e+01 4.967e+01 -1.546 0.1828
## macrop.L2 -1.147e+00 1.325e+00 -0.866 0.4261
## macrop.L3 -4.234e+00 2.297e+00 -1.843 0.1246
## paw.L1 2.871e-01 3.533e-01 0.813 0.4534
## paw.L2 4.063e-01 4.013e-01 1.012 0.3579
## paw.L3 -2.454e-01 3.143e-01 -0.781 0.4703
## poros.L1 NA NA NA NA
## poros.L2 NA NA NA NA
## poros.L3 NA NA NA NA
## rp.L1 5.246e-01 2.939e-01 1.785 0.1344
## rp.L2 3.114e-01 1.826e-01 1.705 0.1489
## rp.L3 5.091e-01 2.845e-01 1.790 0.1335
## sand.L1 8.419e-03 7.942e-03 1.060 0.3376
## sand.L2 1.340e-03 7.949e-03 0.169 0.8727
## sand.L3 -7.378e-03 8.636e-03 -0.854 0.4320
## S.L1 -2.616e+01 1.158e+01 -2.260 0.0734 .
## S.L2 7.479e+01 2.761e+01 2.709 0.0423 *
## S.L3 -1.407e+01 5.076e+01 -0.277 0.7928
## tha.L1 -5.555e+00 1.820e+00 -3.053 0.0283 *
## tha.L2 -1.404e+00 1.353e+00 -1.038 0.3469
## tha.L3 -8.996e-01 4.417e+00 -0.204 0.8466
## thn.L1 2.723e-01 5.571e-01 0.489 0.6457
## thn.L2 -3.006e+00 1.120e+00 -2.683 0.0437 *
## thn.L3 -4.510e-01 1.607e+00 -0.281 0.7902
## thr.L1 7.470e+00 3.469e+01 0.215 0.8380
## thr.L2 4.184e+01 5.354e+01 0.781 0.4699
## thr.L3 -4.687e+01 5.244e+01 -0.894 0.4124
## ths.L1 -6.362e-01 1.290e+01 -0.049 0.9626
## ths.L2 -7.624e+01 3.220e+01 -2.367 0.0642 .
## ths.L3 6.143e+01 3.907e+01 1.572 0.1767
## MO2 -1.098e-02 4.441e-02 -0.247 0.8145
## FLF 1.398e-03 1.754e-02 0.080 0.9395
## OLF 9.339e-02 3.626e-02 2.576 0.0497 *
## HF NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2215 on 5 degrees of freedom
## Multiple R-squared: 0.9589, Adjusted R-squared: 0.5886
## F-statistic: 2.59 on 45 and 5 DF, p-value: 0.1443# Selection.
m1 <- step(m0,
direction = "both",
trace = FALSE,
k = c(AIC = 2, BIC = log(nrow(dd)))[2])
summary(m1)##
## Call:
## lm(formula = prod ~ clay.L2 + clay.L3 + dp.L1 + dp.L2 + dp.L3 +
## ds.L1 + ds.L2 + ds.L3 + fwc.L1 + fwc.L2 + fwc.L3 + I.L1 +
## I.L3 + macrop.L1 + macrop.L2 + macrop.L3 + paw.L1 + paw.L2 +
## paw.L3 + rp.L1 + rp.L2 + rp.L3 + sand.L1 + sand.L3 + S.L1 +
## S.L2 + S.L3 + tha.L1 + tha.L2 + tha.L3 + thn.L1 + thn.L2 +
## thr.L2 + thr.L3 + ths.L2 + ths.L3 + MO2 + OLF, data = dd[,
## -1])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.218658 -0.034847 -0.005472 0.035689 0.162656
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.159839 17.381405 4.036 0.001650 **
## clay.L2 0.006012 0.001938 3.103 0.009143 **
## clay.L3 -0.006653 0.002123 -3.134 0.008630 **
## dp.L1 12.616134 4.606870 2.739 0.017980 *
## dp.L2 1.239317 0.629916 1.967 0.072683 .
## dp.L3 -1.530817 0.608616 -2.515 0.027145 *
## ds.L1 -30.189640 9.862083 -3.061 0.009877 **
## ds.L2 -14.567815 2.995455 -4.863 0.000389 ***
## ds.L3 6.815874 3.244148 2.101 0.057445 .
## fwc.L1 -74.650193 26.527110 -2.814 0.015630 *
## fwc.L2 36.606425 14.334169 2.554 0.025289 *
## fwc.L3 -49.546296 13.545621 -3.658 0.003280 **
## I.L1 -6.304119 0.734342 -8.585 1.81e-06 ***
## I.L3 -1.814630 1.534268 -1.183 0.259819
## macrop.L1 -71.996387 25.426623 -2.832 0.015133 *
## macrop.L2 -1.433614 0.558451 -2.567 0.024676 *
## macrop.L3 -4.777030 0.826577 -5.779 8.75e-05 ***
## paw.L1 0.211563 0.024465 8.647 1.68e-06 ***
## paw.L2 0.275495 0.154461 1.784 0.099780 .
## paw.L3 -0.192812 0.140354 -1.374 0.194635
## rp.L1 0.557220 0.130369 4.274 0.001080 **
## rp.L2 0.333454 0.095552 3.490 0.004465 **
## rp.L3 0.498221 0.139102 3.582 0.003771 **
## sand.L1 0.009320 0.002633 3.540 0.004074 **
## sand.L3 -0.006937 0.003169 -2.189 0.049110 *
## S.L1 -25.528630 4.210821 -6.063 5.65e-05 ***
## S.L2 73.729809 13.367072 5.516 0.000133 ***
## S.L3 -29.914968 7.664192 -3.903 0.002098 **
## tha.L1 -5.215967 0.642006 -8.124 3.21e-06 ***
## tha.L2 -1.015303 0.285043 -3.562 0.003910 **
## tha.L3 -2.136361 1.560249 -1.369 0.196005
## thn.L1 0.298845 0.170507 1.753 0.105142
## thn.L2 -2.824535 0.466093 -6.060 5.67e-05 ***
## thr.L2 26.510417 21.185762 1.251 0.234659
## thr.L3 -40.846364 25.450243 -1.605 0.134484
## ths.L2 -68.173368 13.601603 -5.012 0.000303 ***
## ths.L3 68.107417 18.589076 3.664 0.003243 **
## MO2 -0.021534 0.010602 -2.031 0.065008 .
## OLF 0.101108 0.017068 5.924 6.99e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1468 on 12 degrees of freedom
## Multiple R-squared: 0.9566, Adjusted R-squared: 0.8192
## F-statistic: 6.963 on 38 and 12 DF, p-value: 0.0004561# Diagnostics.
par(mfrow = c(2, 2))
plot(m0)
layout(1)6 Additional relations between soil variables
#-----------------------------------------------------------------------
dd <- list()
dd[[1]] <- db %>%
filter(cam == "L1") %>%
select(faz, tal, MO2, FLF, OLF, HF)
dd[[2]] <- db %>%
select(faz, tal, cam, rp, paw, S, clay, sand)
# Join the tables.
dd <- left_join(dd[[2]], dd[[1]])
# Replace an outlier by a missing value.
i <- which(dd$paw > 140)
dd$paw[i] <- NA
i <- which(dd$S > 0.3)
dd$S[i] <- NA
# Does stack two times to create all pairs.
de <- dd %>%
gather(key = "xvar", value = "xval", rp:sand) %>%
gather(key = "yvar", value = "yval", MO2:HF) %>%
filter(is.finite(yval), is.finite(xval))
ggplot(data = de,
mapping = aes(x = xval,
y = yval,
colour = cam)) +
geom_point() +
geom_smooth(method = lm) +
facet_grid(facets = yvar ~ xvar,
scales = "free")
# unique(de$xvar)
# levels(factor(de$yvar))
# levels(factor(de$xvar))
de <- de %>%
mutate(
xvar = factor(xvar,
levels = c("clay", "sand", "paw", "S", "rp"),
labels = c("Clay", "Sand", "PAW", "S-index", "SRP")),
yvar = factor(yvar,
levels = c("HF", "FLF", "OLF", "MO2"),
labels = c("HF of SOC", "FLF of SOC", "OLF of SOC", "SOM"))
)
dem <- de %>%
group_by(yvar, xvar) %>%
summarise(ym = mean(yval))
ggplot(data = filter(de, cam == "L1"),
mapping = aes(x = xval,
y = yval)) +
geom_point() +
geom_smooth(method = lm, color = "black") +
# geom_smooth(method = lm, formula = y ~ 1, se = FALSE,
# linetype = 2, col = "black", size = 1) +
geom_hline(data = dem,
mapping = aes(yintercept = ym),
linetype = 2) +
facet_grid(facets = yvar ~ xvar,
scales = "free",
labeller = labeller(variable = label_legend)) +
xlab("Values of the regressor variable") +
ylab("Values of the response variable")
#-----------------------------------------------------------------------
# Function to ease analyse the relations.
do_fit <- function(y, x, xlab, ylab, lay) {
# To store the results.
res <- list()
# Get the data.
dd$x <- dd[[x]]
dd$y <- dd[[y]]
# Exploratory plot.
gg0 <- ggplot(data = dd,
mapping = aes(x = x, y = y)) +
geom_point() +
geom_smooth(span = 1) +
geom_hline(yintercept = mp, lty = 2) +
facet_wrap(facets = ~cam)
res$gg0 <- gg0
# Reshape data.
de <- dd %>%
select(tal, cam, x, y) %>%
spread(key = "cam", value = "x") %>%
na.omit()
# Response mean.
mp <- mean(de$y)
# # Second order polynomial regression.
# m0 <- lm(y ~
# poly(L1, degree = 2) +
# poly(L2, degree = 2) +
# poly(L3, degree = 2),
# data = de)
# res$m0_drop <- drop1(m0, test = "F")
# res$m0_summary <- summary(m0)
# # First order regression.
# m1 <- update(m0, . ~ L1 + L2 + L3)
# res$m1_summary <- summary(m1)
# # Only one layer.
# m2 <- update(m0, as.formula(sprintf(". ~ %s", lay)))
m2 <- lm(as.formula(sprintf("y ~ %s", lay)),
data = de)
res$m2_summary <- summary(m2)
# Final plot.
gg <- ggplot(data = de,
mapping = aes(x = eval(parse(text = lay)),
y = y)) +
geom_point() +
geom_smooth(method = lm, formula = y ~ x, col = "black") +
geom_hline(yintercept = mp, lty = 2) +
xlab(label = eval(parse(text = xlab))) +
ylab(label = eval(parse(text = ylab)))
res$gg <- gg
# Return the results.
return(res)
}6.1 HF x clay
fit <- do_fit(x = "clay",
y = "HF",
xlab = leg$clay,
ylab = leg$HF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.9953 -3.6769 -0.5548 3.6460 14.8072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.217949 3.630914 1.988 0.0522 .
## L1 0.012539 0.006031 2.079 0.0427 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.739 on 51 degrees of freedom
## Multiple R-squared: 0.07813, Adjusted R-squared: 0.06005
## F-statistic: 4.322 on 1 and 51 DF, p-value: 0.04267tail(fit, n = 1)## $gg
6.2 HF x sand
fit <- do_fit(x = "sand",
y = "HF",
xlab = leg$sand,
ylab = leg$HF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.1913 -3.3716 -0.4141 3.6855 14.1077
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.930111 1.699697 10.549 2.02e-14 ***
## L1 -0.011073 0.004994 -2.217 0.0311 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.708 on 51 degrees of freedom
## Multiple R-squared: 0.08792, Adjusted R-squared: 0.07004
## F-statistic: 4.916 on 1 and 51 DF, p-value: 0.03109tail(fit, n = 1)## $gg
6.3 HF x paw
fit <- do_fit(x = "paw",
y = "HF",
xlab = leg$paw,
ylab = leg$HF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7308 -5.3082 0.5582 3.8358 15.4102
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.59968 2.14137 6.351 6.22e-08 ***
## L1 0.05752 0.10151 0.567 0.574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.936 on 50 degrees of freedom
## Multiple R-squared: 0.00638, Adjusted R-squared: -0.01349
## F-statistic: 0.3211 on 1 and 50 DF, p-value: 0.5735tail(fit, n = 1)## $gg
6.4 HF x S
fit <- do_fit(x = "S",
y = "HF",
xlab = leg$S,
ylab = leg$HF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.3661 -5.0753 0.6003 3.8263 15.7411
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.978 2.479 6.042 1.88e-07 ***
## L1 -4.752 22.977 -0.207 0.837
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.996 on 50 degrees of freedom
## Multiple R-squared: 0.0008546, Adjusted R-squared: -0.01913
## F-statistic: 0.04277 on 1 and 50 DF, p-value: 0.837tail(fit, n = 1)## $gg
6.5 HF x rp
fit <- do_fit(x = "rp",
y = "HF",
xlab = leg$rp,
ylab = leg$HF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.4121 -5.0275 0.8329 3.8905 15.6272
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.3256 2.0329 7.047 4.59e-09 ***
## L1 0.1636 1.1665 0.140 0.889
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.976 on 51 degrees of freedom
## Multiple R-squared: 0.0003854, Adjusted R-squared: -0.01921
## F-statistic: 0.01966 on 1 and 51 DF, p-value: 0.889tail(fit, n = 1)## $gg
6.6 FLF x clay
fit <- do_fit(x = "clay",
y = "FLF",
xlab = leg$clay,
ylab = leg$FLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6045 -2.4339 0.2239 1.8550 8.8188
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.908134 2.538758 4.297 7.8e-05 ***
## L1 0.007080 0.004217 1.679 0.0993 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.013 on 51 degrees of freedom
## Multiple R-squared: 0.05237, Adjusted R-squared: 0.03379
## F-statistic: 2.819 on 1 and 51 DF, p-value: 0.0993tail(fit, n = 1)## $gg
6.7 FLF x sand
fit <- do_fit(x = "sand",
y = "FLF",
xlab = leg$sand,
ylab = leg$FLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.4869 -2.6018 0.3722 2.0311 8.5677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.090630 1.185165 14.420 <2e-16 ***
## L1 -0.006696 0.003482 -1.923 0.0601 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.98 on 51 degrees of freedom
## Multiple R-squared: 0.0676, Adjusted R-squared: 0.04932
## F-statistic: 3.697 on 1 and 51 DF, p-value: 0.06009tail(fit, n = 1)## $gg
6.8 FLF x paw
fit <- do_fit(x = "paw",
y = "FLF",
xlab = leg$paw,
ylab = leg$FLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.6793 -2.4994 0.2182 2.1308 9.4356
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.42799 1.49629 9.643 5.39e-13 ***
## L1 0.03118 0.07093 0.440 0.662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.148 on 50 degrees of freedom
## Multiple R-squared: 0.00385, Adjusted R-squared: -0.01607
## F-statistic: 0.1932 on 1 and 50 DF, p-value: 0.6621tail(fit, n = 1)## $gg
6.9 FLF x S
fit <- do_fit(x = "S",
y = "FLF",
xlab = leg$S,
ylab = leg$FLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.311 -2.862 -0.083 2.678 7.837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.30 1.67 7.369 1.59e-09 ***
## L1 26.98 15.48 1.743 0.0875 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.039 on 50 degrees of freedom
## Multiple R-squared: 0.05727, Adjusted R-squared: 0.03842
## F-statistic: 3.037 on 1 and 50 DF, p-value: 0.08751tail(fit, n = 1)## $gg
6.10 FLF x rp
fit <- do_fit(x = "rp",
y = "FLF",
xlab = leg$rp,
ylab = leg$FLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.8204 -2.2952 0.1938 1.7182 9.5008
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.6825 1.4010 10.480 2.54e-14 ***
## L1 0.2422 0.8039 0.301 0.764
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.118 on 51 degrees of freedom
## Multiple R-squared: 0.001777, Adjusted R-squared: -0.0178
## F-statistic: 0.0908 on 1 and 51 DF, p-value: 0.7644tail(fit, n = 1)## $gg
6.11 OLF x clay
fit <- do_fit(x = "clay",
y = "OLF",
xlab = leg$clay,
ylab = leg$OLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.8274 -1.8128 -0.7500 0.9472 8.2848
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.532949 1.878914 1.348 0.1836
## L1 0.005752 0.003121 1.843 0.0712 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.97 on 51 degrees of freedom
## Multiple R-squared: 0.06244, Adjusted R-squared: 0.04405
## F-statistic: 3.396 on 1 and 51 DF, p-value: 0.07116tail(fit, n = 1)## $gg
6.12 OLF x sand
fit <- do_fit(x = "sand",
y = "OLF",
xlab = leg$sand,
ylab = leg$OLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7341 -1.8034 -0.6824 0.7680 8.1288
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.576438 0.874685 8.662 1.36e-11 ***
## L1 -0.005509 0.002570 -2.143 0.0369 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.938 on 51 degrees of freedom
## Multiple R-squared: 0.08264, Adjusted R-squared: 0.06465
## F-statistic: 4.594 on 1 and 51 DF, p-value: 0.03687tail(fit, n = 1)## $gg
6.13 OLF x paw
fit <- do_fit(x = "paw",
y = "OLF",
xlab = leg$paw,
ylab = leg$OLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1992 -1.6686 -0.6862 0.7878 9.1723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.14225 1.03357 5.943 2.69e-07 ***
## L1 -0.01984 0.04900 -0.405 0.687
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.865 on 50 degrees of freedom
## Multiple R-squared: 0.003267, Adjusted R-squared: -0.01667
## F-statistic: 0.1639 on 1 and 50 DF, p-value: 0.6873tail(fit, n = 1)## $gg
6.14 OLF x S
fit <- do_fit(x = "S",
y = "OLF",
xlab = leg$S,
ylab = leg$OLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1807 -1.8089 -0.7099 0.5465 8.9778
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.460 1.276 4.278 8.51e-05 ***
## L1 4.196 11.831 0.355 0.724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.087 on 50 degrees of freedom
## Multiple R-squared: 0.002509, Adjusted R-squared: -0.01744
## F-statistic: 0.1258 on 1 and 50 DF, p-value: 0.7244tail(fit, n = 1)## $gg
6.15 OLF x rp
fit <- do_fit(x = "rp",
y = "OLF",
xlab = leg$rp,
ylab = leg$OLF,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4654 -1.8912 -0.3885 0.5056 9.4250
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.3409 1.0152 4.276 8.35e-05 ***
## L1 0.9861 0.5825 1.693 0.0966 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.984 on 51 degrees of freedom
## Multiple R-squared: 0.0532, Adjusted R-squared: 0.03463
## F-statistic: 2.865 on 1 and 51 DF, p-value: 0.0966tail(fit, n = 1)## $gg
6.16 MO2 x clay
fit <- do_fit(x = "clay",
y = "MO2",
xlab = leg$clay,
ylab = leg$MO2,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.9769 -3.4351 -0.6393 3.2453 8.2925
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.659032 2.674910 7.723 3.94e-10 ***
## L1 0.025370 0.004443 5.710 5.82e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.228 on 51 degrees of freedom
## Multiple R-squared: 0.39, Adjusted R-squared: 0.378
## F-statistic: 32.6 on 1 and 51 DF, p-value: 5.82e-07tail(fit, n = 1)## $gg
6.17 MO2 x sand
fit <- do_fit(x = "sand",
y = "MO2",
xlab = leg$sand,
ylab = leg$MO2,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.3472 -2.8860 -0.5997 2.4401 7.1677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.597179 1.169268 36.431 < 2e-16 ***
## L1 -0.023278 0.003436 -6.776 1.23e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.927 on 51 degrees of freedom
## Multiple R-squared: 0.4737, Adjusted R-squared: 0.4634
## F-statistic: 45.91 on 1 and 51 DF, p-value: 1.231e-08tail(fit, n = 1)## $gg
6.18 MO2 x paw
fit <- do_fit(x = "paw",
y = "MO2",
xlab = leg$paw,
ylab = leg$MO2,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.247 -3.771 1.043 3.013 10.425
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.16992 1.95553 17.473 <2e-16 ***
## L1 0.06886 0.09270 0.743 0.461
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.421 on 50 degrees of freedom
## Multiple R-squared: 0.01092, Adjusted R-squared: -0.008866
## F-statistic: 0.5518 on 1 and 50 DF, p-value: 0.4611tail(fit, n = 1)## $gg
6.19 MO2 x S
fit <- do_fit(x = "S",
y = "MO2",
xlab = leg$S,
ylab = leg$MO2,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.136 -2.307 1.319 2.951 10.996
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.743 2.182 15.005 <2e-16 ***
## L1 26.419 20.227 1.306 0.197
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.278 on 50 degrees of freedom
## Multiple R-squared: 0.03299, Adjusted R-squared: 0.01365
## F-statistic: 1.706 on 1 and 50 DF, p-value: 0.1975tail(fit, n = 1)## $gg
6.20 MO2 x rp
fit <- do_fit(x = "rp",
y = "MO2",
xlab = leg$rp,
ylab = leg$MO2,
lay = "L1")
head(fit[-1], n = -1)## $m2_summary
##
## Call:
## lm(formula = as.formula(sprintf("y ~ %s", lay)), data = de)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.2115 -3.6442 0.5151 3.5726 10.6883
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.349 1.810 18.43 <2e-16 ***
## L1 1.392 1.038 1.34 0.186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.32 on 51 degrees of freedom
## Multiple R-squared: 0.03403, Adjusted R-squared: 0.01508
## F-statistic: 1.796 on 1 and 51 DF, p-value: 0.1861tail(fit, n = 1)## $gg