Öffnen Sie die Daten http://www.farys.org/daten/sams98.dta. Es handelt sich um einen Auszug aus dem Schweizer Arbeitsmarktsurvey 1998. Variablen sind:
Damit wir Ereignisdatenmodelle verwenden können, müssen die Daten etwas vorbereitet werden. Nutzen Sie hierfür die folgenden Code-Blöcke.
use http://www.farys.org/daten/sams98.dta, clear
generate geburt = a5j * 12 + a5m
generate heirat = h9j * 12 + cond(h9m>0, h9m, 7)
generate haushalt = h10j * 12 + cond(h10m>0, h10m, 7)
generate scheidung = h15j * 12 + cond(h15m>0, h15m, 7) if h15j>0
generate todpartner = h17j * 12 + cond(h17m>0, h17m, 7) if h17j>0
generate ersteskind = h19_1j * 12 + cond(h19_1m>0, h19_1m, 7) if h19_1j>0
generate interview = 98 * 12 + monat
recode h1 (1 2 = 1) (3 = 3) (4 = 2), into(ehestatus)
Zudem betrachten wir der Einfachheit halber nur die erste Ehe. Daher muss für Verheiratete und Geschiedene geprüft werden ob sie verheiratet/geschieden in erster oder späterer Ehe sind und ggf. die Zustandsvariable angepasst werden (z.B. ist “verheiratet in zweiter Ehe” eigentlich ein “geschieden in erster Ehe”)
replace ehestatus = 2 if h14==3
replace ehestatus = 3 if h14==4
label define ehestatus 1 "verheiratet" 2 "geschieden" 3 "verwitwet"
label values ehestatus ehestatus
label variable ehestatus "Status (erste) Ehe"
gen double ehedauer = interview + 1 - heirat if ehestatus==1
replace ehedauer = scheidung + 1 - heirat if ehestatus==2
replace ehedauer = todpartner + 1 - heirat if ehestatus==3
replace ehedauer = . if ehedauer<=0
Nicht in allen Fällen sind die Angaben konsistent, denn für die Ehedauer resultieren zum Teil negative Werte. Es wäre nun sinnvoll, diesen Fällen genauer nachzugehen. Um den Rahmen der Übung nicht zu sprengen, wollen wird diese Fälle hier aber einfach ausschliessen:
fre ehedauer if ehedauer<=0
replace ehedauer = . if ehedauer<=0
generate byte geschieden = ehestatus==2 if ehestatus<.
sts
). Stellen Sie die Ergebnisse grafisch dar. Definieren Sie hierfür zunächst den Datensatz als Ereignisdaten (stset
).
Die meisten Stata-Befehle für Ereignisdaten verlangen, dass die Daten mit dem stset
-Befehl zuerst entsprechend deklariert werden (ltable
wäre eine Ausnahme). Im vorliegenden Fall:
. stset ehedauer, failure(ehestatus==2) scale(12)
failure event: ehestatus == 2
obs. time interval: (0, ehedauer]
exit on or before: failure
t for analysis: time/12
------------------------------------------------------------------------------
2047 total obs.
4 event time missing (ehedauer>=.) PROBABLE ERROR
------------------------------------------------------------------------------
2043 obs. remaining, representing
347 failures in single record/single failure data
38756.58 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 48.66667
Die Option scale(12)
führt dazu, dass die Prozesszeit jeweils in Jahren berichtet wird anstatt in Monaten. Der Product-Limit-Schätzer von Überlebensfunktion und Hazardrate kann dann mit dem sts
-Befehl berechnet werden:
. sts graph, ci
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
. sts graph, hazard ci
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Testen Sie folgende Hypothesen:
Wir erzeugen einen Dummy ob vorher zusammen gelebt wurde (Monat erster gemeinsamer Haushalt kleiner Monat erste Ehe)
. generate byte kohab = haushalt < heirat if !missing(haushalt, heirat)
. fre kohab
kohab
-----------------------------------------------------------
| Freq. Percent Valid Cum.
--------------+--------------------------------------------
Valid 0 | 1128 55.11 55.11 55.11
1 | 919 44.89 44.89 100.00
Total | 2047 100.00 100.00
-----------------------------------------------------------
. sts graph, by(kohab) ci
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
. graph export surv_kohab.png, replace
(file surv_kohab.png written in PNG format)
. sts graph, by(kohab) hazard ci
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
. graph export hazard_kohab.png, replace
(file hazard_kohab.png written in PNG format)
Log-Rank und Wilcoxon-Test, ob die Überlebensfunktionen beider Gruppen gleich sind:
. sts test kohab, logrank
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Log-rank test for equality of survivor functions
| Events Events
kohab | observed expected
------+-------------------------
0 | 200 246.51
1 | 147 100.49
------+-------------------------
Total | 347 347.00
chi2(1) = 32.34
Pr>chi2 = 0.0000
. sts test kohab, wilcoxon
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Wilcoxon (Breslow) test for equality of survivor functions
| Events Events Sum of
kohab | observed expected ranks
------+--------------------------------------
0 | 200 246.51 -61631
1 | 147 100.49 61631
------+--------------------------------------
Total | 347 347.00 0
chi2(1) = 25.43
Pr>chi2 = 0.0000
. sts test kohab, tware
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Tarone-Ware test for equality of survivor functions
| Events Events Sum of
kohab | observed expected ranks
------+--------------------------------------
0 | 200 246.51 -1670.046
1 | 147 100.49 1670.046
------+--------------------------------------
Total | 347 347.00 0
chi2(1) = 28.85
Pr>chi2 = 0.0000
Die Unterschiede sind sehr deutlich. Ehen, denen eine Kohabitationsphase vorangegangen ist, werden schneller geschieden. Die Gründe könnten sein, dass dies die Ehen sind, bei denen sich die Partner von Anfang an nicht so sicher waren (deshalb die Probephase). Weiterhin könnte man argumentieren, dass man bei diesen Fällen die gesamte Dauer betrachten sollte, also ab Gründung eines gemeinsamen Haushalts und nicht erst ab Heirat.
. generate byte kind = (ersteskind - 6) < (heirat + ehedauer) if !missing(heirat, ehedauer)
(4 missing values generated)
. fre kind
kind
-----------------------------------------------------------
| Freq. Percent Valid Cum.
--------------+--------------------------------------------
Valid 0 | 301 14.70 14.73 14.73
1 | 1742 85.10 85.27 100.00
Total | 2043 99.80 100.00
Missing . | 4 0.20
Total | 2047 100.00
-----------------------------------------------------------
. sts graph, by(kind) ci
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
. graph export surv_kind.png, replace
(file surv_kind.png written in PNG format)
. sts graph, by(kind) hazard ci
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
. graph export hazard_kind.png, replace
(file hazard_kind.png written in PNG format)
Tests:
. sts test kind, logrank
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Log-rank test for equality of survivor functions
| Events Events
kind | observed expected
------+-------------------------
0 | 89 25.20
1 | 258 321.80
------+-------------------------
Total | 347 347.00
chi2(1) = 176.16
Pr>chi2 = 0.0000
. sts test kind, wilcoxon
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Wilcoxon (Breslow) test for equality of survivor functions
| Events Events Sum of
kind | observed expected ranks
------+--------------------------------------
0 | 89 25.20 110417
1 | 258 321.80 -110417
------+--------------------------------------
Total | 347 347.00 0
chi2(1) = 229.39
Pr>chi2 = 0.0000
. sts test kind, tware
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Tarone-Ware test for equality of survivor functions
| Events Events Sum of
kind | observed expected ranks
------+--------------------------------------
0 | 89 25.20 2646.9584
1 | 258 321.80 -2646.9584
------+--------------------------------------
Total | 347 347.00 0
chi2(1) = 208.15
Pr>chi2 = 0.0000
Bezüglich Hypothese 2 sind die Resultate sogar noch deutlicher. Kinderlose Ehepaare haben demnach ein viel höheres Scheidungsrisiko.
streg
) die beide Kovariaten (Kohabitation und Kind) beinhalten.
. streg kohab kind, distribution(exp)
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Iteration 0: log likelihood = -1208.5214
Iteration 1: log likelihood = -1147.2636
Iteration 2: log likelihood = -1142.8179
Iteration 3: log likelihood = -1142.811
Iteration 4: log likelihood = -1142.811
Exponential regression -- log relative-hazard form
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(2) = 131.42
Log likelihood = -1142.811 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.705317 .1878204 4.85 0.000 1.374217 2.116191
kind | .2533167 .0315697 -11.02 0.000 .1984186 .3234038
_cons | .0239558 .0029404 -30.40 0.000 .0188335 .0304713
------------------------------------------------------------------------------
Kohabitation erhöht das Scheidungsrisiko um 70,5%, oder anders: vergleicht man Per-sonen mit und ohne vorheriger Kohabitation, so ist die durchschnittliche Dauer bis zur Scheidung gut 40% kürzer:
. di 1-(1/1.705)
.41348974
Ein Kind senkt das Scheidungsrisiko um fast 75%. Ist ein Kind vorhanden, dauert es bis zur Scheidung fast 4 mal so lang als wenn kein Kind vorhanden wäre.
. streg kohab kind, distribution(wei)
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
Fitting constant-only model:
Iteration 0: log likelihood = -1208.5214
Iteration 1: log likelihood = -1208.1549
Iteration 2: log likelihood = -1208.1548
Fitting full model:
Iteration 0: log likelihood = -1208.1548
Iteration 1: log likelihood = -1141.4184
Iteration 2: log likelihood = -1135.2612
Iteration 3: log likelihood = -1135.2493
Iteration 4: log likelihood = -1135.2493
Weibull regression -- log relative-hazard form
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(2) = 145.81
Log likelihood = -1135.2493 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.955513 .2263935 5.79 0.000 1.558531 2.453612
kind | .2341921 .0295016 -11.52 0.000 .1829552 .2997779
_cons | .0129751 .0026957 -20.91 0.000 .008635 .0194965
-------------+----------------------------------------------------------------
/ln_p | .1840672 .0454319 4.05 0.000 .0950223 .2731122
-------------+----------------------------------------------------------------
p | 1.202097 .0546136 1.099683 1.314048
1/p | .8318799 .0377939 .7610074 .9093527
------------------------------------------------------------------------------
Die Ergebnisse sind ähnlich. Das Weibull-Modell nimmt allerdings keine konstante Ausfallrate (Scheidungsrate) an sondern besitzt einen zusätzlichen Parameter p der angibt, wie sich die Ausfallrate über die Zeit verändert. p=1 entspricht dem Exponentialmodell, p>1 deutet auf eine mit der Zeit steigende Ausfallrate hin. Alternativ könnte auch ein piecewise exponential Model geschätzt werden, um die restriktive Annahme für die Ausfallrate zu lockern.
stsplit
).
Zunächst muss auch id()
per stset
gesetzt werden (damit nachher klar ist, welche Teilstücke zur gleichen Episode gehören):
. stset ehedauer, failure(ehestatus==2) scale(12) id(id)
id: id
failure event: ehestatus == 2
obs. time interval: (ehedauer[_n-1], ehedauer]
exit on or before: failure
t for analysis: time/12
------------------------------------------------------------------------------
2047 total obs.
4 event time missing (ehedauer>=.) PROBABLE ERROR
------------------------------------------------------------------------------
2043 obs. remaining, representing
2043 subjects
347 failures in single failure-per-subject data
38756.58 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 48.66667
Dann teilen wir (recht arbiträr) die Zeitspanne z.B. in 5 Teile (je 10 Jahre) und verwenden diese Kategorien als Dummies:
. preserve
. stsplit tp, at(0(10)50)
(2885 observations (episodes) created)
. tab tp, gen(tp)
tp | Freq. Percent Cum.
------------+-----------------------------------
0 | 2,043 41.46 41.46
10 | 1,403 28.47 69.93
20 | 852 17.29 87.22
30 | 478 9.70 96.92
40 | 152 3.08 100.00
------------+-----------------------------------
Total | 4,928 100.00
. streg kohab kind tp1-tp4, dist(exp)
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
Iteration 0: log likelihood = -1208.5214
Iteration 1: log likelihood = -1139.0096
Iteration 2: log likelihood = -1129.2683
Iteration 3: log likelihood = -1128.8877
Iteration 4: log likelihood = -1128.8869
Iteration 5: log likelihood = -1128.8869
Exponential regression -- log relative-hazard form
No. of subjects = 2043 Number of obs = 4928
No. of failures = 347
Time at risk = 38756.58333
LR chi2(6) = 159.27
Log likelihood = -1128.8869 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.708687 .1965501 4.66 0.000 1.363794 2.1408
kind | .2435845 .0306 -11.24 0.000 .1904226 .3115881
tp1 | 3.37265 3.390327 1.21 0.227 .4702278 24.18991
tp2 | 5.289844 5.314779 1.66 0.097 .7382934 37.90153
tp3 | 4.10503 4.144934 1.40 0.162 .5673364 29.70244
tp4 | 1.483598 1.573642 0.37 0.710 .1855465 11.86259
_cons | .0064057 .0064404 -5.02 0.000 .0008928 .0459598
------------------------------------------------------------------------------
. restore
Die Koeffizienten für kind
und kohab
sind nach wie vor sehr robust.
stcox
).
. stcox kohab kind
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
Iteration 0: log likelihood = -2470.0547
Iteration 1: log likelihood = -2459.6095
Iteration 2: log likelihood = -2402.8313
Iteration 3: log likelihood = -2402.5817
Iteration 4: log likelihood = -2402.5816
Refining estimates:
Iteration 0: log likelihood = -2402.5816
Cox regression -- Breslow method for ties
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(2) = 134.95
Log likelihood = -2402.5816 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.740051 .1998736 4.82 0.000 1.389273 2.179398
kind | .2331076 .0292499 -11.61 0.000 .1822843 .2981011
------------------------------------------------------------------------------
Die Ergebnisse sind recht ähnlich zu den parametrischen Modellen.
. generate byte koho4660 = inrange(h9j, 46, 60)
. generate byte koho6170 = inrange(h9j, 61, 70)
. generate byte koho7180 = inrange(h9j, 71, 80)
. generate byte koho8190 = inrange(h9j, 81, 90)
. generate byte koho9198 = inrange(h9j, 91, 98)
. generate heiratsalter = trunc((heirat - geburt)/12)
. generate byte r_katholisch = i4==1
. generate byte r_reformiert = i4==2
. generate byte r_andere = inlist(i4, 3, 4, 5)
. generate byte r_keine = i4==6
. summarize koho4660 koho6170 koho7180 koho8190 koho9198 heiratsalter r*
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
koho4660 | 2047 .1431363 .3502974 0 1
koho6170 | 2047 .2344895 .4237829 0 1
koho7180 | 2047 .198339 .3988463 0 1
koho8190 | 2047 .2510992 .4337514 0 1
koho9198 | 2047 .172936 .3782843 0 1
-------------+--------------------------------------------------------
heiratsalter | 2047 25.47386 4.951591 14 63
r_katholisch | 2047 .4626282 .4987232 0 1
r_reformiert | 2047 .3873962 .4872744 0 1
r_andere | 2047 .0517831 .221643 0 1
r_keine | 2047 .0981925 .2976474 0 1
.
. stcox kohab kind koho* heiratsalter r*
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
note: koho9198 omitted because of collinearity
note: r_keine omitted because of collinearity
Iteration 0: log likelihood = -2470.0547
Iteration 1: log likelihood = -2422.2824
Iteration 2: log likelihood = -2352.7027
Iteration 3: log likelihood = -2352.3646
Iteration 4: log likelihood = -2352.3645
Refining estimates:
Iteration 0: log likelihood = -2352.3645
Cox regression -- Breslow method for ties
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(10) = 235.38
Log likelihood = -2352.3645 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.631181 .217383 3.67 0.000 1.256218 2.118066
kind | .1981133 .0255421 -12.56 0.000 .1538761 .2550682
koho4660 | .4220232 .1465991 -2.48 0.013 .2136253 .8337196
koho6170 | .9555926 .3008426 -0.14 0.885 .5155786 1.771131
koho7180 | .9794043 .2962249 -0.07 0.945 .5413929 1.771787
koho8190 | 1.039867 .3100269 0.13 0.896 .5796955 1.865331
koho9198 | 1 (omitted)
heiratsalter | .9141727 .0129838 -6.32 0.000 .8890759 .9399779
r_katholisch | .3933232 .0610333 -6.01 0.000 .2901787 .5331304
r_reformiert | .5403167 .0818558 -4.06 0.000 .4015085 .7271131
r_andere | .4770719 .1395917 -2.53 0.011 .2688571 .8465374
r_keine | 1 (omitted)
------------------------------------------------------------------------------
Der Effekt des Zusammenlebens wird etwas abgeschwächt, was wohl mit den Kohorten zusammenhängt. Der Effekt der Kinder zeigt sich noch etwas stärker. Die Effekte der Kontrollvariablen machen mehr oder weniger Sinn: In der Referenzkohorte (1946-1960) ist die Scheidungsrate am geringsten (dass der Unterschied zur Heiratskohorte 1991-1998 weniger deutlich signifikant ist, hat u.a. damit zu tun, dass es in dieser Kohorte sehr viele zensierte Fälle gibt). Das Heiratsalter wirkt sich negativ auf das Scheidungsrisiko aus. Katholiken lassen sich am wenigsten scheiden, Reformierte und andere Religion lassen sich ebenfalls weniger scheiden als Konfessionslose.
estat phtest, detail
). Betrachten Sie zudem für kritische Kovariate stcox
diagnostic plots (z.B. stphplot
).
. stcox kohab kind koho* heiratsalter r*
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
note: koho9198 omitted because of collinearity
note: r_keine omitted because of collinearity
Iteration 0: log likelihood = -2470.0547
Iteration 1: log likelihood = -2422.2824
Iteration 2: log likelihood = -2352.7027
Iteration 3: log likelihood = -2352.3646
Iteration 4: log likelihood = -2352.3645
Refining estimates:
Iteration 0: log likelihood = -2352.3645
Cox regression -- Breslow method for ties
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(10) = 235.38
Log likelihood = -2352.3645 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.631181 .217383 3.67 0.000 1.256218 2.118066
kind | .1981133 .0255421 -12.56 0.000 .1538761 .2550682
koho4660 | .4220232 .1465991 -2.48 0.013 .2136253 .8337196
koho6170 | .9555926 .3008426 -0.14 0.885 .5155786 1.771131
koho7180 | .9794043 .2962249 -0.07 0.945 .5413929 1.771787
koho8190 | 1.039867 .3100269 0.13 0.896 .5796955 1.865331
koho9198 | 1 (omitted)
heiratsalter | .9141727 .0129838 -6.32 0.000 .8890759 .9399779
r_katholisch | .3933232 .0610333 -6.01 0.000 .2901787 .5331304
r_reformiert | .5403167 .0818558 -4.06 0.000 .4015085 .7271131
r_andere | .4770719 .1395917 -2.53 0.011 .2688571 .8465374
r_keine | 1 (omitted)
------------------------------------------------------------------------------
. estat phtest, log detail
Test of proportional-hazards assumption
Time: Log(t)
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
kohab | 0.09957 3.93 1 0.0476
kind | 0.30259 30.25 1 0.0000
koho4660 | -0.02967 0.29 1 0.5924
koho6170 | -0.03251 0.35 1 0.5537
koho7180 | -0.06554 1.46 1 0.2270
koho8190 | -0.04671 0.75 1 0.3854
o.koho9198 | . . 1 .
heiratsalter| 0.00075 0.00 1 0.9878
r_katholisch| 0.05537 1.10 1 0.2946
r_reformiert| 0.05265 0.96 1 0.3277
r_andere | 0.00069 0.00 1 0.9895
o.r_keine | . . 1 .
------------+---------------------------------------------------
global test | 39.19 10 0.0000
----------------------------------------------------------------
. stphplot, by(kind)
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
Weiterhin kann die Proportionalitätsannahme auch getestet werden, indem Interaktionen mit der Prozesszeit explizit ins Modell aufgenommen werden. Im Cox-Modell kann man dazu die Optionen tvc()
und texp()
verwenden:
. stcox kohab kind koho* heiratsalter r*, tvc(kind) texp(ln(_t))
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
note: koho9198 omitted because of collinearity
note: r_keine omitted because of collinearity
Iteration 0: log likelihood = -2470.0547
Iteration 1: log likelihood = -2388.8119
Iteration 2: log likelihood = -2361.4186
Iteration 3: log likelihood = -2336.4658
Iteration 4: log likelihood = -2336.1582
Iteration 5: log likelihood = -2336.1579
Refining estimates:
Iteration 0: log likelihood = -2336.1579
Cox regression -- Breslow method for ties
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(11) = 267.79
Log likelihood = -2336.1579 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
main |
kohab | 1.57455 .2100683 3.40 0.001 1.212254 2.045122
kind | .0283504 .0113925 -8.87 0.000 .0128976 .0623176
koho4660 | .5279324 .1853749 -1.82 0.069 .2652728 1.050664
koho6170 | 1.180672 .3759522 0.52 0.602 .6325397 2.203793
koho7180 | 1.25611 .3846006 0.74 0.456 .6892926 2.289032
koho8190 | 1.270023 .3819405 0.79 0.427 .7044102 2.289801
koho9198 | 1 (omitted)
heiratsalter | .9159835 .0128134 -6.27 0.000 .8912108 .9414448
r_katholisch | .3996302 .0620866 -5.90 0.000 .2947241 .5418774
r_reformiert | .5523031 .0838282 -3.91 0.000 .4101878 .7436562
r_andere | .4903119 .1433904 -2.44 0.015 .2764017 .8697694
r_keine | 1 (omitted)
-------------+----------------------------------------------------------------
tvc |
kind | 2.620743 .4999735 5.05 0.000 1.803168 3.809016
------------------------------------------------------------------------------
Note: variables in tvc equation interacted with ln(_t)
Ein signifikanter Interaktionsterm weist auf eine Verletzung der Proportionalitätsannahme hin. Die Aufnahme von Zeit-Interaktionen ist letztlich auch die Methode, mit der man die Proportionalitätsannahme umgehen kann, d.h. man modelliert ganz einfach die Nicht-Proportionalität. Die Effekte sind dann allerdings etwas schwieriger zu interpretieren. Im vorliegenden Fall ist der Effekt zu Beginn stark negativ und schwächt sich dann über die Zeit ab. Zum Zeitpunkt t = 1 beträgt der Effekt 0.028 ⋅ 2.62ln(1) = 0.028 (Verringerung der Hazardrate um 97%). Nach zehn Ehejahren ist der Effekt 0.028 ⋅ 2.62ln(10) = 0.257 (Verringerung der Hazardrate um 74%). Zumindest bei kategorialen Variablen mit relativ wenigen Ausprägungen hat man auch die Möglichkeit, ein stratifiziertes Modell zu schätzen. In einem stratifizierten Modell wird für jede Kategorie (jedes Stratum) eine eigene Baseline-Hazardfunktion verwendet, wodurch sich die Proportionalitätsannahme erübrigt (allerdings nur bezüglich der Variable, nach der stratifiziert wird). Beispiel:
. stcox kohab koho* heiratsalter r*, strata(kind)
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
note: koho9198 omitted because of collinearity
note: r_keine omitted because of collinearity
Iteration 0: log likelihood = -2238.6625
Iteration 1: log likelihood = -2182.9919
Iteration 2: log likelihood = -2178.9808
Iteration 3: log likelihood = -2178.9597
Iteration 4: log likelihood = -2178.9597
Refining estimates:
Iteration 0: log likelihood = -2178.9597
Stratified Cox regr. -- Breslow method for ties
No. of subjects = 2043 Number of obs = 2043
No. of failures = 347
Time at risk = 38756.58333
LR chi2(9) = 119.41
Log likelihood = -2178.9597 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
kohab | 1.567688 .2091894 3.37 0.001 1.206915 2.036302
koho4660 | .5185118 .1820508 -1.87 0.061 .2605551 1.031853
koho6170 | 1.15775 .3686864 0.46 0.646 .6202244 2.161128
koho7180 | 1.241498 .3800053 0.71 0.480 .6814045 2.261971
koho8190 | 1.252356 .3767576 0.75 0.454 .6944692 2.258408
koho9198 | 1 (omitted)
heiratsalter | .9153559 .0127771 -6.34 0.000 .8906527 .9407442
r_katholisch | .3992674 .0620478 -5.91 0.000 .2944311 .5414322
r_reformiert | .5528388 .0839698 -3.90 0.000 .4104978 .7445368
r_andere | .49246 .1440684 -2.42 0.015 .2775577 .8737531
r_keine | 1 (omitted)
------------------------------------------------------------------------------
Stratified by kind
Der Nachteil der Stratifizierung ist, dass man für das betreffende Merkmal keinen expliziten Effekt erhält. Der Effekt ist nur implizit als Unterschied zwischen den Baseline-Raten im Modell enthalten. Zu bemerken ist bezüglich des vorliegenden Beispiels, dass die Effekte der anderen Kovariaten durch die Verletzung der Proportionalitätsannahme durch die Variable kind
kaum tangiert werden. Die Effekte sind im stratifizierten und im nicht stratifizierten Modell etwa gleich.
stsplit
). Unterscheiden Sie zudem, ob es sich um ein voreheliches Kind handelt oder nicht. Prüfen Sie auch hier die Proportionalitätsannahme.
Der zeitveränderliche Charakter des Vorhandenseins gemeinsamer Kinder kann mit Hilfe eines Episodensplits berücksichtigt werden. Die Ehe-Episoden werden dabei am Geburtsdatum des ersten Kindes in zwei Teile zerlegt. Im vorliegenden Beispiel werden vom Geburtsdatum noch 6 Monate abgezogen, da wohl nicht erst die Geburt das relevante Ereignis ist, sondern bereits die (gesicherte) Schwangerschaft. Zuerst einige Hilfsvariablen:
. generate t_kind = (ersteskind - 6) - heirat
(275 missing values generated)
. replace t_kind = 9999 if t_kind>=.
(275 real changes made)
. generate byte kindvorehe = t_kind<=0
. summarize t_kind kindvorehe
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
t_kind | 2047 1355.18 3406.179 -246 9999
kindvorehe | 2047 .2535418 .4351446 0 1
. stset ehedauer, failure(ehestatus==2) scale(12) id(id)
id: id
failure event: ehestatus == 2
obs. time interval: (ehedauer[_n-1], ehedauer]
exit on or before: failure
t for analysis: time/12
------------------------------------------------------------------------------
2047 total obs.
4 event time missing (ehedauer>=.) PROBABLE ERROR
------------------------------------------------------------------------------
2043 obs. remaining, representing
2043 subjects
347 failures in single failure-per-subject data
38756.58 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 48.66667
. stsplit tvc_kind, at(0) after(t_kind)
(1224 observations (episodes) created)
. tab tvc_kind
tvc_kind | Freq. Percent Cum.
------------+-----------------------------------
-1 | 1,525 46.68 46.68
0 | 1,742 53.32 100.00
------------+-----------------------------------
Total | 3,267 100.00
. replace tvc_kind = tvc_kind + 1
(3267 real changes made)
. fre tvc_kind
tvc_kind
-----------------------------------------------------------
| Freq. Percent Valid Cum.
--------------+--------------------------------------------
Valid 0 | 1525 46.62 46.68 46.68
1 | 1742 53.26 53.32 100.00
Total | 3267 99.88 100.00
Missing . | 4 0.12
Total | 3271 100.00
-----------------------------------------------------------
Die zeitveränderliche Kovariate tvc_kind
kann nun wie jede andere Variable ins Modell aufgenommen werden:
. stcox tvc_kind kindvorehe kohab koho* heiratsalter r*
failure _d: ehestatus == 2
analysis time _t: ehedauer/12
id: id
note: koho9198 omitted because of collinearity
note: r_keine omitted because of collinearity
Iteration 0: log likelihood = -2470.0547
Iteration 1: log likelihood = -2391.2127
Iteration 2: log likelihood = -2380.3406
Iteration 3: log likelihood = -2380.2702
Iteration 4: log likelihood = -2380.2702
Refining estimates:
Iteration 0: log likelihood = -2380.2702
Cox regression -- Breslow method for ties
No. of subjects = 2043 Number of obs = 3267
No. of failures = 347
Time at risk = 38756.58333
LR chi2(11) = 179.57
Log likelihood = -2380.2702 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
tvc_kind | .2899442 .0414785 -8.65 0.000 .2190505 .3837819
kindvorehe | 1.279262 .1683038 1.87 0.061 .9884911 1.655566
kohab | 1.604148 .2141291 3.54 0.000 1.234873 2.083852
koho4660 | .3387934 .117113 -3.13 0.002 .1720657 .6670763
koho6170 | .7389865 .2309183 -0.97 0.333 .4005469 1.363389
koho7180 | .7842027 .2355496 -0.81 0.418 .4352657 1.41287
koho8190 | .8478646 .2515764 -0.56 0.578 .4739799 1.516677
koho9198 | 1 (omitted)
heiratsalter | .9250254 .0131537 -5.48 0.000 .8996007 .9511687
r_katholisch | .3847717 .059753 -6.15 0.000 .2838024 .5216631
r_reformiert | .5228438 .0789956 -4.29 0.000 .3888349 .7030377
r_andere | .4699521 .1378435 -2.57 0.010 .2644748 .8350699
r_keine | 1 (omitted)
------------------------------------------------------------------------------
Auch in dieser Modellierung zeigt sich ein stark negativer Effekt der Geburt des ersten Kindes auf das Scheidungsrisiko (wobei es aber natürlich auch sein kann, dass ein Kind gerade darum zur Welt kommt, weil die Ehe stabil ist). Wenn ein Kind schon vor der Eheschliessung geboren wurde, ist das Scheidungsrisiko hingegen leicht erhöht. Für eine Verletzung der Proportionalitätsannahme finden sich in dieser Modellierung übrigens nur noch schwache Hinweise (man beachte, dass der globale Test nicht signifikant ist):
. qui stcox tvc_kind kindvorehe kohab kind koho* heiratsalter r*
. estat phtest, log detail
Test of proportional-hazards assumption
Time: Log(t)
----------------------------------------------------------------
| rho chi2 df Prob>chi2
------------+---------------------------------------------------
tvc_kind | . . 1 .
kindvorehe | -0.02903 0.31 1 0.5774
kohab | 0.10109 4.06 1 0.0440
kind | 0.25526 21.20 1 0.0000
koho4660 | -0.03047 0.30 1 0.5816
koho6170 | -0.03300 0.36 1 0.5464
koho7180 | -0.06780 1.58 1 0.2094
koho8190 | -0.04850 0.82 1 0.3665
o.koho9198 | . . 1 .
heiratsalter| -0.00185 0.00 1 0.9702
r_katholisch| 0.05053 0.93 1 0.3347
r_reformiert| 0.04959 0.85 1 0.3552
r_andere | -0.00230 0.00 1 0.9649
o.r_keine | . . 1 .
------------+---------------------------------------------------
global test | 31.21 11 0.0010
----------------------------------------------------------------
.
Dies macht Sinn, denn die Nichtproportionalität im vorherigen Modell war ein Artefakt der Definition der Variable kind
. Die Scheidungsraten für Fälle mit kind==1
werden in diesem Modell zu Beginn der Prozesszeit unterschätzt, weil eine Ehe mindestens eine gewisse Zeit überdauern muss (nämlich bis zur Geburt des ersten Kindes), damit die Variable kind
überhaupt den Wert eins annehmen kann. Nur mit einer zeitveränderlichen Variable ist also eine adäquate Modellierung möglich. Das Episodensplitting kann übrigens rückgängig gemacht werden, indem man die zeitveränderlichen Variablen eliminiert und dann stjoin
anwendet:
. drop tvc_kind
. stjoin
(1224 obs. eliminated)
Laden Sie die Daten http://www.farys.org/daten/dropout.dta. Es handelt sich um Daten im long-Format (mehrere Zeitpunkte pro id
) Die Variable eventtype
enthält die Info, ob die jeweilige Person in der betrachteten Periode immatrikuliert war oder abgeschlossen hat. Personen die 2014 weder immatrikuliert sind noch ein Diplom (vorher) erworben haben sind ein Dropout (Abbruch des Studiums).
Es müssen drei Fälle unterschieden werden:
diplom
(0 oder 1) ob in der betrachteten Periode ein Abschluss erlangt wurde.
. use http://www.farys.org/daten/dropout.dta, clear
.
. * Abschluss
. gen byte diplom = eventtyp==2
.
Es gibt noch ein paar inkonsistente Fälle, die noch Beobachtungen nach ihrem Abschluss haben.
. * Fälle löschen, die nach Diplom weitere Einträge haben
. * wichtig: Sortierung innerhalb der Person nach Beobachtungsdatum
. sort id e_date
. capture noisily by id: assert(diplom==0) if _n<_N
6 contradictions in 1838 observations
assertion is false
. by id: drop if diplom[_n-1]==1
(6 observations deleted)
dropout
, der 1 wird, wenn die letzte Beobachtung einer Person vor 2014 war und kein Abschluss erreicht wurde.
. * Abhängige variable: dropout (letzte Beobachtung vor 2014 und kein Diplom)
. by id: gen byte dropout = (_n==_N) & (jahr<2014) & (diplom==0)
Die Beobachtungen der einzelnen Personen beziehen sich auf unterschiedliche Zeitspannen. Zur Ereignisdatenanalyse müssen alle Personen am Startdatum ausgerichtet werden. Erzeugen Sie dafür eine neue Variable studienjahr
, die pro id abbildet in welchem Studienjahr sich die Person zum jeweiligen Zeitpunkt befunden hat.
. * Studienjahr
. by id: gen studienjahr = jahr - jahr[1] + 1
. fre studienjahr
studienjahr
-----------------------------------------------------------
| Freq. Percent Valid Cum.
--------------+--------------------------------------------
Valid 1 | 775 29.73 29.73 29.73
2 | 605 23.21 23.21 52.93
3 | 490 18.80 18.80 71.73
4 | 417 16.00 16.00 87.73
5 | 217 8.32 8.32 96.05
6 | 69 2.65 2.65 98.70
7 | 22 0.84 0.84 99.54
8 | 7 0.27 0.27 99.81
9 | 4 0.15 0.15 99.96
10 | 1 0.04 0.04 100.00
Total | 2607 100.00 100.00
-----------------------------------------------------------
logit
)
cloglog
)
reg
)
. * discrete-time logistic hazard model
. * 6 oder mehr zusammenfassen
. gen studienjahr2 = min(studienjahr, 6)
. logit dropout frau alter i.studienjahr2, or
Iteration 0: log likelihood = -452.72111
Iteration 1: log likelihood = -428.55191
Iteration 2: log likelihood = -422.63342
Iteration 3: log likelihood = -422.41983
Iteration 4: log likelihood = -422.41855
Iteration 5: log likelihood = -422.41855
Logistic regression Number of obs = 2607
LR chi2(7) = 60.61
Prob > chi2 = 0.0000
Log likelihood = -422.41855 Pseudo R2 = 0.0669
------------------------------------------------------------------------------
dropout | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
frau | 1.08568 .2205899 0.40 0.686 .7290442 1.616777
alter | .996716 .0296837 -0.11 0.912 .9402023 1.056627
|
studienjahr2 |
2 | .4803576 .1173349 -3.00 0.003 .2976085 .7753255
3 | .3288643 .1010742 -3.62 0.000 .1800553 .6006586
4 | .0540521 .0392591 -4.02 0.000 .0130188 .2244173
5 | .1042663 .0762847 -3.09 0.002 .0248529 .4374328
6 | .2239986 .1685786 -1.99 0.047 .0512438 .9791501
|
_cons | .0917428 .059248 -3.70 0.000 .0258739 .3252986
------------------------------------------------------------------------------
.
. * discrete-time proportional hazards model
. cloglog dropout frau alter i.studienjahr2, eform
Iteration 0: log likelihood = -422.48564
Iteration 1: log likelihood = -422.41072
Iteration 2: log likelihood = -422.41069
Complementary log-log regression Number of obs = 2607
Zero outcomes = 2498
Nonzero outcomes = 109
LR chi2(7) = 60.62
Log likelihood = -422.41069 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
dropout | exp(b) Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
frau | 1.086742 .2139647 0.42 0.673 .7388159 1.598514
alter | .9967213 .0286665 -0.11 0.909 .9420904 1.05452
|
studienjahr2 |
2 | .4909819 .1166969 -2.99 0.003 .3081424 .7823113
3 | .3383908 .1019241 -3.60 0.000 .1875144 .6106642
4 | .0563091 .0407302 -3.98 0.000 .013642 .2324234
5 | .1083594 .0787663 -3.06 0.002 .0260693 .4504057
6 | .2315006 .171935 -1.97 0.049 .0539966 .9925166
|
_cons | .0877971 .0547216 -3.90 0.000 .0258791 .2978598
------------------------------------------------------------------------------
.
. * discrete-time LPM (zeigt direkt die Abbruchquoten nach Jahr)
. regress dropout frau alter ibn.studienjahr2, nocons robust
Linear regression Number of obs = 2607
F( 8, 2599) = 14.50
Prob > F = 0.0000
R-squared = 0.0623
Root MSE = .19831
------------------------------------------------------------------------------
| Robust
dropout | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
frau | .0032372 .0079402 0.41 0.684 -.0123325 .0188068
alter | -.0001301 .0009824 -0.13 0.895 -.0020564 .0017962
|
studienjahr2 |
1 | .0833708 .0226169 3.69 0.000 .0390218 .1277198
2 | .0422387 .0230961 1.83 0.068 -.0030499 .0875274
3 | .0296206 .0240188 1.23 0.218 -.0174773 .0767185
4 | .0059685 .0236551 0.25 0.801 -.0404161 .0523532
5 | .010361 .0254953 0.41 0.684 -.0396322 .0603542
6 | .0209215 .0254045 0.82 0.410 -.0288936 .0707366
------------------------------------------------------------------------------
Referenzkategorie ist der Abbruch direkt nach dem ersten Jahr. D.h. die Chance abzubrechen sinkt erstmal wenn das erste Jahr “überstanden” ist. Die Chance steigt jedoch für ältere Semester tendenziell wieder an (dünne Datenlage). Zwischen Männern und Frauen gibt es keinen Unterschied in der Abbruchwahrscheinlichkeit. Auch das Alter hat keinen Einfluss. Mit dem LPM lassen sich direkt die Abbruchwahrscheinlichkeiten je Studienjahr darstellen.