/* Stata-Syntax zu den Folien "Zusammenhangsanalysen nominaler Daten"
   (Nominale_Daten.pdf) und "weitere Korrelationskoeffizienten"
   (Korrelationskoeffizienten.pdf) (Statistik-1 für Kriminologen) */

* ------------------------------------------------------------------------------
* Präliminarien:

/* Bitte die Pfadangabe in der folgenden Zeile auf die eigenen Bedingungen
   anpassen: */
global pfad = "d:\lehre\statistik1\ws1112\"

* URL zum Herunterladen von Dateien:
global BaseURL = "http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Lehre/StatIKrim/"

version 11.2
set more off

* Installiere falls nötig fehlende .ado-Files:
cap which fre
if _rc ssc install fre
cap which moments2
if _rc ssc install moments2
cap which polychoric
if _rc net install polychoric.pkg, from("http://web.missouri.edu/~kolenikovs/stata/")
cap which somersd
if _rc ssc install somersd

* Generiere eine Textdatei für die Ausgabe:
cap log close
log using "${pfad}Korrelationen.log", text replace

* ==============================================================================
* A) Programme für Zusammenhangsanalysen nominaler Daten:

* ------------------------------------------------------------------------------
* Kontingenzkoeffizient C:

cap program drop cont_c
program define cont_c, rclass
  version 8.2
  tempvar N chi2 c
  syntax varlist(min=2 max=2 numeric) [if] [in]
  marksample touse
  qui tabulate `varlist' if `touse', chi2
  scalar `N' = r(N)
  scalar `chi2' = r(chi2)
  scalar `c' = sqrt(`chi2'/(`N'+`chi2'))
  di as txt "(obs=" `N' ")" _n _n ///
     as txt "Contingency coefficient C = " as res `c'
  return scalar C = `c'
  return scalar chi2 = `chi2'
  return scalar N = `N'
end

* ------------------------------------------------------------------------------
* phi-Koeffizient:

/* Man beachte, dass der phi-Koeffizient nur für dichotome Daten berechnet
   werden kann! */

cap program drop phi
program define phi, rclass
  version 8.2
  tempvar N chi2 phi
  syntax varlist(min=2 max=2 numeric) [if] [in]
  marksample touse
  qui tabulate `varlist' if `touse', chi2 V
  if (r(r) != 2 | r(c) != 2) {
    di as err "Error: All variables must be dichotomous!"
    error 420
  }
  scalar `N' = r(N)
  scalar `chi2' = r(chi2)
  scalar `phi' = r(CramersV)
  di as txt "(obs=" `N' ")" _n _n ///
     as txt "phi = " as res `phi'
  return scalar phi = `phi'
  return scalar chi2 = `chi2'
  return scalar N = `N'
end

* ------------------------------------------------------------------------------
* B) Programme für "weitere Korrelationskoeffizienten":

cap program drop sens_spec
program define sens_spec, rclass
  version 8.2
  tempname N cells p00 p10 p01 p11 p0_ p1_ sens spec
  syntax varlist(min=2 max=2 numeric) [if] [in]
  marksample touse
  qui tabulate `varlist' if `touse', matcell(`cells')
  if (r(r) != 2 | r(c) != 2) {
    di as err "Error: All variables must be dichotomous!"
    error 420
  }
  scalar `N' = r(N)
  scalar `p00' = `cells'[1,1]
  scalar `p01' = `cells'[1,2]
  scalar `p10' = `cells'[2,1]
  scalar `p11' = `cells'[2,2]
  scalar `sens' = `p11'/(`p01'+`p11')
  scalar `spec' = `p00'/(`p00'+`p10')
  di as txt "(obs=" `N' ")" _n _n ///
     as txt "Sensitivity = " as res `sens' _n ///
     as txt "Specificity = " as res `spec'
  return scalar spec = `spec'
  return scalar sens = `sens'
  return scalar N = `N'
end

* ------------------------------------------------------------------------------
* Odds Ratio:

cap program drop oddsratio
program define oddsratio, rclass
  tempname N or
  syntax varlist(min=2 max=2 numeric) [if] [in]
  marksample touse
  qui logit `varlist' if `touse'
  local iv : word 2 of `varlist'
  scalar `N' = e(N)
  scalar `or' = exp(_b[`iv'])
  di as txt "(obs=" `N' ")" _n _n ///
     as txt "Odds ratio = " as res `or'
  return scalar or = `or'
  return scalar N = `N'
end

* ------------------------------------------------------------------------------
* Yules Y:

cap program drop yules_y
program define yules_y, rclass
  tempname N or yulesy
  syntax varlist(min=2 max=2 numeric) [if] [in]
  marksample touse
  qui logit `varlist' if `touse'
  local iv : word 2 of `varlist'
  scalar `N' = e(N)
  scalar `or' = exp(_b[`iv'])
  scalar `yulesy' = (sqrt(`or')-1)/(sqrt(`or')+1)
  di as txt "(obs=" `N' ")" _n _n ///
     as txt "Yules Y = " as res `yulesy'
  return scalar YulesY = `yulesy'
  return scalar N = `N'
end

* ------------------------------------------------------------------------------
* nu-Koeffizient:

cap program drop nucoef
program define nucoef, rclass
  version 9.0
  tempname N cells p00 p10 p01 p11 p0_ p1_ delta nu
  syntax varlist(min=2 max=2 numeric) [if] [in]
  marksample touse
  qui tabulate `varlist' if `touse', matcell(`cells')
  if (r(r) != 2 | r(c) != 2) {
    di as err "Error: All variables must be dichotomous!"
    error 420
  }
  scalar `N' = r(N)
  scalar `p00' = `cells'[1,1]/`N'
  scalar `p01' = `cells'[1,2]/`N'
  scalar `p10' = `cells'[2,1]/`N'
  scalar `p11' = `cells'[2,2]/`N'
  scalar `p0_' = `p00'+`p01'
  scalar `p1_' = `p10'+`p11'
  scalar `delta' = invnormal(`p00'/`p0_')-invnormal(`p10'/`p1_')
  scalar `nu' = `delta'/sqrt(`delta'^2 + 1/(`p1_'*(1-`p1_')))
  local vnat : word 1 of `varlist'
  local vart : word 2 of `varlist'
  di as txt "(obs=" `N' ")" _n _n ///
     as txt "   naturally dichotomous: " as inp "`vnat'" _n ///
     as txt "artificially dichotomous: " as inp "`vart'" _n _n ///
     as txt "nu(`vnat',`vart') = " as res `nu'
  return scalar nu = `nu'
  return scalar N = `N'
end

* ==============================================================================
/* A) Zusammenhangsanalysen nominaler Daten:

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

clear
input maennl fach anzahl
 0 1 18
 1 1 27
 0 2 36
 1 2 54
 0 3 18
 1 3 27
end
label define geschlecht 0 "weiblich" 1 "männlich"
label define fach 1 "Sprach- & Kult.wiss" ///
                  2 "Rechts-,Wi-,Soz.wiss." ///
				  3 "Mathem. & Nat.wiss."
label values maennl geschlecht
label values fach fach
label variable maennl "Geschlecht"
label variable fach "Studienfach"
expand anzahl
drop anzahl

tabulate maennl fach, exp /* Normalerweise würde man die potenziell abhängige
                             Variable Zeilenvariable spezifizieren, hier ist
							 es umgekehrt, um die Tabelle im Skript zu reprodu-
							 zieren. */

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

clear
input maennl fach anzahl
 0 1 28
 1 1 17
 0 2 31
 1 2 59
 0 3 13
 1 3 32
end
label define geschlecht 0 "weiblich" 1 "männlich"
label define fach 1 "Sprach- & Kult.wiss" ///
                  2 "Rechts-,Wi-,Soz.wiss." ///
				  3 "Mathem. & Nat.wiss."
label values maennl geschlecht
label values fach fach
label variable maennl "Geschlecht"
label variable fach "Studienfach"
expand anzahl
drop anzahl

tabulate maennl fach, exp chi2 V
cont_c maennl fach

* ------------------------------------------------------------------------------
*  3. Verurteilung und Geschlecht:

clear
input verurt maennl anzahl
 0 0 1869739
 1 0    1461
 0 1 1883749
 1 1   91251
end
label define geschlecht 0 "weiblich" 1 "männlich"
label define neinja 0 "nein" 1 "ja"
label values maennl geschlecht
label values verurt neinja
label variable maennl "Geschlecht"
label variable verurt "verurteilt"
expand anzahl
drop anzahl

tabulate verurt maennl
phi verurt maennl

* Prozent weiblicher und männlicher Jugendlicher verurteilt:
tabulate verurt maennl, nof co

* ==============================================================================
* B) "weitere Korrelationskoeffizienten":

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

clear
input test krank anzahl
 0 0 518
 1 0 182
 0 1 117
 1 1 183
end
label define test 0 "negativ" 1 "positiv"
label define neinja 0 "nein" 1 "ja"
label values test test
label values krank neinja
label variable test "Testergebnis"
label variable krank "Krankheit"
expand anzahl
drop anzahl

tabulate test krank
phi test krank
sens_spec test krank

* Basisrate Männer (% krank):
fre krank

* 2a) Odds-Ratio Krankheit bei positivem vs. negativem Test (Männer):
oddsratio test krank

* 3a) Yules-Y für Krankheit und Testergebnis (Männer):
yules_y test krank

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* 1b) Zusammenhang zwischen Test und Gesundheit (Frauen):

clear
input test krank anzahl
 0 0 666
 1 0 234
 0 1  39
 1 1  61
end
label define test 0 "negativ" 1 "positiv"
label define neinja 0 "nein" 1 "ja"
label values test test
label values krank neinja
label variable test "Testergebnis"
label variable krank "Krankheit"
expand anzahl
drop anzahl

tabulate test krank
phi test krank
sens_spec test krank

* Basisrate Frauen (% krank):
fre krank

* 2a) Odds-Ratio Krankheit bei positivem vs. negativem Test (Frauen):
oddsratio test krank

* 3a) Yules-Y für Krankheit und Testergebnis (Frauen):
yules_y test krank

* ------------------------------------------------------------------------------
/* 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): */

twoway (function normalden(invnormal(x))/sqrt(x*(1-x)), range(0 1) n(1000) lcolor(blue)), ///
       ytitle("Proportionale Reduktion von r") ytitle(, size(small)) yscale(range(0 1)) ///
       ylabel(0(0.2)1, grid glwidth(vthin) glcolor(gs12) glpattern(dash) labsize(small)) ///
       xtitle("Anteil der Population unter dem Teilungspunkt") xtitle(, size(small)) ///
       xlabel(, grid glwidth(vthin) glcolor(gs12) glpattern(dash) labsize(small)) ///
       title("Korrelation und Dichotomisierung kontinuierlicher Variablen", size(medsmall)) ///
       name("ReduktR",replace)

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

clear
set seed 24
set obs 1000000
gen double x = rnormal()
gen double y = rnormal()
mata
  X = J(1000000,2,.)
  X[.,1] = st_data(.,1)
  X[.,2] = st_data(.,2)
  sigma = J(2,2,1)
  sigma[1,2] = .6937
  sigma[2,1]=sigma[1,2]
  U = cholesky(sigma)
  X = X*U
  st_store(., "x", X[.,1])
  st_store(., "y", X[.,2])
  mata clear
end

* Korrelation der kontinuierlichen Variablen:
cor x y

* Dichotomisierung der Y-Variable beim Quantil 1/3:
local p100_3 = 100/3
qui centile y, c(`p100_3')
local cut_3 = r(c_1)
recode y (min/`r(c_1)'=0 "unteres Drittel") ///
         (`r(c_1)'/max=1 "obere 2 Drittel"), gen(yd)
fre yd

/* Korrelation der kontinuierlichen X- mit der dichotomisierten Y-Variablen
   (Produkt-Moment-Korrelation = punkt-biseriale Korrelation): */
cor x yd

* Bestimmung des Reduktionsfaktors:
local p = 1/3
local z = invnormal(`p')
local h = normalden(`z')
local red_fakt = `h'/sqrt(`p'*(1-`p'))

/* Schätzung der Korrelation von X mit der latenten kontinuierlichen Variablen Y
   (polyseriale Korrelation, die Korrelation der kontinuierlichen Variablen
   sollte reproduziert werden):
*/
di as txt "geschätzte Korrelation = " as res `r(rho)'/`red_fakt'

* Überprüfung, ob die Funktion polychor() das gleiche Ergebnis liefert:
polychoric x yd, verbose  /* bitte warten!! */

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

local dpx = invnorm(1/3)
local dpxt : display round(`dpx',.0001)
local dpy = normalden(invnorm(1/3))
local dpyt : display round(`dpy',.0001)
twoway (function normalden(x), range(-3.5 3.5) n(1000) lcolor(blue)) ///
       (scatteri 0 `dpx' `dpy' `dpx', ///
        recast(line) lcolor(black) lwidth(vthin) lpattern(dash)) ///
       (scatteri `dpy' `dpx' `dpy' -4, ///
        recast(line) lcolor(black) lwidth(vthin) lpattern(dash)), ///
       text(0.05 -1.64 "p= 0.3333", place(e) size(small)) ///
       text(0.05 -0.1 "q = (1-p) = 0.6667", place(e) size(small)) ///
       text(0 `dpx' "z = `dpxt'", place(e) size(small) col(red)) ///
       text(`dpy' -3.7 "h = `dpyt'", place(n) size(small) col(red)) ///
       ytitle("Dichte (h)") ytitle(, size(small)) ylabel(, labsize(small)) ///
       xtitle("z", size(small)) xlabel(, labsize(small)) ///
       title("Standardnormalverteilung", size(medsmall)) ///
       legend(off) name("StandNorm",replace)

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

use "${BaseURL}JugGew99.dta", clear

* Statistische Kennwerte für GE und SK:
moments2 GE SK

* Produkt-Moment-Korrelation GE und SK:
cor GE SK

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* Dichotomisierung der Selbstkontrolle (niedrig: <= 50; hoch: > 50):
recode SK (min/50=0 "SK <= 50") (50/max=1 "SK > 50"), gen(SK_d50)

* Anzahl der Fälle mit niedriger und hoher Selbstkontrolle:
fre SK_d50

* Streudiagramm von GE gegen SK mit Trennlinie für SK_d50:
twoway (scatter GE SK, msize(small)), ///
       yscale(range(0 100)) xscale(range(0 100)) ///
       ytitle(, size(small)) ylabel(, labsize(small)) ///
       xtitle(, size(small)) xlabel(, labsize(small)) ///
       xlabel(0(20)100) aspectratio(1) legend(off) ///
       xline(50, lwidth(thin) lpattern(solid) lcolor(blue)) ///
       name("GE_SKd50", replace)

* Produkt-Moment- (punkt-biseriale) Korrelation von GE und SK_d50:
cor GE SK_d50

* Polyseriale (biseriale) Korrelation von GE mit SK_d50:
polychoric GE SK_d50

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* Dichotomisierung der Gewalteinstellung (niedrig: < 50; hoch >= 50):
recode GE (min/50=0 "GE <= 50") (50/max=1 "GE > 50"), gen(GE_d50)

/* Kreuztabelle von dichotomisierter Selbstkontrolle und dichotomisierter
   Gewalteinstellung (ACHTUNG: die Reihenfolge von niedriger und hoher Gewalt-
   einstellung wird hier temporär so abgeändert, dass die Lage der Häufigkeiten
   in der Kreuztabelle mit der der Lage der Fälle in den Quadranten des Streu-
   diagramms korrespondieren): */
recode GE_d50 (1=0 "GE > 50") (0=1 "GE <= 50"), gen(GE_d50_tmp)
label variable GE_d50_tmp "RECODE OF GE (Gewalteinstellung)"
tabulate GE_d50_tmp SK_d50
drop GE_d50_tmp

* Streudiagramm mit Trennlinien für GE_d50 und SK_d50:
twoway (scatter GE SK, msize(small)), ///
       yscale(range(0 100)) xscale(range(0 100)) ///
       ytitle(, size(small)) ylabel(, labsize(small)) ///
       xtitle(, size(small)) xlabel(, labsize(small)) ///
       xlabel(0(20)100) aspectratio(1) legend(off) ///
       yline(50, lwidth(thin) lpattern(solid) lcolor(blue)) ///
       xline(50, lwidth(thin) lpattern(solid) lcolor(blue)) ///
       name("GEd50_SKd50", replace)

* Produkt-Moment-Korrelation (phi) von GE_d50 und SK_d50:
cor GE_d50 SK_d50

* Polychorische (tetrachorische) Korrelation von GE_d50 mit SK_d50:
polychoric GE_d50 SK_d50

* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
* Dichtomisierung der Selbstkontrolle für SK > 33.33:
local d33 = 100/3
recode SK (min/`d33'=0 "SK <= 33.33") (`d33'/max=1 "SK > 33.33"), gen(SK_d33)

/* Kreuztabelle von dichotomisierter Selbstkontrolle und dichotomisierter
   Gewalteinstellung (ACHTUNG: die Reihenfolge von niedriger und hoher Gewalt-
   einstellung wird hier temporär so abgeändert, dass die Lage der Häufigkeiten
   in der Kreuztabelle mit der der Lage der Fälle in den Quadranten des Streu-
   diagramms korrespondieren): */
recode GE_d50 (1=0 "GE > 50") (0=1 "GE <= 50"), gen(GE_d50_tmp)
label variable GE_d50_tmp "RECODE OF GE (Gewalteinstellung)"
tabulate GE_d50_tmp SK_d33
drop GE_d50_tmp

* Streudiagramm mit Trennlinien für GE_d50 und SK_d33:
local d33 = 100/3
twoway (scatter GE SK, msize(small)), ///
       yscale(range(0 100)) xscale(range(0 100)) ///
       ytitle(, size(small)) ylabel(, labsize(small)) ///
       xtitle(, size(small)) xlabel(, labsize(small)) ///
       xlabel(0(20)100) aspectratio(1) legend(off) ///
       yline(50, lwidth(thin) lpattern(dash) lcolor(red)) ///
       xline(`d33', lwidth(thin) lpattern(dash) lcolor(red)) ///
       name("GEd50_SKd33", replace)

* Produkt-Moment-Korrelation (phi) von GE_d50 und SK_d33:
cor GE_d50 SK_d33

* Polychorische (tetrachorische) Korrelation von GE_d50 mit SK_d33:
polychoric GE_d50 SK_d33

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

clear
input x
   5
   9
   9
   9
  12
  17
  17
  25
end

egen xrank = rank(x)
list x xrank, sep(0)

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

clear
input x y
 10 11
 14 17
  2  5
  7  1
 19 18
 15 16
  9 13
end
gen id = _n

* Rohwerte:
list id x y, sep(0)

* Produkt-Moment-Korrelation von x und y:
cor x y

* Tabelle der rangtransformierten Werte sortiert nach x:
egen xrank = rank(x)
egen yrank = rank(y)
sort x
list id xrank yrank, sep(0)

/* Produkt-Moment-Korrelation der rangtransformierten Werte (= Spearmans rho)
   (identisch mit Spearmans rho der untransformierten Werte): */
cor xrank yrank
spearman x y

* Kendalls tau:
ktau x y

* Polychorische Korrelation:
polychoric x y

* ------------------------------------------------------------------------------
* 10. Nu-Koeffizient:

* Annahme: Beispieltabelle basiert auf 100 Fällen:

clear
input maennl risiko anzahl
 0 0 15
 1 0 10
 0 1 25
 1 1 50
end
label define geschlecht 0 "weiblich" 1 "männlich"
label define neinja 0 "nein" 1 "ja"
label values maennl geschlecht
label values risiko neinja
label variable maennl "Geschlecht"
label variable risiko "Risiko (dichotomisiert)"
expand anzahl
drop anzahl

tabulate risiko maennl

* Produkt-Moment-Korrelation (= phi):
cor maennl risiko

/* Nu-Koeffizient (Achtung: zuerst muss die natürliche, dann die künstlich di-
   chotome Variable spezifiziert werden!): */
nucoef maennl 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 einlesen:
clear
input x y
  9  0
  5  1
  7  0
  8  0
  1  1
  5  1
  9  0
  4  1
  5  0
  9  0
  9  1
  7  1
end

/* Biseriale Rangkorrelation (= Somers D)
   (Achtung: y muss die dichotome Variable darstellen!): */
somersd y x

* ==============================================================================
log close