Logistic Regression
Jonathan Harrington
Befehle: logistic.txt
1. Logistic Regression: allgemeine Einführung
Dalgaard, P. (2002) Introductory Statistics with R.
Insbesondere Kap. 11
D. Cook, P. Dixon, W. M. Duckworth, M.S. Kaiser, K. Koehler, W.
Q. Meeker and W. R. Stephenson. Binary Response and Logistic Regression Analysis.
http://www.faculty.sbc.edu/bkirk/Biostatistics/course%20documents%20for
%202006/Logistic%20Regression%20Analysis.doc
Literatur
Baayen, R.H. Analyzing Linguistic Data: A practical introduction to Statistics. S. 213-234
Johnson, Keith (in press). Quantitative Methods in Linguistics.
Blackwell. Kapitel 5.
Verzani, J. (2005). Using R for Introductory Statistics
(Ebook ueber die LMU UB). Kapitel 12
1. Logistic Regression: allgemeine Einführung
Mit logistic Regression wird eine Regressionslinie an Proportionen angepasst.
Aus verschiedenen Gründen kann jedoch die lineare (least-squares) Regression nicht auf Proportionen angewandt werden.
Vor allem liegen Proportionen zwischen 0 und 1 während lineare Regression keine solchen Grenzen kennt (und
daher könnte ein lineares Regressionsmodell Proportionen unter 0 oder über 1 vorhersagen).
Außerdem wird in der linearen Regression eine konstante
Varianz angenommen; jedoch kann bewiesen werden,
dass je höher der Proportionsdurchschnitt, umso größer
die Varianz.
Diese (und andere) Probleme können überwunden werden:
1. wenn log-odds statt Proportionen modelliert werden
2. Durch Einsetzung von 'maximum likelihood' anstatt 'least squares'.
1. Logistic Regression: allgemeine Einführung
y = mx + b
logodds(y) = mx + b
Least-squares Regression Logistic Regression
Es wird nicht angenommen, dass die Werte Stichproben aus einer Normalverteilung sind.
Ein Vorteil von logistic Regression:
Einige Daten
high low 1950 30 5 1960 18 21 1971 15 26 1980 13 20 1993 4 32 2005 2 34
In 1950 produzierten 30 Sprecher /lo:st/ und 5 /lɔst/.
jahr = as.numeric(rownames(lost)) jahr = jahr - 1950
lost
Log-odds
p q=1-p Odds = p/q Log-Odds = log(p/q)
0.8 0.2 log(4) = 1.39
0.5 0.5 1 0
4
bedeutet 4:1 (wie im Pferderennen). Die Wahrscheinlichkeit vom Erfolg (p) ist 4 Mal so groß wie Scheitern (q)
p: Proportion 'Erfolg'.
lo:st lɔst
32 8
n 40
p 0.8 (prop. lOst)
(prop. lo:st)
0.0 0.2 0.4 0.6 0.8 1.0
-4-2024
Die Logit-Funktion
p (Erfolg)
Log-odds log(p/(1-p))
Log-odds also log (p/q) als Funktion von p Log-odds haben Werte zwischen ±∞
Log-odds
2. Anwendung der logistic Regression in R: glm()
Das Ziel: nach der Anwendung von logistic Regression geben wir einen beliebigen Jahrgang ein, und das Modell soll uns die Proportion von /lo:st/ vorhersagen. zB Eingabe 1962, Proportion (lo:st) = ?
Jahr ist daher in diesem Fall die unabhängige Variable, Proportion von /lo:st/ die abhängige Variable.
Mit logistic Regression ist die abhängige Variable immer ein kategorialer Wert von 2 Möglichkeiten:
ja oder nein, rot oder grün, 0 oder 1, weiblich oder
männlich, wach oder eingeschlafen, /lo:st/ oder
/lɔst/, Erfolg oder Scheitern, usw.
g = glm(lost ~ jahr, binomial) Ergebnis: ein Log-
Odd pro Jahr
unabhängige Variable (der Jahrgang)
bedeutet:
logistic
Regression ('binomial' weil wie in der
binomialen Verteilung wir mit 2 Werten (ja/nein,
Erfolg/Scheitern zu tun haben).
wird modelliert durch Abhängige Variable
Eine 2-spaltige Matrix: Anzahl von 'ja' und 'nein' (hier /lo:st/ und /lOst/)
lost
high low 1950 30 5 1960 18 21 1971 15 26 1980 13 20 1993 4 32 2005 2 34
Da die Ausgabe der Regression in log-odds ist,
müssen wir die Proportionen ebenfalls umwandeln, wenn wir die Regressionslinie sehen wollen.
Die vorhergesagten Werte überlagern
text(jahr, predict(g), "x", col=3)
Eine Abbildung der Daten in diesem Raum:
# Proportion von /lo:st/ berechnen p = lost[,1]/apply(lost, 1, sum)
# log-odds
lodd = log(p/(1-p))
plot(jahr, lodd, type="b")
# Regressionslinie überlagern abline(g, col=2)
0 10 20 30 40 50
-3-2-101
jahr
lodd
x
x
x
x
x
x
3. Abbildung der Regressionslinie
neuerwert = data.frame(jahr = 12)
ergebnis = predict(g, neuerwert, se.fit=T) ergebnis$fit
Vorhersage: Wert für 1962
abline(h=ergebnis$fit, lty=2, col="blue")
Wir können durch die Transformation (2) die Regressionslinie auch in einem Raum von Jahr x Proportionen abbilden.
Von Proportionen in log-odds
Von log-odds zurück in Proportionen
1 )
log( p L p
L = log(p/(1-p))
L L
e p e
1 p = exp(L)/(1+exp(L)) [1] 1.386294
[1] 0.8
Abbildung der Regression
(1)
(2)
p = 0.8
curve(exp(m*x + k)/(1+ exp(m*x+k)), xlim=c(0, 60), add=T, col=2)
# Regression überlagern
# Proportionen von /lo:st/ berechnen p = lost[,1]/apply(lost, 1, sum)
plot(jahr,p)
# Abbildung Jahr x Proportionen
Die Regression
) (
) (
1
mx kk mx
e
y e
coef(g)
(Intercept) jahr 1.10432397 -0.07026313
k = coef(g)[1]
m = coef(g)[2]
Abbildung: Jahr x Proportionen
und die vorhergesagten Werte liegen wieder auf der Linie:
vorher = predict(g)
text(jahr, exp(vorher)/(1+exp(vorher)), "x", col=3)
0 10 20 30 40 50
0.20.40.60.8
jahr
p
x
x
x
x
x
x
Abbildung Jahr x Proportionen
4. Anwendung der logistic Regression auf die Berechung einer perzeptiven Grenzen zwischen Kategorien
Experiment. Anhand der Sprachsynthese wurde ein F2- Kontinuum in 11 Schritten synthetisiert. 5 Vpn. (L1-
Englisch) mussten zu jedem Stimulus mit "I" oder "U"
antworten. Bei welchem F2-Wert liegt die Grenze
zwischen den Vokalen? Die Anzahle der Bewertungen ist hier:
u i 2311 0 5 2176 0 5 2023 0 5 1885 0 5 1770 0 5 1667 0 5 1548 2 3 1437 4 1 1351 5 0 1269 5 0
Bei F2 = 1437 Hz gab es 4 Urteil für "U", ein Urteil für "I"
ui
Ein Vektor von Proportionen p = ui[,1]/apply(ui, 1, sum)
p
2311 2176 2023 1885 1770 1667 1548 1437 1351 1269 0.0 0.0 0.0 0.0 0.0 0.0 0.4 0.8 1.0 1.0
Bei 1437 Hz waren 80% der Urteile "U" (und daher 20% "I")
Eine Abbildung von F2 als Funktion dieser Proportionen f2werte = as.numeric(rownames(ui))
plot(f2werte, p, ylab="Proportion /u/ Urteile", xlab="F2 (Hz)")
Eine logistische Regression an diese Werte anpassen
logui = glm(ui ~ f2werte, family=binomial) Die Urteile aus den F2-Werten vorhersagen
Mit der Methode auf S. 13-14 die logistische Regressionskurve überlagern
m = coef(logui)[2]
k = coef(logui)[1]
curve(exp(m*x + k)/(1+ exp(m*x+k)), xlim=c(1200, 2400), add=T, col=2)
Die Koeffiziente
Die logistische Regressionskurve
Die 50% Grenze (Umkipppunkt) Die 50% Grenze (Umkipppunkt)
= zu welchem F2-Wert, ist ein Urteil für "I" genauso wahrscheinlich wie ein Urteil für "U"?
) (
) (
1
mx kk mx
e
y e
Es kann bewiesen werden, dass dies mit
-k/m in dieser Formel
gegeben wird (in diesem Beispiel ist y die Proportion, p, und x ist F2werte)
-k/m
1516.723
abline(v=-k/m, lty=2, col="blue")
Lineare Regression: R
2oder adjusted R
2und ein F-test Logistic Regression: G
2und ein
2-test.
G
2= Null deviance – residual deviance
Für ein signifikantes Ergebnis wollen wir daher, dass Null deviance hoch und Residual deviance klein ist.
wenn dieser Wert 0 wäre, dann wären alle Proportionen in allen Jahren gleich (und die Regressionslinie wäre horizontal)
je höher dieser Wert, umso unwahrscheinlicher ist es, dass die Werte überhaupt durch die Regression
modelliert werden können.
5. Signifikanz-Test
Was ist die Wahrscheinlichkeit, dass die Proportion von
/lo:st-lɔst/ durch den Jahrgang vorhergesagt werden kann?
69.3634 - 8.2422
[1] 61.1212
G
2= Null deviance - residual deviance summary(g)
Null deviance: 69.3634 on 5 degrees of freedom
Residual deviance: 8.2422 on 4 degrees of freedom
…
Der Test mit anova() ist ob G
2signifikant von 0 abweicht:
anova(g, test="Chisq")
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 5 69.363 jahr 1 61.121 4 8.242 5.367e-15
Die Proportionen folgen einem Trend (
2(1)=61.2, p < 0.001)
g = glm(lost ~ jahr, binomial)
6. Zwei unabhängige Variablen.
Hier sind genau dieselben Daten aber zusätzlich nach männlich-weiblich aufgeteilt.
female lost 1950 1960 1971 1980 1993 2005
n n 16 9 8 8 4 1
y n 14 9 7 5 0 1
n y 0 6 10 7 10 15
y y 5 15 16 13 22 19
lost
high low 1950 30 5 1960 18 21 1971 15 26 1980 13 20 1993 4 32 2005 2 34
In 1971 waren 26 Tokens [lost] und 15 [lo:st]
von diesen 26 waren 10 von Männern und 16 von Frauen
erzeugt.
8M, 7F
(a) Gibt es einen Trend?
Also weniger [lo:st] in
späteren Jahren? (b) Ist die Proportion [lost]/[lo:st] in
M und F unterschiedlich verteilt?
female lost 1950 1960 1971 1980 1993 2005
n n 16 9 8 8 4 1
y n 14 9 7 5 0 1
n y 0 6 10 7 10 15
y y 5 15 16 13 22 19
Dies ist ein Problem der mehrfachen Logistic Regression:
logodds (lo:st) = b
0+ b
1year + b
2Geschlecht (also in diesem Fall eine Linie im 3D-Raum)
Geschlecht
Year logodds(lo:st)
Und eine gerade Linie in einem 3D-Raum
(b
0ist das Intercept,
b
1und b
2die Neigungen)
pfad = "das Verzeichnis wo ich lost2.txt gespeichert habe"
lost2 = as.matrix(read.table(paste(pfad, "lost2.txt", sep="/")))
high low 0.0 16 0 10.0 9 6 21.0 8 10 30.0 8 7 43.0 4 10 55.0 1 15 0.1 14 5 10.1 9 15 21.1 7 16 30.1 5 13 43.1 0 22 55.1 1 19
M
W
1950 1960 1971 1980 1993 2005 1950 1960 1971 1980 1993 2005
} }
high = Spalte 1 = /lo:st/
low = Spalte 2 = /lOst/
Daten-Vorbereitung
J = c(jahr, jahr)
G = c(rep(0, 6), rep(1, 6)) J
G
[1] 0 10 21 30 43 55 0 10 21 30 43 55
[1] 0 0 0 0 0 0 1 1 1 1 1 1
Zuerst eine Abbildung…
p = lost2[,1]/apply(lost2, 1, sum) interaction.plot(J, G, p)
0.00.40.8
J
mean of p
0 10 21 30 43 55
G m f
Nimmt die Proportion von
/lo:st/ in späteren Jahren ab?
(Die Unterschiede zwischen m und f ignorieren).
Ja Nein Vielleicht
Unterscheiden sich m und f in der Proportion von /lo:st/?
(Die Unterschiede in den Jahrgängen ignorieren).
Ja Nein Vielleicht
Modell berechnen…
mehrg = glm(lost2 ~ J + G, binomial)
g2 = glm(lost2 ~ J, binomial) anova(g2, test="Chisq")
Analysis of Deviance Table
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 11 89.557 year 1 61.121 10 28.436 5.367e-15
Wenn wir übrigens G weglassen, dann müssten wir
trotz der anderen Aufteilung der Daten das gleiche
Ergebnis wir vorhin bekommen:
Coefficients:
(Intercept) J Gm
1.87754 -0.07524 1.20282
Degrees of Freedom: 11 Total (i.e. Null); 9 Residual Null Deviance: 89.56
Residual Deviance: 15.61 AIC: 51.51
logodds(lo:st) = 1.87754 - 0.07524J+ 1.20282G mehrg = glm(lost2 ~ J + G, binomial)
anova(mehrg, test="Chisq")
Df Deviance Resid.Df Resid. Dev P(>|Chi|) NULL 11 89.557 J 1 61.121 10 28.436 5.367e-15 G 1 12.822 9 15.613 3.425e-04
M und F unterscheiden sich in der Proportion von lo:st/lOst,
2(1) = 12.82, p < 0.001
Die Proportion von 'lo:st' nimmt in späteren Jahren ab,
2(1) = 61.12, p < 0.001.
mehrg
Mit 2 oder mehr Variablen soll auch geprüft werden, ob sie miteinander interagieren.
Eine Interaktion zwischen den unabhängigen Variablen – in diesem Fall Geschlecht und Jahrgang – liegt vor, wenn sie eine unterschiedliche Wirkung auf die abhängige
Variable ausüben wie in 1 und 2, aber nicht in 3 und 4
7. Die Interaktion zwischen 2 Variablen
1950 2000
prop(lo:st)
1950 2000 1950 2000 1950
2000 prop(lo:st) prop(lo:st) prop(lo:st)
1 2 3 4
m
f
Wenn eine Interaktion vorliegt, dann können signifikante Ergebnisse in einer der unabhängigen Variablen nicht uneingeschränkt akzeptiert werden.
zB wenn eine Interaktion vorkommt, gibt es vielleicht eine Wirkung von Jahrgang auf die Proportion von /lo:st/ nur in Männern aber nicht in Frauen usw.
Die Interaktion zwischen 2 Variablen
0.00.20.40.60.81.0
year
mean of p
0 10 21 30 43 55
geschl m f
dies scheint aber hier nicht
der Fall zu sein.
Die Interaktion zwischen 2 unabhängigen Variablen, A und B, kann in R mit A:B geprüft werden.
Daher in diesem Fall
g = glm(lost2 ~ J + G + J:G, binomial)
Eine Abkürzung dafür (und mit genau demselben Ergebnis) g = glm(lost2 ~ J * G, binomial)
Die Interaktion zwischen 2 Variablen
anova(g, test="Chisq")
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 11 89.557
J 1 61.121 10 28.436 5.367e-15 G 1 12.822 9 15.613 3.425e-04 J:G 1 0.017 8 15.596 0.896
d.h. die Interaktion ist nicht signifikant und J:G kann
aus dem Regressionsmodell weggelassen werden.
Start: AIC= 53.49 lost2 ~ J * G
Df Deviance AIC - J:G 1 15.613 51.506
<none> 15.596 53.489 Df Deviance AIC
<none> 15.613 51.506 - G 1 28.436 62.328 - J 1 80.018 113.910
Wir bleiben also bei
Call: glm(formula = lost2 ~ J + G, family = binomial) Residual Deviance: 15.61 AIC: 51.51
library(MASS) stepAIC(g)
Dies wird auch durch stepAIC() bestätigt:
AIC wird kleiner wenn wir
J:G weglassen
Weitere Folien zum Durchlesen...
Aus dem vorigen Beispiel wird auch klar, dass ähnlich wie
2Logistic Regression angewandt werden kann, auch wenn die Gruppe nur aus 2 Ebenen besteht.
Gibt es einen signifikanten Unterschied zwischen M und F?
gmf = glm(lost2 ~ G, "binomial") anova(gmf, test="Chisq")
Df Deviance Resid. Df Resid. Dev P(>|Chi|)
NULL 11 89.557 G 1 9.539 10 80.018 0.002
M und F unterscheiden sich in der Proportion von lo:st/lOst (
2(1) = 9.5, p < 0.002).
8. Logistic Regression und zwei Ebenen
high low
0.0 16 0 10.0 9 6 21.0 8 10 30.0 8 7 43.0 4 10 55.0 1 15 0.1 14 5 10.1 9 15 21.1 7 16 30.1 5 13 43.1 0 22 55.1 1 19
m = apply(lost2[1:6,], 2, sum) f = apply(lost2[7:12,], 2, sum) mf = rbind(m, f)
rownames(mf) = c(0, 1)
colnames(mf) = c("high", "low") mf
Wir bekommen dasselbe Ergebnis wenn Logistic Regression auf die entsprechende Tabelle
angewandt wird:
lost2
l.mf = c(0,1)
gmf2 = glm(mf ~ l.mf, "binomial") anova(gmf2, test="Chisq")
high low 0 46 48 1 36 90
= (kodiert nur nach M und F)
und man bekommt dann fast das gleiche Ergebnis mit einem
2-Test, der direkt auf die Tabelle angewandt wird:
chisq.test(mf)
Pearson's Chi-squared test with Yates' continuity correction
data: mf
X-squared = 8.6985, df = 1, p-value = 0.003185
Ein 2-Test kann jedoch nicht verwendet
werden, bei einer Gruppenanzahl von > 2 …
Haben (a) Alter und (b) Geschlecht einen Einfluss auf die Proportion von /lo:st/?
Hier sind dieselben Daten aufgeteilt in 2 Altersgruppen sowie M/F
Gruppe 1 = Vokal = high/low
Gruppe 2 = Geschl = M/F (=0/1) Gruppe 3 = Alter = alt/jung
3 Gruppen jeweils 2 Ebenen
high low alt.0 43 35 alt.1 30 15 jung.0 3 13 jung.1 6 75
lost3
Zuerst eine Abbildung high low alt.0 43 35 alt.1 30 15 jung.0 3 13 jung.1 6 75
# Alter kodieren A = c(0, 0, 1, 1)
# Geschlecht kodieren G = c(0, 1, 0, 1)
prop = lost3[,1]/apply(lost3, 1, sum)
interaction.plot(A, G, prop)
im Geschlecht?
ja nein
Signifikanter Einfluss auf lo:st/lOst?
im Alter? vielleicht
ja nein vielleicht Interaktion zwischen A und G? ja nein vielleicht
0.10.20.30.40.50.6
A
mean of prop
0 1
G 01
g = glm(lost3 ~ A * G, binomial) anova(g, test="Chisq")
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 3 67.758 A 1 64.452 2 3.306 9.893e-16 G 1 0.398 1 2.908 0.528 A:G 1 2.908 0 -4.441e-16 0.088