# Tag Archives: bootstrap

On Thursday, March 9nd, I will give the second lecture of the PhD course on advanced tools for econometrics, on simulation techniques (and bootstrap). Slides are available online.

The first part is this Thurdays, on Nonlinearities in Econometric models.

I will give a short graduate course for PhD students, in Rennes, on Thurday mornings, in March (2nd, 9th, 23rd and 30th). The agenda will be

1. Nonlinear Regression Models and Smoothing Techniques

2. Bootstrapping and Regression

3. Penalized Regression Models and LASSO

4. Quantile Regression and Expectiles

There will be slides available by the end of February.

# Confidence Regions for Parameters in the Simplex

Consider here the case where, in some parametric inference problem, parameter $\boldsymbol{\theta}$ is a point in the Simplex,

$\mathcal{S} =\{\boldsymbol{\theta}=(\theta_1 , \cdots,\theta_k) | \theta_i \ge 0, 1 \le i \le k, \sum_{i=1}^{k} \theta_i=1\}$For instance, consider some regression, on compositional data,

> library(compositions)
>  data(DiagnosticProb)
>  Y=DiagnosticProb[,"type"]-1
>  X=DiagnosticProb[,c("A","B","C")]
>  model = glm(Y~ilr(X),family=binomial)
>  b = ilrInv(coef(model)[-1],orig=X)
>  as.numeric(b)
[1] 0.3447106 0.2374977 0.4177917

We can visualize that estimator on the simplex, using

>  tripoint=function(s){
+    p=s/sum(s)
+    abc2xy(matrix(p,1,3))
+  }

>  lab=LETTERS[1:3]
>  xl=c(-.1,1.25)
>  yl=c(-.1,1.15)
>  library(trifield)
>  A=abc2xy(matrix(c(1,0,0),1,3))
>  B=abc2xy(matrix(c(0,1,0),1,3))
>  C=abc2xy(matrix(c(0,0,1),1,3))
>  plot(0:1,0:1,col="white",
+  xlim=xl,ylim=yl,xlab="",ylab="",axes=FALSE)
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

If we want to compute a ‘confidence region’, we can either use Bayesian models (with a Dirichlet distribution as prior distribution), or use bootstrap. We will use here the second idea

>  MB=matrix(NA,1e4,2)
>  for(sim in 1:1e4){
+    idx=sample(1:nrow(DiagnosticProb),
+    size=nrow(DiagnosticProb),replace=TRUE)
+  Y=DiagnosticProb[idx,"type"]-1
+  X=DiagnosticProb[idx,c("A","B","C")]
+  model = glm(Y~ilr(X),family=binomial)
+  MB[sim,]=tripoint(as.numeric(
+    ilrInv(coef(model)[-1],orig=X)))}

To get some ‘confidence region’, we can then use the bagplot, to get either a region where 50% of the boostraped estimators are, or 95%,

>  library(aplpack)
> P1=bagplot(MB[,1],MB[,2], factor =1.96, cex=.9,
+ dkmethod=2,show.baghull=TRUE)
> P2=bagplot(MB[,1],MB[,2], factor =0.67, cex=.9,
+ dkmethod=2,show.baghull=TRUE) 

Then we can easily plot those two regions,

>  plot(0:1,0:1,col="white")
>  polygon(rbind(A,B,C),col="light yellow")
>  text(B[1],-.05,lab[2])
>  text(A[1],1.05,lab[1])
>  text(C[1],-.05,lab[3])
>  polygon(P1$hull.loop,col="yellow",border=NA) > polygon(P2$hull.loop,col="orange",border=NA)
>  segments((A[1]+C[1])/2,(A[2]+C[2])/2,B[1],B[2],col="grey",lty=2)
>  segments((A[1]+B[1])/2,(A[2]+B[2])/2,C[1],C[2],col="grey",lty=2)
>  segments((B[1]+C[1])/2,(B[2]+C[2])/2,A[1],A[2],col="grey",lty=2)
>  points(tripoint(b),pch=19,cex=2,col="red")

# Some thoughts on Economics, Mathematics, Econometrics, Statistics, Machine Learning, etc

There were a lot of posts, recently, related to those topics, starting with Noah Smith ‘s piece entitled “Economics has a Math Problem” and more recently “Econometrics, Math, and Machine Learning…what?” by Matt Bogard. I don’t have (yet) a clear mind on those issues, but there are still a few thoughts that I wanted to share. I did not really want to, but I’ve been asked, on Twitter, and I thought it might be good to write them down, to clarify some ideas I have, but also (probably, hopefully) to get interesting feedbacks.

# An Attempt to Understand Boosting Algorithm(s)

Last tuesday, at the annual meeting of the French Economic Association, I was having lunch with Alfred, and while we were chatting about modeling issues (econometric models against machine learning prediction), he asked me what boosting was. Since I could not be very specific, we’ve been looking at wikipedia webpage.

Boosting is a machine learning ensemble meta-algorithm for reducing bias primarily and also variance in supervised learning, and a family of machine learning algorithms which convert weak learners to strong ones

One should admit that it is not very informative. But at least, there is the idea that ‘weak learners’ can be used to provide a good predictor. Now, to be honest, I guess I understand the concept. But I still can’t reproduce what I got with standard ‘boosting’ packages.

There are a lot of publications about the concept of ‘boosting’. In 1988, Michael Kearns published Thoughts on Hypothesis Boosting, which is probably the oldest one. About the algorithms, it is possible to find some references. Consider for instance Improving Regressors using Boosting Techniques, by Harris Drucker. Or The Boosting Approach to Machine Learning An Overview by Robert Schapire, among many others. In order to illustrate the use of boosting in the context of regression (and not classification, since I believe it provides a better visualisation) consider the section in Dong-Sheng Cao’s In The boosting: A new idea of building models.

# Normalité des estimateurs dans une régression

Cette semaine en cours, on a évoqué la normalité des estimateurs dans une régression linéaire. En supposant la normalité des résidus, on a vu en cours que était un estimateur Gaussien. En particulier, chacun des estimateurs est alors Gaussien, au sens où , ce qui peut se visualiser sur le graphique suivant (la constante est en abscisse, et la pente en ordonnée), avec un intervalle de confiance à 95%,

En fait, la variance étant inconnue (mais pouvant être estimée), si on remplace la variance des résidus par leur estimateur, la loi de l’estimateur est une loi de Student. De même pour l’estimateur de la pente de la droite de régression,

Le fait que le couple soit également Gaussien fait que l’on peut construire non plus un intervalle de confiance (on n’est plus en dimension 1) mais une ellipse de confiance. On retient une forme elliptique car c’est la plus petite région dans laquelle on se trouvera avec une probabilité de 95% (comme discuté dans un vieux billet).

Mais pour mieux comprendre cette notion d’ellipse de confiance, le plus simple est d’avoir recours à du rééchantillonnage dans la base de données. En effet, comme on ne dispose que d’un jeu de données, on a un estimateur. Et la discussion sur la loi de notre estimateur est purement théorique. Pour visualiser la loi de notre estimateur, il faudrait des centaines, voire des milliers de bases similaires. Que l’on n’a pas. La solution est alors de tirer au hasard des points de notre échantillon (avec remise), i.e. de faire du bootstrap,

set.seed(1)
COEF=matrix(NA,10000,2)
for(s in 1:nrow(COEF)){
I=sample(1:nrow(cars),nrow(cars),replace=TRUE)
COEF[s,]=lm(dist~speed,data=cars[I,])$coefficients } Si on regarde ce qui se passe, tirage après tirage, on génère un paquet d’échantillon, et pour chaque échantillon, on ajuste une droite de régression. Les distributions – sur tous ces échantillons – de la constante et de la pente semblent effectivement Gaussiennes, hist(COEF[,1],col="light blue",prob=TRUE) u=seq(min(COEF[,1]),max(COEF[,1]),length=500) v=dnorm(u,mean(COEF[,1]),sd(COEF[,1])) lines(u,v,lwd=3,col="red") hist(COEF[,2],col="light blue",prob=TRUE) u=seq(min(COEF[,2]),max(COEF[,2]),length=500) v=dnorm(u,mean(COEF[,2]),sd(COEF[,2])) lines(u,v,lwd=3,col="red") pour la distribution de l’estimateur de la pente, alors que pour l’estimateur de la constante, on a la distribution suivante, Si on regarde maintenant la loi jointe, on retrouve une forme elliptique pour le nuage des points (et une forte corrélation négative entre les deux estimateurs). Mais si on creuse un peu library(ellipse) reg=lm(dist~speed,data=cars) e=ellipse(reg) plot(e,type="l",lwd=2) polygon(e,col="light blue") points(COEF,cex=.5) l’ellipse ne coïncide pas (tout à fait) avec celle obtenu théoriquement. Ce qui laisse à penser que l’hypothèse de normalité des résidus est peut-être à revoir.., # Quels résidus et quelle loi simuler ? Suite à une question sur les résidus, je vais essayer de prendre deux minutes pour reprendre un point que je n’avais pas abordé en cours, faute de temps. Les résidus les plus intéressants en provisionnement (et dans les GLM) sont les résidus de Pearson, Mais il existe également des résidus dits ajustés, permettant de prendre en compte le fait que le nombre de paramètres est ici relativement grand, ramené au nombre d’observations. On pose alors A première vue, utiliser l’un ou l’autre devrait donner la même chose si l’on génère des pseudo-triangles, puisque dans un cas, on utiliserait et dans l’autre Bref, multiplier par une constante pour ensuite diviser par la même constante, ça revient au même. Sauf qu’il peut sembler légitime de poser (on reviendra là dessus tout à l’heure). En attendant, voilà nos deux résidus > (E=residuals(regp,"pearson")) 1 2 3 4 9.488238e-01 2.404895e-02 1.168421e-01 -1.082940e+00 5 6 7 8 1.302749e-01 -1.007348e-13 -1.128013e+00 2.773332e-01 9 10 11 13 5.669707e-02 8.919633e-01 -2.110748e-01 -1.533031e+00 14 15 16 19 -2.213449e+00 -1.024162e+00 4.237393e+00 -4.899687e-01 20 21 25 26 7.929194e-01 -2.972380e-01 -4.275912e-01 4.140426e-01 31 -6.202125e-15 > n=sum(is.na(Y)==FALSE) > k=ncol(PAID)+nrow(PAID)-1 > (R=residuals(regp,"pearson")*sqrt(n/(n-k))) 1 2 3 4 1.374976e+00 3.485024e-02 1.693203e-01 -1.569329e+00 5 6 7 8 1.887862e-01 -1.459787e-13 -1.634646e+00 4.018940e-01 9 10 11 13 8.216186e-02 1.292578e+00 -3.058764e-01 -2.221573e+00 14 15 16 19 -3.207593e+00 -1.484151e+00 6.140566e+00 -7.100321e-01 20 21 25 26 1.149049e+00 -4.307387e-01 -6.196386e-01 6.000048e-01 31 -8.987734e-15 L’autre point sur lequel je voulais revenir sur les simulations, est que la loi de Poisson n’est peut-être pas adaptée pour simuler des scénarios de paiements futurs. En effet, si on fait une régression quasipoisson, on voit que le paramètre de surdispersion n’est pas négligeable, > regqp=glm(Y~as.factor(D)+as.factor(A), + data=base,family=quasipoisson(link="log")) > summary(regqp) Call: glm(formula = Y ~ as.factor(D) + as.factor(A), family = quasipoisson(link = "log"), data = base) Deviance Residuals: Min 1Q Median 3Q Max -2.3426 -0.4996 0.0000 0.2770 3.9355 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.05697 0.02769 290.995 < 2e-16 *** as.factor(D)2 -0.96513 0.02427 -39.772 2.41e-12 *** as.factor(D)3 -4.14853 0.11805 -35.142 8.26e-12 *** as.factor(D)4 -5.10499 0.22548 -22.641 6.36e-10 *** as.factor(D)5 -5.94962 0.43338 -13.728 8.17e-08 *** as.factor(D)6 -5.01244 0.39050 -12.836 1.55e-07 *** as.factor(A)2002 0.06440 0.03731 1.726 0.115054 as.factor(A)2003 0.20242 0.03615 5.599 0.000228 *** as.factor(A)2004 0.31175 0.03535 8.820 4.96e-06 *** as.factor(A)2005 0.44407 0.03451 12.869 1.51e-07 *** as.factor(A)2006 0.50271 0.03711 13.546 9.28e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Dispersion parameter for quasipoisson family taken to be 3.18623 Null deviance: 46695.269 on 20 degrees of freedom Residual deviance: 30.214 on 10 degrees of freedom (15 observations deleted due to missingness) AIC: NA Number of Fisher Scoring iterations: 4 On peut alors envisager de simuler non pas des lois de Poisson, mais des lois quasipoissons, en reprenant un vieux code, > rqpois = function(n, lambda, phi, roundvalue = TRUE) { + b = phi + a = lambda/phi + r = rgamma(n, shape = a, scale = b) + if(roundvalue){r=round(r)} + return(r) + } A partir de ces deux remarques, on peut reprendre le code qui permettait de générer des scénarios de paiements pour les années futures, > Yp=predict(regp,type="response",newdata=base) > Rs=rep(NA,20000) > for(s in 1:20000){ + serreur=sample(erreur, + size=36,replace=TRUE) + E=matrix(serreur,6,6) + sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6)) + sbase=data.frame(sY=as.vector(sY),D,A) + sbase$sY[is.na(Y)==TRUE]=NA
+ sreg=glm(sY~as.factor(D)+as.factor(A),
+ sYp=predict(sreg,type="response",
+ newdata=sbase)
+ sYpscenario=rqpois(36,sYp,phi=3.18623)
+ Rs[s]=sum(sYpscenario[is.na(Y)==TRUE])
+ }

si on regarde la densité de la distribution des paiements futurs, on obtient,

plot(density(Rs))

avec en bleu la distribution obtenue sur les résidus de Pearson brut, et en simulant une loi de Poisson, et en rouge la distribution des provisions avec les deux modifications. On observe clairement que les quantiles (par exemple) ont fortement changé. On peut le voir encore plus précisément sur un box-plot

Si on revient cinq minutes sur ce qu’on vient de faire. Dans le premier rappelons que l’on a besoin d’utiliser, pour prédire les paiements futurs, que

(avec la filtration correspondant à la partie supérieure du triangle disponible), i.e.

Comme on utilise un modèle de Poisson, on suppose aussi que

i.e.

c’est à dire que

Il faut avoir des résidus centrés et de variance unitaire. A première vue, c’est le cas pour nos résidus de Pearson (à peu près),

> mean(residuals(regp,"pearson"))
[1] -0.02462518
> sd(residuals(regp,"pearson"))
[1] 1.261934

sauf que l’estimateur de la variance est biaisé. Certes, classiquement, l’estimateur est basé sur une normalisation par un facteur  (classique en statistique quand on possède  observations)

> sqrt(sum((residuals(regp,"pearson")-
+ mean(residuals(regp,"pearson")))^2)/
+ (length(residuals(regp,"pearson"))-1))
[1] 1.261934

Sauf qu’ici on ne perd pas un degré de liberté (remplacer l’espérance des observations par leur moyenne empirique): en régression, il faut corriger par le nombre de variables explicatives. Et donc, pour avoir un estimateur sans biais de la variance de nos résidus, il faut corriger par un facteur . D’où l’utilisation des résidus dits ajustés afin d’avoir des résidus de variance (vraiment) unitaire.
Pour le second point, on utilise le fait qu’en pratique, la dispersion des paiements est plus grande que ce que donnerait un modèle de Poisson. Et il convient d’en tenir compte dans nos simulations. En l’occurrence, on veut que  , mais aussi que. On a vu que le paramètre de surdispersion ne pouvait être supposé unitaire. On a alors le choix, entre simuler une loi Gamma de paramètres

et

ou bien simuler une loi binomiale négative, de moyenne , et de paramètre de surdispersion  telle que la variance s’écrive

Aussi, on prend  alors que

Avec ces méthodes, on obtient quelque chose… qui est programmé sous R,

> BootChainLadder(PAID,10000,"od.pois")
BootChainLadder(Triangle = PAID, R = 10000,
process.distr = "od.pois")

Latest Mean Ultimate Mean IBNR SD IBNR IBNR 75% IBNR 95%
1  4,456         4,456       0.0     0.0        0        0
2  4,730         4,752      22.2    12.2       28       44
3  5,420         5,456      35.5    15.2       43       64
4  6,020         6,086      65.7    19.8       78      102
5  6,794         6,946     151.9    28.7      170      202
6  5,217         7,364   2,147.4   111.0    2,218    2,339

Totals
Latest:         32,637
Mean Ultimate:  35,060
Mean IBNR:       2,423
SD IBNR:           132
Total IBNR 75%:  2,506
Total IBNR 95%:  2,653

> quantile(Rs,c(.75,.95))
75%  95%
2509 2653

et qui donne les mêmes résultats que ce qu’on vient de reprogrammer…

# Simulations et incréments négatifs

Pour faire suite (rapidement) à mon dernier billet, le cas le plus fréquent pour observer des incréments négatifs est lorsque l’on fait des simulations pour obtenir la distribution du montant de paiements futurs. Rappelons que pour générer un pseudo-triangle, on utilise les résidus de Pearson, et on pose

On voit que l’on risque d’avoir un incrément négatif si le résidu est négatif, et plus grand (en valeur absolue) que la racine carrée de la prédiction (par notre modèle Poissonnien, en l’occurrence). Or dans le triangle vu en cours, on a

> source("https://perso.univ-rennes1.fr/
arthur.charpentier/bases.R")
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> Y=as.vector(INC)
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
> Yp=predict(reg,type="response",
+ newdata=base)
> erreurs=residuals(reg,"pearson")
> min(sqrt(Yp[is.na(Y)==FALSE]))
[1] 2.868171
> min(erreurs)
[1] -2.213449

autrement dit, on n’aura jamais d’incrément négatif lors de nos boucle, si l’on génère les résidus par bootstrap. Sinon, avec des lois paramétriques non-bornées inférieurement (loi normale, loi de Student), il est tout a fait possible d’obtenir des incréments négatifs. La méthode pour éviter le problème des incréments négatifs dans les boucles… est probablement de ne pas prendre en compte les scénarios où des incréments négatifs ont été obtenus, i.e.

R=rep(NA,10000)
for(s in 1:10000){
serreur=sample(erreurs,
size=36,replace=TRUE)
E=matrix(serreur,6,6)
sY=matrix(Yp,6,6)+E*sqrt(matrix(Yp,6,6))
if(min(sY[is.na(Y)==FALSE])>=0){
sbase=data.frame(sY=as.vector(sY),D,A)
sbase$sY[is.na(Y)==TRUE]=NA sreg=glm(sY~as.factor(D)+as.factor(A), data=sbase,family=poisson(link="log")) sYp=predict(sreg,type="response", newdata=sbase) R[s]=sum(sYp[is.na(Y)==TRUE])} } Une stratégie un peu plus propre pourrait être de mettre des 0 dès qu’on obtient un incrément négatif… Mais encore une fois, c’est de la cuisine. # Confidence interval for predictions with GLMs Consider a (simple) Poisson regression . Given a sample where , the goal is to derive a 95% confidence interval for given , where is the prediction. Hence, we want to derive a confidence interval for the prediction, not the potential observation, i.e. the dot on the graph below > r=glm(dist~speed,data=cars,family=poisson) > P=predict(r,type="response", + newdata=data.frame(speed=seq(-1,35,by=.2))) > plot(cars,xlim=c(0,31),ylim=c(0,170)) > abline(v=30,lty=2) > lines(seq(-1,35,by=.2),P,lwd=2,col="red") > P0=predict(r,type="response",se.fit=TRUE, + newdata=data.frame(speed=30)) > points(30,P1$fit,pch=4,lwd=3)

i.e.

Let denote the maximum likelihood estimator of . Then

where is Fisher information of (from standard maximum likelihood theory). Recall that

where computation of those values is based on the following calculations

In the case of the log-Poisson regression

Let us get back to our initial problem.

• confidence interval for the linear combination

A first idea to get a confidence interval for is to get a confidence interval for (by taking exponential values of bounds, since the exponential is a monotone function). Asymptotically, we know that

thus, an approximation for the variance matrix of will be based on , obtained by plugging estimators of the parameters.
Then, since as an asymptotic multivariate distribution, any linear combination of the parameters will also be normal, i.e.
has a normal distribution, centered on , with variance where is the variance of . All those quantities can be easily computed. First, we can get the variance of the estimators

> i1=sum(predict(reg,type="response"))
> i2=sum(cars$speed*predict(reg,type="response")) > i3=sum(cars$speed^2*predict(reg,type="response"))
> I=matrix(c(i1,i2,i2,i3),2,2)
> V=solve(I)

Hence, if we compare with the output of the regression,

> summary(reg)$cov.unscaled (Intercept) speed (Intercept) 0.0066870446 -3.474479e-04 speed -0.0003474479 1.940302e-05 > V [,1] [,2] [1,] 0.0066871228 -3.474515e-04 [2,] -0.0003474515 1.940318e-05 Based on those values, it is easy to derive the standard deviation for the linear combination, > x=30 > P2=predict(r,type="link",se.fit=TRUE, + newdata=data.frame(speed=x)) > P2$fit
1
5.046034

$se.fit [1] 0.05747075$residual.scale
[1] 1

> sqrt(V[1,1]+2*x*V[2,1]+x^2*V[2,2])
[1] 0.05747084
> sqrt(t(c(1,x))%*%V%*%c(1,x))
[,1]
[1,] 0.05747084

And once we have the standard deviation, and normality (at least asymptotically), confidence intervals are derived, and then, taking the exponential of the bounds, we get confidence interval

> segments(30,exp(P2$fit-1.96*P2$se.fit),
+ 30,exp(P2$fit+1.96*P2$se.fit),col="blue",lwd=3)

Based on that technique, confidence intervals are no longer centered on the prediction. But who cares ?

• delta method

Actually, those who like to use “more or less” expressions for confidence intervals will not like non centered intervals. So, an alternative is to use the delta method. Instead of writing (again) something on the theory, we can use a package which computes that method,

> estmean=t(c(1,x))%*%coef(reg)
> var=t(c(1,x))%*%summary(reg)$cov.unscaled%*%c(1,x) > library(msm) > deltamethod (~ exp(x1), estmean, var) [1] 8.931232 > P1=predict(r,type="response",se.fit=TRUE, + newdata=data.frame(speed=30)) > P1$fit
1
155.4048

$se.fit 1 8.931232$residual.scale
[1] 1

The delta method gives us (asymptotic) normality, so once we have a standard deviation, we get the confidence interval.

> segments(30,P1$fit-1.96*P1$se.fit,30,
+ P1$fit+1.96*P1$se.fit,col="blue",lwd=3)

Note that those quantities – obtained with two different approaches – are rather close here

> exp(P2$fit-1.96*P2$se.fit)
1
138.8495
> P1$fit-1.96*P1$se.fit
1
137.8996
> exp(P2$fit+1.96*P2$se.fit)
1
173.9341
> P1$fit+1.96*P1$se.fit
1
172.9101
• bootstrap techniques

And a third method (but far from what I expect to teach on that course) is to use bootstrap techniques to about those results based on asymptotic normality (we have only 50 observations). The idea is to sample from out dataset, and to run a log-Poisson regression on those new samples, and to repeat a lot of time,