Actuarial Severity Fitting for Limits and Deductibles

Greg Sollenberger

What this package assumes

actuarialfitdist is designed around a claim-level ground-up severity likelihood.

This assumption is important enough to state repeatedly:

For an exact observation, the row likelihood is proportional to f(loss) / S(deductible). For a limit-censored observation it is S(limit) / S(deductible).

If a row has likelihood weight w, its log-likelihood contribution is w times the ordinary row log-likelihood. An integer weight of 2 therefore has exactly the same likelihood effect as repeating that claim row twice.

1. Auto property-damage limit analysis

Suppose an insurer is analyzing property-damage severity. There is no deductible in this example, but policies have different limits. A claim that reaches its policy limit is recorded at the limit and is right-censored.

n <- 1200
true_loss <- rlnorm(n, meanlog = log(9000), sdlog = 0.95)
policy_limit <- sample(c(10000, 25000, 50000, 100000), n, replace = TRUE,
                       prob = c(0.15, 0.35, 0.35, 0.15))

pd <- data.frame(
  loss = pmin(true_loss, policy_limit),
  deductible = 0,
  limit = policy_limit
)

mean(pd$loss >= pd$limit)
#> [1] 0.1341667

Fit several plausible severity families:

pd_candidates <- fit_severity_candidates(
  pd,
  loss = "loss",
  limit = "limit",
  distributions = c("lognormal", "weibull", "gamma", "burr", "loglogistic")
)

pd_candidates$comparison
#>   distribution converged    logLik      AIC      BIC parameters    n n_censored
#> 1    lognormal      TRUE -10843.45 21690.90 21701.08          2 1200        161
#> 2      weibull      TRUE -10915.10 21834.20 21844.38          2 1200        161
#> 3        gamma      TRUE -10897.90 21799.79 21809.98          2 1200        161
#> 4         burr      TRUE -10841.99 21689.99 21705.26          3 1200        161
#> 5  loglogistic      TRUE -10843.38 21690.75 21700.93          2 1200        161
#>   error elapsed_seconds function_evaluations gradient_evaluations
#> 1  <NA>            0.03                   20                    6
#> 2  <NA>            0.06                   35                   11
#> 3  <NA>            0.18                  104                   26
#> 4  <NA>            0.29                  139                   39
#> 5  <NA>            0.05                   22                    6
#>   convergence_code
#> 1                0
#> 2                0
#> 3                0
#> 4                0
#> 5                0

The package does not automatically choose one distribution. The statistical fit can be compared with the business output directly:

pd_factors <- severity_factors(
  pd_candidates,
  deductible = 0,
  limit = c(10000, 25000, 50000, 100000, 250000),
  base_deductible = 0,
  base_limit = 25000
)

head(pd_factors)
#>             distribution  model scale_value deductible  limit valuation_factor
#> lognormal.1    lognormal single          NA          0  10000                1
#> lognormal.2    lognormal single          NA          0  25000                1
#> lognormal.3    lognormal single          NA          0  50000                1
#> lognormal.4    lognormal single          NA          0 100000                1
#> lognormal.5    lognormal single          NA          0 250000                1
#> weibull.1        weibull single          NA          0  10000                1
#>             expected_payment limited_expected_severity unlimited_mean LER
#> lognormal.1         7403.710                  7403.710       14275.64   0
#> lognormal.2        11349.075                 11349.075       14275.64   0
#> lognormal.3        13230.395                 13230.395       14275.64   0
#> lognormal.4        14023.143                 14023.143       14275.64   0
#> lognormal.5        14255.733                 14255.733       14275.64   0
#> weibull.1           7436.863                  7436.863       13367.07   0
#>                   ILF rating_factor base_deductible base_limit
#> lognormal.1 0.6523624     0.6523624               0      25000
#> lognormal.2 1.0000000     1.0000000               0      25000
#> lognormal.3 1.1657685     1.1657685               0      25000
#> lognormal.4 1.2356199     1.2356199               0      25000
#> lognormal.5 1.2561140     1.2561140               0      25000
#> weibull.1   0.6289921     0.6289921               0      25000
#>             base_valuation_factor
#> lognormal.1                     1
#> lognormal.2                     1
#> lognormal.3                     1
#> lognormal.4                     1
#> lognormal.5                     1
#> weibull.1                       1

The factor table is intentionally long-format. compare_factors() provides a side-by-side view when desired.

compare_factors(
  pd_candidates,
  metric = "rating_factor",
  deductible = 0,
  limit = c(10000, 25000, 50000, 100000),
  base_deductible = 0,
  base_limit = 25000
)
#>             scale_value deductible limit valuation_factor base_deductible
#> lognormal.1          NA          0 10000                1               0
#>             base_limit base_valuation_factor value.lognormal value.weibull
#> lognormal.1      25000                     1       0.6523624     0.6289921
#>             value.gamma value.burr value.loglogistic
#> lognormal.1   0.6394561  0.6637499         0.6652174

2. Auto collision deductible analysis

Now consider collision severity with varying deductibles and no finite policy limit. Ground-up losses below the deductible never appear in the claim sample.

n <- 1500
collision_ground_up <- rgamma(n, shape = 2.2, scale = 3200)
collision_ded <- sample(c(250, 500, 1000, 2000), n, replace = TRUE,
                        prob = c(0.10, 0.35, 0.40, 0.15))
keep <- collision_ground_up > collision_ded

collision <- data.frame(
  loss = collision_ground_up[keep],
  deductible = collision_ded[keep],
  limit = Inf
)

The ground-up amount is retained. We do not subtract the deductible before fitting.

collision_fit <- fit_severity(
  collision,
  loss = "loss",
  deductible = "deductible",
  limit = "limit",
  distribution = "gamma"
)

summary(collision_fit)
#> Actuarial severity model summary
#> 
#> Model: single
#> Distribution: gamma
#> Observations: 1468 (effective weighted count: 1468)
#> Right-censored: 0
#> LogLik: -14225.1  AIC: 28454.2  BIC: 28464.78
#> 
#>          estimate    std.error
#> scale 3232.648535 140.23679450
#> shape    2.206463   0.09288173
#> 
#> The internal scale is the Gamma scale parameter; with common shape, mean severity is proportional to this scale.

Generate deductible factors relative to a $500 deductible:

severity_factors(
  collision_fit,
  deductible = c(250, 500, 1000, 2000, 5000),
  limit = Inf,
  base_deductible = 500,
  base_limit = Inf
)
#>   distribution  model scale_value deductible limit valuation_factor
#> 1        gamma single          NA        250   Inf                1
#> 2        gamma single          NA        500   Inf                1
#> 3        gamma single          NA       1000   Inf                1
#> 4        gamma single          NA       2000   Inf                1
#> 5        gamma single          NA       5000   Inf                1
#>   expected_payment limited_expected_severity unlimited_mean        LER ILF
#> 1         6882.829                  7132.721       7132.721 0.03503457   1
#> 2         6633.681                  7132.721       7132.721 0.06996494   1
#> 3         6140.903                  7132.721       7132.721 0.13905188   1
#> 4         5197.377                  7132.721       7132.721 0.27133320   1
#> 5         2919.909                  7132.721       7132.721 0.59063188   1
#>   rating_factor base_deductible base_limit base_valuation_factor
#> 1     1.0375581             500        Inf                     1
#> 2     1.0000000             500        Inf                     1
#> 3     0.9257158             500        Inf                     1
#> 4     0.7834831             500        Inf                     1
#> 5     0.4401642             500        Inf                     1

3. Homeowners deductible analysis with varying Coverage A

For homeowners, claim severity can vary systematically with the insured amount. The optional scale_by column lets the fitted severity scale vary linearly with an analyst-supplied numeric variable.

Here Coverage A is the scaling variable. The same mechanism could instead use TIV, vehicle value, or a user-created transformation in another line of business.

n <- 1600
coverage_a <- sample(seq(150000, 900000, by = 25000), n, replace = TRUE)
home_ded <- sample(c(1000, 2500, 5000, 10000), n, replace = TRUE)
shape <- 1.35
scale_coef <- 0.075
latent_scale <- scale_coef * coverage_a
home_ground_up <- rweibull(n, shape = shape, scale = latent_scale)
home_limit <- coverage_a
keep <- home_ground_up > home_ded

home <- data.frame(
  loss = pmin(home_ground_up[keep], home_limit[keep]),
  deductible = home_ded[keep],
  limit = home_limit[keep],
  coverage_a = coverage_a[keep]
)

Fit a Weibull model with zero intercept:

home_fit <- fit_severity(
  home,
  loss = "loss",
  deductible = "deductible",
  limit = "limit",
  distribution = "weibull",
  scale_by = "coverage_a",
  scale_intercept = FALSE
)

summary(home_fit)
#> Actuarial severity model summary
#> 
#> Model: single
#> Distribution: weibull
#> Observations: 1475 (effective weighted count: 1475)
#> Right-censored: 0
#> LogLik: -16679.47  AIC: 33362.93  BIC: 33373.53
#> 
#>              estimate   std.error
#> scale_coef 0.07416593 0.001680669
#> shape      1.33880873 0.033190229
#> 
#> The internal scale is the Weibull scale parameter; with common shape, mean severity is proportional to this scale.

When scale_intercept = FALSE, the intercept is exactly zero, not merely shrunk toward zero or omitted from printed output.

For the Weibull family, the native scale is scale_coef * coverage_a, so the fitted mean is that quantity times Gamma(1 + 1 / shape). Other supported families use parameterizations that preserve the same intended proportional-mean relationship when their mean is finite.

Factor output can be produced for arbitrary Coverage A values:

severity_factors(
  home_fit,
  scale_value = c(200000, 400000, 600000, 800000),
  deductible = c(1000, 2500, 5000),
  limit = c(100000, 250000, 500000),
  base_deductible = 1000,
  base_limit = 250000
)
#>    distribution  model scale_value deductible  limit valuation_factor
#> 1       weibull single       2e+05       1000 100000                1
#> 2       weibull single       4e+05       1000 100000                1
#> 3       weibull single       6e+05       1000 100000                1
#> 4       weibull single       8e+05       1000 100000                1
#> 5       weibull single       2e+05       2500 100000                1
#> 6       weibull single       4e+05       2500 100000                1
#> 7       weibull single       6e+05       2500 100000                1
#> 8       weibull single       8e+05       2500 100000                1
#> 9       weibull single       2e+05       5000 100000                1
#> 10      weibull single       4e+05       5000 100000                1
#> 11      weibull single       6e+05       5000 100000                1
#> 12      weibull single       8e+05       5000 100000                1
#> 13      weibull single       2e+05       1000 250000                1
#> 14      weibull single       4e+05       1000 250000                1
#> 15      weibull single       6e+05       1000 250000                1
#> 16      weibull single       8e+05       1000 250000                1
#> 17      weibull single       2e+05       2500 250000                1
#> 18      weibull single       4e+05       2500 250000                1
#> 19      weibull single       6e+05       2500 250000                1
#> 20      weibull single       8e+05       2500 250000                1
#> 21      weibull single       2e+05       5000 250000                1
#> 22      weibull single       4e+05       5000 250000                1
#> 23      weibull single       6e+05       5000 250000                1
#> 24      weibull single       8e+05       5000 250000                1
#> 25      weibull single       2e+05       1000 500000                1
#> 26      weibull single       4e+05       1000 500000                1
#> 27      weibull single       6e+05       1000 500000                1
#> 28      weibull single       8e+05       1000 500000                1
#> 29      weibull single       2e+05       2500 500000                1
#> 30      weibull single       4e+05       2500 500000                1
#> 31      weibull single       6e+05       2500 500000                1
#> 32      weibull single       8e+05       2500 500000                1
#> 33      weibull single       2e+05       5000 500000                1
#> 34      weibull single       4e+05       5000 500000                1
#> 35      weibull single       6e+05       5000 500000                1
#> 36      weibull single       8e+05       5000 500000                1
#>    expected_payment limited_expected_severity unlimited_mean        LER
#> 1         12633.778                  13622.32       13622.33 0.07256755
#> 2         26168.293                  27157.77       27244.66 0.03653727
#> 3         38690.383                  39636.75       40866.99 0.02440477
#> 4         49076.014                  49942.21       54489.32 0.01831911
#> 5         11218.044                  13622.32       13622.33 0.17649523
#> 6         24710.469                  27157.77       27244.66 0.09034749
#> 7         37282.994                  39636.75       40866.99 0.06062379
#> 8         47780.975                  49942.21       54489.32 0.04559916
#> 9          9085.913                  13622.32       13622.33 0.33301285
#> 10        22375.348                  27157.77       27244.66 0.17649523
#> 11        34985.114                  39636.75       40866.99 0.11959269
#> 12        45644.876                  49942.21       54489.32 0.09034749
#> 13        12633.791                  13622.33       13622.33 0.07256755
#> 14        26249.214                  27244.66       27244.66 0.03653727
#> 15        39868.923                  40866.23       40866.99 0.02440477
#> 16        53464.485                  54461.65       54489.32 0.01831911
#> 17        11218.054                  13622.33       13622.33 0.17649523
#> 18        24783.173                  27244.66       27244.66 0.09034749
#> 19        38388.818                  40866.23       40866.99 0.06062379
#> 20        51979.490                  54461.65       54489.32 0.04559916
#> 21         9085.919                  13622.33       13622.33 0.33301285
#> 22        22436.107                  27244.66       27244.66 0.17649523
#> 23        35979.022                  40866.23       40866.99 0.11959269
#> 24        49543.471                  54461.65       54489.32 0.09034749
#> 25        12633.791                  13622.33       13622.33 0.07256755
#> 26        26249.214                  27244.66       27244.66 0.03653727
#> 27        39869.640                  40866.99       40866.99 0.02440477
#> 28        53491.123                  54489.32       54489.32 0.01831911
#> 29        11218.054                  13622.33       13622.33 0.17649523
#> 30        24783.173                  27244.66       27244.66 0.09034749
#> 31        38389.478                  40866.99       40866.99 0.06062379
#> 32        52004.651                  54489.32       54489.32 0.04559916
#> 33         9085.919                  13622.33       13622.33 0.33301285
#> 34        22436.107                  27244.66       27244.66 0.17649523
#> 35        35979.596                  40866.99       40866.99 0.11959269
#> 36        49566.345                  54489.32       54489.32 0.09034749
#>          ILF rating_factor base_deductible base_limit base_valuation_factor
#> 1  0.9999990     0.9999990            1000     250000                     1
#> 2  0.9969172     0.9969172            1000     250000                     1
#> 3  0.9704396     0.9704396            1000     250000                     1
#> 4  0.9179180     0.9179180            1000     250000                     1
#> 5  0.9999992     0.8879397            1000     250000                     1
#> 6  0.9970664     0.9413794            1000     250000                     1
#> 7  0.9711941     0.9351392            1000     250000                     1
#> 8  0.9192275     0.8936956            1000     250000                     1
#> 9  0.9999993     0.7191755            1000     250000                     1
#> 10 0.9972919     0.8524197            1000     250000                     1
#> 11 0.9723753     0.8775034            1000     250000                     1
#> 12 0.9213096     0.8537420            1000     250000                     1
#> 13 1.0000000     1.0000000            1000     250000                     1
#> 14 1.0000000     1.0000000            1000     250000                     1
#> 15 1.0000000     1.0000000            1000     250000                     1
#> 16 1.0000000     1.0000000            1000     250000                     1
#> 17 1.0000000     0.8879404            1000     250000                     1
#> 18 1.0000000     0.9441491            1000     250000                     1
#> 19 1.0000000     0.9628757            1000     250000                     1
#> 20 1.0000000     0.9722246            1000     250000                     1
#> 21 1.0000000     0.7191760            1000     250000                     1
#> 22 1.0000000     0.8547344            1000     250000                     1
#> 23 1.0000000     0.9024328            1000     250000                     1
#> 24 1.0000000     0.9266613            1000     250000                     1
#> 25 1.0000000     1.0000000            1000     250000                     1
#> 26 1.0000000     1.0000000            1000     250000                     1
#> 27 1.0000180     1.0000180            1000     250000                     1
#> 28 1.0004982     1.0004982            1000     250000                     1
#> 29 1.0000000     0.8879404            1000     250000                     1
#> 30 1.0000000     0.9441492            1000     250000                     1
#> 31 1.0000172     0.9628923            1000     250000                     1
#> 32 1.0004841     0.9726953            1000     250000                     1
#> 33 1.0000000     0.7191760            1000     250000                     1
#> 34 1.0000000     0.8547344            1000     250000                     1
#> 35 1.0000160     0.9024471            1000     250000                     1
#> 36 1.0004617     0.9270892            1000     250000                     1

Scaling diagnostics

A severity histogram across the entire Coverage A range is not a comparison to one fitted distribution: each Coverage A value implies a different fitted severity distribution. fit_plot_data() can therefore restrict the comparison to a meaningful scale range while still retaining the exact row-level scale, deductible, and limit values inside that range.

plot_fit(
  home_fit,
  type = "density",
  scale_range = c(300000, 450000)
)

A separate scaling diagnostic asks whether the assumed relationship across Coverage A is reasonable.

plot_scaling_fit(home_fit, bins = 10, nsim = 3, seed = 12)

The actual and modeled observed means are directly comparable because the modeled series reproduces each claim row’s deductible and censoring structure. The fitted ground-up mean series is also shown to make the underlying severity-scaling relationship visible.

4. Spliced body-tail fitting

A single family may fit the body of a severity distribution reasonably while missing the upper tail. fit_spliced_severity() supports a body family below a threshold and a separate tail family above it.

The two densities are normalized on their respective sides of the threshold. A body probability is estimated. Density continuity is not required.

The threshold may be fixed by the analyst or searched over candidate values.

n <- 1500
threshold_true <- 30000
is_body <- runif(n) < 0.90

# Draw body values from a lognormal conditional on being below the threshold.
Fb <- plnorm(threshold_true, log(9000), 0.75)
body_draw <- qlnorm(runif(n, 0, Fb), log(9000), 0.75)

# Draw Pareto-II tail values conditional on exceeding the threshold.
shape_tail <- 2.1
scale_tail <- 18000
Ft <- 1 - (1 + threshold_true / scale_tail)^(-shape_tail)
u_tail <- runif(n, Ft, 1)
tail_draw <- scale_tail * ((1 - u_tail)^(-1 / shape_tail) - 1)

splice_ground_up <- ifelse(is_body, body_draw, tail_draw)
splice_ded <- sample(c(0, 1000, 2500), n, replace = TRUE)
splice_lim <- sample(c(50000, 100000, 250000, Inf), n, replace = TRUE)
keep <- splice_ground_up > splice_ded

splice_dat <- data.frame(
  loss = pmin(splice_ground_up[keep], splice_lim[keep]),
  deductible = splice_ded[keep],
  limit = splice_lim[keep]
)

Search for a threshold:

splice_fit <- fit_spliced_severity(
  splice_dat,
  loss = "loss",
  deductible = "deductible",
  limit = "limit",
  body = "lognormal",
  tail = "pareto2",
  threshold_method = "search",
  threshold_probs = c(0.60, 0.95),
  threshold_grid = 4
)
#> Warning: Spliced-fit diagnostic: Extreme tail scale and shape: possible weak
#> identifiability

splice_fit
#> Actuarial severity fit
#> Model:        spliced
#> Distribution: lognormal+pareto2
#> Threshold:     50000 (search)
#> Ground-up observations: 1476
#> Right-censored:          23
#> Log-likelihood:          -15243.51
#> AIC:                     30499.02
#> Converged:               yes
#> 
#> Coefficients:
#>   body_scale   body_sdlog   tail_scale   tail_shape  splice_prob 
#> 9.400326e+03 8.157510e-01 2.381546e+07 5.053511e+02 9.494472e-01 
#> 
#> Note: Spliced densities are normalized on each side of the selected threshold; density continuity at the threshold is not imposed.
splice_fit$threshold_search
#>   threshold n_body n_exact_tail n_censored    logLik      AIC      BIC
#> 1  12192.85    886          567         23 -15249.63 30509.27 30535.75
#> 2  16508.30   1058          395         23 -15246.43 30502.86 30529.35
#> 3  22961.79   1230          223         23 -15248.76 30507.53 30534.01
#> 4  50000.00   1414           53          9 -15243.51 30497.02 30523.50
#>   converged reason                                                tail_warning
#> 1      TRUE   <NA>                                                        <NA>
#> 2      TRUE   <NA>                                                        <NA>
#> 3      TRUE   <NA>                                                        <NA>
#> 4      TRUE   <NA> Extreme tail scale and shape: possible weak identifiability

The same factor-generation interface works for the spliced model:

severity_factors(
  splice_fit,
  deductible = c(0, 1000, 5000),
  limit = c(25000, 50000, 100000, 250000),
  base_deductible = 0,
  base_limit = 50000
)
#>         distribution   model scale_value deductible  limit valuation_factor
#> 1  lognormal+pareto2 spliced          NA          0  25000                1
#> 2  lognormal+pareto2 spliced          NA       1000  25000                1
#> 3  lognormal+pareto2 spliced          NA       5000  25000                1
#> 4  lognormal+pareto2 spliced          NA          0  50000                1
#> 5  lognormal+pareto2 spliced          NA       1000  50000                1
#> 6  lognormal+pareto2 spliced          NA       5000  50000                1
#> 7  lognormal+pareto2 spliced          NA          0 100000                1
#> 8  lognormal+pareto2 spliced          NA       1000 100000                1
#> 9  lognormal+pareto2 spliced          NA       5000 100000                1
#> 10 lognormal+pareto2 spliced          NA          0 250000                1
#> 11 lognormal+pareto2 spliced          NA       1000 250000                1
#> 12 lognormal+pareto2 spliced          NA       5000 250000                1
#>    expected_payment limited_expected_severity unlimited_mean        LER
#> 1         11814.759                  11814.76       16243.18 0.00000000
#> 2         10953.495                  11814.76       16243.18 0.06152844
#> 3          7781.570                  11814.76       16243.18 0.28611659
#> 4         13851.063                  13851.06       16243.18 0.00000000
#> 5         12901.667                  13851.06       16243.18 0.06152844
#> 6          9443.465                  13851.06       16243.18 0.28611659
#> 7         15410.717                  15410.72       16243.18 0.00000000
#> 8         14428.672                  15410.72       16243.18 0.06152844
#> 9         10846.576                  15410.72       16243.18 0.28611659
#> 10        16207.626                  16207.63       16243.18 0.00000000
#> 11        15208.946                  16207.63       16243.18 0.06152844
#> 12        11563.720                  16207.63       16243.18 0.28611659
#>          ILF rating_factor base_deductible base_limit base_valuation_factor
#> 1  0.8529858     0.8529858               0      50000                     1
#> 2  0.8489984     0.7908054               0      50000                     1
#> 3  0.8240165     0.5618031               0      50000                     1
#> 4  1.0000000     1.0000000               0      50000                     1
#> 5  1.0000000     0.9314568               0      50000                     1
#> 6  1.0000000     0.6817863               0      50000                     1
#> 7  1.1126018     1.1126018               0      50000                     1
#> 8  1.1183572     1.0417014               0      50000                     1
#> 9  1.1485802     0.7830862               0      50000                     1
#> 10 1.1701359     1.1701359               0      50000                     1
#> 11 1.1788358     1.0980346               0      50000                     1
#> 12 1.2245210     0.8348616               0      50000                     1

Valuation or depreciation at the payment stage

The fitting layer does not decide whether historical claims are RCV, ACV, or some other valuation basis. The analyst should first express the experience on a consistent ground-up basis.

At factor generation, a multiplicative valuation factor can then represent an ACV-type payment transformation. If the ground-up replacement-cost severity is X, the package evaluates payments of the form

min(max(valuation_factor * X - deductible, 0), limit)

For example:

severity_factors(
  home_fit,
  scale_value = 400000,
  deductible = 2500,
  limit = 25000,
  valuation_factor = c(1.00, 0.80, 0.60, 0.40, 0.20),
  base_deductible = 2500,
  base_limit = 25000,
  base_valuation_factor = 1
)
#>   distribution  model scale_value deductible limit valuation_factor
#> 1      weibull single       4e+05       2500 25000              1.0
#> 2      weibull single       4e+05       2500 25000              0.8
#> 3      weibull single       4e+05       2500 25000              0.6
#> 4      weibull single       4e+05       2500 25000              0.4
#> 5      weibull single       4e+05       2500 25000              0.2
#>   expected_payment limited_expected_severity unlimited_mean        LER ILF
#> 1        16897.766                  18289.18       27244.66 0.09034749   1
#> 2        14954.852                  18289.18       27244.66 0.11232866   1
#> 3        12183.568                  18289.18       27244.66 0.14831883   1
#> 4         8238.316                  18289.18       27244.66 0.21767889   1
#> 5         3252.768                  18289.18       27244.66 0.40285079   1
#>   rating_factor base_deductible base_limit base_valuation_factor
#> 1     1.0000000            2500      25000                     1
#> 2     0.8850194            2500      25000                     1
#> 3     0.7210165            2500      25000                     1
#> 4     0.4875388            2500      25000                     1
#> 5     0.1924970            2500      25000                     1

Weights and fixed parameters

A likelihood weight multiplies the row log-likelihood. Integer weights therefore behave as frequency weights.

weighted_collision <- collision[1:1000, ]
weighted_collision$w <- sample(1:3, nrow(weighted_collision), replace = TRUE)

weighted_fit <- fit_severity(
  weighted_collision,
  loss = "loss",
  deductible = "deductible",
  limit = "limit",
  distribution = "gamma",
  weights = "w",
  fixed = list(shape = 2.2)
)

coef(weighted_fit)
#>   scale   shape 
#> 3223.05    2.20

Hessian and bootstrap uncertainty

vcov() and confint() use a numerical Hessian of the optimized negative log-likelihood. These are large-sample Wald approximations; the Hessian is evaluated numerically rather than supplied as hand-derived second derivatives for every family.

vcov(collision_fit)
#>             scale         shape
#> scale 19666.35853 -11.885324597
#> shape   -11.88532   0.008627015
confint(collision_fit)
#>             lower       upper
#> scale 2957.789469 3507.507602
#> shape    2.024419    2.388508

A lightweight parametric or case bootstrap is also available. Parametric bootstrap retains the fitted portfolio’s deductible, limit, scaling, and weight structure.

boot <- bootstrap_fit(collision_fit, R = 5, type = "parametric", seed = 99)
head(boot$estimates)
boot$success_rate

For a searched spliced fit, refit_threshold = FALSE reuses the selected threshold in each bootstrap replicate by default. Set it to TRUE to repeat threshold search, understanding that this is much more computationally expensive.

Custom ggplot2 and Plotly graphics

The plot-building data are part of the public interface:

d <- fit_plot_data(home_fit, scale_range = c(300000, 450000), nsim = 3, seed = 5)
head(d)
#>   row source simulation observed_loss censored deductible  limit scale_value
#> 1   1 actual          0      17733.91    FALSE       2500 450000      450000
#> 2   2 actual          0      26398.36    FALSE      10000 350000      350000
#> 3   3 actual          0      48168.44    FALSE      10000 375000      375000
#> 4   4 actual          0      19954.84    FALSE      10000 400000      400000
#> 5   5 actual          0      34167.79    FALSE       1000 400000      400000
#> 6   6 actual          0      15006.80    FALSE       2500 350000      350000
#>   weight
#> 1      1
#> 2      1
#> 3      1
#> 4      1
#> 5      1
#> 6      1

This allows analysts to use their own graphics conventions without reproducing the statistically important observation logic. See inst/examples/08_plotly_diagnostics.R for Plotly examples.

Optional terms, anomaly audits, and log-scale plots

The policy terms are optional. A sample with no deductible or censoring limit can be fitted by naming only the loss column and family. In particular, Coverage A is not a policy censoring limit merely because it is an insured value.

set.seed(20261008)
no_terms <- data.frame(loss = rlnorm(250, log(12000), .8))
fit_no_terms <- fit_severity(no_terms, loss = "loss",
  distribution = "lognormal", hessian = FALSE)
print(fit_no_terms)
#> Actuarial severity fit
#> Model:        single
#> Distribution: lognormal
#> Ground-up observations: 250
#> Right-censored:          0
#> Log-likelihood:          -2624.113
#> AIC:                     5252.225
#> Converged:               yes
#> 
#> Coefficients:
#>        scale        sdlog 
#> 1.152856e+04 7.595185e-01 
#> 
#> Note: The internal scale is exp(meanlog). With scale_by, mean severity is proportional to scale_coef * scale_by + scale_intercept.
plot_fit(fit_no_terms, x_scale = "log10",
  xlim = c(1000, 100000))

If a supplied row is inconsistent with the observation mechanism, the default is to exclude the row after warning; the original positional row number is reported. This is an audit aid, not a substitute for understanding how the original ground-up amounts were collected.

problem <- no_terms
problem$deductible <- 0
problem$deductible[3] <- problem$loss[3] + 10
clean_fit <- suppressWarnings(fit_severity(problem, loss = "loss",
  deductible = "deductible", distribution = "lognormal", hessian = FALSE))
clean_fit$excluded_rows
#> [1] 3
clean_fit$excluded_details
#>   row                   reason
#> 1   3 loss at/below deductible

For single-fit plots and multi-candidate plots, use x_scale = "log10" by default; xlim zooms without discarding data during statistic calculation. The histogram density is computed in log-severity coordinates when log-scaled.

plot_candidate_fits(pd_candidates, type = "density",
  x_scale = "log10", xlim = c(200, 150000))