Sample-based estimation of mean electricity consumption curves for small domains
Section 4. Model-based estimation methods

• Previous
• Next

In this section, we use the model from Valliant, Dorfman and Royall (2000), in which curves $Y_i$ are considered to be random and we propose four innovative approaches for responding to our problem estimating curves with a mean demand for small domains. Assuming that $Y_i\left(t\right)$ and the auxiliary information vector $X_i$ are available for each individual $i$ in the domain $d$ and that the mean ${\overline{X}}_d={\displaystyle {\sum }_{i\in U_d}}{X}_i/N_d$ is also known.

We assume that the auxiliary variables are related to the demand curves according to a functional model of superpopulation on all of the population that is generally expressed as:

$$\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)$$

with $f_d,$ an unknown regression function to be estimated, which can vary from one domain to another, and ${ϵ}_i$ a process of zero expectation noise, zero covariance for different individuals and non-null regarding time.

If the size of domain $N_d$ is large, then the mean ${\mu }_d$ will be estimated by

$${\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{,}$$

where ${\widehat{Y}}_i\left(t\right)\mathrm{<mo>=</mo>}{\widehat{f}}_d\left(X_i\mathrm{,}t\right)$ is the prediction of $Y_i\left(t\right).$ Otherwise, the mean ${\mu }_d$ is estimated by (see Valliant et al., 2000):

$${\widehat{\mu }}_d\left(t\right)=\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)$$

The quality of our estimates thus depends on the quality of our model: if the model is false, that may lead to biases in the estimates.

4.1  Functional linear model

The simplest model of form (4.1) is the functional linear regression model from 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)$$

where the residuals ${\epsilon }_i\left(t\right)$ are independent for $i\ne j,$ distributed based on a law of means of 0 and of variance of ${\sigma }_i^2\left(t\right)\mathrm{.}$ If the size of domain $N_d$ is large, then the mean of $Y$ in domain $d$ is estimated by

$${\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{,}$$

where ${\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)$ is the best linear unbiased (BLU) estimator of $\beta$ that does not depend on domain $d\mathrm{.}$ Estimator ${\widehat{\mu }}_d^{\text{blu}}\left(t\right)$ can be expressed as a weighted sum of $Y_i\left(t\right):$

$${\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{,}$$

where the weight $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)$ now dependant on time $t\mathrm{.}$ If $n_d/N_d$ is not insignificant, then the mean ${\mu }_d$ is estimated using (4.2) by:

$${\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{.}$$

This estimator can again be expressed as a weighted sum of $Y_i\left(t\right)$ with weights that will still depend on time $t\mathrm{.}$ The variance function (based on the model) of ${\widehat{\mu }}_d^{\text{blu}}\left(t\right)$ can be derived using Rao and Molina (2015), Chapter 7. The variance function ${\sigma }_i^2\left(t\right)$ is unknown and can be estimated by following Rao and Molina (2015). By replacing ${\sigma }_i^2\left(t\right)$ with ${\widehat{\sigma }}_i^2\left(t\right)\mathrm{,}$ we will obtain the empirical best linear unbiased predictor (EBLUP) of ${\mu }_d$ and its variance can be obtained using the method set out by Rao and Molina (2015). This EBLUP estimator does not use the sample weight $d_i$ and therefore is not consistent in terms of the sampling plan (unless the sample weights are constant for units in the same domain $d).$ A modified estimator, also referred to as a pseudo-EBLUP, can be constructed using the new approach described by Rao and Molina (2015, Chapter 7), equal in this case to the estimator set out in (3.4).

If $Y_i\left(t\right)$ is unknown for the units in domain $d\mathrm{,}$ the following indirect estimator can be used:

$${\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)$$

with $\widehat{\beta }\left(t\right)$ given in (3.5) and the weights ${\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$ are not dependant on time $t,$ unlike $w_{id}^{\text{blu}}\text{}\mathrm{.}$ Thus, the estimators proposed in this section have the benefit of being able to be used for unsampled domains.

4.2  Unit-level linear mixed models for functional data

The unit-level linear mixed models proposed by Battese, Harter and Fuller (1988) are very useful in estimating total actual variables for domains. As we will see later in more detail, they can translate both the effect of auxiliary information on the interest variable (by fixed effects), and the specifics of the domains (by random effects).

In this section, we thus attempt to adapt those models to the context of functional data. To that end, we will project curves in a space of defined dimensions and in that way, transform our functional problem into several problems in estimating total or mean real uncorrelated variables for small domains, which we will then resolve using the usual methods. The use of projection bases thus makes it possible to preserve the temporal correlation structure of our data while arriving at several unrelated subproblems in estimating real variables, which we treat independently using the usual methods.

4.2.1  Estimation of curves using unit-level linear mixed models applied to PCA scores

Like PCA in finite dimensions, functional PCA is a dimension-reduction method that makes it possible to summarize information contained in a data set. It was proposed by Deville (1974), its theoretical properties were studied in Dauxois, Pousse and Romain (1982) or Hall, Mülller and Wang (2006) and, finally, it was adapted to surveys by Cardot et al. (2010).

More formally, the curves $Y_i\mathrm{<mo>=</mo>}{\left(Y_i\left(t\right)\right)}_{t\in \left[\mathrm{0,}T\right]}$ are functions of time $t$ and we assume that they belong to $L^2\left[\mathrm{0,}T\right]\mathrm{,}$ the space of square-integrable functions in the interval $\left[\mathrm{0,}T\right]\mathrm{.}$ That space is equipped with the usual scalar product $\mathrm{<mo><</mo>}f\mathrm{,}g\mathrm{<mo>></mo>\; <mo>=</mo>}{\displaystyle {\int }_0^T}f\left(t\right)g\left(t\right)dt$ and the standard $\Vert f\Vert \mathrm{<mo>=</mo>}{\left({\displaystyle {\int }_0^T}{f}^2(t)dt\right)}^{1/2}\text{}\mathrm{.}$ The variance covariance function $v\left(s\mathrm{,}t\right)$ defined by:

$$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)$$

with $\mu \mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>\; 1}}^N}{Y}_i/N$ the mean curve of $Y$ on the population $U\mathrm{.}$

Let ${\left({\lambda }_k\right)}_{k\mathrm{<mo>=</mo>}1}^N$the eigen values of $v$ with ${\lambda }_1\ge {\lambda }_2\dots \ge {\lambda }_N\ge 0$ and ${\left({\xi }_k\right)}_{k\mathrm{<mo>=</mo>}1}^N$the related orthonormal eigen vectors, $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{.}$

The best approximation of $Y$ in a dimensional space $K$ smaller than $N$ is given by the projection of $Y$ in the space created by the first eigen vectors ${\xi }_k\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}q$ (Ramsay and 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)$$

where $f_{ik}$ is the projection (or score) of $Y_i$ on the component ${\xi }_k$ and $R_i\left(t\right),$ the rest representing the difference between curve $i$ and its projection. The score $f_{ik}$ is independent of the domain and can be calculated as the scalar product between ${\xi }_k$ and $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{.}$ The decomposition given in (4.6) is also known as Karhunen-Loève.

Using (4.6), the mean ${\mu }_d$ on the domain $d$ can be approximated by

$${\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)$$

The unknown mean $\mu$ is estimated using the Horvitz-Thompson estimator

$$\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)$$

and the ${\xi }_k\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K$ are estimated by ${\widehat{\xi }}_k\mathrm{,}$ the eigen vectors in the estimated variance-covariance function $\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 al., 2010).

Thus, to estimate ${\mu }_d\mathrm{,}$ we must estimate the mean scores on the principal components for the domain $d\mathrm{,}$ i.e., ${\overline{f}}_{dk}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in U_d}}{f}_{ik}/N_d\text{}\mathrm{.}$ To that end, for each component $f_{ik}\text{}\mathrm{,}$ $k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}K,$ we consider a unit-level linear mixed model, also known as a nested error regression model (Battese et al., 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)$$

with ${\beta }_k^{\prime }{X}_i$ the fixed effect of auxiliary information, ${\nu }_{dk}$ the random effect of the domain $d$ and ${ϵ}_{ik}$ the residual of unit $i.$ We assume that the random effects of the domains are independent and follow a common law of means of 0 and of variance of ${\sigma }_{\nu k}^2.$ The residuals are also independent, distributed based on a law of means of 0 and of variance of ${\sigma }_{ϵk}^2\mathrm{.}$ In addition, the random effects and residuals are also assumed to be independent. The parameter $\beta$ in the model can be estimated by ${\tilde{\beta }}_k\text{}\mathrm{,}$ the best linear unbiased estimator (BLUP) (Rao and Molina, 2015, Chapter 4.7) and the BLUP estimator ${\overline{f}}_{dk}$ is thus expressed as a composite estimator (see Rao and 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)$$

with ${\gamma }_k\mathrm{<mo>=</mo>}{\sigma }_{\nu k}^2/\left({\sigma }_{\nu k}^2+{\sigma }_{ϵk}^2/n_d\right)$ and ${\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$ the respective means of the vectors $X_i$ and the scores ${\widehat{f}}_{ik}$ on $s_d.$ Finally, the mean ${\mu }_d$ is estimated by

$${\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)$$

with $\widehat{\mu }$ and ${\widehat{\xi }}_k$ the estimates of $\mu$ and $k^{\text{th}}$ principal component ${\xi }_k$ given previously.

The variances ${\sigma }_{\nu k}^2$ and ${\sigma }_{ϵk}^2$ for $k\mathrm{<mo>=</mo>\; 1,}\dots ,K$ are unknown and are estimated by ${\widehat{\sigma }}_{\nu k}^2$ and ${\widehat{\sigma }}_{ϵk}^2$ obtained, for example, by restricted maximum likelihood (Rao and Molina, 2015). The estimator for ${\overline{f}}_{dk}$ is obtained by replacing ${\gamma }_k$ in (4.10) with ${\widehat{\gamma }}_k\mathrm{<mo>=</mo>}{\widehat{\sigma }}_{\nu k}^2/\left({\widehat{\sigma }}_{\nu k}^2+{\widehat{\sigma }}_{ϵk}^2/n_d\right)$ and is known as empirical best linear unbiased prediction (EBLUP). Nonetheless, the calculation of the variance function (based on the model) of ${\widehat{\mu }}_d^{\text{BHF}}\left(t\right)$ is more complicated in this case because of the estimators ${\widehat{\xi }}_k$ of the principal components ${\xi }_k$ and will be examined elsewhere.

We note that a simpler model, without the random effects, could have been considered for the 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)$$

with ${ϵ}_{ik}$ a null mean residual ${\sigma }_k^2.$ In this case, the parameter ${\beta }_k$ is estimated by ${\widehat{\beta }}_k\mathrm{,}$ the BLUP estimator and the mean score on the domain $d$ by ${\widehat{\overline{f}}}_{dk}\mathrm{<mo>=</mo>}{\widehat{\beta }}_k^{\prime }{\overline{X}}_d\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots ,K\mathrm{.}$

If the rate of $n_d/N_d$ is not insignificant, then ${\widehat{\mu }}_d$ is obtained using the procedure described in Section 4.1.

Note 1: Here, the PCA is not used as a dimension-reduction method, but to decompose our problem into several unrelated subproblems in the estimation of total real variables, which we know how to resolve. We thus keep a number $K$ of principle components as high as possible, i.e., equal to the minimum number of instants of discretization and the number of individuals in the sample.

Note 2: When certain explanatory variables in vector $X_i$ are categorical, our method, defined in the case of real variables, must be adapted: to that end, we propose transforming each categorical variable into a series of one hot encoding indicators $0-1.$ As well, when the number of explanatory variables $p$ is large, it may also be relevant to introduce penalties, RIDGE-style for example, in the regression problem.

Note 3: Other projection bases can be considered, such as wavelets (see Mallat, 1999), as they are particularly suited to irregular curves. Finally, another solution would be to apply the functional linear mixed models for curve values to the instants of discretization; however, that method would not allow for consideration of temporal correlations in the problem, unlike the previous ones.

4.2.2  Estimating variance by parametric bootstrap

To estimate the accuracy (variance based on the model) of mean curve estimators, we propose declining the parametric bootstrap method proposed by González-Manteiga, Lombarda, Molina, Morales and Santamara (2008) and then reiterated by Molina and Rao (2010). This is a replicate method that consists of generating a large number $B$ of replicates $s^{*\left(b\right)}\text{}\mathrm{,}b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B$ of size $n$ by simple random sampling with replacement in $s$ and randomly generating the random and fixed effects in the estimated superpopulation model. Note ${\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)$ the estimated projection residual for the unit $i\in s$ (see also the formula in (4.6)). For $b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B\mathrm{,}$ we proceed as follows for each $t\in \left[\mathrm{0,}T\right]:$

1. Generate the random bootstrap effects of each domain, for each principal component:

$${\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. and, independent of those random effects, generate the individual bootstrap errors for each unit $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}n$ and for each principal component:

$${ϵ}_{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. Calculate the $n$ bootstrap replicates of the projection residuals ${\widehat{R}}_i^{*\left(b\right)}\left(t\right)\mathrm{,}$ for $i\in s^{*\left(b\right)}$ (this means selection with replacement of $n$ projection residuals among the $n$ residuals ${\widehat{R}}_i(t)).$
3. Calculate the bootstrap replicates $Y_i^{*\left(b\right)}\left(t\right)$ conditional to the explanatory variables $X_i$ using the estimated model:

$$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. We see that $f_{k\text{}\mathrm{,}i}^{*\left(b\right)}\text{},$ the simulated score for the unit $i,$ is obtained using the same

4. For each domain $d,$ calculate the bootstrap replicate ${\widehat{\mu }}_d^{*\left(b\right)}$ on the replicate $s^{*\left(b\right)}$ by declining the entire process: PCA and estimation of linear mixed models on principal components by means of EBLUP.
5. Estimate the estimator’s variance ${\widehat{\mu }}_d\left(t\right)$ by the empirical variance of the $B$ replicates ${\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{.}$$

This approach will also be the one used to approximate the variance of the functional linear regression (omitting step 1 of generating random effects ${\nu }_{k\text{}\mathrm{,}d}^{*}\text{}).$

4.3  Non-parametric estimation using regression trees and random forests for small curve domains

In this section, to obtain individual predictions ${\widehat{Y}}_i\left(t\right),$ we use non-parametric models, regression trees adapted to functional data, and random forests, which no longer require a linear form in the relation between auxiliary information and interest variable and allow more flexibility in modelling. In fact, regression trees for functional data are frequently used at EDF and are known to give satisfactory results on electricity consumption curves. As well, in literature, regression trees have been adapted to surveys by Toth and Eltinge (2011), but not for estimating totals in small domains.

In this section and the next section, we are therefore seeking to estimate a specific case of the general model (4.1) in which the function $f$ does not depend on the domain of the 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  Regression trees for functional data

The classification and regression tree (CART) proposed by Breiman, Friedman, Stone and Olshen (1984) is a very popular non-parametric statistical technique. Its goal is to predict the value of a real or categorical target variable $Y$ based on a vector of real or categorical explanatory variables $X\mathrm{<mo>=</mo>}\left(X_1\mathrm{,}\dots \mathrm{,}{X}_j\mathrm{,}\dots ,X_p\right).$ To that end, we determine a partitioning of the space of $X$ by repeatedly splitting the data set in two, using the decision rule (split criteria) involving a single explanatory variable. Thus, our sample $s$ is the first node $\lambda$ in a tree (its “root”) that we seek to subdivide into two separate nodes ${\lambda }_l$ and ${\lambda }_r$ such that the values of the real target variable $Y_i$ are as homogenous as possible in each node. The inertia criterion $\kappa \left(\lambda \right)$ used to quantify the homogeneity of a node is frequently the sum of the squares of residuals between the values of $Y_i$ for units $i$ in node $\lambda$ and the mean of those values in the node: $\kappa \left(\lambda \right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in \lambda }}{\left(Y_i-{\overline{Y}}_{\lambda }\right)}^2$ where ${\overline{Y}}_{\lambda }$ is the mean of $Y_i$ in node $\lambda .$

For the variables $X_j$ which are quantitative, the decision rules are expressed as

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

with $c$ a cut-off point to be optimized among all possible values of $X_j.$ For qualitative variables, they consist of dividing into two separate subsets of modalities. The search for the optimal split criterion is a matter of resolving the optimization problem

$$\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)$$

Each of these nodes will then be subdivided in turn into two child nodes and the partitioning process continues until a minimal node size is obtained, until the value of the target variable is the same for all units of the node, or until a given maximum depth is attained. The final partition of the space is then made up of the final nodes in the tree, also called leaves. A summary of each of those leaves (often the mean for a quantitative target variable) then becomes the predicted variable for all units assigned to the leaf. The various parameters (minimum node size and depth) can be selected by cross-validation.

When the variable $Y$ to be predicted is not a real variable, but a dimension vector $m\mathrm{<mo>></mo>\; 1},$ the regression tree principle extends very naturally: the tree construction algorithm and the choice of cross-validation parameters remain unchanged, but the inertia criterion is modified. Thus, the minimization problem is still in the form of (4.15), but this time the criterion is in the form of $\kappa \left(\lambda \right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in \lambda }}{\Vert Y_i-{\overline{Y}}_{\lambda }\Vert }^2$ where $\Vert \cdot \Vert$ is, for example, the Euclidean norm or the Mahalanobis distance norm. Multivariate regression trees were used, for example, by De’Ath (2002) in an ecological application.

Finally, when the variable to be predicted $Y$ is a curve, the algorithm for construction of the tree and for choosing the parameters is the same, but this time a functional inertia criterion $\kappa$ must be used. There are many possible choices. We chose to follow the “Courbotree” approach described in Stéphan and Cogordan (2009) and frequently used at EDF for building segmentations of data sets of electricity consumption curves based on explanatory variables. In that approach, we apply the method presented in the previous paragraph for multivariate $Y$ on vectors $Y_i\mathrm{<mo>=</mo>}\left(Y_i\left(t_1\right)\mathrm{,}\dots \mathrm{,}{Y}_i\left(t_L\right)\right)$ of curve values at the instants of discretization, with the Euclidian distance. The Euclidian distance on instants of discretization can thus be seen as an approximation of the norm $L^2\left[\mathrm{0,}T\right].$ More formally, the functional criterion is thus expressed as

$$\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)$$

with ${\overline{Y}}_{\lambda }\left(t_l\right)\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\in \lambda }}{Y}_i\left(t_l\right)/n_{\lambda }$ where $n_{\lambda }$ is the number of units in the sample that belong to the node $\lambda .$

In practice, when working on electricity consumption data, the curves considered are at extremely similar levels, and the Courbotree algorithm based on the Euclidian distance may not work well when applied to raw data. Often, the Courbotree algorithm is therefore only used on the curve forms, i.e., the normalized curves ${\tilde{Y}}_i\left(t\right)\mathrm{<mo>=</mo>}{Y}_i\left(t\right)/{\overline{Y}}_i$ where ${\overline{Y}}_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{\mathcal{l}\mathrm{<mo>=</mo>\; 1}}^L}{Y}_i\left(t_{\mathcal{l}}\right)/L$ is the mean of $Y_i\left(t\right)$ (or the level) on all measurement instants $t_1\mathrm{,}\dots \mathrm{,}{t}_L$ (method also known as normalized Courbotree). We then calculate the prediction ${\overline{Y}}_i$ using a linear regression and finally obtain the prediction of $Y_i\left(t\right)$ by obtaining the product between the prediction of ${\tilde{Y}}_i\left(t\right)$ and that of ${\overline{Y}}_i\mathrm{.}$

4.3.2  Variance estimation

To estimate the variance under the model of our estimators for mean curves by domain, we will follow a bootstrap approach very similar to the one proposed for parametric models. Here, our superpopulation model is expressed as

$$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)$$

Let $\widehat{f}\left(X_i\mathrm{,}t\right)\mathrm{,}$ for all $i\in s$ the predicted value for the unit $i$ by regression tree, and ${\widehat{ϵ}}_i\left(t\right)\mathrm{<mo>=</mo>}{Y}_i\left(t\right)-\widehat{f}\left(X_i\mathrm{,}t\right)\mathrm{,}$ for all $i\in s$ the estimated residual for that unit. The idea of our accuracy approximation method is, as with linear mixed models, to generate a large number $B$ of replicates $s^{*\left(b\right)}\text{}\mathrm{,}b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B$ of size $n$ by simple random sampling with replacement in $s,$ and calculate the estimator of the mean curve by domain on each replicate and, finally, deduct the variance from the estimator by the variability between replicates. The bootstrap method used here is also known as residual bootstrap in linear model cases.

More specifically, for $b\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}B\mathrm{,}$ we proceed as follows for each $t\in \left[\mathrm{0,}T\right]:$

1. Calculate the bootstrap replications of the estimated residuals ${\widehat{ϵ}}_i^{*\left(b\right)}\left(t\right)$ for $i\in s^{*\left(b\right)}\text{}.$
2. Calculate the bootstrap replications for $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. and recalculate, for each domain $d,$ the mean estimator ${\widehat{\mu }}_d^{*\left(b\right)}\left(t\right)$ on this replicate.

3. Estimate the variance by the empirical variance of the $B$ bootstrap replicates ${\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{.}$

The process is identical if we estimate the function $f$ by random forests rather than regression trees.

4.4  Aggregation of predictions by random forests for curves

The literature often highlights the mediocre predictive performances of regression trees compared to other techniques such as SVMs (see, for example, Cristianini and Shawe-Taylor, 2000). Regression trees can be unstable and highly dependant on the sample on which they were built. To resolve that default, Breiman (2001) proposed the random forest algorithm. This is a set technique that, as its name suggests, consists of aggregating predictions resulting from different regression trees. The fact that the aggregation of unstable predictors leads to a reduction in variance was particularly shown by Breiman (1998). For a quantitative target variable, the aggregation of predictions is performed by taking the mean predictions for each of the trees.

To reduce the variance of the aggregate prediction, the objective is to build trees that are very different from each other. The Breiman algorithm introduces variability in the construction of the trees on the one hand by means of replication (simple random sampling with replacement) of units and, on the other hand, by means of random selection, for each “split” in the tree, of a subset of candidate explanatory variables. For a regression tree, there are therefore two additional parameters to be adjusted for a random forest: the number of trees and the number of candidate explanatory variables in each split.

When the interest variable is multivariate (or functional), the algorithm proposed by Breiman adapts easily, by aggregating the multivariate (or functional) regression trees presented in the previous paragraph. Multivariate random forests have, for example, been studied by Segal and Xiao (2011).

The algorithm that we are proposing here, called “Courboforest,” simply consists of aggregating the functional regression trees constructed using the “Courbotree” approach, i.e., the multivariate regression trees applied to the vectors $\left(Y_i\mathrm{<mo>=</mo>}{Y}_i\left(t_1\right)\mathrm{,}\dots \mathrm{,}{Y}_i\left(t_L\right)\right)$ of the values of the curves at the instants of discretization, with the split criterion being the inertia based on the Euclidean distance defined by equation (4.16).

• Previous
• Article main page
• Next