2
$\begingroup$

This question already has an answer here:

I am aware of a variety of methods for simultaneously predicting multiple outcomes known sometimes as multivariate regression/analysis. However, my situation is a little more special. I am trying to predict a vector of values and where each part ranges from 0-1, and the sum of the vector must be equal to 1. A typical example of this would be population fractions of exclusive groups.

The most simple approach to this I can think of is to use OLS and rescale the predictions so they do not violate data structure. Here's an example using R and US state population data.

> #packages
> suppressPackageStartupMessages(library(tidyverse))

> suppressPackageStartupMessages(library(poliscidata))

> suppressPackageStartupMessages(library(rms))

> suppressPackageStartupMessages(library(magrittr))

> d = states

> #recode
> d$black = d$blkpct10 / 100

> d$hispanic = d$hispanic10 / 100

> d$white = 1 - d$black - d$hispanic

> #test sums
> rowSums(d %>% select(black, hispanic, white))
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 
 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
46 47 48 49 50 
 1  1  1  1  1 

> #regression
> #using 4 variables
> #OLS
> (ols_black = ols(black ~ abortion_rank12 + ba_or_more + cig_tax12 + conserv_advantage, data = d))
Linear Regression Model

 ols(formula = black ~ abortion_rank12 + ba_or_more + cig_tax12 + 
     conserv_advantage, data = d)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
 Obs      50    LR chi2      4.81    R2       0.092    
 sigma0.0958    d.f.            4    R2 adj   0.011    
 d.f.     45    Pr(> chi2) 0.3068    g        0.033    

 Residuals

      Min       1Q   Median       3Q      Max 
 -0.11768 -0.06635 -0.03056  0.05120  0.25075 


                   Coef    S.E.   t     Pr(>|t|)
 Intercept          0.2077 0.1505  1.38 0.1744  
 abortion_rank12   -0.0016 0.0013 -1.28 0.2063  
 ba_or_more         0.0000 0.0042  0.00 0.9970  
 cig_tax12         -0.0209 0.0211 -0.99 0.3272  
 conserv_advantage -0.0014 0.0024 -0.56 0.5803  


> (ols_hispanic = ols(hispanic ~ abortion_rank12 + ba_or_more + cig_tax12 + conserv_advantage, data = d))
Linear Regression Model

 ols(formula = hispanic ~ abortion_rank12 + ba_or_more + cig_tax12 + 
     conserv_advantage, data = d)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
 Obs      50    LR chi2      3.11    R2       0.060    
 sigma0.1010    d.f.            4    R2 adj  -0.023    
 d.f.     45    Pr(> chi2) 0.5400    g        0.028    

 Residuals

      Min       1Q   Median       3Q      Max 
 -0.13098 -0.05541 -0.02833  0.02041  0.34087 


                   Coef    S.E.   t     Pr(>|t|)
 Intercept          0.1112 0.1586  0.70 0.4869  
 abortion_rank12    0.0011 0.0013  0.87 0.3904  
 ba_or_more         0.0001 0.0044  0.01 0.9899  
 cig_tax12         -0.0069 0.0222 -0.31 0.7558  
 conserv_advantage -0.0014 0.0026 -0.56 0.5764  


> (ols_white = ols(white ~ abortion_rank12 + ba_or_more + cig_tax12 + conserv_advantage, data = d))
Linear Regression Model

 ols(formula = white ~ abortion_rank12 + ba_or_more + cig_tax12 + 
     conserv_advantage, data = d)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
 Obs      50    LR chi2      1.29    R2       0.025    
 sigma0.1344    d.f.            4    R2 adj  -0.061    
 d.f.     45    Pr(> chi2) 0.8635    g        0.024    

 Residuals

      Min       1Q   Median       3Q      Max 
 -0.29197 -0.11366  0.02988  0.08802  0.19163 


                   Coef    S.E.   t     Pr(>|t|)
 Intercept          0.6811 0.2111  3.23 0.0023  
 abortion_rank12    0.0005 0.0018  0.26 0.7938  
 ba_or_more        -0.0001 0.0059 -0.01 0.9903  
 cig_tax12          0.0278 0.0296  0.94 0.3517  
 conserv_advantage  0.0028 0.0034  0.82 0.4167  


> #get predictions
> d$ols_black = predict(ols_black)

> d$ols_hispanic = predict(ols_hispanic)

> d$ols_white = predict(ols_white)

> #inspect predictions
> d %>% select(starts_with("ols")) %>% mutate(sum = rowSums(.))
    ols_black ols_hispanic ols_white sum
1  0.08129250   0.10802760 0.8106799   1
2  0.11823321   0.08031043 0.8014564   1
3  0.14136889   0.07023172 0.7883994   1
4  0.13189772   0.07624972 0.7918526   1
5  0.10280396   0.15374533 0.7434507   1
6  0.12930893   0.11325306 0.7574380   1
7  0.06261715   0.13842026 0.7989626   1
8  0.11806119   0.12687147 0.7550673   1
9  0.11508647   0.10817106 0.7767425   1
10 0.15097533   0.08315063 0.7658740   1
11 0.06825851   0.13257364 0.7991678   1
12 0.09290831   0.11626503 0.7908267   1
13 0.12003177   0.09022122 0.7897470   1
14 0.09482377   0.12513629 0.7800399   1
15 0.14432113   0.07970917 0.7759697   1
16 0.14472842   0.08870156 0.7665700   1
17 0.14109725   0.09872039 0.7601824   1
18 0.15769332   0.06708056 0.7752261   1
19 0.09469234   0.14480432 0.7605033   1
20 0.09043738   0.14094802 0.7686146   1
21 0.10129247   0.11791403 0.7807935   1
22 0.12064080   0.10132318 0.7780360   1
23 0.11602388   0.11910634 0.7648698   1
24 0.16377700   0.09077838 0.7454446   1
25 0.12524587   0.07730750 0.7974466   1
26 0.07060243   0.10921933 0.8201782   1
27 0.12703322   0.11021081 0.7627560   1
28 0.13367696   0.07430030 0.7920227   1
29 0.15029178   0.07798233 0.7717259   1
30 0.10685246   0.12199943 0.7711481   1
31 0.07058313   0.13901865 0.7903982   1
32 0.09142246   0.12212662 0.7864509   1
33 0.10757628   0.12890212 0.7635216   1
34 0.03493501   0.13189804 0.8331670   1
35 0.13176324   0.09598102 0.7722557   1
36 0.14334473   0.06494498 0.7917103   1
37 0.10788090   0.14958898 0.7425301   1
38 0.14930675   0.08299583 0.7676974   1
39 0.09010154   0.12273524 0.7871632   1
40 0.12960586   0.09186749 0.7785267   1
41 0.12586666   0.08545994 0.7886734   1
42 0.12252536   0.09645709 0.7810176   1
43 0.12473591   0.08529285 0.7899712   1
44 0.09099896   0.08246437 0.8265367   1
45 0.16097493   0.09583502 0.7431901   1
46 0.07564659   0.14598479 0.7783686   1
47 0.05854270   0.14244206 0.7990152   1
48 0.09670796   0.09239404 0.8108980   1
49 0.10933095   0.10965351 0.7810155   1
50 0.09307565   0.09622426 0.8107001   1

> #fit accuracies
> d %>% select(starts_with("ols"), black, hispanic, white) %>% cor() %>% .[1:3, 4:6]
                  black    hispanic       white
ols_black     0.3029909 -0.17497676 -0.08992821
ols_hispanic -0.2159748  0.24547482 -0.02824165
ols_white    -0.1708902 -0.04347992  0.15944403

> #multivariate OLS
> (ols_joint = ols(cbind(black, hispanic, white) ~ abortion_rank12 + ba_or_more, data = d))
Linear Regression Model

 ols(formula = cbind(black, hispanic, white) ~ abortion_rank12 + 
     ba_or_more, data = d)

                Model Likelihood     Discrimination    
                   Ratio Test           Indexes        
 Obs     150    LR chi2    342.38    R2       0.898    
 sigma0.1911    d.f.            8    R2 adj   0.892    
 d.f.    141    Pr(> chi2) 0.0000    g        0.314    

 Residuals

           black hispanic    white
 Min    -0.12063 -0.12726 -0.28038
 1Q     -0.06388 -0.05599 -0.10255
 Median -0.03327 -0.02423  0.04949
 3Q      0.05821  0.01528  0.09565
 Max     0.24543  0.34242  0.18706


       Coef    S.E.   t     Pr(>|t|)
  [1,]  0.1552 0.1636  0.95 0.3442  
  [2,] -0.0018 0.0022 -0.82 0.4133  
  [3,]  0.0001 0.0067  0.02 0.9868  
  [4,]  0.0384 0.1636  0.23 0.8146  
  [5,]  0.0013 0.0022  0.62 0.5392  
  [6,]  0.0012 0.0067  0.18 0.8549  
  [7,]  0.8063 0.1636  4.93 <0.0001 
  [8,]  0.0004 0.0022  0.21 0.8378  
  [9,] -0.0013 0.0067 -0.20 0.8419  


> #predictions
> d$joint_ols_black = predict(ols_joint)[, 1]

> d$joint_ols_hispanic = predict(ols_joint)[, 2]

> d$joint_ols_white = predict(ols_joint)[, 3]

> #inspect predictions
> d %>% select(starts_with("joint_ols")) %>% mutate(sum = rowSums(.))
   joint_ols_black joint_ols_hispanic joint_ols_white sum
1       0.09555303         0.11815034       0.7862966   1
2       0.12188779         0.09234947       0.7857627   1
3       0.15017948         0.06705244       0.7827681   1
4       0.14913609         0.07664226       0.7742216   1
5       0.07086325         0.14100841       0.7881283   1
6       0.11448727         0.11617229       0.7693404   1
7       0.07865753         0.14265444       0.7786880   1
8       0.10473606         0.11402244       0.7812415   1
9       0.11151654         0.10446698       0.7840165   1
10      0.14218850         0.08435132       0.7734602   1
11      0.08335854         0.13124113       0.7854003   1
12      0.09180628         0.11898908       0.7892046   1
13      0.11851983         0.09737339       0.7841068   1
14      0.09420884         0.12441664       0.7813745   1
15      0.14521110         0.07551151       0.7792774   1
16      0.13883168         0.08949833       0.7716700   1
17      0.12714583         0.08709083       0.7857633   1
18      0.15582745         0.06610206       0.7780705   1
19      0.08789627         0.13914201       0.7729617   1
20      0.08224830         0.14009239       0.7776593   1
21      0.10274571         0.11314933       0.7841050   1
22      0.12575708         0.09286476       0.7813782   1
23      0.10862763         0.11478390       0.7765885   1
24      0.14372208         0.08017760       0.7761003   1
25      0.13056949         0.08268238       0.7867481   1
26      0.08490326         0.12719051       0.7879062   1
27      0.10986041         0.10728667       0.7828529   1
28      0.13662966         0.08628645       0.7770839   1
29      0.14754681         0.08020051       0.7722527   1
30      0.10152407         0.12076966       0.7777063   1
31      0.07674517         0.14264299       0.7806118   1
32      0.09003875         0.12057784       0.7893834   1
33      0.08785901         0.11761214       0.7945288   1
34      0.07293158         0.14274317       0.7843252   1
35      0.12928101         0.08956415       0.7811548   1
36      0.15418246         0.06904484       0.7767727   1
37      0.07973434         0.13343390       0.7868318   1
38      0.15280485         0.07494186       0.7722533   1
39      0.10672641         0.11489554       0.7783781   1
40      0.12393384         0.09383805       0.7822281   1
41      0.13476186         0.08676737       0.7784708   1
42      0.11662975         0.09760812       0.7857621   1
43      0.13301661         0.08860231       0.7783811   1
44      0.11545267         0.10572082       0.7788265   1
45      0.14112283         0.09369495       0.7651822   1
46      0.07479938         0.14226225       0.7829384   1
47      0.06919598         0.14370501       0.7870990   1
48      0.12051018         0.09824650       0.7812433   1
49      0.09809657         0.10401752       0.7978859   1
50      0.09703090         0.11336115       0.7896079   1

> #accuracies
> d %>% select(starts_with("joint_ols"), black, hispanic, white) %>% cor() %>% .[1:3, 4:6]
                        black    hispanic       white
joint_ols_black     0.2679728 -0.22508994 -0.02574528
joint_ols_hispanic -0.2605606  0.23149306  0.01537504
joint_ols_white    -0.1416356  0.07306988  0.04870974

Thus, in my case, as far as I understand, the OLS does not actually know the values must range from 0-1, and sum to 1, yet this somehow happens with the data when doing regressions one by one, as well as in the joint OLS (i.e. multivariate regression).

My questions are:

  1. How does one generally model outcomes that have a specified range, such as 0-1? This is somewhat different from modeling binary data because while in that case the predictions must be within 0-1 (when transformed from logits), the training data is not binary in this case.
  2. In general, how does one force outcome constraints such as the sum=1 in the above case?
  3. How does OLS know in the above code not to output inappropriate values?

I imagine the above issues can be approached with explicit Bayesian approaches where the priors reflect the outcome ranges, though not sure how one would handle the sum constraint. Is there a frequentist approach as well?

ETA

Stephan Kolassa points to this prior question, which also concerns the issue of proportional data, also called compositional data. However, it differs from the present question in that the question concerns time series data, whereas the present does not, and that I'm explicitly interested in retaining modeling on the proportional data, whereas the approach in the answer question is to transform the data. I am essentially looking for e.g. links to R packages (or maybe Python) that enable modeling on proportional data and predictions on the same scale.

$\endgroup$

marked as duplicate by Stephan Kolassa, Siong Thye Goh, mdewey, Michael Chernick, user158565 Apr 17 at 21:16

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    $\begingroup$ The paper I cite in my answer to the proposed duplicate is concerned with time series forecasting, but the transformation it uses is applicable for more general compositional-data $\endgroup$ – Stephan Kolassa Apr 14 at 5:51
  • $\begingroup$ Thanks, Stephan. That's not entirely what I was looking for, but good to know it's called compositional data. I'll be googling around for that term. $\endgroup$ – Deleet Apr 14 at 6:49
  • 2
    $\begingroup$ I do not think you will be successful in terms of finding a way to model the original proportional data, but having predictions respect the compositional constraints. Could you explain just why you would prefer to model the original proportions? Modeling transformed data is very common, e.g., in logistic regression. Yes, interpretation becomes more difficult, but it's better to have a model that is hard to interpret, rather than no model at all... $\endgroup$ – Stephan Kolassa Apr 14 at 6:58
  • $\begingroup$ Sometimes the models become uninterpretable. For instance, I would want to keep things in original metric so that I can make predictions in that metric. Let's say we have a dataset with missing compositional data, then it would be useful to be able to impute these, but that requires a model that can output appropriate predictions. I did skim the paper you posted, but it doesn't seem to be able to predict values in proper scale (I may be wrong). $\endgroup$ – Deleet Apr 14 at 7:10
  • 2
    $\begingroup$ Well, that paper is explicitly about forecasting, and of course on the original scale. You need to back-transform the predictions from the transformed scale. (You may need some bias-correction in the back-transformation, similar to what you need in Box-Cox transformations, and I don't recall whether that paper goes into that.) $\endgroup$ – Stephan Kolassa Apr 14 at 7:13
1
$\begingroup$

In response to you showing and asking about the results of OLS models, in general, OLS is not appropriate for truncated data, as it will not adjust the estimates of the coefficients to take into account the effect of truncation, resulting in biased coefficients. It is similar to estimating a logistic regression using OLS. But just as with logistic regression, a linear model will sometimes fit just as well, or almost as well, as a non-linear model. (That's why there is something called a linear probability model that models continuous probability values between 0 and 1 using OLS and that is sometimes used by practitioners instead of logistic regression.) The fact that you obtained "appropriate" predicted values is simply due to your data allowing for a good linear approximation of the data generating process.

You'll probably find good examples of modeling compositional data, but for reference, since you asked for general information and are using R, the censReg package allows for estimation of a censored regression model with both lower and upper limits of the dependent variable that can be any numbers, which is a generalization of the the standard Tobit model with a dependent variable left-censored at zero:

https://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf

The model will be estimated using Maximum Likelihood, as is usually done for censored/truncated/limited dependent variable models, not using OLS.

$\endgroup$
  • $\begingroup$ But the values are not approximate in the above case, that's what's so odd. They are exactly summed to 1 (7 digits), but I don't see why that would be the case. In fact, when I wrote the example code, I wanted it to illustrate the problem with using the linear probability model approach, as you also allude to. $\endgroup$ – Deleet Apr 14 at 7:13
  • $\begingroup$ Also, if you decide to use separate models, I recommend you check out Seemingly Unrelated Regression. It's appropriate (and not equivalent to OLS) when your separate regressions have correlated error terms. With this method, the equations are treated as independent from one another, except that the errors are modeled as jointly normally distributed. All equations are estimated simultaneously. $\endgroup$ – AlexK Apr 14 at 7:34
  • $\begingroup$ Note: "SUR is in fact equivalent to OLS ... when each equation contains exactly the same set of regressors on the right-hand-side." $\endgroup$ – nanoman Apr 14 at 10:20
  • $\begingroup$ @nanoman, sure, with one small caveat, as described in the comment to this answer $\endgroup$ – AlexK Apr 14 at 19:00
1
$\begingroup$

First, your OLS models automatically respect the sum-to-1 constraint (up to roundoff error) because both the models and the constraint are linear in the observations. That is, the sum of the predictions for each input is equivalent to the prediction of an OLS model fit to the sum of observations (a constant of 1). What is not guaranteed is that the OLS predictions will be between 0 and 1, although they happen to be for your data.

The multinomial logit model provides one natural way to transform between unconstrained linear-regression predictions and "probabilities" (positive fractions that sum to 1). This is the transformation $y_{it}\mapsto \log(y_{it}/y_{0t})$ in the linked answer (where here you would take $t$ as a sample index, not time).

To make a proper fit systematically, you would need to decide on your error model, which determines the likelihood function. For example, do you imagine that the model can precisely predict probabilities from which the individuals are sampled, and the only error in the population fractions comes from the multinomial distribution based on the population size? If so, you could directly use the likelihood from the multinomial logit model, casting your observations as counts of individuals. But this is implausible for the US state example, because factors not in the model have a much greater effect than the tiny multinomial fluctuations in a state population of millions.

More likely, you just want a reasonable error model with a tractable likelihood, which could be a Gaussian error with unknown variance on the linear predictor functions (log-odds-ratios) of the multinomial logit -- i.e., OLS applied to the transformed observations.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.