Statistik II

Lineare Algebra, Multivariate Statistik, ALM

Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.

Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Psychologische Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.

Versionshistory:

• v1.11: minimale Korrekturen (3.6.2024)
• v1.1: Ergänzung Variablen- vs. Personenraum (Abschnitt 10.1.4) (27.5.2024)
• v1.0: erste online-gestellte Version (3.3.2024)

10 Einführung in die Lineare Algebra

Dieser Teil befasst sich – auf den ersten Blick – mit einem gar nicht so statistisch anmutendem Thema, nämlich der Linearen Algebra. Allerdings werden in den folgenden Teilen Vektoren und Matrizen sich als äußerst nützliche und kompakte Schreibweisen erweisen, um bestimmte Themen darstellen zu können. Daher werden wir hier Grundlagen von Vektoren und Matrizen und deren Rechenregeln behandeln, einige Eigenschaften diskutieren und am Ende etwas ausführlicher die Determinante einer Matrix sowie Eigenwerte und Eigenvektoren einer Matrix behandeln.

10.1 Vektoren

Ein Vektor ist eine Zusammenfassung von \(n\) Zahlen \(x_i\) (auch Elemente genannt). Wir beschränken uns dabei auf den Raum der reellen Zahlen, d.h. jedes Element des Vektors ist eine reelle Zahl, also \(x_i \in \mathbb{R}\). Im Folgenden sind Vektoren immer klein und fett gedruckt und die Notation lautet dann:

\[\begin{align*} \bf x &= \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} \end{align*}\]

In Kurzform schreiben wir \(\mathbf{x}\in \mathbb{R}^n\), wobei \(n\) die Anzahl der Elemente darstellt. Sofern nicht anders vermerkt sind Vektoren immer Spaltenvektoren, während Zeilenvektoren dargestellt werden als:

\[\begin{align*} \bf x' &= \begin{pmatrix} x_{1} \phantom{x} x_{2} \phantom{x} \dots \phantom{x} x_{n} \phantom{x} \end{pmatrix} \end{align*}\]

Dabei zeigt das Apostroph an, dass der Vektor transponiert wurde, d.h. der Vektor nun einer Zeile entspricht.

Geometrisch sind Vektoren Geraden in einem \(p\)-dimensionalen Raum, die hinsichtlich ihrer Länge und Richtung, nicht aber ihrer Lage, spezifiziert sind. In einem 2-dimensionalen Raum kann der Vektor \[ \mathbf{x}=\begin{pmatrix} 2\\ 3 \end{pmatrix} \] wie folgt gezeichnet werden:

Wenn der Vektor seinen Ausgang im Ursprung hat, das heißt bei den Koordinaten \((0,0\)), dann spricht man von der Standardposition eines Vektors und in diesem Fall entsprechen die Werte eines Vektors den Koordinaten der Vektorspitze.

10.1.1 Rechenregeln für Vektoren

Die meisten Rechenregeln sind vermutlich intuitiv und werden – ohne großes Nachdenken – in R verwendet. Dennoch gibt es auch einige komplexere Fälle und wir beginnen daher mit einer formalen Einführung in die Rechenregeln für Vektoren.

Im Folgenden seien und zwei reelle Vektoren mit \(n\) Elementen und \(a\in\mathbb{R}\) eine beliebige reelle Zahl.

1. Skalarmultiplikation: Bei der Multiplikation mit einem Skalar \(a\) wird jedes Element einzeln mit \(a\) multipliziert. \[\begin{align*} a \cdot \bf x &= a \cdot \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} = \begin{pmatrix} a\cdot x_{1} \\ a \cdot x_{2} \\ \vdots \\ a \cdot x_{n} \end{pmatrix} \end{align*}\] Beispiel: \[\begin{align*} 5 \cdot \begin{pmatrix} 2 \\ 1 \\ 2.5 \end{pmatrix} = \begin{pmatrix} 5\cdot 2 \\ 5 \cdot 1 \\ 5 \cdot 2.5 \end{pmatrix} = \begin{pmatrix} 10 \\ 5 \\ 12.5 \end{pmatrix} \end{align*}\] Die Skalarmultiplikation kann geometrisch so verstanden werden, dass der Vektor neu skaliert wird, das heißt, in seiner Länge, nicht aber der Richtung verändert wird. Nehmen wir den Beispielvektor von gerade eben und skalieren ihn mit \(a=2\), so wird ersichtlich, dass die Länge nun verdoppelt wird (linke Abbildung). Ist der Skalar negativ, also z.B. \(a=-1.5\), wird der Vektor um 180° gedreht (rechte Abbildung). \[ 2\cdot\mathbf{x}=\begin{pmatrix} 2\\ 3 \end{pmatrix} =\begin{pmatrix} 4\\ 6 \end{pmatrix} \quad\text{und}\quad -1.5\cdot\mathbf{x}=\begin{pmatrix} 2\\ 3 \end{pmatrix} =\begin{pmatrix} -3\\ -4.5 \end{pmatrix} \]

2. Vektoraddition: Bei der Vektoraddition werden die Elemente beider Vektoren elementweise addiert. \[\begin{align*} \bf x + \bf y = \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{pmatrix} + \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{pmatrix} = \begin{pmatrix} x_{1} + y_{1} \\ x_{2} + y_{2} \\ \vdots \\ x_{n} + y_{n} \end{pmatrix} \end{align*}\] Beispiel: \[\begin{align*} \begin{pmatrix} 1 \\ 2.5 \\ 3 \end{pmatrix} + \begin{pmatrix} 4 \\ 3 \\ 1.2 \end{pmatrix} = \begin{pmatrix} 1 + 4 \\ 2.5 + 3 \\ 3 + 1.2 \end{pmatrix} = \begin{pmatrix} 5 \\ 5.5 \\ 4.2 \end{pmatrix} \end{align*}\] Geometrisch bedeutet die Vektoraddition, dass der Ausgangort des zweiten Vektors an die Spitze des ersten Vektors gelegt wird und der Summenvektor vom Ausgangsort des ersten zur Spitze des zweiten Vektors geht (die beiden schwarzen Pfeile in der linken Abbildung); das Ergebnis ist dann der blaue Pfeil (die roten Pfeile illustrieren die beschriebene Verknüpfung). Die Subtraktion zweier Vektoren kann so verstanden werden, dass der zweite Vektor zunächst mit \(-1\) skaliert und dann zum ersten Vektor addiert wird (roter Pfeil in der rechten Abbildung). Das Resultat (blauer Pfeil) ist das gleiche, als wenn man von der Spitze des zweiten Vektors zur Spitze des ersten Vektors gehen würde (gestrichelter blauer Pfeil): \[ \begin{pmatrix} 4\\ 2 \end{pmatrix} + \begin{pmatrix} 1\\ 2 \end{pmatrix} = \begin{pmatrix} 5\\ 4 \end{pmatrix} \quad\text{und}\quad \begin{pmatrix} 4\\ 2 \end{pmatrix} - \begin{pmatrix} 1\\ 2 \end{pmatrix} = \begin{pmatrix} 4\\ 2 \end{pmatrix} +(-1)\cdot \begin{pmatrix} 1\\ 2 \end{pmatrix} =\begin{pmatrix} 3\\ 0 \end{pmatrix} \]

3. Skalarprodukt: Das Skalarprodukt ist die Multiplikation zweier Vektoren. Dabei werden die Elemente an der gleichen Position multipliziert und aus allen Produkten wird die Summe gebildet. Das Skalarprodukt ist also eine (reelle) Zahl. \[\begin{align*} <\bf x, \bf y > &= \bf x' \cdot \bf y = \begin{pmatrix} x_{1}\phantom{x} x_{2} \phantom{x}\dots\phantom{x} x_{n} \end{pmatrix} \cdot \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \\ \end{pmatrix} = \it x_1 \cdot y_1 + x_2 \cdot y_2 + \dots + x_n \cdot y_n \\ &= \sum^n_{i = 1}x_i\cdot y_i \end{align*}\]

Beispiel: \[\begin{align*} <\begin{pmatrix} 1 \\ 2.5 \\ 3\end{pmatrix} , \begin{pmatrix}4 \\3 \\1.2\end{pmatrix}> &= \begin{pmatrix} 1 & 2.5 & 3\end{pmatrix} \cdot \begin{pmatrix}4 \\3 \\1.2\end{pmatrix}= 4 + 7.5 + 3.6 = 15.1 \end{align*}\]

4. Norm (bzw. Länge) eines Vektors: Die Norm ist die Länge eines Vektors und ergibt sich, indem man erst alle einzelnen Elemente quadriert, dann summiert und anschließend die Wurzel der Summe berechnet. In anderen Worten: Die Norm ist die Wurzel aus dem Skalarprodukt eines Vektors mit sich selber.

\[\begin{align*} \left\lVert\bf x\right\rVert = \sqrt{<\mathbf{x},\mathbf{x}>} = \sqrt{\sum^n_{i=1} x^2_i} \end{align*}\]

Beispiel: \[\begin{align*} \left\lVert\begin{pmatrix} 1 & 2.5 & 3\end{pmatrix}\right\rVert = \sqrt{1^2 + 2.5^2 + 3^2}=4.03 \end{align*}\]

10.1.2 Euklidische Distanz zweier Vektoren

Die euklidische Distanz zweier Vektoren \(d(\mathbf{x},\mathbf{y})\) kann leicht berechnet werden, wenn Subtraktion von Vektoren und die Norm gemeinsam betrachtet werden: Sie ist nichts anderes als die Norm des Differenzvektors.

\[ d(\mathbf{x},\mathbf{y}) = \left\lVert\mathbf{x}-\mathbf{y}\right\rVert = \sqrt{<\mathbf{x}-\mathbf{y},\mathbf{x}-\mathbf{y}>} = \sqrt{\sum^n_{i=1} (x_i-y_i)^2} \]

10.1.3 Winkel zwischen zwei Vektoren und Orthogonalität

Der Kosinus des Winkels \(\theta\), also \(\cos(\theta)\), zweier Vektoren bestimmt sich als:

\[ \cos(\theta) = \frac{<\bf x \cdot y>}{\bf \left\lVert x\right\rVert \cdot \left\lVert y\right\rVert} \]

Hieraus folgt, dass der Winkel \(\theta\) sich errechnet, indem man den Arkuskosinus auf den Bruch anwendet:

\[ \theta = \arccos\left(\frac{<\bf x \cdot y>}{\bf \left\lVert x\right\rVert \cdot \left\lVert y\right\rVert}\right) \]

Nun ergibt sich eine interessante Eigenschaft von Vektoren, wenn diese orthogonal, also rechtwinklig mit einem Winkel von 90, zueinander stehen. Denn der obige Ausdruck ergibt genau dann \(\theta = 90\), wenn \(\arccos(0)\) gilt, der Bruch innerhalb der Klammern des \(\arccos\) also insgesamt den Wert 0 annimmt.

Allgemein gilt somit: Ist das Skalarprodukt zweier Vektoren (der Zähler des Bruchs) 0, so sind die Vektoren auch gleichzeitig orthogonal zueinander.

Beispiel: Gegeben seien die orthogonal zueinanderstehenden Vektoren \(\mathbf{x} = (-3\ 1\ 5)\) und \(\mathbf{y} = (1\ 3\ 0)\). Der Winkel zwischen beiden Vektoren ergibt sich durch:

\[ \theta = \arccos\left(\frac{-3\cdot 1 + 1 \cdot 3 + 5 \cdot 0}{\sqrt{-3^2 + 1^2 + 5^2} \cdot \sqrt{1^2 + 3^2 + 0^2}}\right) = \arccos\left(\frac{0}{\sqrt{35} \cdot \sqrt{10}}\right) = \arccos(0) = \frac{\pi}{2}(=90^\circ) \] Ein Winkel wird zwar oft in Grad angegeben, er kann aber auch als Bogenmaß mit der Einheit Radiant (rad) angegeben werden. Da diese am Einheitskreis definiert sind, entspricht der volle Umfang eines Kreises von 360° dann genau \(2\cdot\pi\) rad. Durch Umformungen gelangt man zu den beiden Formeln zur Umrechnung, nämlich \(x\text{ rad}=x\cdot\frac{180}{\pi}\text{ Grad}\) bzw. \(y\text{ Grad}=y\cdot\frac{\pi}{180}\text{ rad}\). Üblicherweise arbeiten Funktionen in R mit der Einheit Radiant, \(\frac{\pi}{2}\) rad entsprechen dann 90\(^\circ\).

10.1.4 Variablen- vs. Personenraum

An dieser Stelle lohnt es sich, wenn wir uns von der bisher üblichen Darstellung von Daten etwas lösen. Dies wird uns in folgenden Teilen die Anschauung etwas vereinfachen. Bisher hatten wir die Variablen eines Datensatzes immer als Achsen eines Koordinatensystems aufgefasst und dann (z.B. für einen Scatterplot) einen Punkt für jede Versuchsperson eingezeichnet. Diese Darstellung wird als Variablenraum bezeichnet. Hier werden insbesondere Muster individueller Beobachtungen dargestellt.

Sollen hingegen die Variablen und ihre Zusammenhänge betont werden, so eignet sich eine andere Darstellung, bei der die Personen die Achsen darstellen und die Variablen als Vektoren in diesem sog. Personenraum aufgefasst werden.

Beide Betrachtungsweisen werden an einem Beispiel aus Wickens (2014) verdeutlicht. Die folgende Tabelle fasst die Werte zweier Personen \(P_1\) und \(P_2\) auf den Variablen \(X\) und \(Y\) zusammen:

\[ \begin{align} \begin{array}{c|cc} & X & Y \\\hline P_1 &-1 & 2 & \\ P_2 & 3 & 1 & \end{array} \end{align} \] Der Variablenraum in diesem Beispiel ist zweidimensional und besteht aus einer Achse \(X\) und einer Achse \(Y\) (Spalten der Matrix). Der Wert von Person 1 besitzt die Koordinaten \(P_1 = (-1, 2)\) und der Wert der Person 2 die Koordinaten \(P_2 = (3, 1)\). Dies ist im linken Teil der folgenden Abbildungen dargestellt. Der Personenraum ist ebenfalls zweidimensional und besteht aus einer Achse \(P_1\) und einer Achse \(P_2\) (Zeilen der Matrix). Der Wert der Variable \(X\) besitzt die Koordinaten \(X = (-1, 3)\) und der Wert der Variable \(Y\) die Koordinaten \(Y = (2,1)\). Ein Vektor vom Ursprung bis zum Punkt \((-1, 3)\) stellt die Variable \(X\) dar, und ein Vektor vom Ursprung bis zum Punkt \((2,1)\) stellt die Variable \(Y\) dar. Diese Situation ist im rechten Teil der folgenden Abbildung visualisiert.

Die Geometrie im zwei- bzw. dreidimensionalen Raum kann genutzt werden, um Prinzipien zu erfassen, die auch in höherdimensionalen Räumen gelten, die nicht leicht visualisierbar sind.

Wichtige Eigenschaften von Variablen ergeben sich daraus, wenn man sie als zentrierte Vektoren im Personenraum auffasst. Zentrierte Variablen erleichtern die Analyse der Beziehungen von Variablen zueinander und werden gewonnen, indem von jedem Element (eines Vektors) der Mittelwert aller Elemente subtrahiert wird. Um deutlich zu machen, dass wir von einem zentrierten Vektor sprechen, markieren wir diesen mit einem hochgestellten \(c\) und berechnen ihn als \[ {\bf x^c} = \bf x - M_{\bf x} \] wobei \(M_{\bf x}\) der Mittelwert aller Elemente ist. Durch diese Zentrierung verändern wir zwar die Lage im Raum, wenn wir die Vektoren zeichnen würden, nicht aber deren (Varianz und) Korrelation untereinander, hier verdeutlicht an einem Beispiel aus Wickens (2014):

# Ausgangsvariablen/-vektoren
X1 = c(1, 1, 3, 3, 4, 4, 5, 6, 6, 7)
X2 = c(4, 7, 9, 12, 11, 12, 17, 13, 18, 17)
cor(X1,X2)   

## [1] 0.9039935

# jetzt zentrieren
X1c = X1 - mean(X1)
X2c = X2 - mean(X2)
cor(X1c, X2c) # bleibt gleich

## [1] 0.9039935

Zudem gilt eine Beziehung zwischen der Länge eines (zentrierten) Vektors (d.h. dessen Norm) und der Standardabweichung der entsprechenden Variablen (welche gleich ist für die Ausgangsvariable und die zentrierte Variable bzw. deren Vektoren): \[ \left\lVert{\bf x^c}\right\rVert = \sqrt{N}\cdot s_{\bf x} \]

# zwei Funktionen die die Norm und die Standardabweichung berechnen
norm = function(x){sqrt(x%*%x)}
stddev = function(x){sqrt((1/length(x))*sum((x-mean(x))^2))}

# Vergleich
norm(X1c)

##          [,1]
## [1,] 6.164414

sqrt(10)*stddev(X1c)

## [1] 6.164414

Da also Norm und Standardabweichung proportional zueinander sind, ist ein (zentrierter) Vektor umso länger, je größer die Varianz der Ausgangsvariablen ist. Zudem sagt der Winkel zwischen zentrierten Vektoren etwas über deren Korrelation aus und es gilt \[ r_{\bf x_1x_2}=r_{\bf x^c_1x^c_2}=\frac{<\bf{\bf x_1^c\cdot \bf x_2^c}>}{\left\lVert\bf x_1^c\right\rVert \left\lVert\bf x_2^c\right\rVert} = \cos[\measuredangle(\bf x_1^c,\bf x_2^c)] \] Damit können wir aus den zentrierten Vektoren direkt die Korrelation (auch der Ausgangsvektoren) berechnen:

# Berechnung der Korrelation aus den zentrierten Vektoren
r = (X1c %*% X2c)/(norm(X1c)*norm(X2c))
r

##           [,1]
## [1,] 0.9039935

# daraus Berechnung des Winkels der zentrierten Vektoren (und Umwandlung in Grad)
acos(r)/pi*180

##          [,1]
## [1,] 25.31193

Insgesamt haben der Winkel zwischen den zentrierten Vektoren und deren Korrelation eine direkte Entsprechung: Im Fall perfekter Korrelationen, sind beide jeweils Linearkombinationen voneinander; sind beide unkorreliert, liegt ein Winkel von 90 Grad vor. Je kleiner der Winkel zwischen ihnen ist, umso größer ist ihre Korrelation.

Anmerkung 1: Die geometrische Anschauung verdeutlicht auch, dass nicht alle beliebigen Korrelationsmuster möglich sind (siehe Wickens, 2014, S. 20-21). Sind z.B. zwei Variablen \(X_1\) und \(X_2\) korreliert mit \(r_{X_1X_2} = 0.8\), dann kann für eine dritte Variable \(X_3\) und deren Korrelation mit \(X_1\) und \(X_2\) nicht gleichzeitig \(r_{X_1X_3} = 0.5\) und \(r_{X_2X_3} = -0.5\) gelten. Die Korrelation \(r_{X_1X_2} = 0.8\) bedeutet, dass die zentrierten Vektoren \(X_1^c\) und \(X_2^c\) einen Winkel von etwa 37 Grad haben. Der zentrierte Vektor der dritten Variablen kann aber nicht so platziert werden, dass er einen Winkel von 60 Grad zu \(X_1^c\) und gleichzeitig einen Winkel von 120 Grad zu \(X_2^c\) hat (die beiden Winkel die für die jeweiligen Korrelation gelten müssten). Dies wird aus den bloßen Zahlen, die die Korrelation beschreiben, nicht deutlich.

Anmerkung 2: Die mitunter benutzte geometrische Intuition, dass zwei Variablen deren Vektoren in einem Winkel von 90 Grad zueinander stehen auch unkorreliert sind, ist nicht ganz korrekt. Dies kann man sich leicht an einem Beispiel verdeutlichen:

X1 = c(4,0,0) # haben in einem 3D-Koordinatensystem einen 90 Grad Winkel zueinander
X2 = c(0,4,0)
cor(X1,X2)    # r = -0.5 (obwohl 90 Grad)

## [1] -0.5

# jetzt die zentrierten Vektoren
X1c = X1 - mean(X1)
X2c = X2 - mean(X2)
cor(X1c, X2c) # natürlich auch r = -0.5

## [1] -0.5

# Korrelation und Winkel mit den Ausgangsvektoren berechnet
( r = (X1 %*% X2)/(norm(X1)*norm(X2)) ) # r = 0 ??

##      [,1]
## [1,]    0

acos(r)/pi*180                          # 90 Grad (wie intuitiv)

##      [,1]
## [1,]   90

# jetzt mit den zentrierten Vektoren
( r = (X1c %*% X2c)/(norm(X1c)*norm(X2c)) ) # passt zu oben

##      [,1]
## [1,] -0.5

acos(r)/pi*180                          # aber kein 90 Grad Winkel!

##      [,1]
## [1,]  120

Die Beziehung zwischen Korrelation und Winkel zweier Variablen gilt nur nachdem die Variablen zentriert wurden. D.h. nur wenn die zentrierten Vektoren einen Winkel von 90 Grad haben, sind die (Ausgangs-)Variablen unkorreliert.

10.1.5 Lineare Unabhängigkeit bzw. Abhängigkeit

Eine Menge von \(p\) Vektoren \(\bf\{ x_1,\ x_2, \ \dots, \ x_p\}\) können linear abhängig oder linear unabhängig sein. Dabei heißen die Vektoren genau dann linear abhängig, wenn sich einer der Vektoren \(\bf x_i\) als Linearkombination der anderen Vektoren darstellen lässt. Mit anderen Worten: Eine Menge von Vektoren ist linear abhängig, wenn es reelle Zahlen \(a_1, \ \dots\ , \ a_{p-1}\) gibt, sodass gilt:

\[\begin{align*} {\bf x_p} = a_1 {\bf x_1} + a_2 {\bf x_2} + \dots + a_{p-1} {\bf x_{p-1}} \end{align*}\]

Existiert keine solche Kombination, sind die Vektoren der Menge per Definition linear unabhängig.

Als Beispiel sei die Menge \(\{(2.5,\ 4.5,\ 1.5)',\ (1,\ 3,\ 1)',\) \((4,\ 6,\ 2)'\}\) gegeben. Die hierin enthaltenen Vektoren sind linear abhängig, da sich der erste Vektor aus den letzten beiden Vektoren darstellen lässt:

\[\begin{align*} \begin{pmatrix} 2.5 \\ 4.5 \\ 1.5 \end{pmatrix} = \frac{1}{2} \cdot \begin{pmatrix} 1 \\ 3 \\ 1 \end{pmatrix} + \frac{1}{2} \cdot \begin{pmatrix} 4 \\ 6 \\ 2 \end{pmatrix} \end{align*}\]

Da lineare Abhängigkeit bedeutet, dass sich ein Vektor linear durch die übrigen Vektoren kombinieren lässt, ergibt sich auch unmittelbar eine geometrische Interpretation: Eine Menge von Vektoren ist genau dann linear abhängig, wenn einer der Vektoren im Raum, welcher von den übrigen Vektoren aufgespannt wird, liegt.

10.1.6 Vektorrechnung in R

In R haben wir bereits viel mit Vektoren gearbeitet, die sich mit Hilfe des c()-Operators erstellen lassen:

x <- c(1, 2, 3, 5)
x

## [1] 1 2 3 5

Wir werden nun einmal die Beispiele der oben eingeführten Rechenregeln in R nachvollziehen

1. Skalarmultiplikation:

x <- c(2, 1, 2.5)
5 * x

## [1] 10.0  5.0 12.5

2. Vektoraddition:

x <- c(1, 2.5, 3)
y <- c(4, 3, 1.2)
x + y

## [1] 5.0 5.5 4.2

3. Skalarprodukt:

Möchte man in R zwei Vektoren i.S. eines Skalarproduktes multiplizieren, so geht das leider nicht einfach mit dem naheliegenden *-Operator. In diesem Falle würde R alle Elemente der Vektoren “elementweise” multiplizieren und es resultiert wieder ein Vektor:

x <- c(1, 2.5, 3)
y <- c(4, 3, 1.2)
x*y

## [1] 4.0 7.5 3.6

Natürlich kann auch einfach die Summe des gerade berechneten Vektors berechnet werden…

sum(x*y)

## [1] 15.1

…aber es gibt für das Skalarprodukt (allgemeiner: für das Matrixprodukt, siehe weiter unten) auch einen eigenen Operator, der geschrieben wird als %*%:

x %*% y

##      [,1]
## [1,] 15.1

Vielleicht wurde bemerkt, dass der Vektor \(\bf x\) nicht transponiert werden musste, so wie oben eigentlich beschrieben wurde. Das ist wichtig, denn R unterscheidet als Default nicht zwischen Spalten- und Zeilenvektoren. Wird ein Vektor verwendet, so formt R ihn einfach passend um (was meistens dann auch passt). Transponieren muss man aber Vektoren (und vor allem Matrizen) dennoch manchmal und dies kann mit der Funktion t() getan werden:

a <- t(x) # dann "Zeilenvektor"
a

##      [,1] [,2] [,3]
## [1,]    1  2.5    3

b <- t(a) # jetzt wieder "Spaltenvektor"
b

##      [,1]
## [1,]  1.0
## [2,]  2.5
## [3,]  3.0

R betrachtet diesen Vektor dann übrigens als Matrix, wobei natürlich ein Vektor auch nur ein Spezialfall einer Matrix mit einer Zeile oder einer Spalte ist:

class(a)

## [1] "matrix" "array"

10.2 Matrizen

Matrizen sind in gewisser Weise eine Verallgemeinerung von Vektoren. Mehrere Spalten- oder Zeilenvektoren lassen sich in Form einer Tabelle mit \(n\) Zeilen und \(q\) Spalten zusammenfassen; einer (\(n\times q\))-Matrix (sprich: “n kreuz q”). Die Notation erfolgt als:

\[\begin{align*} \bf X = \begin{pmatrix} x_{11} & x_{12} & \dots & x_{1q} \\ x_{21} & x_{22} & \dots & x_{2q} \\ \vdots & \vdots & \dots & \vdots \\ x_{n1} & x_{n2} & \dots & x_{nq}\end{pmatrix} \end{align*}\]

Im Folgenden werden wir Matrizen großschreiben und fett drucken. Ein einzelnes Element indiziert man mit \(x_{ij}\) mit \(i \in \{1, \dots, n\}\) und \(j \in \{1, \dots, q\}\). Wenn eine Matrix genauso viele Zeilen wie Spalten besitzt (\(n\times n\)), dann bezeichnet man sie als quadratisch. Eine quadratische Matrix, die auf der Hauptdiagonalen Zahlen \(a_i\in\mathbb{R}\) enthält (die ungleich von Null sind), aber ansonsten nur Nullen außerhalb der Diagonalen, nennt man eine Diagonalmatrix. Hier ein Beispiel für eine \(3 \times 3\) Diagonalmatrix:

\[\begin{align*} \bf D = \begin{pmatrix} a_1 & 0 & 0 \\ 0 & a_2 & 0 \\ 0 & 0 & a_3 \end{pmatrix} \end{align*}\]

Sind zusätzlich all diese \(a_i=1\), so redet man von der Einheitsmatrix \(\bf I\). Manchmal wird als Index noch die Anzahl der Zeilen und Spalten dazugeschrieben, wie hier im Beispiel: \[\begin{align*} \bf I_3 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix} \end{align*}\]

10.2.1 Rechenregeln für Matrizen

Da Vektoren “Spezialfälle” von Matrizen sind, also \((n\times 1)\)- bzw. \((1\times q)\)-Matrizen sind, ähneln viele der Rechenregeln für Matrizen den Rechenregeln für Vektoren. Im Folgenden seien \(\bf X\) und \(\bf Z\) zwei \((n\times q)\)-Matrizen sowie \(a \in \mathbb{R}\) eine reelle Zahl.

1. Skalarmultiplikation: Jedes Element der Matrix wird mit \(a\) multipliziert.

\[\begin{align*} c\cdot \bf X =\begin{pmatrix} c\cdot x_{11} & c\cdot x_{12} & \dots & c\cdot x_{1q} \\ c\cdot x_{21} & c\cdot x_{22} & \dots & c\cdot x_{2q} \\ \vdots & \vdots & \vdots & \vdots \\ c\cdot x_{n1} & c\cdot x_{n2} & \dots & c\cdot x_{nq}\end{pmatrix} \end{align*}\] Beispiel: \[\begin{align*} 5 \cdot \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{pmatrix} = \begin{pmatrix} 5 & 10 \\ 15 & 20 \\ 25 & 30 \end{pmatrix} \end{align*}\]

2. Matrizenaddition: Die Elemente der Matrizen werden paarweise addiert.

\[\begin{align*} \bf X + Z &= \begin{pmatrix} x_{11} & x_{12} & \dots & x_{1q} \\ x_{21} & x_{22} & \dots & x_{2q} \\ \vdots & \vdots & \vdots & \vdots \\ x_{n1} & x_{n2} & \dots & x_{nq} \end{pmatrix}+ \begin{pmatrix} z_{11} & z_{12} & \dots & z_{1q} \\ z_{21} & z_{22} & \dots & z_{2q} \\ \vdots & \vdots & \vdots & \vdots \\ z_{n1} & z_{n2} & \dots & z_{nq} \end{pmatrix} \\ \\ &= \begin{pmatrix} x_{11} + z_{11} & x_{12} + z_{12} & \dots & x_{1q} + z_{1q} \\ x_{21} + z_{21} & x_{22} + z_{21} & \dots & x_{2q} + z_{2q} \\ \vdots & \vdots & \vdots & \vdots \\ x_{n1} + z_{n1} & x_{n2} + z_{n2} & \dots & x_{nq} + z_{nq}\end{pmatrix} \end{align*}\]

3. Transposition: Die Transposition ist auch eine Erweiterung dessen was wir für Vektoren bereits eingeführt haben. Im Wesentlichen ist gemeint, dass Zeilen und Spalten vertauscht werden. Mit anderen Worten: die erste Zeile wird zur ersten Spalte, die zweite Zeile wird die zweite Spalte usw.

Beispiel: \[\begin{align*} \begin{pmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{pmatrix}' = \begin{pmatrix} 1 & 3 & 5 \\ 2 & 4 & 6 \end{pmatrix} \end{align*}\]

Nun kommen wir zu zwei Operationen, die etwas komplizierter sind und mehr Verallgemeinerung benötigen:

4. Matrizenmultiplikation: Das Produkt zweier Matrizen ist wieder eine Matrix, deren Elemente die Skalarprodukte der Zeilen der ersten und der Spalten der zweiten Matrix sind. Sei nun \(\bf Y\) eine \((q\times m)\)-Matrix, dann lässt sich das Matrixprodukt aus \(\bf X\) und \(\bf Y\) formal darstellen als:

\[\begin{align*} \bf X \cdot Y = \begin{pmatrix} \sum_{k=1}^q x_{1k} y_{k1} & \dots & \sum_{k=1}^q x_{1k}y_{km} \\ \vdots & & \vdots \\ \sum_{k=1}^q x_{nk}y_{k1} & \dots & \sum_{k=1}^q x_{nk}y_{km} \end{pmatrix} \end{align*}\]

Die “Rechenregel” lautet dabei also in etwa “Zeile mal Spalte” (vgl. das unten folgende Beispiel zum einfacheren Verständnis). Dabei ist die Matrizenplutimikation nur dann möglich, wenn die Matrix \(\bf X\) genauso viele Spalten besitzt, wie die Matrix \(\bf Y\) Zeilen. Weiterhin gilt, dass das Ergebnis des Produkts einer \((n\times q)\)-Matrix mit einer \((q\times m)\)-Matrix eine Matrix mit \(n\) Zeilen und \(m\) Spalten ergibt (etwas informell: \((n\times q)\cdot (q\times m) =(n\times m)\)). Daraus folgt, dass die Matrizenmultiplikation, so wie sie definiert ist, nicht kommutativ ist, d.h. die beiden Matrizen nicht “ihre Reihenfolgen tauschen können” (wie dies bei der normalen Multiplikation möglich ist).

Im folgenden Beispiel ist zu sehen, wie aus einer \((3\times 2)\)-Matrix und einer \((2\times 3)\)-Matrix eine \((3\times 3)\)-Matrix resultiert. Die roten Felder zeigen an, welche Zeilen und Spalten miteinander verrechnet werden und dann ein Element des Resultats ergeben: \[\begin{align*} \begin{pmatrix} 1 & 2 \\ {\color{Red} 3} & {\color{Red} 4} \\ 5 & 6 \end{pmatrix} \cdot \begin{pmatrix} 5 & {\color{Red} 3} & 5 \\ 1 & {\color{Red} 2} & 2 \end{pmatrix} = \begin{pmatrix} 1\cdot5 + 2\cdot1 & 1 \cdot 3 + 2 \cdot 2 & 1\cdot5+2\cdot2 \\ 3\cdot5 + 4\cdot1 & \color{Red}{ 3\cdot 3 + 4\cdot 2} & 3\cdot5 + 4\cdot2 \\ 5\cdot5 + 6\cdot1 & 5\cdot3 + 6\cdot2 & 5\cdot5 + 6\cdot2 \end{pmatrix} \end{align*}\]

5. Matrixinversion: Die oben kurz erwähnte Einheitsmatrix spielt schließlich eine sehr wichtige Rolle für die Inverse einer Matrix (die wir später brauchen werden). Sei nun \(\bf A\) eine quadratische Matrix. Die Inverse zu \(\bf A\) ist diejenige Matrix \(\bf A^{-1}\), für die gilt: \[\begin{align*} \bf A \cdot A^{-1} = A^{-1} \cdot A = I \end{align*}\]

Das heißt, multipliziert man eine Matrix mit ihrer Inversen, ergibt sich die Einheitsmatrix. Dabei besitzt eine quadratische Matrix \(\bf A\) nur dann eine Inverse, wenn ihre Spalten- (oder auch Zeilenvektoren) linear unabhängig sind. Man sagt dann, die Matrix sei von vollem Rang. Eine Möglichkeit die Inverse einer Matrix zu bestimmen ist das Gauß-Jordan-Verfahren. Wir werden das Verfahren hier nicht formal einführen, sondern uns auf die Berechnung mit R beschränken.

10.2.2 Matrizenrechnung in R

Typischerweise erstellen wir eine Matrix indem wir die matrix()-Funktion nutzen oder mehrere Vektoren miteinander zeilen- oder spaltenweise verbinden. Dies wird im folgenden Code schrittweise demonstriert:

# 1. Elemente mit matrix() in Matrix überführen
elemente <- c(1, 2, 3, 4, 5, 6) # ein Vektor
mat <- matrix(elemente, 
              nrow = 2,         # ergibt eine 2 x... 
              ncol = 3,         # ...3 Matrix
              byrow = FALSE)    # die Elemente werden spaltenweise in die Matrix geschrieben
                                # Probieren Sie auch byrow = TRUE!
mat 

##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6

# 2. Vektoren zu Matrix zusammenfassen
a <- c(1, 2, 3)
b <- c(4, 5, 6)
mat <- rbind(a, b) # zeilenweise ("row-bind")
mat # man beachte, dass die Label "a" und "b" beibehalten werden

##   [,1] [,2] [,3]
## a    1    2    3
## b    4    5    6

mat <- cbind(a, b) # spaltenweise ("column-bind")
mat # man beachte, dass die Label "a" und "b" beibehalten werden

##      a b
## [1,] 1 4
## [2,] 2 5
## [3,] 3 6

1. Skalarmultiplikation: Funktioniert genau wie bei einem Vektor.

5 * mat

##       a  b
## [1,]  5 20
## [2,] 10 25
## [3,] 15 30

2. Matrizenaddition: R kann einerseits einen Skalar zu jedem Element der Matrix addieren…

mat_new <- mat + 2 # erhöht jedes Element um 2
mat_new

##      a b
## [1,] 3 6
## [2,] 4 7
## [3,] 5 8

…und natürlich zwei Matrizen mit gleich vielen Spalten und Zeilen elementweise addieren:

summe <- mat_new + mat
summe

##      a  b
## [1,] 4 10
## [2,] 6 12
## [3,] 8 14

3. Transposition: Funktioniert ganz ähnlich wie bei Vektoren, hier am Beispiel der zuletzt erstellten Matrix.

t(summe)

##   [,1] [,2] [,3]
## a    4    6    8
## b   10   12   14

4. Matrizenmultiplikation: Zwei Matrizen werden multipliziert, indem – wie auch bei Vektoren – der %*%- Operator verwendet wird.

elemente <- c(1, 2, 3,
              4, 5, 6) 
mat.X <- matrix(elemente,
                nrow = 2,
                ncol = 3,
                byrow = FALSE) 
mat.X # 2 x 3 - Matrix

##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6

mat.Y <- matrix(elemente,
                nrow = 3,
                ncol = 2, 
                byrow = FALSE) 
mat.Y # 3 x 2 - Matrix

##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6

mat.X %*% mat.Y # ...es resultiert eine 2 x 2 Matrix

##      [,1] [,2]
## [1,]   22   49
## [2,]   28   64

mat.Y %*% mat.X # ...es resultiert eine 3 x 3 Matrix

##      [,1] [,2] [,3]
## [1,]    9   19   29
## [2,]   12   26   40
## [3,]   15   33   51

5. Matrixinversion: Die Inverse einer Matrix kann mit der solve()-Funktion bestimmt werden.

# Aufpassen, Matrix muss von vollem Rang sein!
quad_mat <- matrix(c(1, 2, 10,
                     4, 5, -2,
                     7, 8, 5),
                   nrow = 3,
                   ncol = 3,
                   byrow = FALSE)
quad_mat

##      [,1] [,2] [,3]
## [1,]    1    4    7
## [2,]    2    5    8
## [3,]   10   -2    5

quad_mat.I <- solve(quad_mat)  # Bestimmung der Inversen
quad_mat.I

##            [,1]       [,2]        [,3]
## [1,] -0.7192982  0.5964912  0.05263158
## [2,] -1.2280702  1.1403509 -0.10526316
## [3,]  0.9473684 -0.7368421  0.05263158

round(quad_mat %*% quad_mat.I, 2) # Matrix mal Inverse ergibt die Einheitsmatrix

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

10.3 Determinanten, Eigenwerte und Eigenvektoren

10.3.1 Determinanten

10.3.1.1 Berechnung im allgemeinen Fall

Die Determinante einer Matrix \(\mathbf{A}\) ist letztlich eine einzige Zahl, in deren Berechnung sämtliche Werte einer Matrix eingehen. Für \((2\times 2)\)-Matrizen gibt es eine einfache Formel zur Berechnung. Wenn also eine Matrix \(\mathbf{A}\) beschrieben wird durch \[ \mathbf{A}= \begin{pmatrix} a & b \\ c & d \end{pmatrix}, \] dann wird ihre Determinante dadurch berechnet, dass das Produkt der Nebendiagonalen von dem der Hauptdiagonalen subtrahiert wird: \[|\mathbf{A}| = ad-bc\]

Nun betrachten wir \((3\times 3)\)-Matrizen der Form

\[ \mathbf{A}= \begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix} \] Die Determinante von \(\mathbf{A}\) wird dann berechnet als \[ |\mathbf{A}|:= \begin{vmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{vmatrix}:= a_{11}\cdot \begin{vmatrix} a_{22} & a_{23}\\ a_{32} & a_{33} \end{vmatrix}- a_{12}\cdot \begin{vmatrix} a_{21} & a_{23}\\ a_{31} & a_{33} \end{vmatrix}+ a_{13}\cdot \begin{vmatrix} a_{21} & a_{22}\\ a_{31} & a_{32} \end{vmatrix} \] Der Aufbau der Produktsumme auf der rechten Seite der Gleichung wird vielleicht deutlicher, wenn man sie in etwas anderer Form aufschreibt (wobei in den Matrizen rechts vom Gleichheitszeichen die nicht-fettgedruckten Einträge weggestrichen zu denken sind):

\[ \begin{vmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{vmatrix} :=a_{11}\cdot \begin{vmatrix} a_{11}&a_{12}&a_{13}\\ a_{21}&{\bf a_{22}}&{\bf a_{23}}\\ a_{31}&{\bf a_{32}}&{\bf a_{33}}\\ \end{vmatrix} -a_{12}\cdot \begin{vmatrix} a_{11}&a_{12}&a_{13}\\ {\bf a_{21}}&a_{22}&{\bf a_{23}}\\ {\bf a_{31}}&a_{32}&{\bf a_{33}}\\ \end{vmatrix} +a_{13}\cdot \begin{vmatrix} a_{11}&a_{12}&a_{13}\\ {\bf a_{21}}&{\bf a_{22}}&a_{23}\\ {\bf a_{31}}&{\bf a_{32}}&a_{33}\\ \end{vmatrix} \] Man sieht, dass beispielsweise die Determinante im ersten Summanden die Determinante jener Matrix ist, die aus \(\mathbf{A}\) durch Streichung der ersten Zeile und der ersten Spalte hervorgeht. Diese Determinante wird multipliziert mit dem Matrixelement, das im Schnittpunkt der gestrichenen Spalte und Zeile steht. Analoges gilt für die beiden anderen Summanden. Auf die spezifische Struktur der Vorzeichen in dieser Produktsumme wird weiter unten noch einmal eingegangen.

Bezeichnet man nun die Matrix, die aus \(\mathbf{A}\) durch Streichung der \(i\)-ten Zeile und \(j\)-ten Spalte hervorgeht, mit \(\mathbf{A_{ij}}\), so lässt sich auch kompakter schreiben: \[ |\mathbf{A}|=a_{11}\cdot|\mathbf{A_{11}}|-a_{12}\cdot|\mathbf{A_{12}}|+a_{13}\cdot|\mathbf{A_{13}}| \] Hier ist ein wichtiger Punkt, dass die Bestimmung von Determinanten für \((3\times 3)\)-Matrizen auf die Bestimmung von Determinanten für \((2\times 2)\)-Matrizen zurückgeführt wird. Rechnet man die Determinanten für die auftretenden \((2\times 2)\)-Matrizen explizit aus, so erhält man \[\begin{equation*} |\mathbf{A}|=a_{11}(a_{22}a_{33}-a_{23}a_{32})-a_{12}(a_{21}a_{33}-a_{23}a_{31}) +a_{13}(a_{21}a_{32}-a_{22}a_{31}) \end{equation*}\] und weiteres Ausrechnen liefert \[ |\mathbf{A}|=a_{11}a_{22}a_{33}+a_{12}a_{23}a_{31}+a_{13}a_{21}a_{32} -a_{11}a_{23}a_{32}-a_{12}a_{21}a_{33}-a_{13}a_{22}a_{31} \] Die letzte Gleichung wird auch als die Regel des Sarrus bezeichnet: Schreibt man sozusagen die erste und die zweite Spalte als vierte und fünfte Spalte der Matrix, ergibt sich die Determinante als die Summe der Produkte der drei Hauptdiagonalen minus die Summe der Produkte der drei Nebendiagonalen:

Die dargestellte Formel “entwickelt”, so sagt man, die Determinante nach der ersten Zeile. Natürlich kann man sie auch nach der zweiten oder dritten Zeile entwickeln, wobei dann die Formeln wie folgt aussehen (man achte auf die Vorzeichen!). Bei der Entwicklung nach der zweiten Zeile erhalten wir \[ |\mathbf{A}|=-a_{21}\cdot|\mathbf{A_{21}}|+a_{22}\cdot|\mathbf{A_{22}}|-a_{23}\cdot|\mathbf{A_{23}}| \] und für die Entwicklung nach der 3. Zeile resultiert: \[ |\mathbf{A}|=+a_{31}\cdot|\mathbf{A_{31}}|-a_{32}\cdot|\mathbf{A_{32}}|+a_{33}\cdot|\mathbf{A_{33}}| \] Dies kann man einfach durch Ausrechnen nachweisen. Rechnet man nämlich die Determinante nach diesen beiden Formeln aus, so kommt man jedes mal zur gleichen Produktsumme wie bei der Entwicklung nach der ersten Zeile (alle entsprechen also der Regel des Sarrus). Egal, welche der Formeln man also verwendet, man erhält jedes mal die Determinante von \(\mathbf{A}\). Die in den drei Formeln jeweils zu verwendenden Vorzeichen kann man sich leicht durch das folgenden Vorzeichenschema merken: \[\begin{equation*} \begin{pmatrix} + & - & +\\ - & + & -\\ + & - & + \end{pmatrix} \end{equation*}\]

Dabei erhält also ein Term der Form \(a_{ij}\cdot|\mathbf{A_{ij}}|\) ein negatives Vorzeichen, falls die Summe \(i+j\) ungerade, und ein positives Vorzeichen, falls sie gerade ist.

Man kann die Determinante \(\mathbf{A}\) auch nach einer Spalte entwickeln. Für die Entwicklung nach der ersten Spalte bekommt man die Formel (vgl. nun die Vorzeichen mit denen in der ersten Spalte) des Vorzeichenschemas):

\[\begin{equation*} |\mathbf{A}|=+a_{11}\cdot|\mathbf{A_{11}}|-a_{21}\cdot|\mathbf{A_{21}}|+a_{31}\cdot|\mathbf{A_{31}}| \end{equation*}\]

Einfaches Ausrechnen zeigt, dass dies zur gleichen Produktsumme wie die vorangegangenen Varianten führt, das heißt die Determinante \(|\mathbf{A}|\) ergibt. Die Formulierung der beiden Berechnungsformeln für die Entwicklung der Determinante nach der zweiten und dritten Spalte und der Nachweis, das sich mit diesen Formeln ebenfalls die Determinante berechnen lässt, wird analog durchgeführt.

Anzumerken bleibt noch, dass die Entscheidung, nach welcher Zeile oder Spalte die Determinante entwickelt werden sollte, von der jeweiligen Struktur der betrachteten Matrix abhängen wird. Besitzt beispielsweise eine der Spalten oder Zeilen besonders viele Nullen, dann wird es rechnerisch vorteilhaft sein, die Determinante nach dieser Spalte oder Zeile zu entwickeln. So entwickelt man die Determinante der Matrix \[ \mathbf{A}= \begin{pmatrix} \phantom{-}3 & 0 & 8 \\ \phantom{-}5 & 0 & 7\\ -1 & 4 & 2 \end{pmatrix} \] sinnvollerweise nach der zweiten Spalte: \[ |\mathbf{A}|= -0\cdot \begin{vmatrix} \phantom{-}5 & 7\\ -1 & 2 \end{vmatrix} + 0\cdot \begin{vmatrix} \phantom{-}3 & 8\\ -1 & 2 \end{vmatrix} -4\cdot \begin{vmatrix} 3 & 8\\ 5 & 7 \end{vmatrix} = -4\cdot \begin{vmatrix} 3 & 8\\ 5 & 7 \end{vmatrix} = 76 \]

Nun kommen wir zum allgemeinen Fall der Determinanten von \((n\times n)\)-Matrizen. Zuletzt hatten wir die Determinante einer \((3\times 3)\)-Matrix z.B. mit folgender Formel berechnet \[\begin{equation*} |\mathbf{A}|=a_{11}|\mathbf{A_{11}}|-a_{12}|\mathbf{A_{12}}|+a_{13}|\mathbf{A_{13}}| \end{equation*}\] Hierbei bezeichnete \(\mathbf{A_{ij}}\) jene Matrix, die aus \(\mathbf{A}\) durch Streichung der \(i\)-ten Zeile und \(j\)-ten Spalte hervorging. Weiter wurde festgehalten, dass das Vorzeichen eines Summanden auf der rechten Seite dieses Ausdrucks negativ ist, wenn die Summe der Indizes von \(\mathbf{A_{ij}}\) ungerade, und positiv ist, wenn die Summe gerade ist. Die Vorzeichen ergeben sich also durch den Ausdruck \[\begin{equation*} (-1)^{i+j} \end{equation*}\] Dadurch kann die Gleichung auch kompakter geschrieben werden als: \[\begin{equation*} |\mathbf{A}|=\sum_{j=1}^3 (-1)^{1+j} a_{1j}|\mathbf{A_{1j}}| \end{equation*}\]

Diese Gleichung zeigt nun den Weg, wie die Determinante allgemein für \((n\times n)\)-Matrizen zu definieren ist. Der Laplace’sche Entwicklungssatz besagt:

Sei \[\begin{equation*} \mathbf{A}= \begin{pmatrix} a_{11} & \dots & a_{1n}\\ \vdots & \ddots & \vdots\\ a_{n1} & \dots & a_{nn} \end{pmatrix} \end{equation*}\] eine \((n\times n)\)-Matrix. Dann ist die Determinante \(|\mathbf{A}|\) die Zahl \[\begin{equation*} |\mathbf{A}|:=\sum_{j=1}^n (-1)^{1+j} a_{1j}|\mathbf{A_{1j}}| \end{equation*}\] Für eine \((1\times1)\)-Matrix gilt: \(|\mathbf{A}|=(a_{11})=a_{11}\)

Diese Formel führt die Bestimmung der Determinanten einer \((n\times n)\)-Matrix \(\mathbf{A}\) auf die Bestimmung der Determinanten von \(n\) vielen \(((n-1)\times(n-1))\)-Matrizen \(\mathbf{A_{1j}}\) zurück. Wie aber bestimmt man die Determinanten dieser \(((n-1)\times(n-1))\)-Matrizen? Indem man den Laplace’schen Entwicklungssatz auf diese Matrizen anwendet. Dies wird sooft getan, bis man bei der Bestimmung der Determinanten von \((1\times 1)\)-Matrizen ankommt. Berechnungsformeln dieser Art nennt man auch rekursive Formeln, da man durch wiederholte Anwendung der Formel auf eine Situation zurückgeführt wird, die man zu behandeln weiß.

Die Determinante einer Matrix zu berechnen kann eine sehr aufwändige Sache sein. Es ist deshalb natürlich auch im allgemeinen Fall ratsam, sich die Matrix genau anzusehen, bevor man deren Determinante zu bestimmen beginnt. Besitzt die Matrix beispielsweise eine Zeile oder Spalte mit vielen Nullen, dann ist es in der Regel von Vorteil, die Determinante nach dieser Zeile bzw. Spalte zu entwickeln, da dann viele rechenaufwändige Summanden wegfallen können.

10.3.1.2 Die geometrische Interpretation

Gegeben seien die beiden Vektoren

\[ v_1 = \begin{pmatrix} 4 \\ 1 \end{pmatrix}\quad\text{und}\quad v_2 = \begin{pmatrix} 1 \\ 3 \end{pmatrix}\] Wir können nun beide Vektoren als Spalten einer Matrix \(\mathbf{X}\) auffassen. Der Betrag der Determinanten der Matrix entspricht dann dem Volumen des Parallelogramms, welches durch diese Vektoren aufgespannt wird, also: \[|\mathbf{X}|=\begin{pmatrix} 4 & 1 \\ 1 & 3 \end{pmatrix} = 11\] Auch mit R ergibt sich leicht das gleiche Ergebnis:

v1 <- c(4, 1)
v2 <- c(1, 3)
X <- cbind(v1,v2)
det(X)

## [1] 11

Das durch die Vektoren aufgespannte Parallelogramm ist in der folgenden Abbildung dargestellt. Wir wollen nun dessen Flächeninhalt einmal beispielhaft “per Hand” berechnen, um die Richtigkeit der obigen Behauptung zum Betrag der Determinanten und der Fläche zu testen.

Die Fläche \(A\) eines Parallelogramms berechnet sich als \[ A = h\cdot a \] wobei \(h\) die Höhe des Parallelogramms ist und \(a\) die dazugehörige Grundseite. In unserem Beispiel nehmen wir \(v_1\) als Grundseite und müssen nun die Höhe \(h\) ermitteln–das ist der etwas kompliziertere Teil.

Zur Berechnung von \(h\) brauchen wir zuerst den dazugehörigen Vektor sowie den Punkt \(b\) (beschrieben durch einen Vektor), an dem \(h\) den Vektor \(v_1\) berührt. Wir suchen also zuerst die orthogonale Projektion von \(v_2\) auf \(v_1\). Diese ergibt sich (ohne Beweis) wie folgt: \[b = \frac{< v_1, v_2 >}{<v_1,v_1>} \cdot v_1\] Berechnen wir \(b\) mit R ergibt sich also:

v1 <- c(4, 1)
v2 <- c(1, 3)
# in den aktuelleren Versionen von R muss das Skalarprodukt explizit in einen
# Vektor konvertiert werden, daher das as.vector() davor (es reicht aber auch c(...))
b <- as.vector(v1 %*% v2) / as.vector(v1 %*% v1) * v1 
b                                                   

## [1] 1.6470588 0.4117647

Und schließlich ergibt sich der Vektor für \(h\) durch die Differenz der Vektoren \(b\) und \(v_2\). Die Länge dieses Vektors wiederum ergibt dann letztlich unsere gesuchte Höhe.

h_vec <- v2 - b # vektor für h
h <- sqrt(sum(h_vec * h_vec)) # Länge des Vektors
h

## [1] 2.667892

Damit haben wir alles zusammen und können die Fläche des Parallelogramms ausrechnen:

a <- sqrt(sum(v1 * v1)) # Länge der Grundseite = Länge des Vektors v1
A <- h * a
A # ergibt tatsächlich 11

## [1] 11

10.3.2 Eigenwerte und Eigenvektoren

Formal werden Eigenwerte und Eigenvektoren einer Matrix wie folgt definiert:

1. Sei \(\mathbf{A}\) eine quadratische \((n\times n)\)-Matrix. Dann heißen alle \(\lambda\), die die charakteristische Gleichung \[|\mathbf{A}-\lambda\cdot \mathbf{I}|=0\] erfüllen, Eigenwerte von \(\mathbf{A}\).

2. Sei \(\mathbf{A}\) eine quadratische \((n\times n)\)-Matrix. Dann heißen alle Vektoren \[\mathbf{x}=\begin{pmatrix}x_1\\\vdots\\x_n\end{pmatrix}\neq\begin{pmatrix}0\\\vdots\\0\end{pmatrix}\] die die Gleichung \[\mathbf{A}\cdot\mathbf{x}=\lambda\cdot\mathbf{x}\] für eine reelle Zahl \(\lambda\) erfüllen, Eigenvektoren von \(\mathbf{A}\) zum Eigenwert \(\lambda\) (\(\mathbf{A}\cdot\mathbf{x}\) und \(\lambda\cdot\mathbf{x}\) sind dabei \(n\)-dimensionale Vektoren). Normalerweise wird jeder Vektor bei Multiplikation mit der Matrix \(\mathbf{A}\) gestreckt und gedreht. Eigenvektoren sind genau diejenigen Vektoren, die bei Multiplikation mit \(\mathbf{A}\) zwar um \(\lambda\) gestreckt, nicht aber gedreht werden.

Bevor wir zu einem Beispiel zur Berechnung kommen, sind im Folgenden wichtige Eigenschaften von Eigenwerten und Eigenvektoren aufgelistet (\(\mathbf{A}\) und \(\mathbf{B}\) seien quadratische, invertierbare (sog. “reguläre”) Matrizen):

• \(\lambda\) ist Eigenwert zu \(\mathbf{A}\Leftrightarrow \lambda^{-1}\) ist Eigenwert zu \(\mathbf{A}^{-1}\)
• \(\lambda\) ist Eigenwert zu \(\mathbf{A}\Leftrightarrow \lambda\) ist Eigenwert zu \(\mathbf{B}^{-1}\mathbf{A}\mathbf{B}\)
• \(\lambda\) ist Eigenwert zu \(\mathbf{A}\mathbf{B}\Leftrightarrow \lambda\) ist Eigenwert zu \(\mathbf{B}\mathbf{A}\)
• \(\lambda\) ist Eigenwert zu \(\mathbf{A}\Leftrightarrow \lambda+1\) ist Eigenwert zu \(\mathbf{A}+\mathbf{I}\)
• Die Eigenwerte einer Diagonalmatrix sind gerade ihre Hauptdiagonalelemente.
• Die Eigenvektoren einer symmetrischen Matrix sind orthogonal zueinander.

Wenngleich Eigenwerte und Eigenvektoren mit der R-Funktion eigen() sehr schnell berechnet werden können, folgt hier zunächst ein kleines Beispiel zur Berechnung im Fall einer \((2\times 2)\)-Matrix. Gesucht sind also die Eigenwerte und Eigenvektoren zur Matrix \[\mathbf{A}=\begin{pmatrix}1 & 0\\1 & 3\end{pmatrix}\] Wir beginnen mit den Eigenwerten: \[\begin{equation*} \begin{split} &|\mathbf{A}-\lambda\cdot\mathbf{I}|=0 \\ \Leftrightarrow & \left|\begin{pmatrix}1 & 0\\1 & 3\end{pmatrix}-\lambda\cdot \begin{pmatrix}1&0\\0&1\end{pmatrix}\right|=0\\ \Leftrightarrow & \left|\begin{pmatrix}1-\lambda&0\\1&3-\lambda\end{pmatrix}\right|=0 \\ \Leftrightarrow & (1-\lambda)(3-\lambda)=0\\ \Rightarrow & \lambda_1=3\quad\text{und}\quad\lambda_2=1 \end{split} \end{equation*}\]

Nun haben wir die zwei Eigenwerte der Matrix bestimmt und können zu jedem Eigenwert einen Eigenvektor berechnen. Wir beginnen mit \(\lambda_1\): \[\begin{equation*}\begin{split} &\mathbf{A}\cdot\mathbf{x}=\lambda_1\cdot\mathbf{x}\\ \Rightarrow&\begin{pmatrix}1 & 0\\1 & 3\end{pmatrix}\begin{pmatrix}x_1\\x_2 \end{pmatrix}=3\cdot\begin{pmatrix}x_1\\x_2 \end{pmatrix}\\ \Rightarrow&\quad\begin{matrix}x_1&&&=&3x_1\\x_1&+&3x_2&=&3x_2\end{matrix}\\ \Rightarrow&\begin{pmatrix}x_1\\x_2\end{pmatrix}=n\cdot \begin{pmatrix}0\\1\end{pmatrix}\quad(n\in\mathbb{R}) \end{split} \end{equation*}\] Auf die gleiche Art berechnen wir den Eigenvektor zu \(\lambda_2\): \[\begin{equation*}\begin{split} &\mathbf{A}\cdot\mathbf{x}=\lambda_2\cdot\mathbf{x}\\ \Rightarrow&\begin{pmatrix}1 & 0\\1 & 3\end{pmatrix}\begin{pmatrix}x_1\\x_2 \end{pmatrix}=1\cdot\begin{pmatrix}x_1\\x_2 \end{pmatrix}\\ \Rightarrow&\quad\begin{matrix}x_1&&&=&x_1\\x_1&+&3x_2&=&x_2\end{matrix}\\ \Rightarrow&\begin{pmatrix}x_1\\x_2\end{pmatrix}=n\cdot \begin{pmatrix}-2\\\phantom{-}1\end{pmatrix}\quad(n\in\mathbb{R}) \end{split} \end{equation*}\]

Hier wird noch deutlich, dass Eigenvektoren nur bis auf einen Skalar \(n\) genau bestimmbar sind, d.h. sie eigentlich nur in der Richtung und nicht in der Länge eindeutig definiert sind.

In R gibt es die Funktion eigen() um Eigenwerte und -vektoren zu berechnen:

A <- rbind(c(1, 0),
           c(1, 3))
eigen.results <- eigen(A)
eigen.results

## eigen() decomposition
## $values
## [1] 3 1
## 
## $vectors
##      [,1]       [,2]
## [1,]    0  0.8944272
## [2,]    1 -0.4472136

Während die Eigenwerte den händisch berechneten Eigenwerten entsprechen, schauen die Eigenvektoren – hier der zweite – auf den ersten Blick anders aus. Ursache dafür ist, dass die Funktion eigen() normierte Eigenvektoren ausgibt, das heißt, Vektoren mit einer Länge von 1. Durch Wahl eines geeigneten Wertes für \(n\) können wir die normierten Vektoren aber quasi auch zurückrechnen in die händisch bestimmten: \[ \begin{aligned} &n\cdot 0.8944272 = -2 \\ \Leftrightarrow &n = \frac{-2}{0.8944272}\\ \Leftrightarrow &n = -2.236068 \end{aligned} \]

n <- -2 / eigen.results$vectors[1,2]
n

## [1] -2.236068

eigen.results$vectors[,2] * n    # entspricht dem händisch berechneten Eigenvektor

## [1] -2  1

Bisher haben wir beliebige (quadratische) Matrizen betrachtet. Eine besondere Form sind symmetrische Matrizen, bei denen das obere Dreieck an Werten und das untere Dreieck an Werten quasi an der Hauptdiagonalen gespiegelt sind. Zwei Beispiele hierfür sind: \[ \mathbf{A}= \begin{pmatrix} 4 & 1 \\ 1 & 2 \end{pmatrix} \quad\text{oder}\quad \mathbf{B}= \begin{pmatrix} 2 & 1 & 3 \\ 1 & 4 & 2 \\ 3 & 2 & 3 \end{pmatrix} \] Hierbei haben die Eigenvektoren der Matrizen noch die Eigenschaften, dass sie orthogonal aufeinander stehen, wie am Kreuzprodukt leicht gesehen werden kann:

# Beispiel 1
A <- rbind(c(4, 1),
           c(1, 2))
(e.vector <- eigen(A)$vectors)

##            [,1]       [,2]
## [1,] -0.9238795  0.3826834
## [2,] -0.3826834 -0.9238795

e.vector[,1] %*% e.vector[,2]

##      [,1]
## [1,]    0

# Beispiel 2
B <- rbind(c(2, 1, 3),
           c(1, 4, 2),
           c(3, 2, 3))
(e.vector <- eigen(B)$vectors)

##            [,1]       [,2]       [,3]
## [1,] -0.4954257  0.4800681  0.7239393
## [2,] -0.5795168 -0.8034930  0.1362322
## [3,] -0.6470809  0.3520421 -0.6762786

round(e.vector[,1] %*% e.vector[,2], 2)

##      [,1]
## [1,]    0

round(e.vector[,1] %*% e.vector[,3], 2)

##      [,1]
## [1,]    0

round(e.vector[,2] %*% e.vector[,3], 2)

##      [,1]
## [1,]    0

Wir werden sehr bald eine typische Matrix kennenlernen und benutzen (der wir aber in Statistik I tatsächlich schon begegnet sind), die immer symmetrisch ist und für die daher die Orthogonalität der Eigenvektoren gilt.

10.4 Literatur

Wickens, T. D. (2014). The geometry of multivariate statistics. Psychology Press.

11 Multivariate Statistik

Bisher haben wir nur Verfahren mit einer abhängigen Variablen betrachtet. Diese Verfahren heißen univariat. Multivariate Verfahren generalisieren diese nun auf mehr als eine abhängige Variable.

# Pakete die in diesem Teil benutzt werden:
library(car)
library(mvtnorm)
library(DescTools)
library(schoRsch)
library(expm)

11.1 Daten, Mittelwertvektor, (Varianz-)Kovarianzmatrix

Im Folgenden werden Datenformate sowie Kenngrößen von Variablen im uni- und im multivariaten Fall einander gegenübergestellt:

11.1.1 Datenformate

Im univariaten Fall wurden Daten als Datenreihe einer Variablen \(X\) angegeben. Im multivariaten Fall wird hierzu eine Datenmatrix verwendet, in der üblicherweise die \(p\)-vielen Variablen in Spalten dargestellt werden und jede der \(n\)-vielen Beobachtungseinheiten (z.B. Versuchspersonen) in Zeilen:

11.1.2 Zentrale Tendenz

Das wichtigste Maß der zentralen Tendenz im univariaten Fall war der Mittelwert \(M_X\). Im multivariaten Fall ist dies ein Vektor \(\mathbf{\bar{x}}\) mit \(p\)-vielen Einträgen, dessen Elemente jeweils die Mittelwerte der \(p\)-vielen Variablen sind. Dieser Vektor wird auch Zentroid genannt:

11.1.3 Datenvariabilität

Das wichtigste Maß der Variabilität war im univariaten Fall die Varianz \(S^2_X\) (sowie die daraus abgeleiteten Größen). Die einfache Generalisierung für den multivariaten Fall ist die (Varianz-)Kovarianzmatrix \(\mathbf{S}\) (im Folgenden nur noch als Kovarianzmatrix bezeichnet).

11.1.4 Beispiel

Für die Datenmatrix \[\mathbf{X}=\begin{pmatrix}2&1\\4&3\\1&1\\5&2\\3&3\end{pmatrix}\] ergibt sich als Mittelwertsvektor \[\bar{\mathbf{x}}=\begin{pmatrix}3\\2\end{pmatrix}\] und als Kovarianzmatrix (\(S_{ij}:=\text{Kov}(X_i,X_j)\)): \[\mathbf{S}=\begin{pmatrix}S_{11}&S_{12}\\S_{21}&S_{22}\end{pmatrix}= \begin{pmatrix}2&0.8\\0.8&0.8\end{pmatrix}\] Die Korrelation zwischen \(X_1\) und \(X_2\) beträgt dann \(r=\frac{\text{Kov}(X_1,X_2)}{S_{X_1}S_{X_2}}=\frac{0.8}{\sqrt{2}\cdot\sqrt{0.8}}=.63\).

11.1.5 Geometrische Darstellung (im Variablenraum)

Im Variablenraum, das heißt, wenn die Variablen die Achsen bilden und wir uns für die individuellen Beobachtungen interessieren, werden Daten in Form einer Punktwolke dargestellt (vgl. Einführung Kovarianz in Teil 5 von Statistik I). Die einzelnen Versuchspersonen werden typischerweise als Punkte dargestellt (die eigentlich eine Vektorspitze repräsentieren). Um diese Punktwolken kann, wie später noch ausgeführt wird, ein Ellipsoid \(\varepsilon\) gelegt werden, der eine solche Form hat, dass bei minimalem Volumen maximal viele Punkte in ihm enthalten sind. Zwei der Parameter, die in dessen Berechnung eingehen, sind der Zentroid \(\mathbf{\bar{x}}\) (bestimmt den Mittelpunkt) sowie die Kovarianzmatrix \(\mathbf{S}\) (bestimmt die Form und Größe). Ein dritter Parameter \(k\), der nicht in den Daten enthalten ist, skaliert die Größe (d.h., den Radius im Sinne der sog. Mahalanobis-Distanz; vgl. Abschnitt 12.2.3).

Die folgende Abbildung zeigt zwei Beispiele von Punktwolken zweier Variablen \(X_1\) und \(X_2\) mit entsprechend eingezeichnetem Ellipsoiden. Im linken Beispiel haben beide die gleiche Varianz und korrelieren positiv; im rechten Beispiel sind beide Variablen unkorreliert, sie haben aber verschiedene Varianzen. Auf die eingezeichneten Linien wird gleich eingegangen.

11.2 Das Ellipsoid und multivariate Analoga zur Varianz

11.2.1 Das Ellipsoid \(\varepsilon\)

Während die analytische Darstellung der Kovarianzstruktur in Form der Matrix \(\mathbf{S}\) erfolgt, kann das bereits oben erwähnte Ellipsoid \(\varepsilon\) zur geometrischen Darstellung der Kovarianzmatrix \(\mathbf{S}\) herangezogen werden. In die exakte Darstellung gehen insbesondere Eigenwerte und Eigenvektoren von \(\mathbf{S}\) mit ein.

Die allgemeine Definition des (Volumens des) Ellipsoiden lautet: \[ \varepsilon(\mathbf{\bar{x}},\mathbf{S}, k^2):= \{\mathbf{x}\in\mathbb{R}^p|(\mathbf{x}-\mathbf{\bar{x}})' \cdot\mathbf{S}^{-1}\cdot (\mathbf{x}-\mathbf{\bar{x}})\leq k^2 \} \]

Umgangssprachlich bedeutet dies in etwa: Die Menge aller Punkte, die innerhalb eines \(p\)-dimensionalen Ellipsoiden mit dem “Radius” \(k\) liegen, wobei üblicherweise das Ellipsoid mit dem Radius \(k=1\) gezeichnet wird (und wir das im Folgenden auch so tun, wenn es nicht anders vermerkt ist).

Gezeichnet werden kann das Ellipsoid relativ leicht, nämlich mithilfe der Eigenvektoren und Eigenwerte der Kovarianzmatrix \(\mathbf{S}\). Wenn die Vektoren \(\mathbf{e_1},\ldots,\mathbf{e_p}\) die Eigenvektoren der Matrix \(\mathbf{S}\) sind, dann sind dies genau die Hauptachsen des Ellipsoiden, und die Wurzeln der Eigenwerte, also \(\sqrt{\lambda_1},\ldots,\sqrt{\lambda_p}\) sind die halben Längen der jeweiligen Hauptachsen (wenn \(k=1\)). Dies ist in der folgenden Abbildung für den Fall \(p=2\) illustriert.

Eine naheliegende Frage ist: Wieviele Datenpunkte liegen innerhalb bzw. außerhalb des Ellipsoiden? Natürlich kann dies nur relativ beantwortet werden, aber es gilt, dass der relative Anteil der Datenpunkte, die außerhalb des um den Faktor \(k\) vergrößerten Ellipsoiden liegen, höchstens \(\frac{p}{k^2}\) beträgt. Anders ausgedrückt: Innerhalb des \(k\)-fachen Ellipsoiden \(\varepsilon(\mathbf{S}, \mathbf{\bar{x}},k^2)\) liegen folglich mindestens \((1-\frac{p}{k^2})\cdot 100\%\) der Datenpunkte. Für manche Konstellationen von \(p\) und \(k\) ergeben sich hier allerdings wenig hilfreiche Werte. So liegen nämlich z.B. für \(p=2\) und \(k=1\) \(-100\%\) der Datenpunkte innerhalb des Ellipsoiden, allerdings ergibt sich mit \(k=p=2\) direkt ein sinnvollerer Wert von \(50\%\).

Schließlich wird für \(p=1\) (also dem univariaten Fall) das Ellipsoid mit Radius \(k=1\) eine Strecke mit der Länge \(2\cdot S_X\) und dem Mittelpunkt \(M_X\). In diesem Fall gilt allgemein \[P(|x-M_X|\geq k\cdot S_X)\leq \frac{1}{k^2}\] Dies ist die – oft nur im univariaten Fall behandelte – Tschebyscheff’sche Ungleichung, die hier als Spezialfall resultieren würde.

11.2.2 Multivariate Analoga zur Varianz

Detaillierte Informationen über die multivariate Variabilität sind also in der Kovarianzmatrix \(\mathbf{S}\) enthalten und können vermöge des Ellipsoiden im Variablenraum visualisiert werden. Die Eigenwerte von \(\mathbf{S}\) geben auch Auskunft über die Form des Ellipsoiden und sie gehen oft – manchmal indirekt über die Determinante oder die sog. Spur – in die Berechnung multivariater Prüfgrößen mit ein.

Ist es hingegen sinnvoll oder erforderlich, einen einzigen Kennwert für die Datenvariabilität anzugeben, gibt es zwei Konzepte.

Die verallgemeinerte Varianz (engl. generalized variance) ist die Determinante der Matrix \(\mathbf{S}\), also \(\det S\). Sie enthält Informationen über “die Kombination der Varianzen und Kovarianzen”, das spezifische Gewicht geht aber verloren. Sie ist…

• …proportional zum Volumen des Ellipsoiden und…
• …proportional zum Volumen des Parallelepipeds, welches durch die Spalten der Kovarianzmatrix \(\mathbf{S}\) aufgespannt wird.

Die verallgemeinerte Varianz wird kleiner, wenn die Kovarianzen (betragsmäßig) größer oder wenn die Varianzen kleiner werden. Kovarianz- und Varianzinformation sind also hier konfundiert. Die folgende Abbildung zeigt vier Beispiele für Kovarianzstrukturen im Variablenraum (obere Reihe) und den Parallelepipeden der Spalten der dazugehörigen Kovarianzmatrix (untere Reihe). Dabei wurden die Daten so konstruiert, dass die folgenden Kovarianzmatrizen \(\mathbf{S}\) resultieren:

\[ \mathbf{S}= \begin{pmatrix} 1 & 0.2 \\ 0.2 & 1 \end{pmatrix}\qquad \mathbf{S}= \begin{pmatrix} 1 & 0.7 \\ 0.7 & 1 \end{pmatrix}\qquad \mathbf{S}= \begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix}\qquad \mathbf{S}= \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \]

Die Gesamtvarianz (engl. total variance) ist die sogenannte Spur der Matrix \[ \text{Gesamtvarianz}:=\sum_{j=1}^p S_{jj}\quad\text{oder}\quad\sum_{j=1}^p \lambda_{j} \] In ihr ist keinerlei Information über Richtung und Struktur der Kovariation mehr enthalten, aber sie kann interpretiert werden als der “mittlere, quadrierte Abstand der Datenpunkte vom Zentroiden”.

11.2.3 Mahalanobis-Transformation und -Distanz

In der univariaten Statistik hatten wir durch die \(z\)-Transformation \[Z=\frac{X-M_X}{S_X}\] eine Datenreihe so transformiert, dass \(M_Z=0\) und \(S_Z=1\) gilt. Die multivariate Verallgemeinerung ist die Mahalanobis-Transformation.

11.2.3.1 Mahalanobis-Transformation

Werden Daten Mahalanobis-transformiert, so erfüllen sie drei Eigenschaften:

• Der neue Zentroid ist der Nullvektor \(\mathbf{0}\).
• Alle Varianzen/Standardabweichungen werden 1, das heißt \(S_{ii}=1\quad\forall i\in\{1,\ldots,p\}\).
• Alle Kovarianzen werden Null, das heißt \(S_{ij}=0\quad\forall i,j\in\{1,\ldots,p\}\text{ mit }i\neq j\).

Aus den letzten beiden Eigenschaften folgt, dass die neue Kovarianzmatrix die Einheitsmatrix \(\mathbf{I_p}\) ist. Dies bedeutet auch, dass die Transformation aus dem Ellipsoiden eine Einheitskugel macht.

Diese drei Eigenschaften werden erfüllt durch die Transformation \[ \mathbf{z}=\mathbf{S^{-\frac{1}{2}}}\cdot (\mathbf{x}-\mathbf{\bar{x}}) \] wobei \(\mathbf{S^{-\frac{1}{2}}}\) diejenige Matrix ist, für die gilt: \[ (\mathbf{S^{-\frac{1}{2}}})'\cdot \mathbf{S^{-\frac{1}{2}}}=\mathbf{S^{-1}} \] Diese Matrix ist also eine Art Quadratwurzel der Inversen der Matrix \(\mathbf{S}\). In R kann sie mit einer Kombination von solve() und der Funktion sqrtm() aus dem Paket expm bestimmt werden (für mehr Informationen, siehe hier).

Wir benutzen im folgenden Code die Daten aus Abschnitt 12.1.5 (linke Abbildung), die bspw. wie folgt erzeugt werden können:

set.seed(6)                  # damit das Ergebnis der Zufallsziehung
                             # reproduzierbar ist
n <- 100                     # Stichprobengröße
mu <- c(0, 0)                # Zentroid der Population
Sigma <- rbind(c(1,   0.5),
               c(0.5, 1))    # Kovarianzmatrix der Population
dat <- rmvnorm(n, mu, Sigma) # Daten ziehen (aus multivariater Normalverteilung)

Zur Visualisierung des Ellipsoiden benutzen wir die Funktion ellipse() aus dem Paket car. Die gestrichelten Linien in der rechten Abbildung illustrieren, dass das Ellipsoid mit Radius = 1 (also \(k=1\) bei mahalanobis-transformierten Daten tatsächlich den Einheitskreis ergibt.

par(mfrow = c(1,2))         # zwei Abbildungen nebeneinander
par(pty = "s")              # wichtig, damit beide Achsen gleich skaliert werden

mean <- colMeans(dat)       # empirischer Mittelwertvektor
Sigma <- (n-1)/n * cov(dat) # empirische Kovarianzmatrix

plot(dat[,1],dat[,2], 
     asp = 1,               # skaliert beide Achsen dann gleich
     xlim = c(-5, 5),
     ylim = c(-5, 5),
     main = "Originaldaten", 
     xlab = expression(italic(X[1])),
     ylab = expression(italic(X[2])), 
     col = "black",
     pch = 21, 
     cex = 0.5,
     cex.main = 0.8,
     bg = "lightblue")      # Datenpunkte

ellipse(mean,
        Sigma,
        1,
        col = "red")

# Funktion, um Mahalanobis-Transformation für eine Person (einen Zeilenvektor)
# durchzuführen
maha_trans <- function(x, mu, Sigma) {
  x <- as.matrix(x) # notwendig, um Matrizenoperationen durchzuführen
  Sigma.m12 <- sqrtm(solve(Sigma))   # Sigma^(-1/2) berechnen
  Sigma.m12 %*% (x - mu)             # Mahalanobis-Transformation
}

z.dat <- apply(X = dat,              # Funktion für die Transformation aufrufen
               MARGIN = 1,           # für jede Zeile, d.h. Person
               FUN = maha_trans, 
               mu = mean,
               Sigma = Sigma)

z.dat <- as.data.frame(t(z.dat))     # als DataFrame (wenn nötig)
names(z.dat) <- c("X1", "X2")        # ggf. noch Variablennamen benennen

plot(z.dat[,1],z.dat[,2], 
     asp = 1,
     xlim = c(-5, 5),
     ylim = c(-5, 5),
     main = "Mahalanobis-transformierte Daten",
     xlab = expression(italic(X[1])),
     ylab = expression(italic(X[2])), 
     col = "black",
     pch = 21, 
     cex = 0.5,
     cex.main = 0.8,
     bg = "lightblue") # Datenpunkte

ellipse(colMeans(z.dat),
        (n-1)/n * cov(z.dat),
        1,
        col = "red")
abline(h = c(1,-1),
       v = c(1,-1),
       lty = 2)

# Zum Check:
round(colMeans(z.dat), 2)       # Mittelwertvektor

## X1 X2 
##  0  0

round((n-1)/n * cov(z.dat), 2)  # (unkorrigierte) Kovarianzmatrix

##    X1 X2
## X1  1  0
## X2  0  1

11.2.3.2 Mahalanobis-Distanz

Gegenüber der euklidischen Distanz relativiert die Mahalanobis-Distanz die Entfernung zweier Punkte durch Ausrechnung der Kovarianzstruktur: In einer längeren Richtung einer Ellipse zählt Abstand weniger, als in einer kürzeren Richtung.

Die Mahalanobis-Distanz zweier Vektoren \(\mathbf{x_1}\) und \(\mathbf{x_2}\) beträgt: \[MD(\mathbf{x_1},\mathbf{x_2}):=\sqrt{(\mathbf{x_2}-\mathbf{x_1})'\cdot \mathbf{S}^{-1}\cdot(\mathbf{x_2}-\mathbf{x_1})}\]

Würde man hier statt \(\mathbf{S}^{-1}\) die Einheitsmatrix \(\mathbf{I_p}\) verwenden, das heißt, würde man die Gleichung ohne \(\mathbf{S}^{-1}\) aufschreiben, so ergibt sich die euklidische Distanz.

Zur Herleitung der Mahalanobis-Distanz seien \(\mathbf{z_1}\) und \(\mathbf{z_2}\) zwei mahalanobis-transformierte Vektoren. Dann ist deren quadrierte, euklidische Distanz \(||\mathbf{z_2}-\mathbf{z_1}||^2\) und es gilt weiter: \[ \begin{aligned} ||\mathbf{z_2}-\mathbf{z_1}||^2&=<\mathbf{z_2}-\mathbf{z_1},\mathbf{z_2}-\mathbf{z_1}>\\ &=(\mathbf{z_2}-\mathbf{z_1})'\cdot(\mathbf{z_2}-\mathbf{z_1})\\ &=\left[ \mathbf{S}^{-\frac{1}{2}}\cdot(\mathbf{x_2}-\bar{\mathbf{x}})-\mathbf{S}^{-\frac{1}{2}}\cdot(\mathbf{x_1}-\bar{\mathbf{x}})\right]'\cdot \left[\mathbf{S}^{-\frac{1}{2}}\cdot(\mathbf{x_2}-\bar{\mathbf{x}})- \mathbf{S}^{-\frac{1}{2}}\cdot(\mathbf{x_1}-\bar{\mathbf{x}})\right]\\ % &=\left[ \mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_2}-\mathbf{S}^{-\frac{1}{2}}\cdot\bar{\mathbf{x}}- \mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_1}+\mathbf{S}^{-\frac{1}{2}}\cdot\bar{\mathbf{x}}\right]'\cdot \left[\mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_2}-\mathbf{S}^{-\frac{1}{2}}\cdot\bar{\mathbf{x}}- \mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_1}+\mathbf{S}^{-\frac{1}{2}}\cdot\bar{\mathbf{x}}\right]\\ % &=\left[\mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_2}-\mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_1}\right]' \cdot\left[\mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_2}-\mathbf{S}^{-\frac{1}{2}}\cdot\mathbf{x_1} \right]\\ &=\left[\mathbf{S}^{-\frac{1}{2}}(\mathbf{x_2}-\mathbf{x_1})\right]'\cdot \left[\mathbf{S}^{-\frac{1}{2}}(\mathbf{x_2}-\mathbf{x_1})\right]\\ &=(\mathbf{x_2}-\mathbf{x_1})'\cdot\left(\mathbf{S}^{-\frac{1}{2}}\right)'\cdot \mathbf{S}^{-\frac{1}{2}}\cdot(\mathbf{x_2}-\mathbf{x_1})\\ &=(\mathbf{x_2}-\mathbf{x_1})'\cdot \mathbf{S}^{-1}\cdot(\mathbf{x_2}-\mathbf{x_1}) \end{aligned} \] Zwei Punkte die den Abstand \(||\mathbf{z_2}-\mathbf{z_1}||^2\) haben, haben im ursprünglichen Koordinatensystem also die Beziehung \((\mathbf{x_2}-\mathbf{x_1})'\cdot \mathbf{S}^{-1}\cdot(\mathbf{x_2}-\mathbf{x_1})\). Diese Beziehung ist die Mahalanobis-Distanz.

11.3 Wahrscheinlichkeitstheoretische Grundlagen

In Statistik I haben wir bereits Zufallsvariablen kennengelernt, die einem Ergebnis eines Zufallsexperimentes \(\omega\in\Omega\) (oft eine Person oder eine Stichprobe, also einem Element von \(\Omega^n\)) eine Zahl aus \(\Omega'\) (oft \(\Omega'\subseteq\mathbb{R}\)) zuordnet. Im multivariaten Fall erweitern wir dieses Verständnis dahingehend, dass den Ergebnissen eines Zufallsexperimentes nicht nur ein Wert, sondern verschiedene Werte zugeordnet werden (die auch strukturell verschieden sein können, wie z.B. Reaktionszeiten und Fehler(raten)).

11.3.1 Multivariate Verteilungen

Unter einem Zufallsvektor versteht man einen Vektor \(\mathbf{X}\), dessen \(p\)-viele Elemente Zufallsvariablen \(\mathbf{X_1},\ldots,\mathbf{X_p}\) sind: \[ \mathbf{X}= \begin{pmatrix} \mathbf{X_1} \\ \vdots\\ \mathbf{X_p} \end{pmatrix} \] Der Übersichtlichkeit halber beschränken wir uns hier auf den Fall \(p=2\). Wir beginnen mit dem Fall, dass beide Elemente von \(\mathbf{X}\) diskrete Zufallsvariablen sind.

Nehmen wir an, wir hätten zwei Münzen geworfen. Dabei besitzt jede Münze ihre eigene Menge möglicher Ergebnisse: \[ \Omega_1 = \{K, Z\}\quad\text{und}\quad\Omega_2 = \{K, Z\} \] Nun bilden wir das Mengenprodukt beider Mengen, also \(\Omega_1\times\Omega_2\), welches sämtliche Kombinationen der beiden Mengen beschreibt. Das Ergebnis ist wiederum eine neue Menge, bei der jedes Element jeweils zwei Ereignisse (= ein 2-Tupel) enthält: \[ \Omega = \{(K,K) ; (K,Z) ; (Z,K) ; (Z,Z)\} \] Ein Zufallsvektor \(\mathbf{X}\) teilt dann jedem Element/Tupel \(\omega\in\Omega\) einen Vektor zu, wobei eine mögliche Zuordnungsvorschrift sein könnte:

• Das erste Element des zugeordneten Vektors enthält die Häufigkeit, mit der das Teilergebnis “Zahl” im Tupel \(\omega\) vorkommt.
• Das zweite Element des zugeordneten Vektors nimmt den Wert 1 an, wenn das Tupel \(\omega\) das Teilergebnis “Zahl” beinhaltet, sonst 0.

Dieser Zufallsvektor kann tabellarisch wie folgt dargestellt werden: Ganz ähnlich können für den diskreten Fall nun Wahrscheinlichkeitsfunktionen \(P(\mathbf{X}=\mathbf{x})\) angegeben werden. Bei Annahme eines Laplace-Experimentes tritt die Realisierung \((1,1)'\) mit \(p_{(1,1)'}=0.5\) auf, während sie für die anderen beiden Realisierungen \(0.25\) ist. Die Summe aller Wahrscheinlichkeiten ergibt, wie im univariaten Fall, wieder \(1\). Eine grafische Illustration ist in der folgenden Abbildung dargestellt.

Im stetigen Fall ist die Wahrscheinlichkeit für eine einzelne Realisierung wieder 0 und Wahrscheinlichkeiten müssen nun als Volumina angegeben werden. Auch hier gilt aber: Integriert man eine multivariate Dichtefunktion über alle Dimensionen gleichzeitig, so ergibt sich der Wert 1. Die wichtigste stetige Dichtefunktion im multivariaten Fall ist die multivariate Normalverteilung, die im nächsten Abschnitt ausführlicher behandelt wird.

11.3.2 Multivariate Normalverteilung

11.3.2.1 Kenngrößen und Schätzer

Ähnlich wie für die Beschreibung der Daten (siehe Abschnitt 1) gibt es auch auf der Populationsebene Kennwerte bzw. Parameter. Die konzeptuelle Erweiterung ist, dass wir im univariaten Fall Zufallsvariablen \[ \mathbf{X}:\Omega \mapsto \Omega' \] betrachtet haben und diese nun auf Zufallsvektoren verallgemeinern, also \[ \mathbf{X}= \begin{pmatrix} \mathbf{X_1} \\ \vdots\\ \mathbf{X_p} \end{pmatrix}: \Omega \mapsto \Omega'^p \] wobei die einzelnen \(\mathbf{X_i}\) Zufallsvariablen sind. Die folgende Tabelle stellt die Parameter der zentralen Tendenz einander gegenüber:

Ganz ähnlich ergeben sich für die Parameter der Datenvariabilität:

Als Schätzer eignen sich…

• …der Mittelwertvektor \(\mathbf{\bar{x}}\) für den Erwartungswertvektor \(\mathbf{\mu}\), sowie…
• …die korrigierte (Varianz-)Kovarianzmatrix \(\mathbf{\hat{S}}=\frac{n}{n-1}\mathbf{S}\) für die Kovarianzmatrix \(\Sigma\).

11.3.2.2 Definition und Dichte

Zur Erinnerung: Die univariate Normalverteilung ist durch ihre Varianz \(\sigma^2\) und ihren Erwartungswert \(\mu\) bestimmt und ihre Dichtefunktion lautet: \[ f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}} \]

Ganz ähnlich ist eine multivariate Normalverteilung durch ihren Erwartungswertvektor \(\mathbf{\mu}\) und ihre Kovarianzmatrix \(\Sigma\) eindeutig charakterisiert. Ist ein Zufallsvektor \(\mathbf{X}\) \(p\)-variat normalverteilt mit \(\mathbf{\mu}\) und \(\Sigma\), schreibt man auch: \[ \mathbf{X}\sim N_p(\mathbf{\mu},\Sigma) \] Ist \(\Sigma\) regulär, besitzt \(N_p(\mathbf{\mu},\Sigma)\) eine Wahrscheinlichkeitsdichte, die einer Realisierung, das heißt einem konkreten Vektor \(\mathbf{x}\), den Wert \[ f(\mathbf{x})=\frac{1}{\sqrt{(2\pi)^p \det\Sigma}}e^{{-\frac{1}{2}} (\mathbf{x}-\mathbf{\mathbf{\mu}})'\Sigma^{-1}(\mathbf{x}-\mathbf{\mathbf{\mu}})} \] zuweist. Man sieht, dass der wesentliche Bestandteil des Exponenten die Mahalanobis-Distanz von \(\mathbf{x}\) und \(\mathbf{\mu}\) ist.

11.3.2.3 Geometrie

Ein Beispiel für die Dichte einer Multinormalverteilung bei \(p=2\) ist in der folgenden Abbildung dargestellt. Die Mengen, auf denen die Dichtefunktion konstant ist (bei nur zwei betrachteten Dimensionen also die ‘’Höhenlinien’’) sind die Ränder des Ellipsoiden \(\varepsilon(\mathbf{\mu},{\mathbf{\Sigma}},k)\) für feste Werte von \(k\). Diese Konturen sind also die Menge aller Punkte, deren quadrierte Mahalanobis-Distanz zum Zentroiden konstant ist. Wahrscheinlichkeiten sind nun nicht mehr Flächen, sondern Volumina und das Gesamtvolumen bei \(p=2\) unter der bivariaten Normalverteilung ist 1. Wird sie gegen eine Fläche bestehend aus der \(y\)-Achse und der \(X_1\)- oder \(X_2\)-Achse projiziert, ergibt sich die jeweilige univariate Normalverteilung mit einer Fläche von 1 unter der Kurve.

Die Konturen (oder Niveaulinien) der bivariaten Normalverteilung können auch bestimmt werden für einen vorgegebenen Wert der Dichte oder eine gegebene Wahrscheinlichkeit. Ergänzende Informationen hierfür finden Sie hier.

11.4 Hotelling’s \(T^2\)-Test

Hotelling’s \(T^2\)-Tests sind Verallgemeinerungen des klassischen \(t\)-Tests mit einer abhängigen Variablen auf mehrere abhängige Variablen. Hierbei gibt es ganz ähnliche Situationen, wie die, denen wir in Teil 12, Statistik I von Statistik I bereits begegnet sind: Einstichprobenfall, zwei unabhängige Stichproben und zwei abhängige Stichproben.

11.4.1 Einstichprobenfall

Bei einer multivariaten Fragestellung sind prinzipiell zwei Vorgehensweisen denkbar:

11.4.1.1 p-viele univariate \(t\)-Tests

In diesem Fall gelten zunächst die gleichen Voraussetzungen wie beim bereits bekannten \(t\)-Test (Teil 12, Statistik I):

• Die jeweilige Variable \(X\) ist in der Population normalverteilt mit unbekanntem Erwartungswert \(\mu\) und unbekannter Varianz \(\sigma^2\): \(X\sim N(\mu,\sigma^2)\)
• Die einzelnen Messungen sind unabhängig voneinander.

Die zu testenden Hypothesen lauten (im ungerichteten Fall):

• \(H_0:\quad\mu=\mu_0\)
• \(H_1:\quad\mu\neq\mu_0\)

Aus den Daten wird als Prüfgröße \[t=\frac{\bar{X}-\mu_0}{\frac{\hat{S}}{\sqrt{n}}}\] berechnet, und für die entsprechende Zufallsvariable \(\mathbf{t}\) gilt bei Annahme der Gültigkeit der Nullhypothese \[\mathbf{t}\sim t_{n-1}.\]

Als Entscheidungsregel wurde damals festgehalten: “Verwirf \(H_0\), falls \(|t|\geq t_{n-1;1-\frac{\alpha}{2}}\)”.

11.4.1.2 Ein multivariater \(T^2\)-Test

Die Voraussetzungen des \(T^2\)-Tests sind Verallgemeinerungen des \(t\)-Tests:

• \(\mathbf{X}=\begin{pmatrix}X_1\\\ldots\\X_p\end{pmatrix}\) ist multivariat-normalverteilt mit unbekanntem Erwartungswertvektor und unbekannter, regulärer Kovarianzmatrix, also \(\mathbf{x}\sim N_p(\mathbf{\mu},\Sigma)\)
• Die Messungen sind unabhängig voneinander.

Entsprechende Generalisierungen sind auch die Hypothesen, die sich hier auf Erwartungswertvektoren beziehen:

• \(H_0:\quad\mathbf{\mu}=\mathbf{\mu_0}\), d.h. \(\begin{pmatrix}\mu_1\\\vdots\\\mu_p\end{pmatrix}=\begin{pmatrix} \mu_{01}\\\vdots\\\mu_{0p}\end{pmatrix}\)
• \(H_1:\quad\mathbf{\mu}\neq\mathbf{\mu_0}\), d.h. \(\begin{pmatrix}\mu_1\\\vdots\\\mu_p\end{pmatrix}\neq\begin{pmatrix} \mu_{01}\\\vdots\\\mu_{0p}\end{pmatrix}\)

Die Prüfgröße nennt sich naheliegenderweise \(T^2\) und berechnet sich als \[T^2=(n-1)\cdot(\mathbf{\bar{x}}-\mathbf{\mu_0})'\cdot \mathbf{S}^{-1}\cdot (\mathbf{\bar{x}}-\mathbf{\mu_0})\]

wobei \(\mathbf{\bar{x}}\) der Mittelwertvektor und \(\mathbf{S}\) die Kovarianzmatrix der Stichprobe sind.

Von einer entsprechenden Zufallsvariable kann nun wieder die theoretische Verteilung bestimmt werden und es gilt unter Annahme der Gültigkeit der Nullhypothese \[\frac{n-p}{p(n-1)}\cdot T^2\sim F_{p,n-p}\] Damit können wir als Entscheidungsregel formulieren: “Verwirf \(H_0\), falls \(\frac{n-p}{p(n-1)}\cdot T^2\geq F_{p,n-p;1-\alpha}\).”

Eine kleine Begründung der multivariaten \(T^2\)-Statistik kann ersichtlich werden, wenn vom univariaten \(t\)-Bruch ausgegangen wird: \[\begin{aligned} t=&\frac{\bar{X}-\mu_0}{\frac{\hat{S}}{\sqrt{n}}}\\ \Rightarrow t^2=&\frac{(\bar{X}-\mu_0)^2\cdot n}{ \hat{S^2}}\\ \Leftrightarrow t^2=&\frac{(\bar{X}-\mu_0)^2\cdot n}{\frac{n}{n-1}S^2}= (n-1)\cdot\frac{(\bar{X}-\mu_0)^2}{S^2}=(n-1)\cdot (\bar{x}-\mu_0) \cdot(S^2)^{-1}\cdot (\bar{X}-\mu_0) \end{aligned}\] Verallgemeinert man den letzten Ausdruck auf den multivariaten Fall, erhält man \[T^2=(n-1)\cdot(\mathbf{\bar{x}}-\mathbf{\mu_0})'\cdot \mathbf{S}^{-1}\cdot (\mathbf{\bar{x}}-\mathbf{\mu_0})\]

11.4.1.3 Ein Beispiel

Es soll geprüft werden, ob ein neuartiges Entspannungsverfahren eine physiologische Wirkung hat. Dazu werden die beiden Variablen “Herzfrequenz” (\(X_1\)) und “diastolischer Blutdruck” (\(X_2\)) nach der Entspannung erhoben und mit den Werten verglichen, die diese Variablen ohne Entspannung üblicherweise haben. Erhoben werden die in der folgenden Tabelle dargestellten Werte: Bekannt seien die Erwartungswerte für beide Variablen: \(\mu_{01}=60\) und \(\mu_{02}=80\). Es soll getestet werden, ob der Erwartungswertvektor der physiologischen Werte nach der Entspannung vom Vektor \(\mathbf{\mu_0}=(60,80)'\) abweicht, wobei wir \(\alpha=.10\) wählen (Wenn wir dann später zum Vergleich zwei univariate \(t\)-Tests berechnen, korrigieren wir dann für multiples Testen und verwenden dort ein Bonferroni-adjustiertes \(\alpha'=\frac{\alpha}{2}=.05\).):

• \(H_0:\quad\mathbf{\mu}=\mathbf{\mu_0}\quad\text{also}\quad \begin{pmatrix}\mu_1\\\mu_2\end{pmatrix}=\begin{pmatrix}60\\80 \end{pmatrix}\)
• \(H_1:\quad\mathbf{\mu}\neq\mathbf{\mu_0}\quad\text{also}\quad \begin{pmatrix}\mu_1\\\mu_2\end{pmatrix}\neq\begin{pmatrix}60\\80 \end{pmatrix}\)

X1 <- c(55, 59, 53, 55, 61, 59)         # Daten X1
X2 <- c(71, 82, 79, 76, 74, 74)         # Daten X2
n <- 6                                  # Anzahl Personen
p <- 2                                  # Anzahl Variablen
mu0 <- c(60,80)                         # Erwartungswertvektor
x <- c(mean(X1), mean(X2))              # Mittelwertvektor
Kov <- (n-1)/n*cov(X1, X2)              # Kovarianz
S2.X1 <- (n-1)/n*cov(X1, X1)            # Varianz X1
S2.X2 <- (n-1)/n*cov(X2, X2)            # Varianz X2
S <- rbind(c(S2.X1, Kov),               # Kovarianzmatrix S
           c(Kov, S2.X2))
round(S, 2)                             # Kovarianzmatrix ausgeben

##       [,1]  [,2]
## [1,]  8.00 -0.33
## [2,] -0.33 13.00

# dann kann daraus die Prüfgröße T^2 berechnet werden
T2 <- (n-1) * t(x-mu0) %*% solve(S) %*% (x-mu0)
T2

##          [,1]
## [1,] 12.17647

Wegen \[\frac{n-p}{p(n-1)}\cdot T^2\sim F_{p,n-p}\] wird noch folgende Transformation durchgeführt:

(n-p) / (p*(n-1)) * T2

##          [,1]
## [1,] 4.870588

Schließlich benötigen noch wir den kritischen Wert \(F_{\text{krit}}\), den wir wie gewohnt bestimmen können:

qf(p = 0.90, 
   df1 = 2,
   df2 = 4)

## [1] 4.324555

Da sich insgesamt also \[4.87\geq F_{\text{krit}}\] ergibt, entscheiden wir uns zugunsten der Alternativhypothese; der \(T^2\)-Test ist somit signifikant.

Eine entsprechende R-Funktion steht mit HotellingsT2Test() im Paket DescTools zur Verfügung. Im hier besprochenen Fall wird der Funktion eine Matrix übergeben, deren Spalten die Variablen \(X_1\) und \(X_2\) sind, sowie der Erwartungsvektor, auf den getestet werden soll:

library(DescTools)
X <- cbind(X1, X2)         # X1 und X2 als Spalten einer Matrix...
HotellingsT2Test(x = X,    # ...die hier übergeben wird
                 mu = mu0)

## 
##  Hotelling's one sample T2-test
## 
## data:  X
## T.2 = 4.8706, df1 = 2, df2 = 4, p-value = 0.08474
## alternative hypothesis: true location is not equal to c(60,80)

Auch mit dem \(p\)-Wert folgt also – für \(\alpha=.10\) – dass der Test signifikant ist.

Zum Vergleich schauen wir uns nun noch die beiden univariaten \(t\)-Tests mit \(\alpha'=.05\):

t_out(t.test(X1, 
             mu = mu0[1])) # nicht signifikant

##                 Test                           Results
## 1 One Sample t-test: t(5) = -2.37, p = .064, d = -0.97
## 
## NOTE: Reporting unadjusted estimate for Cohen's d.

t_out(t.test(X2,
             mu = mu0[2])) # nicht signifikant

##                 Test                           Results
## 1 One Sample t-test: t(5) = -2.48, p = .056, d = -1.01
## 
## NOTE: Reporting unadjusted estimate for Cohen's d.

Der \(T^2\)-Test ist also insgesamt offenbar sensibler (d.h. er verfügt über eine größere Power), zumindest dann, wenn mit für multiples Testen korrigierten \(\alpha\)-Werten gearbeitet wird.

11.4.2 Multivariate Konfidenzbereiche

Zunächst berechnen wir – wie in Teil 13 von Statistik I – univariate Vertrauensbereiche (d.h. Konfidenzintervalle).

Dass der \(t\)-Test für eine der Variablen \(X_1\) oder \(X_2\) signifikant wird, bedeutet allgemein, dass \(|t|\geq t_{n-1;1-\frac{\alpha'}{2}}\) sein muss. Mit anderen Worten: \[ \begin{aligned} \frac{M_X-\mu_0}{\frac{\hat{S}}{\sqrt{n}}}\geq t_{n-1;1-\frac{\alpha'}{2}} \quad&\text{oder}\quad \frac{M_X-\mu_0}{\frac{\hat{S}}{\sqrt{n}}}\leq -t_{n-1;1-\frac{\alpha'}{2}}\\ \Leftrightarrow \mu_0\leq M_X-\frac{\hat{S}}{\sqrt{n}}t_{n-1;1-\frac{\alpha'}{2}} \quad&\text{oder}\quad \mu_0\geq M_X+\frac{\hat{S}}{\sqrt{n}}t_{n-1;1-\frac{\alpha'}{2}} \end{aligned} \] Anders ausgedrückt wird der \(t\)-Test nicht signifikant, wenn der Testwert \(\mu_0\) innerhalb des \((1-\alpha')\cdot 100\%\)-Konfidenzintervalls um \(M_X\) liegt, also wenn \[\mu_0\in \left[M_X-\frac{\hat{S}}{\sqrt{n}}t_{n-1;1-\frac{\alpha'}{2}}; M_X+\frac{\hat{S}}{\sqrt{n}}t_{n-1;1-\frac{\alpha'}{2}}\right]\]

Für die Beispieldaten können wir die Konfidenzintervalle leicht mit der Funktion t.test() ausgegeben lassen:

t.test(X1, mu = mu0[1])$conf.int

## [1] 53.74844 60.25156
## attr(,"conf.level")
## [1] 0.95

t.test(X2, mu = mu0[2])$conf.int

## [1] 71.85506 80.14494
## attr(,"conf.level")
## [1] 0.95

Daraus ergeben sich also die folgenden beiden Beziehungen:

• \(t\)-Test für \(X_1\) wird nicht signifikant \(\Leftrightarrow \mu_{01}\in [53.75;60.25]\)
• \(t\)-Test für \(X_2\) wird nicht signifikant \(\Leftrightarrow \mu_{02}\in [71.85;80.15]\)

Die nachfolgende Abbildung illustriert die Bereiche von \(\mu_0\) in denen keiner, einer oder beide der \(t\)-Tests signifikant werden würden.

Analog zu diesem Vorgehen, bilden wir nun einen bivariaten Vertrauensbereich. Dass der \(T^2\)-Test signifikant wird, bedeutet, dass \[ \begin{aligned} &\frac{n-p}{p(n-1)}\cdot T^2\geq F_{p,n-p;1-\alpha}\\ \Leftrightarrow&\frac{n-p}{p(n-1)}\cdot (n-1)\cdot(\mathbf{\bar{x}}-\mathbf{\mu_0})'\cdot \mathbf{S}^{-1}\cdot (\mathbf{\bar{x}}-\mathbf{\mu_0})\geq F_{p,n-p;1-\alpha}\\ \Leftrightarrow& (\mathbf{\bar{x}}-\mathbf{\mu_0})'\cdot \mathbf{S}^{-1}\cdot (\mathbf{\bar{x}}-\mathbf{\mu_0})\geq \frac{p}{n-p}\cdot F_{p,n-p;1-\alpha}\\ \end{aligned} \] Anders ausgedrückt wird der \(T^2\)-Test nicht signifikant, wenn \(\mathbf{\mu_0}\) innerhalb eines Ellipsoiden mit Mittelpunkt \(\mathbf{\bar{x}}\) und Radius \(\frac{p}{n-p}\cdot F_{p,n-p;1-\alpha}\) liegt, oder kurz wenn \[ \mathbf{\mu_0}\in\varepsilon \left(\mathbf{\bar{x}},\mathbf{S}, \frac{p}{n-p}\cdot F_{p,n-p;1-\alpha}\right) \]

Der bivariate Vertrauensbereich ist also ein Ellipsoid und der rote Bereich in der folgenden Abbildung ist entsprechend der Bereich, in dem der \(T^2\)-Test für \(H_0: \mathbf{\mu}=\mathbf{\mu_0}\) nicht signifikant wird (die gestrichelten Linien sind die univariaten Konfidenzintervalle).

Man sieht, dass bei einem normalen \(t\)-Test sich ein Kreuz zweier Konfidenzintervalle ergibt, während bei dem multivariaten Vertrauensbereich sich aufgrund der Kovarianzmatrix eine Ellipse ergibt. In Richtung der kürzeren Hauptachse ist ein \(T^2\)-Test einzelnen \(t\)-Tests überlegen. Insbesondere wird der Unterschied größer, je höher die Variablen miteinander korreliert sind: Je größer die Korrelation, desto wahrscheinlicher ist ein Vorteil für den multivariaten Test.

Zusammenfassend also: \[ \varepsilon\left(\mathbf{\bar{x}},\mathbf{S}, \frac{p}{n-p}\cdot F_{p,n-p;1-\alpha}\right)\text{ ist der $(1-\alpha)\cdot 100\%$ ''bivariate'' Vertrauensbereich für } \mathbf{\mu} \]

Wegen der eingehenden Kovarianzstruktur spielt hierbei also nicht mehr nur die Lage, sondern auch die Form des Vertrauensbereiches eine Rolle.

11.4.3 Zwei unabhängige Stichproben

Der hier betrachtete Fall ist nun die Generalisierung des \(t\)-Tests für zwei unabhängige Stichproben auf mehr als eine abhängige Variable.

Seien \[\mathbf{X_1}= \begin{pmatrix} X_{11}\\ \vdots\\ X_{1p} \end{pmatrix}\quad\text{und}\quad \mathbf{X_2}= \begin{pmatrix} X_{21}\\ \vdots\\X_{2p} \end{pmatrix}\]

zwei \(p\)-dimensionale Zufallvektoren, wobei \(\mathbf{X_1}\sim N_p(\mathbf{\mu_1},\Sigma_1)\) und \(\mathbf{X_2}\sim N_p(\mathbf{\mu_2},\Sigma_2)\) verteilt sind. Für die Kovarianzmatrizen gilt weiter, dass \(\Sigma_1,\) und \(\Sigma_2\) regulär und homogen sind, also \(\Sigma_1=\Sigma_2\).

Die Hypothesen lauten:

• \(H_0:\quad\mathbf{\mu_1}=\mathbf{\mu_2}\text{, d.h. } \begin{pmatrix}\mu_{11}\\\vdots\\\mu_{1p}\end{pmatrix}= \begin{pmatrix}\mu_{21}\\\vdots\\\mu_{2p}\end{pmatrix}\)
• \(H_1:\quad\mathbf{\mu_1}\neq\mathbf{\mu_2}\text{, d.h. } \begin{pmatrix}\mu_{11}\\\vdots\\\mu_{1p}\end{pmatrix}\neq \begin{pmatrix}\mu_{21}\\\vdots\\\mu_{2p}\end{pmatrix}\)

Mit zwei unabhängig voneinander gezogenen Stichproben der Umfänge \(n_1\) bzw. \(n_2\) ergeben sich als Mittelwertsvektoren dieser Stichproben \(\mathbf{\bar{x}_1}\) und \(\mathbf{\bar{x}_2}\) sowie \(\mathbf{S_1}\) und \(\mathbf{S_2}\) als Kovarianzmatrizen.

Die Teststatistik lautet dann \[T^2=\frac{n_1n_2(n_1+n_2-2)}{n_1+n_2}\cdot (\mathbf{\bar{x}_1-\bar{x}_2})'\cdot (n_1\mathbf{S_1}+n_2\mathbf{S_2})^{-1}\cdot (\mathbf{\bar{x}_1 -\bar{x}_2}) \] woraus sich als Entscheidungsregel ergibt: “Verwirf \(H_0\), falls \(\frac{n_1+n_2-p-1}{(n_1+n_2-2)p}\cdot T^2\geq F_{p,n_1+n_2-p-1;1-\alpha}\).”

11.4.4 Zwei abhängige Stichproben

Wir betrachten nun noch den Fall der Generalisierung des \(t\)-Tests für zwei abhängige Stichproben auf mehr als eine abhängige Variable. Ganz ähnlich wie im univariaten Fall ist das Vorgehen, die Situation auf einen Einstichproben-\(T^2\)-Test der Differenzvariablen zurückzuführen.

Seien \[\mathbf{X_1}= \begin{pmatrix} X_{11}\\ \vdots\\ X_{1p} \end{pmatrix}\quad\text{und}\quad \mathbf{X_2}= \begin{pmatrix} X_{21}\\ \vdots\\X_{2p} \end{pmatrix}\]

zwei \(p\)-dimensionale Zufallvektoren, die die Erhebung von \(p\)-vielen Werten an gepaarten (normalerweise also identischen) Versuchspersonen zu zwei verschiedenen Zeitpunkten beschreiben. Der Zufallsvektor \((\mathbf{X_1},\mathbf{X_2})'\) sei in der Population \(p\)-variat normalverteilt mit regulärer Kovarianzmatrix \(\Sigma\).

Die Hypothesen lauten:

• \(H_0:\quad\mathbf{\mu_1}=\mathbf{\mu_2}\text{, d.h. }\begin{pmatrix}\mu_{11}\\\vdots\\\mu_{1p}\end{pmatrix}=\begin{pmatrix}\mu_{21}\\\vdots\\\mu_{2p}\end{pmatrix}\)
• \(H_1:\quad\mathbf{\mu_1}\neq\mathbf{\mu_2}\text{, d.h. }\begin{pmatrix}\mu_{11}\\\vdots\\\mu_{1p}\end{pmatrix}\neq\begin{pmatrix}\mu_{21}\\\vdots\\\mu_{2p}\end{pmatrix}\)

Die zwei Stichproben vom Umfang \(n\) ergeben insgesamt zwei \((n\times p)\)-Datenmatrizen \[ \mathbf{X_1} = \begin{pmatrix} x_{111} & x_{112} & \cdots & x_{11p} \\ x_{121} & x_{122} & \cdots & x_{12p} \\ \vdots & \vdots & \vdots & \vdots \\ x_{1n1} & x_{1n2} & \cdots & x_{1np} \\ \end{pmatrix} \quad\text{und}\quad \mathbf{X_2} = \begin{pmatrix} x_{211} & x_{212} & \cdots & x_{21p} \\ x_{221} & x_{222} & \cdots & x_{22p} \\ \vdots & \vdots & \vdots & \vdots \\ x_{2n1} & x_{2n2} & \cdots & x_{2np} \\ \end{pmatrix} \]

aus denen die Differenzmatrix \[ \mathbf{D} = \mathbf{X_1}-\mathbf{X_2} = \begin{pmatrix} x_{111}-x_{211} & x_{112}-x_{212} & \cdots & x_{11p}-x_{21p} \\ x_{121}-x_{221} & x_{122}-x_{222} & \cdots & x_{12p}-x_{22p} \\ \vdots & \vdots & \vdots & \vdots \\ x_{1n1}-x_{2n1} & x_{1n2}-x_{2n2} & \cdots & x_{1np}-x_{2np} \\ \end{pmatrix} \]

berechnet wird. Im Weiteren wird dann der Mittelwertvektor von \(\mathbf{D}\), \(\mathbf{\bar{d}}\), sowie die Kovarianzmatrix von \(\mathbf{D}\), \(\mathbf{S_D}\), berechnet. Getestet wird dann die \(H_0:\mathbf{\mu_D}=\mathbf{0}\) wobei \(\mathbf{0}\) der Vektor ist, dessen \(p\)-viele Elemente 0 sind.

Die Teststatistik lautet dann \[T^2=(n-1)\cdot \mathbf{\bar{d}}'\cdot \mathbf{S_D^{-1}}\cdot \mathbf{\bar{d}}\] woraus sich als Entscheidungsregel ergibt: “Verwirf \(H_0\), falls \(\frac{n-p}{p(n-1)}\cdot T^2\geq F_{p,n-p;1-\alpha}\).”

11.5 Ausblick

Wir haben hier nur einen kleinen Einblick in multivariate Verfahren geben können. Natürlich gibt es weitere Verfahren, die wir auch bereits kennen, und die auch auf den Fall mehrerer abhängiger Variablen verallgemeinert werden können. Insbesondere trifft dies auf die Varianzanalyse zu.

11.5.1 Multivariate Varianzanalyse

Im Rahmen der Varianzanalyse (ANOVA) haben wir verschiedene Fälle betrachtet, denen aber gemeinsam war, dass es nur eine abhängige Variable gab. Es gibt aber auch die Multivariate Varianzanalyse (MANOVA) mit mehreren abhängigen Variablen; gewissermaßen die Verallgemeinerung des \(T^2\)-Tests auf z.B. mehr als zwei Gruppen, Bedingungen, usw.

Ein Ansatz zur Durchführung ist, ähnlich wie in der univariaten ANOVA, die Zerlegung der Beobachtungen (die hier ja Vektoren sind) in mehrere Bestandteile: Im einfaktoriellen Fall in einen Grand Mean Vektor, einen Effektvektor und einen Residuenvektor. Aus diesen Abweichungen werden Kreuzprodukte gebildet und über Gruppen und/oder Personen aufsummiert. Auch ähnlich zum univariaten Fall resultieren dann eine Kreuzproduktsumme within \(\mathbf{W}\), eine Kreuzproduktsumme between \(\mathbf{B}\) und eine eine Kreuzproduktsumme total \(\mathbf{T}\) mit \[\mathbf{W}+\mathbf{B}=\mathbf{T}\] Als Prüfgröße ergibt sich Wilk’s Lambda \[ \Lambda=\frac{\det\mathbf{W}}{\det\mathbf{T}} \] und die entsprechende Zufallsvariable ist bei Annahme der Nullhypothese Wilk’s Lambda-verteilt, also \[\Lambda\sim \Lambda_{N-L,L-1}\] wobei \(N\) die Gesamtzahl der Versuchspersonen ist und \(L\) die Anzahl der Gruppen.(\(\Lambda\) entspricht zudem dem Produkt der Eigenwerte von \(\mathbf{WT}^{-1}\).) Es gilt weiter \(0\leq \Lambda\leq 1\); allerdings sprechen hier kleinere Werte gegen die \(H_0\) und für die \(H_1\).

Für größere Werte von \(N\) (z.B. \(N>8\)) gilt außerdem, dass das Bartlett-transformierte \(\Lambda\) approximativ \(\chi^2\)-verteilt ist, also \[ \left( N-1-\frac{L+p}{2}\right)\cdot(-\ln\Lambda)\sim\chi^2_{p(L-1)} \] wobei \(p\) die Anzahl der Variablen ist.

Während \(\Lambda\) als Likelihood-Ratio-Test konstruiert wurde, gibt es auch andere Möglichkeiten, Prüfgrößen zu konstruieren. Im Fall der MANOVA führt dies tatsächlich zu anderen Prüfgrößen (ohne hier auf die Verteilungsformen einzugehen):

• Roy’s größter Eigenwert oder Roy’s Maximalwurzel ist der größte Eigenwert \(\lambda_1\) von \(\mathbf{W}^{-1}\mathbf{B}\).
• Pillai’s Spurkriterium ist die Spur der Matrix \(\mathbf{T}^{-1}\mathbf{B}\).
• Die Hotelling-Lawley Spur ist die Spur der Matrix \(\mathbf{W}^{-1}\mathbf{B}\).

Im Detail haben die Prüfgrößen ihre (leichten) Vor- und Nachteile und es kann passieren, dass ein Test bei einer Prüfgröße signifikant wird, bei einer anderen nicht (siehe Olson, 1976). Grundlage für alle Prüfgrößen sind aber die Matrizen \(\mathbf{W}\) und \(\mathbf{B}\) bzw. genauer die Eigenwerte von \(\mathbf{W}^{-1}\) und \(\mathbf{B}\).

Wenn \(p=1\) (also univariater Fall) oder \(J=2\) (also zwei Gruppen) ist, dann sind alle vier Prüfgrößen äquivalent und im Fall \(J=2\) sind alle vier auch äquivalent zu Hotelling’s \(T^2\).

11.5.2 Multivariater Ansatz zur Messwiederholungs-Varianzanalyse

In Teil 5 von Statistik II hatten wir die Messwiederholungs-ANOVA eingeführt, die einen univariaten Zugang darstellt. Jede Variable (d.h. die Variablen zu den Messzeitpunkten) wird für sich als normalverteilt aufgefasst, die Zusammenhänge zwischen den Variablen werden nicht beschrieben und Spärizität war eine wichtige Voraussetzung.

Es gibt für die Problemstellung der Messwiederholungs-ANOVA zudem einen multivariaten Ansatz, bei dem die \(p\)-vielen Messzeitpunkte als \(p\)-viele Variablen und in der Folge als \(p\)-variat normalverteilte (\(p\)-dimensionale) Variable aufgefasst werden. Die Idee ist dann, eine Kontrastmatrix zu definieren, die alle paarweisen Differenzen von Erwartungswerten testet und die Nullhypothese lautet, dass sie alle Null sind. Diese Vorgehensweise benötigt die Annahme der Sphärizität nicht, hat aber eine geringere Power als die univariate Herangehensweise.

11.6 Literatur

Olson, C.L. (1976). On choosing a test statistic in multivariate analysis of variance. Psychological Bulletin, 83, 579-586.

12 Einführung in das Allgemeine Lineare Modell (ALM)

Wir haben bei der multiplen linearen Regression festgehalten, dass die Lösung des Minimierungsproblems nur bei orthogonalen Prädiktoren problemlos funktioniert, auch wenn Statistikprogramme wie R keinerlei Probleme haben und wir auch mit der Funktion lm() Regressionen verschiedener Art bereits berechnet haben. Dabei haben wir auch schon auf das sog. Allgemeine Lineare Modell (ALM) hingewiesen, welches den formalen Überbau über Regressionen und Varianzanalysen bietet. Dass beide Verfahren zu gleichen Ergebnissen führen und Varianzanalysen auch wie Regressionen behandelt werden können, haben wir auch bereits angemerkt. In diesem Teil wollen wir daher das ALM in seinen Grundzügen kennenlernen.

Mit dem Vorwissen aus Teil 11 zu Grundlagen der Vektoren- und Matrixrechnung, können wir uns dem ALM leichter zuwenden, da wird relativ viel Gebrauch von den dort eingeführten Notationen machen werden. In den Teilen 7 und 8 haben wir die multiple Regression ausführlich kennengelernt. Nun gehen wir den Schritt zur Notation im ALM und leiten daran – etwas oberflächlich – eine wichtige Beziehung her. Im Anschluss werden wir mit R ausprobieren, ob diese Rechnung tatsächlich richtig ist und formulieren danach exemplarisch noch den \(t\)-Test für zwei unabhängige Stichproben im ALM.

12.1 Von der Regression zum ALM

Wir beginnen mit dem Fall dreier Versuchspersonen und zwei Prädiktoren und generalisieren anschließend diesen Fall. Gibt es nur drei Versuchspersonen und zwei Prädiktoren, kann die Regressiongleichung wie folgt für jede Versuchsperson (aus)geschrieben werden: \[\begin{equation*} \begin{aligned} y_1 &= a + b_1\cdot x_{11} + b_2\cdot x_{12} + \epsilon_1\\ y_2 &= a + b_1\cdot x_{21} + b_2\cdot x_{22} + \epsilon_2\\ y_3 &= a + b_1\cdot x_{31} + b_2\cdot x_{32} + \epsilon_3\\ \end{aligned} \end{equation*}\] Hier können wir nun unsere Kenntnisse der Matrizenrechnungen direkt einbringen und diese Gleichungen in Form mehrerer Matrizen bzw. Vektoren schreiben. Dazu setzen wir folgende Matrizen fest: \[\begin{align*} \bf y = \begin{pmatrix} y_1\\ y_2 \\ y_3 \end{pmatrix} \qquad \bf X = \begin{pmatrix} 1 & x_{11} & x_{12}\\ 1 & x_{21} & x_{22}\\ 1 & x_{31} & x_{32} \end{pmatrix} \qquad \bf b = \begin{pmatrix} a\\ b_1 \\ b_2 \end{pmatrix} \qquad \mathbf{\epsilon} = \begin{pmatrix} \epsilon_1\\ \epsilon_2 \\ \epsilon_3 \end{pmatrix} \end{align*}\] Wir nennen dabei \(\bf b\) den Parametervektor und \(\bf X\) die Designmatrix. Mithilfe dieser Matrizen bzw. Vektoren können die Gleichungen der drei Versuchspersonen auch geschrieben werden als: \[\bf y=\bf X\cdot \bf b + \mathbf{\epsilon}\] Nun nehmen wir an, es gäbe \(n\) Versuchspersonen und \(q\) Prädiktoren. Dann schreiben wir für jede Beobachtung \(i\) die Regressionsgleichung auch zunächst einmal explizit aus:

\[\begin{equation*} \begin{aligned} y_1 &= a + b_1\cdot x_{11} + b_2\cdot x_{12} + \ldots + b_q\cdot x_{1q} + \epsilon_1\\ y_2 &= a + b_1\cdot x_{21} + b_2\cdot x_{22} + \ldots + b_q\cdot x_{2q} + \epsilon_2\\ \vdots &\\ y_i &= a + b_1\cdot x_{i1} + b_2\cdot x_{i2} + \ldots + b_q\cdot x_{iq} + \epsilon_i\\ \vdots &\\ y_n &= a + b_1\cdot x_{n1} + b_2\cdot x_{n2} + \ldots + b_q\cdot x_{nq} + \epsilon_n\\ \end{aligned} \end{equation*}\]

Ganz ähnlich wie gerade, können wir diese Gleichungen in Form von Matrizen schreiben:

\[\begin{align*} \bf y = \begin{pmatrix} y_1\\ y_2 \\ \vdots\\ y_n \end{pmatrix} \qquad \bf X = \begin{pmatrix} 1 & x_{11} & x_{12} & \dots & x_{1q}\\ 1 & x_{21} & x_{22} & \dots & x_{2q}\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} & \dots & x_{nq} \end{pmatrix} \qquad \bf b = \begin{pmatrix} a\\ b_1 \\ b_2 \\ \vdots\\ b_q \end{pmatrix} \qquad \mathbf{\epsilon} = \begin{pmatrix} \epsilon_1\\ \epsilon_2 \\ \vdots\\ \epsilon_n \end{pmatrix} \end{align*}\]

Auch hier können wir die Gleichungen mit diesen Matrizen kurz schreiben als: \[\bf y=\bf X\cdot \bf b + \mathbf{\epsilon}\] Der schwierige Teil ist nun: Wie bestimmen bzw. schätzen wir \(\bf b\), damit diese Gleichung bei einem Datensatz einen möglichst geringen Fehler ergibt? Dies kann geometrisch begründet werden. Im Prinzip sagen wir ja wieder Werte vorher und wir nennen den Vektor der vorhergesagten Werte \(\mathbf{\hat{y}}\) und den Vektor der die Gleichung dann erfüllt nennen wir \(\mathbf{\hat{b}}\), den Schätzer des Parametervektors. Es gilt dann: \[\mathbf{\hat{y}}=\bf X\cdot \mathbf{\hat{b}}\] Wir müssen nun einen Weg finden, \(\mathbf{\hat{b}}\) durch bekannte Größen schreiben zu können, also durch die Matrix \(\bf X\) und den Vektor \(\bf y\). Der Vektor \(\mathbf{\hat{y}}\) ergibt sich nach der obigen Gleichung durch die Linearkombination unserer Prädiktoren (die in \(\bf X\) stecken) mit dem Schätzer des Parametervektors (\(\mathbf{\hat{b}}\)). Hieraus folgt, in den Worten der Linearen Algebra, dass sich \(\mathbf{\hat{y}}\) nur in demjenigen Unterraum bewegen kann, der durch die Prädiktoren \(\bf x_j\) der Matrix \(\bf X\) aufgespannt wird. Haben wir beispielsweise zwei Prädiktoren \(\bf x_1\) und \(\bf x_2\), so muss \(\mathbf{\hat{y}}\) innerhalb der Ebene liegen, die durch \(\bf x_1\) und \(\bf x_2\) aufgespannt wird:

Da wir wieder eine Minimierung der Abweichung der empirischen von den vorhergesagten Werten anstreben, suchen wir nun die kürzeste Distanz zwischen \(\bf y\) und \(\mathbf{\hat{y}}\). Dies ist genau dann der Fall, wenn \(\mathbf{\hat{y}}\) die orthogonale (“senkrechte”) Projektion von \(\bf y\) auf den von den Prädiktoren aufgespannten Unterraum \(L(\mathbf{x_1},\mathbf{x_2},\ldots,\mathbf{x_q})\) ist. Die folgende Abbildung illustriert diese Situation für das Beispiel mit zwei Prädiktoren. Die gestrichelte Linie stellt hierbei die Abweichung der empirischen von den vorhergesagten Werten ab, die minimiert werden soll (je kürzer diese Linie ist, desto besser ist unsere Vorhersage):

Wichtig ist nun die Einsicht, dass die orthogonale Projektion von \(\bf y\) auf den Unterraum (also \(\mathbf{y}-\mathbf{\hat{y}}\)) auch orthogonal zu jedem einzelnen Vektor \(\bf x_j\) steht. Da das Skalarprodukt orthogonaler Vektoren aber immer Null ergibt, folgt daraus:

\[\begin{equation*}\begin{split} <\mathbf{y}-\mathbf{\hat{y}},\mathbf{x_1}>=0&\quad\Leftrightarrow \quad<\mathbf{y},\mathbf{x_1}>-<\mathbf{\hat{y}}, \mathbf{x_1}>=0\\ <\mathbf{y}-\mathbf{\hat{y}},\mathbf{x_2}>=0&\quad\Leftrightarrow \quad<\mathbf{y},\mathbf{x_2}>-<\mathbf{\hat{y}}, \mathbf{x_2}>=0\\ \vdots&\quad\Leftrightarrow\quad\vdots\\ <\mathbf{y}-\mathbf{\hat{y}},\mathbf{x_q}>=0&\quad\Leftrightarrow \quad<\mathbf{y},\mathbf{x_q}>-<\mathbf{\hat{y}}, \mathbf{x_q}>=0\\ \end{split} \end{equation*}\] Unter Ausnutzung der Matrizenschreibweise kann dieses System von \(q\)-vielen Gleichungen auch wie folgt geschrieben werden: \[ \mathbf{X}'\cdot (\mathbf{y}-\mathbf{\hat{y}})=\mathbf{0} \Leftrightarrow \mathbf{X}'\cdot\\\mathbf{y}-\mathbf{X}'\cdot \mathbf{\hat{y}}=\mathbf{0} \] Wegen \(\mathbf{\hat{y}}=\mathbf{X}\cdot\mathbf{\hat{b}}\), folgt daraus \[\begin{equation*}\begin{split} &\mathbf{X}'\cdot\mathbf{y}-\mathbf{X}'\cdot\mathbf{X}\cdot \mathbf{\hat{b}}=\mathbf{0}\\ \Leftrightarrow&\mathbf{X}'\cdot\mathbf{y}=(\mathbf{X}' \cdot\mathbf{X})\cdot\mathbf{\hat{b}} \end{split} \end{equation*}\]

Das Produkt \(\mathbf{X}'\cdot\mathbf{X}\) ist nun eine symmetrische Matrix, die, voller Rang vorausgesetzt, zudem invertierbar ist. Unter diesen Umständen existiert dann also auch \[(\mathbf{X}'\cdot \mathbf{X})^{-1}\] und Multiplikation von links ergibt:
\[ \boxed{(\mathbf{X}'\cdot\mathbf{X})^{-1}\cdot\mathbf{X}'\cdot \mathbf{y}=\mathbf{\hat{b}}} \]

Die letzte Beziehung ist die vielleicht wichtigste Ableitung aus der vorangegangenen Betrachtung. Während wir in Teil 7 und 8 zwar Multiple Regressionen mit R berechnet haben, aber kein Werkzeug hatten, um die Koeffizienten quasi ohne R bzw. die Funktion lm() zu bestimmen, haben wir nun die nötige Gleichung dafür.

12.2 Anwendung in R

Um die Richtigkeit der Gleichung zu demonstrieren, benutzen wir wieder einen Datensatz, den wir in den vorangegangenen Teilen bereits verwendet haben. Wir laden den Datensatz, berechnen die Regression und geben uns die Koeffizienten aus:

daten <- read.table("./Daten/Album_Verkaeufe.dat",
                    header = TRUE)
modell <- lm(Verkaeufe ~ Werbung + Frequenz + Attraktivitaet, 
             data = daten)
coef(modell)

##    (Intercept)        Werbung       Frequenz Attraktivitaet 
##   -26.61295836     0.08488483     3.36742517    11.08633520

Nun bestimmen wir die Matrizen, die wir benötigen, nämlich \(\bf y\) und \(\bf X\). Der Vektor \(\bf y\) ist sehr einfach zu bestimmen, da dieser ja nur die Werte der Personen auf dem Kriterium beinhaltet, also auf der Variablen Verkaeufe. Wir nehmen daher einfach die entsprechende Spalte des DataFrames daten und wandeln diese direkt mit as.matrix in eine Matrix um:

y <- as.matrix(daten$Verkaeufe)

Die Designmatrix \(\bf X\) hat als erste Spalte 200-mal eine 1 und danach folgen die Spalten 1, 3 und 4 aus dem DataFrame daten, also die Prädiktoren. Diese Matrix können wir wie folgt zusammenstellen:

einsen <- rep(1, 200)                       # 200-mal eine 1
praediktoren <- as.matrix(daten[,c(1,3,4)]) # die drei Prädiktoren
X <- cbind(einsen, praediktoren)
head(X)                                     # die ersten 6 Zeilen ausgeben

##      einsen  Werbung Frequenz Attraktivitaet
## [1,]      1   10.256       43             10
## [2,]      1  985.685       28              7
## [3,]      1 1445.563       35              7
## [4,]      1 1188.193       33              7
## [5,]      1  574.513       44              5
## [6,]      1  568.954       19              5

Nun haben wir alle Bestandteile, um \(\mathbf{\hat{b}}\) gemäß der Formel \[ (\mathbf{X}'\cdot\mathbf{X})^{-1}\cdot\mathbf{X}'\cdot \mathbf{y}=\mathbf{\hat{b}} \] zu bestimmen. Der Übersichtlichkeit halber machen wir das hier schrittweise:

teil_1 <- t(X) %*% X       # X' * X
teil_2 <- solve(teil_1)    # (X' * X)^-1 <-- die Inverse bestimmen...
teil_3 <- teil_2 %*% t(X)  # ...das dann mit X' multiplizieren...
b.dach <- teil_3 %*% y     # ...und das alles noch mit y multiplizieren 
b.dach                     # Tadaaaaa.... genau die Koeffizienten wie mit lm() bestimmt

##                        [,1]
## einsen         -26.61295836
## Werbung          0.08488483
## Frequenz         3.36742517
## Attraktivitaet  11.08633520

coef(modell)               # zum Vergleich

##    (Intercept)        Werbung       Frequenz Attraktivitaet 
##   -26.61295836     0.08488483     3.36742517    11.08633520

Dass R intern das übrigens genauso löst – lm() steht ja auch für linear model – sieht man auch daran, wenn man mit der Funktion model.matrix() die Designmatrix eines berechneten lm()-Objektes ausgibt. Wir tun das hier für die ersten sechs Zeilen und sehen, dass dies genau \(\bf X\) ist:

designmatrix <- model.matrix(modell)
head(designmatrix)

##   (Intercept)  Werbung Frequenz Attraktivitaet
## 1           1   10.256       43             10
## 2           1  985.685       28              7
## 3           1 1445.563       35              7
## 4           1 1188.193       33              7
## 5           1  574.513       44              5
## 6           1  568.954       19              5

Lediglich ist die Spalte einsen überschrieben mit (Intercept), da durch diese Spalte der Achsenabschnitt a in das Modell einfließt.

12.3 Ein weiteres Beispiel: der \(t\)-Test für zwei unabhängige Stichproben

Wir können nun – quasi per Analogie – die bisher verwendeten Formeln noch etwas verallgemeinern. Insbesondere den Parametervektor verändern wir leicht und nennen ihn fortan \(\mathbf{\beta}\). Die Parameter in der multiplen Regression waren der Achsenabschnitt \(a\) und die Koeffizienten \(b_1,\ldots,b_q\). Im allgemeinen Fall nennen wir die Parameter nun einfach \(\beta_1,\ldots,\beta_p\): \[ \mathbf{\beta}= \begin{pmatrix} \beta_1\\ \beta_2\\ \vdots\\ \beta_p \end{pmatrix} \] Mit \(\mathbf{\hat{\beta}}\) meinen wir nun den Vektor, der \(\mathbf{\beta}\) schätzt. Da wir nur einen anderen Namen benutzen, ändert sich nichts an der Schätzung und es gilt immer noch: \[ \boxed{\mathbf{\hat{\beta}}=(\mathbf{X}'\cdot\mathbf{X})^{-1}\cdot\mathbf{X}'\cdot \mathbf{y}} \] Die “Kunst” ist nun sich zwei Dinge zu überlegen:

1. Was sind die Parameter für den vorliegenden Spezialfall (also \(t\)-Test, Varianzanalyse, einfache lineare Regression, multiple lineare Regression, …)?
2. Wie sieht die Designmatrix \(\bf X\) in diesem Fall aus? (Oft gibt es übrigens verschiedene Möglichkeiten dafür.)

Wir schauen nun nur noch kurz den \(t\)-Test für zwei unabhängige Stichproben an. Hier sind die Parameter, die geschätzt werden müssen \(\mu_A\) und \(\mu_B\). Ein Vektor \(\bf a\) enthalte nun die Werte der Gruppe \(A\) und ein Vektor \(\bf b\) die Werte der Gruppe \(B\). Der Vektor \(\bf y\) soll nun “untereinander” die Werte von \(\bf a\) und \(\bf b\) enthalten. Als Parametervektor setzen wir \[ \mathbf{\beta}= \begin{pmatrix} \mu_A\\ \mu_B \end{pmatrix} \] und eine geeignete Designmatrix sieht wie folgt aus: \[ \mathbf{X}= \begin{pmatrix} 1 & 0\\ \vdots & \vdots\\ 1 & 0\\ 0 & 1\\ \vdots & \vdots\\ 0 & 1\\ \end{pmatrix} \] Wir probieren direkt in R aus, ob der Parameterschätzvektor \(\mathbf{\hat{\beta}}\) die beiden Schätzer für \(\mu_A\) und \(\mu_B\) enthält, also die beiden Mittelwerte der Stichproben \(A\) und \(B\):

a <- c(1, 2, 1, 3)
b <- c(5, 6, 4, 7)
y <- c(a, b)
X1 <- rep(c(1,0),
          each = 4)
X2 <- rep(c(0,1),
          each = 4)
X <- cbind(X1, X2) # Designmatrix erstellen...
X                  # ...und anzeigen

##      X1 X2
## [1,]  1  0
## [2,]  1  0
## [3,]  1  0
## [4,]  1  0
## [5,]  0  1
## [6,]  0  1
## [7,]  0  1
## [8,]  0  1

Zur Information hier noch die beiden Mittelwerte:

mean(a)

## [1] 1.75

mean(b)

## [1] 5.5

Nun benutzen wir genau den gleichen Code wie oben (nur wir haben statt b.dach direkt beta.dach geschrieben, weil dies die allgemeinere Form ist):

teil_1 <- t(X) %*% X          # X' * X
teil_2 <- solve(teil_1)       # (X' * X)^-1 <-- die Inverse bestimme n...
teil_3 <- teil_2 %*% t(X)     # ...das dann mit X' multiplizieren...
beta.dach <- teil_3 %*% y     # ...und das alles noch mit y multiplizieren 
beta.dach                     # Tadaaaaa.... genau die Mittelwerte

##    [,1]
## X1 1.75
## X2 5.50

12.4 Ausblick

Wir haben hier einen Einblick in das ALM bekommen, welches eine Art Dach für viele spezielle Verfahren darstellt. (Tatsächlich ist selbst das ALM übrigens nur ein Spezialfall des sog. Verallgemeinerten Linearen Modells, welches wir in Teil 9 schon kennengelernt haben). Statistikprogramme handhaben diese speziellen Verfahren oft auch einfach als ALM, ohne dass wir dies merken und gesagt bekommen. Natürlich haben wir hier die Herleitung ein wenig oberflächlich betrieben und es fehlen natürlich auch noch viele Aspekte, wie z.B. Signifikanztests der Parameter.