Skip to contents

Profile-likelihood confidence intervals for fitted VG parameters

Source: R/confint_fit_swrc.R
confint.fit_swrc.Rd

Computes confidence intervals for the Van Genuchten parameters returned by fit_swrc() using the profile likelihood ratio (PLR) approach. For each free parameter, one dimension of the parameter space is scanned while the remaining free parameters are re-optimised at each step. The interval endpoints are the values at which the profile F-statistic exceeds the critical value:

Usage

# S3 method for class 'fit_swrc'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

object

An object of class fit_swrc, as returned by fit_swrc().

parm

Character vector of parameter names for which to compute intervals. Defaults to all free parameters (those not fixed in the original call to fit_swrc()).

level

Confidence level. Defaults to 0.95.

...

Not used; present for S3 compatibility.

Value

A tibble with columns:

  • Group keys (if the original fit was grouped)

  • param — parameter name

  • estimate — point estimate (MLE)

  • lower — lower confidence bound

  • upper — upper confidence bound

  • level — confidence level

Details

$$F(p) = \frac{RSS(p) - RSS_0}{RSS_0 / (n - k)} \leq F_{1,\, n-k}(\text{level})$$

where RSS₀ is the minimum (fitted) residual sum of squares, n is the number of observations, and k is the number of free parameters.

PLR intervals are generally more reliable than normal-approximation intervals (i.e., estimate ± z · SE) for the non-linear VG model, especially when sample sizes are small or parameter distributions are skewed.

See also

fit_swrc(), fit_swrc_hcc()

Examples

library(tibble)

observed <- tibble(
  h     = c(0, 10, 100, 500, 1000, 5000, 15000),
  theta = c(0.44, 0.40, 0.32, 0.25, 0.20, 0.12, 0.08)
)

fit <- fit_swrc(observed, theta_col = theta, h_col = h)
confint(fit)
#> # A tibble: 4 × 5
#>   param   estimate   lower upper level
#>   <chr>      <dbl>   <dbl> <dbl> <dbl>
#> 1 theta_r   0      0       0.119  0.95
#> 2 theta_s   0.425  0.374   0.484  0.95
#> 3 alpha     0.0182 0.00335 0.173  0.95
#> 4 n         1.27   1.16    1.82   0.95

# 90% intervals for alpha and n only
confint(fit, parm = c("alpha", "n"), level = 0.90)
#> # A tibble: 2 × 5
#>   param estimate   lower  upper level
#>   <chr>    <dbl>   <dbl>  <dbl> <dbl>
#> 1 alpha   0.0182 0.00527 0.0885   0.9
#> 2 n       1.27   1.18    1.58     0.9

# Grouped: returns CIs for each group
library(dplyr)
#> 
#> Attaching package: ‘dplyr’
#> The following objects are masked from ‘package:stats’:
#> 
#>     filter, lag
#> The following objects are masked from ‘package:base’:
#> 
#>     intersect, setdiff, setequal, union
tibble(
  horizon = rep(c("A", "B"), each = 7),
  h       = rep(c(0, 10, 100, 500, 1000, 5000, 15000), 2),
  theta   = c(0.44, 0.40, 0.32, 0.25, 0.20, 0.12, 0.08,
              0.38, 0.33, 0.27, 0.21, 0.16, 0.10, 0.07)
) |>
  group_by(horizon) |>
  fit_swrc(theta_col = theta, h_col = h) |>
  confint()
#> # A tibble: 8 × 6
#>   horizon param   estimate   lower upper level
#>   <chr>   <chr>      <dbl>   <dbl> <dbl> <dbl>
#> 1 A       theta_r   0      0       0.119  0.95
#> 2 A       theta_s   0.425  0.374   0.484  0.95
#> 3 A       alpha     0.0182 0.00335 0.173  0.95
#> 4 A       n         1.27   1.16    1.82   0.95
#> 5 B       theta_r   0      0       0.116  0.95
#> 6 B       theta_s   0.361  0.309   0.430  0.95
#> 7 B       alpha     0.0215 0.00256 0.459  0.95
#> 8 B       n         1.26   1.14    2.20   0.95