5. Linear regression of histogram variables

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.

get.MatH.main.info(Agronomique)$varnames
#> [1] "Y_TSC"   "X_Wheat" "X_Pig"   "X_Cmilk"
pars <- WH.regression.two.components(Agronomique, Yvar = 1, Xvars = 2:4)
pars
#> (AV_Intercept)     AV_X_Wheat       AV_X_Pig     AV_X_Cmilk    CEN_X_Wheat 
#>   1.564758e+04   7.207630e-01   7.356636e-01   2.791063e-01   8.342455e-01 
#>      CEN_X_Pig    CEN_X_Cmilk 
#>   5.866346e-01   2.996154e-01

Predictions and goodness of fit

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

pred <- WH.regression.two.components.predict(Agronomique[, 2:4], parameters = pars)
pred[1:3, ]
#> a matrix of distributions 
#>  1  variables  3  rows 
#>  each distibution in the cell is represented by the mean and the standard deviation 
#>                       Predicted       
#>  Ile-de-France [m= 44988  ,s= 17713 ]
#>  Champagne-Ardenne [m= 28290  ,s= 15313 ]
#>  Picardie [m= 47427  ,s= 24750 ]

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

gof <- WH.regression.GOF(observed = Agronomique[, 1], predicted = pred)
gof$RMSE_W
#> [1] 7056.517
gof$OMEGA
#> [1] 0.8617325
gof$PSEUDOR2$index
#> [1] 0.6784146

Plotting observed and predicted histograms

plotPredVsObs(pred[1:8, ], Agronomique[1:8, 1], type = "HISTO", ncolu = 2)

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:

plot_errors(pred, Agronomique[, 1], type = "HISTO_QUA")

#> [1] 7056.517

Simplified computations

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

pars20 <- WH.regression.two.components(Agronomique, Yvar = 1, Xvars = 2:4,
                                       simplify = TRUE, qua = 20)
round(rbind(exact = pars, simplified = pars20), 4)
#>            (AV_Intercept) AV_X_Wheat AV_X_Pig AV_X_Cmilk CEN_X_Wheat CEN_X_Pig
#> exact            15647.58     0.7208   0.7357     0.2791      0.8342    0.5866
#> simplified       15652.29     0.7203   0.7310     0.2790      0.8317    0.5864
#>            CEN_X_Cmilk
#> exact           0.2996
#> simplified      0.3016

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