DFFITs für die Beta-Regression
Ich versuche, DFFITS für GLM zu berechnen, wobei die Antworten einer Beta-Verteilung folgen. Mit dem betaregR-Paket. Aber ich denke, dieses Paket wird nicht unterstützt, influence.measures()weil durch die Verwendung von 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
es gibt nach
Fehler in if (Modell $ rank == 0) {: Argument hat die Länge Null
Ich bin ein Neuling in R. Ich weiß nicht, wie ich dieses Problem lösen soll. Bitte helfen Sie, dieses Problem zu lösen, oder geben Sie mir einige Tipps über den R-Code, wie DFFITs manuell berechnet werden. Grüße
Antworten
dffitssind nicht für "betareg"Objekte implementiert , aber Sie können versuchen, sie manuell zu berechnen.
Entsprechend dieser Frage / Antwort zum Stapelüberlauf könnten wir diese Funktion schreiben:
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)
})
}
Es funktioniert gut für lm:
fit <- lm(mpg ~ hp, mtcars)
dffits1(fit)
stopifnot(all.equal(dffits1(fit), dffits(fit)))
Versuchen wir jetzt 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
Sieht nicht schlecht aus.
Anmerkungen:
- Auch wenn dies im Code funktioniert, sollten Sie überprüfen, ob es Ihren statistischen Anforderungen entspricht!
- Ich habe
suppressWarningsin Zeilen5:6von verwendetdffits1.predict(bfit, ReadingSkills)lässt dascontrastsirgendwie fallen, wohingegenpredict(bfit)nicht (sollte praktisch gleich sein). Die Ergebnisse sind jedoch identisch:all.equal(predict(bfit, ReadingSkills), predict(bfit))Daher ist das Ignorieren der Warnungen sicher.