# R-Syntax zu den Folien "Zusammenhangsanalysen nominaler Daten"
# (Nominale_Daten.pdf) und "weitere Korrelationskoeffizienten"
# (Korrelationskoeffizienten.pdf) (Statistik für Kriminologen, WS 2006/2007)

library(polycor) # Paket zur Berechnung polyserialer und polychorischer
                 # Korrelationskoeffizienten (setzt Pakt 'mvtnorm' voraus;
                 # beide Pakete müssen zuvor extra installiert worden sein)

BaseURL="http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Lehre/StatIKrim/Uebungen/"

# Folgende Pfadangabe bitte anpassen (dabei den Vorwärts-Slash (/) als letztes
# Zeichen nicht vergessen!):
Pfad=''

# Nur bei Internetanschluss verwenden:
download.file(paste(BaseURL,'crosstabs.r',sep=""),paste(Pfad,'crosstabs.r',sep=""))
download.file(paste(BaseURL,'JugGew99.RData',sep=""),paste(Pfad,'JugGew99.RData',sep=""),mode="wb")
download.file(paste(BaseURL,'kurtosis.r',sep=""),paste(Pfad,'kurtosis.r',sep=""))
download.file(paste(BaseURL,'skewness.r',sep=""),paste(Pfad,'skewness.r',sep=""))

# Die Dateien 'crosstabs.r','JugGew99.RData','kurtosis.r' und 'skewness.r'
# müssen auf der Festplatte vorliegen:
source(paste(Pfad,'crosstabs.r',sep=''))
load(paste(Pfad,'JugGew99.RData',sep=''))
source(paste(Pfad,'kurtosis.r',sep=''))
source(paste(Pfad,'skewness.r',sep=''))

# A) Funktionen für Zusammenhangsanalysen nominaler Daten:
# ------------------------------------------------------------------------------

cont.C = function(x,y=0) # x=table(x,y) or x,y=vectors
{
   if (sum(as.numeric(y))==0) {
   chi = chisq.test(x,correct=F)
   }
   else {
   chi = chisq.test(x,y,correct=F)
   }
   n = sum(chi$obs)
   C = sqrt(chi$stat/(n+chi$stat))
   names(C)="Continguency coefficient C"
#   cat('Continguency coefficient C =',round(V,3),'\n')
   C
}

# ------------------------------------------------------------------------------
# Im Fall von 2*2-Kreuztabellen ist Cramérs V gleich dem Absolutbetrag von phi
# und gleich dem Absolutbetrag der Produkt-Moment-Korrelation:

Cramers.V = function(x,y=0) # x=table(x,y) or x,y=vectors
{
   if (sum(as.numeric(y))==0) {
   chi = chisq.test(x,correct=F)
   }
   else {
   chi = chisq.test(x,y,correct=F)
   }
   n = sum(chi$obs)
   k = min(dim(chi$obs))
   V = sqrt(chi$stat/(n*(k-1)))
   names(V)="Cramér's V"
#   cat('Cramérs V =',round(V,3),'\n')
   V
}

# ------------------------------------------------------------------------------
# Man beachte, dass phi = Produkt-Moment-Korrelation für den Fall, dass die
# Daten in Vektorform vorliegen. Falls die Daten als Tabelle vorliegen, wird
# nur der Absolutbetrag des phi-Koeffizienten ausgegeben:

phi = function(x,y=0) # x=table(x,y) or x,y=vectors
{
   if (sum(as.numeric(y))==0) {
     chi = chisq.test(x,correct=F)
     n = sum(chi$obs)
     phi = sqrt(chi$stat/n)
     names(phi)="|phi|"
#     cat('|phi| =',round(phi,3),'\n')
   }
   else {
     phi = cor(x,y)
     names(phi)="phi"
#     cat('phi =',round(phi,3),'\n')
   }
   phi
}

# B) Funktionen für "weitere Korrelationskoeffizienten":
# ------------------------------------------------------------------------------

sens.spec = function(x,y=0) # x=group, y=test; x=table(x,y) or x,y=vectors
{
   if (sum(as.numeric(y))==0) {
     a = x[1]
     c = x[2]
     b = x[3]
     d = x[4]
   }
   else {
     tab = table(y,x)
     a = tab[1]
     c = tab[2]
     b = tab[3]
     d = tab[4]
   }
   sensitivity = d/(b+d)
   specificity = a/(a+c)
#   cat('Sensitiviy = ',round(sensitiviy,3),
#       ', Specificity = ',round(specificity,3),'\n',sep='')
   list(Sensitivity=sensitivity,Specificity=specificity)
}

# ------------------------------------------------------------------------------

OR = function(x,y=0) # x=table(x,y) or x,y=vectors
{
   if (sum(as.numeric(y))==0) {
     a = x[1]
     c = x[2]
     b = x[3]
     d = x[4]
   }
   else {
     tab = table(y,x)
     a = tab[1]
     c = tab[2]
     b = tab[3]
     d = tab[4]
   }
   OddsRatio = (a*d)/(b*c)
#   cat('Odds-Ratio =',round(OddsRatio,3),'\n')
   OddsRatio
}

# ------------------------------------------------------------------------------

Yules.Y = function(x,y=0) # x=table(x,y) or x,y=vectors
{
   if (sum(as.numeric(y))==0) {
     a = x[1]
     c = x[2]
     b = x[3]
     d = x[4]
   }
   else {
     tab = table(y,x)
     a = tab[1]
     c = tab[2]
     b = tab[3]
     d = tab[4]
   }
   OR = (a*d)/(b*c)
   Y = (sqrt(OR)-1)/(sqrt(OR)+1)
#   cat('Yules Y =',round(Y,3),'\n')
   Y
}

# ------------------------------------------------------------------------------

nu = function(x,y=0) # x=table(x,y) or x,y=vectors
                     # x = naturally dichotomous, y = artificially dichotomous
{
   if (sum(as.numeric(y))==0) {
     p00 = x[1]/sum(x)
     p10 = x[2]/sum(x)
     p01 = x[3]/sum(x)
     p11 = x[4]/sum(x)
   }
   else {
     tab = table(x,y)
     p00 = tab[1]/sum(tab)
     p10 = tab[2]/sum(tab)
     p01 = tab[3]/sum(tab)
     p11 = tab[4]/sum(tab)
   }
   p0. = p00+p01
   p1. = p10+p11
   delta = qnorm(p00/p0.)-qnorm(p10/p1.)
   nu.coefficient = delta/sqrt(delta^2 + 1/(p1.*(1-p1.)))
#   cat('nu-coefficient =',round(nu.coefficient,3),'\n')
   nu.coefficient
}

# ------------------------------------------------------------------------------
# Achtung: Keine Fehlerchecks für die Funktion rank.bis.corr!

rank.bis.corr = function(x,y)   # x = ordinal, y = dichotomous
{
  xrank = rank(x)
  n = length(x)
  rank.b.r = 2/n*(mean(xrank[y==max(y)])-mean(xrank[y==min(y)]))
#  cat('rank-biserial r = ',round(rank.b.r,3),'\n')
  rank.b.r
}

# ==============================================================================
# A) Zusammenhangsanalysen nominaler Daten:

# 1. Hypothetisches Beispiel: Perfekte Unabhängigkeit zwischen Studienfach-
#    präferenz und Geschlecht

Tabelle.1 = as.table(matrix(c(18,27,36,54,18,27),nrow=2,ncol=3,
                     dimnames=list(Geschlecht =c('weiblich','männlich'),
                                   Studienfach=c('Sprach- & Kult.wiss.',
                                                 'Rechts-,Wi-,Soz.wiss.',
                                                 'Mathem. & Nat.wiss.'))))
crosstabs(Tabelle.1)

# ------------------------------------------------------------------------------
# 2. Fiktives Beispiel: Zusammenhang zw. Studienfachpräferenz und Geschlecht:

Tabelle.2 = as.table(matrix(c(28,17,31,59,13,32),nrow=2,ncol=3,
                     dimnames=list(Geschlecht =c('weiblich','männlich'),
                                   Studienfach=c('Sprach- & Kult.wiss.',
                                                 'Rechts-,Wi-,Soz.wiss.',
                                                 'Mathem. & Nat.wiss.'))))
crosstabs(Tabelle.2,expected=T,chisq=T)
cont.C(Tabelle.2)
Cramers.V(Tabelle.2)

# ------------------------------------------------------------------------------
# 3. Verurteilung und Geschlecht:

Tabelle.3 = as.table(matrix(c(1869739,1461,1883749,91251),nrow=2,ncol=2,
                     dimnames=list(verurteilt=c('nein','ja'),
                                   Geschlecht=c('weiblich','männlich'))))
crosstabs(Tabelle.3)
phi(Tabelle.3)

#Prozent weiblicher Jugendlicher verurteilt:
100*Tabelle.3[2]/(Tabelle.3[1]+Tabelle.3[2])

#Prozent männlicher Jugendlicher verurteilt:
100*Tabelle.3[4]/(Tabelle.3[3]+Tabelle.3[4])

# ==============================================================================
# B) "weitere Korrelationskoeffizienten":

# 1a) Zusammenhang zwischen Test und Gesundheit (Männer):

Tabelle.4 = as.table(matrix(c(518,182,117,183),nrow=2,ncol=2,
                     dimnames=list(Test      = c('negativ','positiv'),
                                   Krankheit = c('gesund','krank'))))
crosstabs(Tabelle.4)
phi(Tabelle.4)
sens.spec(Tabelle.4)

# 1b) Zusammenhang zwischen Test und Gesundheit (Frauen):

Tabelle.5 = as.table(matrix(c(666,234,39,61),nrow=2,ncol=2,
                     dimnames=list(Test      = c('negativ','positiv'),
                                   Krankheit = c('gesund','krank'))))
crosstabs(Tabelle.5)
phi(Tabelle.5)
sens.spec(Tabelle.5)

# Basisrate Männer:
(Tabelle.4[3]+Tabelle.4[4])/sum(Tabelle.4)

# Basisrate Frauen:
(Tabelle.5[3]+Tabelle.5[4])/sum(Tabelle.5)

# ------------------------------------------------------------------------------
# 2a) Odds-Ratio Krankheit bei positivem vs. negativem Test (Männer):

OR(Tabelle.4)

# 2b) Odds-Ratio Krankheit bei positivem vs. negativem Test (Frauen):

OR(Tabelle.5)

# ------------------------------------------------------------------------------
# 3a) Yules-Y für Krankheit und Testergebnis (Männer):

Yules.Y(Tabelle.4)

# 3b) Yules-Y für Krankheit und Testergebnis (Frauen):

Yules.Y(Tabelle.5)

# ------------------------------------------------------------------------------
# 4) Grafik des Anteils der Korrelation in Abhängigkeit vom Teilungspunkt (eine
#    kontinuierliche und eine dichotome Variable unter der Annahme normalver-
#    teilter manifester und latenter Variablen):

windows(width=9.5,record=T)

curve(dnorm(qnorm(x))/sqrt(x*(1-x)),from=0,to=1,n=10001,ylim=c(0,1),col='blue',
      ylab='Proportionale Reduktion von r',
      xlab='Anteil der Population unter dem Teilungspunkt',
      main='Korrelation und Dichotomisierung kontinuierlicher Variablen')
grid()

# ------------------------------------------------------------------------------
# 5) Erzeugung normalverteilter Zufallszahlen für X und Y mit einer Korrelation
#    von .57 (1 Millionen Fälle, dem Computer bitte Zeit lassen!):

data = rmvnorm(1000000, c(0, 0), matrix(c(1, .57, .57, 1), 2, 2))
x = data[,1]
y = data[,2]

# Korrelation der kontinuierlichen Variablen:
cor(x,y)

# Dichotomisierung der Y-Variable beim Quantil 1/3:
p = 1/3
y.d = cut(y,c(-Inf,quantile(y,p),+Inf))

# Korrelation der kontinuierlichen X- mit der dichotomisierten Y-Variablen
# (Produkt-Moment-Korrelation = punkt-biseriale Korrelation):
cor(x,y.d)

# Bestimmung des Reduktionsfaktors:
z = qnorm(p)
h = dnorm(z)
red.fakt = h/sqrt(p*(1-p))

# Schätzung der Korrelation von X mit der latenten kontinuierlichen Variablen Y
# (die Korrelation der kontinuierlichen Variablen sollte reproduziert werden):
cor(x,y.d)/red.fakt

# Überprüfung, ob die Funktion polyserial() das gleiche Ergebnis liefert (das
# Paket 'polycor' muss zuvor installiert und aktiviert worden sein):
polyserial(x,y.d)

# ------------------------------------------------------------------------------
# 6) Grafik zur Bestimmung von h (Ordinate der Standardnormalverteilung) für
#    die Teilung der Population am Quantil 1/3:

par(las=0)
curve(dnorm(x),-3,3,ylab='Dichte (h)',xlab='z',
      main='Standardnormalverteilung',col='blue',lwd=2,font=2)
lines(x=c(qnorm(1/3),qnorm(1/3)),y=c(-1,dnorm(qnorm(1/3))),
      col='red',lty=1,lwd=2)
lines(x=c(qnorm(1/3),-2.8),y=c(dnorm(qnorm(1/3)),dnorm(qnorm(1/3))),
      col='red',lty=2,lwd=2)
text(-0.8,0.05,'p = 0.3333',pos=2,font=2,cex=1.2)
text(qnorm(1/3)+abs(-0.8-qnorm(1/3)),0.05,'q = (1-p) = 0.6667',
     pos=4,font=2,cex=1.2)
mtext(paste('z =',round(qnorm(1/3),3)),side=1,at=qnorm(1/3),
      adj=0.23,col='red',cex=1.2,font=2,padj=0.1)
par(las=1)
mtext(paste('h = ',round(dnorm(qnorm(1/3)),3),' ',sep=''),side=2,
      at=dnorm(qnorm(1/3)),col='red',cex=1.2,font=2,adj=0.34)
par(las=0)

# ------------------------------------------------------------------------------
# 7) Beispiel: Korrelation Gewalteinstellung (GE) und Selbstkontrolle (SK)
#    (da hier nicht exakt die gleiche Stichprobe vorliegt wie im Skriptbeispiel,
#    fallen die Ergebnisse geringfügig anders aus):

load(paste(Pfad,'JugGew99.RData',sep=''))
attach(JugGew99)

# Statistische Kennwerte für GE und SK:
mean(GE);sd(GE)
kurtosis(GE);skewness(GE)

mean(SK);sd(SK)
kurtosis(SK);skewness(SK)

# Produkt-Moment-Korrelation GE und SK:
cor(GE,SK)

# Dichotomisierung der Selbstkontrolle (niedrig: <= 50; hoch: > 50):
SK.d50 = factor(SK > 50,labels=c('SK <= 50','SK > 50'))

# Anzahl der Fälle mit niedriger und hoher Selbstkontrolle:
table(SK.d50)

# Streudiagramm von GE gegen SK mit Trennlinie für SK.d50:
plot(GE~SK,ylim=c(0,100),xlim=c(0,100),
     ylab='Gewalteinstellung',xlab='Selbstkontrolle')
abline(v = 50, col='blue')

# Produkt-Moment- (punkt-biseriale) Korrelation von GE und SK.d50:
cor(GE,as.numeric(SK.d50))

# Polyseriale (biseriale) Korrelation von GE mit SK.d50 (Achtung: die kontinu-
# ierliche Variable muss als erste spezifiziert werden!):
polyserial(GE,SK.d50)

# Dichotomisierung der Gewalteinstellung (niedrig: < 50; hoch >= 50):
GE.d50 = factor(GE >= 50,labels=c('GE < 50','GE >= 50'))

# Kreuztabelle von dichotomisierter Selbstkontrolle und dichotomisierter
# Gewalteinstellung (Achtung: die Reihenfolge von niedriger und hoher Gewalt-
# einstellung ist umgekehrt dargestellt als im Streudiagramm und im Skript):
table(GE.d50,SK.d50)

# Streudiagramm mit Trennlinien für GE.d50 und SK.d50:
plot(GE~SK,ylim=c(0,100),xlim=c(0,100),
     ylab='Gewalteinstellung',xlab='Selbstkontrolle')
abline(v = 50, col='blue', lty=2)
abline(h = 50, col='blue', lty=2)

# Produkt-Moment-Korrelation (phi) von GE.d50 und SK.d50:
cor(as.numeric(GE.d50),as.numeric(SK.d50))

# Polychorische (tetrachorische) Korrelation von GE.d50 mit SK.d50:
polychor(GE.d50,SK.d50)

# Dichtomisierung der Selbstkontrolle für SK > 33.33:
SK.d33 = factor(SK > 100/3,labels=c('SK <= 33.33','SK > 33.33'))

# Kreuztabelle von dichotomisierter Selbstkontrolle und dichotomisierter
# Gewalteinstellung (Achtung: die Reihenfolge von niedriger und hoher Gewalt-
# einstellung ist umgekehrt dargestellt als im Streudiagramm und im Skript):
table(GE.d50,SK.d33)

# Streudiagramm mit Trennlinien für GE.d50 und SK.d33:
plot(GE~SK,ylim=c(0,100),xlim=c(0,100),
     ylab='Gewalteinstellung',xlab='Selbstkontrolle')
abline(v = 100/3, col='red', lty=2)
abline(h = 50, col='red', lty=2)

# Produkt-Moment-Korrelation (phi) von GE.d50 und SK.d33:
cor(as.numeric(GE.d50),as.numeric(SK.d33))

# Polychorische (tetrachorische) Korrelation von GE.d50 mit SK.d33:
polychor(GE.d50,SK.d33)

# ------------------------------------------------------------------------------
# 8. Rangtransformation der Werte 5,9,9,9,12,17,17,25:

x = c(5,9,9,9,12,17,17,25)
rank(x)

# ------------------------------------------------------------------------------
# 9. Produkt-Moment-Korrelation, Spearmans rho, Kendalls tau und polychorische
# Korrelation:

id = 1:7
x  = c(10,14,2,7,19,15, 9)
y  = c(11,17,5,1,18,16,13)

# Tabelle der Rohwerte:
rbind(id,x,y)

# Produkt-Moment-Korrelation von X und Y:
cor(x,y)

# Tabelle der rangtransformierten Werte sortiert nach X:
rbind(id,'x rangtr.='=rank(x),'y rangtr.'=rank(y))[,order(rank(x))]

# Produkt-Moment-Korrelation der rangtransformierten Werte (= Spearmans rho)
# (identisch mit Spearmans rho der untransformierten Werte):
cor(rank(x),rank(y))
cor(x,y,method='spearman')

# Kendalls tau:
cor(x,y,method='kendall')

# Polychorische Korrelation:
polychor(x,y)

# ------------------------------------------------------------------------------
# 10. Nue-Koeffizient:

# Annahme: Beispieltabelle basiert auf 100 Fällen (Achtung: Die Zeilen-Variable
# muss natürlich und die Spaltenvariable künstlich dichotom sein! Im Skriptbei-
# spiel definiert aus anschaulichen Gründen die X-Variable die Spalten, hier
# müssen es die Zeilen sein):
Tabelle.6 = as.table(matrix(c(15,10,25,50),nrow=2,ncol=2,
                     dimnames=list(Geschlecht = c('weiblich','männlich'),
                                   Risiko     = c('niedrig' ,'hoch'))))
Tabelle.6

# Nue-Koeffizient:
nu(Tabelle.6)

# Betrag von Phi (= Betrag der Produkt-Moment-Korrelation):
phi(Tabelle.6)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Das gleiche nicht mit Tabellendaten sondern mit X- und Y-Variablen (Vektoren)
# (Achtung: Die X-Variable muss natürlich und die Y-Variable künstlich dichotom
# sein! Im Skriptbeispiel definiert aus anschaulichen Gründen die X-Variable die
# Spalten, hier müssen es die Zeilen sein):
x = c(rep(0,40),rep(1,60))
y = c(rep(0,15),rep(1,25),rep(0,10),rep(1,50))
table(x,y)

# x- und y lassen sich zu sogenannten Faktor-Variablen umdefinieren (für das
# Beispiel nicht unbedingt nötig, aber eine hilfreiche Demonstration):

Geschlecht = factor(x,labels=c('weiblich','maennlich'))
Risiko = factor(y,labels=c('niedrig','hoch'))
table(Geschlecht,Risiko)

# Nue-Koeffizient (Achtung: zuerst muss die natürliche, dann die künstlich di-
# chotome Variable spezifiziert werden!):
nu(Geschlecht,Risiko)

# Produkt-Moment-Korrelation (= phi):
cor(as.numeric(Geschlecht),as.numeric(Risiko))

# ==============================================================================
# C) Ergänzung
#
# In den Folien fehlt ein Beipiel zur der biserialen Rangkorrelation. Die
# Berechnung ist einfach: Die biseriale Rangkorrelation ist die Differenz der
# durchschnittlichen Ränge für die beiden Gruppen der dichotomen Variablen
# multipliziert mit zwei und relativiert an der Stichprobengröße:
#
# r(rank-biserial) = 2/n * [mean(rank2) - mean(rank1)]
#
# Sie kann auch als Effektgrößenmaß für den Mann-Whitney Test interpretiert
# (und auch anhand der U-Werte des Mann-Whitney Tests bestimmt) werden.
#
# Beispiel: X sei eine ordinal skalierte Variable, y eine natürlich dichotome
# Variable. Die Werte seien:
#
#  Fall |  x-Wert | y-Wert |   Rang(x)  | Rang für y(1) | Rang für y(2)
# ------|---------|--------|------------|---------------|--------------
#     1 |      9  |     0  |      10.5  |         10.5  |
#     2 |      5  |     1  |       4.0  |               |         4.0
#     3 |      7  |     0  |       6.5  |          6.5  |
#     4 |      8  |     0  |       8.0  |          8.0  |
#     5 |      1  |     1  |       1.0  |               |         1.0
#     6 |      5  |     1  |       4.0  |               |         4.0
#     7 |      9  |     0  |      10.5  |         10.5  |
#     8 |      4  |     1  |       2.0  |               |         2.0
#     9 |      5  |     0  |       4.0  |          4.0  |
#    10 |      9  |     0  |      10.5  |         10.5  |
#    11 |      9  |     1  |      10.5  |               |        10.5
#    12 |      7  |     1  |       6.5  |               |         6.5
# ------|---------|--------|------------|---------------|--------------
#  mean |         |        |            |  R1 =  8.333  |  R2 = 4.667
#
# r(rank-biserial) = 2/12 * (4.667 - 8.333) = -0.611
#
# ------------------------------------------------------------------------------
#
# Daten:
x = c(9,5,7,8,1,5,9,4,5,9,9,7)
y = c(0,1,0,0,1,1,0,1,0,0,1,1)

# Biseriale Rangkorrelation
# (Achtung: y muss die dichotome Variable darstellen!):
rank.bis.corr(x,y)