# Simulate a distribution from a fitted beta-regression model for a density plot in R [duplicate]

I have produced the following figure by simulating some values from a fitted gamma regression with a low AIC value that provides the closest approximation of my raw data out of all of my models, and the code used to create it.

> gam_age_tiss <- glm(value ~ age*tissue, data = SEAB_pivot, family = Gamma(link = "log"))
> gam_age_tiss

Call:  glm(formula = value ~ age * tissue, family = Gamma(link = "log"),
data = SEAB_pivot)

Coefficients:
(Intercept)                 ageHY                 ageJU
0.24801              -1.28749              -0.34595
tissuefeathers  ageHY:tissuefeathers  ageJU:tissuefeathers
0.36543               1.05740              -0.07576

Degrees of Freedom: 71999 Total (i.e. Null);  71994 Residual
Null Deviance:      22020
Residual Deviance: 7987     AIC: 60260

> summary(gam_age_tiss)

Call:
glm(formula = value ~ age * tissue, family = Gamma(link = "log"),
data = SEAB_pivot)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)           0.248006   0.002545   97.43   <2e-16 ***
ageHY                -1.287488   0.004409 -292.03   <2e-16 ***
ageJU                -0.345954   0.004409  -78.47   <2e-16 ***
tissuefeathers        0.365432   0.003600  101.52   <2e-16 ***
ageHY:tissuefeathers  1.057399   0.006235  169.59   <2e-16 ***
ageJU:tissuefeathers -0.075761   0.006235  -12.15   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.1166247)

Null deviance: 22016.1  on 71999  degrees of freedom
Residual deviance:  7987.2  on 71994  degrees of freedom
AIC: 60259

Number of Fisher Scoring iterations: 4

> sims1 <- simulate(gam_age_tiss, nsim = 50)
> plot(density(SEAB_pivot$value),main = "gam_age_tiss") > for (i in 1:50) lines(density(sims1[[i]]), col = "red")  I would like to create a similar figure for a fitted beta regression, but the simulate function from the {stats} r package that I used to sample from the gamma distribution apparently only works with glm objects. I have been trying to manually extract the the two shape parameters (alpha and beta) from my fitted beta regression, so that I can just use a function like sample or rbeta, without any luck, and am not understanding moment matching for mu specifically from the betareg output. Can someone help me either: a) simulate samples from the distribution of the fitted betareg as done with the gamma regression, b) extract the alpha and beta parameters from the betareg output, or c) help me use the coefficients in the betareg to manually calculate alpha and beta? Here is the density plot of raw data and the code I used for the betareg. > plot(density(beta_merge$eccentricity),main = "Beta Regression on Year (lowest AIC)")

> beta_year <- betareg(eccentricity ~ year, data = maxmin)
> beta_year

Call:
betareg(formula = eccentricity ~ year, data = maxmin)

Coefficients (mean model with logit link):
(Intercept)     year2021
0.9243      -0.6400

Phi coefficients (precision model with identity link):
(phi)
4.361

> byear <- summary(beta_year)
> byear

Call:
betareg(formula = eccentricity ~ year, data = maxmin)

Standardized weighted residuals 2:
Min      1Q  Median      3Q     Max
-2.8116 -0.6589  0.0613  0.7658  3.1518

Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.92432    0.04596   20.11   <2e-16 ***
year2021    -0.64002    0.06238  -10.26   <2e-16 ***

Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi)   4.3613     0.2011   21.68   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type of estimator: ML (maximum likelihood)
Log-likelihood: 210.4 on 3 Df
Pseudo R-squared: 0.126
Number of iterations: 13 (BFGS) + 1 (Fisher scoring)


See Cribari-Neto, F., & Zeileis, A. (2010). Beta regression in R. Journal of Statistical Software, 34(2), Article 1. https://doi.org/10.18637/jss.v034.i02 for formula and definition in betareg. It specifies that $$\alpha = p = \mu \phi$$, $$\beta = q = (1 - \mu) \phi$$, $$\mu = \alpha / (\alpha + \beta)$$, $$\phi = \alpha + \beta$$, where $$\alpha, \beta, p, q$$ are shape parameters, $$\mu$$ is the predicted mean, and $$\phi$$ is the predicted precision parameter. According to the function arguments of

predict(object, newdata = NULL,
type = c("response", "link", "precision", "variance", "quantile"),
na.action = na.pass, at = 0.5, ...)


We can extract $$\alpha$$ and $$\beta$$ by

predict(beta_year, type = "response") * predict(beta_year, type = "precision")
(1 - predict(beta_year, type = "response")) * predict(beta_year, type = "precision")


You estimated a fixed dispersion beta regressions. betareg reports $$\phi = \alpha + \beta$$ as the precision parameter, which is specified as a constant in your model. $$1/\phi = 1/(\alpha + \beta)$$ is known as the dispersion parameter. You can also have a variable dispersion beta regression that associates the precision parameter with certain predictors. The default link function of the precision parameter is "log" to ensure a positive precision parameter, unless formula is of type y ~ x where the default is "identity".

betareg(eccentricity ~ year | year, link = "logit", link.phi = "log", data = maxmin)

• The development version of betareg on R-Forge also provides a simulate() method. You can install it via install.packages("betareg", repos = "R-Forge.R-project.org"). Commented Feb 18 at 1:15
• @AchimZeileis I removed my old betareg package and installed the development version, but when I try to run it I get this error: > betareg::simulate(beta_year, nsim = 50) Error: 'simulate' is not an exported object from 'namespace:betareg' Commented Feb 19 at 17:04
• @DrTerryTao I used the predict function to generate vectors of alpha and beta, which I plugged into rbeta and was able to plot. Thank you for your help - I think I was just expecting a single value for each of the parameters which was what threw me off. Commented Feb 19 at 17:26
• Maybe this depends on the R version you are using. Adding type = "source" might help: install.packages("betareg", repos = "https://R-Forge.R-project.org", type = "source"). Afterwards running library("betareg") and then example(betareg) and simulate(gy, nsim = 3) should work. Commented Feb 20 at 22:35