Beta Regression

1What the method is for

When the response is a proportion, keep the fitted values from leaving the interval by construction.

Let the response \(y_s\) be closed inside \((0,1)\) — say, the proportion of vocabulary lost in settlement \(s\). Fit a linear regression \(y_s = \mathbf{x}_s^\top \boldsymbol{\beta} + \varepsilon_s\) with a normal distribution and the linear predictor ranges over the whole real line, so fitted values fall below 0 or rise above 1. The moment a value appears that cannot be read as a proportion, the model contradicts the domain of its own response.

Two workarounds are commonly reached for, and neither is adequate. One is to dichotomise at a threshold — everything at or above 0.5 counts as 1 — and hand the result to logistic regression, which throws away a quantity that was measured on a continuum. The other is to apply logistic regression directly, but the binomial distribution counts trials that come out 0 or 1; it is not a distribution for a continuous proportion.

Beta regression puts the beta distribution, continuous on \((0,1)\), under the response, and connects its mean to the linear predictor through a link function. The problem of the domain is dissolved by the choice of distribution itself (Ferrari & Cribari-Neto 2004).

2Origin — who made it, and in which field

It came out of a gap: there was no framework for regressing on proportion data.

The formulation is due to Silvia Ferrari and Francisco Cribari-Neto, published in the Journal of Applied Statistics in 2004. The work of two Brazilian statisticians, and the motivation was plain. Rates and proportions turn up in every field, yet no standard regression framework took them as the response. The inadequacy of the two workarounds above is exactly where the paper starts.

It spread quickly, into every area where proportion data are handled. What deserves notice, though, is that the paper invented neither a new distribution nor a new principle of estimation. The beta distribution had long been known, and the generalized linear model framework had been in place since 1972. What Ferrari and Cribari-Neto contributed was assembling existing parts into something usable as a regression model. The lineage in the next section shows this.

3Where it comes from mathematically

Exponential family → generalized linear models → beta regression: a chain of three.

The starting point is the discovery that distributions which look unrelated share a common form. The exponential family, formulated independently by Fisher (1934) and others, is the name for those whose density can be written as follows.

Nelder & Wedderburn (1972) gave a framework in which any distribution from this family can serve as the response under the same fitting procedure. The linear predictor \(\mathbf{x}^\top\boldsymbol{\beta}\) and the mean \(\mu\) of the response distribution are joined by a link function \(g\).

4Definition and assumptions

Rewrite the beta distribution with a mean and a precision, then attach the predictor through a logit link.

The beta distribution is usually written \(\mathrm{Beta}(a,b)\) with two shape parameters, but for regression it helps to handle the mean directly. Setting \(\mu = a/(a+b)\) and \(\phi = a+b\) gives the following form.

That \(\mu\) is attached to the linear predictor through a logit. For the rate of vocabulary loss in settlement \(s\), the main model of the Maibara survey takes the following form.

5Handling the endpoints — when 0 and 1 occur exactly

The beta distribution lives on an open interval. The endpoints carry no density.

The domain is \((0,1)\), not \([0,1]\). If \(y=0\) or \(y=1\) appears in the data the likelihood is undefined and the fit fails. Endpoints are not rare in proportion data, so in practice some treatment is required.

The standard remedy is the linear squeeze given by Smithson & Verkuilen (2006), which pushes values inside the interval relative to the sample size \(n\).

6Maximum likelihood in origin, Bayesian in implementation

The distinction changes how the results are read.

The formulation of Ferrari & Cribari-Neto (2004) is maximum likelihood estimation. The generalized linear model framework itself was assembled as a procedure for obtaining maximum likelihood estimates by iteratively reweighted least squares (Nelder & Wedderburn 1972). In this lineage, a parameter is a fixed unknown value and it is the sample that varies at random.

The settled model of this study, however, is fitted by Bayesian estimation in Stan through brms. Priors are placed on the parameters and the posterior is obtained by MCMC. The intervals reported are credible intervals, not confidence intervals, and \(P(\beta > 0)\) can be read as a probability about the parameter.

7How it was used in this research

Does population growth push up the rate of vocabulary loss, once sampling effort and spatial structure are controlled?

1. Constructing the response　For item \(j\) in settlement \(s\), write \(g_{sj}\) for the proportion of the older generation (aged 60 and over) using it and \(m_{sj}\) for the proportion among middle-school students, and define the loss rate as \(\mathrm{loss}_s = (1/J)\sum_j \max(g_{sj}-m_{sj},\,0)\). The gain rate \(\mathrm{gain}_s\) reverses the direction of the difference.
2. Frequentist baseline　Fit a beta regression with a logit link in betareg, with standardised population size, population growth rate and area as predictors. A first sketch taken on a small sample (n=25).
3. Moving to the hierarchical Bayesian model　Fit a hierarchical beta regression with a random intercept for middle-school district in brms with cmdstanr/Stan. The logs of the numbers surveyed in each generation enter as control terms for sampling effort.
4. Priors　Weakly informative: \(\mathrm{Normal}(0,1)\) on the slopes, \(\mathrm{Normal}(0,2)\) on the intercept, \(\mathrm{Gamma}(0.1,0.1)\) on the precision \(\phi\).
5. Sampling　Four chains, 4,000 iterations each (2,000 as warm-up), adapt_delta = 0.97, max_treedepth = 12. The adopted model has \(\hat{R}=1.00\), bulk-ESS > 3,000, and at most one divergent transition.
6. Checking the spatial term　Build a binary adjacency matrix treating settlements within 3 km as neighbours, and compare models with and without an ICAR spatial random effect by approximate leave-one-out cross-validation.

Hierarchical Bayesian beta regression (n=25, random intercept for middle-school district). Coefficients are on the logit scale and predictors are standardised. Figures are quoted from loss_growth_coefs.csv in the repository and from run results recorded in the working log.

8Limits — what the method cannot say

What the choice of distribution resolves is the problem of the domain, and nothing else.

9Sources and data

A secondary analysis of a three-generation dialect vocabulary survey covering five middle-school districts from Maibara in Shiga Prefecture to Sekigahara in Gifu Prefecture. The questionnaire elicited 106 dialect vocabulary variants across 40 question concepts, classified in advance into 45 Kansai/Ōmi forms, 31 Tokyo/standard forms and 26 Mino/Owari forms, with four items of youth speech and new dialect. The original survey was conducted in 2015 by Habiro-Net, and this research analyses the anonymised data. Every stage of the analysis is scripted.