Tag Archives: GLM

Talk at Desjardins General Insurance

This afternoon, I will give a talk at the seminar of the R&D department at Desjardins General Insurance, on correlation in claims reserving. A lot of interesting papers have been published recently on that topic. On multivariate Chain Ladder, some interesting articles have been published, e.g. the one by Carsten Prohl and Klaus Schmidt (here) or the one by Michael Merz and Mario Wuthrich (there).

But I think another interesting perspective (so far, not in claims reserving, but one should find some time to look at it) should be about multivariate regression (multivariate GLM’s), e.g.

All that will be mentioned in the talk. Slides can be downloaded here,

The dataset used in the example can be obtained with the code below

> P.corp=read.table("http://freakonometrics.blog.free.fr/public/data/auto-corporel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.corp=as.matrix(P.corp) > n=nrow(P.corp) > P.mat =read.table("http://freakonometrics.blog.free.fr/public/data/auto-materiel.csv", +        header=FALSE,sep=";",na.strings = "NA",dec=",") > P.mat=as.matrix(P.mat) > P.mat=P.mat[1:n,1:n] >  P.mat = P.mat[2:10,1:9] >  P.corp= P.corp[2:10,1:9] > n=9 > P.tot = P.mat + P.corp

Too large datasets for regression ? What about subsampling….

recently, a classmate working in an insurance company told me he had too large datasets to run simple regressions (GLM, which involves optimization issues), and that they were thinking of a reward for the one who will write the best R-code (at least the fastest). My first idea was to use subsampling techniques, saying that 10 regressions on 100,000 observations can take less time than a regression on 1,000,000 observations. And perhaps provide also better results…

  • Time to run a regression, as a function of the number of observations

Here, I generate a dataset as follows


and we fit


where http://freakonometrics.hypotheses.org/files/2016/11/largesamp03.png is a spline function (just to make it as general as possible, since in insurance ratemaking, we include continuous variates that do not influence claims frequency linearly in the score). Yes, there might be also useless variables, including one of them which is strongly correlated with one that has an impact in the regression. The code to generate the dataset is simply

> n=10000
> X1=rexp(n)
> X2=sample(c("A","B","C"),size=n,replace=TRUE)
> X3=runif(n)
> Z=rmnorm(n,c(0,0),matrix(c(1,0.8,.8,1),2,2))
> X4=Z[,1]
> X5=Z[,2]
> X6=X1^2
> E=runif(n)
> lambda=.2*X5-4*dbeta(X3,2,5)+X1+
> Y=rpois(n,exp(lambda))
> base=data.frame(Y,X1,X2,X3,X4,X5,X6,E)

We would like the study the time it takes to run a regression, as a function of the size (i.e. the number of lines http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png) of the dataset.

> system.time( glm(Y~bs(X1)+X2+X3+X4+
+ X5+X6+offset(log(E)),family=poisson,
+ data=base) )
utilisateur     système      écoulé
0.25        0.00        0.25

Here, the time I look at is the last one. But so far, it was rather simple, but it is not the best model I can get. Let us use a stepwise (backward) variable selection,

> system.time( step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=poisson,
+ data=base)) )
Start:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
Step:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + offset(log(E))
Df Deviance    AIC
<none>        2236.0 2882.1
- X5      1   2240.1 2884.2
- X4      1   2244.1 2888.2
- X3      1   4783.2 5427.3
- X2      2   5311.4 5953.5
- bs(X1)  3   6273.7 6913.8
utilisateur     système      écoulé
1.82        0.03        1.86

Finally, from the first regression, we have points in black (based on 200 simulated datasets), and with a stepwise procedure, we have the points in red.

i.e. it might look linear (proportional), but if it was linear, then on a log-log scale, we should have also straigh lines, with slope 1,

Actually, it looks like a convex function.

The interpretation of that convexity might lead to misinterpretation. On the graph below on the left, on a dataset two times bigger than the previous one (black point) will be less than two times longer to run, while on the right, it will be more than two timess longer,

Convexity can simply be interpreted as “too large datasets take time, and too small too…”. Which is a first step: it should be interesting, in some cases, to run several regressions on smaller datasets….

  • Running 100 regressions on 100 lines, or running 1 regression on 10,000 lines ?

Here, we have datasets with http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=200,000 lines. The questions is how long will it take if we subdived into http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png subsamples (of equal size), and run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regressions ?

> nk=trunc(n/k)rep(1:k,each=nk); nt=nk*k
> base=data.frame(Y[1:nt],X1[1:nt],
+ X2[1:nt],X3[1:nt],X4[1:nt],X5[1:nt],
+ X6[1:nt],E[1:nt],classe)
> system.time( for(j in 1:k){
+  glm(Y~bs(X1)+X2+X3+X4+X5+
+ X6+offset(log(E)),family=poisson
+ ,data=base,subset=classe==j) })
utilisateur     système      écoulé
1.31        0.00        1.31
> system.time( for(j in 1:k){
+      step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=
+ poisson,data=base,subset=classe==j)) })
Start:  AIC=183.97
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))


  Df Deviance    AIC
<none>        117.15 213.04
- X2      2   250.15 342.04
- X3      1   251.00 344.89
- X4      1   420.63 514.53
- bs(X1)  3   626.84 716.74
utilisateur     système      écoulé
11.97        0.03       12.31

On the graph below, we have the time (y-axis, here on a log scale) it took to run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regression on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png, as function of http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png (x-axis), including the time it took to run the regression on a dataset of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png which is the concentration of dots on the left (i.e. http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png=1), both on the 6 regressors – in black – and with a strepwise procedure – in red. One has to keep in mind that I did not remove the printing option in the stepwise procedure, so it might be difficult to compare the two clouds (black vs. red). Nevertheless, we clearly see that if we run http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png regression on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png, when http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png is not too large, i.e. less than 10 or 15, it is not longer than the regression on http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=200,000 lines.

So here we see that running 100 regressions on 2,000 lines is longer than running 1 regression on 200,000 lines… But maybe we are not comparing things that are actually comparable: what if it takes a bit longer, but we strongely improve the quality of our estimators ?

  • What about the quality of the output ?

Here, we consider only one dataset, with http://freakonometrics.hypotheses.org/files/2016/11/largesamp04.png=100,000 lines (just to make it run a bit faster). And http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png=20 subsets. Recall that the generated dataset is from


and we fit


Here, we plot here http://freakonometrics.hypotheses.org/files/2016/11/largesamp07.png and a confidence interval, defined as


The lightblue segment is the initial estimator, while the blue one is obtained from the stepwise procedure. The grey area represent the estimation on the overall sample, while the http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png segments on the right are the http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png estimators (each on samples of size http://freakonometrics.hypotheses.org/files/2016/11/largesamp06.png).

We can see that we have much more volatility on those http://freakonometrics.hypotheses.org/files/2016/11/largesamp05.png estimators, but the average (horizontal doted lines) are not so bad… The true value (i.e. the one used to generate the dataset is the dotter black horizontal line).
And if we repeat that on 1,000 simulated dataset, we obtaind the following distribution for http://freakonometrics.hypotheses.org/files/2016/11/largesamp07.png (blue line), so we have an unbiased estimator of our parameter (the verticular line being here the true value), here including a stepwise procedure,

But if we add the the red curve is the average of the http://freakonometrics.hypotheses.org/files/2016/11/largesamp09.png the previous one being now the clear blue line in the back, we see that taking average of estimators on subsamples is not bad at all, on the contrary,

and for those who think that the stepwise procedure is a mistake, here is what we get without it,

So what we can see is that running 20 regressions can take (a little) more time (from what we’ve seen earlier) than running only one on the whole dataset…. but it provides better estimates. So the tradeoff is not that simple, and maybe running several regressions on huge datasets can be a proper alternative.

La tarification avec SAS

En tarification, il est possible d’utiliser d’autres logiciels que R, en particulier, il semble que l’on puisse faire deux ou trois choses avec SAS…. J’en parle un peu car il semble  que, paradoxalement, les asssureurs préfèrent encore SAS à R (par exemple). Et comme plusieurs étudiants m’avaient demandé “et comment on fait avec SAS ?“. Bon, par contre je ne mets que les choses de base, parce que SAS est assez limité sur ce qu’il peut faire….

Pour suivre un peu le plan du cours, la première étape est de définir une variable d’exposition dans la table,

DATA contrats;
SET lib.contrats;
lnexpo = log(expo);

Pour faire une régression de Poisson, ce n’est pas forcément compliqué,

ODS OUTPUT ParameterEstimates=Genmod1_Param
MODEL nbsin = ageconducteur /
                  dist = poisson   
                  link = log   
                  offset = lnexpo 

La sortie SAS a alors l’allure suivante

                                  The GENMOD Procedure
                    Critère pour évaluer la qualité de l'ajustement
              Critère                   DF          Valeur       Valeur/DF
              Deviance                63E3      26872.5334          0.4237
              Scaled Deviance         63E3      26872.5334          0.4237
              Pearson Chi-Square      63E3      73275.5362          1.1553
              Scaled Pearson X2       63E3      73275.5362          1.1553
              Log Likelihood                   -18474.2667

       Algorithm converged.
                       Analyse des résultats estimés de paramètres

                              Erreur      Wald 95Limites
Paramètre    DF   Estimation   standard      de confiance %    Khi 2   Pr > Khi 2
Intercept     1      -3.5164     0.0851    -3.6832  -3.3496   1708.02       <.0001
ageconducteur 1       0.0168     0.0014     0.0141   0.0195    146.73       <.0001
Scale         0       1.0000     0.0000     1.0000   1.0000
NOTE: The scale parameter was held fixed.

                         Statistiques LR pour Analyse de Type 3
                      Source           DF      Khi 2    Pr > Khi 2
                      ageconducteur     1     148.72        <.0001

Il est aussi possible de faire des GAM (i.e. du lissage de la variable explicative – continue – avec des fonctions splines)

MODEL nbsin = spline(ageconducteur) / dist = Poisson;

et on peut faire des prédictions avec ce modèle (la sortie n’apporte pas grand chose, en pratique),

DATA gam;
SET gam;
pred_nbsin_gam = exp(P_nbsin);
KEEP ageconducteur pred_nbsin_gam;

Enfin, on peut tenter de faire un joli graphique. Pour cela, on calcule les prédictions de trois modèles, le premier étant des nombres moyens de sinistres par âge

PROC SORT DATA = base; BY ageconducteur; RUN; QUIT;
BY ageconducteur;
VAR nbsin;
WEIGHT expo;
OUTPUT OUT = nbsin_age (DROP = _TYPE_ _FREQ_) MEAN=mo

ensuite, on fait un modèle GLM, et  un modèle GAM, et on récupère les sorties

PROC SORT DATA = nbsin_age; BY age_cond; RUN; QUIT;
PROC SORT DATA = gam; BY age_cond; RUN; QUIT;
DATA nbsin_age;
MERGE nbsin_age
BY age_cond;

On essaye de faire le dessin (je passe les lignes de commande, il y en a une vingtaine)

Pour faire une régression quasiPoisson, le code a l’allure suivante,s

ODS OUTPUT ParameterEstimates=Genmod1bis_Param
MODEL nbsin = ageconducteur /
                 dist = poisson  
                 link = log      
                 offset = lnexpo 
                 scale = deviance;

La sortie donne alors l’estimation du paramètre de surdispersion (ou sur cet exemple de sousdispersion)

                      Analyse des résultats estimés de paramètres

                                   Erreur    Wald 95Limites
Paramètre      DF   Estimation   standard    de confiance %     Khi 2   Pr > Khi 2

Intercept        1     -3.5164     0.0554  -3.6249  -3.4078   4031.29       <.0001
ageconducteur    1      0.0168     0.0009   0.0150   0.0186    346.32       <.0001
Scale            0      0.6509     0.0000   0.6509   0.6509

On notera que pour calculer le critère d’Akaike, ça n’est pas forcément trivial,

%MACRO CALCUL_AIC_BIC(infomodel=, param=);
    DATA _null_;
    SET &infomodel.;
    IF Criterion = "Log Likelihood" THEN CALL SYMPUT("Loglike", Value);
    IF Criterion = "Deviance" THEN CALL SYMPUT("n_etoile", Df);
    DATA _null_;
    SET &param.  end=fin;
    RETAIN nb_df 0;
    nb_df = nb_df + df;
    IF fin THEN CALL SYMPUT("k", nb_df);
    DATA genmod_aic_bic;
    SET &param.;
    FORMAT Loglike 12.2 K 10. N 10. AIC_CALC 12.2 BIC_CALC 12.2;
    Loglike = 0; K = 0; N = 0; AIC_CALC = 0; BIC_CALC = 0;
    IF Parameter = "Intercept";
    DATA genmod_aic_bic;
    SET genmod_aic_bic;
    Loglike = &loglike.;
    K = &k.;
    N = %eval(&n_etoile. + &k.);
    AIC_CALC = 2 * Loglike + 2 * K;
    BIC_CALC = 2 * Loglike + K * log(N);
    PROC PRINT DATA = genmod_aic_bic;
    RUN; QUIT;
%CALCUL_AIC_BIC(infomodel=Genmod2_InfoModele, param=Genmod2_Param);