Ehescheidung

Öffnen Sie die Daten http://www.farys.org/daten/sams98.dta. Es handelt sich um einen Auszug aus dem Schweizer Arbeitsmarktsurvey 1998. Variablen sind:

Datenaufbereitung

Damit wir Ereignisdatenmodelle verwenden können, müssen die Daten etwas vorbereitet werden. Nutzen Sie hierfür die folgenden Code-Blöcke.

  1. Brauchen wir eine geeignete Skala in der alle Zeitpunkte gemessen sind. Z.B. Jahrhundertmonate (1 entspricht dem Januar 1900):
  2. 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
  3. Brauchen wir eine Zustandsvariable, ob die Person geschieden ist:
  4. Es macht Sinn, mehrere Situationen zu unterscheiden:
    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"
  5. Aus Heiratsmonat, Interviewmonat, Todesmonat und Scheidungsmonat lässt sich die Ehedauer für diverse Szenarien berechnen:
  6. 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
  7. Abschliessend bauen wir noch die Variable die das Ereignis (Scheidung) anzeigt:
  8. generate byte geschieden = ehestatus==2 if ehestatus<.

Übungsaufgaben

  1. Berechnen Sie die Überlebens- und Hazardfunktion mit Hilfe des Product-Limit-Schätzers (Befehl sts). Stellen Sie die Ergebnisse grafisch dar. Definieren Sie hierfür zunächst den Datensatz als Ereignisdaten (stset).
  2. 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

  3. Testen Sie folgende Hypothesen:

    Berechnen Sie dazu den Product-Limit-Schätzer für die Überlebens- und Hazardfunktionen und führen Sie Gruppenvergleiche durch. Kohabitation liegt vor, wenn das Datum der Gründung eines gemeinsamen Haushalts vor dem Datum der Eheschliessung liegt. Nehmen Sie weiterhin vereinfachend an, dass gemeinsame Kinder vorhanden sind, wenn das Geburtsdatum des ersten Kindes vor dem Scheidungsdatum liegt.
  4. Hypothese 1: Voreheliches Zusammenleben erhöht Scheidungsrisiko

    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.

    Hypothese 2: Kinder verringern Scheidungsrisiko

    . 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.

  5. Schätzen Sie für verschiedene Verteilungen parametrische Modelle (streg) die beide Kovariaten (Kohabitation und Kind) beinhalten.
  6. Erweitern Sie das Exponentialmodell nun zu einem “piecewise exponential model” (siehe stsplit).
  7. 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.

  8. Prüfen Sie erneut die beiden Hypothesen:
    1. Voreheliches Zusammenleben (Kohabitation) erhöht das Scheidungsrisiko (Selektions-effekt).
    2. Das Scheidungsrisiko wird durch das Vorhandensein gemeinsamer Kinder verringert.
    Verwenden Sie hierfür diesmal jedoch ein semi-parametrisches Modell (Cox-Regression; in Stata stcox).
  9. . 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.

  10. Bonus: Erzeugen Sie einige Kontrollvariablen und erweitern Sie das Modell um diese. Geben Sie auch eine kurze Interpretation der Befunde aus den Kontrollvariablen.
  11. Bilden der Kontrollvariablen

    . 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
    
    . 

    Cox-Regression mit Kontrollvariablen

    . 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.

  12. Prüfen Sie die Proportionalitätsannahme (estat phtest, detail). Betrachten Sie zudem für kritische Kovariate stcox diagnostic plots (z.B. stphplot).
  13. . 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.

  14. Um die zweite Hypothese adäquat zu prüfen, sollte der zeitveränderliche Charakter des Vorhandenseins gemeinsamer Kinder berücksichtigt werden. Führen Sie die Berechnungen (Cox-Regression ohne und ggf. mit Drittvariablen) mit einer zeitveränderlichen Variable aufgrund des Geburtsdatums des ersten Kindes durch (siehe stsplit). Unterscheiden Sie zudem, ob es sich um ein voreheliches Kind handelt oder nicht. Prüfen Sie auch hier die Proportionalitätsannahme.
  15. 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)

Modelle für diskrete Zeit

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:

  1. Erzeugen Sie eine Variable diplom (0 oder 1) ob in der betrachteten Periode ein Abschluss erlangt wurde.
  2. . 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)
  3. Erzeugen Sie nun einen Dummy dropout, der 1 wird, wenn die letzte Beobachtung einer Person vor 2014 war und kein Abschluss erreicht wurde.
  4. . * Abhängige variable: dropout (letzte Beobachtung vor 2014 und kein Diplom)
    
    . by id: gen byte dropout = (_n==_N) & (jahr<2014) & (diplom==0)
  5. 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           
    -----------------------------------------------------------
  6. Fassen Sie Studienjahre grösser gleich 6 zusammen in eine Kategorie. Schätzen Sie dann, inkl. Periodendummies, Geschlecht und Alter mehrere Modelle. Was lässt sich über die Schätzer (Geschlecht/Alter/Studienjahr) sagen?
  7. . * 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.