# Delta method and quantile estimation

An alternative (to profile likelihood techniques) to derive confidence intervals is to use the delta method. Consider an estimator such that

then for any differentiable transformation

The proof of that result is based on Taylor’s expansion (see here or there for more details on the theory, or even on this blog – here, in French or there in English – for some codes in R). This can be used to derive an asymptotic confidence interval for a quantile. Consider the following dataset

```> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
> library(evir)
> X=base1\$Loss.in.DKM```

It is possible to fit a Generalized Pareto distribution on observations above a given threshold,

In that case, if exceed the threshold out of a sample of size , the estimator of the quantile

i.e. . Then

while

i.e. it is now possible to implement the delta-method to derive the asymptotic variance of the quantile, and also (asymptotic) confidence intervals.

```> u=5
> GPD=gpd(X,u)
> theta=GPD\$par.ests
> sigma=GPD\$varcov
> k=GPD\$n.exceed
> n=length(X)
> p=.975
> Q=u+theta[2]/theta[1]*((n*(1-p)/k)^(-theta[1])-1)
> nabla=c(-theta[2]/theta[1]^2*((1-p)^(-theta[1])-1)-
+ theta[2]/theta[1]*(1-p)^(-theta[1]*log(1-p)),
+ 1/theta[1]*((1-p)^(-theta[1])-1))
> variance=t(nabla)%*%sigma%*%nabla```

Based on the assumption of normality, it is possible to derive confidence intervals, and to compare them with the one obtained in R,

```> c(Q-1.96/sqrt(k)*sqrt(variance),
+   Q+1.96/sqrt(k)*sqrt(variance))
[1] 13.11562 16.82852
+ tailplot(gpd(X,5))
+ gpd.q(tailplot(gpd(X,5)), .975, ci.type =
+ "likelihood", ci.p = 0.95,like.num = 50)
Lower CI Estimate Upper CI
13.33329 14.97207 17.18094```

# 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
+ 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,

# Statistiques, de la théorie à la pratique (partie 2)

Il y a une quinzaine de jours (ici) on avait vu comment construire des estimateurs pour l’indice de queue de la loi de Pareto. On avait 4 valeurs (différentes). Maintenant, on peut se demander lequel devrait être le “meilleur“.

• propriétés de l’estimateur du maximum de vraisemblance

Maintenant que l’on a un peu avancé dans le cours, on peut continuer l’étude de la distribution des couts de sinistres (étude commencée ici).

i.e.

de telle sorte que

et

Aussi l’information de Fisher est .
La variance aymptotique de  est alors

• propriétés de l’estimateur de la méthode des moments

Nous avions noté que

i.e. , alors . Aussi, comme le théorème central limite limite

la delta method permet d’écrire, en posant ,

i.e.

La variance asymptotique est alors

• propriétés de l’estimateur de la méthode de la médiane

Pour le dernier estimateur, le théorème de Glivenko Cantelli sur les fonctions quantiles permet d’écrire

i.e. pour la médiane,

où la variance asymptotique (pour la médiane) peut se simplifier sous la forme

Là encore, en utilisant la delta method (je peux renvoyer ici pour une application numérique de cette méthode), avec

i.e. , on en déduit que

La variance asymptotique de notre estimateur est alors (en enlevant le facteur au dénominateur – car il est commun à tous),

Pour trois de nos estimateurs (je n’ai pas fait le quatrième, ça serait hors programme), on a pu obtenir facilement des variances asymptotiques. Graphiquement, on a les évolutions suivantes pour les variances asymptotiques de nos trois estimateurs, en fonction de la valeur de , avec le maximum de vraisemblance (en rouge), celui de la méthode des moments (en bleu) et celui basé sur la médiane (en mauve)

Manifestement le premier estimateur (le maximum de vraisemblance) a une variance asymptotique plus faible que les deux autres (mais on l’avait vu en cours, c’est ce que garantit l’inégalité de Cramér-Rao). Autrement dit, l’estimateur du maximum de vraisemblance est probablement le plus fiable….
En utilisant ces résultats mentionnés, on peut aussi construire des intervalles de confiance asymptotiques pour la valeur du paramètre , à l’aide de nos trois estimateurs, i.e. pour un intervalle de confiance à 90%,

où les estimations des variances asymptotiques sont obtenues en substituant à   dans les formules au dessus.
Si on considère que l’estimateur du maximum de vraisemblance est le meilleur que l’on puisse utiliser, on obtient comme intervalle de confiance pour ,

```2.845093-1.64/sqrt(100)*2.845093
[1] 2.378498

2.845093+1.64/sqrt(100)*2.845093
[1] 3.311688```

L’étape d’après sera d’introduire des tests. On pourra ainsi se demander si l’on peut accepter l’hypothèse . Mais c’est une autre histoire… à suivre donc….

# Delta method avec R, ou transformation des paramètres

Un petit billet aujourd’hui sur une question qui m’a été posée (et qui m’a permis de découvrir des nouvelles fonctions sous R). On suppose que l’on a un modèle économétrique de la forme

que l’on peut encore écrire

Supposons que ce soit la formulation en terme de ces derniers coefficients qui nous intéresse, et que l’on cherche un intervalle de confiance.
La Delta Method nous dit que si on dispose d’un estimateur asymptotiquement gaussien (mais c’est le cas pour la plupart des estimateurs que l’on construit en pratique), tout transformation bijective de cet estimateur reste asymptotiquement Gaussien. Bref, on va pouvoir utiliser ces résultats pour en déduire un intervalle de confiance (asymptotique) pour le nouveau paramètre. En effet, l’estimateur par moindre carrés vérifie la normalité asymptotique (sous quelques hypothèses)

ce qui donne avec les notations précédentes

Autrement dit, on peut facilement en déduire un intervalle de confiance (asymptotique). Pour cela, on a trois méthodes (au moins) avec R. Pour la suite, je vais considerer la régression sur la base de données cars (comme toujours), et on s’intéresse à l’estimation de

(ne me demandez pas pourquoi, il me fallait un exemple numérique, j’ai pris un truc nonlinéaire pour qu’au moins ça puisse avoir un vague intérêt).

• En faisant un peu de calcul formel

Pour les calculs symboliques, il y a la fonction D dans le package de base,

```> f <- expression(sin(cos(x + y^2)))
> ( Dfx <- D(f, "x") )
-(cos(cos(x + y^2)) * sin(x + y^2))```

ce qui pourra alors être utilisé pour faire les calculs de la variance asymptotique.

```> x=c(0,1,2); y=1
> eval(D.sc)
[1] -0.7216061 -0.8316919 -0.0774320```

On peut également utiliser cette fonction pour faire plein de jolis dessins,

```> x <-seq(-5, 5, length=100)
> fx=x^2
> Df <-D(expression(x^2), "x")
> Dfx<-eval(Df)
> plot(x, fx, type = "l", las = 1,ylim=c(-2,2))
> for (i in 1:length(x)) {
+ abline(a = fx[i] - Dfx[i] * x[i], b = Dfx[i],
+ col = rainbow(100)[i])
+ }```

Mais bon, c’était juste pour rigoler un peu… On peut utiliser cette fonction pour calculer l’écart type de l’estimateur de la grandeur qui nous intéresse.

```> reg1 = lm(dist~speed, data = cars)
> estmean = coef(reg1)
> estvar = summary(reg1)\$cov.unscaled * summary(reg1)\$sigma^2
> f = function(x,y){
+ -y/(x+y)
+ }
> X=f(estmean[1],estmean[2])
> g=expression(-y/(x+y))
> x<-as.numeric(estmean[1]); y<-as.numeric(estmean[2])
> V=eval(c(eval(D(g, "x")),eval(D(g, "y"))))
> sqrt(t(V) %*% estvar %*% V)
[,1]
[1,] 0.1063255```

Voilà pour la première méthode. Petite remarque: pour faire le calcul symbolique de la dérivée, j’utilise expression et non pas function.

• En utilisant la fonction deltamethod de la librarie msm

Ici on commence comme avant, sauf que la phase de calcul est immédiat…

```> reg2 <- lm(dist~speed, data = cars)
> estmean = coef(reg2)
> estvar = summary(reg2)\$cov.unscaled * summary(reg2)\$sigma^2
> deltamethod (~ -x2 / (x1 + x2), estmean, estvar)
[1] 0.1063255```

et d’ailleurs, ça peut s’appliquer pour un paramètre en dimension 2

```> deltamethod (list(~ -x2 / (x1 + x2), ~ x1+x2), estmean, estvar)
[1] 0.1063255 6.3664368```

En revanche, on obtient des estimateur des écart-types, et pas de la matrice de variance-covariance des estimateurs. Autrement dit on ne peut pas en déduire l’ellipse de confiance, seulement des intervalles de confiance marginaux (en utilisant la normalité asymptotique).

• En utilisant la fonction delta.method de librarie alr3

Pour finir, une dernière fonction. La seule chose qui change c’est la synthaxe. Mais la sortie est plus élégante,

```> reg3 = lm(dist~speed, data = cars)
> delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255```

En fait, l’écart-type indiqué ici est précisément celui qui est utilisé pour écrire l’intervalle de confiance. Dans cet exemple, l’intervalle de confiance approché (car on suppose que la normalité – en théorie seulement asymptotique – est vérifiée) et estimé (car les paramètres sont inconnus et estimés) est de la forme suivant, pour un seuil de confiance à 95%

```> DM=delta.method(reg3, "-b1/(b1+b0)")
Functions of parameters:  expression(-b1/(b1+b0))
Estimate = 0.2881585 with se = 0.1063255
> DM\$estimate+qt(.975,df=nrow(cars)-1)*DM\$se
[1] 0.5018277
> DM\$estimate+qt(.025,df=nrow(cars)-1)*DM\$se
[1] 0.07448936```

C’est à dire qu’on a un intervalle de confiance de la forme [0.074;0.502].