Regression (zweiter Teil)

Regression (zweiter Teil)

Jonathan Harrington

1. Regression und die Locus-Theorie

2. Voraussetzungen für die Durchführung einer Regression

3. Mehrfache (multiple) Regression: zwei oder mehrere Regressoren

Neigungen von Regressionslinien miteinander vergleichen

4. Polynomiale Regression

// vor hinteren Vokalen: keine einheitliche Locus-Frequenz

F2-Transitionenen und die Locus-Theorie

Die F2-Transitionen, die die Artikulationsstelle auditiv

vermitteln, sind durch eine Kombination vom Locus und der Vokalzielposition voraussagbar

Frequenz

Dauer

i

b F2-Lokus

/b/ ca. 700 Hz d

/d/ ca. 1800 Hz



// vor vorderen Vokalen: 3000 Hz

In der gesprochenen Sprache weichen F2-Onsets (F2-

Offsets) von der theoretischen Locusfrequenz ab, teilweise wegen antizipatorischer oder perzeveratorischer V-auf-K Koartikulation.

F2 (Hz)

d i:

o:

d d

d

0 50 150 250

50010002000

Time (ms)

i:

o:

Max. V auf K Koartikulation

Dauer 500

2000

Dauer 500

2000

bb 800

b

Keine V-auf-K Koartikulation

F2-Locus und V-auf-K Koartikulation

• daher ist der F2 Target vom F2 Onset nicht vorhersagbar

• F2 Target ist vom F2 Onset vorhersagbar

• Kein Locus • Locus ist vom Vokal unabhängig

(Regressionlinie im Raum von F2 Target x F2 Onset)

• Die Neigung liegt zwischen 0 und 1

• Je steiler (näher an 1) die Neigung, umso bedeutender die V auf K Koartikulation

Locusgleichung

0 500 1500 2500

050015002500

bb b

F2(Onset)=Locus,

Regressionsneigung = 0

Dauer 500

2000

F2 Onset 800

bb b

F2 Target

0 500 1500 2500

050015002500

bb b

F2(Target) = F2(Onset) Regressionsneigung =1

Dauer 500

2000

F2 Target

F2 Onset

V-auf-K Koartikulation in /dV/ im Vgl. zu /Vd/

Ein Sprecher (australisch Englisch) erzeugte /dVd/ Silben, fuer verschiedene Vokale.

In welchen Reihenfolgen, initiale /dV/, oder finale /Vd/, ist die V-auf-K Koartikulation wahrscheinlich am

bedeutendsten: initial oder final?

F2 (Hz)

o:

d

d

0 50 150 250

50010002000

Dauer (ms)

V-target = 600 Hz, F2in=1800 Hz, F2off=1300 Hz

600 1000 1400

10001600

F2-Target

F2-initial do:

600 1000 1400

10001600

F2-Target

F2 offset

o:d

/dV/ /Vd/

F2-Target F2-Target

F2-Onset F2-Offset

Welche Regressionslinie müsste steiler sein?

dV Vd

0 100 200 300 400

5001500

time (ms)

1

-400 -300 -200 -100 0

5001500

time (ms)

1

F2-verschiedene Vokale

pfad = "Das Verzeichnis, wo die Daten gespeichert ist"

fdat = read.table(paste(pfad, "dform.txt", sep="/"), header=T) attach(fdat) names(fdat)

"lab" "F2targ" "F2in" "F2fin"

F2 /dVd/-Daten einlesen

1000 1500 2000

140016001800

F2-targ

F2-initial

A i:

I

a:V

u:

Ee:

o:

U @:

O

I

plot(F2targ, F2in, type="n", xlab="F2-targ", ylab="F2-initial") text(F2targ, F2in, lab)

Regressionslinie berechnen

k mF

F ^{^} 2

 2



regin = lm(F2in ~ F2targ)

überlagerte Regressionslinie abline(regin)

Dasselbe für F2-Offset

k mF

F ^{^} 2

 2



regfin = lm(F2fin ~ F2targ) Koeffiziente F2-Initial:

coef(regin)

(Intercept) F2targ

1220.2502132 0.2710194

coef(regfin)

Intercept) F2targ

829.4801710 0.4641336 Koeffiziente F2-Final:

Nebenbei: Vorhersage-Intervall

lmint(F2targ, F2in, level=.95)

1000 1500 2000

140016001800

x

y

Die Wahrscheinlichkeit, dass Werte innerhalb der blauen Linie fallen = 95%

pfad = "Das Verzeichnis, wo die Daten gespeichert ist"

source(paste(pfad, "lmint.S", sep="/"))

möglicher Ausreißer

Die Neigungen vergleichen

Hilfeseite: siehe Slope.test.pdf in unserer Webseite.

pfad = "Das Verzeichnis, wo die Daten gespeichert ist"

source(paste(pfad, "Slope.test.S", sep="/")) args(Slope.test) # siehe unten, Hilfseite

function (...) … bedeutet: beliebig viele Argumente.

jedes Argument soll eine zweispaltige Matrix sein. Die y-Werte (F2in oder F2fin) in Spalte 1, die x-Werte (F2-targ) in Spalte 2.

Slope.test(mat1, mat2)

die Regressionslinie für die Daten in mat1 wird berechnet;

die Regressionslinie für die Daten in mat2 wird berechnet;

die Neigungen dieser beiden R-Linien werden statistisch miteinander verglichen.

$separate

r-sq F ratio df df prob. line fits data intercept slope first.out 0.8632720 69.45167 1 11 0.9999956 1220.2502 0.2710194 first.out 0.7437726 31.93061 1 11 0.9998513 829.4802 0.4641336

$combined

F ratio Probability of them being DIFFERENT df df intercept 5.971148 0.9773694 1 23 slope 4.778655 0.9602575 1 22

Die Neigungen der Regressionslinien unterscheiden sich F=4.78, df = 1, 22, p < 0.05 oder F(1,22) = 4.78, p < 0.05 onset = cbind(F2in, F2targ)

offset = cbind(F2fin, F2targ) Slope.test(onset, offset)

# 2-spaltige Matrix der F2in-Werte

# 2-spaltige Matrix der F2fin-Werte

Da sich die Regressionslinien signifikant

unterscheiden, ist die V-auf-K Koartikulation unterschiedlich in /Vd/ im Vergleich zu /dV/.

2. Voraussetzung für die

Durchführung der Regression

siehe vor allem http://www.duke.edu/~rnau/testing.htm und E-Buch in der LMU Autor = Verzani Kap. 10)

qqnorm(resid(regin))

qqline(resid(regin)) plot(regin, 2)

-1.5 -0.5 0.5 1.5

-1.5-0.50.51.5

Theoretical Quantiles

Standardized residuals

lm(F2in ~ F2targ) Normal Q-Q

10 1 12

2.1. Folgen die 'residuals' einer Normalverteilung?

oder

(siehe regression.ppt, vorige Woche)

shapiro.test(resid(regin))

Shapiro-Wilk normality test data: resid(regin)

W = 0.9274, p-value = 0.3148

Zusätzlicher Test

Die Verteilung der

residuals weicht nicht signifikant von einer Normalverteilung ab Die Werte sollen nicht

allzusehr von der

geraden Linie abweichen

2.2. Ist Linearität erkennbar?

Sind wir sicher, dass die y/x Werte wirklich mit einer Linie modelliert werden können?

0.5 1.0 1.5 2.0 2.5 3.0

10002000

x

y

0.5 1.0 1.5 2.0 2.5 3.0

10002000

x

y

oder vielleicht nicht besser mit einer Parabel?

oder

plot(predict(regin), resid(regin)) plot(regin, 1) Die Werte in einer Abbildung der residuals als

Funktion der eingeschätzen Werte, y, sollten mehr oder weniger auf einer randomisierten Weise um die Linie residuals = 0 verteilt sein.

2.2. Test für Linearität

1400 1600 1800

-100-50050

Fitted values

Residuals

lm(F2in ~ F2targ) Residuals vs Fitted

10

12 1

^

2. 3. Keine Autokorrelation

(die aufeinanderfolgenden Werte in der x oder y- Achse müssen voneinander unabhängig sein).

16460 16470 16480 16490 16500

-600004000

times

data[, k]

Ein gutes Beispiel von autokorrelierten Daten: ein stimmfahfes Sprachsignal

Man soll überhaupt sehr vorsichtig sein, eine Regression auf Daten mit einer Zeit-Achse anzuwenden (da Werte als Funktion der Zeit sich oft wiederholen, also sie sind oft

miteinander korreliert).

2. 3. Test für Autokorrelation

Die Autokorrelation der residuals berechnen

acf(resid(regin))

Autokorrelation lag k

inwiefern sind die ersten n-k Werte des Signals mit den letzten n-k Werte korreliert?

zB Autokorrelation lag 1 bei einem 4-Punkt Signal:

die ersten 3 Werte werden mit den letzten 3 Werte verglichen Autokorrelation lag 0: kann ignoriert werden: sie

hat immer den Wert 1.

Alle anderen lag-Werte variieren zwischen -1 und 1

95% Vertrauensintervall um 0 n = length(F2targ)

2/sqrt(n)

Insbesondere die Werte bei lag 1 und 2 beobachten: diese

sollten in jedem Fall innerhalb des Vertauensintervalls liegen.

Wenn die meisten ACF-Werte innerhalb der blauen Linien liegen, gibt es keine Autokorrelation.

2. 3. Test für Autokorrelation

0 2 4 6 8 10

-0.50.00.51.0

Lag

ACF

2. 4. Konstante Varianz oder

'homoscedasticity' der Residuals?

Insbesondere sollten die residuals nicht wesentlich größer am Anfang/Ende sein, sondern auf eine

randomisierte Weise um die 0 Linie verteilt sein.

2 4 6 8 10

-50050

Index

resid(regin)

plot(resid(regin))

2.5. Ausreißer

Ein einziger Ausreißer vor allem vom Mittelpunkt weit entfernt kann die Regressionslinie deutlich beeinflussen.

Ausreißer die eventuell (aber nicht unbedingt) die Regressionslinie stark beeinflussen, können mit dem sog. Cookes Distance aufgedeckt werden

2 4 6 8 10 12

0.00.10.20.30.4

Obs. number

Cook's distance

lm(F2in ~ F2targ) Cook's distance

1

12

3

plot(regin, 4)

Die Werte 1 und 12 sind vielleicht solche Ausreißer, also

onset[c(1,12),]

plot(F2targ, F2in, cex = 10*sqrt(cooks.distance(regin))) text(F2targ, F2in, as.character(1:length(F2in)))

Die Cookes-Entfernungen und daher Ausreißer könnten auch mit einem sogenannten 'bubble plot' (siehe Verzani) abgebildet werden

2.5. Ausreißer

1000 1500 2000

14001500160017001800

F2targ

F2in

2 1 3

5 4

6

78

9

11 10

12

13

Ausreißer

Man könnte dann prüfen, ob sich die Regressionlinien signifikant mit/ohne Ausreißer unterscheiden:

onsetohne = onset[-c(1,12),]

Slope.test(onset, onsetohne)

$combined

F ratio Probability of them being DIFFERENT df df intercept 0.38475966 0.45825958 1 21 slope 0.01596279 0.09927878 1 20

detach(fdat)

2.5. Ausreißer

3. Mehrfache Regression

Mehrfache Regression

k bx

y  

Einfache Regression ^

k x

b x

b

y  _{1 1}  _{2 2} 

^

x₁

x₂ y

Und eine gerade Linie in einem 3D-Raum

In diesem Fall: 2 Regressors (x₁, x₂) , 2 Neigungen (b₁, b₂), ein Intercept, k

Es können auch mehrere Regressoren sein…

k x

b x

b x

b x

b

y  _{1 1}  _{2 2}  _{3 3}  ... _{n n} 

^

Eine gerade Linie in einem n-dimensionalen Raum n verschiedene Neigungen, ein Intercept

Mehrfache Regression

Einige Daten

pfad = "Das Verzeichnis, wo die Daten gespeichert ist"

ydata = read.table(paste(pfad, "ydata.txt", sep="/"), header=T) attach(ydata)

names(ydata)

[1] "F2" "DORSY" "DORSX" "LIPY" "LIPX"

[y] Vokale, alle Werte zum zeitlichen Mittelpunkt

DORSX, DORSY (horizontale und vertikale Position des Zungendorsums)

LIPX, LIPY (horizontale Verlagerung und vertikale Position der Unterlippe)

Wir wollen versuchen, den F2-Wert aus diesen artikulatorischen Parametern vorherzusagen.

Einige informellen Vorhersagen F2 steigt wenn:

DORSY steigt

DORSX nach vorne

LIPY

LIPX

steigt/fällt

nach vorne/hinten fällt

(die Lippen werden offener) nach hinten

(Lippen sind nicht so gerundet)

steigt/fällt

nach vorne/hinten

F2 = b^ ₁DORSX +b₂DORSY + b₃LIPX + b₄ LIPY + k

Ein mehrfaches Regressionsmodell

Mehrfache Regression in R

Festlegung von b₁, b₂, b₃, b₄, k

Sind alle diese Parameter notwendig? Wenn nicht, welches Parameter oder Parameterkombination hat die deutlichste lineare Beziehung zu F2?

hat die Bedeutung

F2 = ein Gewicht mal die horizontale Position der Zunge + ein anderes Gewicht mal die vertikale Position der Zunge +

…. + ein Intercept

^

pairs(ydata)

F2

0.8 1.0 1.2 1.4 -1.8 -1.6 -1.4

145016001750

0.81.11.4

DORSY

DORSX

3.43.84.2

-1.8-1.5

LIPY

1450 1600 1750 3.4 3.8 4.2 -1.4 -1.0 -0.6

-1.4-1.0-0.6

LIPX

hoch tief

vorne hinten*

oben unten

hinten vorne

* Richtung Glottis

F2 = b^ ₁DORSX +b₂DORSY + b₃LIPX + b₄ LIPY + k

regm = lm(F2 ~ DORSX+DORSY+LIPX+LIPY)

coef(regm)

Normalverteilung der Residuals?

plot(regm, 2) shapiro.test(resid(regm))

usw. wie für einfache Regression Koeffiziente

summary(regm)

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) 1355.05 822.63 1.647 0.113 DORSX -38.33 157.21 -0.244 0.810 DORSY 63.08 254.27 0.248 0.806 LIPX -67.36 110.90 -0.607 0.550 LIPY -164.08 205.04 -0.800 0.432

Residual standard error: 83 on 23 degrees of freedom

Multiple R-Squared: 0.3939, Adjusted R-squared: 0.2884 F-statistic: 3.736 on 4 and 23 DF, p-value: 0.01747

F2 kann mit einer multidimensionalen Regressionslinie aus diesen artikulatorischen Parametern modelliert werden:

Adjusted R² = 0.29, F (4, 23) = 3.74, p < 0.05.

Adjusted R

R²: (siehe vorige Vorlesung) ist die Proportion der Varianz, die durch die Regressionslinie erklärt werden kann (variiert

zwischen 0 und 1; 1 bedeutet alle Werte liegen auf der Linie)

Daher muss in der Berechnung von R² für die Anzahl der Regressoren kompensiert werden, wenn wir - wie in diesem Fall –Regressionslinien mit unterschiedlichen Anzahlen von Regressoren miteinander vergleichen wollen.

R² wird mit einer zunehmenden Anzahl von Regressoren größer.

Adjusted R

1 ) 1

1 (

1 ²





 

 n k

R n

Adjusted R² =

n ist die Anzahl der Stichproben, k ist die Anzahl der Regressoren.

1-(1-0.3939) * ( (n-1)/(n-4-1) ) n = length(F2)

[1] 0.2884913 Für diesen Fall

kann auch negativ sein ist weniger als R²

Modell-Prüfung durch AIC (Akaike's Information Criterion)

Mit der stepAIC() Funktion in library(MASS) wird

geprüft, ob für die Regression wirklich alle (in diesem Fall 4) Regressoren benötigt werden.

Je kleiner AIC, umso nützlicher die Kombination für die Regression (umso höher adjusted R²)

library(MASS) stepAIC(regm)

Start: AIC= 251.95

F2 ~ DORSX + DORSY + LIPX + LIPY Df Sum of Sq RSS AIC - DORSX 1 410 158861 250 - DORSY 1 424 158875 250 - LIPX 1 2541 160993 250 - LIPY 1 4412 162863 251

<none> 158451 252

sortiert nach AIC. Dies ist der AIC-Wert, wenn aus diesem Modell DORSX weggelassen wäre.

Daher wird im nächsten Modell DORSX weggelassen.

Vor allem ist dieser AIC Wert weniger als AIC mit allen Parametern (= 251.95).

Step: AIC= 250.02

F2 ~ LIPX + LIPY + DORSY

Df Sum of Sq RSS AIC - DORSY 1 1311 160172 248 - LIPY 1 4241 163102 249

<none> 158861 250 - LIPX 1 16377 175238 251

AIC ist am kleinsten, wenn aus diesem Modell DORSY weggelassen wird. Und dieser Wert ohne DORSY ist kleiner als derjenige mit LIPX+LIPY+DORSY

zusammen.

Daher wird DORSY weggelassen…

Step: AIC= 248.25 F2 ~ LIPX + LIPY

Df Sum of Sq RSS AIC

<none> 160172 248 - LIPX 1 25225 185397 250 - LIPY 1 50955 211127 254

Wenn wir entweder LIPX oder LIPY weggelassen, dann wird AIC höher im Vergleich zu AIC mit beiden Parametern

zusammen.

Daher bleiben wir bei F2 ~ LIPX + LIPY

Dieses Modell F2 ~ LIPX + LIPY müsste auch den höchsten adjusted R² haben. Prüfen, zB:

summary(regm) summary(lm(F2~LIPX+LIPY)) Adjusted R-squared: 0.3383 Adjusted R-squared: 0.2884

Also wird die Variation in F2 in [y] am meisten durch die horizontale und vertikale Position der Unterlippe erklärt.

4. Polynomiale Regression

k bx

y  

^

ein Regressor, 2 Koeffiziente

k x

b x

b

y^{^}  ₁  ₂ ² 

ein Regressor ein Koeffizient

bestimmt die Krümmung;

b₂ ist negativ b₂ ist positiv b₂ ist näher an 0

k x

b x

b x

b x

b x

b

y^{^}  ₁  ₂ ²  ₃ ³  ₄ ⁴  ... _n1 ⁿ 

ein Regressor, n Koeffiziente

In allen Fällen handelt es sich um

Abbildung/Beziehungen im 2D-Raum (da wir mit einem Regressor zu tun haben).

detach(ydata)

attach(edat) plot(COG, F2)

0.5 1.0 1.5 2.0 2.5 3.0

1000150020002500

COG

F2

pfad = "Das Verzeichnis, wo die Daten gespeichert ist"

edat = read.table(paste(pfad, "epg.txt", sep="/"))

k COG

b COG

b

F^{^}2  ₁  ₂ ² 

regp = lm(F2 ~ COG + I(COG^2)) coef(regp)

(Intercept) COG I(COG^2) -294.3732 2047.8403 -393.5154

294.3732 393.5

2047.8

2  COG  COG² 

F^{^}

plot(COG, F2)

Die Parabel überlagern

curve(-294.3732+2047.8403*x -393.5154*x^2, add=T)

Coefficients:

Estimate Std. Error t value Pr(>|t|) (Intercept) -294.37 139.95 -2.103 0.0415 * COG 2047.84 192.83 10.620 1.81e-13 ***

I(COG^2) -393.52 54.17 -7.264 6.10e-09 ***

---

Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 202.1 on 42 degrees of freedom

Multiple R-Squared: 0.9092, Adjusted R-squared: 0.9049 F-statistic: 210.4 on 2 and 42 DF, p-value: < 2.2e-16

summary(regp)

Beide Komponente, COG und COG² der Parabel scheinen notwendig zu sein, um die F2-COG

Beziehungen zu modellieren.

mit stepAIC() kann wieder festgestellt werden, ob wir den COG² Parameter wirklich benötigen:

stepAIC(regp)

Start: AIC= 480.69 F2 ~ COG + I(COG^2)

Df Sum of Sq RSS AIC

<none> 1715550 481 - I(COG^2) 1 2155469 3871019 515 - COG 1 4606873 6322423 537 Call:

lm(formula = F2 ~ COG + I(COG^2)) Coefficients:

(Intercept) COG I(COG^2) -294.4 2047.8 -393.5

Scheinbar ja (da AIC höher wird, wenn entweder COG oder COG² aus dem Modell weggelassen werden).