Modeling loss aversion with extended-support beta regression

The recently-proposed extended-support beta regression model in R package betareg is illustrated by simultaneously modeling the occurrence and extent of loss aversion in a behavioral economics experiment.

Motivation

To illustrate the benefits of extended-support 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ätzle-Rü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 high-school students in a risky lottery with positive expected payouts. Glätzle-Rü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ätzle-Rützler et al. by employing a similar model for the mean investments but additionally exploring distributional specifications that allow for a probabilistic, rather than mean-only, 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 6-8 or upper grades 10-12?
• 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.

1. 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() or glm() function.

   data("LossAversion", package = "betareg")
   la_ols <- glm(invest ~ grade * (arrangement + age) + male, data = LossAversion)
   summary(la_ols)
   

2. Heteroscedastic censored normal model (CN), also known as heteroscedastic two-limit 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)
   

3. Beta regression (B) after ad-hoc 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)
   

4. Extended-support beta mixture model (XBX) with the same specification as B but adding an extra exceedance parameter to be estimated (instead of the ad-hoc scaling). This can also be fitted with betareg since version 3.2-0 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 log-likelihood, 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 so-called 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 object-oriented implementation of rootograms is available, along with other tools for working with probabilistic models, in the topmodels package on R-Forge (hopefully soon to be submitted to CRAN). You can install it from R-Forge or R-universe and then create the rootograms for the four models:

install.packages("topmodels", repos = "https://zeileis.R-universe.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 10-12, 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 model-based 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 model-based effects, as shown above, the procast() function (for probabilistic forecasts) from the topmodels package can be used. This is again an object-oriented 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 = "10-12")
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 three-part hurdle models for “zero-and-one-inflated” beta regression. However, these models do not work well here: unappealing interpretation, too many parameters, quasi-complete separation of boundary and non-boundary observations. Hence, we do not show the details in this post.