Fitting hydraulic vulnerability curves

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
)

## Warning: `fit_many()` was deprecated in photosynthesis 2.1.3.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# 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)

References

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.

Joseph R. Stinziano and Christopher D. Muir

2023-08-15

References