Skip to contents

Fit epidemic models using nonlinear regression

Source: R/fit_nlin.R
fit_nlin.Rd

Fit exponential, monomolecular, logistic, and Gompertz models to disease progress data using nonlinear regression.

Usage

fit_nlin(
  time,
  y,
  starting_par = list(y0 = 0.01, r = 0.03),
  maxiter = 50,
  weights = NULL,
  weight_method = c("none", "binomial", "mean", "cv", "power"),
  weight_eps = 0.01,
  weight_power = 1
)

Arguments

time

Numeric vector of assessment times.

y

Numeric vector of disease intensity values.

starting_par

Named list with starting values for `y0` and `r`. When omitted or partially specified, `epifitter` supplies data-driven fallback values.

maxiter

Maximum number of iterations. Must be a positive number.

weights

Optional numeric vector of positive weights, or a function that receives a data frame with columns `time`, `y`, `predicted`, and `model` and returns positive weights. When supplied, these weights are used directly and `weight_method` must be `"none"`.

weight_method

Weighting strategy for nonlinear least squares. Use `"none"` for ordinary nonlinear least squares, `"binomial"` for \(1 / (\hat{p}(1 - \hat{p}) + \epsilon)\), `"mean"` for \(1 / (\hat{p} + \epsilon)\), `"cv"` for \(1 / (\hat{p}^2 + \epsilon)\), or `"power"` for \(1 / (|\hat{p}|^{2\theta} + \epsilon)\).

weight_eps

Small positive constant added to fitted-value variance approximations to keep weights finite near 0 and 1.

weight_power

Non-negative power \(\theta\) used only when `weight_method = "power"`.

Value

A list with fit statistics, parameter estimates, and prediction data.

Details

Weighted fits use weighted nonlinear least squares, not a binomial or beta likelihood. `weight_method` options other than `"none"` use a two-step working approximation: first fit the model without weights, derive weights from the fitted values, and then refit the model. Report the selected weighting rule as an assumption.

Examples

set.seed(1)
epi <- sim_logistic(N = 30, y0 = 0.01, dt = 5, r = 0.3, alpha = 0.2, n = 4)
fit_nlin(time = epi$time, y = epi$random_y, starting_par = list(y0 = 0.01, r = 0.03))
#> Warning: The following nonlinear model(s) did not converge and will be returned as NA: Monomolecular.
#> Results of fitting population models 
#> 
#> Stats:
#>                  CCC r_squared    RSE
#> Gompertz      0.9986    0.9981 0.0217
#> Logistic      0.9986    0.9974 0.0220
#> Exponential   0.9038    0.8493 0.1665
#> Monomolecular     NA        NA     NA
#> 
#>  Infection rate:
#>                 Estimate   Std.error      Lower     Upper
#> Gompertz      0.22500335 0.006307726 0.21203764 0.2379691
#> Logistic      0.31571569 0.009477990 0.29623341 0.3351980
#> Exponential   0.06632831 0.007801830 0.05029142 0.0823652
#> Monomolecular         NA          NA         NA        NA
#> 
#>  Initial inoculum:
#>                   Estimate    Std.error         Lower        Upper
#> Gompertz      1.279015e-08 2.122676e-08 -3.084208e-08 5.642237e-08
#> Logistic      8.926801e-03 1.286286e-03  6.282802e-03 1.157080e-02
#> Exponential   1.560581e-01 3.149510e-02  9.131899e-02 2.207972e-01
#> Monomolecular           NA           NA            NA           NA