library(ggplot2) library(dplyr) source(file.path(pfadu, "sig.fn.R")) ovokal = read.table(file.path(pfadu, "ovokal.txt")) pvp = read.table(file.path(pfadu, "pvp.txt")) sz = read.table(file.path(pfadu, "sz.txt")) ########################################################################### ###### 0. Einleitung: Logistische Regression ########################################################################### # Mit der logistischen Regression wird geprueft, ob # Proportionen (abhaengige Variable) von einem (oder mehreren) unabhaengigen # Faktoren beeinflusst werden. # # Die abhaengige Variable ist immer kategorial und immer binaer. # Die unabhaengige Variable # kann numerisch oder kategorial (auch mehrstufig) sein. # # Beispiele: # 1. Inwiefern wird die Vokalieriung # von einem final /l/ im Englischen (feel vs. 'feeu') # vom Dialekt beeinflusst? # Abhaengige Variable: Vokalisiert (kategorial, 2 Stufen: ja, nein) # Unabhaengige Variable: Dialekt (kategorial: 2 oder mehrere Stufen) # # 2. Wird 'passt' in Augsburg im Vgl. zu Muenchen eher mit /ʃ/ produziert? # Abhaengige Variable: Frikativ (kategorial, 2 Stufen: /s/, /ʃ/) # Unabhaengige Variable: Dialekt (kategorial, 2 Stufen: Augsburg, Muenchen) # # 3. Ein offener Vokal in /lVm/ wird mit unterschiedlichen Dauern synthetisiert. # Nimmt die Wahrnehmung von 'lahm' vs. 'Lamm' mit zunehmender Dauer zu? # Abhaengige Variable: Vokal (kategorial, 2 Stufen: /a/, /a:/) # Unabhaengige Variable: Dauer (kontinuierlich) # In der linearen (least-squares) Regression wurde geprueft, inwiefern eine lineare Beziehung zwischen zwei Variablen, x und y vorliegt. # Um dies zu tun, wurde eine Regressionslinie mit Steigung m und Intercept k an die Stichproben angepasst. # yhut in: # # yhut = m x + k # # sind dann die eingeschaetzen Werte, die auf der Regressionslinie liegen # # Wir koennen aber eine solche Linie nicht an Proportionen anpassen, # weil Proportionen zwischen 0 und 1 begrenzt sind (und die lineare # Regression erwartet Werte zwischen ±unendlich). # Stattdessen wird eine gerade Linie an sogenannte Log(Odd)s angepasst # # log(Odds) = m x + k # # # Odds: P/Q. (Odds = Gewinnchancen). # P ist 'Erfolg' - z.B. in Bsp. 3 Haeufigkeit der 'lahm'-Urteile # Q ist 'Misserfolg':Haeufigkeit der 'lamm'-Urteile # log(Odds) ist dann einfach log(P/Q) # log(Odds) haben Werte zwischen -∞ und +∞ (zwischen - und + unendlich) ########################################################################### ###### 1. Erfolg (P), Misserfolg (Q), Log-Odds log(P/Q), und Proportionen (p) ########################################################################### head(ovokal) # Zwischen 1950 und 2005 sollen Woerter wie 'lost' in # einer aristokratischen Form der Standardaussprache von England # immer weniger mit einem hohen Vokal /lo:st/ # und zunehmend mit einem tiefen Vokal /lɔst/ produziert. # Ist dies Der Fall = # Wird der Vokal (hoch vs. tief) = abhaengige Variable # vom Jahr (1950... 2005) = unabhaengige numerische Variable beeinflusst? ################################ P ist 'Erfolg' # Wir nehmen den zweiten Wert von levels(ovokal$Vokal) # als P P = ovokal$Vokal == "tief" ################################ Q ist 'Misserfolg' (= nicht Erfolg) Q = !P # P, Q in den Data-Frame einbinden ovokal = cbind(ovokal, P, Q) ################################# Die Summe von P, Q pro Jahr berechnen ### hierbei Vpn ignorieren (=ueber alle Antworten aufsummieren) ovokal.m = ovokal%>% group_by(Jahr)%>% summarise(P=sum(P),Q=sum(Q)) ovokal.m # z.B. Reihe 1 bedeutet: es gab 5 Mal Erfolg (tief) und 30 Mal # Misserfolg (hoch) in 1950 # Nun die Proportionen mit P/(P+Q) berechnen; # Proportionen berechnen = die Proportion vom Erfolg (0 <= p <= 1) # neue Spalte p mit mutate() erstellen: ovokal.m = ovokal.m %>% mutate(p = P/(P+Q)) ovokal.m # Log-odds ("lodd") berechnen: log(P/Q) ovokal.m = ovokal.m %>% mutate(lodd = log(P/Q)) ########################################################################### ###### 2. Die logistische Regressionslinie ########################################################################### # Abbildung im Raum Log-Odds x unabhaengige Variable plot(lodd ~ Jahr, data = ovokal.m, ylab="Log-Odds") # Eine Regressionslinie in diesem Raum wird mit # glm(family=binomial) berechnet. # glm: generalised linear model # Entweder Anwendung auf den urspruenglichen Data-Frame head(ovokal) lreg = glm(Vokal ~ Jahr, family=binomial, data = ovokal) # Oder Anwendung auf P, Q in dem gemittelten Data-Frame lreg = glm(cbind(P, Q) ~ Jahr, family=binomial, data = ovokal.m) # Das Intercept und die Steigung sind hier coef(lreg) # Und diese kann auf die Daten im Log-Odds Raum ueberlagert werden abline(lreg) ########################################################################### ###### 3. Die Pruefstatistik ########################################################################### # Die Wahrscheinlichkeit, dass die Steigung von Null abweicht # (weil wenn die Steigung Null waere, dann waeren alle log-odds gleich # und wir haetten eine gerade linie) anova(lreg, test="Chisq") # Df Deviance Resid. Df Resid. Dev Pr(>Chi) #NULL 5 69.363 #Jahr 1 61.121 4 8.242 5.367e-15 *** # Jahr hatte einen signifikanten Einfluss auf die # Proportion von lost mit tiefem/hohem Vokal # (X^2[1] = 61.1, p < 0.001). ########################################################################### ###### 4. Die Sigmoid Funktion und der Umkipppunkt ########################################################################### # Die Regressionslinie wird berechnet im Raum # Log-Odds x unabhaengige Variable (Jahr). # Dieselben (k, m) Werte koennen verwendet werden, um # statt Log-Odds auf der y-Achse Proportionen abzubilden. # In dem Fall wandelt sich die gerade Linie in eine # sogenannte Sigmoid-Funktion um. # Eine Sigmoid-Funktion. Die mathematische Formel dafuer: # f(x) = e^(mx + k)/(1 + e^(mx + k)) # e ist die Exponentialfunktion, m und k sind die Steigung und Intercept # Umsetzung in R # sig() # m ist die Steigung: je groeßer m, umso steiler kippt # die S-Kurve um # Steigungen von 1, .5, , .25 (schwarz, rot, gruen) sig(m = c(1, .5, .25)) # Daher: wenn die Steigung 0 (Null) ist, bekommt # man eine gerade Linie, um den 0.5 Wert (wenn k = 0) sig(m = 0) # Hoehere/tiefere k-Werte verschieben die Linie um den 0.5 Wert sig(k = c(0, 1, -1), m = 0) # Der Umkipppunkt is der Punkt, zu dem (i) der Sigmoid # am steilsten ist. An diesem Punkt ist die Proportion # (auf der vertikalen Achse) immer 0.5. # Den Umkipppunkt bekommt man mit -k/m sig(k = 4, m = .8) # Hier ist m = 0.8, k = 4 # Daher u = -k/m = -4/.8 = -5 #vertikale Linie: abline(v = -5, lty=2) #horizontale Linie: abline(h = .5) ########################################################################### ###### 5. Proportionen abbilden ########################################################################### plot(p ~ Jahr, data = ovokal.m, ylab = "Proportion 'tief'") # von vorher lreg = glm(Vokal ~ Jahr, family=binomial, data = ovokal) # Intercept lreg.k = coef(lreg)[1] # Steigung lreg.m = coef(lreg)[2] # Angepasste Sigmoid sig(lreg.k, lreg.m, add=T) # verifizieren, dass es wirklich ein Sigmoid ist! plot(p ~ Jahr, data = ovokal.m, xlim = c(1920, 2020), ylim = c(0, 1), ylab = "Proportion 'tief'") sig(lreg.k, lreg.m, add=T) # Der Umkipppunkt (das Jahr, zu dem sich laut dem Modell # die Entscheidungen von hoch auf tief umkippt) -lreg.k/lreg.m # 1965.717 abline(v = -lreg.k/lreg.m, lty = 2) abline(h = .5, lty=2) ########################################################################### ###### 6. Umkipppunkte in einem synthetischen Kontinuum ########################################################################### # Ein 11-stufiges Kontinuum wurde synthesisert zwischen /pUp/ und /pYp/. # Die Stimuli wurden 10 Mal einem Hoerer einzeln praesentiert. # Der Hoerer musste pro Stimulus entscheiden: PUPP oder PUEPP? # Zu welchem F2-Wert kommt der Umkipppunkt vor? # (= zu welchem F2-Wert kippt die Entscheidung um von PUPP auf PUEPP?) # P, Q, Proportionen, logodds berechnen levels(pvp$Urteil) # Erfolg P = pvp$Urteil == "Y" # Misserfolg Q = !P # in den Data-Frame einbinden pvp = cbind(pvp, P, Q) # summieren pro F2-Wert pvp.m = pvp %>% group_by(F2) %>% summarise(P = sum(P),Q = sum(Q)) # Proportionen pvp.m = pvp.m %>% mutate(p = P/(P+Q)) plot(p ~ F2, data = pvp.m, ylab = "Proportion /Y/-Urteile") # (k,m) der Sigmoid berechnen pvp.glm = glm(Urteil ~ F2, family=binomial, data = pvp) # oder pvp.glm = glm(cbind(P, Q) ~ F2, family=binomial, data = pvp.m) # Koeffiziente pvp.k = coef(pvp.glm)[1] pvp.m = coef(pvp.glm)[2] sig(pvp.k, pvp.m, add=T) # Umkipppunkte u = -pvp.k/pvp.m abline(v = u, lty=2) # Die Wahrscheinlichkeit, dass die Urteile # durch F2-aenderungen beeinflusst werden: anova(pvp.glm, test = "Chisq") # Df Deviance Resid. Df Resid. Dev Pr(>Chi) #NULL 10 111.59 #F2 1 108.96 9 2.63 < 2.2e-16 *** # Die Proportion von pUp/pYp-Antworten wurde signifikant # von F2 beeinflusst (X^2[1] = 109.0, p < 0.001). ########################################################################### ###### 7. Der unabhaengige Faktor ist kategorial ########################################################################### # Die logistische Regression kann auf eine aehnliche Weise # verwendet werden, wenn der unabhaengige Faktor kategorial ist. # (Der wesentliche Unterschied: man braucht nicht einen Sigmoid abzubilden, # und es wird kein Umkipppunkt berechnet) head(sz) # 20 Vpn. 9 aus Bayern, 11 aus Schleswig-Holstein produzierten 'Sonne'. # Der initiale Frikativ wurde als /z/ oder /s/ wahrgenommen. # Wird die Stimmhaftigkeit vom Dialekt beeinflusst? # Die Abbildung: beide Variablen sind kategorial, daher geom_bar() ggplot(sz) + aes(fill = Frikativ, x = Dialekt) + geom_bar(position="fill") # Test sz.glm = glm(Frikativ ~ Dialekt, family=binomial, data = sz) anova(sz.glm, test = "Chisq") # Df Deviance Resid. Df Resid. Dev Pr(>Chi) #NULL 19 27.726 #Dialekt 1 5.3002 18 22.426 0.02132 * # Die s/z Verteilung wurde signifikant vom Dialekt # beeinflusst (X^2[1] = 5.3, p < 0.05)