# How do I compute beta regression response using coefficients?

For this example assume we have a data Reading skills from a betareg package.

data("ReadingSkills", package = "betareg")
##   accuracy dyslexia     iq
## 1  0.88386       no  0.827
## 2  0.76524       no  0.590
## 3  0.91508       no  0.471
## 4  0.98376       no  1.144
## 5  0.88386       no -0.676
## 6  0.70905       no -0.795


And I have a beta regression model from the example:

rs_beta <- betareg(accuracy ~ dyslexia * iq | dyslexia + iq,
data = ReadingSkills, hessian = TRUE)



So I have a coefficients :

coef(rs_beta)
##       (Intercept)          dyslexia                iq       dyslexia:iq
##         1.1232251        -0.7416450         0.4863696        -0.5812569
## (phi)_(Intercept)    (phi)_dyslexia          (phi)_iq
##         3.3044312         1.7465642         1.2290731


If I compute the response using the predict() function it gives me:

predict(rs_beta, newdata = ReadingSkills[1,], type = "response")
## 1
## 0.9397888


I know that by default it is using logit link function. So I am trying to compute this manually for the first observation with the coefficients above and it gives me:

$$\begin{eqnarray*} \mathrm{logit}(r) & = & 1.1232251 + 0.4863696 \cdot 0.827 = 1.5245 \\ r & = & 0.82134 \end{eqnarray*}$$

Where do I make a mistake?

The reason is that dyslexia uses sum contrasts rather than the default treatment contrasts, see ?ReadingSkills. This means that it is coded -1/1 rather than 0/1. The reason for this is to match the results from the original Smithson & Verkuilen (2006, Psychological Methods) paper.

Thus, the first row of the model matrix is:

model.matrix(rs_beta)[1,]
## (Intercept)    dyslexia          iq dyslexia:iq
##       1.000      -1.000       0.827      -0.827


Consequently, the prediction is (using a little bit of rounding for simplicity):

$$\begin{eqnarray*} \mathrm{logit}(\mu) & = & 1.1232 - 0.7417 \cdot (-1) + 0.4864 \cdot 0.827 - 0.5813 \cdot (-0.827) = 2.7479 \\ \mu & = & 0.9398 \end{eqnarray*}$$

The exact result can be obtained by:

plogis(model.matrix(rs_beta)[1,] %*% coef(rs_beta, model = "mean"))
##           [,1]
## [1,] 0.9397888

• Indeed thanks a lot, just a moment ago I have noticed it when I read a first example from the pdf below (I attach , maybe it will help someone) where instead of classical dummy coding the -1/1 coding was implemented for both variables. And thank you for your clarification! researchgate.net/publication/… Commented Jun 17, 2020 at 20:56