In general linear models verwendet man kein R^2. Das VGAM-Paket gibt beispielsweise gar kein R^2 bei Modellen aus: Ein guter Hinweis, dass das vielleicht auch keine gute Idee ist!
Wenn ich nun aber dennoch einen Wert angeben will, wie "gut" mein Modell ist? Dann gibt es diverse "Pseudo-R^2".
Eine einfache Überlegung ist, eine Korrelation zwischen den echten und den vorhergesagten Werten zu brechnen. Das Cox-Shell-pseudo-R^2 hingegen berechnet einen Wert auf Basis der Loglikelihood des Modells und der Loglikelihood eines Modells mit derselben Verteilung, aber ohne Prediktoren (intercept only model).
Wie schneiden diese beiden Werte ab?
Zunächst baue ich mir einen Testdatensatz aus normalverteilten Zufallszahlen als abhängige Variable Y, einem Fehlerwert e, der so berechnet ist, dass er im Prinzip 10 Prozent von Y abweicht (hierfür verwende ich den mvrnorm-Befehl aus dem MASS-Paket), und einer unabhängigen Variable X = Y + e.
cfit sollte identisch mit c sein.
Wie sieht es jetzt mit einer Gumbel-Verteilung aus? Ich ändere die Y-Werte, nicht aber den Fehlerterm e.
Wenn ich nun aber dennoch einen Wert angeben will, wie "gut" mein Modell ist? Dann gibt es diverse "Pseudo-R^2".
Eine einfache Überlegung ist, eine Korrelation zwischen den echten und den vorhergesagten Werten zu brechnen. Das Cox-Shell-pseudo-R^2 hingegen berechnet einen Wert auf Basis der Loglikelihood des Modells und der Loglikelihood eines Modells mit derselben Verteilung, aber ohne Prediktoren (intercept only model).
Wie schneiden diese beiden Werte ab?
Zunächst baue ich mir einen Testdatensatz aus normalverteilten Zufallszahlen als abhängige Variable Y, einem Fehlerwert e, der so berechnet ist, dass er im Prinzip 10 Prozent von Y abweicht (hierfür verwende ich den mvrnorm-Befehl aus dem MASS-Paket), und einer unabhängigen Variable X = Y + e.
Normalverteilte Werte
Da ich gefittete Werte mit den echten Werten vergleichen will, nehme ich gleich schon den vglm-Befehl aus dem VGMA-Paket, da ich später eine Gumbel-Verteilung testen will.# requires(MASS, VGMA)cn sollte jetzt einen Wert von ca. 0.95.
library(MASS)
library(VGMA)
Y = rnorm(10000, 0, 1)
e = mvrnorm(10000, 0, .1)
X = Y+e
cn = cor(Y,X)
fit = vglm(Y~X, family=gaussianff)
P = predict(fit)
cfit = cor(Y,P)
cfit sollte identisch mit c sein.
Gumbelverteilte Werte
Wie sieht es jetzt mit einer Gumbel-Verteilung aus? Ich ändere die Y-Werte, nicht aber den Fehlerterm e. Y = rgumbel(10000, 0, 1)Anschließend berechne ich noch das Cox-Shell-pseudo-R^2:
X = Y+e
cg = cor(Y,X)
fitg = vglm(Y~X, family=gumbel)
P = predict(fitg)
cfitg = cor(Y,P)[1]
# Cox and Shell pseudo R^2Während der Wert für die bloße Korrelation über den Werten der Normalverteilung liegt (eben aufgrund der Gumbel-Verteilung), obwohl sich der Fehlerterm nicht geändert hat, liegt das Cox-Shell-pseudo-R^2 deutlich drunter. Pseudo-R^2 ist also nicht gleich Pseudo-R^2...
inter = vglm(Y~1, family=gumbel)
N <- 10000
m.ll <- logLik(fitg)[1]
n.ll <- logLik(inter)[1]
cs.R2 <- 1 - exp(-2/N ∗ (m.ll - n.ll))
cs.R2
print(c(cn, cfit, cg, cfitg, cs.R2))
> (c(cn, cfit, cg, cfitg, cs.R2)) [1] 0.9555202 0.9555202 0.9714456 0.9714456 0.9213320
Liebe Studierende,
AntwortenLöschennur falls sich jemand wundert: Dieser Post hat nichts mit dem Seminar zu tun. Wer nur die seminarrelevanten Einträge sehen will, sollte den Link "Datamining für SozialwissenschaftlerInnen" anklicken oder das Label "Seminar".
Grüße
SH