Daten & Packages laden

Laden Sie die folgenden Packages und Data-Frames:

library(tidyverse)
library(afex)
library(emmeans)
library(broom)
urla = "https://www.phonetik.uni-muenchen.de/studium_lehre/"
urlb = "lehrmaterialien/R_speech_processing/Rdf"
url = paste0(urla, urlb)
preasp = read.table(file.path(url, "preasp.txt"), 
                    stringsAsFactors = T)
ital.df = preasp %>%
  filter(vtype == "e" & cplace == "kk") %>% 
  mutate(spk = factor(spk))
rownames(ital.df) = 1:nrow(ital.df)
nv.df = read.table(file.path(url, "nv.df.txt"), 
                    stringsAsFactors = T)
nv.df = nv.df %>% 
  mutate(Context = 
           factor(Context, levels = c("a", "na", "an")))

Varianzanalyse

Q1: Der Dataframe nv.df besteht aus folgenden Spalten:
Speaker: Sprecher
Sex: M oder F
Context: a, na, an. Die Bedeutung davon: /a/ ist ein oraler Vokal (in z.B. ‘pad’), na ein Vokal nach einem nasalen Konsonant (in z.B. ‘mad’), und an ein Vokal vor einem nasalen Konsonanten (in z.B. ‘pan’).
F1: Frequenz von F1

Prüfen Sie durch eine Abbildung und statistischen Test inwiefern F1 von Sex und oder Context beeinflusst wird.

Test für within und between: Sex ist between, Context ist within.

nv.df %>%
  mutate(Fac = interaction(Sex, Context)) %>%
  select(Speaker, Fac) %>%
  table()

##        Fac
## Speaker F.a M.a F.na M.na F.an M.an
##     S1    0   1    0    1    0    1
##     S10   1   0    1    0    1    0
##     S2    0   1    0    1    0    1
##     S3    0   1    0    1    0    1
##     S4    0   1    0    1    0    1
##     S5    0   1    0    1    0    1
##     S6    1   0    1    0    1    0
##     S7    1   0    1    0    1    0
##     S8    1   0    1    0    1    0
##     S9    1   0    1    0    1    0

Die Abbildung zeigt, dass Kontext einen erheblichen Einfluss auf F1 hat, sowohl in Männern als auch in Frauen. (Und wie zu erwarten: F1 ist höher in Frauen als in Männern). Wir erwarten eine Interaktion zwischen diesen Faktoren, da der F-M Unterschied in na viel größer ist als in den anderen Kontexten. Darüber hinaus erwarten wir eine Interaktion, weil der Unterschied zwischen an und den anderen beiden Kontexten viel größer ist als zwischen a und na.

nv.df %>%
  ggplot +
  aes(y = F1, x = Context, col = Sex) +
  geom_boxplot()

Statistischer Test: Es gab einen signifikanten Einfluss auf F1 von Sex (F[1,8] = 11.08, p < 0.05) und von Context (F[1.62, 12.94] = 56.38, p < 0.001) und eine signifikante (F[1.62, 12.94] = 5.51, p < 0.05) Interaktion zwischen diesen beiden Faktoren.

nv.aov = nv.df %>%
  aov_4(F1 ~ Sex * Context + (Context|Speaker), .)

## Contrasts set to contr.sum for the following variables: Sex

nv.aov

## Anova Table (Type 3 tests)
## 
## Response: F1
##        Effect          df     MSE         F  ges p.value
## 1         Sex        1, 8 3778.87   11.08 * .376    .010
## 2     Context 1.62, 12.94 3046.82 56.38 *** .800   <.001
## 3 Sex:Context 1.62, 12.94 3046.82    5.51 * .281    .023
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG

Post-hoc Tests zeigten signifikante Unterschiede zwischen Männern und Frauen nur in an (p < 0.01) aber nicht in den anderen beiden Kontexten. Es gab weiterhin signifikante Unterschiede zwischen a und an (Frauen: p < 0.001; Männer: p < 0.05), und zwischen na und an (Frauen: p < 0.001; Männer: p < 0.05) jedoch nicht zwischen a und na.

nv.aov %>%
  emmeans(., ~ Sex * Context) %>%
  pairs(., simple = "each", combine=T)

##  Context Sex contrast estimate   SE df t.ratio p.value
##  a       .   F - M        29.8 38.2  8   0.780  1.0000
##  na      .   F - M        34.6 38.2  8   0.907  1.0000
##  an      .   F - M       159.8 23.8  8   6.723  0.0013
##  .       F   a - na      -46.0 28.4  8  -1.622  1.0000
##  .       F   a - an     -287.4 38.2  8  -7.521  0.0006
##  .       F   na - an    -241.4 26.3  8  -9.182  0.0001
##  .       M   a - na      -41.2 28.4  8  -1.453  1.0000
##  .       M   a - an     -157.4 38.2  8  -4.119  0.0301
##  .       M   na - an    -116.2 26.3  8  -4.420  0.0200
## 
## P value adjustment: bonferroni method for 9 tests

Regression

Q2: Prüfen Sie durch eine Abbildung und statistischen Test für den Dataframe trees, inwiefern Height aus Volume vorhersagbar ist. Schätzen Sie Height ein für ein Volumen von 110. Führen Sie zwei diagnostische Tests durch, um die Gültigkeit des statistischen Tests zu prüfen.

trees %>%
  ggplot() +
  aes(y = Height, x = Volume) +
  geom_point() +
  geom_smooth(method = "lm", color = "blue")

## `geom_smooth()` using formula = 'y ~ x'

Es gibt zumindest eine Tendenz, dass sie miteinander linear assoziiert sind.

trees %>%
  lm(Height ~ Volume, .) %>%
  summary()

## 
## Call:
## lm(formula = Height ~ Volume, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.7777  -2.9722  -0.1515   2.0804  10.6426 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 69.00336    1.97443  34.949  < 2e-16 ***
## Volume       0.23190    0.05768   4.021 0.000378 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.193 on 29 degrees of freedom
## Multiple R-squared:  0.3579, Adjusted R-squared:  0.3358 
## F-statistic: 16.16 on 1 and 29 DF,  p-value: 0.0003784

Es gibt eine signifikante (\(R^2\) = 0.3579, F[1,29] = 16.16, p < 0.001) lineare Beziehung zwischen Height und Volume.

trees %>%
  lm(Height ~ Volume, .) %>%
  predict(., data.frame(Volume = 110))

##        1 
## 94.51234

Ist der Mittelwert der Residuals nah an 0 (Null)? Ja.

trees %>%
  lm(Height ~ Volume, .) %>%
  augment() %>%
  summarise(mean(.resid))

## # A tibble: 1 × 1
##   `mean(.resid)`
##            <dbl>
## 1      -2.43e-14

Ist es wahrscheinlich, dass die Residuals normalverteilt sind? Ja.

# Ja
trees %>%
  lm(Height ~ Volume, .) %>%
  augment() %>%
  pull(.resid) %>%
  shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.9822, p-value = 0.8707

Q3: Für diese Daten wurde F2-F1 (Hz) in einem Vokal gemittelt über mehrere Sprecher zwischen 1910 und 1997 gemessen. Prüfen Sie durch eine Abbildung, ob sich F2-F1 mit der Zeit ändert. Schätzen Sie den Wert von F2-F1 ein im Jahr 2012. Führen Sie vier Tests durch, um die Gültigkeit des statistischen Tests zu prüfen.

jahr = c(1910, 1920, 1930, 1940, 1950, 1959, 1969, 1978, 1987, 1997)
fdat = c(139, 149, 157, 175, 216, 303, 390, 449, 462, 487)

# Dataframe bauen
f.df = data.frame(Jahr = jahr, Fdiff = fdat)

# Bild
f.df %>%
  ggplot() +
  aes(y = Fdiff, x = Jahr) +
  geom_point() +
  geom_smooth(method = "lm", se = F, colour = "blue")

## `geom_smooth()` using formula = 'y ~ x'

Es gibt eine Beziehung, aber eventuell nicht linear (sondern S-formig)

f.df %>%
  lm(Fdiff ~ Jahr, .) %>%
  predict(., data.frame(Jahr = 2012))

##       1 
## 566.341

Test-Gültigkeit:

# Das ist OK
f.df %>%
  lm(Fdiff ~ Jahr, .) %>%
  augment %>%
  summarise(mean(.resid))

## # A tibble: 1 × 1
##   `mean(.resid)`
##            <dbl>
## 1       9.99e-13

# Das ist OK
f.df %>%
  lm(Fdiff ~ Jahr, .) %>%
  augment %>%
  pull(.resid) %>%
  shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.95893, p-value = 0.7736

# Das ist eventuell problematisch: Die Residuals folgen einer S-Kurve
f.df %>%
  lm(Fdiff ~ Jahr, .) %>%
  plot(which = 1)

# Aber das ist OK.
f.df %>%
  lm(Fdiff ~ Jahr, .) %>%
  augment() %>%
  pull(.resid) %>%
  acf()

Q4: Für den Dataframe ital.df prüfen Sie durch eine Abbildung und statistischen Test, ob eine lineare Beziehung zwischen Vokaldauer (vdur) und Dauer der Verschlusslösung (clorel) vorliegt. Prüfen Sie die Gültigkeit des statistsichen Tests.

Es gibt eine leicht positive Beziehung zwischen den Variablen.

ital.df %>%
  ggplot +
  aes(y = vdur, x = clorel) +
  geom_point() +
  geom_smooth(method = "lm", col="blue")

## `geom_smooth()` using formula = 'y ~ x'

Test ob die Steigung signifikant von 0 (Null) abweicht. Ja das ist der Fall: Es gibt eine signifikante, lineare Beziehung zwischen Vokaldauer und Verschlusslösung (\(R^2\) = 0.3579, F[1,117] = 7.49, p < 0.01)

ital.df %>%
  lm(vdur ~ clorel, .) %>%
  summary()

## 
## Call:
## lm(formula = vdur ~ clorel, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.034642 -0.008564 -0.000372  0.007199  0.054463 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.039848   0.007764   5.132 1.15e-06 ***
## clorel      0.111661   0.040798   2.737  0.00717 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01209 on 117 degrees of freedom
## Multiple R-squared:  0.06017,    Adjusted R-squared:  0.05214 
## F-statistic: 7.491 on 1 and 117 DF,  p-value: 0.007169

Gültigkeit von dem Test.

# Gültigkeittest 1: Mittelwert. Alles OK
ital.df %>%
  lm(vdur ~ clorel, .) %>%
  augment() %>%
  summarise(mean(.resid))

## # A tibble: 1 × 1
##   `mean(.resid)`
##            <dbl>
## 1      -1.47e-17

# Gültigkeittest 2: konstante Varianz in den Residuals? 
# Ja - die Linie ist mehr oder weniger flach
ital.df %>%
  lm(vdur ~ clorel, .) %>%
  plot(which = 1)

# Gültigkeittest 3: Autokorrelation in den Residuals?
# Alles OK - fast alle Lagwerte sind zwischen ±0.2
ital.df %>%
  lm(vdur ~ clorel, .) %>%
  augment() %>%
  pull(.resid) %>%
  acf()

# Gültigkeittest 4: normalverteilt?
# Nein: p < 0.05
ital.df %>%
  lm(clorel ~ vdur, .) %>%
  augment() %>%
  pull(.resid) %>%
  shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.97585, p-value = 0.03058

Den QQplot anschauen. Es gibt Ausreißer: Beobachtungen 32, 55, 88

ital.df %>%
  lm(clorel ~ vdur, .) %>%
  plot(which = 2)

# Diese entfernen, und noch einmal prüfen, ob
# die Residuals dann normalverteilt sind. Ja
# jetzt ist es OK
ital.df %>%
  slice(-c(32, 55, 88)) %>%
  lm(clorel ~ vdur, .) %>%
  augment() %>%
  pull(.resid) %>%
  shapiro.test()

## 
##  Shapiro-Wilk normality test
## 
## data:  .
## W = 0.98788, p-value = 0.3899

Daher den Test noch einmal durchführen, ohne Ausreißer.

# Bild
ital.df %>%
  slice(-c(32, 55, 88)) %>%
  ggplot +
  aes(y = vdur, x = clorel) +
  geom_point() +
  geom_smooth(method = "lm", col="blue")

## `geom_smooth()` using formula = 'y ~ x'

# Test
ital.df %>%
  slice(-c(32, 55, 88)) %>%
  lm(clorel ~ vdur, .) %>%
  summary()

## 
## Call:
## lm(formula = clorel ~ vdur, data = .)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.052764 -0.014475  0.001263  0.016269  0.047366 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.15669    0.01125   13.93  < 2e-16 ***
## vdur         0.48618    0.18144    2.68  0.00846 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02329 on 114 degrees of freedom
## Multiple R-squared:  0.05925,    Adjusted R-squared:  0.051 
## F-statistic:  7.18 on 1 and 114 DF,  p-value: 0.008462

Es gab eine signifikante (\(R^2\) = 0.05925, F[1,114] = 7.18, p < 0.01) lineare Beziehung zwischen der Vokal- und Lösungsdauer.

Übung 8

Jonathan Harrington / Johanna Cronenberg

Daten & Packages laden

Varianzanalyse

Regression