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.)