library(ggplot2) source(file.path(pfadu, "proben.R")) form = read.table(file.path(pfadu, "bet.txt")) e.df = read.table(file.path(pfadu, "e.txt")) ################################################################## ### 1. SE (Standard Error) und Konfidenzintervall ################################################################## # (Siehe Vorlesung, Normalverteilung) ####################### SE # Der SE ist die Populationsstandardabweichung von Mittelwerten. # z.B. Wir werfen 12 Wuerfel zusammen, und berechnen den Mittelwert # der Zahlen. # Der SE setzt sich zusammen aus σ/√k # σ ist die Populationsstandardabweichung = sd(1:6) * sqrt(5/6) # k ist die Anzahl der Wuerfel = 12 SE = sd(1:6) * sqrt(5/6) / sqrt(12) # Bedeutung: wenn wir 12 Wuerfel zusammenwerfen, # davon den Mittelwert berechnen, # diesen Vorgang unendlich viel Mal wiederholen, # sodass wir unendlich viele Mittelwerte haetten, # dann davon die Standardabweichung berechnen, # dann waere diese Standardabweichung genau # SE = sd(1:6) * sqrt(5/6) / sqrt(12) # Wir muessten ziemlich nah an diesen SE mit z.B. # 50000 Mittelwerten kommen: o = proben(k=12, N = 50000) sd(o) ######################## Konfidenzintervall # Wir wollen zwei Werte a und b auf eine solche Weise berechnen, # sodass der Mittelwert zwischen a und b # mit einer Wahrscheinlichkeit von 0.95 liegt. # Hier benoetigen wir den Populations-Mittelwert: # mu = mean(1:6) # und SE SE = sd(1:6) * sqrt(5/6) / sqrt(12) # a und b a = qnorm(0.025, mu, SE) b = qnorm(0.975, mu, SE) # d.h. wenn wir 12 Wuerfel zusammen werfen, # davon den Mittelwert berechnen, # dann faellt der Mittelwert zwischen # (a) 2.533725 und (b) 4.466275 # mit einer Wahrscheinlichkeit von 0.95 # (= ein Mittelwert von weniger # als 2.533725 oder mehr als 4.466275 wird meistens nur # in 5/100 Faellen vorkommen). # Pruefen: 12 Wuerfel werfen, davon den Mittelwert berechnen, diesen # Vorgang 100 Mal wiederholen (daher 100 Mittelwerte) m = proben(k = 12, N = 100) # Wieviele dieser Mittelwerte sind weniger # als 2.486562 oder groeßer als 4.513438? sum(m < 2.533725 | m > 4.466275) ################################################################## ### 2. Konfidenzintervall einer Stichprobe ################################################################## # Hier sind 12 Dauerwerte von einem /a:/ Vokal d = c(119, 111, 105, 130, 133, 122, 124, 129, 95, 100, 109, 111) # Wir wollen aufgrund der Stichprobe ein Konfidenzintervall # fuer den Dauer-Mittelwert von /a:/ erstellen (= dieser Mittelwert # faellt zwischen a und b mit einer Wahrscheinlichkeit von 95%). # Zu diesem Zweck muessen: # i. (mu, SE) aufgrund # der Stichprobe wie folgt eingeschaetzt werden, angenommen dass es # sich um eine randomisiert ausgewaehlte Stichprobe handelt: # beste Einschaetzung von mu mu = mean(d) # beste Einschaetzung von SE SE = sd(d)/sqrt(12) # ii. Wenn (mu, SE) nicht aus theoretischen Ueberlegungen (wie beim Wuerfel-Spiel) berechnet werden, dann wird das Konfidenzintervall nicht mit der Normalverteilung sondern durch die t-Verteilung # mit einer gewissen Anzahl von Freiheitsgraden erstellt. # # Die t-Verteilung ist der Normalverteilung sehr aehnlich, # und naehert sich der Normalverteilung # je hoeher die Anzahl der Freiheitsgrade # z.B. t-Verteilung # mit 3, 6, 20 Freiheitsgrade auf die Normalverteilung ueberlagern: ylab = "Wahrscheinlichkeitsdichte" curve(dnorm(x, 0, 1), xlim = c(-4, 4), ylab = ylab) # t-Verteilung mit 3 df curve(dt(x, 3), add = T, col = "red") # t-Verteilung mit 6 df curve(dt(x, 6), add = T, col = "blue") # t-Verteilung mit 20 df curve(dt(x, 20), add = T, col = "green") # Die Anzahl der Freiheitsgrade = die Anzahl der Stichproben minus 1 df = 11 # Das Konfidenzintervall fuer die obige Stichprobe wird mit qt() # statt qnorm() berechnet: a = mu + SE * qt(0.025, df) b1 = mu - SE * qt(0.025, df) #oder b2 = mu + SE * qt(0.975, df) a b1 b2 # Bedeutung: aufgrund dieser Stichprobe faellt der Dauermittelwert # fuer /a/ zwischen 107.8 ms und 123.5 ms mit einer # Wahrscheinlichkeit von 0.95. ################################################################## ### 3. Konfidenzintervall fuer den Unterschied zwischen 2 Stichproben ################################################################## # Die Dauerwerte wenn 12 Sprecher (Vpn1, Vpn2... Vpn12) einen betonten /a/ produzierten sind wie folgt: bet = c(119, 111, 105, 130, 133, 122, 124, 129, 95, 100, 109, 111) # Die Dauerwerte wenn die selben 12 Sprecher (Vpn1, Vpn2... Vpn12) einen unbetonten /a/ produzierten sind: un = c(110, 95, 108, 80, 120, 110, 120, 95, 72, 83, 90, 95) # Wir wollen einen 95% Konfidenzintervall erstellen fuer den # Dauerunterschied zwichen betont und unbetont mu = mean(bet - un) SE = sd(bet - un)/sqrt(12) # 95% Konfidenzintervall mu + SE * qt(0.025, 11) mu + SE * qt(0.975, 11) # Das 95% Konfidenzintervall fuer den Unterschied zwischen den # Mittelwerten ist 8.7 ms ≤ mu ≤ 26.3 ms # Bedeutung: der Unterschied zwischen # dem Dauer-Mittelwert von einem betontem /a/ # und dem Dauer-Mittelwert von einem unbetonten /a/ # liegt zwischen # 8.7 ms und 26.3 ms mit einer Wahrscheinlichkeit von 0.95. ################################################################## ### 4. Pruefen, ob sich Mittelwerte signifikant unterscheiden ################################################################## # H0 (die Null Hypothese). # Betonung hat keinen Einfluss auf die Dauer. # Bedeutung: der Unterschied zwischen den Mittelwerten ist Null (0). # H1 (die Alternativ-Hypothese). Betonung beeinflusst die Dauer. # Bedeutung: der Unterschied zwischen den Mittelwerten weicht von Null ab. # α-Wert (alpha-Wert) # 0.05 ist der hier gewaehlte, sogenannte α-Wert (alpha-Wert) # bei dem wir H0 verwerfen (Ein Konfidenzintervall von 0.95 # bedeutet einen α-Wert von 1 - 0.95 = 0.05). # Pruefen: wenn 0 außerhalb des Konfidenzintervalls # 8.7 ms ≤ mu ≤ 26.3 ms faellt, verwerfen wir H0 # und akzepierten H1. # Schlussfolgerung fuer das obige Beispiel: # Wir verwerfen H0 und akzeptieren H1. # Berichten: # # Die Dauer wird signifikant von der Betonung beeinflusst (p < 0.05) # Bedeutung (i): # die Wahrscheinlichkeit, dass die Dauer NICHT von der # Betonung beeinflusst wird, ist weniger als 0.05 (weniger als 5%). # Bedeutung (ii): die Wahrscheinlichkeit, dass der Unterschied # zwischen dem Dauer-Mittelwert von betontem /a/ und # dem Dauer-Mittelwert von unbetontem /a/ Null sein koennte, ist # weniger als 0.05. ################################################################## ### 5. Der t-test ################################################################## # Wir bekommen dieselbe Auswertung einfacher durch den sogenannten t-test: # Manuell noch einmal: mu = mean(bet - un) SE = sd(bet - un)/sqrt(12) # 95% Konfidenzintervall mu + SE * qt(0.025, 11) mu + SE * qt(0.975, 11) # mit dem t-test # t.test(bet - un) # Identifizieren: # mu (17.5) # # Freiheitgrade (11) # # 95% Konfidenzintervall (8.739938 26.260062) # # t = 4.3969 (die t-Statistik). # Der Abstand zwischen mu und 0 (Null) in SE-Einheiten. (mu - 0)/SE # # p-value = 0.001069: # Bedeutung (i): Die Wahrscheinlichkeit dass # der Unterschied zwischen den Mittelwerten Null sein koennte. # Bedeutung (ii): Die Wahrscheinlichkeit, dass H0 zutrifft. # Bedeutung (iii): die Wahrscheinlichkeit, dass der Unterschied zwischen # den Mittelwerten außerhalb des Konfidenzintervalls faellt, auch: (1- pt(4.396914, 11)) * 2 # Ergebnis berichten. # Wir waehlen immer drei α-Werte aus: p < 0.05, p < 0.01, p < 0.001 # und waehlen den α-Wert, der am naechsten ueber dem p-Wert liegt # Hier waehlen wir p < 0.01, da der p-Wert (0.001069) ueber 0.001 aber unter 0.01 liegt. # Berichten # Entweder: # Betonung hatte einen signifikanten Einfluss auf die # Dauer (t[11] = 4.4, p < 0.01) # Oder: # Die Betonung wurde signifikant von der Dauer beeinflusst # (t[11] = 4.4, p < 0.01) # # Sollte p > 0.05 sein, dann ist das Ergebnis nicht signifikant # (Wir verwerfen nicht H0), und schreiben: # Die Betonung hatte keinen signifikanten Einfluss auf die Dauer # oder # Die Dauer wurde nicht signifikant von der Dauer beinflusst. ################################################################## ### 6. Gepaarter t-test ################################################################## # Das obige Beispiel ist ein gepaarter t-test, weil eine # Differenz - ob betont oder unbetont - pro Paar berechnet wird # (daher 12 Paare in dem obigen Beispiel). # Ein gepaarter t-test kommt in der Phonetik meistens # vor, wenn Stichproben-Paare pro Versuchsperson verglichen werden. # # Zweites Beispiel: wenn (wie ueblich) # die Werte in einem Data-Frame stecken. # 12 Versuchspersonen produzierten jeweils einen betonten # unbetonten /i/. Unterscheiden sich der betonte und unbetonte /i/ # in F2? dim(form) head(form) form # Die Frage immer umstellen (Wird y von x beeinflusst?) # Wurde F2 (abhaengige Variable) von der Betonung # (unabhaengiger Faktor mit 2 Stufen: betont/unbetont) beeinflusst? # Der Test ist gepaart: es gibt ein paar von Werten pro Versuchsperson. # 1. Boxplot # 1. Unterschiede pro Paar (hier Versuchsperson) berechnen unterschied = aggregate(F2 ~ Vpn, diff, data = form) # oder library(dplyr) unterschied = form%>% group_by(Vpn)%>% summarise(F2=diff(F2)) unterschied # 2. Boxplot der Unterschiede # am einfachsten ohne ggplot2: boxplot(unterschied$F2) # mit ggplot2 (x einfach weglassen): ggplot(unterschied) + aes(y = F2) + geom_boxplot() #oder: ggplot(unterschied) + aes(x = "betont-unbetont", y = F2) + geom_boxplot() + xlab("Unterschiede in F2") # 3. t.test: Prueft, ob der Mittelwert der Unterschiede signifikant abweicht von 0 oder nicht t.test(unterschied$F2) # F2 wurde signifikant von der Betonung beeinflusst # (t[11] = 4.4, p < 0.01). # Alternativ # y ~ x unter der Beruechsichtigung, dass x gepaarte Werte hat t.test(F2 ~ Bet, paired = T, data = form) #dazu passend eine alternative Abbildung der selben Daten: ggplot(form) + aes(x = Bet, y = F2,col = Vpn, group = Vpn) + geom_point() + geom_line() ################################################################## ### 7. Ungepaarter t-test ################################################################## # Wenn nicht Paare von Stichproben sondern # zwei Gruppen verglichen werden. # In der Phonetik: meistens wenn zwei verschiedene # Sprecher-Gruppen (maennlich/weiblich; Bayern/Hessen; englisch/deutsch) # verglichen werden. # Vorgang: wie oben aber # die Stichproben werden nicht paarweise voneinander # subtrahiert (da sie nicht gepaart sind). # # Unterscheiden sich deutsch und englisch in F2 von /e/? # = Wird F2 (abhaengige Variable) von der Sprache (unabhaengige Variable # mit 2 Stufen: englisch/deutsch) beeinflusst? head(e.df) dim(e.df) # 1. Boxplot oder Densitplot der Unterschiede ggplot(e.df) + aes(y = F2, x = Sprache) + geom_boxplot() # density plot mit geom_density() ggplot(e.df) + aes(x = F2, col = Sprache) + geom_density() # Hier pruefen wir, ob signifikante Unterschiede # zwischen den Mittelwerten der beiden Gruppen vorliegen. # (NB: nicht gepaart!) t.test(F2 ~ Sprache, data = e.df) # "F2 wurde signifikant (t[21.1] = 2.3, p < 0.05) von der Sprache beeinflusst." # (In einem ungepaarten Welch Two Sample t-Test bekommt man # immer Gleitkommazahlen als Freiheitsgrade, # da die Freiheitsgrade mittels einer Formel (der sogenannten # Welch-Satterthwaite-Formel) verändert (“approximiert”) werden. # Durch diese Änderung müssen wir uns weniger Gedanken darüber machen, # wann ein solcher Test erlaubt ist oder nicht.). # 8. Voraussetzung für den t-Test prüfen: shapiro.test() # zur Überprüfung von Verteilungen # Der shapiro.test() überprüft, ob die Verteilungen der Daten # signifikant von einer Normalverteilung abweichen. # Sollte dies der Fall sein, dürfen wir den t-Test nicht anwenden, #sondern müssen statt dessen den Wilcoxon-Vorzeichen-Rang-Test # # anwenden (wilcox.test()). #Ergibt der shapiro.test() keine signifikante Abweichung # der Verteilungen unserer Stichproben von einer Normalverteilung, # bleiben wir beim t-Test. #Beispiel für den gepaarten Fall: shapiro.test(unterschied$F2) #daher wilcox.test(unterschied$F2) #“F2 wurde signifikant von der Betonung # beeinflusst, wie ein Wilcoxon-Test zeigte (V=3,p<0.01).” #ungepaart: with(e.df, tapply(F2, Sprache, shapiro.test)) #da beide nicht signifikant --> t.test(F2 ~ Sprache, data = e.df) # “F2 wurde signifikant (t[21.1]=2.3,p<0.05) # von der Sprache beeinflusst.”