    * Survival 11 years test of standard errors.sps.
    * Produce file of survival times for teeth filled in period.
    * Adjusted for time dependency.
    * Test standard errors by using random censoring as alternative to weighted adjustment.
    * Written by PSKL on 24/09/02.

    * Work from file of observed restoration lives.

    get file='D:\Longitudinal Data\Observed Restoration Lives 1991 to 2001.sav'.

    * Restrict to restorations placed before the end of the observation period 31/12/2001.

    select if (doplacn<yrmoda(2001,12,31)-yrmoda(1899,12,31)).

    * Check whether time is in acceptable range, and set reint flag.

    compute reint=0.
    if not (missing(time)) reint=1.
    if (time le 0) reint=0.
    if (doplacn+time>yrmoda(2001,12,31)-yrmoda(1899,12,31)) reint=0.

    execute.

    * Create weights to adjust for time dependent censoring.

    * Double up the cases with no re-intervention and weight by probability of re-intervention.
    compute caseid=$casenum.
    execute.

    * Look up probability of reattendance after last visit.
    sort cases by interval.
    match files file=*
     /table='D:\Longitudinal Data\reattendance probability.sav'
     /by interval.
    compute recid=1.
    sort cases by caseid.
    save outfile='d:\temp1.sav'.

    select if (reint=0).
    compute recid=2.
    sort cases by caseid.
    save outfile='d:\temp2.sav'.

    add files file='d:\temp1.sav'
     /file='d:\temp2.sav'
     /by caseid .

    Compute weight = 1000.
    Do if (reint=0).
    if (recid=1) weight = 1000-trunc(1000*prreatt+0.5).
    if (recid=2) weight = trunc(1000*prreatt+0.5).
    end if.

    save outfile='d:\temp3.sav'.

    get file='d:\temp3.sav'.

    ***************************************************.
    * Compare and calculate survival function and statistics, and standard errors.

    * Assume censored cases are lost at one year (365 days) after last visit.

    * recalculate time if reint=0.
    * Impute effective time of censoring.
    Do if (reint=0).
    If (recid=1) time=min(lastvn+365,yrmoda(2001,12,31)-yrmoda(1899,12,31))- doplacn.
    if (recid=2) time=yrmoda(2001,12,31)-yrmoda(1899,12,31)- doplacn.
    end if.

    save outfile= 'd:\temp4.sav'.
    get file='d:\temp4.sav'.

    *  Test 1 Adjusted KM, but increased weights give incorrect confidence intervals.

    weight by weight.

    KM
      time  /STATUS=reint(1)
        /PRINT MEAN
        /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95)
        /Save survival (survadj).

# Kaplan-Meier

Notes
Output Created | 19-NOV-2002 12:33:21
Comments |
Input | Data | d:\temp4.sav
Filter | <none>
Weight | WEIGHT
Split File | <none>
N of Rows in Working Data File | 843599
Syntax | KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95) /Save survival (survadj).
Resources | Elapsed Time | 0:01:36.44

     Survival Analysis for TIME

              Survival Time    Standard Error   95% Confidence Interval

     Mean:      2537.38                 .09     (  2537.21,   2537.55 )
     (Limited to   4016.0 )
     Median:    3003.00                 .45     (  3002.12,   3003.88 )

                                                                Percentiles

                         45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

     Value             3761.00      3003.00      2430.00      1964.00      1562.00      1230.00       949.00       709.00       500.00
     Standard Error        .            .45          .31          .24          .19          .14          .11          .09          .06

                                                                Percentiles

                         90.00        95.00

     Value              321.00       168.00
     Standard Error        .05          .

    >Warning # 3211
    >On at least one case, the value of the weight variable was zero, negative,
    >or missing.  Such cases are invisible to statistical procedures and graphs
    >which need positively weighted cases, but remain on the file and are
    >processed by non-statistical facilities such as LIST and SAVE.

    * Test 2 Unadjusted KM, for confidence intervals.

    weight off.

    * Censor all non-returners.
    select if (recid=1).

    KM
      time  /STATUS=reint(1)
       /PRINT MEAN
       /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95)
       /Save survival (survunad).

# Kaplan-Meier

Notes
Output Created | 19-NOV-2002 12:35:01
Comments |
Input | Data | d:\temp4.sav
Filter | <none>
Weight | <none>
Split File | <none>
N of Rows in Working Data File | 503965
Syntax | KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95) /Save survival (survunad).
Resources | Elapsed Time | 0:01:30.12

     Survival Analysis for TIME

              Survival Time    Standard Error   95% Confidence Interval

     Mean:      2483.62                2.83     (  2478.07,   2489.17 )
     (Limited to   4016.0 )
     Median:    2766.00               11.51     (  2743.43,   2788.57 )

                                                                Percentiles

                         45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

     Value             3362.00      2766.00      2276.00      1857.00      1495.00      1188.00       924.00       696.00       496.00
     Standard Error        .          11.51         8.73         6.88         5.45         4.29         3.36         2.62         1.99

                                                                Percentiles

                         90.00        95.00

     Value              321.00       168.00
     Standard Error       1.40          .

    * Test 3 Censoring using randomisation instead of weighting.

    * Censor the proportion prreatt of all non-returners at the end of the observation period (instead of earlier).

    Set seed = 10102002.
    compute random=uniform(1).
    if (reint=0 and random<prreatt) time=yrmoda(2001,12,31)-yrmoda(1899,12,31)- doplacn.

    KM
      time  /STATUS=reint(1)
       /PRINT MEAN
       /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95)
       /Save survival (survrand).

# Kaplan-Meier

Notes
Output Created | 19-NOV-2002 12:36:32
Comments |
Input | Data | d:\temp4.sav
Filter | <none>
Weight | <none>
Split File | <none>
N of Rows in Working Data File | 503965
Syntax | KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95) /Save survival (survrand).
Resources | Elapsed Time | 0:02:31.14

     Survival Analysis for TIME

              Survival Time    Standard Error   95% Confidence Interval

     Mean:      2537.29                2.77     (  2531.87,   2542.71 )
     (Limited to   4016.0 )
     Median:    3001.00               14.27     (  2973.02,   3028.98 )

                                                                Percentiles

                         45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

     Value             3761.00      3001.00      2429.00      1964.00      1562.00      1230.00       949.00       709.00       500.00
     Standard Error        .          14.27         9.89         7.57         5.91         4.55         3.53         2.72         2.04

                                                                Percentiles

                         90.00        95.00

     Value              321.00       168.00
     Standard Error       1.42          .

    save outfile='D:\temp5.sav'.

    * Rerun to test robustness.

    get file='d:\temp4.sav'.
    select if (recid=1).
    compute random=uniform(1).
    if (reint=0 and random<prreatt) time=yrmoda(2001,12,31)-yrmoda(1899,12,31)- doplacn.

    KM
      time  /STATUS=reint(1)
       /PRINT MEAN
       /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95).

# Kaplan-Meier

Notes
Output Created | 19-NOV-2002 12:40:16
Comments |
Input | Data | d:\temp4.sav
Filter | <none>
Weight | <none>
Split File | <none>
N of Rows in Working Data File | 503965
Syntax | KM time /STATUS=reint(1) /PRINT MEAN /PERCENTILES = (45,50,55,60,65,70,75,80,85,90,95).
Resources | Elapsed Time | 0:00:40.91

     Survival Analysis for TIME

              Survival Time    Standard Error   95% Confidence Interval

     Mean:      2537.56                2.76     (  2532.14,   2542.98 )
     (Limited to   4016.0 )
     Median:    3003.00               14.22     (  2975.13,   3030.87 )

                                                                Percentiles

                         45.00        50.00        55.00        60.00        65.00        70.00        75.00        80.00        85.00

     Value             3758.00      3003.00      2431.00      1965.00      1564.00      1231.00       949.00       709.00       500.00
     Standard Error        .          14.22         9.90         7.58         5.92         4.56         3.53         2.72         2.04

                                                                Percentiles

                         90.00        95.00

     Value              321.00       168.00
     Standard Error       1.42          .

    **************************************************************************************.

    * Create an Excel file for comparing survivor functions.
    get file='d:\temp5.sav'.

    weight off.

    select if (reint=1).
    sort cases by time.
    aggregate outfile=*
     /presorted
     /break time
     / survadj survunad survrand
    = first( survadj survunad survrand ).

    SAVE TRANSLATE OUTFILE='I:\Research Projects\longevity\1991 to 2001 adj unadj and rand.xls'
      /TYPE=XLS /MAP /REPLACE /FIELDNAMES.

    Data written to I:\Research Projects\longevity\1991 to 2001 adj unadj and rand.xls.
    4 variables and 3783 cases written to range: SPSS.
    Variable: TIME       Type: Number   Width:  8   Dec: 2
    Variable: SURVADJ    Type: Number   Width: 10   Dec: 5
    Variable: SURVUNAD   Type: Number   Width: 10   Dec: 5
    Variable: SURVRAND   Type: Number   Width: 10   Dec: 5

    * Create an Excel file for calculation of standard errors.

    get file='d:\temp4.sav'.

    weight by weight.

    sort cases by time reint.

    save outfile = 'd:\temp6.sav'.
    get file='d:\temp6.sav'.

    compute censored=0.
    if (reint=0) censored=1.

    aggregate outfile=*
     /presorted
     /break time
     /reints censored =sum(reint censored ).

    >Warning # 3211
    >On at least one case, the value of the weight variable was zero, negative,
    >or missing.  Such cases are invisible to statistical procedures and graphs
    >which need positively weighted cases, but remain on the file and are
    >processed by non-statistical facilities such as LIST and SAVE.

    compute reints=reints/1000.
    compute censored = censored/1000.

    SAVE TRANSLATE OUTFILE='I:\Research Projects\longevity\1991 to 2001 times for se.xls'
      /TYPE=XLS /MAP /REPLACE /FIELDNAMES.

    Data written to I:\Research Projects\longevity\1991 to 2001 times for se.xls.
    3 variables and 4012 cases written to range: SPSS.
    Variable: TIME       Type: Number   Width:  8   Dec: 2
    Variable: REINTS     Type: Number   Width:  8   Dec: 2
    Variable: CENSORED   Type: Number   Width:  8   Dec: 2