A few years ago, I was reviewing a risk model written by a junior analyst who had a background in general data science but had never touched actuarial work. He had built a model to predict insurance claim costs for a fleet of delivery vans using a standard linear regression. When I looked at the residuals, I started laughing—not because it was funny, but because it was a classic mistake. The model was predicting negative claim amounts for some vehicles, and it was completely ignoring the "black swan" events—those three or four massive accidents that accounted for 80% of the total loss. He was treating loss data like it followed a bell curve, but loss data is almost never Gaussian; it's skewed, it's non-negative, and it has a heavy tail.

Loss modeling is essentially the art of breaking a problem into two distinct pieces: frequency and severity. You don't just ask "how much money will we lose?" Instead, you ask "how often will a loss occur?" (frequency) and "given that a loss occurred, how big will it be?" (severity). Once you have these two components, you can combine them to find the aggregate loss. In R, this means moving away from lm() and embracing Generalized Linear Models (GLMs) with distributions that actually fit the reality of risk.

Modeling Frequency with Poisson and Negative Binomial

When we talk about frequency, we're dealing with counts. You can't have 2.5 accidents; you have 2 or 3. This is why we typically start with the Poisson distribution. In R, you'll use glm() with family = poisson. It's straightforward, but there's a catch: the Poisson distribution assumes the mean is equal to the variance. In the real world, insurance data is often "overdispersed," meaning the variance is much higher than the mean because some drivers are just significantly more prone to accidents than others.

If you find that your Poisson model is underestimating the variance, I recommend switching to a Negative Binomial distribution. While R's base glm() doesn't handle Negative Binomial out of the box, the MASS package provides glm.nb(), which is the industry standard for handling that extra volatility in frequency data.

# Quick look at a Poisson frequency model
freq_model <- glm(claim_count ~ driver_age + vehicle_type, 
                  data = claims_data, 
                  family = poisson(link = "log"))

# If overdispersion is present, we shift to Negative Binomial
library(MASS)
freq_nb <- glm.nb(claim_count ~ driver_age + vehicle_type, data = claims_data)

Estimating Severity with Gamma and Lognormal

Severity is where things get tricky. Loss amounts are always positive and usually heavily right-skewed. If you use a normal distribution, you're effectively telling the model that a $10,000 claim is just as likely to be "offset" by a -$10,000 claim, which is physically impossible. To handle this, we use the Gamma distribution or the Lognormal distribution.

The Gamma distribution is particularly useful because it's flexible and works well within the glm() framework. When using family = Gamma(link = "log"), you're modeling the expected cost per claim. I prefer the log link because it ensures the predicted losses remain positive, regardless of the coefficients. You'll notice that your predictions will now properly account for those few high-cost claims without letting them pull the entire average into an unrealistic space.

# Modeling the cost per claim (severity)
# We filter for only those who actually had a claim
severity_data <- subset(claims_data, claim_count > 0)

sev_model <- glm(claim_amount ~ vehicle_age + coverage_level, 
                 data = severity_data, 
                 family = Gamma(link = "log"))

Calculating the Aggregate Loss

The real magic happens when you combine these two. The aggregate loss is the product of the predicted frequency and the predicted severity. If your frequency model says a driver is expected to have 0.15 claims per year, and your severity model says the average claim cost is $2,000, your expected annual loss for that driver is $300.

In a production pipeline, you'll typically generate predictions for both models across your entire portfolio and multiply them element-wise. This approach allows you to see exactly where your risk is coming from. Is the loss driven by a high volume of small claims (a frequency problem) or a low volume of catastrophic claims (a severity problem)? Being able to distinguish between the two is what separates a software engineer who knows R from a true risk modeler.

📋 Practical Task

Simulating an Aggregate Loss Distribution for a Small Fleet

You have been tasked with estimating the potential total loss for a fleet of 100 delivery vans over the next year. You don't have a full dataset, but your actuarial team has given you the following parameters:

• Frequency: The number of claims per van follows a Poisson distribution with a mean ($\lambda$) of 0.2.
• Severity: The cost of each individual claim follows a Gamma distribution with a shape ($\alpha$) of 2 and a scale ($\theta$) of 1500.

Write an R script that:

1. Simulates the number of claims for each of the 100 vans using rpois().
2. For every claim generated, simulates a cost using rgamma().
3. Calculates the total aggregate loss for the entire fleet.
4. Repeats this simulation 10,000 times to create a distribution of total annual losses.
5. Calculates the 95th percentile of this distribution (the Value at Risk) to determine how much capital the company needs to hold to cover losses in 95% of years.