---
title: "Traditional-style reserving using GLMs"
author: "Grainne McGuire"
date: '2019-11-07'
description: "This post provides a worked example in R of fitting a GLM to some non-life claims reserving data"
categories: ["GLM", "R", "Reserving"]
image: "heatmap.png"
resources: 
  - "_traditional_style_reserving_glms.qmd" 
  - "_glms_meyersshi.csv"
---


## April 2023 update

The article below is my original article working through fitting a GLM to aggregate triangular data using the methods and data from the CAS Monograph [Stochastic Loss Reserving using Generalized Linear Models](https://www.casact.org/monograph/cas-monograph-no-3).

Since then some or all of the article has re-published with some changes at:

- [Reserving with GLMs on the MLRWP blog](https://institute-and-faculty-of-actuaries.github.io/mlr-blog/post/foundations/r-glms/)
- [Reserving with GLMs in Python on the MLRWP blog](https://institute-and-faculty-of-actuaries.github.io/mlr-blog/post/foundations/python-glms/)


## Introduction

This post provides a worked example in R of fitting a GLM to some non-life claims reserving data. 
The example and data are drawn from the CAS Monograph [Stochastic Loss Reserving using Generalized Linear Models](/publication/2016_monograph).
Here we follow through the application of the cross-classified model from that monograph to the data (Chapter 3), and follow through with the additional work to firstly simplify the model and secondly to improve the model fit through the use of some interaction terms (Chapter 7).

- The data used in this article is available [here](_glms_meyersshi.csv)
- You can download this notebook as a quarto file to run the code yourself from [here](_traditional_style_reserving_glms.qmd)


## Setup

Before we begin, we first set up the R session by loading in the various packages that we need.

* [here](https://cran.rstudio.com/web/packages/here/index.html): I first heard about here from reading [Jenny Bryan's](https://github.com/jennybc/here_here) article on it and have been a fan of `here` and the R project structure ever since. It really helps with portability of code.
    + Basically it allows you to use relative rather than absolute file paths.
    + If you want to run this notebook and don't want to use `here` then all you need to do is put an appropriate pathname in for loading in the data from a CSV file.
      Location is not used anywhere else.
      
* [data.table](https://github.com/Rdatatable/data.table) I really like `data.table` - its power, speed and terseness. At some point though I may replace the `data.table` code with base R to reduce dependencies. For now though, there isn't a huge amount of `data.table` code.
    + Even if you don't like `data.table` syntax, the `fread` and `fwrite` functions can be very useful for reading and writing CSV files.

* [ggplot2](https://ggplot2.tidyverse.org/): create nice graphs easily
    + [viridis](https://cran.r-project.org/web/packages/viridis/vignettes/intro-to-viridis.html) nice colour palettes that are tested for common forms of colour-blindness
    + [cowplot](https://cran.r-project.org/web/packages/cowplot/index.html) - an easy way of grouping graphs into a single figure
    + [patchwork]https://cran.r-project.org/web/packages/patchwork/index.html) - an easy way of grouping graphs into a single figure


* [knitr](https://cran.r-project.org/web/packages/knitr/index.html) The notebook uses `kable` from the `knitr` package. If you're using RStudio to run this code in notebook format, then you should already have it. Otherwise you can install it, or you can simply replace all the `kable` commands with print statements.

If you don't have any of these packages you will need to install them via `install.packages` or, if using RStudio, via the Install buttom in the packages pane. 


```{r setup, message=FALSE, collapse=TRUE}

library(here)
library(data.table)
library(ggplot2)
library(viridis)
library(patchwork)
# I have elected not to attach knitr, so we need to use knitr::kable() below

options(scipen = 99)   # get rid of scientific notation

# use this theme in all plots
ggplot2::theme_set(theme_classic())


```



## Data

The data are presented in Section 1.3 of the monograph. Ultimately the data were sourced from the Meyers and Shi (2011) database, and are the workers compensation triangle of the New Jersey Manufacturers Group.

For ease of use, I have created a CSV file with the data which is loaded in this code chunk. As noted above I use relative path names with the `here` package. If you don't want to have a setup that works with `here`, just ensure the full pathname to the file is included in the `fread` statement below.

Once the data is loaded in, have a look at the start of it.


```{r}

msdata <- data.table::fread(here::here("post/2019-11-07-traditional-style-reserving-using-glms/_glms_meyersshi.csv"))
# if needed replace here::here("data/glms_meyershi.csv") with
# the correct path and filename of where you put the data

setDT(msdata)

print(msdata[1:6,])

```

So we have four columns:

* acc_year: accident year, numbered from 1 to 10

* dev_year: development year, also numbered from 1 to 10

* cumulative: cumulative payments to date

* incremental: incremental payments for that accident year, development year combination.



Let's look at the data visually.

First we plot the cumulative amounts in each accident year



```{r, fig.height=6, fig.width=8}

ggplot(data=msdata, aes(x=dev_year, y=cumulative, colour=as.factor(acc_year))) +
	geom_line(size=1) +
	scale_color_viridis_d(begin=0.9, end=0) + 
	ggtitle("Cumulative amounts by development year") + 
	theme(legend.position = "right", legend.title=element_blank(), legend.text=element_text(size=8))


```

Now look at the incremental amounts

```{r, fig.height=6, fig.width=8}

ggplot(data=msdata, aes(x=dev_year, y=incremental, colour=as.factor(acc_year))) +
	geom_line(size=1) +
	scale_color_viridis_d(begin=0.9, end=0) + 
	ggtitle("Incremental amounts by development year") + 
	theme(legend.position = "right", legend.title=element_blank(), legend.text=element_text(size=8))


```


The data look quite well behaved - each year seems to have a similar development pattern.



## Modelling

### Initial model

The first model applied here is the Over-dispersed Poisson (ODP) cross classified (cc) model (Sections 3.3.2 and 3.3.3 of the monograph).
This model has been shown to give the same results as the chain ladder algorithm.

To apply the model, we will use the `glm` function from the base R `stats` package. The cross-classified model requires separate levels for each of accident and development year so we first make a factor version of these variates.



```{r}

msdata[, acc_year_factor := as.factor(acc_year)
	   ][, dev_year_factor := as.factor(dev_year)
	     ][, cal_year := acc_year + dev_year - 1]

```


Now we fit the model and look at the results via `summary`.

* The family is the *quasipoisson* - this is how we fit an ODP model with `glm`.
* The link is log
* The formula is simply "incremental ~ 0 + acc_year_factor + dev_year_factor"
    + The 0 tells `glm` to fit a model without an intercept
    + We choose to do that here because then we can more easily compare the results to those in the monograph.


```{r}

glm_fit1 <- glm(data = msdata, 
	family = quasipoisson(link = "log"),
	formula = "incremental ~ 0 + acc_year_factor + dev_year_factor")


summary(glm_fit1)

```


We now extract the coefficient table in a more convenient way and append it onto the `glm_fit1` object for later use.

We will also print the table again in a nicer form to make it easier to compare to the first 3 columns of Table 3-5 of the monograph.

* If you do this, you should see that the results match.



```{r}

# save the data for later use
glm_fit1$coeff_table <- data.table(parameter = names(glm_fit1$coefficients), 
								   coeff_glm_fit1 = glm_fit1$coefficients)

# print out the table so we can compare against Table 3-5.
glm_fit1$coeff_table |> 
	knitr::kable(digits=c(0,3))

```


### Loss reserve

* We will calculate the loss reserve for this model
* This should give the same answers as the chain ladder algorithm

```{r}

# first make the lower triangle data set
ay <- NULL
dy <- NULL


for(i in 2:10){
	ay <- c(ay, rep(i, times=(i-1)))
	dy <- c(dy, (10-i+2):10)
}

futdata <- data.table(acc_year = ay, dev_year = dy)

# make factors
futdata[, cal_year := acc_year + dev_year
		][, acc_year_factor := as.factor(acc_year)
		  ][, dev_year_factor := as.factor(dev_year)]

# make the prediction and sum by acc_year
x <- predict(glm_fit1, newdata = futdata, type="response")
futdata[, incremental := x]


ocl_year <- futdata[,  lapply(.SD, sum), .SDcols=c("incremental"), by="acc_year"]

ocl_year |> 
	knitr::kable(digits=c(0, 0))


```

* As expected, this matches the results in Table 3-2 of the monograph.



### Model diagnostics

It's always important to check that a model fits the data well, so here we look at the following:

* Residual Scatterplots
   + by linear predictor
   + by accident, development and calendar years
   
* Heat map of actual vs fitted
   + In this we get the actual/fitted ratio in each (acc, dev) cell [subject to lower and upper bounds of (0.5, 2)] and then plot the colour-coded triangle of the actual/fitted values
   + heat maps are helpful to check for model fit and may help to identify missing interactions.


**Note on residuals with `glm`**

* The residuals in a glm object accessed with `$residuals` are residuals used in the model fitting algorithm.
* For diagnostic purposes, the standardised deviance residuals are usually preferable to work with.
  + These are the signed square roots of the contribution of the i*th* observation to the deviance, divided by hat matrix values.
  + The `stats::rstandard()` function may be used with glm objects to extract the standardised deviance residuals.

<br><br>

**Generating the diagnostics**

* First we prepare the data by adding the fitted values and residuals. 
  +  Because this model has a lot of parameters, there are two observations where the fitted is exactly equal to the actual -- (acc_year=1, dev_year=10) and (acc_year=10, dev_year=0).
  This is because these observations have a unique parameter.
  + The deviance calculations below return NaN (not a number) for these points, but the residual should really be 0 so this adjustment is made manually. 
  
* Also add actual/fitted ratios and the log of these (restricted to the range [log(0.5), log(2)]) - these will be used for a heatmap later. 
  + The restricted range is used to generate easier to read shadings in the heat-map, while the conversion to log means that the shading scales will be similar intensity for $x$% and $1/x$ % 


```{r}
msdata[, residuals1 := rstandard(glm_fit1)
	   ][, fitted1 := glm_fit1$fitted.values
	     ][, linear_predictor1 := log(fitted1)
	       ][, AvsF1 := incremental / fitted1
	         ][, AvsF_restricted1 := log(pmax(0.5, pmin(2,AvsF1)))]

# check for NaN residuals
msdata[is.nan(residuals1),]

# these occur where we expect them so so replace with 0
msdata[is.nan(residuals1), residuals1 := 0]


```

Look at the first 10 rows of msdata

```{r}
head(msdata, 10)
```


Now let's look at the residual scatterplots - here I use the `cowplot` package to combine all 4 graphs into one plot.

In the linear predictor scatterplot, the points are colour coded so that the lighter points belong to the earlier development years, and the darker points belong to the later ones.


```{r, fig.height=6, fig.width=9}

p1 <- ggplot(data=msdata, aes(x=linear_predictor1, y=residuals1, colour=dev_year)) +
	geom_point(size=2) +
	scale_colour_viridis(begin=0.9, end=0) +
	theme(legend.position = "none") +
	ggtitle("Linear predictor")


p2 <- ggplot(data=msdata, aes(x=acc_year, y=residuals1)) +
	geom_point(size=2, colour="#2d708eff") +
	ggtitle("Accident year")

p3 <- ggplot(data=msdata, aes(x=dev_year, y=residuals1)) +
	geom_point(size=2, colour="#2d708eff") +
	ggtitle("Development year")

p4 <- ggplot(data=msdata, aes(x=cal_year, y=residuals1)) +
	geom_point(size=2, colour="#2d708eff") +
	ggtitle("Calendar year")

# combine plots with patchwork
p <- (p1 + p2) / (p3 + p4)

p

```


Now construct and draw the heat map. Note that the colours are:

* blue (A/F = 50%)
* white (A/F = 100%)
* red (A/F = 200%)


```{r}

# heatmap code
# to get the correct shading I've plotted the log of the restricted A/F values

p_hm <- ggplot(data=msdata, aes(x=dev_year, y=acc_year)) + 
	geom_tile(aes(fill = AvsF_restricted1))+scale_y_reverse()+
	scale_fill_gradient2(name="AvF_min", low="royalblue", mid="white", high="red", midpoint=0, space="Lab", na.value="grey50", guide="colourbar")+
	labs(x="Development year", y="Accident year")+
	theme(legend.position = "none")+
	theme(axis.title.x = element_text(size=8), axis.text.x  = element_text(size=7))+
	theme(axis.title.y = element_text(size=8), axis.text.y  = element_text(size=7))+
	theme(panel.background = element_rect(fill = "grey", colour = "grey", size = 2, linetype = "solid"),
		  panel.grid = element_line(colour="grey")) + 
	NULL

print(p_hm)


```



## Refining the model

We could stop here - and just use the results from this model, which match those produced by the chain ladder. The diagnostics suggest that the model fits quite well.
However can we:

* identify simplifications to the model to make it more parsinomious?
* identify any areas of poorer fit than may suggest missing model terms including interactions?


### Simplifying the model

First we consider if we can use a parametric shape for the accident and development year parameters.


### Accident year

First plot the accident year parameters.


```{r}

# extract the data
dt_acc_year <- glm_fit1$coeff_table[grepl("acc_year", parameter),  
									][, acc_year := as.integer(gsub("acc_year_factor", "", parameter))]


# plot
ggplot(data=dt_acc_year, aes(x=acc_year, y=coeff_glm_fit1)) +
	geom_line(size=2, colour="#440154ff") +
	geom_point(size=4, colour="#440154ff") + 
	ggtitle("Accident year parameters")


```

* Note that their shape closely resembles that of a parabola.
* This suggests that we can replace the 10 accident year parameters by
    + the overall intercept
    + an acc_year term
    + an acc_year squarted term

* So refit the model on this basis.
    + Drop the 0 from the glm_fit1 formula to allow the model to have an intecept
    + Replace the acc_year_factor term with the parabola terms.


```{r}

# add an x and x^2 term
msdata[, acc_year_2 := acc_year^2]

glm_fit2 <- glm(data = msdata, 
	family = quasipoisson(link = "log"),
	formula = "incremental ~ acc_year + acc_year_2 + dev_year_factor")


summary(glm_fit2)


```

We see in the coefficient table part of the summary that the two acc_year terms are highly significant.

<br>

Now extract the coefficients and compare the previous and current fits.

* Remember that the intercept must be included in these calculations.
* Again, save the coefficient table in this extracted form with the glm_fit2 object for later use.


```{r}
# extract the coefficient table
glm_fit2$coeff_table <- data.table(parameter = names(glm_fit2$coefficients), coeff_glm_fit2 = glm_fit2$coefficients)
print(glm_fit2$coeff_table)

```


Now compare the past and current parameter estimates for accident year.

```{r}

# pull out the acc year coefficinents only
dt_acc_year[, coeff_glm_fit2 := glm_fit2$coeff_table[parameter == "acc_year", coeff_glm_fit2]*acc_year + 
				glm_fit2$coeff_table[parameter == "acc_year_2", coeff_glm_fit2]*acc_year^2 + 
				glm_fit2$coeff_table[parameter == "(Intercept)", coeff_glm_fit2]]

# make long for ggplot
dt_acc_year_plot <- melt(dt_acc_year, id.vars = "acc_year", measure.vars = c("coeff_glm_fit1", "coeff_glm_fit2"), variable.name="model", value = "estimate")

# remove the coeff_ from the model names
dt_acc_year_plot[, model := gsub("coeff_", "", model, fixed=TRUE)]

ggplot(data=dt_acc_year_plot, aes(x=acc_year, y=estimate, colour=model)) +
	geom_line(size=2) +
	geom_point(size=4) +
	scale_colour_viridis_d(begin=0, end=0.5) + 
	ggtitle("Accident year parameters")



```

* This looks very good - the fit is very similar, but we have 7 fewer parameters.


### Development year

* Now we do the same thing for development year
* Note that the glm_fit2 model (and the glm_fit1 model too) do not have a parameter for dev_year = 1 as this is the base level.
    + This means that the parameter is really 0, so we must remember to include this.


```{r}

# extract the data
dt_dev_year <- glm_fit2$coeff_table[grepl("dev_year", parameter),  
									][, dev_year := as.integer(gsub("dev_year_factor", "", parameter))][]   # known data.table printing bug

# add year 1
dt_dev_year <- rbind(dt_dev_year, data.table(parameter="dev_year_factor1", coeff_glm_fit2=0, dev_year=1))
setorder(dt_dev_year, dev_year)


# plot
ggplot(data=dt_dev_year, aes(x=dev_year, y=coeff_glm_fit2)) +
	geom_line(size=2, colour="#440154ff") +
	geom_point(size=4, colour="#440154ff") +
	ggtitle("Development year parameters")


```


* Looking at this plot, it appears that a straight line would fit quite well
* This fit would be improved by allowing the straight line to bend (have a knot) at dev_year = 7
    + So let's try this below
    + note we actually fit dev_year - 1 rather than dev_year
        - this means that the parameter estimate at dev_year = 1 is 0, just as it is in the glm_fit2 model, so it makes the results comparable
        - if we fit dev_year, then the parameter estimate at dev_year=1 would be non-zero, so the two fits would be shifted relative to each other and we would need to adjust for that.


```{r}

# add dev-1 and dev-7 terms
msdata[, dev_year_m1 := dev_year - 1]
msdata[, dev_year_ge_7 := pmax(dev_year-7.5, 0)]

# fit the model
glm_fit3 <- glm(data = msdata, 
	family = quasipoisson(link = "log"),
	formula = "incremental ~ acc_year + acc_year_2 + dev_year_m1 + dev_year_ge_7")

# extract and save the coefficient table
glm_fit3$coeff_table <- data.table(parameter = names(glm_fit3$coefficients), coeff_glm_fit3 = glm_fit3$coefficients)

# display a summary of the model
summary(glm_fit3)


```


* Assuming the fit is satisfactory, our original model with 19 parmaeters has now been simplified to 5 parameters - much more parsimonious and robust.
* Let's check the fit by dev_year to see.


```{r}

# get the dev_year fit under the new model and add to the data.table containing the factor level parameters
p1 <- glm_fit3$coeff_table[parameter == "dev_year_m1", coeff_glm_fit3]
p2 <- glm_fit3$coeff_table[parameter == "dev_year_ge_7", coeff_glm_fit3]
dt_dev_year[, coeff_glm_fit3 := p1*(dev_year-1) + p2*pmax(0, dev_year-7.5) ]


# make long for ggplot
dt_dev_year_plot <- melt(dt_dev_year, id.vars = "dev_year", measure.vars = c("coeff_glm_fit2", "coeff_glm_fit3"), variable.name="model", value = "estimate")

# remove the coeff_ from the model names
dt_dev_year_plot[, model := gsub("coeff_", "", model, fixed=TRUE)]


ggplot(data=dt_dev_year_plot, aes(x=dev_year, y=estimate, colour=model)) +
	geom_line(size=2) +
	geom_point(size=4) +
	scale_colour_viridis_d(begin=0, end=0.5) +
	ggtitle("Accident year parameters")



```


* This looks good.
* However dev_year = 2 is a bit underfit in the latest model, so we can add something to improve this fit
* So refit and replot.



```{r}

msdata[, dev_year_eq_2 := as.integer(dev_year == 2)]

glm_fit4 <- glm(data = msdata, 
	family = quasipoisson(link = "log"),
	formula = "incremental ~ acc_year + acc_year_2 + dev_year_m1 + dev_year_ge_7 + dev_year_eq_2")


glm_fit4$coeff_table <- data.table(parameter = names(glm_fit4$coefficients), coeff_glm_fit4 = glm_fit4$coefficients)


p1 <- glm_fit4$coeff_table[parameter == "dev_year_m1", coeff_glm_fit4]
p2 <- glm_fit4$coeff_table[parameter == "dev_year_ge_7", coeff_glm_fit4]
p3 <- glm_fit4$coeff_table[parameter == "dev_year_eq_2", coeff_glm_fit4]
dt_dev_year[, coeff_glm_fit4 := p1*(dev_year-1) + p2*pmax(0, dev_year-7.5) + p3*(dev_year == 2) ]


# make long for ggplot
dt_dev_year_plot <- melt(dt_dev_year, id.vars = "dev_year", measure.vars = c("coeff_glm_fit2", "coeff_glm_fit4"), variable.name="model", value = "estimate")

# remove the coeff_ from the model names
dt_dev_year_plot[, model := gsub("coeff_", "", model, fixed=TRUE)]


ggplot(data=dt_dev_year_plot, aes(x=dev_year, y=estimate, colour=model)) +
	geom_line(size=2) +
	geom_point(size=4) +
	scale_colour_viridis_d(begin=0, end=0.5) +
	ggtitle("Accident year parameters")



```

* Looks good!


### Identifying missing structure

* The second part of the model refining process involves checking for missing structure.
* Let's have a better look at the heat map.


```{r}
msdata[, residuals4 := rstandard(glm_fit4)
	   ][, fitted4 := glm_fit4$fitted.values
	     ][, linear_predictor4 := log(fitted4)
	       ][, AvsF4 := incremental / fitted4
	         ][, AvsF_restricted4 := log(pmax(0.5, pmin(2,AvsF4)))]


p_hm <- ggplot(data=msdata, aes(x=dev_year, y=acc_year)) + 
	geom_tile(aes(fill = AvsF_restricted4))+scale_y_reverse()+
	scale_fill_gradient2(name="AvF_min", low="royalblue", mid="white", high="red", midpoint=0, space="Lab", na.value="grey50", guide="colourbar")+
	labs(x="Development year", y="Accident year")+
	theme(legend.position = "none")+
	theme(axis.title.x = element_text(size=8), axis.text.x  = element_text(size=7))+
	theme(axis.title.y = element_text(size=8), axis.text.y  = element_text(size=7))+
	theme(panel.background = element_rect(fill = "grey", colour = "grey", size = 2, linetype = "solid"),
		  panel.grid = element_line(colour="grey")) + 
	NULL

print(p_hm)


```



Let's look at the heatmap again, with some annotations


```{r}

p_hm + 
	annotate(geom="rect", xmin= 0.5, xmax=1.5, ymin=0.5, ymax=6.5, colour="darkblue", alpha=0.1, size=1.5) +
	annotate(geom="rect", xmin= 0.5, xmax=1.5, ymin=6.5, ymax=10.5, colour="darkred", alpha=0.1, size=1.5) +
	annotate(geom="rect", xmin= 1.5, xmax=2.5, ymin=0.5, ymax=6.5, colour="darkred", alpha=0.1, size=1.5) +
	annotate(geom="rect", xmin= 1.5, xmax=2.5, ymin=6.5, ymax=9.5, colour="darkblue", alpha=0.1, size=1.5) +
	annotate(geom="segment", x=3, xend=3, y=1, yend=8, arrow=arrow(), colour="darkblue", size=2) +
	annotate(geom="rect", xmin= 3.5, xmax=4.5, ymin=0.5, ymax=7.5, colour="darkred", alpha=0.1, size=1.5) 



```


We see:

* development year 1, a distinct area of blue in the earlier accident years (A < F), followed by red (A > F)
* development year 2, a distinct area of red in the earlier accident years (A > F), followed by blue (A < F)
* development year 3, a possible progression from red to blue with increasing accident year (F increasing relative to A)
* development year 4, nearly all red (A > F)

This suggests the payment pattern has altered and can be accommodated by (mostly) interaction terms within the GLM. Consider adding the following terms:

* (development year = 1) * (accident year is between 1 and 6)
* (development year = 2) * (accident year is between 1 and 6)
* (development year = 3) * (accident year linear trend)
* (development year = 4) 


So, let's refit the model with terms to capture these and have a look at the heat map again


```{r}

# add the new terms
msdata[, dev_year_eq_1 := as.integer(dev_year == 1)]
msdata[, dev_year_eq_3 := as.integer(dev_year == 3)]
msdata[, dev_year_eq_4 := as.integer(dev_year == 4)]
msdata[, acc_year_1_6 := as.integer(acc_year >= 1 & acc_year <= 6)]


glm_fit5 <- glm(data = msdata, 
	family = quasipoisson(link = "log"),
	formula = "incremental ~ acc_year + acc_year_2 + dev_year_m1 + dev_year_ge_7 + dev_year_eq_2 + dev_year_eq_4 +
	dev_year_eq_1:acc_year_1_6 +  dev_year_eq_2:acc_year_1_6 + dev_year_eq_3:acc_year ")



glm_fit5$coeff_table <- data.table(parameter = names(glm_fit5$coefficients), coeff_glm_fit5 = glm_fit5$coefficients)

# print the coefficient table

glm_fit5$coeff_table |> 
	knitr::kable(digits=c(0, 4))

```


* This model should match that displayed in Table 7-5 of the monograph - and indeed it does (some very minor differences in parameter values - the model in the monograph was fitted in SAS).

* Look at the heat map again with annotations - has the model resolved the identified issues?


```{r}
# attach fitteds and residuals
msdata[, residuals5 := rstandard(glm_fit5)
	   ][, fitted5 := glm_fit5$fitted.values
	     ][, linear_predictor5 := log(fitted5)
	       ][, AvsF5 := incremental / fitted5
	         ][, AvsF_restricted5 := log(pmax(0.5, pmin(2,AvsF5)))]



p_hm <- ggplot(data=msdata, aes(x=dev_year, y=acc_year)) + 
	geom_tile(aes(fill = AvsF_restricted5))+scale_y_reverse()+
	scale_fill_gradient2(name="AvF_min", low="royalblue", mid="white", high="red", midpoint=0, space="Lab", na.value="grey50", guide="colourbar")+
	labs(x="Development year", y="Accident year")+
	theme(legend.position = "none")+
	theme(axis.title.x = element_text(size=8), axis.text.x  = element_text(size=7))+
	theme(axis.title.y = element_text(size=8), axis.text.y  = element_text(size=7))+
	theme(panel.background = element_rect(fill = "grey", colour = "grey", size = 2, linetype = "solid"),
		  panel.grid = element_line(colour="grey")) + 
	annotate(geom="rect", xmin= 0.5, xmax=1.5, ymin=0.5, ymax=6.5, colour="darkblue", alpha=0.1, size=1.5) +
	annotate(geom="rect", xmin= 0.5, xmax=1.5, ymin=6.5, ymax=10.5, colour="darkred", alpha=0.1, size=1.5) +
	annotate(geom="rect", xmin= 1.5, xmax=2.5, ymin=0.5, ymax=6.5, colour="darkred", alpha=0.1, size=1.5) +
	annotate(geom="rect", xmin= 1.5, xmax=2.5, ymin=6.5, ymax=9.5, colour="darkblue", alpha=0.1, size=1.5) +
	annotate(geom="segment", x=3, xend=3, y=1, yend=8, arrow=arrow(), colour="darkblue", size=2) +
	annotate(geom="rect", xmin= 3.5, xmax=4.5, ymin=0.5, ymax=7.5, colour="darkred", alpha=0.1, size=1.5) 


print(p_hm)



```


* This looks much better.

* We should also look at the residual plots again


```{r, fig.height=6, fig.width=9}

p1 <- ggplot(data=msdata, aes(x=linear_predictor5, y=residuals5, colour=dev_year)) +
	geom_point(size=2) +
	scale_colour_viridis(begin=0.9, end=0) +
	theme_bw() + 
	theme(legend.position = "none") +
	ggtitle("Linear predictor")


p2 <- ggplot(data=msdata, aes(x=acc_year, y=residuals5)) +
	geom_point(size=2, colour="#2d708eff") +
	theme_bw() + 
	ggtitle("Accident year")

p3 <- ggplot(data=msdata, aes(x=dev_year, y=residuals5)) +
	geom_point(size=2, colour="#2d708eff") +
	theme_bw() + 
	ggtitle("Development year")

p4 <- ggplot(data=msdata, aes(x=cal_year, y=residuals5)) +
	geom_point(size=2, colour="#2d708eff") +
	theme_bw() + 
	ggtitle("Calendar year")

p <- (p1 + p2) / (p3 + p4)

p

```



## Loss reserve

* Now that we have a model, let's produce the estimate of the outstanding claims by accident year and in total.
    + Take the lower triangle data [futdata] created above
    + Add on the new variates we created
    + Score the model on this data
    + Summarise the results

<br>
Create the data and score using `predict`


```{r}

# add all model variates
futdata[, acc_year_2 := acc_year^2
		][, dev_year_m1 := dev_year - 1
		  ][, dev_year_ge_7 := pmax(0, dev_year - 7.5)
		    ][, dev_year_eq_1 := as.integer(dev_year == 1)
		      ][, dev_year_eq_2 := as.integer(dev_year == 2)
		        ][, dev_year_eq_3 := as.integer(dev_year == 3)
		          ][, dev_year_eq_4 := as.integer(dev_year == 4)
		            ][, acc_year_1_6 := as.integer(acc_year>=1 & acc_year <=6)]


x <- predict(glm_fit5, newdata = futdata, type="response")
futdata[, incremental := x]

head(futdata)


```


Get reserves by accident year and in total 

```{r}

ocl_year <- futdata[,  lapply(.SD, sum), .SDcols=c("incremental"), by="acc_year"]
ocl_total <- ocl_year[, sum(incremental)]


ocl_year |> 
	knitr::kable(digits=c(0, 0))

```


The total reserve is

```{r}

ocl_total |> round(0)

```


* These results are similar, though not identical, to the results given in Table 7-6 of the monograph.
* This is because the *forecast* column of the monograph contains bootstrapped means rather than the model mean.


## Conclusion

The aim of this article has been to demonstrate fitting a GLM to a loss reserve following the example used in the CAS monograph.
We started with the chain ladder equivalent - the cross classified model with an over-dispersed Poisson distribution, then first simplified it and second, extended it to include some interactions.
We also cover how to create some of the plots discussed in the monograph in R, in particular residual scatter plots and the heat maps.



## Session information

To assist with reproducibility, here are details of my R session.


```{r}
sessionInfo()  
```




