I am working on beta regression model with two grouping variables (farm and years). Climatic variables are my predictors. Response variable is male proportion. I standardized all variables prior to analyses. When I used years as grouping variable, the phi coefficient (intercept) becomes -84. Can this be possible ?

model <- betareg(male_proportion ~ tmin_M +  ppt_M + vpdmin_M + vpdmax_M | years, data = Data, type = "ML")

Standardized weighted residuals 2:
    Min      1Q  Median      3Q     Max 
  -2.945 -0.7879  0.0711  0.5999  2.0158 

Coefficients (mean model with log link):
             Estimate Std. Error  z value Pr(>|z|)    
(Intercept) -0.857335   0.004008 -141.237  < 2e-16 ***
tmin_M      -0.075446   0.006367  -11.193  < 2e-16 ***
ppt_M        0.093383   0.004434   14.343  < 2e-16 ***
vpdmin_M     0.054695   0.006637    7.471  < 2e-16 ***
vpdmax_M    -0.057208   0.006190   -5.566 4.18e-11 ***

Phi coefficients (precision model with log link):
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -84.286687   8.912758  -7.402   <2e-16 ***
      years   0.053939   0.005953   9.871   <2e-16 ***

I have used 4 variables from prism database (minimum temperature, precipitation, vapor pressure deficit minimum and vapor pressure deficit maximum. All these variables are standardized. Male proportion varies within 0-1. Years vary from 1962-2021.

  • 1
    $\begingroup$ Why not? After, -84 is utterly meaningless by itself, just as (say) 17 is utterly meaningless as a temperature (what scale are you using?). Tell us the units of measurement and describe your predictor variables; only then could we supply more specific information. (I see you state you standardized the variables, which is a helpful clue, but it's still not dispositive, because standardized variables can have arbitrarily large ranges, depending on the amount of data, and assuming you mean "standardizing" to zero mean and unit variance.) $\endgroup$
    – whuber
    Apr 15 at 21:20
  • $\begingroup$ @whuber I have added the code that I used. I am confused with -ve value of intercept Phi coefficient. Can that be negative ??? $\endgroup$
    – Rahul
    Apr 16 at 1:36
  • 3
    $\begingroup$ The estimate of phi itself is not negative. The intercept in a model for it is. Take a look at the model for phi in the software you're using $\endgroup$
    – Glen_b
    Apr 16 at 3:50
  • 1
    $\begingroup$ To add, if your first year is 1962 the lowest actually fitted value is about $-84+105$ which is decidedly not negative. Even so you're using a log link, so negative values are certainly possible in the link scale; they are simply less than your log's base in the original scale (and the log link restricts the original scale to be positive). $\endgroup$
    – PBulls
    Apr 16 at 6:43
  • 1
    $\begingroup$ I suspect you are confusing two conventions for formulas. Please read the documentation for betareg, which explains years is not a random effect: it is the variable (along with an implicit intercept) used for the regression of the precision parameter $\phi.$ $\endgroup$
    – whuber
    Apr 17 at 15:03

1 Answer 1


Super-short answer: The phi submodel uses a log-link which assures that all parameter estimates for $\phi$ are positive.

Side remark: There is a serious misunderstanding regarding the | operator in betareg. The terms after | are not random effects but the terms for the linear predictor in the phi submodel. Thus, you fit a log-linear trend for the $\phi$ parameter here. Random effects are not supported by betareg. If you need that, you can consider using glmmTMB or possibly mgcv if using GAM-style random-effects is ok for you.

Details: The model you fitted here assumes $y_i \sim \text{Beta}(\mu_i, \phi_i)$ which you estimated as:

\begin{eqnarray*} \log(\hat \mu_i) & = & x_i^\top \hat \beta \\ & = & -0.86 - 0.075 \cdot \text{tmin}_i + 0.093 \cdot \text{ppt}_i + 0.055 \cdot \text{vpdmin}_i - 0.057 \cdot \text{vpdmax}_i \\ \log(\hat \phi_i) & = & z_i^\top \hat \gamma \\ & = & -84.3 + 0.054 \cdot \text{years}_i \end{eqnarray*}

Thus, for the year 2021 we would get a predicted $\hat \phi = \exp(-84.3 + 0.054 \cdot 2021) \approx 6 \cdot 10^{10}$. This seems to be fishy to me...but maybe years is standardized somehow.

Also, there is a discrepancy between the betareg(...) call and the summary output you posted. The former would imply a logit link for the mu submodel whereas the output claims to use a log link.

Summary: Please review the model specification and its properties carefully, e.g., using vignette("betareg", package = "betareg"). Also re-think the way you want to handle the explanatory variables.

  • $\begingroup$ +1. That's excellent documentation, btw -- many thanks. $\endgroup$
    – whuber
    Apr 17 at 15:04
  • 1
    $\begingroup$ :-) Thank you for your tireless support here on StackExchange!! $\endgroup$ Apr 17 at 15:13
  • $\begingroup$ Thanks so much Achim Zelieis for the nice explanation $\endgroup$
    – Rahul
    Apr 18 at 17:01
  • 1
    $\begingroup$ Please accept the answer by clicking on the check mark on the left. Then it is flagged as resolved here on SO. $\endgroup$ Apr 18 at 18:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.