Estimation en régression de Poisson

Stéphane Ivanoff

12 octobre 2012

Sous la direction de Vincent Rivoirard

Table des matières

1 Introduction 2

2 Présentation du modèle 2
2.1 Régression de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.2 Le Lasso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.3 Estimateur de Dantzig . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3 Choix des λj 4
3.1 Processus de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.2 Un théorème pour la calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.3 Inégalité oracle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4 Perspectives 7
4.1 De meilleures inégalités oracles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
4.2 Group-Lasso et Group-Dantzig . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1
1 Introduction
En statistiques paramétriques, on cherche, de manière générale, à estimer un nombre p de para-
mètres à partir d’un nombre n ≥ p d’observations. Typiquement, en régression linéaire gaussienne,
la variable expliquée Y dépend linéairement de paramètres X1 , ..., Xp selon des coefficients inconnus
β1 , ..., βp , avec de plus un terme d’erreur  qu’on suppose gaussien. On observe n représentations
indépendantes de ce modèle, d’où, en notant Xij le j-ième paramètre de l’observation Yi , le modèle :

Y = Xβ + .

Un résultat classique affirme

 que le meilleur estimateur de β est la solution des moindres carrés
β̂ = argmin kY − Xβk22 , qui vaut β̂ = (X ∗ X)−1 X ∗ Y .
Cependant, que se passe-t-il lorsque p > n (ce que l’on appelle les statistiques en grande di-
mension) ? Le résultat précédent n’est plus applicable car la matrice X ∗ X n’est plus inversible. Cela
traduit en fait un problème plus général : on n’a pas de solution unique à la minimisation de kY −Xβk22
(intuitivement, on peut comparer cette situation à celle d’un système de n équations linéaires à p in-
connues). Il est donc nécessaire de développer d’autres méthodes pour l’estimation de β. Généralement,
ces méthodes ont recours à la pénalisation pour effectuer un choix. L’exemple le plus classique est
celui du Lasso, défini ainsi :
n o
β̂ L = argmin kY − Xβk22 + λkβk1 .

Cet estimateur est une référence en statistiques en grande dimension. Nous reviendrons plus loin sur ces
propriétés ; le lecteur intéressé par le Lasso et les statistiques en grandes dimensions en général pourra
consulter [3]. Ce problème, loin d’être purement théorique, est au contraire au coeur de nombreuses
problématiques modernes. Dans l’étude du génome humain, par exemple, on peut chercher à déterminer
l’impact de milliers de gènes sur un phénomène particulier, alors qu’on ne dispose que de quelques
dizaines d’observations.

2 Présentation du modèle
2.1 Régression de Poisson
Pour la suite de cet article, nous allons considérer non plus la régression linéaire gaussienne mais
la régression de Poisson, qui s’écrit :

Yi |Xi ∼ P(f0 (Xi )),

où les Xi peuvent être aléatoires ou déterministes. On les supposera dans [0, 1].
En d’autres termes, on connait les Xi et les Yi et on suppose que Yi suit une loi de Poisson dont le
paramètre est une fonction déterministe inconnue f0 de Xi . Cette fonction étant strictement positive,
on peut écrire f0 = exp(g0 ) et on va supposer que g0 se décompose sur un dictionnaire, c’est-à-dire
qu’il existe une famille {ϕ1 , ..., ϕp } connue de fonctions telle que

X
p
g0 (x) = β0j ϕj (x).
j=1

2
Estimer f0 revient à estimer les coefficients β0j . On cherchera donc un estimateur de f0 sous la forme
X
p 
fβ (x) = exp βj ϕj (x) , β ∈ Rp .
j=1

Les βj sont appelés coefficients de régression.

2.2 Le Lasso
On s’intéresse à la log-vraisemblance du modèle. La vraisemblance vaut :
Y
n
e−fβ (Xi ) fβ (Xi )Yi
L(β) =
Yi !
i=1

d’où on déduit la log-vraisemblance :

X
n
l(β) = (Yi log fβ (Xi ) − fβ (Xi ))
i=1

à une constante près.

On définit C(β) = −l(β). Maximiser la log-vraisemblance revient donc à minimiser C. Cependant,

les solutions sont multiples à cause de l’hypothèse p > n, comme établi précédemment. On ajoute
P
donc une pénalité sur les coefficients de β : on définit pen(β) = pj=1 λj |βj | et l’estimateur suivant :
n o
β̂ L = argmin C(β) + pen(β) .

Les λj sont des constantes strictement positives laissées arbitraires pour le moment.
Cet estimateur est connu sous le nom de Lasso (Least Absolute Shrinkage and Selection Operator),
introduit par Tibshirani en 1996. L’ajout de la norme L1 pondérée dans la quantité à minimiser permet
d’assurer l’existence d’une unique solution.

2.3 Estimateur de Dantzig

On définit maintenant un autre estimateur pour les grandes dimensions, l’estimateur de Dantzig.
Il vient de l’approche suivante : en posant Aij = ϕj (Xi ), on peut réécrire la log-vraisemblance ainsi :
X
p X
n X
p


∗
l(β) = βj (A Y )j − exp βj Aij .
j=1 i=1 j=1

Ses dérivées partielles s’écrivent alors :

∂l(β)


= A∗j Y − exp Aβ ,
∂βj
où Y = (Y1 , ..., Yn )∗ et exp(Aβ) = (exp((Aβ)1 ), ..., exp((Aβ)n ))∗ . La condition qu’on voudrait imposer
est ∂l(β) ∗
∂βj = 0 mais c’est impossible. On relaxe donc cette condition en |Aj (Y − exp(Aβ))| ≤ λj , et on
définit : n o
β̂ D = argmin kβk1 : ∀j, |X ∗ (Y − exp(Aβ))| ≤ λj .

3
 Il est intéressant de comparer ces deux estimateurs, d’autant plus qu’ils font intervenir les mêmes
λj . L’estimateur de Dantzig est plus difficile à manipuler que le Lasso, mais d’un point de vue algo-
rithmique, il est plus rapide à mettre en place. Tous deux présentent cependant l’avantage d’estimer
une partie des coefficients précisément à zéro, à cause de la géométrie de la norme L1 .

{|A∗j (Y − exp(Aβ))| ≤ λj ∀j}

ligne de niveau de kβk1

β̂ D

Cette propriété est utile en pratique. En effet, dans l’exemple de l’étude du génome, de nombreux
gènes sont considéres dans le modèle mais n’ont en fait aucune influence sur le phénomène considéré.
Il est d’ailleurs classique de supposer que le vecteur β0 est S-sparse, c’est-à-dire que seuls S de ses
coefficients sont non nuls (typiquement, S est bien plus petit que p).
L’estimateur Lasso comme l’estimateur de Dantzig permettent ainsi de faire le tri entre les gènes
influents (coefficients non nuls) et les autres ; on parle de sélection de modèles.

3 Choix des λj
Dans cette section, on cherche comment choisir les λj de manière optimale : on parle de calibra-
tion.
Cette calibration des λj fait intervenir la notion de processus de Poisson, qu’on rappelle ou présente
ici.

3.1 Processus de Poisson

Soit X un ensemble mesurable. Un processus ponctuel N sur X est un ensemble aléatoire au
plus dénombrable de points de X.
Pour A ⊂ X mesurable, on note NA le nombre de points de N dans A.

On dit alors que N est un processus de Poisson de mesure moyenne µ si :

– ∀n, ∀A1 , ..., An ⊂ X disjoints mesurables, NA1 , ..., NAn sont des variables aléatoires indépendantes.
– ∀A ⊂ X mesurable, NA ∼ P(µ(A)).  
Remarque. µ ne peut avoir d’atomes. En effet, si µ({x}) > 0, alors P N{x} ≥ 2 > 0, or il ne
peut pas y avoir de répétition dans N car c’est un ensemble.

4
 Supposons maintenant X = R+ . Alors dµ = λ(x)dx et on appelle λ l’intensité du processus. Si λ
est une constante, on est dans le cas plus simple d’un processus de Poisson homogène sur R+ , qu’on
peut concevoir, et construire, comme une suite de variables aléatoires indépendantes exponentielles de
paramètre λ :
Soit (τn )n∈N une suite de variables aléatoires i.i.d. ∼ E(λ). On pose Ti = τ1 + ... + τi et le processus
ponctuel {T1 , ..., Tn , ...} est un processus de Poisson d’intensité λ.
P
Dans le cas plus général d’une fonction λ non constante, on définit la mesure ponctuelle dN = t∈N δt ,
de sorte que pour toute fonction f ,
Z X
f (x)dNx = f (t).
t∈N

On a alors les formules suivantes sur l’espérance et la variance, qui lient le processus de Poisson N à
sa mesure moyenne µ :
Z Z  Z
si f dµ < ∞, E f dN = f dµ.
Z Z  Z
si f 2 dµ < ∞, Var f dN = f 2 dµ.

On dispose alors d’inégalités de concentrations. On présente la suivante, de “type” Bernstein :

Lemme 1. Soit f une fonction bornée. Alors, ∀x > 0,

Z !
√ kf k∞ x
P f (dN − dµ) ≥ 2vx + ≤ e−x , (1)
3
R  R
où v = Var f dN = f 2 dµ.

Ces notions, et en particulier ce lemme, permettent d’établir un théorème qui énonce le meilleur
choix de λj .

3.2 Un théorème pour la calibration

Les λj doivent être choisis en fonction des données du problème. On va voir ici comment réaliser
ce choix. Pour le moment, on va s’intéresser uniquement à l’estimateur de Dantzig.
On voudrait que β0 satisfasse à la condition de Dantzig : ∀j, |X ∗ (Y − exp(Aβ))| ≤ λj .
Comme exp(Aβ0 ) = E[Y ], cela revient à chercher les λj tels que pour tout j,

|A∗j (Y − E[Y ])| ≤ λj

avec forte probabilité.

λj ne doit donc pas être pris trop petit. Cependant, s’il est pris trop grand, on risque d’estimer trop
de paramètres à zéro. Un compromis est donc nécessaire.
La calibration se fait via le résultat suivant. On définit
X
n
Vj = Var(A∗j Y ) = f0 (Xi )ϕ2j (Xi )
i=1

5
et son estimateur
X
n
V̂j = Yi ϕ2j (Xi ).
i=1

Théorème 1. Soit γ > 1. On pose

q
Ṽj = V̂j + 2γ(log p)kϕj k2n,∞ V̂j + 3γ(log p)kϕj k2n,∞ V̂j

où kϕj kn,∞ = maxi∈{1,...,n} |ϕj (Xi )|.

Alors
  3
P |A∗j (Y − E[Y ])| > λj ≤ γ , (2)
p
q
où λj = 2γ log pṼj + γ log p
3 kϕ̃j k∞ .

γ est un paramètre laissé libre. En pratique, on obtient de bons résultats en le choisissant proche
de 1.
La preuve se fait en ramenant la quantité à contrôler |A∗j (Y −E[Y ])| à une quantité faisant intervenir
un processus de Poisson, de la manière suivante. On revient à notre modèle de régression de Poisson.
Quitte à les renuméroter, on peut supposer que Xi < Xi+1 pour tout i. On définit alors la fonction h
constante par morceaux qui vaut
f0 (Xi )
h(t) =
Xi+1 − Xi
pour t ∈ [Xi , Xi+1 [.
Soit N un processus de Poisson d’intensité h, donc de mesure moyenne µ = h(t)dt. Des calculs simples
permettent de lier la quantité à contrôler |A∗j (Y − E[Y ])| au processus de Poisson N :
Z
|A∗j (Y − E[Y ])| = ϕ̃j (t)(dN (t) − h(t))dt .

Le lemme 1 permet de contrôler cette quantité en fonction de la variance Vj , puis de contrôler Vj par
la quantité Ṽj . On obtient ainsi une inégalité valable avec probabilité 1 − 3e−x . Le choix x = γ log p
permet d’aboutir à (2).

3.3 Inégalité oracle

Il reste encore à montrer l’efficacité de l’estimateur de Dantzig pour les λj choisis. On construit
pour cela une inégalité oracle.
Il nous faut d’abord définir certains coefficients. Soit X une matrice de taille n × p. Pour tout en-
semble T ⊂ {1, ..., p}, on note XT la matrice de taille n × |T | obtenue en prenant les colonnes de X
correspondant aux indices de l’ensemble T .
Pour tout S ≤ p, on définit δSX comme étant la plus grande quantité telle que, pour tout ensemble
T ⊂ {1, ..., p} de taille |T | ≤ S et tout vecteur c de longueur |T |,

kXT ck22 ≥ δSX kck22 .

6
Pour tous S et S 0 tels que S + S 0 ≤ p, on définit de plus θS,S X
0 comme la plus petite quantité telle que,

pour tous sous-ensembles disjoints T et T de {1, ..., p} et tous vecteurs c et c0 de longueurs respectives
0

|T | ≤ S et |T 0 | ≤ S 0 ,
|(XT c)∗ (XT 0 c0 )| ≤ θS,S
X 0
0 kck2 kc k2 .
p p
Soit à la matrice de terme général Ãij = f0 (Xi )Aij = f0 (Xi )ϕj (Xi ). Dans la suite, on note
δS = δSA et δ̃S = δSÃ , θS,S 0 = θS,S
A Ã
0 et θ̃S,S 0 = θS,S 0 .

On a alors le théorème suivant :

Théorème 2. On suppose :
– β0 est S-sparse : l’ensemble T0 = {j/ β0j 6= 0} est de taille S ;
– 3S ≤ p ;
– θ̃S,2S < δ̃2S .
Alors, en choisissant les λj comme dans la section précédente, on a, avec probabilité 1 − O( pγ−1
1
),
P
j∈T01 V̂j
kβ̂ −
D
β0 k22 ≤ 256γ log p ,
(δ̃2S − θ̃S,2S )2

où T01 = T0 ∪ T1 avec T1 l’ensemble des S indices j ∈

/ T0 correspondant aux S plus grandes valeurs de
β̂.

Remarque. L’hypothèse θ̃S,2S < δ̃2S n’est pas vérifiable car elle porte sur la matrice à inconnue.
Si on suppose f0 continue, il existe alors des constantes m et M strictement positives telles que
m ≤ f0 ≤ M sur [0, 1] et on peut remplacer notre hypothèse par celle-ci, plus forte :

M θS,2S < mδ2S .

On remplace également δ̃2S − θ̃S,2S par mδ2S − M θS,2S dans l’inégalité.

4 Perspectives
Il reste encore beaucoup à faire sur l’estimation en régression de Poisson, d’autant plus que c’est
un modèle qui pourrait devenir assez actuel maintenant que certains phénomènes biologiques, jusqu’ici
modélisés de manière gaussienne, sont mieux représentés par des processus de Poisson. Voici quelques
perspectives pour l’avenir.

4.1 De meilleures inégalités oracles

Le théorème 2 n’est pas optimal. Il porte sur le carré de l’erreur L2 d’estimation de β0 , là où on
souhaiterait plutôt un résultat sur l’estimation de f0 .
De plus, pour être véritablement applicable, il faut supposer non seulement que f0 est bornée (ce qui
n’est pas une grosse contrainte), mais aussi qu’on a des informations sur ses bornes, ce qui est bien
plus fort en pratique.
Il est probablement nécessaire de poser des conditions plus précises sur le dictionnaire (ϕj )j∈{1,...,p} .
Une meilleure majoration peut être trouvée si on prend pour dictionnaire une base orthonormée clas-
sique (base de Haar, ondelettes...).

7
 P
Une autre piste, à base de l’estimateur Lasso,  une inégalité oracle plus propre. On
permet d’obtenir
p
rappelle que pour tout β, on note fβ (x) = exp j=1 βj ϕj (x) . On peut alors définir la divergence
de Kullback :

K̃n (β0 , β) = EPβ0 [l(f0 ) − l(fβ )]

Xn h   i
= f0 (Xi ) log f0 (Xi ) − log fβ (Xi ) − (f0 (Xi ) − fβ (Xi )) .
i=1

Alors un résultat (de [6]) affirme que, pour tout ζ > 0, on peut poser a0 = 3 + ζ4 et alors, sous
l’hypothèse RE (s, a0 ) qui demande essentiellement que la matrice Ã∗ Ã vérifie une forme de “défini-
positivité” restreinte aux vecteurs b vérifiant kbJ c k1 ≤ a0 kbJ k1 pour tout sous-ensemble J de {1, ..., p}
de taille inférieure ou égale à s, on a :
( )
K̃n (β0 , β̂) ≤ (1 + ζ) inf K̃n (β0 , β) + C(ζ, s)|S(β)|(max λj )2 .
β∈Rp ,|S(β)|≤s j

Il s’agit ici d’une véritable inégalité oracle, qui affirme que l’estimateur β̂ approche β0 aussi bien
que le meilleur β, à un terme additif proportionnel à la taille du support de β, et à un facteur (1 + ζ)
près.
Malheureusement, ce théorème demande pour l’instant une hypothèse supplémentaire trop forte (l’équation
(12) de [6]). Il semble néanmoins possible d’améliorer cette hypothèse.

4.2 Group-Lasso et Group-Dantzig

Toujours dans l’exemple de l’étude du génome, il s’avère en réalité que les gènes n’agissent pas seuls,
mais en groupe. On cherche donc non plus à sélectionner des gènes, mais des groupes indissociables
de gènes.
On définit donc dans cette optique les estimateurs Group-Lasso et Group-Dantzig. Soit {G1 , ..., GM }
une partition de {1, ..., p}. On pose alors
( )
X
M
gL
β̂ = argmin C(β) + λk kβGk k2
k=1

et ( )
X
M
β̂ gD = argmin kβGk k2 : ∀k ∈ {1, ..., M }, kA∗Gk (Y − exp(Aβ))k2 ≤ λk .
k=1

On remarque que si M = p, alors pour tout k, quitte à effectuer une permutation, Gk = k, kβGk k2 =
|βk | et on retrouve simplement les estimateurs Lasso et de Dantzig classiques.

8
 Je remercie Vincent Rivoirard pour son aide et sa patience tout au long de mon stage, ainsi que
Patricia Reynaud-Bouret pour son accueil chaleureux à l’université de Nice.

Références
[1] Antoniadis, A., Fryzlewicz, P. and Letué, F. (2010) The Dantzig selector in Cox’s propor-
tional hazards model. Scand. J. Stat. 37(4), 531–552.
[2] Bertin, K., Le Pennec, E. and Rivoirard, V. (2011) Adaptive Dantzig density estimation,
Annales de l’Institut Henri Poincaré Probabilités et Statistiques, 47(1), 43–74.
[3] Bühlmann, P. and van de Geer, S. (2011) Statistics for high-dimensional data. Springer
Series in Statistics. Springer, Heidelberg.
[4] Candès, E. and Tao, T. (2007) The Dantzig selector : statistical estimation when p is much
larger than n, Ann. Statist. 35(6), 2313–2351.
[5] James, G.M. and Radchenko, P. (2009) A generalized Dantzig selector with shrinkage tun-
ing. Biometrika 96(2), 323–337.
[6] Lemler, S. (2012) Oracle inequalities for the Lasso for the conditional hazard rate in a
high-dimensional setting.
[7] Kingman, J.F.C. (1993) Poisson processes. Oxford studies in Probabilites.
[8] Reynaud-Bouret, P. and Rivoirard, V. (2010) Near optimal thresholding estimation of a
Poisson intensity on the real line. Electronic Journal of Statistics, 4, 172–238.