2
$\begingroup$

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

data("ReadingSkills", package = "betareg") 
head(ReadingSkills)
##   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?

$\endgroup$

1 Answer 1

1
$\begingroup$

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
$\endgroup$
1
  • 1
    $\begingroup$ 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/… $\endgroup$
    – huberttt
    Commented Jun 17, 2020 at 20:56

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.