Why are my beta regression results biased?

Question

I would like to use the betareg package, and started with some simulations to ensure I understand how it works.

I seem to be getting biased coefficient estimates in my simulation. I have made my simulation so that it roughly matches the type of data I plan to use (i.e., I have a continuous covariate, but data for evenly-spaced values of it).

# Covariate
theta <- rep(c(0, 1, 2, 3, 4), 200)

# Set a slope of 1, with no intercept 
beta <- 1 

# Define the linear predictor
eta <- theta*beta

# Mean mu equals logistic function of the linear predictor
mu <- exp(eta)/(1 + exp(eta))

# If we set the first shape parameter to 1, then the following 
# is the relationship between mu and the second shape parameter 
b <- 1/mu - 1

# Draw responses from beta 
resp <- rbeta(length(mu), 1, b) 

# There will be some 1s, so squeeze the data using the strategy
# from the betareg documentation 
resp.squeeze <- (resp*(length(resp) - 1) + 0.5)/length(resp) 

# Run beta regression
br <- betareg(resp.squeeze ~ theta | theta)
summary(br)

My summary call outputs these coefficients for the conditional $\mu$ model:

Coefficients (mean model with logit link):
             Estimate  Std. Error z value Pr(>|z|)    
(Intercept)  0.19401    0.06850   2.832  0.00462 ** 
theta        0.77837    0.03605  21.593  < 2e-16 ***

So, it is estimating a nonzero intercept when it was designed to have a zero intercept, and a slope that is also incorrect.

I can recover the correct estimates if I change the example so that the linear predictor does not go too high, so that there are fewer responses epsilon away from one. Is the problem purely a function of needing to squeeze the data? I am just worried about using it with my actual data, given that the squeeze in the simulation is quite small, the model is correctly specified -- but the coefficients are still off. I have a lot of data in my data set (so the squeeze will also be small), but obviously I do not know that the data was actually generated through a process that follows a conditional beta.

Any advice would be much appreciated. I would love to avoid the 0-1 inflated models, if I can. (The betareg package has a ton of great functionality I would like to use.)

Thanks!

Answer 1

The parameterization you use is not the same as in betareg hence you cannot recover your original parameters. The rbeta() function has parameters shape1 and shape2 while betareg() employs mu and phi. The mapping is mu = shape1/(shape1 + shape2) and phi = shape1 + shape2 and conversely shape1 = mu * phi and shape2 = (1 - mu) * phi.

The problem is that in your example you correctly set up mu so that it is linked to a linear predictor in theta. However, you do not specify phi but instead shape1 directly. The consequence is that the implicit log(phi) is not a linear function of theta as fitted in betareg().

To avoid this problem, set up a phi parameter using a log-link and say an intercept of 1 but a slope of zero. Additionally, I shift the regressor theta somewhat to avoid the mu values extremely close to 1.

theta <- rep(-2:2, 200)
mu <- exp(0 + 1 * theta) / (1 + exp(0 + 1 * theta))
phi <- exp(1 + 0 * theta)

And then you can simulate the response

set.seed(1)
y <- rbeta(length(mu), shape1 = mu * phi, shape2 = (1 - mu) * phi)

And then betareg() essentially recovers your parameters:

betareg(y ~ theta | theta)
## Call:
## betareg(formula = y ~ theta | theta)
## 
## Coefficients (mean model with logit link):
## (Intercept)        theta  
##    -0.01135      1.01680  
## 
## Phi coefficients (precision model with log link):
## (Intercept)        theta  
##     1.01608     -0.01465  

Answer 2

I believe the problem is either not convex or not locally quadratic, so a more robust optimizer than BFGS, like simulated annealing ("SANN") is needed.

The likelihood function I wrote below (zoib.ll) prints out the intercept and slope at each iteration, so you can see it jumping out of local minima. Using "Nelder-Mead" as the optimizer is much faster but only works if the problem is not too hard (e.g.,, rnorm(1000, 3, 1)).

library(extraDistr) # provides a parametrization of beta in terms of the mean
# and the precision
theta <- rnorm(10000, 5, 1)
beta <- 1
rho <- 1
eta.mu <- theta*beta
eta.phi <- theta*rho
mu <- exp(eta.mu)/(1 + exp(eta.mu))
phi <- exp(eta.phi)
util <- rprop(length(mu), size = phi, mean = mu)
data <- data.frame(util, theta)

zoib.ll <- function(par, data){
  x <- data$theta
  y <- data$util
  i.zero <- ifelse(y == 0, 1, 0)
  i.one <- ifelse(y == 1, 1, 0)
  i.gr <- ifelse(y > 0, 1, 0)
  i.btwn <- ifelse((y > 0) & (y < 1), 1, 0) 

  eta.mu <- par[1] + par[2]*x
  eta.zero <- par[3] + par[4]*x
  eta.one <- par[5] + par[6]*x
  eta.v <- par[7] + par[8]*x

  print(par[1])
  print(par[2])

  p.zero <- exp(eta.zero)/(1 + exp(eta.zero))
  p.one <- exp(eta.one)/(1 + exp(eta.one))  
  mu <- exp(eta.mu)/(1 + exp(eta.mu))
  v <- exp(eta.v)

  lik.zo <- i.zero*log(p.zero) + i.gr*log(1-p.zero) + i.one*log(p.one) + i.btwn*log(1-p.one)
  lik.b <- i.btwn*(lgamma(v) - lgamma(v*mu) - lgamma(v*(1-mu)) + (v*mu - 1)*log(y) + (v*(1-mu) - 1)*log(1-y))
  lik.b[is.nan(lik.b)] <- 0
  ll <- sum(lik.zo + lik.b)
  return(-ll)
}

inits <- c(0, 0, 0, 0, 0, 0, 1, 1)
opt1 <- optim(par = inits, zoib.ll, data = data, method = "SANN")
opt1$par