Modeling loss aversion with extendedsupport beta regression
Motivation
To illustrate the benefits of extendedsupport beta regression models, suggested in a recent arXiv paper with Ioannis Kosmidis, we revisit the analysis of a behavioral economics experiment conducted and published by GlätzleRützler et al. (2015, Journal of Economic Behavior & Organization, doi:10.1016/j.jebo.2014.12.021). The outcome variable is the proportion of tokens invested by highschool students in a risky lottery with positive expected payouts. GlätzleRützler et al. focused on the effects of several experimental factors on the mean investments, which reflect the players’ willingness to take risks. In their study they employed linear regression models, estimated by ordinary least squares (OLS) with standard errors adjusted for potential clustering and heteroscedasticity.
Here, we extend the analysis from GlätzleRützler et al. by employing a similar model for the mean investments but additionally exploring distributional specifications that allow for a probabilistic, rather than meanonly, interpretation of the effects. From an economic perspective this is of interest because it allows to interpret both the mean willingness to take risks in this experiment, and the probability to behave like a rational Homo oeconomicus, who would invest (almost) all tokens in this lottery because it has positive expected payouts.
The full replication code for the analyses from the arXiv paper is available in lossaversion.R with some auxiliary functions in beta01.R. Below we only provide the most important R snippets to provide a feeling for the workflow in R. The rest of the discussion here highlights the main insights from the analysis.
An aggregated version of the data from all nine rounds of the experiment is available as LossAversion
in the betareg
package. Interest is in linking the variable invest
with the proportion of the total tokens invested in all nine rounds to explanatory information:
grade
: Is the player from lower grades 68 or upper grades 1012?arrangement
: Is the player an individual or a team of two?male
: Is (at least one of) the player(s) male?age
: (Average) age of the player(s) in years.
Models
We compare four different models for invest
which all employ the same equation for the mean submodel.
And all except the OLS reference model employ the main effects of the three experimental factors for the dispersion submodel.

Normal linear model (N) with constant variance, corresponding to the OLS approach from the original study. In R, this can be fitted with the base
lm()
orglm()
function.data("LossAversion", package = "betareg") la_ols < glm(invest ~ grade * (arrangement + age) + male, data = LossAversion) summary(la_ols)

Heteroscedastic censored normal model (CN), also known as heteroscedastic twolimit tobit model in econometrics. This can be fitted with the crch package (for censored regression with conditional heteroscedasticity).
library("crch") la_htobit < crch(invest ~ grade * (arrangement + age) + male  arrangement + male + grade, data = LossAversion, left = 0, right = 1) summary(la_htobit)

Beta regression (B) after adhoc scaling of the investments to the open unit interval (to avoid the boundary observations). This can be fitted with the betareg package.
library("betareg") LossAversion$invests < (LossAversion$invest * (nrow(LossAversion)  1) + 0.5)/ nrow(LossAversion) la_beta < betareg(invests ~ grade * (arrangement + age) + male  arrangement + male + grade, data = LossAversion) summary(la_beta)

Extendedsupport beta mixture model (XBX) with the same specification as B but adding an extra exceedance parameter to be estimated (instead of the adhoc scaling). This can also be fitted with
betareg
since version 3.20 with XBX regression being automatically selected in case of boundary observations in the response.la_xbx < betareg(invest ~ grade * (arrangement + age) + male  arrangement + male + grade, data = LossAversion) summary(la_xbx)
If you run the code and compare the model summaries, note that the coefficients from N and CN use an identity link for the mean parameter whereas B and XBX use a logit link. In addition, the loglikelihood, and, hence, AIC and BIC, are comparable only between CN and XBX because those two models have the same support for the response variable, that is the unit interval with point masses at 0 and 1. See also the accompanying arXiv paper for the full summary tables. By and large, the mean parameters for the N and CN models are rather similar and those for B and XBX are rather similar, only with some differences that occur due to using models with point masses on the boundaries (CN and XBX) or not (N and B).
However, instead of studying the individual estimated coefficients in more detail we rather assess the models graphically by visualizing their goodness of fit and different types of fitted effects.
Goodness of fit
To illustrate how different the fitted probability distributions of the four models are, we employ socalled hanging rootograms. These compare the empirical marginal distributions of the response variable (proportion of tokens invested) to the aggregated fitted distributions from the models. The quality of the fit can be judged by the deviation of the hanging bars from the zero reference line.
An objectoriented implementation of rootograms is available, along with other tools for working with probabilistic models, in the topmodels package on RForge (hopefully soon to be submitted to CRAN). You can install it from RForge or Runiverse and then create the rootograms for the four models:
install.packages("topmodels", repos = "https://zeileis.Runiverse.dev")
library("topmodels")
rootogram(la_ols, breaks = 6:16 / 10, main = "N")
rootogram(la_htobit, main = "CN")
rootogram(la_beta, main = "B")
rootogram(la_xbx, main = "XBX")
A more refined version of the plots is shown below. See the full replication script linked above for the code details.
The square root of the expected frequencies are shown as red dots and the square root of the observed frequencies are hanging from the points as gray bars. The dashed lines are the Tukey warning limits at +/ 1. These plots show that models N and B fit poorly in the tails. In contrast, models CN and XBX fit very well with almost all bars hanging close to the zero reference line. The fit for XBX appears to be slightly better than for CN.
Effects
While models XBX and CN provide a much better probabilistic fit than their uncensored counterparts B and N, it turns out that the predicted mean investments from all four models are still very similar. But XBX and CN allow for interpretations beyond the mean, including economically relevant interpretations of probability effects.
For illustration, we focus on the team arrangement effect for a subsample with a large share of very rational subjects: male players or teams with at least one male, in grades 1012, and between 15 and 17 years of age. The figure below shows the estimated arrangement effect for the mean E(Y), i.e., the expected proportion of tokens invested, and for the probability to behave very rationally and invest almost everything, i.e., P(Y > 0.95). The empirical quantities are shown in black for the subsample between 15 and 17 years of age while the modelbased effects are shown at an age of 16.
The graphic shows that all models do a reasonable job in estimating E(Y) but the censored models XBX and CN are much better at estimating the probability to behave very rationally. Similarly, it could be shown that the fit for the probability P(Y < 0.05) is also much better for XBX and CN than for B and N but this is not done here, because that probability does not have such an appealing economic interpretation like P(Y > 0.95).
For obtaining the modelbased effects, as shown above, the procast()
function (for probabilistic forecasts) from the topmodels
package can be used. This is again an objectoriented implementation that facilitates obtaining not only moments (such as means and variances) but also entire probability distributions (as S3 objects) and corresponding probabilities, densities, and quantiles.
Here, we only briefly show the code for the fitted XBX model but the same function calls can be applied to the other fitted model objects. First, we set up the new data that only varies arrangement
from single to team but keeps all other variables fixed. Then, both kinds of effects are computed with procast()
.
la_nd < data.frame(arrangement = c("single", "team"), male = "yes", age = 16, grade = "1012")
procast(la_xbx, newdata = la_nd, type = "mean")
## mean
## 1 0.4713
## 2 0.6861
procast(la_xbx, newdata = la_nd, type = "cdf", at = 0.95, lower.tail = FALSE)
## probability
## 1 0.07161
## 2 0.18501
Thus, the mean invested proportion goes up from 47.1% to 68.6% for teams vs. single players in this setting, while the probability to behave almost fully rationally increases from 7.2% to 18.5%.
Again, the full code for creating the figure and underlying table is provided in the replication script linked above.
The script also inlucdes some further illustrations, e.g., the comparison with threepart hurdle models for “zeroandoneinflated” beta regression. However, these models do not work well here: unappealing interpretation, too many parameters, quasicomplete separation of boundary and nonboundary observations. Hence, we do not show the details in this post.