---
title: "5. Linear regression of histogram variables"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{5. Linear regression of histogram variables}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 5, dpi = 72, out.width = "100%")
```

```{r setup}
library(HistDAWass)
```

## The two-component model

HistDAWass implements the linear regression of a histogram variable $Y$ on
histogram variables $X_1, \dots, X_p$ proposed by Irpino and Verde (2015). Each
quantile function is decomposed into its mean and its centred quantile
function, $Q(p) = \bar{y} + Q^c(p)$, and the model is
$$Q_{y_i}(p) = \beta_0 + \sum_{j} \beta_j \bar{x}_{ij} +
\sum_{j} \gamma_j Q^c_{x_{ij}}(p) + \varepsilon_i(p).$$
The $\beta$ parameters (prefix `AV_`) describe the relation between the
means, as in an ordinary regression; the $\gamma$ parameters (prefix `CEN_`)
describe how the variability and the shape of the explanatory histograms
transfer to the response. The $\gamma$ are estimated by non-negative least
squares, so that the predictions are always valid quantile functions. The
parameters minimise the sum of the squared L2 Wasserstein distances between
observed and predicted histograms.

## Fitting the model

The `Agronomique` data describe 22 French regions by the histograms of the
total costs of the farms (`Y_TSC`) and of the costs for wheat, pigs and cow
milk.

```{r}
get.MatH.main.info(Agronomique)$varnames
pars <- WH.regression.two.components(Agronomique, Yvar = 1, Xvars = 2:4)
pars
```

## Predictions and goodness of fit

The predicted histograms are obtained from the explanatory variables (in the
same order used for the estimation):

```{r}
pred <- WH.regression.two.components.predict(Agronomique[, 2:4], parameters = pars)
pred[1:3, ]
```

`WH.regression.GOF()` returns the root mean squared error based on the L2
Wasserstein distance (`RMSE_W`), the `OMEGA` index and a pseudo-$R^2$ (both
between 0 and 1, the higher the better):

```{r}
gof <- WH.regression.GOF(observed = Agronomique[, 1], predicted = pred)
gof$RMSE_W
gof$OMEGA
gof$PSEUDOR2$index
```

## Plotting observed and predicted histograms

```{r, fig.height = 8}
plotPredVsObs(pred[1:8, ], Agronomique[1:8, 1], type = "HISTO", ncolu = 2)
```

```{r, fig.height = 4}
plotPredVsObs(pred[1:4, ], Agronomique[1:4, 1], type = "CDF", ncolu = 2)
```

`plot_errors()` shows, for each individual, the difference between the
observed and the predicted quantile functions (or densities), and prints the
RMSE:

```{r, fig.height = 4}
plot_errors(pred, Agronomique[, 1], type = "HISTO_QUA")
```

## Simplified computations

With many individuals or histograms with many bins, `simplify = TRUE`
approximates the histograms with `qua` quantiles:

```{r}
pars20 <- WH.regression.two.components(Agronomique, Yvar = 1, Xvars = 2:4,
                                       simplify = TRUE, qua = 20)
round(rbind(exact = pars, simplified = pars20), 4)
```

## Reference

Irpino, A., Verde, R. (2015) Linear regression for numeric symbolic variables:
a least squares approach based on Wasserstein Distance. *Advances in Data
Analysis and Classification*, 9, 81-106.
<https://doi.org/10.1007/s11634-015-0197-7>