Estimation de courbes moyennes de consommation électrique pour des petits domaines à partir d’échantillons
Section 4. Méthodes d’estimation dans l’approche basée sur un modèle

• Précédent
• Suivant

Dans cette section, nous nous plaçons dans l’approche basée sur le modèle Valliant, Dorfman et Royall (2000) dans laquelle les courbes $Y_i$ sont considérées aléatoires et nous proposons quatre approches innovantes permettant de répondre à notre problématique d’estimation de courbes de charge moyennes de petits domaines. Supposons que $Y_i\left(t\right)$ et le vecteur d’information auxiliaire $X_i$ sont disponibles pour chaque individu $i$ du domaine $d$ et en plus, que la moyenne ${\overline{X}}_d={\displaystyle {\sum }_{i\in U_d}}{X}_i/N_d$ est connue.

On suppose que les variables auxiliaires sont liées aux courbes de charge selon un modèle fonctionnel de superpopulation sur l’ensemble de la population qui s’écrit de manière générale :

$$\xi \mathrm{<mo>:</mo>}{Y}_i\left(t\right)\mathrm{<mo>=</mo>}{f}_d\left(X_i\mathrm{,}t\right)+{ϵ}_i\left(t\right)\mathrm{,}i\in U_d\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}(4.1)$$

avec $f_d$ une fonction de régression inconnue à estimer, qui peut varier d’un domaine à l’autre et ${ϵ}_i$ un processus de bruit d’espérance nulle, de covariance nulle pour des individus différents et non-nulle par rapport au temps.

Si la taille du domaine $N_d$ est grande, alors la moyenne ${\mu }_d$ sera estimée par

$${\widehat{\mu }}_d\left(t\right)\mathrm{<mo>=</mo>}\frac{1}{N_d}{\displaystyle \sum _{i\in U_d}}{\widehat{Y}}_i\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}$$

où ${\widehat{Y}}_i\left(t\right)\mathrm{<mo>=</mo>}{\widehat{f}}_d\left(X_i\mathrm{,}t\right)$ est la prédiction de $Y_i\left(t\right).$ Dans le cas contraire, la moyenne ${\mu }_d$ est estimée par (voir Valliant et coll., 2000) :

$${\widehat{\mu }}_d\left(t\right)\mathrm{<mo>=</mo>}\frac{1}{N_d}\left({\displaystyle \sum _{i\in s_d}}{Y}_i\left(t\right)+{\displaystyle \sum _{i\in U_d-s_d}}{\widehat{Y}}_i\left(t\right)\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{.}(4.2)$$

La qualité de nos estimations dépend alors de la qualité de notre modèle : si le modèle est faux, cela pourra conduire à des biais dans les estimations.

4.1  Modèle linéaire fonctionnel

Le plus simple modèle de la forme (4.1) est le modèle de regression linéaire fonctionnelle Faraway (1997) :

$$Y_i\left(t\right)=X_i^{\prime }\beta \left(t\right)+{\epsilon }_i\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}i\in U_d\mathrm{.}(4.3)$$

où les résidus ${\epsilon }_i\left(t\right)$ sont indépendants pour $i\ne j,$ distribués selon une loi de moyenne 0 et de variance ${\sigma }_i^2\left(t\right)\mathrm{.}$ Si la taille du domaine $N_d$ est large, alors la moyenne de $Y$ dans le domaine $d$ est estimée par

$${\widehat{\mu }}_d^{\text{blu}}\left(t\right)\mathrm{<mo>=</mo>}{\overline{X}}_d^{\prime }{\widehat{\beta }}_{\text{BLU}}\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}$$

où ${\widehat{\beta }}_{\text{BLU}}\left(t\right)\mathrm{<mo>=</mo>}{\left({\displaystyle {\sum }_{i\in s}}{X}_i{X}_i^{\prime }/{\sigma }_i^2\left(t\right)\right)}^{-1}{\displaystyle {\sum }_{i\in s}}{X}_i{Y}_i\left(t\right)/{\sigma }_i^2\left(t\right)$ est le «Best Linear Unbiased» (BLU) estimateur de $\beta$ qui ne dépend pas du domaine $d\mathrm{.}$ L’estimateur ${\widehat{\mu }}_d^{\text{blu}}\left(t\right)$ peut s’écrire comme une somme pondérée de $Y_i\left(t\right)\mathrm{<mo>:</mo>}$

$${\widehat{\mu }}_d^{\text{blu}}\left(t\right)\mathrm{<mo>=</mo>}\frac{1}{N_d}{\displaystyle \sum _{i\in s}}{w}_{id}^{\text{blu}}\left(t\right)Y_i\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}$$

où les poids $w_{id}^{\text{blu}}\left(t\right)\mathrm{<mo>=</mo>}\left({\displaystyle {\sum }_{j\in U_d}}{X}_j^{\prime }\right){\left({\displaystyle {\sum }_{j\in s}}{X}_j{X}_j^{\prime }/{\sigma }_j^2\left(t\right)\right)}^{-1}{X}_i/{\sigma }_i^2\left(t\right)$ dependent maintenant du temps $t\mathrm{.}$ Si $n_d/N_d$ n’est pas négligeable, alors la moyenne ${\mu }_d$ est estimée en utilisant (4.2) par :

$${\widehat{\mu }}_d^{\text{blu}}\left(t\right)\mathrm{<mo>=</mo>}\frac{1}{N_d}{\displaystyle \sum _{i\in s_d}}\left(Y_i\left(t\right)-X_i^T{\widehat{\beta }}_{\text{BLU}}\left(t\right)\right)+{\overline{X}}_d^{\prime }{\widehat{\beta }}_{\text{BLU}}\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{.}$$

Cet estimateur peut s’écrire encore comme une somme pondérée des $Y_i\left(t\right)$ avec des poids qui dépendront toujours du temps $t\mathrm{.}$ La fonction de variance (sous le modèle) de ${\widehat{\mu }}_d^{\text{blu}}\left(t\right)$ peut être dérivée en utilisant Rao et Molina (2015), chapitre 7. La fonction de variance ${\sigma }_i^2\left(t\right)$ est inconnue et elle peut-être estimée en suivant Rao et Molina (2015). En remplaçant ${\sigma }_i^2\left(t\right)$ par ${\widehat{\sigma }}_i^2\left(t\right)\mathrm{,}$ nous obtiendrons l’estimateur EBLUP « Empirical Best Linear Unbiased Predictor » de ${\mu }_d$ et sa variance peut-être obtenue en suivant la méthode donnée dans Rao et Molina (2015). Cet estimateur EBLUP n’utilise pas les poids de sondage $d_i$ et par conséquent, il n’est pas consistant par rapport au plan de sondage (sauf si les poids de sondage sont constants pour les unités appartenant au même domaine $d).$ Un estimateur modifié, appelé aussi pseudo-EBLUP, peut-être construit en suivant de nouveau l’approche décrite en (Rao et Molina, 2015, chapitre 7) et qui sera égal dans ce cas à l’estimateur donné dans (3.4).

Si $Y_i\left(t\right)$ n’est pas connu pour les unités du domaine $d\mathrm{,}$ l’estimateur indirect suivant peut être utilisé :

$${\widehat{\mu }}_d^{\text{ind}}\left(t\right)\mathrm{<mo>=</mo>}{\overline{X}}_d\widehat{\beta }\left(t\right)\mathrm{<mo>=</mo>}\frac{1}{N_d}{\displaystyle \sum _{i\in s}}{\tilde{w}}_{id}^{\text{ind}}{Y}_i\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}(4.4)$$

avec $\widehat{\beta }\left(t\right)$ donné dans (3.5) et les poids ${\tilde{w}}_{id}^{\text{ind}}\mathrm{<mo>=</mo>}\left({\displaystyle {\sum }_{j\in U_d}}{X}_j^{\prime }\right){\left({\displaystyle {\sum }_{i\in s}}{d}_i{X}_i{X}_i^{\prime }\right)}^{-1}{\displaystyle {\sum }_{i\in s}}{d}_i{X}_i$ ne dépendent pas de temps $t$ contrairement à $w_{id}^{\text{blu}}\text{}\mathrm{.}$ Ainsi, les estimateurs proposés dans cette section ont l’avantage de pouvoir être utilisés pour les domaines non-échantillonnés.

4.2  Modèles linéaires mixtes au niveau unité pour des données fonctionnelles

Les modèles linéaires mixtes au niveau unité proposés par Battese, Harter et Fuller (1988) sont très utilisés dans le cadre de l’estimation de totaux de variables réelles pour des domaines. En effet, il permettent, comme nous le verrons plus en détail par la suite, de traduire à la fois l’effet de l’information auxiliaire sur la variable d’intérêt (par les effets fixes), et les spécificités des domaines (par les effets aléatoires).

Dans cette partie, nous cherchons donc à adapter ces modèles au contexte des données fonctionnelles. Pour cela, nous allons projeter les courbes dans un espace de dimension finie et transformer de cette façcon notre problème fonctionnel en plusieurs problèmes d’estimation de totaux ou de moyennes de variables réelles non-corrélées pour des petits domaines, que nous résoudrons ensuite par des méthodes usuelles. L’utilisation de bases de projection permet donc de préserver la structure de corrélation temporelle de nos données tout en se ramenant à plusieurs sous-problèmes décorrélés d’estimation de moyennes de variables réelles que l’on traite indépendamment par les méthodes usuelles.

4.2.1  Estimation de courbes moyennes par modèles linéaires mixtes au niveau unité appliqués aux scores de l’ACP

Tout comme l’EM en dimension finie, l’EM fonctionnelle est une méthode de réduction de dimension permettant de résumer l’information contenue dans un jeu de données. Elle a été proposée par Deville (1974), ses propriétés théoriques ont été étudiées dans Dauxois, Pousse et Romain (1982) ou Hall, Müller et Wang (2006) et enfin elle a été adaptée dans le cadre des sondages par Cardot et coll. (2010).

Plus formellement, les courbes $Y_i\mathrm{<mo>=</mo>}{\left(Y_i\left(t\right)\right)}_{t\in \left[\mathrm{0,}T\right]}$ sont des fonctions du temps $t$ et nous supposons qu’elles appartiennent à $L^2\left[\mathrm{0,}T\right]\mathrm{,}$ l’espace des fonctions de carré intégrable sur l’intervalle $\left[\mathrm{0,}T\right]\mathrm{.}$ Cet espace est équipé avec le produit scalaire usuel $\mathrm{<mo><</mo>}f\mathrm{,}g\mathrm{<mo>></mo>\; <mo>=</mo>}{\displaystyle {\int }_0^T}f\left(t\right)g\left(t\right)dt$ et la norme $\Vert f\Vert \mathrm{<mo>=</mo>}{\left({\displaystyle {\int }_0^T}{f}^2\left(t\right)dt\right)}^{1/2}\text{}\mathrm{.}$ Soit la fonction de variance covariance $v\left(s\mathrm{,}t\right)$ définie par

$$v\left(s\mathrm{,}t\right)\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^N}\left(Y_i\left(s\right)-\mu \left(s\right)\right)\left(Y_i\left(t\right)-\mu \left(t\right)\right)\mathrm{,}s\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}(4.5)$$

avec $\mu \mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^N}{Y}_i/N$ la courbe moyenne de $Y$ sur la population $U\mathrm{.}$

Soient ${\left({\lambda }_k\right)}_{k\mathrm{<mo>=</mo>}1}^N$ les valeurs propres de $v$ avec ${\lambda }_1\ge {\lambda }_2\dots \ge {\lambda }_N\ge 0$ et ${\left({\xi }_k\right)}_{k\mathrm{<mo>=</mo>}1}^N$ les vecteurs propres orthonormés associés, $v\left(s\mathrm{,}t\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^N}{\lambda }_k{\xi }_k\left(s\right){\xi }_k\left(t\right)\mathrm{.}$

La meilleure approximation de $Y$ dans un espace de dimension $K$ plus petite que $N$ est donnée par la projection de $Y$ dans l’espace engendré par les premiers vecteurs propres ${\xi }_k\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}q$ (Ramsay et Silverman, 2005):

$$Y_i\left(t\right)\mathrm{<mo>=</mo>}\mu \left(t\right)+{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}{f}_{ik}{\xi }_k\left(t\right)+R_i\left(t\right)\mathrm{,}i\in U\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}(4.6)$$

où $f_{ik}$ est la projection (où le score) de $Y_i$ sur la composante ${\xi }_k$ et $R_i\left(t\right)$ le reste représentant l’écart entre la courbe $i$ et sa projection. Le score $f_{ik}$ est indépendant du domaine et peut être calculé comme le produit scalaire entre ${\xi }_k$ et $Y_i-\mu \mathrm{,}$ $f_{ik}\mathrm{<mo>=</mo>\; <mo><</mo>}{Y}_i-\mu \mathrm{,}{\xi }_k\mathrm{<mo>></mo>\; <mo>=</mo>}{\displaystyle {\int }_0^T}\left(Y_i-\mu \right)\left(t\right){\xi }_k\left(t\right)dt\mathrm{.}$ La décomposition donnée dans (4.6) est connue aussi sous le nom de Karhunen-Loève.

En utilisant (4.6), la moyenne ${\mu }_d$ sur le domaine $d$ peut être approximée par

$${\mu }_d\left(t\right)\simeq \mu \left(t\right)+{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}\left(\frac{1}{N_d}{\displaystyle \sum _{i\in U_d}}{f}_{ik}\right){\xi }_k\left(t\right)\mathrm{,}d\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}D\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{.}(4.7)$$

La moyenne inconnue $\mu$ est estimée par l’estimateur d’Horvitz-Thompson

$$\widehat{\mu }\left(t\right)\mathrm{<mo>=</mo>}\frac{1}{N}{\displaystyle \sum _{i\in s}}{d}_i{Y}_i\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right](4.8)$$

et les ${\xi }_k\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K$ sont estimés par ${\widehat{\xi }}_k\mathrm{,}$ les vecteurs propres de la fonction de variance-covariance estimée $\widehat{v}\left(s\mathrm{,}t\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in s}}{d}_i\left(Y_i\left(s\right)-\widehat{\mu }\left(s\right)\right)\left(Y_i\left(t\right)-\widehat{\mu }\left(t\right)\right)/N$ (Cardot et coll., 2010).

Donc, afin d’estimer ${\mu }_d\mathrm{,}$ il nous faut estimer la moyenne des scores sur les composantes principales pour le domaine $d\mathrm{,}$ i.e., ${\overline{f}}_{dk}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in U_d}}{f}_{ik}/N_d\text{}\mathrm{.}$ Pour cela, nous considérons pour chaque composante $f_{ik}\text{}\mathrm{,}$ $k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K$ un modèle linéaire mixte au niveau unité aussi appelé nested error regression model (Battese et coll., 1988) :

$$f_{ik}\mathrm{<mo>=</mo>}{\beta }_k^{\prime }{X}_i+{\nu }_{dk}+{ϵ}_{ik}\mathrm{,}i\in U_d\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K\mathrm{,}(4.9)$$

avec ${\beta }_k^{\prime }{X}_i$ l’effet fixe des informations auxiliaires, ${\nu }_{dk}$ l’effet aléatoire du domaine $d$ et ${ϵ}_{ik}$ le résidu de l’unité $i.$ On suppose que les effets aléatoires des domaines sont indépendants, et suivent une loi commune de moyenne 0 et de variance ${\sigma }_{\nu k}^2.$ Les résidus sont également indépendants, distribués selon une loi de moyenne 0 et de variance ${\sigma }_{ϵk}^2\mathrm{.}$ En outre, les effets aléatoires et les résidus sont également supposés indépendants. Le paramètre $\beta$ du modèle peut être estimé par ${\tilde{\beta }}_k\text{}\mathrm{,}$ l’estimateur «Best Linear Unbiased Estimator» (BLUP) (Rao et Molina, 2015, chapitre 4.7) et l’estimateur BLUP de ${\overline{f}}_{dk}$ s’écrit alors comme un estimateur composite (voir Rao et Molina, 2015):

$${\tilde{\overline{f}}}_{dk}\mathrm{<mo>=</mo>}{\gamma }_k\left({\overline{f}}_{dk\text{}\mathrm{,}s}-{\left({\overline{X}}_{d\text{}\mathrm{,}s}-{\overline{X}}_d\right)}^{\prime }{\tilde{\beta }}_k\right)+\left(1-{\gamma }_k\right){\overline{X}}_d^{\prime }{\tilde{\beta }}_k\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots ,K(4.10)$$

avec ${\gamma }_k\mathrm{<mo>=</mo>}{\sigma }_{\nu k}^2/\left({\sigma }_{\nu k}^2+{\sigma }_{ϵk}^2/n_d\right)$ et ${\overline{X}}_{d\text{}\mathrm{,}s}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in s_d}}{X}_i/n_d\text{}\mathrm{,}$ ${\overline{f}}_{dk\text{}\mathrm{,}s}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in s_d}}{f}_{ik}/n_d$ les moyennes respectives des vecteurs $X_i$ et des scores ${\widehat{f}}_{ik}$ sur $s_d.$ Finalement, la moyenne ${\mu }_d$ est estimée par

$${\widehat{\mu }}_d^{\text{BHF}}\left(t\right)\mathrm{<mo>=</mo>}\widehat{\mu }\left(t\right)+{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}{\widehat{\overline{f}}}_{dk}{\widehat{\xi }}_k\left(t\right)\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{,}(4.11)$$

avec $\widehat{\mu }$ et ${\widehat{\xi }}_k$ les estimations du $\mu$ et de la $k^{\text{e}}$ composante principale ${\xi }_k$ données auparavant.

Les variances ${\sigma }_{\nu k}^2$ et ${\sigma }_{ϵk}^2$ pour $k\mathrm{<mo>=</mo>\; 1,}\dots ,K$ sont inconnues et elles sont estimées par ${\widehat{\sigma }}_{\nu k}^2$ et ${\widehat{\sigma }}_{ϵk}^2$ obtenues par exemple par maximum de vraisemblance restreint (Rao et Molina, 2015). L’estimateur du ${\overline{f}}_{dk}$ est obtenu en remplaçant dans (4.10) ${\gamma }_k$ avec ${\widehat{\gamma }}_k\mathrm{<mo>=</mo>}{\widehat{\sigma }}_{\nu k}^2/\left({\widehat{\sigma }}_{\nu k}^2+{\widehat{\sigma }}_{ϵk}^2/n_d\right)$ et appelé EBLUP (Empirical Best Linear Unbiased Prediction). Néanmoins, le calcul de la fonction de variance (sous le modèle) de ${\widehat{\mu }}_d^{\text{BHF}}\left(t\right)$ est plus compliqué dans ce cas à cause des estimateurs ${\widehat{\xi }}_k$ de composantes principales ${\xi }_k$ et elle sera traitée ailleurs.

Nous remarquons qu’un modèle plus simple, sans les effets aléatoires, aurait pu être considéré pour les scores $f_{ik}\mathrm{<mo>:</mo>}$

$$f_{ik}\mathrm{<mo>=</mo>}{\beta }_k^{\prime }{X}_i+{ϵ}_{ik}\mathrm{,}i\in U\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K\mathrm{,}(4.12)$$

avec ${ϵ}_{ik}$ un résidu de moyenne nulle et de variance ${\sigma }_k^2.$ Dans ce cas, le paramètre ${\beta }_k$ est estimé par ${\widehat{\beta }}_k\mathrm{,}$ l’estimateur de type BLUP et le score moyen sur le domaine $d$ par ${\widehat{\overline{f}}}_{dk}\mathrm{<mo>=</mo>}{\widehat{\beta }}_k^{\prime }{\overline{X}}_d\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots ,K\mathrm{.}$

Si le taux $n_d/N_d$ n’est pas négligeable, alors ${\widehat{\mu }}_d$ est obtenu en suivant le procédé décrit dans la section 4.1.

Remarque 1. Ici, l’EM n’est pas utilisée en tant que méthode de réduction de dimension mais dans le but de décomposer notre problème en plusieurs sous-problèmes non corrélés d’estimation de totaux de variables réelles, que l’on sait résoudre. On garde donc un nombre $K$ de composantes principales aussi élevé que possible, c’est-à-dire égal au minimum du nombre d’instants de discrétisation et du nombre d’individus de l’échantillon.

Remarque 2. Lorsque certaines des variables explicatives du vecteur $X_i$ sont catégorielles, notre méthode, définie dans le cas de variables réelles, doit être adaptée : pour cela, on propose de transformer chaque variable catégorielle en un ensemble d’indicatrices par codage $0-1$ « one hot encoding ». De plus, lorsque le nombre de variables explicatives $p$ est grand, il pourrait également être pertinent d’introduire des pénalisations, de type RIDGE par exemple, dans le problème de régression.

Remarque 3. D’autres bases de projection peuvent être considérées, comme par exemple les ondelettes (voir Mallat, 1999), celles-ci étant particulièrement adaptées aux courbes irrégulières. Une autre solution consisterait enfin à appliquer les modèles linéaires mixtes fonctionnels sur les valeurs des courbes aux instants de discrétisation; néanmoins cette façon de faire ne permettrait pas de prendre en compte les corrélations temporelles de la problématique contrairement aux précédentes.

4.2.2  Estimation de variance par booststrap paramétrique

Pour estimer la précision (variance sous le modèle) des estimateurs de courbes moyennes, on propose de décliner l’approche par bootstrap paramétrique proposée par González-Manteiga, Lombarda, Molina, Morales et Santamara (2008) et reprise ensuite par Molina et Rao (2010). Il s’agit d’une méthode de rééchantillonnage qui consiste à générer un grand nombre $B$ de rééchantillons $s^{*\left(b\right)}\text{}\mathrm{,}b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B$ de taille $n$ par sondage aléatoire simple avec remise dans $s$ et de générer de façon aléatoire les effets aléatoires et fixes dans le modèle de superpopulation estimé. Notons ${\widehat{R}}_i\left(t\right)\mathrm{<mo>=</mo>}{Y}_i\left(t\right)-\widehat{\mu }\left(t\right)+{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>}1}^K}{f}_{k\text{}\mathrm{,}i}{\widehat{\xi }}_k\left(t\right)$ le résidu de projection estimé pour l’unité $i\in s$ (voir aussi formule (4.6)). Pour $b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B\mathrm{,}$ on procède de la manière suivante pour chaque $t\in \left[\mathrm{0,}T\right]:$

1. Générer les effets aléatoires bootstrap de chacun des domaines, pour chaque composante principale :

$${\nu }_{k\text{}\mathrm{,}d}^{*\left(b\right)}\sim \mathcal{N}\left(\mathrm{0,}{\widehat{\sigma }}_{k\text{}\mathrm{,}\nu }^2\right)\mathrm{,}d\mathrm{<mo>=</mo>\; 1,}\dots ,D\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K$$

1. et générer indépendamment de ces effets aléatoires les erreurs individuelles bootstrap pour chaque unité $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}n$ et pour chaque composante principale :

$${ϵ}_{k\text{}\mathrm{,}i}^{*\left(b\right)}\sim \mathcal{N}\left(\mathrm{0,}{\widehat{\sigma }}_{k\text{}\mathrm{,}ϵ}^2\right)\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}n\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K\mathrm{.}$$

2. Calculer les $n$ réplications bootstrap des résidus de projection ${\widehat{R}}_i^{*\left(b\right)}\left(t\right)\mathrm{,}$ pour $i\in s^{*\left(b\right)}$ (cela revient à sélectionner avec remise $n$ résidus de projection parmi les $n$ résidus ${\widehat{R}}_i(t)).$
3. Calculer les réplications bootstrap $Y_i^{*\left(b\right)}\left(t\right)$ conditionnellement aux variables explicatives $X_i$ en suivant le modèle estimé :

$$Y_i^{*\left(b\right)}\left(t\right)=\widehat{\mu }\left(t\right)+{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^K}\underset{f_{k\text{}\mathrm{,}i}^{*\left(b\right)}}{\underbrace{\left(X_i^{\prime }{\widehat{\beta }}_k+{\nu }_{k\text{}\mathrm{,}d}^{*\left(b\right)}+{ϵ}_{k\text{}\mathrm{,}i}^{*\left(b\right)}\right)}}{\widehat{\xi }}_k\left(t\right)+{\widehat{R}}_i^{*\left(b\right)}\left(t\right)\mathrm{,}\forall i\in s_d^{*\left(b\right)}\mathrm{<mo>=</mo>}{s}^{*\left(b\right)}\cap s_d\mathrm{.}$$

1. On remarque que $f_{k\text{}\mathrm{,}i}^{*\left(b\right)}\text{},$ le score simulé pour l’unité $i,$ est obtenu selon la même approche que dans González-Manteiga et coll. (2008).

4. Calculer pour chaque domaine $d,$ la réplication bootstrap ${\widehat{\mu }}_d^{*\left(b\right)}$ sur ce rééchantillon $s^{*\left(b\right)}$ en déclinant l’ensemble de la démarche : EM puis estimation de modèles linéaires mixtes sur les composantes principales par EBLUP.
5. Estimer la variance de l’estimateur ${\widehat{\mu }}_d\left(t\right)$ par la variance empirique des $B$ réplications ${\widehat{\mu }}_d^{*\left(b\right)}\mathrm{<mo>:</mo>}$

$$\widehat{V}\left({\widehat{\mu }}_d\left(t\right)\right)\mathrm{<mo>=</mo>}\frac{1}{B-1}{\displaystyle \sum _{b\mathrm{<mo>=</mo>\; 1}}^B}{\left({\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)-\frac{1}{B}{\displaystyle \sum _{b\mathrm{<mo>=</mo>\; 1}}^B}{\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)\right)}^2\text{}\mathrm{.}$$

Cette approche sera également celle suivie pour approximer la variance de la régression linéaire fonctionnelle (en omettant l’étape 1 de génération des effets aléatoires ${\nu }_{k\text{}\mathrm{,}d}^{*}\text{}).$

4.3  Estimation non-paramétrique par arbres de régression et forêts aléatoires sur des petits-domaines de courbes

Pour obtenir les prédictions individuelles ${\widehat{Y}}_i\left(t\right)$ nous utilisons dans cette section des modèles non paramétriques, des arbres de régression adaptés aux données fonctionnelles et des forêts aléatoires, qui n’imposent plus une forme linéaire à la relation entre informations auxiliaires et variable d’intérêt et permettent plus de souplesse dans la modélisation. En effet, les arbres de régression pour données fonctionnelles sont fréquemment utilisés à EDF et sont connus pour donner des résultats satisfaisants sur les courbes de consommation électrique. Par ailleurs, dans la littérature, les arbres de régression ont été adaptés au cadre des sondages par Toth et Eltinge (2011) mais pas dans une optique d’estimation de totaux sur des petits domaines.

Dans cette section et la suivante, on cherche donc à estimer un cas particulier du modèle général (4.1) dans lequel la fonction $f$ ne dépend pas du domaine auquel appartient l’unité $i,$

$$Y_i\left(t\right)\mathrm{<mo>=</mo>}f\left(X_i\mathrm{,}t\right)+{ϵ}_i\left(t\right)\forall i\in U\mathrm{,}t\in \left[\mathrm{0,}T\right]\mathrm{.}(4.13)$$

4.3.1  Arbres de régression pour des données fonctionnelles

L’arbre de régression et de classification (CART) proposé par Breiman, Friedman, Stone et Olshen (1984) est une technique de statistique non paramétrique très populaire. Son objectif est de prédire la valeur d’une variable cible $Y$ réelle ou catégorielle en fonction d’un vecteur des variables explicatives $X\mathrm{<mo>=</mo>}\left(X_1\mathrm{,}\dots \mathrm{,}{X}_j\mathrm{,}\dots ,X_p\right)$ réelles ou catégorielles. Pour cela, on détermine un partitionnement de l’espace des $X$ en séparant en deux itérativement le jeu de données, selon une règle de décision (critère de « split ») impliquant une unique variable explicative. Ainsi, notre échantillon $s$ constitue le premier nœud $\lambda$ d’un arbre (sa « racine ») que l’on cherche à subdiviser en deux nœuds disjoints ${\lambda }_l$ et ${\lambda }_r$ de façon à ce que les valeurs de la variable cible $Y_i$ réelle soit les plus homogènes possible dans chacun des nœuds. Le critère d’inertie $\kappa \left(\lambda \right)$ utilisé pour quantifier l’homogénéité d’un nœud est fréquemment la somme des carrés des résidus entre les valeurs de $Y_i$ pour les unités $i$ du nœud $\lambda$ et la moyenne de ces valeurs dans le nœud : $\kappa \left(\lambda \right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in \lambda }}{\left(Y_i-{\overline{Y}}_{\lambda }\right)}^2$ où ${\overline{Y}}_{\lambda }$ est la moyenne des $Y_i$ dans le nœud $\lambda .$

Pour les variables $X_j$ quantitatives, les règles de décision sont de la forme

$$\{\begin{array}{lll}i\in {\lambda }_l\hfill & \text{si}\hfill & X_{ji}\mathrm{<mo><</mo>}c\hfill \\ i\in {\lambda }_r\hfill & \text{sinon}\mathrm{,}\hfill & \hfill \end{array}(4.14)$$

avec $c$ un point de coupure à optimiser parmi l’ensemble des valeurs possibles de $X_j.$ Pour les variables qualitatives, elles consistent en un découpage en deux sous-ensembles disjoints de modalités. La recherche du critère de split optimal revient à résoudre le problème d’optimisation

$$\underset{{\lambda }_l\text{}\mathrm{,}{\lambda }_r}{\text{arg}\text{max}}\left(\kappa \left(\lambda \right)-\kappa \left({\lambda }_l\right)-\kappa \left({\lambda }_r\right)\right)\mathrm{.}(4.15)$$

Chacun de ces nœuds sera ensuite à son tour subdivisé en deux nœuds fils et le processus de partitionnement se poursuit jusqu’à atteindre une taille minimale de nœud, jusqu’à ce que la valeur de la variable cible soit la même pour l’ensemble des unités du nœud, ou encore jusqu’à atteindre une profondeur maximale donnée. La partition finale de l’espace est alors constituée par les nœuds finaux de l’arbre, aussi appelés des feuilles. Un résumé de chacune de ces feuilles (très souvent la moyenne pour une variable cible quantitative) devient alors la variable prédite pour l’ensemble des unités affectées à la feuille. Les différents paramètres (taille minimale de nœud et profondeur) peuvent être choisis par validation croisée.

Lorsque la variable $Y$ à prédire n’est plus une variable réelle mais un vecteur de dimension $m\mathrm{<mo>></mo>\; 1},$ le principe de l’arbre de régression s’étend très naturellement : l’algorithme de construction de l’arbre et de choix des paramètres par validation croisée reste inchangé mais le critère d’inertie est modifié. Ainsi, le problème de minimisation s’écrit toujours sous la forme (4.15) mais cette fois le critère est de la forme $\kappa \left(\lambda \right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in \lambda }}{\Vert Y_i-{\overline{Y}}_{\lambda }\Vert }^2$ où $\Vert \cdot \Vert$ est par exemple la norme euclidienne ou la norme associé à la distance de Mahalanobis. Les arbres de régression multivariés ont été utilisés par exemple par De’Ath (2002) dans le cadre d’une application à l’écologie.

Enfin, lorsque la variable à prédire $Y$ est une courbe, l’algorithme de construction de l’arbre et de choix des paramètres est identique mais cette fois, on doit utiliser un critère d’inertie $\kappa$ fonctionnel. De nombreux choix sont possibles. Nous avons choisi de suivre l’approche dite du « Courbotree », décrite dans Stéphan et Cogordan (2009) et fréquemment employée à EDF pour construire des segmentations de jeux de données de courbes de consommation électrique en fonction de variables explicatives. Dans cette approche, on applique la méthode présentée dans le paragraphe précédent pour $Y$ multivariée sur les vecteurs $Y_i\mathrm{<mo>=</mo>}\left(Y_i\left(t_1\right)\mathrm{,}\dots \mathrm{,}{Y}_i\left(t_L\right)\right)$ des valeurs des courbes aux instants de discrétisation, avec la distance euclidienne. La distance euclidienne sur les instants de discrétisation peut alors être vue comme une approximation de la norme $L^2\left[\mathrm{0,}T\right].$ Plus formellement, le critère fonctionnel s’écrit alors

$$\kappa \left(\lambda \right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\in \lambda }}{\displaystyle \sum _{l\mathrm{<mo>=</mo>\; 1}}^L}{\left(Y_i\left(t_l\right)-{\overline{Y}}_{\lambda }\left(t_l\right)\right)}^2\mathrm{,}(4.16)$$

avec ${\overline{Y}}_{\lambda }\left(t_l\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in \lambda }}{Y}_i\left(t_l\right)/n_{\lambda }$ où $n_{\lambda }$ est le nombre d’unités de l’échantillon appartenant au nœud $\lambda .$

En pratique, lorsque l’on travaille sur des données de consommation électrique, les courbes considérées ont des niveaux extrêmement hétérogènes, et l’algorithme du Courbotree basé sur la distance euclidienne peut mal fonctionner lorsqu’il est appliqué sur les données brutes. Fréquemment, on n’utilise donc l’algorithme Courbotree que sur les formes des courbes, c’est-à-dire sur les courbes normalisées ${\tilde{Y}}_i\left(t\right)\mathrm{<mo>=</mo>}{Y}_i\left(t\right)/{\overline{Y}}_i$ où ${\overline{Y}}_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{\mathcal{l}\mathrm{<mo>=</mo>\; 1}}^L}{Y}_i\left(t_{\mathcal{l}}\right)/L$ est la moyenne de $Y_i\left(t\right)$ (ou le niveau) sur tous les instants de mesure $t_1\mathrm{,}\dots \mathrm{,}{t}_L$ (méthode appelée aussi courbotree normalisée). Ensuite, on calcule la prédiction ${\overline{Y}}_i$ à l’aide d’une régression linéaire et on obtient finalement la prédiction de $Y_i\left(t\right)$ en faisant le produit entre la prédiction de ${\tilde{Y}}_i\left(t\right)$ et celle de ${\overline{Y}}_i\mathrm{.}$

4.3.2  Estimation de la variance

Afin d’estimer la variance sous le modèle de nos estimateurs de courbes moyennes par domaine, on va suivre une approche par bootstrap très similaire à celle proposée pour les modèles paramétriques. Ici, notre modèle de superpopulation s’écrit sous la forme

$$Y_i\left(t\right)\mathrm{<mo>=</mo>}f\left(X_i\mathrm{,}t\right)+{ϵ}_i\left(t\right)\mathrm{,}\forall i\in U\mathrm{.}(4.17)$$

Soient $\widehat{f}\left(X_i\mathrm{,}t\right)\mathrm{,}$ pour tout $i\in s$ la valeur prédite pour l’unité $i$ par arbre de régression, et ${\widehat{ϵ}}_i\left(t\right)\mathrm{<mo>=</mo>}{Y}_i\left(t\right)-\widehat{f}\left(X_i\mathrm{,}t\right)\mathrm{,}$ pour tout $i\in s$ le résidu estimé pour cette unité. L’idée de notre méthode d’approximation de précision est, comme dans le cas des modèles linéaires mixtes, de générer un grand nombre $B$ de rééchantillons $s^{*\left(b\right)}\text{}\mathrm{,}b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B$ de taille $n$ par sondage aléatoire simple avec remise dans $s,$ puis de calculer l’estimateur de la courbe moyenne par domaine sur chaque rééchantillon et enfin de déduire la variance de l’estimateur par la variabilité entre rééchantillons. La méthode de bootstrap employée ici est connue aussi sous le nom de bootstrap des résidus dans le cas du modèle linéaire.

Plus précisément, pour $b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B\mathrm{,}$ on procède de la manière suivante pour chaque $t\in \left[\mathrm{0,}T\right]:$

1. Calculer les réplications bootstrap des résidus estimés ${\widehat{ϵ}}_i^{*\left(b\right)}\left(t\right)$ pour $i\in s^{*\left(b\right)}\text{}.$
2. Calculer les réplications bootstrap de $Y_i\left(t\right):$

$$Y_i^{*\left(b\right)}\left(t\right)\mathrm{<mo>=</mo>}\widehat{f}\left(X_i\mathrm{,}t\right)+{\widehat{ϵ}}_i^{*\left(b\right)}\left(t\right)\mathrm{,}\forall i\in s^{*\left(b\right)}$$

1. et recalculer, pour chaque domaine $d,$ l’estimateur de moyenne ${\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)$ sur ce rééchantillon.

3. Estimer la variance par la variance empirique des $B$ réplications bootstrap ${\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)\mathrm{,}$ $\widehat{V}\left({\widehat{\mu }}_d\left(t\right)\right)\mathrm{<mo>=</mo>}\frac{1}{B-1}{\displaystyle {\sum }_{b\mathrm{<mo>=</mo>\; 1}}^B}{\left({\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)-\frac{1}{B}{\displaystyle {\sum }_{b\mathrm{<mo>=</mo>\; 1}}^B}{\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)\right)}^2\text{}\mathrm{.}$

La démarche sera identique si on estime la fonction $f$ par des forêts aléatoires plutôt que des arbres de régression.

4.4  Agrégation de prédictions par forêts aléatoires pour des courbes

La littérature met souvent en évidence les médiocres performances prédictives des arbres de régression en comparaison d’autres techniques telles que les SVM (voir par exemple Cristianini et Shawe-Taylor, 2000). En effet, les arbres de régression peuvent être instables et très dépendants de l’échantillon sur lequel ils ont été construits. Pour remédier à ce défaut, Breiman (2001) a proposé l’algorithme des forêts aléatoires. Il s’agit d’une technique ensembliste qui, comme son nom l’indique, consiste à agréger les prédictions issues de différents arbres de régression. Le fait que l’agrégation de prédicteurs instables induise une réduction de variance a été montré notamment dans Breiman (1998). Pour une variable cible quantitative, l’agrégation des prédictions est réalisée en prenant la moyenne des prédictions de chacun des arbres.

Afin de diminuer la variance de la prédiction agrégée, l’objectif est de construire des arbres très différents les uns des autres. L’algorithme de Breiman introduit de la variabilité dans la construction des arbres d’une part en réalisant un rééchantillonnage (tirage aléatoire simple avec remise) des unités et d’autre part en sélectionnant aléatoirement, pour chaque « split » de l’arbre, un sous-ensemble de variables explicatives candidates. Par rapport à un arbre de régression, il y a donc deux paramètres supplémentaires à ajuster pour une forêt aléatoire : le nombre d’arbres et le nombre de variables explicatives candidates à chaque split.

Lorsque la variable d’intérêt est multivariée (ou fonctionnelle), l’algorithme proposé par Breiman s’adapte aisément, en agrégeant les arbres de régression multivariés (ou fonctionnels) présentés dans le paragraphe précédent. Les forêts aléatoires multivariées ont par exemple été étudiées par Segal et Xiao (2011).

L’algorithme que nous proposons ici, appelé « CourboForest », consiste simplement à agréger des arbres de régression fonctionnels construits selon l’approche « Courbotree », c’est-à-dire des arbres de régression multivariés appliqués sur les vecteurs $\left(Y_i\mathrm{<mo>=</mo>}{Y}_i\left(t_1\right)\mathrm{,}\dots \mathrm{,}{Y}_i\left(t_L\right)\right)$ des valeurs des courbes aux instants de discrétisation, avec pour critère de split l’inertie basée sur la distance euclidienne définie par l’équation (4.16).

• Précédent
• Page principale de l'article
• Suivant