Current approach fits a sigmoidal model and calculates hydraulic parameters from the curve fit (Pammenter & Van der Willigen 1998; Ogle et al. 2009).
library(dplyr)
library(photosynthesis)
# Read in data
dat = system.file("extdata", "hydraulic_vulnerability.csv", package = "photosynthesis") |>
read.csv() |>
mutate(ID = paste(Plot, Tree, sep = "_")) |>
rename(psi = P)
# Fit hydraulic vulnerability curve
fit = fit_hydra_vuln_curve(
filter(dat, Tree == 5, Plot == "Irrigation"),
start_weibull = list(a = 2, b = 1),
title = "Irrigation 5"
)
# Return Sigmoidal model summary
summary(fit[[1]]) ##
## Call:
## lm(formula = H_log ~ psi, data = data[data$H_log < Inf, ])
##
## Residuals:
## 2 3 4 5 6
## 0.40236 -0.63441 0.01791 0.09292 0.12121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1700 0.5344 9.675 0.00234 **
## psi -1.0884 0.1212 -8.982 0.00291 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4427 on 3 degrees of freedom
## Multiple R-squared: 0.9642, Adjusted R-squared: 0.9522
## F-statistic: 80.68 on 1 and 3 DF, p-value: 0.002912# Return Weibull model summary
summary(fit[[4]]) #expecting a = 4.99, b = 3.22##
## Formula: K.Kmax ~ exp(-((psi/a)^b))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 5.3160 0.0902 58.93 4.96e-07 ***
## b 2.7778 0.2393 11.61 0.000315 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02867 on 4 degrees of freedom
##
## Number of iterations to convergence: 8
## Achieved convergence tolerance: 1.49e-08# Return model parameters with 95% confidence intervals
fit[[2]] ## Value Parameter Curve
## b...1 4.749922 b Sigmoidal
## a...2 -1.088445 a Sigmoidal
## b...3 2.777799 b Weibull
## a...4 5.315979 a Weibull# Return hydraulic parameters
fit[[3]] ## P25 P50 P88 P95 S50 Pe Pmax DSI
## 1 3.740581 4.749922 6.580451 7.455102 27.21113 2.912439 6.587406 3.674967
## 2 3.394637 4.658873 6.967591 7.890836 20.66405 2.239211 7.078534 4.839322
## Curve
## 1 Sigmoidal
## 2 Weibull# Return graph
# fit[[5]]
# Fit many curves
fits = fit_many(
data = dat,
group = "ID",
start_weibull = list(a = 4, b = 2),
funct = fit_hydra_vuln_curve,
progress = FALSE
)
# Return model summary
summary(fits[[1]][[1]])##
## Call:
## lm(formula = H_log ~ psi, data = data[data$H_log < Inf, ])
##
## Residuals:
## 44 45 46 47 48
## -0.15427 0.06136 0.23623 0.20568 -0.34900
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.72729 0.34666 16.52 0.000483 ***
## psi -1.41591 0.07861 -18.01 0.000373 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2872 on 3 degrees of freedom
## Multiple R-squared: 0.9908, Adjusted R-squared: 0.9878
## F-statistic: 324.4 on 1 and 3 DF, p-value: 0.0003733# Return sigmoidal model output
fits[[1]][[2]] ## Value Parameter Curve
## b...1 4.044964 b Sigmoidal
## a...2 -1.415905 a Sigmoidal
## b...3 3.905565 b Weibull
## a...4 4.580836 a Weibull# Return hydraulic parameters
fits[[1]][[3]] ## P25 P50 P88 P95 S50 Pe Pmax DSI
## 1 3.269056 4.044964 5.452141 6.124509 35.39764 2.632440 5.457488 2.825047
## 2 3.329677 4.170508 5.552841 6.066678 32.45566 2.629945 5.711071 3.081127
## Curve
## 1 Sigmoidal
## 2 Weibull# Return graph
# fits[[1]][[5]]
# Compile parameter outputs
pars = compile_data(data = fits, output_type = "dataframe", list_element = 3)
# Compile graphs
graphs = compile_data(data = fits, output_type = "list", list_element = 5)Ogle K, Barber JJ, Willson C, Thompson B. 2009. Hierarchical statistical modeling of xylem vulnerability to cavitation. New Phytologist 182:541-554.
Pammenter NW, Van der Willigen CV. 1998. A mathematical and statistical analysis of the curves illustrating vulnerability of xylem to cavitation. Tree Physiology 18:589-593.