Vergleich der Korrelationen metrisch und ordinalskalierter Daten: Ein neuer Ansatz

From Teachwiki
Jump to: navigation, search

Einführung[edit]

Beim Datenreduktionsverfahren Faktorenanalyse werden latente Variablen für eine große Anzahl von beobachteten Variablen ermittelt, die untereinander stark korreliert sind. Die Wahl der geeigneten Berechnung der Korrelation ist abhängig vom Skalenniveau. Liegen ordinalskalierte Variablen vor, hat man die Wahl zwischen den Korrelationskoeffizienten z.B. Spearman, Somers’ D, Tau-b. Nimmt man an, dass den ordinal gemessenen Variablen latente metrische und normalverteilte Merkmale zugrunde liegen, verwendet man die polychorische Korrelation.

Masterarbeit von Frau Pari-Schatz[edit]

In dieser Arbeit werden die Ansätze zum Thema Korrelationsmatrizen für metrisch- und ordinalskalierte Variablen aus der Masterarbeit von Frau Pari-Schatz weitergeführt. In ihrer Masterarbeit weist sie darauf hin, dass die Anwendung der Bravais-Pearson-Korrelation auf ordinale Daten zu Unterschätzung der Korrelation führt. Dies wird auch in den von ihr durchgeführten Simulationen deutlich. Hierfür hat sie zunächst bivariate standardnormalverteilte Variablen simuliert. Die Korrelationen zwischen den marginalen Verteilungen wurden als Vergleichswert festgelegt. Die metrischen Daten hat sie anschließend in zwei bis zehn gleich große und symmetrisch zu Null liegende Kategorien eingeteilt, wobei die Quantile der Standardnormalverteilung als Klassengrenzen dienten. Von diesen ordinalskalierten Daten wurden die Korrelationen berechnet und mit den wahren festgelegten Korrelationen der metrisch skalierten Daten verglichen.

Die polychorische Korrelation fordert, dass die underlying variable der zugrunde liegenden Verteilung standardnormalverteilt ist. Frau Pari-Schatz simulierte in einer weiteren Untersuchung bivariate \chi^2-verteilte Daten. Diese wurden ebenfalls in gleich große Klassen eingeteilt und die Korrelation der marginalen Verteilungen mit der polychorischen Korrelation verglichen. Bei einer bivariaten \chi^2-Verteilung mit jeweils vier Freiheitsgraden in den marginalen Verteilungen zeigte die polychorische Korrelation relativ gute Ergebnisse. Dahingegen bei einer bivariaten \chi^2-Verteilung mit vier und 100 Freiheitsgraden weichen die Ergebnisse für die polychorische Korrelation stark ab.

Weiterführung von Anja Weiß und Melanie Reichelt[edit]

Die Ergebnisse dieser Arbeit führten zu der Frage, wie die polychorische Korrelation auf die Verletzung der Normalverteilungsannahme reagiert bzw. wie entwickeln sich die Ergebnisse, wenn die polychorische Korrelation auf eine schiefe Verteilung angewendet wird.

In der Arbeit mit dem Thema "Vergleich der Korrelationen schiefer und gewölbter metrisch und ordinalskalierter Daten" von Frau Weiß und mir wurde zur Simulation der Verteilungen der Ansatz von Devroye (1991) verwendet. Mit Hilfe der Momente werden die Koeffizienten

a_j=\frac{2j+1}{2}\sum^{j/2}_{k=0}(-1)^k2^{-j}{j\choose k}{2j-2k\choose j}\mu_{j-2k}

bestimmt und mit Hilfe einer Rekursionsformel

(j+1)\varphi_{j+1}(x)-(2j+1)x\varphi_j(x)+j\varphi_{j-1}(x)=0

die Funktionen ermittelt. Diese ergeben dann die Verteilung

f(x)=\sum_{j=0}^{\infty} a_j\varphi_j(x).

Im Gegensatz zur eben beschriebenen Arbeit, verwende ich in dieser Arbeit einen anderen Ansatz, um Daten zu erzeugen, die einer schiefen Verteilung folgen. Hierdurch sollen die Ergebnisse verifiziert werden, da wir uns in der vorangegangenen Arbeit nicht sicher waren, ob der theoretische Ansatz von Devroye richtig umgesetzt wurde. Durch die Summation zweier gewichteter Normalverteilungen werden Mischverteilungen erzeugt, die in der englischen Literatur unter den Begriffen Mixture Models, Mixture Distributions oder Compound Distributions zu finden sind.

Simulation der Daten[edit]

Vorgehensweise[edit]

Die p-dimensionale Normalverteilung besitzt folgende Dichtefunktion

f(x)=|2\pi\Sigma|^{-1/2}\textrm{exp}\left\{-\frac{1}{2}(x-\mu)^{\top}\Sigma^{-1}(x-\mu)\right\}

mit dem Erwartungswert \mu^{\top}=(\mu_1,\mu_2,\ldots,\mu_p) und der Kovarianzmatrix \Sigma>0. Kurz schreibt man X\sim N_{p}(\mu,\Sigma)(siehe Härdle,2003).

Ein Spezialfall der multivariaten Normalverteilung ist die bivariate Normalverteilung mit \mu^{\top}=(\mu_1,\mu_2) und der Kovarianzmatrix

\mathbf{\Sigma} =
\begin{pmatrix} \sigma_{11} & \sigma_{12}\\ \sigma_{21} & \sigma_{22}\end{pmatrix}
=\begin{pmatrix} \sigma^{2}_{1} & \rho\sigma_{1}\sigma_{2}\\ \rho\sigma_{1}\sigma_{2} & \sigma^{2}_{2}\end{pmatrix}.

Die bivariate Normalverteilung hängt von fünf Parametern ab, den beiden Erwartungswerten \mu_1 und \mu_2, den beiden Varianzen \sigma_1 und \sigma_2 und dem Korrelationskoeffizienten \rho.

Abbildung 1: Mischverteilung zweier univariaten Normalverteilungen

Eine Mischverteilung g(x) besteht aus einer Linearkombination von verschiedenen Funktionen f_{j}(x)

g(x)=\sum^{J}_{j=1}r_{j}f_{j}(x)

wobei die f_{j}(x) als Basisfunktion und die r_j als Mischungsparameter bezeichnet werden. Seien nun die Basisfunktionen f_{1}(x) und f_{2}(x) zwei Normalverteilungen, d.h., J=2 und sei r=r_{1} und 1-r=r_{2}, dann hat die Mischverteilung die Form

g(x)=r_{1}f_{1}(x)+(1-r_{1})f_{2}(x).

Abbildung 1 zeigt die Mischfunktionen zweier univariater Normalverteilungen. Die linke Funktion ist bimodal und die rechte Funktion ist unimodal. In Abbildung 2 sieht man die Mischfunktionen zweier bivariater Normalverteilungen. Die linke Funktion ist wiederum bimodal und die rechte Funktion ist unimodal.

Abbildung 2: Mischverteilung zweier bivariaten Normalverteilungen
Bimodale Mischverteilung
Unimodale Mischverteilung

Umsetzung in R[edit]

Die simulierten Daten werden im ersten Schritt in 2 bis 5 gleich große Kategorien eingeteilt. Die Quantile, die als Klassengrenzen dienen, werden geschätzt. In dieser Arbeit verwende ich den Harrell-Davis-Schätzer, da in Simulationsstudien u.a. von Handl (1985) die Effizienz der Quantilschätzer hinsichtlich des mittleren quadratischen Fehlers für eine Vielzahl von Verteilungen verglichen wurden und der Harrell-Davis-Schätzer gute Ergebnisse lieferte, wenn nicht zu extreme Quantile geschätzt wurden. Bei dem Schätzer von Harrell und Davis wird x_p geschätzt, indem \mathop{\mbox{E}}[X_{((n+a)p)}] geschätzt wird. Sie gehen hierbei von einer geordneten Stichprobe x_{1},\ldots,x_{n} aus und die dazugehörigen Zufallsvariablen werden mit X_i,\;i=1,\ldots,n bezeichnet. Es gilt

\lim_{n\to\infty}\mathop{\mbox{E}}[X_{((n+a)p)}]=x_p

Als Schätzer für \mathop{\mbox{E}}[X_{((n+a)p)}] ist gegeben

\hat{\mathop{\mbox{E}}}[X_{((n+a)p)}]	=\sum^{n}_{i=1}w_{n,i}x_{(i)}

mit

w_{n,i}=\frac{\int^{i/n}_{(i-1)/n}y^{(n+1)p-1}(1-y)^{(n+1)(1-p)-1}dy} {\beta((n+1)p,(n+1)(1-p))}

und

\beta(a,b)=\int^{1}_{0}z^{a-1}(1-z)^{b-1}dz.

Den Harrell-Davis-Schätzer gewinnt man in R durch die Funktion hdquantile aus dem Paket Hmisc. So ergibt sich für 2 Kategorien das Quantil x_{0.5} als Grenze. Alle metrischen Werte, die \leq x_{0.5} sind, werden mit 1 kodiert und Werte, die  >x_{0.5} sind, werden mit 2 kodiert. Analog erfolgt die Einteilung in Klassen für die anderen Kategorien.

Von den metrischen Daten werden die Korrelationskoeffizienten berechnet und von den kategoriellen Daten werden die polychorischen Korrelationen berechnet.

Die polychorische Korrelation wird bei geordneten kategoriellen Daten verwendet. Liegen nur zwei Kategorien vor, so spricht man von der tetrachorischen Korrelation (Pearson, 1901). Man nimmt an, dass für jede ordinale Variable z eine kontinuierliche underlying variable  z^* existiert. Die Verbindung zwischen z mit m Kategorien und z^* ist

z=i\ \Leftrightarrow\ \tau_{i-1}<z^{*}<\tau_{i}\  i=1,2,\ldots,m

die Parameter

-\infty=\tau_{0}<\tau_{1}<\tau_{2}<\ldots<\tau_{m-1}<\tau_{m}=+\infty

werden als Schwellenwerte bezeichnet. Wenn die Daten in m Kategorien vorliegen, so existieren m-1 Schwellenwerte. Für die Variable z^{*} kann man im Prinzip jede stetige Verteilung wählen. Günstigerweise wählt man die Standardnormalverteilung mit der Dichtefunktion \phi(u) und der Verteilungsfunktion \Phi(u) für z^{*}. Somit ergibt sich für die geschätzten Schwellenwerte

\hat{\tau_{i}}=\Phi^{-1}(p_{1}+p_{2}+\ldots+p_{i})\;i=1,\ldots,m-1

mit \Phi^{-1} als Inverse der Standardnormalverteilung. Die polychorische Korrelation ist die Korrelation in der zweidimensionalen Normalverteilung der underlying variable z^{*}_{1} und z^{*}_{2}.

In R wird die polychorische Korrelation mit Hilfe der Funktion polychor aus dem Paket polycor bestimmt. Die Funktion polychor berechnet die polychorische Korrelation zwischen zwei ordinalskalierten Variablen unter der Annahme, dass diese bivariat normalverteilt sind. Man hat die Möglichkeit, zwischen der Maximum-Likelihood-Methode oder der Maximum-Likelihood- "Zwei-Schritt" Approximation zu wählen. Bei der "Zwei-Schritt" Approximation werden im ersten Schritt die Schwellenwerte anhand der univariaten marginalen Verteilung geschätzt. Im zweiten Schritt wird mit Hilfe der Schwellenwerte die polychorische Korrelation anhand der bivariaten marginalen Verteilungen geschätzt. Bei der Maximum-Likelihood-Methode werden die Schwellenwerte und die polychorische Korrelation in einem Schritt geschätzt. Man nimmt aber meistens die "Zwei-Schritt" Approximation, da die Schätzungen bei beiden Methoden annähernd gleich sind, die "Zwei-Schritt" Approximation aber schneller ist.

Der Vorgang der Simulation der 1000 metrischen Beobachtungen, die Einteilung in Klassen und das Berechnen der Korrelationen wird 1000 mal wiederholt. In jedem Durchgang wird das Verhältnis der beiden Korrelationen bestimmt und gespeichert. Am Ende des Programms wird dann das arithmetische Mittel über die Verhältnisse gebildet. Des Weiteren wid die Standardabweichung ausgegeben und ein Test auf \mu=1 zum Signifikanzniveau \alpha=0.05 gemacht. Ein Programmtext für R für 2 Kategorien befindet sich im Anhang. Im nächsten Schritt werden die Daten in verschieden große Klassen eingeteilt. Die Vorgehensweise anschließend ist analog.

Ergebnisse[edit]

Ergebnisse gleich großer Klassen[edit]

Um die Korrelationen der metrischen Daten mit denen der ordinalskalierten Daten zu vergleichen, wurde jeweils das Verhältnis zwischen beiden Werten gebildet. Dabei erhält man, einen Verhältniswert von 1, wenn beide Korrelationen annähernd gleich groß sind. Wenn das Verhältnis < 1, dann wird die Korrelation der metrischen Daten unterschätzt, d.h., die Korrelation der metrischen Daten ist größer als die der kategoriellen Daten. Wenn der Verhältniswert > 1, dann wird die Korrelation der metrischen Daten überschätzt. In Abbildung 4 werden die verschiedenen Kategorien miteinander verglichen. Dabei sind auf der x-Achse die unterschiedlichen Mischungsparameter r und auf der y-Achse die Mittelwerte der Verhältnisse aus den Korrelationen abgetragen.

Abbildung 3: Mittelwert bei verschiedenen Mischungsparametern , bimodale Verteilung, gleich große Klassen, 1000 Beobachtungen, 1000 Durchläufe

Man sieht, dass die Ergebnisse um 1 schwanken, d.h., die Korrelation der metrischen Daten und die polychorische Korrelation liefern annähernd gleiche Ergebnisse. Zwar weicht die Kategorie 2 ein wenig von den anderen Kategorien ab, aber insgesamt liegt die Abweichung bei weniger als 1%. Der Test auf \mu=1 bzw. die Konfidenzintervalle zeigen, dass bei Kategorie 2 und 3 die 1 meistens nicht enthalten ist, dahingegen bei Kategorie 4 und 5 immer enthalten ist (siehe Tabelle 1). Abbildung 5 ist analog zu Abbildung 4, nur dass die Daten einer unimodalen Verteilung entstammen.

Abbildung 4: Mittelwert bei verschiedenen Mischungsparametern , unimodale Verteilung, gleich große Klassen, 1000 Beobachtungen, 1000 Durchläufe

Bei dieser Abbildung weicht die Kategorie 2 ebenfalls von den anderen Ergebnissen ab, aber alle Mittelwerte sind hier größer als eins, d.h. die polychorische Korrelation ist kleiner als die metrische Korrelation. Des Weiteren wird, wie in Tabelle 1 zu sehen, die Nullhypothese \mu=1 fast immer abgelehnt.

Tabelle 1: p-Werte (gerundet) des Tests auf \mu=1, grün=größer als Signifikanzniveau, rot=kleiner als Signifikanzniveau
  Bimodal Unimodal
r 2 3 4 5 2 3 4 5
0.1 0.05 0.04 0.85 0.04 0.00 0.00 0.00 0.03
0.2 0.00 0.09 0.52 0.48 0.04 0.01 0.00 0.05
0.3 0.02 0.02 0.20 0.11 0.00 0.11 0.02 0.01
0.4 0.21 0.03 0.07 0.31 0.00 0.00 0.00 0.00

Ergebnisse verschiedener Klassengrößen[edit]

Analog zu den Daten, die in gleich große Klassen kategorisiert wurden, habe ich für die Korrelationen der Daten, die in unterschiedlich große Klassen kategorisiert wurden, den Mittelwert über die Verhältnise der Korrelationen gebildet. Die Abbildung 6 zeigt die Ergebnisse, in der die Korrelation auf 2 künstlich kodierte Variablen, denen jeweils eine unimodale Mischverteilung zugrunde liegt, mit 2 Kategorien angewandt wurde. In der Legende in der Abbildung 6 ist angegeben, welche Quantilsgrenzen gewählt wurden, um die Klassengrößen zu bilden.

Abbildung 5: Mittelwert bei verschiedenen Mischungsparametern , unimodale Verteilung, verschieden große Klassen, 1000 Beobachtungen, 1000 Durchläufe

Bei dieser Abbildung wird sichtbar, dass die Abweichung zwischen der Korrelation der metrischen Daten und der polychorischen Korrelation größer ist, als bei der Einteilung in gleich große Klassen. Je näher sich die Einteilung der symmetrischen Form nähert, um so besser werden die Ergebnisse. Die 1 wird aber nie in den Konfidenzintervalle eingeschlossen. Abschließend sieht man in der Abbildung 7 die Ergebnisse der Simulation für 3 Kategorien und asymmetrischer Klassenaufteilung.

Abbildung 6: Mittelwert bei verschiedenen Mischungsparametern , unimodale Verteilung, verschieden große Klassen, 1000 Beobachtungen, 1000 Durchläufe

Diese Abbildung ähnelt der Abbildung 6. Die Einteilung 80%-10%-10% weicht von den anderen Ergebnissen leicht ab, aber alle Mittelwerte sind größer als eins, d.h. die polychorische Korrelation ist kleiner als die metrische Korrelation. Auch in dieser Simulation wird die 1 nie in die Konfidenzintervalle eingeschlossen.

Zusammenfassung[edit]

Die Ergebnisse dieser Arbeit stehen im Gegensatz zu den erwarteten Ergebnissen. Eigentlich hat man erwartet, dass die polychorische Korrelation stärker auf die Abweichung von der Normalverteilung reagiert. Dies hatte sich teilweise auch in den vorangegangenen Arbeiten gezeigt. Eine mögliche Erklärung hierfür wäre, dass die erzeugten Verteilungen eine zu geringe Bandbreite an Schiefe aufweisen, d.h., zu nah an der Normalverteilung liegen. Um dies zu überprüfen, wurde auf die simulierten Daten der Shapiro-Wilk-Test für multivariate Normalverteilung, auch Royston’s H-Test genannt, angewandt. Der Shapiro-Wilk-Test testet die Hypothese, ob die Daten normalverteilt sind. Die Nullhypothese besagt, dass eine Normalverteilung vorliegt und die Alternativhypothese besagt, dass keine Normalverteilung vorliegt. Die Ergebnisse dieser Arbeit stehen im Gegensatz zu den erwarteten Ergebnissen. Eigentlich hat man erwartet, dass die polychorische Korrelation stärker auf die Abweichung von der Normalverteilung reagiert. Dies hatte sich teilweise auch in den vorangegangenen Arbeiten gezeigt. Eine mögliche Erklärung hierfür wäre, dass die erzeugten Verteilungen eine zu geringe Bandbreite an Schiefe aufweisen, d.h., zu nah an der Normalverteilung liegen. Um dies zu überprüfen, wurde auf die simulierten Daten der Shapiro-Wilk-Test für multivariate Normalverteilung, auch Royston's H-Test genannt, angewandt. Der Shapiro-Wilk-Test testet die Hypothese, ob die Daten normalverteilt sind. Die Nullhypothese besagt, dass eine Normalverteilung vorliegt und die Alternativhypothese besagt, dass keine Normalverteilung vorliegt. Betrachtet werden die Ordnungsstatistiken x_{(1)},\ldots,x_{(n)} einer geordneten unvariaten Stichprobe. Die Teststatistik von Shapiro und Wilk(1965) lautet:

\tilde{T}(x_{1},\ldots,x_{n})=\frac{(\sum^{n}_{i=1}a_{i}x_{(i)})^2}{\sum^{n}_{i=1}(x_{i}-\bar{x})^2}

Die Gewichte a_i sind gebeben durch

(a_{1},\ldots,a_{n})=\frac{m^{\top}V^{-1}}{(m^{\top}V^{-1}V^{-1}m)^{1/2}}

mit m=(m_{1},\ldots,m_{n})^{\top} als den Erwartungswerten der Ränge und V als Kovarianzmatrix der Ränge. Vorausgesetzt x_{1},\ldots,x_{n} stammen aus einer Normalverteilung, schlägt Royston folgende Normaltransformation vor

Z={(1-W)^{\lambda}-\mu}/\sigma,

wobei die Aprroximation durch Polynome von \lambda, \mu und \sigma bei Royston (1983) näher beschrieben wird. Diese Idee erweitert er dann auf den multivariaten Fall. W_j ist dann die univariate Shapiro-Wilk-Teststatistik für die j-te Variable in einer p-dimensionalen Verteilung. Wird die Transformation auf jedes W_j angewendet, ergeben sich die Z_j. Mit Hilfe der Z_j wird dann

K_{j}=[\Phi^{-1}{\Phi(-Z_{j})/2}]^2

berechnet, wobei \Phi(\cdot) das Quantil der Standardnormalverteilung ist. Wenn die Daten multivariat normalverteilt sind, dann ist die Teststatistik

H=e\sum^{p}_{j=1}K_{j}

approximativ \chi^{2}_{e}-verteilt. e wird equivalent degrees of freedom genannt und kann durch folgende Formel geschätzt werden

\hat{e}=p/{1+(p-1)\bar{c}}

mit \bar{c} als Schätzer für die mittlere Korrelation zwischen den K_{j}. In R kann der Test mit Hilfe der Funktion mshapiro.test() aus dem Paket mvnormtest ausgeführt werden. Wenn der ausgegebene p-Wert <0.05 ist, sind die Daten statistisch signifikant nicht normalverteilt (auf dem Signifikanzniveau \alpha=0.05). Es zeigt sich bei den von mir simulierten Daten, dass in nur ca. 10% der Fälle die Nullhypothese abgelehnt werden kann. Ein weiteres Problem ist die bereits angesprochene geringe Spannweite der Schiefe. In einer weiteren Arbeit müsste z.B. durch die Addition mehrer Normalverteilungen oder durch einen neuen Ansatz eine größere Spannweite der Schiefe erreicht werden und dann die Simulationen wiederholt werden.


Literatur[edit]

  • Mischverteilungen, Kapitel 9
  • Fox, J.: The polycor Package, R Software
  • Handl, A.: Verteilungsfreie Quantilschätzer - ein Vergleich. Diskussionsarbeit 1985 Nr. 6/85, Institut für Quantitative Ökonomik und Statistik, Freie Universität Berlin
  • Handl, A.; Quantile schätzen in R
  • Härdle, W. und Simar, L.: Applied Multivariate Statistical Analysis. Springer, 2003
  • Royston, J.P.: Some Techniques for Assessing Multivariate Normality Based on the Shapiro-Wilk W, Applied Statistics, vol. 32, 1983
  • Prof. Dr. Schmidt, V.: Statistik II. Universität Ulm, 2003, Skript
  • Svantesson, T. und Wallac, J.W.: Tests for assessing multivariate normality and the covariance structure of MIMO data Paper

Anhang[edit]

simulate <- function(n,r1)
{
r1 <- r1
r2<-1-r1
data1 <- rmvnorm(n, c(0, 0), matrix(c(2, .5, .5, 2), 2, 2))
data2 <- rmvnorm(n, c(1.5, 1.5), matrix(c(2, .5, .5, 2), 2, 2))
dat<<-r1*data1+r2*data2

}
# !!! Ende Makro Simulate !!!!

vecp2<<- NULL 
k<<-NULL
j=1
while(j<=1000)
{
n <- 1000
simulate(n,0.4)
C <- t(dat[1:1000,1:2])
st<-mshapiro.test(C)
s<-st2
if(s < 0.05)
	{
		k <- c(k,1)
	} else
	{
		k <- c(k,0)
	}

x <- dat[,1]
y <- dat[,2]
correlation <-cor(x, y)  # sample correlation
q <- hdquantile(x,0.5)
q2<- hdquantile(y,0.5)
x <- cut(x, c(-Inf, q1, Inf))
y <- cut(y, c(-Inf, q21, Inf))
polycor<-polychor(x, y)  # 2-step estimate
pc<- correlation/polycor
vecp2 <- c(vecp2,pc)
j <- j+1
}
vecp2
m <- mean(vecp2)
v <-sqrt(var(vecp2))
sum(k)
m
v