--- title: "Statistik 5" subtitle: "Übung 4" format: revealjs: chalkboard: remember-playback: true create-playback: true slide-number: true #embed-resources: true css: - webex.css - statistics_exercises_reveal.css include-after-body: - webex.js - webex-pdf-persist.js height: 1200 width: 1800 webexercises: storage-key: "Uebung_Statistik5_4" --- ```{r setup} #| include: false library("exams2forms") knitr::opts_knit$set("out.format" = NULL, "rmarkdown.pandoc.to" = NULL) ``` ## Einleitung und Ziele der Übung \ Herzlich Willkommen zur Übung Statistik 5. Bei dieser Übung wollen wir *t*-Tests für unabhängige Stichproben mit passenden Konfidenzintervallen und *U*-Tests vergleichen. ## Agenda \ - Konfidenzintervall Mittelwertsunterschiede - *t*-Test für unabhängige Stichproben - *U*-Test bzw. Wilcoxon-Rangsummentest ## 1. Konfidenzintervall Mittelwertsunterschiede - Einführung \ In der Vorlesung haben wir uns bereits die Idee des Signifikanztests angeschaut. In dieser Übung soll geschaut werden, wie sich *t*-Test, Konfidenzintervalle und *U*-Test verhalten. ## 1. Konfidenzintervall Mittelwertsunterschiede - Übung 1 \ Dafür sollen jeweils 100'000 Stichproben aus zwei Populationen (A, B) gezogen werden. Mit diesen werden wir durchweg arbeiten: - Jede Population besteht aus je zwei Subpopulationen 1 und 2, wobei für jede Teilpopulation in jeder Stichprobe $n=64$ gelten soll. - Die Populationsparameter für A1 seien: $\mu_{A1} = 100,\quad\sigma^2_{A1}=144$ - Die Populationsparameter für A2 seien: $\mu_{A2} = 100,\quad\sigma^2_{A2}=144$ - Die Populationsparameter für B1 seien: $\mu_{B1} = 100,\quad\sigma^2_{B1}=144$ - Die Populationsparameter für B2 seien: $\mu_{B2} = 106,\quad\sigma^2_{B2}=144$ Benutzen Sie den Seed 1000. Simulieren Sie die die 100'000 Stichproben je Subpopulation und speichern Sie das Ergebnis jeweils in einer Variable nach dem Schema: Population_Subpopulation, z.B. A_1 Wenn Sie alles richtig gemacht haben, sollten Sie etwa folgende Werte für die jeweils erste Person ihrer Stichprobe haben: 94.7, 110.0, 99.6, 89.4 ```{r} #| eval: false #| include: false #| echo: false Stichproben <- 100000 n <- 64 mu_A1 <- 100 mu_A2 <- 100 mu_B1 <- 100 mu_B2 <- 106 sigma2_A1 <- 144 sigma2_A2 <- 144 sigma2_B1 <- 144 sigma2_B2 <- 144 set.seed(1000) A_1 <- replicate(Stichproben,rnorm(n,mu_A1,sqrt(sigma2_A1))) A_2 <- replicate(Stichproben,rnorm(n,mu_A2,sqrt(sigma2_A2))) B_1 <- replicate(Stichproben,rnorm(n,mu_B1,sqrt(sigma2_B1))) B_2 <- replicate(Stichproben,rnorm(n,mu_B2,sqrt(sigma2_B2))) ``` ## 1. Konfidenzintervall Mittelwertsunterschiede - Übung 2 \ Berechnen sie für jede der Stichprobe und jede Subpopulation die folgenden Stichprobenkennwerte und speichern Sie diese in jeweils in Variable. Sie können dafür wieder apply verwenden. - $\overline x$ (z.B. A_1_M als Variable) - $s$ (z.B. A_1_SD als Variable) Haben Sie alles richtig gemacht, sollten Sie jetzt 8 Variablen mit jeweils 100'000 Werten haben. ```{r} #| eval: false #| include: false #| echo: false A_1_M <- apply(A_1, 2, mean) A_2_M <- apply(A_2, 2, mean) B_1_M <- apply(B_1, 2, mean) B_2_M <- apply(B_2, 2, mean) A_1_SD <- apply(A_1, 2, sd) A_2_SD <- apply(A_2, 2, sd) B_1_SD <- apply(B_1, 2, sd) B_2_SD <- apply(B_2, 2, sd) ``` Berechnen Sie jetzt folgende Grössen für den Vergleich der Subgruppen 1 und 2 wieder getrennt für die Populationen A und B. Sie sollten vier neue Variablen (A_sigma2_inn, B_sigma2_inn, MDiff_Sigma_A, MDiff_Sigma_B) mit jeweils 100'000 Werten haben. Sie brauchen diesmal kein apply () oder Schleifen. $$\widehat\sigma_{inn}^{2}=\frac{(n_{1}−1)*s_{1}^{2}+(n_{2}−1)*s_{2}^{2}}{n_{1}+n_{2}-2}$$ $$\widehat\sigma_{\overline x_1 - \overline x_2}=\sqrt{\frac{\widehat\sigma_{inn}^{2}}{n_{1}}+\frac{\widehat\sigma_{inn}^{2}}{n_{2}}}$$ ```{r} #| eval: false #| include: false #| echo: false A_sigma2_inn <- ((n-1)*A_1_SD^2+(n-1)*A_2_SD^2)/(n+n-2) B_sigma2_inn <- ((n-1)*B_1_SD^2+(n-1)*B_2_SD^2)/(n+n-2) MDiff_Sigma_A <- sqrt(A_sigma2_inn/n+A_sigma2_inn/n) MDiff_Sigma_B <- sqrt(B_sigma2_inn/n+B_sigma2_inn/n) ``` ## 1. Konfidenzintervall Mittelwertsunterschiede - Übung 3 Bestimmen Sie den benötigten *t*-Wert für ein 95%-Konfidenzintervall für Mittelwertsunterschiede und speichern Sie ihn in der Variable t. Prüfen Sie Ihren Wert hier (2 Nachkommastellen):`r forms_string(c("1.98","1.97"), width = 7, obfuscate = FALSE)` ```{r} #| eval: false #| include: false #| echo: false t <- qt(.975,64+64-2) ``` \ Bestimmen Sie jetzt die Ober- und Untergrenzen ihrer Konfidenzintervalle für jede Stichprobe und speichern Sie die Ergebnisse in den Variablen A_ug, A_og, B_ug, B_og. ```{r} #| eval: false #| include: false #| echo: false A_ug <- (A_1_M - A_2_M) - t*MDiff_Sigma_A A_og <- (A_1_M - A_2_M) + t*MDiff_Sigma_A B_ug <- (B_1_M - B_2_M) - t*MDiff_Sigma_B B_og <- (B_1_M - B_2_M) + t*MDiff_Sigma_B ``` \ Prüfen Sie für jedes Konfidenzintervall ob die 0 ausserhalb der Ober- und Untergrenzen liegt und speichern Sie das Ergebnis in den Variablen A_Ausserhalb, B_Ausserhalb. Sie benötigen keine if-Bedingung aber Vergleiche (< >) und den oder-Operator |. ```{r} #| eval: false #| include: false #| echo: false A_Ausserhalb <- (0 < A_ug) | (0 > A_og) # ich bin gar nicht sicher, ob die Klammern notwendig sind, schaden tun sie nicht B_Ausserhalb <- (0 < B_ug) | (0 > B_og) # ich bin gar nicht sicher, ob die Klammern notwendig sind, schaden tun sie nicht table(A_Ausserhalb) table(B_Ausserhalb) ``` \ Beantworten Sie folgende Fragen: - Wie oft lag für Stichprobe A die 0 ausserhalb der Konfidenzintervalle? `r forms_string("5002", width = 7, obfuscate = FALSE)` - Wie oft lag für Stichprobe B die 0 ausserhalb der Konfidenzintervalle? `r forms_string("79923", width = 7, obfuscate = FALSE)` - Wie oft lag für Stichprobe A die 0 innerhalb der Konfidenzintervalle? `r forms_string("94998", width = 7, obfuscate = FALSE)` - Wie oft lag für Stichprobe B die 0 innerhalb der Konfidenzintervalle? `r forms_string("20077", width = 7, obfuscate = FALSE)` ## 1. Konfidenzintervall Mittelwertsunterschiede - Übung 4 \ Was bedeutet es, wenn die 0 innerhalb oder ausserhalb des Konfidenzintervalls liegt? Wir hatten eine Konfidenz von 95% gewollt. Passen unsere Ergebnisse zu diesem Wert? Begründen Sie Ihre Antwort. Population B wurde so gewählt, dass unsere Stichprobe mit einer gewissen Power/Teststärke ein signifikantes Ergebnis liefern. Wie gross sind Fehler 2. Art und die Teststärke für die Stichproben von Population B? ## 2. *t*-Test für unabhängige Stichproben - Einführung Die Statistik für den *t*-Test für unabhängige Stichproben ist definiert als: $$ t=\frac{\overline x_1 - \overline x_2}{\widehat\sigma_{\overline x_1 - \overline x_2}}$$ \ Sie gibt also an, wie gross das Verhältnis zwischen der Mittelwertsdifferenz zur geschätzten Streuung der Mittelwerte (über viele theoretische Stichproben) ist. Je grösser das Verhältnis, desto mehr können wir annehmen, dass der beobachtete Unterschied nicht nur Zufall war. \ Erkenn Sie auch, dass diese Formel der umgestellten Formel für das Konfidenzintervall entspricht, ausser das hier die empirische *t*-Wert berechnet wird, während in die Gleichung des Konfidenzintervalls der kritische *t*-Wert eingesetzt wird. ## 2. *t*-Test für unabhängige Stichproben - Übung 1 \ Sie haben eigentlich schon alles, um die empirischen *t*-Werte zu bestimmen. Bitte bestimmen Sie jetzt die empirischen *t*-Werte für beide Populationen und speichern Sie die Ergebnisse in A_t und B_t. ```{r} #| eval: false #| include: false #| echo: false A_t <- (A_1_M-A_2_M)/MDiff_Sigma_A B_t <- (B_1_M-B_2_M)/MDiff_Sigma_B ``` \ Für die Prüfung der statistischen Signifikanz könnten Sie jetzt einfach die empirischen *t*-Werte mit dem kritischen *t*-Wert von zuvor vergleichen. Wir wollen einen anderen Weg gehen. Bestimmen Sie für jeden *t*-Wert den zugehörigen *p*-Wert mit der Funktion pt(). Probieren Sie die Funktion am besten mit dem kritischen *t*-Wert von zuvor aus. Sie werden das erste, das zweite und das vierte Argument der Funktion brauchen. Denken Sie auch daran, dass der *p*-Wert für eine ungerichtete Hypothese wäre. Ein letzter Hinweis, Sie müssen (um die Richtung ignorieren zu können) die Beträge der *t*-Werte mittels abs() nehmen. Sobald Sie die Funktion verstanden haben, bestimmen sie die *p*-Werte und speichern Sie die Ergebnisse in A_p und B_p. ```{r} #| eval: false #| include: false #| echo: false 2*pt(t,n+n-2, lower.tail = FALSE) A_p <- 2*pt(abs(A_t), n+n-2, lower.tail = FALSE) B_p <- 2*pt(abs(B_t), n+n-2, lower.tail = FALSE) ``` \ Beantworten Sie folgende Fragen: - Wie oft war der *p*-Wert für Stichprobe A kleiner als .05? `r forms_string("5002", width = 7, obfuscate = FALSE)` - Wie oft war der *p*-Wert für Stichprobe B kleiner als .05? `r forms_string("79923", width = 7, obfuscate = FALSE)` - Wie oft war der *p*-Wert für Stichprobe A grösser als .05? `r forms_string("94998", width = 7, obfuscate = FALSE)` - Wie oft war der *p*-Wert für Stichprobe B grösser als .05? `r forms_string("20077", width = 7, obfuscate = FALSE)` ```{r} #| eval: false #| include: false #| echo: false #Logischerweise Sind die Häufigkeiten identisch zu den Häufigkeiten der Konfidenzinzervalle table(A_p < .05) table(B_p < .05) ``` ## 2. *t*-Test für unabhängige Stichproben - Übung 2 \ Wir wissen, dass eigentlich Subpopulation B2 einen höheren Mittelwert in der Population aufweist, als Subpopulation B1. Wie oft haben Sie ein signifikantes Ergebnis erzielt, wobei in der Stichprobe die falsche Gruppe den höheren Mittelwert hatte? Sie können diese Frage mit den Konfidenzintervallen oder den *t*-Tests lösen. `r forms_string("0", width = 7, obfuscate = FALSE)` ```{r} #| eval: false #| include: false #| echo: false #Beide Varianten table(00)) #Wir sehen, dass ist nicht einmal passiert, der Effekt ist auch recht stark in der Bevölkerung mit d = .50 und relativ grossen Stichproben. Es könnte theoretisch aber passieren, nur recht unwahrscheinlich bei diesen Werte :-D ``` \ Wir haben jetzt gesehen, dass Konfidenzintervalle und die zugehörigen parametrischen Tests exakt die gleichen Ergebnisse erzielen. Jetzt wollen wir uns noch die *p*-Werte für die nicht signifikanten Ergebnisse anschauen. Bitte erstellen sie 6 Histogramme für die *p*-Werte von - Population A - Population B - Population A, aber nur die nicht signifikanten *p*-Werte - Population B, aber nur die nicht signifikanten *p*-Werte - Population A, aber nur die signifikanten *p*-Werte - Population B, aber nur die signifikanten *p*-Werte ```{r} #| eval: false #| include: false #| echo: false hist(A_p) hist(B_p) hist(A_p[A_p>.05]) hist(B_p[B_p>.05]) hist(A_p[A_p<.05]) hist(B_p[B_p<.05]) # Was die Studis sehen sollten, p-Werte von Populationen wo kein Effekt existiert sind gleichverteilt, Wenn ein Populationseffekt existiert, wird die Verteilung immer linkssteiler/rechtsschiefer je grösser der Populationseffekt ist. Man kann sich dies auch zunutze machen um Verdachtsfälle von QRPs/p-Hacking in der Forschung zu identifizieren, Stichwort p-Curve. ``` ## 2. *t*-Test für unabhängige Stichproben - Übung 3 \ Was ist Ihnen bei den vorherigen Histogrammen aufgefallen? \ Schauen wir uns noch die Mittelwertsunterschiede in Bezug auf die zugehörigen *p*-Werte an. Erstellen Sie dazu passende Streudiagramme für - Population A - Population B - Population A, aber nur Stichproben mit nicht signifikanten *p*-Werten - Population B, aber nur nur Stichproben mit nicht signifikanten *p*-Werten - Population A, aber nur Stichproben mit signifikanten *p*-Werten - Population B, aber nur Stichproben mit signifikanten *p*-Werten ```{r} #| eval: false #| include: false #| echo: false plot(A_p,A_1_M - A_2_M) plot(B_p,B_1_M - B_2_M) plot(A_p[A_p>.05],A_1_M[A_p>.05] - A_2_M[A_p>.05]) plot(B_p[B_p>.05],B_1_M[B_p>.05] - B_2_M[B_p>.05]) plot(A_p[A_p<.05],A_1_M[A_p<.05] - A_2_M[A_p<.05]) plot(B_p[B_p<.05],B_1_M[B_p<.05] - B_2_M[B_p<.05]) # Bitte betont nochmal dass wir zwar sehen, dass p-Werte über viele Stichproben mit den Populationseffekte im Zusammenhang stehen (Muster zeigen), allerdings einzelne p-Werte nichts über den Populationseffekt aussagen (Ein p-Wert ist ist immer vom Stichprobenfehler und damit von n und der Varianz in der Population abhängig), daher brauchen wir Effektstärken und sollten nicht mehr in die Höhe von p-Werten interpretieren als, dass sie ein Signifikanz-kriterium sind. Wenn wir über Effektstärken reden wollen sind, Konfidenzintervalle und insbesondere Effektstärken definitiv die bessere Wahl. Insbesondere finden wir in beiden Populationen alle möglichen p-Werte ``` \ Was fällt Ihnen auf? ## 3. *U*-Test - Einführung \ Jetzt wollen wir noch den *U*-Test mit dem *t*-Test vergleichen. Der *U*-Test wurde wie folgt durchgeführt: 1. Ränge zuweisen über beide Gruppen hinweg 2. Rangsummen für beide Gruppen berechnen 3. Für jede Gruppe einen *U*-Wert berechnen 4. Das Minimum der beiden U-Werte ist die Teststatistik (der empirische *U*-Wert) 5. Da R die Wilcoxon-Verteilung statt der *U*-Verteilung verwendet, nehmen wir stattdessen die Rangsumme von Gruppe 1 und ziehen die minimal mögliche Rangsumme für diese Gruppe ab.Die beiden Tests sind aber äquivalent. $$U_1 = n_1*n_2+\frac{n_1*(n_1+1)}{2}-R_1$$ $$U_2 = n_1*n_2+\frac{n_2*(n_2+1)}{2}-R_2$$ $$W = R_1-\frac{n_1*(n_1+1)}{2}$$ ## 3. *U*-Test - Übung 1 \ Da die Rangbestimmung etwas kompliziert ist nutzen Sie bitte zunächst folgenden Code, um die Ränge für die jeweilige Subpopulation 1 zu ermitteln: ```{r} #| include: true #| eval: false #| echo: true n <- 64 A_1_Rang <- apply(rbind(A_1, A_2),2,rank)[1:n,] B_1_Rang <- apply(rbind(B_1, B_2),2,rank)[1:n,] ``` \ Nachdem Sie die Ränge wie oben ermittelt haben, bestimmen Sie die Rangsumme je Stichprobe für Gruppe 1 und ziehen Sie davon die theoretisch niedrigste mögliche Rangsumme ab. Speichern Sie Ihre Ergebnisse in A_W und B_W. ```{r} #| eval: false #| include: false #| echo: false A_W <- apply(A_1_Rang,2,sum)-(n*(n+1))/2 B_W <- apply(B_1_Rang,2,sum)-(n*(n+1))/2 ``` ## 3. *U*-Test - Übung 2 \ Der Umgang mit pwilcox ist etwas kompliziert, da die Verteilung nicht um 0 symmetrisch ist sondern der Erwartungswert unter der Nullhypothese bei $(n_1*n_2)/2$ liegt. Daher hier der Code um die *p*-Werte zu bestimmen ```{r} #| include: true #| eval: false #| echo: true # Da n1 = n2 ist gilt: Erwartungswert <- 64^2/2 # normalerweise bekommt pwilcox, als Argumente, W, n_1, n_2, bei uns sind die letzten beiden gleich A_W_p <- 2*(pwilcox(abs(A_W-Erwartungswert)+Erwartungswert-1, n, n, lower.tail = FALSE) ) B_W_p <- 2*(pwilcox(abs(B_W-Erwartungswert)+Erwartungswert-1, n, n, lower.tail = FALSE) ) ``` Beantworten Sie nun folgende Fragen: - Wie oft war der *p*-Wert für Stichprobe A kleiner als .05? `r forms_string("4976", width = 7, obfuscate = FALSE)` - Wie oft war der *p*-Wert für Stichprobe B kleiner als .05? `r forms_string("77962", width = 7, obfuscate = FALSE)` - Wie oft war der *p*-Wert für Stichprobe A grösser als .05? `r forms_string("95024", width = 7, obfuscate = FALSE)` - Wie oft war der *p*-Wert für Stichprobe B grösser als .05? `r forms_string("22038", width = 7, obfuscate = FALSE)` - Im Vergleich zu unseren Ergebnissen für Konfidenzintervalle und parametrische Signifikanztests, was können Sie schlussfolgern? ```{r} #| eval: false #| include: false #| echo: false # Wir sehen, dass der Alphafehler eigentlich verschwindend gefallen ist (reiner Zufall), der Betafehler deutlich mehr (über 2-Prozentpunkt). Damit ist auch die Teststärke gesunken (da 1-Beta = Teststärke). Wir sehen nonparametrische Tests haben tendenziell eine geringere Power/Teststärke als parametrische Tests, wenn die Voraussetzungen der parametrischen Tests erfüllt sind. Die Alphafehler sind gleich. table(A_W_p < .05) table(B_W_p < .05) ```