Open Access
Open Peer Review

This article has Open Peer Review reports available.

How does Open Peer Review work?

# A Monte Carlo simulation study comparing linear regression, beta regression, variable-dispersion beta regression and fractional logit regression at recovering average difference measures in a two sample design

BMC Medical Research Methodology201414:14

DOI: 10.1186/1471-2288-14-14

Received: 30 August 2013

Accepted: 21 January 2014

Published: 24 January 2014

The Erratum to this article has been published in BMC Medical Research Methodology 2016 16:152

## Abstract

### Background

In biomedical research, response variables are often encountered which have bounded support on the open unit interval - (0,1). Traditionally, researchers have attempted to estimate covariate effects on these types of response data using linear regression. Alternative modelling strategies may include: beta regression, variable-dispersion beta regression, and fractional logit regression models. This study employs a Monte Carlo simulation design to compare the statistical properties of the linear regression model to that of the more novel beta regression, variable-dispersion beta regression, and fractional logit regression models.

### Methods

In the Monte Carlo experiment we assume a simple two sample design. We assume observations are realizations of independent draws from their respective probability models. The randomly simulated draws from the various probability models are chosen to emulate average proportion/percentage/rate differences of pre-specified magnitudes. Following simulation of the experimental data we estimate average proportion/percentage/rate differences. We compare the estimators in terms of bias, variance, type-1 error and power. Estimates of Monte Carlo error associated with these quantities are provided.

### Results

If response data are beta distributed with constant dispersion parameters across the two samples, then all models are unbiased and have reasonable type-1 error rates and power profiles. If the response data in the two samples have different dispersion parameters, then the simple beta regression model is biased. When the sample size is small (N0 = N1 = 25) linear regression has superior type-1 error rates compared to the other models. Small sample type-1 error rates can be improved in beta regression models using bias correction/reduction methods. In the power experiments, variable-dispersion beta regression and fractional logit regression models have slightly elevated power compared to linear regression models. Similar results were observed if the response data are generated from a discrete multinomial distribution with support on (0,1).

### Conclusions

The linear regression model, the variable-dispersion beta regression model and the fractional logit regression model all perform well across the simulation experiments under consideration. When employing beta regression to estimate covariate effects on (0,1) response data, researchers should ensure their dispersion sub-model is properly specified, else inferential errors could arise.

### Keywords

Regression modelling Linear regression Beta regression Variable-dispersion beta regression Fractional Logit regression Beta distribution Multinomial distribution Monte Carlo simulation

## Background

In biomedical research it is common to encounter response variables which have support on the interval (0,1). These types of response variables may arise in the form of proportions/percentages, or certain types of fractions and rates. The traditional approach to analyzing these types of response data – across virtually all scientific disciplines - is via linear regression. If desired, the response variable can be transformed prior to estimation of the linear regression parameters. This transformed linear model may improve diagnostic performance; however, this may render interpretation of estimated regression parameters challenging. Alternatively, the beta distribution allows specification of a probability model for continuous random variables with support over the interval (0,1). For many years statisticians have exploited the flexibility of the beta distribution in theoretical modelling exercises; however, its use in applied research settings has not garnered equal attention. Johnson et al. [1] cite numerous instances where the beta distribution has been used in theory/practice and champion increased use of the beta distribution in applied research settings. Gupta et al. [2] also cite numerous applications where the beta distribution provides a useful probability generating model for continuous data with support on the interval (0,1). However, neither of these extensive resources on the beta distribution cites a regression modelling framework for estimating covariate effects on beta distributed response variables. Recent developments by Paulino [3], Ferrari and Cribrari-Neto [4], Smithson and Verkuilen [5] and others have resulted in a more general purpose beta regression machinery. The variable-dispersion beta regression model [5] will be used extensively in our simulation experiments, as it is particularly useful for modelling covariate effects on response variables which are assumed to follow a beta distribution. The beta regression model extends on ideas of generalized linear models [6] both in terms of their specification and estimation. Use of the beta regression model has been increasing in recent years. In slides from an unpublished presentation given by Ferrari [7], the author suggests over 100 instances where beta regression has been used in theoretical and applied research settings. Some application areas include: medicine, veterinary science, pharmacology, odontology, hydrobiology, nutritional science, forest science, waste management, education, political science, economics and finance. Clearly, embedding the beta distribution within a more general regression modelling framework has enhanced its uptake in applied research settings. A final model which we consider for estimating the average proportion/percentage/rate difference in our two-sample model is the fractional logit regression model. The fractional logit model is a popular model for fractional response variables in econometrics and was proposed (independently) by Papke and Wooldridge [8] and by Cox [9]. The fractional logit model is similar to generalized linear regression models [6]; however, it does not make any fully parametric assumptions regarding the distributional form for the response variable. Rather the fractional logit model only specifies a parametric form for the conditional mean and conditional variance of the response. The form of the conditional mean and variance functions are chosen to ensure admissible predictions/fitted-values from such models. In this case, the model specification is chosen to ensure predictions/fitted values from the fractional logit model fall in the interval (0,1). The estimator proposed by Papke and Wooldridge [8] is the one we pursue in this manuscript as they specify forms for robust variance estimators which have more desirable coverage/power properties than the more traditional quasi-likelihood models proposed by Cox [9].

Given the recent popularity of the beta regression model, especially in biomedical research, we thought it prudent to compare linear regression, beta regression, variable-dispersion beta regression and fractional logit regression models for estimating covariate effects on a response variable which lives on the interval (0,1). To accomplish this goal we conducted a Monte Carlo simulation experiment where we generated response variables following different (parametric) probability generating models. First, we considered simulating response data from the continuous beta distribution with support on (0,1). This experiment allows us to compare models when we know the beta regression model is properly specified given the response data. Specifically, we can investigate efficiency gains which may be observed from specifying an appropriate statistical model to observed response data. Additionally, we simulate response data from the discrete multinomial distribution with probability mass observed only on a finite number of points in (0,1). This experiment allows us to investigate model performance when response data is non-continuous. In this case, all models are incorrectly specified given the response data. This experiment allows us to investigate whether estimated regression models are robust to non-continuous response data. In all scenarios, we fit linear regression, beta regression, variable-dispersion beta regression, and fractional logit regression models to these randomly generated response data and compared the finite sample statistical properties of the respective estimators. We are particularly interested in the ability of each estimator to recover the average proportion/percentage/rate differences from a simple two sample design. In terms of statistical properties we will compare the respective estimators in terms of: bias, variance, type-1 error and power. Understanding the performance of these models on simulated datasets (where population parameters are known) is important for applied researchers who must discern whether to estimate covariate effects on (0,1) response data using the traditional linear regression model or more novel regression models, such as beta regression, variable-dispersion beta regression and fractional logit regression models.

## Methods

### Statistical methods

#### The linear regression model

The linear regression model is a workhorse of applied statisticians. It is used to model the effect of continuous/categorical covariates on a scalar response (assumed to be generated according to a Gaussian probability model). Thorough introductions to the linear regression model are given in Weisberg [10], McCullagh and Nelder [6], and White [11].

In this study we consider a simple two sample problem, re-cast under a regression framework, such that our response variable is modelled as a function of a single intercept parameter and a single slope parameter. The linear model and its conditional mean function look as follows:
${Y}_{i}={\mathit{\beta }}_{0}+{\mathit{\beta }}_{1}{X}_{i1}+{\mathit{ϵ}}_{i}$
The notation above suggests that we observe a vector of response variables, Y1…Yn. Further, we have information on a single binary covariate, Xi {0,1}, again for i = 1…n. The regression coefficients β0 and β1 are estimated from the data. Estimation and inferential procedures are justified given that the following assumptions are satisfied [11]:
1. 1.

The model is properly specified

2. 2.

X is a non-stochastic and finite dimensional (n by p) matrix with n ≥ p

3. 3.

(XTX) is non-singular and hence invertible

4. 4.

E(ϵi) = 0 (i = 1…n)

5. 5.

ϵi ~ Normal(0, σ2) (i = 1…n)

In our experiment, we are interested in the ability of the linear regression estimator to recover the average proportion/percentage/rate difference given our simple two sample design. Taking linear combinations of the estimated model parameters we arrive at the following estimator:
$\mathit{\Delta }=E\left({Y}_{i}\phantom{\rule{0.25em}{0ex}}|\phantom{\rule{0.25em}{0ex}}{X}_{i1}=1\right)-E\left({Y}_{i}\phantom{\rule{0.25em}{0ex}}|\phantom{\rule{0.25em}{0ex}}{X}_{i1}=0\right)={\mathit{\beta }}_{1}$

Therefore a test of ∆ = 0 is equivalent to a test of β1 = 0. In this simulation we carry out such a test using a Wald statistic, W, which follows an asymptotic standard normal distribution. We reject the null hypothesis in instances where |W| > 1.96 (corresponding to an α = 5% significance threshold).

#### The beta regression model (and some extensions)

The beta regression model was proposed by Paulino [3], Ferrari and Cribrari-Neto [4], and Smithson and Verkuilen [5] for modelling covariate effects on a continuous response variable which assumes support on the interval (0,1).

The beta distribution is thoroughly described in Johnson et al. [1] and Gupta [2]. The beta density is a very flexible density, assuming support on the interval (0,1). The most common parameterization of the beta density is in terms of its two shape parameters {p,q}:
$f\left(y;p,q\right)=\frac{\mathit{\Gamma }\left(p+q\right)}{\mathit{\Gamma }\left(p\right)\mathit{\Gamma }\left(q\right)}{y}^{p-1}{\left(1-y\right)}^{q-1}$
In the parameterization given above we assume p > 0, q > 0, y (0,1) and use (·) to denote the gamma function (a generalization of the factorial function to non-integer arguments). The density assumes probability mass on the interval (0,1) and is zero elsewhere. Further, under this parameterization we define the mean and variance of the random variable, Y, as follows:
$E\left(Y\right)=\frac{p}{p+q}$
$\mathit{\text{VAR}}\left(Y\right)=\frac{\mathit{pq}}{{\left(p+q\right)}^{2}\left(p+q+1\right)}$
Above E(·) and VAR(·) denote the expectation and variance operators, with respect to the given beta distribution. In regression modelling, it is more common to parameterize the density in terms of a mean (μ) and dispersion parameter (ϕ) instead of two shape parameters, {p,q}. In this parameterization we have the following relationships:
$\mathit{\mu }=\frac{p}{p+q}$
$\mathit{\varphi }=p+q$

This implies: p = μϕ and q = (1 – μ)ϕ.

Given the above relationships we can derive the mean and the variance of the beta density in terms of a mean and dispersion parameter as follows:
$E\left(Y\right)=\mathit{\mu }$
$\mathit{\text{VAR}}\left(Y\right)=\frac{V\left(\mathit{\mu }\right)}{1+\mathit{\varphi }}=\frac{\mathit{\mu }\left(1-\mathit{\mu }\right)}{1+\mathit{\varphi }}$
Given a fixed value for the mean, the larger the value of ϕ the smaller the variance of the response variable, Y (and vice-versa). Under this new parameterization, in terms of a mean and dispersion parameter, the density of Y looks as follows:
$f\left(y;\mathit{\mu },\mathit{\varphi }\right)=\frac{\mathit{\Gamma }\left(\mathit{\varphi }\right)}{\mathit{\Gamma }\left(\mathit{\mu },\mathit{\varphi }\right)\mathit{\Gamma }\left(\left(1-\mathit{\mu }\right)\mathit{\varphi }\right)}{y}^{\mathit{\mu }\mathit{\varphi }-1}{\left(1-y\right)}^{\left(1-\mathit{\mu }\right)\mathit{\varphi }-1}$
In Figure 1, we graphically represent some of the forms the beta density can take on for different values of {p,q}, or alternatively, {μ, ϕ}.
The beta regression model, and the variable-dispersion extensions which we will discuss in this study are being increasingly utilized to model covariate effects on response variables observed on the interval (0,1). The beta regression model is an obvious choice for modelling response data which follow a beta distribution. Consider the scenario where we observe response data Y1…Yn on the interval (0,1). The beta regression model assumes that the mean of these random variables, can be represented in the following form:
$g\left({\mathit{\mu }}_{\mathrm{i}}\right)={\mathit{\eta }}_{\mathrm{i}}={\mathit{\beta }}_{0}+{\mathit{\beta }}_{1}{X}_{i1}$

Our link function g(·) can be any function which is strictly monotone, twice differentiable, and maps the response variable observed on the interval (0,1) to the real line. The most commonly used link function in beta regression is the logit link. Alternative link functions include: the probit, the complementary log-log, the log-log and the Cauchy link. In general, any inverse cumulative distribution function will be an appropriate link function in a beta regression framework as they act to map the interval (0,1) to the real line.

The components of the basic beta regression model can be summarized as:
1. 1.

A response variable from a beta distribution

2. 2.

A linear predictor, ηi

3. 3.

A suitable link function, such that: E(Y i |X i ) = g(μ i) = η i

Given the above components, the log-likelihood of the beta regression model can be written as follows:
$\begin{array}{ll}\phantom{\rule{1em}{0ex}}\mathit{LL}\left({\mathit{\mu }}_{i},\mathit{\varphi }\right)=& \mathit{\text{log}}\mathit{\Gamma }\left(\mathit{\varphi }\right)-\mathit{\text{log}}\mathit{\Gamma }\left({\mathit{\mu }}_{i}\mathit{\varphi }\right)-\text{log}\left(\left(1-{\mathit{\mu }}_{i}\right)\mathit{\varphi }\right)\\ +\left({\mathit{\mu }}_{i}\mathit{\varphi }-1\right)\text{log}\left({y}_{i}\right)\\ +\left\{\left(1-{\mathit{\mu }}_{i}\right)\mathit{\varphi }-1\right\}\text{log}\left(1-{y}_{i}\right)\end{array}$

The log-likelihood function can be maximized numerically as described in Ferrari and Cribrari-Neto [4]. The mean and dispersion parameter estimates are known to be biased, especially in small samples. Kosmidis and Firth [12] discuss the issue of finite sample bias in beta regression. The authors propose a general purpose algorithm for producing bias-reduced and bias-corrected parameter estimates via adjustments to the score function. In our simulation experiment we estimate parameters from the beta regression model via standard maximum likelihood (ML) methods, as well as the bias-reduced (BR) and bias-corrected (BC) methods. In our simulation experiments we employ the simple ML estimators; however, we note that BC/BR methods may improve type-1 error rates in small sample situations.

The beta regression model proposed above assumes that the dispersion parameter is constant for all individuals under consideration. In many biomedical applications this may be an unrealistic assumption (especially if one expects a non-zero mean difference across categorical groups). As its name implies, the variable-dispersion beta regression model [5] allows the value of the dispersion parameter to vary across individuals. Further, the value of the dispersion parameter can actually be modelled as a function of covariates. The variable-dispersion beta regression model is a type of double-index regression model [13], as it contains two regression equations, one modelling the mean as a function of covariates and the other modelling the dispersion as a function of covariates.

Again, we consider the scenario where we observe response data Y1…Yn on the interval (0,1). The variable-dispersion beta regression model assumes that the mean and dispersion of these random variables can be represented in the following form:
$g\left({\mathit{\mu }}_{\mathrm{i}}\right)={\mathit{\eta }}_{\mathrm{i}}={\mathit{\beta }}_{0}+{\mathit{\beta }}_{1}{X}_{i1}$
$h\left({\mathit{\varphi }}_{\mathrm{i}}\right)={\zeta }_{\mathrm{i}}={\gamma }_{0}+{\gamma }_{1}{X}_{i1}$

Once again, we assume that both g(·) and h(·) are strictly monotonic, twice differentiable functions which act to map the mean, μi, and the dispersion, ϕ i , to the real line. Once again, suitable choices of g(·) include the following link functions: logit, probit, complementary log-log, log-log, Cauchy or any other inverse cumulative distribution function. The link function for h(·) is typically chosen to be the log link. The identity link can also be used; however, it has the undesirable property of possibly suggesting non-positive values of ϕ i .

The log-likelihood function for the variable-dispersion beta regression model can be numerically maximized and is subject to similar finite sample biases as the basic beta regression model. Below, we illustrate the log-likelihood function for this model:
$\begin{array}{ll}\phantom{\rule{1em}{0ex}}\mathit{LL}\left({\mathit{\mu }}_{i},\mathit{\varphi }\right)=& \mathit{\text{log}}\mathit{\Gamma }\left({\mathit{\varphi }}_{i}\right)-\mathit{\text{log}}\mathit{\Gamma }\left({\mathit{\mu }}_{i}{\mathit{\varphi }}_{i}\right)-\text{log}\left(\left(1-{\mathit{\mu }}_{i}\right){\mathit{\varphi }}_{i}\right)\\ +\left({\mathit{\mu }}_{i}{\mathit{\varphi }}_{i}-1\right)\text{log}\left({y}_{i}\right)\\ +\left\{\left(1-{\mathit{\mu }}_{i}\right){\mathit{\varphi }}_{i}-1\right\}\text{log}\left(1-{y}_{i}\right)\end{array}$
In the case of both the beta regression model and the variable-dispersion (double-index) beta regression models, estimates of mean and dispersion parameters {β,γ} are achieved by numerically solving the likelihood equations given above. The resulting parameter estimates are asymptotically normally distributed and take the following form:
$\left(\begin{array}{c}\stackrel{̂}{\mathit{\beta }}\\ \stackrel{̂}{\gamma }\end{array}\right)~\mathit{\text{MVN}}\left(\left(\begin{array}{c}\mathit{\beta }\\ \gamma \end{array}\right),{C}^{-1}\right)$
Further,
${C}^{-1}={C}^{-1}\left(\mathit{\beta },\gamma \right)=\left(\begin{array}{cc}{C}^{\mathit{\beta }\mathit{\beta }}& {C}^{\mathit{\beta }\mathit{\gamma }}\\ {C}^{\mathit{\gamma }\mathit{\beta }}& {C}^{\mathit{\gamma }\mathit{\gamma }}\end{array}\right)$

For our purposes it suffices to realize that the estimators of the mean and dispersion parameters are consistent estimators of their target parameters and are distributed according to a multivariate normal distribution, with variance-covariance matrix C-1. Detailed derivations of these formulas (particularly pertaining to the forms of the C-1 matrix) are given in Ferrari and Cribrari-Neto [4].

Again, we are interested in the ability of the (variable-dispersion) beta regression estimator to recover the average proportion/percentage/rate difference given our two sample design. In all of our simulation experiments we assume a logit link for the mean function. The (default) identity link is used in the beta regression modelling context and the log link is used in the variable-dispersion beta regression context. In all scenarios, our target of inference is the average proportion/percentage/rate difference and we view the terms in the dispersion sub-model as a nuisance. A point estimator of the proportion/percentage/rate difference from the beta regression model is:
$\mathit{\Delta }=\mathit{\text{log}}\left(\frac{1}{1+\text{exp}\left(-{\mathit{\beta }}_{0}-{\mathit{\beta }}_{1}{X}_{i1}\right)}\right)-\mathit{\text{log}}\left(\frac{1}{1+\text{exp}\left(-{\mathit{\beta }}_{0}\right)}\right)$

We use the delta method to estimate the variance and standard error of this estimator, respectively. We construct a Wald style test of the null hypothesis that ∆ = 0. The Wald statistic, W, is computed as the ratio of the difference in proportion/percentage/rates over the estimated standard error. The test statistic is presumed to follow an asymptotic standard normal distribution. Again, we use a 5% critical threshold for rejecting the null hypothesis (this corresponds to rejection of H0 if |W| > 1.96)

#### The fractional logit regression model

The final methodology we consider for estimating average proportion/percentage/rate differences in our two-sample design is the fractional logit regression model [8, 9]. The fractional logit regression model is most commonly encountered in the econometrics literature and has been demonstrated as being an effective means for estimating covariate effects on a response variable which lives on (0,1). Hence we consider it in this manuscript – as a result, introducing health services researchers to yet another plausible strategy for modelling proportions/percentages/fractions/rates.

The fractional logit regression model is considered a quasi-parametric regression model. In other words, the fractional logit regression model does not make any parametric assumption regarding the distribution of the response variable being modelled; rather, it makes assumptions regarding only the first two conditional moments of the response variable – the conditional mean and the conditional variance. The choice of the conditional mean and conditional variance function are typically made to ensure that predictions/fitted-values from the specified model are admissible. In our case, this implies that the predictions/fitted-values fall in the interval (0,1).

As mentioned above, quasi-likelihood models typically only make assumptions regarding the first two conditional moments of the response variable [6, 9]. The conditional variance is assumed to be a known function of the mean (up to a scale parameter) and the conditional mean function therein is assumed to be a function of unknown model parameters:

In the first Equation V(·) denotes the variance operator, σ2 is a scale parameter which is estimated from observed data. ν(·) is a known variance function, and μi is the mean function. In the second Equation E(·) denotes the variance operator, g(·) is a known link function and βj represent the unknown mean function parameters which must be estimated from the data. Yi and Xi represent the response variable and observed covariates, respectively.

In describing the fractional logit model we adopt the terminology of Papke and Wooldridge [8]. Our chief assumption relates to the specification of the conditional mean function, namely:
$E\left({Y}_{i}\right)=h\left({\mathit{\mu }}_{i}\right)$
Generally, h(·) is a known function which maps our real valued linear predictor into the interval (0,1). Again, their exist many plausible function which could accomplish this goal, in this manuscript we choose h(·) to be the logistic function and arrive at the fractional logit model. That is:
$h\left({\mathit{\mu }}_{i}\right)=\frac{\mathit{\text{exp}}\left({\mathit{\mu }}_{i}\right)}{1+\mathit{\text{exp}}\left({\mathit{\mu }}_{i}\right)}=\frac{1}{1+\mathit{\text{exp}}\left(-{\mathit{\mu }}_{i}\right)}$
Further, the conditional variance of the response variable is assumed to be:

Papke and Wooldridge [8] argue that this conditional variance assumption is too restrictive for modelling response data with support over (0,1). Therefore, in their manuscript they offer two alternative strategies: first, using robust/sandwich estimators of the variance-covariance matrix and second, adjusting the estimated variance-covariance matrix by the Pearson scale adjustment factor. We considered both approaches; however, noted little difference in performance between the two estimators of the variance-covariance matrix. Hence, we report on only the fractional logit model with sandwich/robust variance-covariance matrix.

Parameter estimation under the fractional logit model proceeds by maximizing the following Bernoulli quasi-likelihood function:
$\mathit{LL}={y}_{i}*\mathit{\text{log}}\left(h\left({\mathit{\mu }}_{i}\right)\right)+\left(1-{y}_{i}\right)*\mathit{\text{log}}\left(h\left({\mathit{\mu }}_{i}\right)\right)$

#### Monte Carlo simulation design

The goal of this simulation experiment is to compare the properties of the linear regression model, the beta regression model, the variable-dispersion beta regression model and the fractional logit regression model at recovering estimates of average proportion/percentage/rate differences from a simple two sample design. In all experiments we simulate data from parametric probability generating models such that the observed response data is on the interval (0,1). Subsequently, we estimate covariate effects on the response variable using one of four regression models: the linear regression model, the beta regression model, the (double-index) variable-dispersion beta regression model and the fractional logit regression model. Given estimates of average proportion/percentage/rate differences from the respective models, we compare statistical properties of the respective estimators, such as: bias, variance, type-1 error and power [14, 15]. We investigate finite sample performance of each of the estimators by varying the sample size within each unique simulation experiment. In all scenarios, the sample size in group 1 is set equal to the sample size in group 2. Group specific sample sizes under consideration in this simulation are: 25, 100, 250, and 750. The total sample size for a given simulation experiment is double the group-specific sample size (as this experiment assumes a 2-sample design). In each instance we consider 20,000 replications of each experiment. We choose 20,000 replicate simulations such that coverage in the type-1 error experiments is based off of approximately 1000 rejections of a true null hypothesis. We present mean estimates of bias, variance, type-1 error and power averaged across the 20,000 replicate simulations. Further we present Monte Carlo error estimates of bias, variance and power. Detailed derivations of Monte Carlo error are described in White [16]. The “seeds” which govern the pseudo-randomness of the various Monte Carlo experiments are given in the attached R/SAS codes.

The first parametric probability model which we consider for generating response data on the interval (0,1) is the beta distribution. Table 1 describes the parameter values used to generate randomly simulated beta response variables. The response variables are generated such that certain mean and dispersion properties are achieved. For example, mean differences of zero are used to assess the type-1 error rates of respective estimators (for both fixed and varying dispersion). Further, non-zero mean differences are used to assess power (again for both fixed and varying dispersion). In this experiment response data are generated as independent draws from the respective beta distributions. That is, observations within and between the two samples are independently distributed. Within the type-1 error and power experiment frameworks, respectively, we have 3 sub-experiments: the first set of experiments consider the scenario where the central tendency of the simulated response distribution is near the center of the support (0.5); the second set of experiments considers the effect of shifting the central tendency to the right such that it is centered near 0.25; and finally, the last experiment considers the effect of shifting the central tendency to the boundary of the support, near 0.05. As sub-scenarios we vary the shape of the beta distribution when data are simulated from the center, right-center and far-right of the support, considering scenarios where the simulated data are symmetric and other scenarios where the simulated data is highly skewed. As the data are beta distributed we expect the beta regression models to perform well in all scenarios; however, we anticipate that the linear model will perform well when data are symmetric and unimodal. That is, we expect the linear model to perform well as the shape/rate parameters both become large and as the ratio of the shape/rate parameters approach 1 (resulting in a symmetric and unimodal beta distribution – which converges to that of a normal distribution).
The next parametric probability model under consideration is the discrete multinomial model which takes probability mass only on a finite number of points on the interval (0,1). More specifically, we assume our response variable Yi can take on the following values:
${Y}_{i}\in \left\{0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95\right\}$
That said, we do not assume the probability of assuming these values is necessarily uniform. Rather, we assign a vector of probabilities to these points, corresponding to the relative likelihood that the response variable assumes that particular value. Table 2 describes the particular probability vectors used to generate response variables for each group in our two sample design. Once again, we vary the expected value of the response to assess differences in type-1 error rates and power across our linear regression, beta regression, (double-index) variable-dispersion beta regression and fractional logit regression models. Again, in this experiment response data are generated as independent draws from the respective multinomial distributions. That is, observations within and between the two samples are independently distributed.

#### Statistical software

This simulation experiment was conducted using R version 3.02 [17] and results were also verified using SAS 9.3 [18].

Simulation of the beta and multinomial response variables were carried out using the rbeta() and rmultinom() functions, respectively. Linear regression modelling was performed using the lm() function. Beta regression was performed using the betareg() function in the betareg library [13]. Fractional logit regression models were estimated using the glm() function and the sandwich() function [19]. Standard errors for the proportion/percentage/rate differences from beta regression and fractional logit regression models were calculated using the deltamethod() function in the msm library [20].

SAS PROC NLMIXED was used to specify the linear regression model, beta regression model, variable-dispersion beta-regression model and fraction logit regression model likelihood equations, respectively, and model parameters were estimated via likelihood methods.

All R and SAS code used to conduct this simulation can be obtained by contacting the corresponding author.

## Results

Detailed results of the Monte Carlo simulation study are given in Tables 3, 4, 5 and 6. Tables 3 and 4 describe the type-1 error and power experiments, respectively, given response data simulated according to independent draws from various parameterizations of the beta distribution. Tables 5 and 6 describe the type-1 error and power experiments, respectively, given response data simulated according to independent draws from various parameterizations of the multinomial distribution.

Table 3 describes the results of the type-1 error experiment (∆ = 0) given response data distributed according to independent draws from a beta distribution. The top half of Table 3 illustrates results when the dispersion parameter is equal across groups; whereas, the bottom half of Table 3 illustrates results when the dispersion parameter varies as a function of group membership. As probability mass moves away from the center of the support (i.e. 0.5) and towards the boundary of the support (0 or 1) we observe that the beta regression model provides biased estimates of the average proportion/percentage/rate difference between the two samples when the dispersion parameters vary as a function of group membership. For example, when μ0 = μ1 = 0.25 and ϕ 0 = 5 and ϕ 1 = 10 we observe biased estimates of effect from the beta regression model (biases range from 2.27E-02 through 2.41E-02). As the dispersion parameters increase (i.e. ϕ 0 = 100 and ϕ 1 = 200) the observed bias in the beta regression model is slightly attenuated (biases range from 1.22E-03 through 1.26E-03). Similar findings are observed when the mean parameters are adjusted, such that μ0 = μ1 = 0.05. Table 4 describes the results of the power experiments (∆ = 0.025) given response data distributed according to independent draws from a beta distribution. Near identical results are observed as were discussed for the type-1 error experiments in Table 3. That is, when the respective means are near the boundary of the support, and the dispersion parameters vary as a function of group membership the beta regression model can yield biased estimates of effect. When the dispersion parameters are small (ϕ 0 = 5 and ϕ 1 = 10) the bias in effect estimates is appreciable (biases range from 2.22E-02 through 2.31E-02). On an absolute scale these biases are meaningful; however, when expressed on a relative scale these biases are even more pronounced. As the dispersion parameters increase in magnitude the magnitude of the bias in the beta regression models is attenuated (biases range from 9.73E-04 through 1.20E-03). Given that the simple beta regression model is biased in certain scenarios we eliminate it from consideration in the results/discussion sections which follow.

In small sample scenarios, when N0 = N1 = 25, the linear regression model had a mean type-1 error of approximately 0.050; whereas, the variable-dispersion beta regression model had mean type-1 error rate of 0.058 and the fractional logit regression model had a mean type-1 error rate of 0.060. As the sample size is increased to 100 per group, 250 per group and 750 per group, respectively, the type-1 error rates of the linear regression model, the variable dispersion beta regression model and the fractional logit regression model became more similar. Further, improvements in the type-1 error rate of the variable-dispersion beta regression model for small samples (N0 = N1 = 25) were observed when we used bias corrected/reduced estimation methods instead of the more traditional ML estimation methods (results not shown; however, can be verified by modifying simulation codes in R).

When considering the power experiments estimates of average bias across the 20,000 replicate experiments were small for the linear regression model, the variable-dispersion beta regression model and the fractional logit regression model (of magnitude 1E-04 through 1E-06 respectively). Further, estimates of average variance across the 20,000 replicate simulations were similar across the linear regression model, the variable-dispersion beta regression model and the fractional logit regression model. These findings imply the estimators have similar average mean squared error. That said, the power for estimated variable-dispersion beta regression models and the fractional logit regression models, respectively, marginally exceeded that of the linear regression model across all simulation experiments considered.

Table 5 describes the results of the type-1 error experiment (∆ = 0) given response data distributed according to independent draws from a multinomial distribution. For the type-1 error experiments all estimators are relatively free of bias. The magnitudes of estimated biases are similar for the linear regression model, variable-dispersion beta regression model and the fractional logit model. Again, average variance across the 20,000 replicate simulations were similar for all models. Type-1 error rates are closest to the desired 5% level for the linear regression model. Again the variable-dispersion beta regression model and the fractional logit regression model have elevated type-1 error rates when sample sizes are small (N0 = N1 = 25). Table 6 describes the results of the power experiment (∆ = 0.10) given response data distributed according to independent draws from a multinomial distribution. When data are simulated according to either a symmetric or asymmetric discrete multinomial distribution we observe that the beta regression model is biased. In the symmetric case biases are attenuated (biases range from 2.32E-03 through 2.55E-03) compared to the asymmetric case (biases range from 1.35E-02 through 1.38E-02). The magnitude of the bias in the linear regression estimator and the fractional logit regression estimator are similar. However, in the case of discrete data we notice that the variable-dispersion beta regression model has slightly elevated mean bias levels. That said, the variable-dispersion beta regression model is slightly more powerful than the linear regression model and the fractional logit regression model (however, this is likely an artifact of the difference in magnitudes of bias in these models). Among models with comparable biases, the fractional logit model is more powerful than the linear regression model when data are generated from a discrete multinomial distribution on (0,1).

## Discussion

The main findings of this Monte Carlo simulation study are summarized in Tables 3, 4, 5 and 6 in the results section. In general, properties of the respective estimators are similar regardless of whether the underlying data generating mechanism is beta distributed (Tables 3 and 4) or multinomial distributed (Tables 5 and 6). Hence we will discuss findings from the type-1 error experiments and the power experiments in general, as results seem to hold irrespective of the probability generating models. We note interesting exceptions where warranted.

Considering the type-1 error experiments (Table 3 and Table 5) we observe that the linear regression model, the variable-dispersion beta regression model and the fractional logit regression model provide unbiased estimates of our population proportion/percentage/rate difference (∆ = 0) under all simulated scenarios. The magnitudes of bias tend to be similar across estimators, ranging from 1E-04 through 1E-06. In many circumstances the simple beta regression model also provide unbiased estimates of our null (∆ = 0) effect. However, in circumstances where the dispersion parameter varied between groups, the simple beta regression model demonstrated fairly substantial bias in its attempt to recover the average population proportion/percentage/rate difference. The impact of non-constant dispersion amongst individuals in this simulation experiment were more pronounced when the dispersion parameters were small (e.g. ϕ 0 = 5 and ϕ 1 = 10) compared to when the dispersion parameters were large (e.g. ϕ 0 = 100 and ϕ 1 = 200). Further, the effects of non-constant dispersion between groups appear more pronounced when the group means are near the boundary of the distributions support (0 or 1) compared to when they are near the center of the support (½). This is demonstrated by observed biases in the beta regression model of about 0.02 units in certain circumstances (Table 3). It is interesting to note that in terms of type-1 error rates the linear regression model performed well regardless of sample size; whereas, the variable-dispersion beta regression model and fractional logit regression model experienced slightly elevated type-1 error rates when the group specific sample sizes were small (N0 = N1 = 25). Another important point is that improvements in the small sample type-1 error rates of the beta regression estimators could be achieved by using the bias corrected/reduced estimation methods in place of the more traditional ML estimators. These BC/BR estimators are easily implemented in the R betareg() procedure [12, 13].

Considering the power experiments (Table 4 and Table 6) we again observe that the linear regression model, the variable-dispersion beta regression model and the fractional logit regression model provide (relatively) unbiased estimates of our proportion/percentage/rate difference (∆ = 0.025 in the beta distributed simulations and ∆ = 0.10 in the multinomial distributed simulations). The magnitude of the average biases across the 20,000 replicate experiments is similar across these three models when the data are beta distributed (Table 4); however, when the data arise from a multinomial distribution the variable-dispersion beta regression model has slightly elevated bias levels compared to the linear regression model and fractional logit regression model. That said, on a relative (or absolute) scale, the observed biases in the variable-dispersion beta regression model are not overly large. Again, in cases where the dispersion parameter varies across groups we observe that the simple beta regression model has trouble recovering the desired epidemiological effect measure. Again, this problem is more pronounced when the dispersion parameters are small and the group means are situated near the boundary of the support. The beta regression model also struggles at recovering the desired difference measure in the multinomial experiment where the response variable is skewed (Table 6); however, demonstrates more comparable performance to the linear regression model, the variable-dispersion beta regression model and the fractional logit model when the response variable is simulated from a symmetric multinomial distribution. In general, the linear regression model, the variable-dispersion beta regression model and the fractional logit regression model perform well in terms of recovering unbiased estimates of the non-zero effect measure. The models have similar power profiles across the continuous beta distributed simulation experiments – with minor power advantages appearing in the variable-dispersion beta regression models and the fractional logit regression models. Further when response data are distributed according to a discrete multinomial distribution minor advantages in power appear for the variable-dispersion beta regression model (at the cost of small magnitude increases in bias) and the fractional logit model compared to the linear regression model (Table 6).

The results of this Monte Carlo simulation study indicate that the linear regression model, the variable-dispersion beta regression model and the fractional logit regression model are capable of producing unbiased estimates of average proportion/percentage/rate differences given response data observed on the interval (0,1) from a two sample design. The simple beta regression model struggles if the dispersion sub-model is incorrectly specified. When sample sizes are small, type-1 error rates appear closer to the nominal 5% level in the linear regression model. The variable-dispersion beta regression model and fractional logit model appear slightly more powerful than the linear regression model when a non-zero difference between groups is present.

A similar study was conducted by Kieschnick and McCullough [21]. In their article they made similar conclusions favouring the (variable-dispersion) beta regression model and the fractional logit regression model for estimating covariate effects on response data observed on (0,1). In their article they dismissed the linear regression model, because in the more complex regression scenarios they were considering it could lead to inadmissible predictions (e.g. predicted values outside of (0,1)). In our simulation experiment we are not necessarily interested in the predictions or fitted values, rather we are interested in the ability of our model to recover the average difference in proportions/percentages/rates across a two-sample design. This is typically viewed as an absolute measure of “effect” in epidemiological research. If one is interested in this measure of effect, rather than in model based predictions then it appears from this simulation that the linear regression model performs similarly well as the more novel variable-dispersion beta regression model and the fractional logit model. That said, if ones interest lies in predicted/fitted-values then it likely behoves the researcher to choose a model which will result in admissible predictions.

Introduction of the beta distribution into a general regression framework has resulted in enhanced attention/use of the beta distribution by both theoretical and applied researchers. Ferrari [7] provides many examples in which the beta regression model has been applied in theoretical/applied modelling exercises. Two biomedical applications where the beta regression framework has been implemented include modelling scales scores, such as SF-6D response data [22] and modelling stroke lesion volumes [23]. Many other biomedical applications of the beta regression model exist as suggested by Ferrari [7].

The purpose of this Monte Carlo simulation experiment was to investigate the properties of the linear regression model, the beta regression model, the variable-dispersion beta regression model and the fractional logit regression model at recovering average proportion/percentage/rate differences from a two sample design. The simplicity of the design aids in interpreting properties of the respective models. In addition, the two sample design is one of the most commonly encountered study designs by epidemiologists and biostatisticians. Hence the simulation study answers a very important question for applied epidemiologists/biostatisticians, namely: given a more traditional linear regression framework for modelling covariate effects on response data observed on the interval (0,1) are there any benefits in fitting a more novel beta regression model, variable-dispersion beta regression model or fractional logit regression model to these same data and using it for inference? Results of this simulation study suggest that the more novel variable-dispersion beta regression model and fractional logit regression model have comparable properties to the traditional linear regression model. While the linear regression model may perform better in terms of type-1 error rates in small samples, the variable-dispersion beta regression model and fractional logit regression model seem slightly more powerful at detecting a true non-zero difference between groups in a two-sample design. Conversely, the simple beta regression model appears to struggle if the dispersion sub-model is incorrectly specified. Given this finding applied researchers should be cautious in fitting off the shelf beta regression models to their (0,1) response data. If a choice is made to fit a beta regression model to observed data practitioners must strive to ensure correct specification of both the mean and dispersion sub-models in order to generate proper inferences.

## Conclusion

The purpose of this Monte Carlo simulation study was to compare the properties of the linear regression model to the more novel beta regression, variable-dispersion beta regression and fractional logit regression models at recovering estimates of average proportion/percentage/rate differences in a two-sample design. We observe that the simple beta regression model is biased if the dispersion sub-model is incorrectly specified. The variable-dispersion beta regression model is unbiased (given proper specification of the dispersion sub-model). The fractional logit regression model is also an unbiased estimator of effect. Moreover, the power and type-1 error profiles are very similar for the linear model, variable-dispersion beta regression model and the fractional logit regression model. These results seem to suggest promise for the beta regression model going forward; however, for the time being applied researchers should be cautious in applying off the shelf beta regression algorithms to their response data observed on the interval (0,1) and should strive to ensure correct specification of both the mean and dispersion sub-models such that proper inferences are generated from the observed data.

## Declarations

### Acknowledgements

The authors would like to acknowledge the contributions of the assigned BMC Medical Research Methodology editor – Monica Taljaard – and two reviewers Jay Verkuilen and Maarten Buis. The reviewers’ comments helped to improve the manuscript greatly.

## Authors’ Affiliations

(1)
Department of Family and Community Medicine, University of Toronto

## References

1. Johnson N, Kotz S, Balakrishnan N: Continuous Univariate Distributions. 1995, Hoboken, New Jersey: Wiley, 2Google Scholar
2. Gupta A, Nadarajah S: Handbook of Beta Distribution and its Applications. 2004, Boca Raton, Florida: CRC Press, 1Google Scholar
3. Paolino P: Maximum likelihood estimation of models with beta distributed dependent variables. Polit Anal. 2001, 9 (4): 325-346. 10.1093/oxfordjournals.pan.a004873.
4. Ferrari S, Cribrari-Neto F: Beta regression for modelling rates and proportions. J Appl Stat. 2004, 10: 1-18.Google Scholar
5. Smithson M, Verkuilen J: A better lemon squeezer? Maximum-likelihood regression with beta distributed dependent variables. Psychol Methods. 2006, 11 (1): 54-71.
6. McCullagh P, Nelder J: Generalized linear models. 1989, Boca Raton: CRC Press, 2
7. Ferrari S: Beta Regression Modelling: Recent Advances and in Theory and Applications. 2013, Unpublished presentation: http://www.ime.usp.br/~sferrari/13EMRslidesSilvia.pdf Google Scholar
8. Papke L, Wooldridge J: Econometric methods for fractional response variables with an application to 401(K) plan participation rates. J Appl Econ. 1996, 11: 619-632. 10.1002/(SICI)1099-1255(199611)11:6<619::AID-JAE418>3.0.CO;2-1.
9. Cox C: Non-linear quasi-likelihood models: applications to continuous proportions. Comput Stat Data Anal. 1996, 21 (4): 449-461. 10.1016/0167-9473(95)00024-0.
10. Weisberg S: Applied Linear Regression. 2005, Hoboken, New Jersey: Wiley, 3
11. White H: Asymptotic Theory for Econometricians. 2000, San Diego, California: Academic PressGoogle Scholar
12. Kosmidis I, Firth D: A generic algorithm for reducing bias in parametric estimation. Electron J Stat. 2010, 4: 1097-1112. 10.1214/10-EJS579.
13. Grun B, Kosmidis I, Zeileis A: Extended beta regression in R: shaken, stirred, mixed and partitioned. J Stat Softw. 2012, 48 (11): 1-25.
14. Knight K: Mathematical Statistics. 2000, Boca Raton, Florida: CRC PressGoogle Scholar
15. Wasserman L: All of Statistics: A Concise Course in Statistical Inference. 2004, New York, New York: Springer
16. White I: SIMSUM: analyses of simulation studies including Monte Carlo Error. Stata J. 10 (3): 369-385.Google Scholar
17. R Development Core Team: R: A language and environment for statistical computing. 2013, Vienna, Austria: R Foundation for Statistical Computing, http://www.R-project.org/, 3-900051-07-0,Google Scholar
18. SAS Institute: 2013, North Carolina, USA, http://www.sas.com/en_us/legal/editorial-guidelines.html,
19. Zeileis A: Econometric computing with HC and HAC covariance matrix estimators. J Stat Softw. 2004, 11 (10): 1-17.
20. Jackson C: Multi-state models for panel data: the msm package for R. J Stat Softw. 2011, 38 (8): 1-29.
21. Kieschnick R, McCullough B: Regression analysis of variates observed on (0,1): percentages, proportions and fractions. Stat Model. 2003, 3 (3): 193-213. 10.1191/1471082X03st053oa.
22. Hunger M, Beaumert J, Holle R: Analysis of SF-6D index data: is beta regression appropriate?. Value Health. 2011, 14: 759-767. 10.1016/j.jval.2010.12.009.
23. Swearingen C, Tilley B, Adams R, Rumboldt Z, Nicholas J, Bandyopadhyay D, Woolson R: Application of beta regression to analyze ischemic stroke volume in NINDS rt-PA clinical trials. Methods in Neuroepidemiology. 2011, 37: 73-82. 10.1159/000330375.
24. ### Pre-publication history

1. The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1471-2288/14/14/prepub