{"word_count": 691, "character_count": 6363}
A few month ago, I did mention a graph, of some so-called Lorenz curves to compare regression models, see e.g. Progressive's slides (thanks Guillaume for the reference)The idea is simple. Consider some model for the pure premium (in insurance, it is the quantity that we like to model), i.e. the conditional expected valeurOn some dataset, we have our predictions, as well as observed quantities, . The curve are obtained simply :sort the observations so thatbased on that ordering (from high risks to low risks, based on our predictions), we plot Lorenz curveThat looks nice, and as mentioned in my previous post, I have two issues regarding that graph. The first one is the 'upper-bound' : I do not understand this "perfect model". The second one is the (pseudo) 'lower-bound', which is the "average model". At least, for the later, I understand how to get it, but I do not understand what it means if we are (sometimes) below that line.In order to understand a bit more, let us consider the following example.we have some covariates given those covariates, we have parameters , with some random effect given those parameters, we have observations Let us generate such a dataset,n = 10s = 1tau = 1X = rnorm(n)Theta = -1+X+rnorm(n)*tauY = exp(Theta+rnorm(n)*s)base = data.frame(X,Theta,S,m)and compare different strategies to get Lorenz curves :the random case : observations are ordered randomly. That simply means that predictions are identical   Ys = rev(base[sample(1:n),"Y"])   L  = c(0,cumsum(Ys))/sum(Ys)the partial-information model : we fit a model using our covariates, and we sort observations base on that model   reg = lm(log(Y)~X,data=base)   base$m = predict(reg)   Ys = rev(base[order(base$m),"Y"])   L  = c(0,cumsum(Ys))/sum(Ys)a second partial-information model : assume that we know that the model should be a function increasing with the covariate. In that case, we can directly use it to sort our sample   Ys = rev(base[order(base$X),"Y"])   L  = c(0,cumsum(Ys))/sum(Ys)the perfect model : here we observe . That's clearly, here, the best we can do (we cannot observe the noise)   Ys = rev(base[order(base$Theta),"Y"])   L  = c(0,cumsum(Ys))/sum(Ys)The empirical curve will be function of the number of observvations, , the variance of the random effect (that cannot be observed) , and the variance of the noise . In order to visualize what's going on, let us generate several samples (say 1,000), and plot some confidence intervals (90% and 50%), as well as the average curve, obtained on those simulated samples. For instance, in the random case, useL1=function(n=10,s=1,ns=1e3,tau=.2){  CY=matrix(NA,ns,n+1)  for(i in 1:ns){    X=rnorm(n)    Theta=-1+X+rnorm(n)/5    S=exp(Theta+rnorm(n)*s)    m=exp(Theta+.5*var(log(S)))    base=data.frame(X,Theta,S,m)    Y=rev(base[sample(1:n),"S"])    CY[i,]=c(0,cumsum(Y))/sum(Y)  }  Qinf1=apply(CY,2,function(x) quantile(x,.05))  Qsup1=apply(CY,2,function(x) quantile(x,.95))  Qinf2=apply(CY,2,function(x) quantile(x,.25))  Qsup2=apply(CY,2,function(x) quantile(x,.75))  Qm=apply(CY,2,mean)  return(data.frame(inf10=Qinf1,inf25=Qinf2,moy=Qm,sup75=Qsup2,sup90=Qsup1))}Then useviz=function(B){  n=nrow(B)-1  u=(0:n)/n  plot(0:1,0:1,col="white",xlab="",ylab="")  polygon(c(u,rev(u)),c(B$inf10,rev(B$sup90)),col=rgb(0,0,1,.2),border=NA)  polygon(c(u,rev(u)),c(B$inf25,rev(B$sup75)),col=rgb(0,0,1,.4),border=NA)  lines(u,B$moy,col="blue",lwd=2)  segments(0,0,1,1)}e.g.viz(L1(s=1,n=100,tau=1))Which is what we were expecting: we are around the diagonal, and on average, we have the diagonal.If we consider now the partial information model,L2=function(n=10,s=1,ns=1e3,tau=.2){  CY=matrix(NA,ns,n+1)  for(i in 1:ns){    X=rnorm(n)    Theta=-1+X+rnorm(n)*tau    S=exp(Theta+rnorm(n)*s)    m=exp(Theta+.5*var(log(S)))    base=data.frame(X,Theta,S,m)    Y=rev(base[order(base$X),"S"])    CY[i,]=c(0,cumsum(Y))/sum(Y)  }  Qinf1=apply(CY,2,function(x) quantile(x,.05))  Qsup1=apply(CY,2,function(x) quantile(x,.95))  Qinf2=apply(CY,2,function(x) quantile(x,.25))  Qsup2=apply(CY,2,function(x) quantile(x,.75))  Qm=apply(CY,2,mean)  return(data.frame(inf10=Qinf1,inf25=Qinf2,moy=Qm,sup75=Qsup2,sup90=Qsup1))} viz(L2(s=1,n=100,tau=1))We have here a curve above the first diagonal, as on the Figure mentioned in the introduction of this post. Note that if we use directly the covariate to sort, instead of  model, we have a rather similar graph,L2b=function(n=10,s=1,ns=1e3,tau=.2){  CY=matrix(NA,ns,n+1)  for(i in 1:ns){    X=rnorm(n)    Theta=-1+X+rnorm(n)*tau    S=exp(Theta+rnorm(n)*s)    m=exp(Theta+.5*var(log(S)))    base=data.frame(X,Theta,S,m)    reg=lm(log(S)~X,data=base)    base$m=predict(reg)    Y=rev(base[order(base$m),"S"])    CY[i,]=c(0,cumsum(Y))/sum(Y)  }  Qinf1=apply(CY,2,function(x) quantile(x,.05))  Qsup1=apply(CY,2,function(x) quantile(x,.95))  Qinf2=apply(CY,2,function(x) quantile(x,.25))  Qsup2=apply(CY,2,function(x) quantile(x,.75))  Qm=apply(CY,2,mean)  return(data.frame(inf10=Qinf1,inf25=Qinf2,moy=Qm,sup75=Qsup2,sup90=Qsup1))} viz(L2b(s=1,n=100,tau=1))And finally, our perfect model would beL3=function(n=10,s=1,ns=1e3,tau=.2){  CY=matrix(NA,ns,n+1)  for(i in 1:ns){    X=rnorm(n)    Theta=-1+X+rnorm(n)*tau    S=exp(Theta+rnorm(n)*s)    m=exp(Theta+.5*var(log(S)))    base=data.frame(X,Theta,S,m)    Y=rev(base[order(base$m),"S"])    CY[i,]=c(0,cumsum(Y))/sum(Y)  }  Qinf1=apply(CY,2,function(x) quantile(x,.05))  Qsup1=apply(CY,2,function(x) quantile(x,.95))  Qinf2=apply(CY,2,function(x) quantile(x,.25))  Qsup2=apply(CY,2,function(x) quantile(x,.75))  Qm=apply(CY,2,mean)  return(data.frame(inf10=Qinf1,inf25=Qinf2,moy=Qm,sup75=Qsup2,sup90=Qsup1))} viz(L3(s=1,n=100,tau=1))The curve is here in the upper left corner, again, as stated in the introductionary Figure.But if we start playing with those functions, we can get weird results, e.g.viz(L3(s=1,n=100,tau=1))Here, we have a large variability for the noise . This might happen in insurance when dealing with liabilty claims, with large variance, and heavy tails. In that case, observe that the perfect model can even be below the diagonal. Which makes me more and more suspicious regarding the use of such graphical tools in insurance, to assess if a model a good, or not. Here, the perfect model might seem worst than a random one (where an identical value for 's is considered, even if the true values of 's are given).