# ------------------------------------------------------------------------------ # -------- Übung zur Vorlesung Generalisierte Regressionsmodelle --------------- # -------- R Code / Blatt 10 --------------- # ------------------------------------------------------------------------------ # -------- Aufgabe 1 ---------------------------------------------------------- # (a) # Daten einlesen knie <- read.table("knie.txt", header=T) attach(knie) str(knie) # Kodierung als ordinale Variable PAIN <- as.ordered(PAIN) # hier nicht unbedingt nötig, da bereits 0/1-kodiert TH <- as.factor(TH) GEN <- as.factor(GEN) # Einzentrieren um 30 (entspricht in etwa dem arithmetischen Mittel) mean(AGE) AGE <- AGE-30 # (d) library(MASS) #help(polr) # Für das kumulative Logit-Modell mit polr: # P(y<=r|x) = exp(\theta{0r}-x'\gamma)/(1+exp(\theta{0r}-x'\gamma)) # bzw. log(P(y<=r|x)/P(y>k|x)) = \theta{0r}-x'\gamma. # ACHTUNG: # In der Vorlesung, bzw. bei Tutz (2000) wird statt # log(P(y<=r|x)/P(y>k|x)) = \theta{0r}-x'\gamma nämlich # log(P(y<=r|x)/P(y>k|x)) = \theta{0r}+x'\gamma verwendet! # Die Definition von polr erleichtert die Interpretation: # positiver Koeffizient -> "positiver" Zusammenhang (in dem Sinne: die Chance auf eine HÖHERE Kategorie steigt) # Modell mit Haupteffekten (falls im Argument Hess=T, wird die beob. Fisher-Matrix gespeichert) ModelAge <- polr(PAIN ~ TH + GEN + AGE) summary(ModelAge) # Kumulatives Logit-Modell: Parameterschätzer für Einflussgrößen sind über alle # Response-Kategorien gleich; man kann Daten in disjunkte Untermengen {1,...,r} # (Y_r=1) und {r+1,...k} (Y_r=0) einteilen und diese (mit zugehörigem Intercept # \theta_{0r}) wie binäre Regressionsmodelle für die abhängige Variable Y_r # interpretieren. exp(-ModelAge$coef[1]) # Interpretation: beim Medikament ist die Chance für eine niedrigere Responsekategorie, # d.h. hier für weniger starke Schmerzen, um den Faktor 2.8 größer als beim Placebo, # bei festgehaltenen übrigen Kovariablen. # Konfidenzintervalle für Parameterschätzer: confint(ModelAge) # Man erkennt: Der Therapie-Effekt ist signifikant von 0 verschieden (negativ), # der Geschlechts-Effekt sowie der Alters-Effekt (linear) offenbar nicht. # (Das gilt wegen Dualität von Konfidenzintervallen und Hypothesentests: die im KI # enthaltenen Werte sind genau diejenigen, die beim Test nicht signifikant verworfen # werden würden.) # (f) # nun: Alter wird quadratisch mit aufgenommen AGE2 <- AGE^2 ModelAge2 <- update(ModelAge, . ~ . + AGE2) summary(ModelAge2) #Konfindenzintervalle: confint(ModelAge2) # Der quadratische Effekt des Alters ist signifikant von 0 verschieden; offenbar # hat das Alter doch einen Einfluss. Er ist nur nicht linear. # Mit einem Likelihood-Ratio Test vergleicht man, ob das Modell mit quadratischem Alter # zum Haupteffekt-Modell reduziert werden kann. anova(ModelAge,ModelAge2) # Signifikanter Unterschied (p-Wert = 0.0149), keine Reduktion möglich. Auch ein # Vergleich des AICs spricht für das Modell mit quadratischem Alter. AIC(ModelAge) AIC(ModelAge2) # Interpretation des Effekts von Alter auf Basis der Koeffizienten schwierig, # da quadratisch. Daher Zeichnung der entsprechenden Funktion alter <- function(a,z,b1,b2) { b1*(a-z) + b2*(a-z)^2 } # Beobachtete Werte des (wahren) Alters summary(knie$AGE) a <- seq(18,52,by=0.1) # Darstellung auf Basis des nicht zentrierten Alters plot(a,alter(a,30,ModelAge2$coef[3],ModelAge2$coef[4]),type="l", xlab = "Alter", ylab = "Effekt auf Schmerzempfinden") # Nach Aussage des Modells sind die (subjektiven) Schmerzen für eine # Person knapp unter 30 am schlimmsten. # Im folgenden überprüfen wir mit Hilfe der Devianzanalyse, ob das Modell # auf andere Art reduziert werden kann. # ModelAge2 ohne Geschlecht ModelAge2_gen <- update(ModelAge2, . ~ . - GEN) anova(ModelAge2,ModelAge2_gen) # keine Signifikanz (p-Wert = 0.824), also kann GEN weggelassen werden. # ModelAge2 ohne Therapie: ModelAge2_th <- update(ModelAge2, . ~ . - TH) anova(ModelAge2,ModelAge2_th) # keine Reduktion möglich (p-Wert = 0.0032). # ohne AGE: Model_age <- update(ModelAge2, . ~ . - AGE - AGE2) anova(ModelAge2,Model_age) # Wie zu erwarten war (vgl. oben): auch hier keine Reduktion möglich # (p-Wert = 0.0173) # Demnach entscheidet man sich für ein Modell mit Kovariablen TH, AGE und # AGE^2 (eine Entscheidung auf Basis des AIC würde hier zum selben Modell führen). # (g) # Verwenden des Sprays newdataTh <- data.frame(TH=factor(1), GEN=factor(0), AGE=22-30, AGE2=(22-30)^2) # Placebo newdataPl <- data.frame(TH=factor(0), GEN=factor(0), AGE=22-30, AGE2=(22-30)^2) # Vorhersage für ModelAge2 ohne GEN (wurde in (c) gewählt): predict(ModelAge2_gen,newdata=newdataTh,type="p") predict(ModelAge2_gen,newdata=newdataPl,type="p") # Die Wahrscheinlichkeit für geringeres Schmerzempfinden ist bei Verwendung des # Sprays deutlich !höher! als bei Verabreichung des Placebos. # Zum Vergleich: Vorhersage für ModelAge2 (mit GEN): predict(ModelAge2,newdata=newdataTh,type="p") predict(ModelAge2,newdata=newdataPl,type="p") # Ignorieren des Geschlechts verändert die vorhergesagten Wahrscheinlichkeiten # kaum. # Die Funktion lrm() aus der library Design kann ebenfalls zum Fit eines # kumulativen Logit-Modells verwendet werden. Achtung: dort ist die # Modellgleichung anders definiert: log(P(y>=r|x)/P(y=r,x) zu fitten, # verwende man als family sratio. # Wir verwenden zunächst den komplementären loglog Link modelvglm2b <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="cloglog",parallel=T)) summary(modelvglm2b) # Vergleich mit kumulativem Modell summary(modelvglm1b) # Die \gamma Koeffizienten sind gleich (vgl. Blatt 9) # Zugriff auf Koeffizienten mit der Function coef(modelvglm2b) # oder auf den entsprechenden Slot mittels @ modelvglm2b@coefficients # Die entsprechenden \theta des kumulativen Modells sind also theta <- modelvglm1b@coefficients[1:3] # Wir berechnen log(exp(theta)[2:3] - exp(theta)[1:2]) # und vergleichen mit modelvglm2b@coefficients[2:3] # Wie auf Blatt 9 gezeigt sind diese (bis auf numerische Ungenauigkeiten) gleich # Als nächstes berechnen wir ein sequentielles Logit- und Probit-Modell modelvglm2a <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="logit",parallel=T)) summary(modelvglm2a) modelvglm2c <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="probit",parallel=T)) summary(modelvglm2c) # Nach der Devianz zu urteilen bieten das komplementäre loglog-Modell # und das sequentielle Logit-Modell die beste Anpassung an die Daten. # (e) # Höhere Flexibilität ergibt sich, wenn man auch über Kategorien variierende # \gamma Vektoren zulässt # Wir betrachten zum Beispiel ein sequentielles Logit-Modell mit veränderlichen # Koeffizienten für Therapie modelvglm3a <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="logit",parallel=F ~ TH)) summary(modelvglm3a) # Modellvergleich über die Devianz 1 - pchisq(deviance(modelvglm2a) - deviance(modelvglm3a), df=2) # Der Vergleich legt veränderliche Koeffizienten nahe. # auch Alterseffekt nicht restringiert modelvglm4a <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="logit",parallel=F)) summary(modelvglm4a) # Modellvergleich 1 - pchisq(deviance(modelvglm3a) - deviance(modelvglm4a), df=4) # Das restringierte Modell sollte beibehalten werden. # Das Ergebnis für das komplementäre loglog-Modell ist ähnlich modelvglm3b <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="cloglog",parallel=F ~ TH)) summary(modelvglm3b) 1 - pchisq(deviance(modelvglm2b) - deviance(modelvglm3b), df=2) modelvglm4b <- vglm(PAIN ~ TH + AGE + AGE2, family=sratio(link="cloglog",parallel=F)) summary(modelvglm4b) 1 - pchisq(deviance(modelvglm3b) - deviance(modelvglm4b), df=4) # Auch das kumulative Logit-Modell kann mit kategorie-spezifischem # Therapie-Effekt stark verbessert werden. modelvglm5a <- vglm(PAIN ~ TH + AGE + AGE2, family=cumulative(link="logit",parallel=F ~ TH)) summary(modelvglm5a) summary(modelvglm1a) 1 - pchisq(deviance(modelvglm1a) - deviance(modelvglm5a), df=2)