DFFIT para la regresión beta

Aug 26 2020

Estoy tratando de calcular DFFITS para GLM, donde las respuestas siguen una distribución Beta. Usando betaregel paquete R. Pero creo que este paquete no es compatible influence.measures()porque al usar dffits() Code

require(betareg)
df<-data("ReadingSkills")
y<-ReadingSkills$accuracy
n<-length(y)

bfit<-betareg(accuracy ~ dyslexia + iq, data = ReadingSkills)
DFFITS<-dffits(bfit, infl=influence(bfit, do.coef = FALSE))
DFFITS

produce

Error en if (model$rank == 0) { : el argumento es de longitud cero

Soy un novato en R. No sé cómo resolver este problema. Por favor, ayúdeme a resolver esto o deme algunos consejos a través del código R sobre cómo calcular los DFFIT manualmente. Saludos

Respuestas

2 jay.sf Aug 26 2020 at 15:12

dffitsno están implementados para "betareg"objetos, pero puede intentar calcularlos manualmente.

De acuerdo con este Stack Overflow Q/A , podríamos escribir esta función:

dffits1 <- function(x1, bres.type="response") {
  stopifnot(class(x1) %in% c("lm", "betareg"))
  sapply(1:length(x1$fitted.values), function(i) {
    x2 <- update(x1, data=x1$model[-i, ]) # leave one out
    h <- hatvalues(x1)
    nm <- rownames(x1$model[i, ])
    num_dffits <- suppressWarnings(predict(x1, x1$model[i, ]) - 
                                     predict(x2, x1$model[i, ]))
    residx <- if (class(x1) == "betareg") {
      betareg:::residuals.betareg(x2, type=bres.type)
    } else {
      x2$residuals
    }
    denom_dffits <- sqrt(c(crossprod(residx)) / x2$df.residual*h[i])
    return(num_dffits / denom_dffits)
  })
}

Funciona bien para lm:

fit <- lm(mpg ~ hp, mtcars)
dffits1(fit)
stopifnot(all.equal(dffits1(fit), dffits(fit)))

Ahora probemos betareg:

library(betareg)
data("ReadingSkills")

bfit <- betareg(accuracy ~ dyslexia + iq, data=ReadingSkills)
dffits1(bfit)
#           1           2           3           4           5           6           7 
# -0.07590185 -0.21862047 -0.03620530  0.07349169 -0.11344968 -0.39255172 -0.25739032 
#           8           9          10          11          12          13          14 
#  0.33722706  0.16606198  0.10427684  0.11949807  0.09932991  0.11545263  0.09889406 
#          15          16          17          18          19          20          21 
#  0.21732090  0.11545263 -0.34296030  0.09850239 -0.36810187  0.09824013  0.01513643 
#          22          23          24          25          26          27          28 
#  0.18635669 -0.31192106 -0.39038732  0.09862045 -0.10859676  0.04362528 -0.28811277 
#          29          30          31          32          33          34          35 
#  0.07951977  0.02734462 -0.08419156 -0.38471945 -0.43879762  0.28583882 -0.12650591 
#          36          37          38          39          40          41          42 
# -0.12072976 -0.01701615  0.38653773 -0.06440176  0.15768684  0.05629040  0.12134228 
#          43          44 
#  0.13347935  0.19670715 

No se ve mal.

Notas:

  • ¡Incluso si esto funciona en código, debe verificar si cumple con sus requisitos estadísticos!
  • He usado suppressWarningsen líneas 5:6de dffits1. predict(bfit, ReadingSkills)cae de contrastsalguna manera, mientras predict(bfit)que no (prácticamente debería ser lo mismo). Sin embargo, los resultados son idénticos: all.equal(predict(bfit, ReadingSkills), predict(bfit)), por lo que ignorar las advertencias es seguro.