¿Cómo puedo calcular la respuesta de regresión beta utilizando los coeficientes?

Question

Para este ejemplo supongamos que tenemos un dato Reading skills de un betareg paquete.

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

Y tengo un modelo de regresión beta del ejemplo:

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

Así que tengo unos coeficientes :

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 

Si calculo la respuesta utilizando el predict() función que me da:

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

Sé que por defecto se utiliza la función de enlace logit. Así que estoy tratando de calcular esto manualmente para la primera observación con los coeficientes anteriores y me da:

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

¿Dónde me equivoco?

Answer 1

La razón es que dyslexia utiliza contrastes de suma en lugar de los contrastes de tratamiento por defecto, véase ?ReadingSkills . Esto significa que se codifica -1/1 en lugar de 0/1. La razón de esto es para que coincida con los resultados del original de Smithson & Verkuilen (2006, Métodos psicológicos ) papel.

Así, la primera fila de la matriz del modelo es:

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

En consecuencia, la predicción es (utilizando un poco de redondeo para simplificar):

$$ \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*} $$

El resultado exacto se puede obtener mediante:

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