Fitting and Comparing Severity Candidates

Workflow and data interpretation

fit_severity_candidates() fits competing parametric ground-up severity families and returns a list with fits, a comparison data frame, and excluded_rows / excluded_details. The observations must be ground-up; net paid losses should not be substituted without reconstructing that basis. A finite policy censoring limit means values at the limit are right-censored; no limit or deductible column is needed if it does not apply.

set.seed(128)
n <- 750
latent <- rlnorm(n, log(14000), 1.0)
limits <- sample(c(20000, 50000, 100000, Inf), n, replace = TRUE)
dat <- data.frame(ground_up_loss = pmin(latent, limits),
                  censor_limit = limits)

Fit, compare, and extract a selected fit

candidates <- fit_severity_candidates(dat, loss = "ground_up_loss",
  limit = "censor_limit", distributions = c("lognormal", "gamma", "weibull"),
  hessian = FALSE)
print(candidates$comparison)
#>   distribution converged    logLik      AIC      BIC parameters   n n_censored
#> 1    lognormal      TRUE -7254.734 14513.47 14522.71          2 750         85
#> 2        gamma      TRUE -7306.483 14616.97 14626.21          2 750         85
#> 3      weibull      TRUE -7311.927 14627.85 14637.09          2 750         85
#>   error elapsed_seconds function_evaluations gradient_evaluations
#> 1  <NA>            0.03                   25                    5
#> 2  <NA>            0.22                  103                   29
#> 3  <NA>            0.06                   35                   11
#>   convergence_code
#> 1                0
#> 2                0
#> 3                0
candidates$comparison[order(candidates$comparison$AIC), ]
#>   distribution converged    logLik      AIC      BIC parameters   n n_censored
#> 1    lognormal      TRUE -7254.734 14513.47 14522.71          2 750         85
#> 2        gamma      TRUE -7306.483 14616.97 14626.21          2 750         85
#> 3      weibull      TRUE -7311.927 14627.85 14637.09          2 750         85
#>   error elapsed_seconds function_evaluations gradient_evaluations
#> 1  <NA>            0.03                   25                    5
#> 2  <NA>            0.22                  103                   29
#> 3  <NA>            0.06                   35                   11
#>   convergence_code
#> 1                0
#> 2                0
#> 3                0
fit_lognormal <- candidates$fits$lognormal
summary(fit_lognormal)
#> Actuarial severity model summary
#> 
#> Model: single
#> Distribution: lognormal
#> Observations: 750 (effective weighted count: 750)
#> Right-censored: 85
#> LogLik: -7254.734  AIC: 14513.47  BIC: 14522.71
#> 
#>           estimate std.error
#> scale 13642.397983        NA
#> sdlog     1.015471        NA
#> 
#> The internal scale is exp(meanlog). With scale_by, mean severity is proportional to scale_coef * scale_by + scale_intercept.

Lower AIC/BIC is better relative to other fits on the same observation model. A higher log likelihood is better. These are not absolute goodness-of-fit measures and can hide tail or layer errors. Not every distribution should be forced onto every loss dataset; inspect convergence and factor stability.

Visual comparison

Modeled portfolio probabilities reproduce row-level censoring limits and deductibles. Thus the histogram comparisons are made on the observed basis. The new analytic default computes expected bin counts without Monte Carlo simulation; method = "simulation" restores the stochastic plot. The default horizontal axis is logarithmic, and specifying xlim changes only the displayed range. Use scale_range = c(lower, upper) when fitting with scale_by.

plot_candidate_fits(candidates, type = "histogram",
  x_scale = "log10", xlim = c(500, 200000),
  x_text_angle = 45, x_text_size = 8)

# Deterministic summaries: the very same bins drawn in the default plot.
head(candidate_plot_data(candidates, output = "plot_ready",
  type = "histogram", xlim = c(500, 200000)))
#>   source         x      value     lower     upper distribution
#> 1 actual  538.8846 0.00000000  500.0000  580.7932    lognormal
#> 2 actual  625.9610 0.02049658  580.7932  674.6414    lognormal
#> 3 actual  727.1077 0.02049658  674.6414  783.6543    lognormal
#> 4 actual  844.5984 0.02049658  783.6543  910.2821    lognormal
#> 5 actual  981.0740 0.06148974  910.2821 1057.3713    lognormal
#> 6 actual 1139.6021 0.02049658 1057.3713 1228.2280    lognormal

# Old row-level simulations remain available when requested.
head(candidate_plot_data(candidates, nsim = 1, seed = 22, output = "raw"))
#>   row source simulation observed_loss censored deductible limit scale_value
#> 1   1 actual          0      20000.00     TRUE          0 2e+04          NA
#> 2   2 actual          0      22684.22    FALSE          0 1e+05          NA
#> 3   3 actual          0      20000.00     TRUE          0 2e+04          NA
#> 4   4 actual          0      13640.75    FALSE          0 1e+05          NA
#> 5   5 actual          0      20000.00     TRUE          0 2e+04          NA
#> 6   6 actual          0      25013.35    FALSE          0   Inf          NA
#>   weight distribution
#> 1      1    lognormal
#> 2      1    lognormal
#> 3      1    lognormal
#> 4      1    lognormal
#> 5      1    lognormal
#> 6      1    lognormal

Convert fits to rating factors and assess sensitivity

limit in severity_factors() is a payment limit after the deductible: E[min(max(v*X-d,0),L)]. This is distinct from the ground-up censoring limit used for fitting. In particular, with valuation_factor = 1, finite payment limits require integrated survival between d and d + L.

factors <- severity_factors(candidates, deductible = c(500, 1000),
  limit = c(25000, 50000, Inf), base_deductible = 500, base_limit = 25000)
head(factors)
#>             distribution  model scale_value deductible limit valuation_factor
#> lognormal.1    lognormal single          NA        500 25000                1
#> lognormal.2    lognormal single          NA       1000 25000                1
#> lognormal.3    lognormal single          NA        500 50000                1
#> lognormal.4    lognormal single          NA       1000 50000                1
#> lognormal.5    lognormal single          NA        500   Inf                1
#> lognormal.6    lognormal single          NA       1000   Inf                1
#>             expected_payment limited_expected_severity unlimited_mean
#> lognormal.1         14234.88                  14598.73       22845.93
#> lognormal.2         13868.96                  14598.73       22845.93
#> lognormal.3         18369.59                  18819.74       22845.93
#> lognormal.4         17919.71                  18819.74       22845.93
#> lognormal.5         22345.99                  22845.93       22845.93
#> lognormal.6         21847.17                  22845.93       22845.93
#>                    LER      ILF rating_factor base_deductible base_limit
#> lognormal.1 0.02188308 1.000000     1.0000000             500      25000
#> lognormal.2 0.04371713 1.000000     0.9742944             500      25000
#> lognormal.3 0.02188308 1.290463     1.2904632             500      25000
#> lognormal.4 0.04371713 1.292073     1.2588593             500      25000
#> lognormal.5 0.02188308 1.569805     1.5698055             500      25000
#> lognormal.6 0.04371713 1.575256     1.5347634             500      25000
#>             base_valuation_factor
#> lognormal.1                     1
#> lognormal.2                     1
#> lognormal.3                     1
#> lognormal.4                     1
#> lognormal.5                     1
#> lognormal.6                     1
compare_factors(candidates, metric = "rating_factor",
  deductible = c(500, 1000), limit = c(25000, 50000),
  base_deductible = 500, base_limit = 25000)
#>             scale_value deductible limit valuation_factor base_deductible
#> lognormal.1          NA        500 25000                1             500
#>             base_limit base_valuation_factor value.lognormal value.gamma
#> lognormal.1      25000                     1               1           1
#>             value.weibull
#> lognormal.1             1

By default, the grid is Cartesian. Use grid = "paired" to match vectors positionally, allowing scalar recycling. This prevents unexpectedly generating 10,000 scenarios from two vectors of length 100.

Excluded claims

Rows with below-deductible or above-censor-limit losses are excluded by default with a warning. The original positions and reasons are retained. This behavior is not a claim that these payments were necessarily erroneous: extra-contractual payments, bad-faith settlements, and mismatched ground-up definitions require specific analyst investigation, not silent censoring.

bad <- dat
bad$ground_up_loss[2] <- bad$censor_limit[2] + 1000
# Use a guaranteed finite limit in the example.
bad$censor_limit[2] <- 10000
bad$ground_up_loss[2] <- 10500
fit <- suppressWarnings(fit_severity(bad, loss = "ground_up_loss",
  limit = "censor_limit", distribution = "lognormal", hessian = FALSE))
fit$excluded_rows
#> [1] 2
fit$excluded_details
#>   row                     reason
#> 1   2 loss above censoring limit

Optimization diagnostics

Inspect the candidate table for outlying runtimes or optimizer evaluation counts. When fitting scale_intercept = TRUE, the optimizer automatically scales its intercept parameter; it does not constrain its value. A caller can still supply control$parscale explicitly. The displayed fit coefficients remain in their original units.

candidates$comparison[, c("distribution", "elapsed_seconds",
  "function_evaluations", "gradient_evaluations", "convergence_code")]
#>   distribution elapsed_seconds function_evaluations gradient_evaluations
#> 1    lognormal            0.03                   25                    5
#> 2        gamma            0.22                  103                   29
#> 3      weibull            0.06                   35                   11
#>   convergence_code
#> 1                0
#> 2                0
#> 3                0