DFFIT per la regressione beta

Aug 26 2020

Sto cercando di calcolare DFFITS per GLM, dove le risposte seguono una distribuzione beta. Utilizzando betaregil pacchetto R. Ma penso che questo pacchetto non supporti influence.measures()perché utilizzando 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

cede

Errore in if (model$rank == 0) { : l'argomento è di lunghezza zero

Sono un principiante in R. Non so come risolvere questo problema. Gentilmente aiutami a risolvere questo problema o dammi alcuni suggerimenti tramite il codice R su come calcolare manualmente i DFFIT. Saluti

Risposte

2 jay.sf Aug 26 2020 at 15:12

dffitsnon sono implementati per "betareg"gli oggetti, ma potresti provare a calcolarli manualmente.

Secondo questo Stack Overflow Q/A potremmo scrivere questa funzione:

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)
  })
}

Funziona bene per lm:

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

Ora proviamo 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 

Non sembra male.

Appunti:

  • Anche se funziona nel codice, dovresti verificare se soddisfa i tuoi requisiti statistici!
  • Ho usato suppressWarningsnelle righe 5:6di dffits1. predict(bfit, ReadingSkills)lascia cadere in contrastsqualche modo, mentre predict(bfit)non (dovrebbe praticamente essere lo stesso). Tuttavia i risultati sono identici: all.equal(predict(bfit, ReadingSkills), predict(bfit)), quindi ignorando gli avvertimenti essere prudenti.