DFFIT per la regressione beta
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
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 righe5:6didffits1.predict(bfit, ReadingSkills)lascia cadere incontrastsqualche modo, mentrepredict(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.