7 Further topics

Guinea-Bissau childhood vaccination study

Table 7.2

Please note, that \(\widehat{\beta}\) will vary with chosen seed, but SD values are stable.

Code show/hide

# Cox model - full cohort
bissau <- read.csv("data/bissau.csv")
bissau$agem <- as.integer(bissau$age/30.4)

library(survival)
options(contrasts=c("contr.treatment", "contr.poly"))
coxph(Surv(fuptime, dead != 0) ~ factor(bcg) + factor(agem), data = bissau)

Call:
coxph(formula = Surv(fuptime, dead != 0) ~ factor(bcg) + factor(agem), 
    data = bissau)

                  coef exp(coef) se(coef)      z      p
factor(bcg)1  -0.34720   0.70667  0.14605 -2.377 0.0174
factor(agem)1  0.11500   1.12187  0.23205  0.496 0.6202
factor(agem)2 -0.25687   0.77347  0.25861 -0.993 0.3206
factor(agem)3  0.19894   1.22011  0.24325  0.818 0.4135
factor(agem)4  0.33252   1.39447  0.24183  1.375 0.1691
factor(agem)5  0.33066   1.39189  0.24957  1.325 0.1852
factor(agem)6 -0.01052   0.98953  0.35340 -0.030 0.9762

Likelihood ratio test=12.37  on 7 df, p=0.08901
n= 5274, number of events= 222

Code show/hide

# 12.5 pct sub-cohort
bissau$s <- rbinom(nrow(bissau), 1, prob = 0.125)
with(bissau, table(dead, s))

    s
dead    0    1
   0 4444  608
   1  189   33

Code show/hide

# Case-cohort data
bis2 <- subset(bissau, s == 0 & dead != 0 | s == 1)
ccfit <- cch(Surv(fuptime, dead!= 0) ~ factor(bcg) + factor(agem), 
             data = bis2, 
             subcoh=bis2$s, id=bis2$id, cohort.size=nrow(bissau))
summary(ccfit)

Case-cohort analysis,x$method, Prentice 
 with subcohort of 641 from cohort of 5274 

Call: cch(formula = Surv(fuptime, dead != 0) ~ factor(bcg) + factor(agem), 
    data = bis2, subcoh = bis2$s, id = bis2$id, cohort.size = nrow(bissau))

Coefficients:
                Coef    HR  (95%   CI)     p
factor(bcg)1  -0.296 0.744 0.531 1.042 0.085
factor(agem)1 -0.042 0.959 0.564 1.630 0.877
factor(agem)2 -0.456 0.634 0.358 1.123 0.118
factor(agem)3  0.284 1.329 0.757 2.331 0.322
factor(agem)4  0.134 1.144 0.656 1.992 0.636
factor(agem)5  0.521 1.684 0.938 3.024 0.081
factor(agem)6  0.166 1.181 0.527 2.643 0.686

Code show/hide

# Nested-case control
library(multipleNCC)
library(Epi)

# Add noise to remove ties

bis2$fuptime_noise <- jitter(bis2$fuptime, factor=1, amount = NULL)

nccdata <- Epi::ccwc(exit=fuptime_noise, 
                     fail=dead!=0, 
                     data=bis2, 
                     include=list(bcg, agem), controls=3, silent=TRUE)
summary(clogit(Fail ~ factor(bcg) + factor(agem), data=nccdata))

Call:
coxph(formula = Surv(rep(1, 888L), Fail) ~ factor(bcg) + factor(agem), 
    data = nccdata, method = "exact")

  n= 888, number of events= 222 

                  coef exp(coef) se(coef)      z Pr(>|z|)  
factor(bcg)1  -0.31264   0.73151  0.17525 -1.784   0.0744 .
factor(agem)1  0.08065   1.08399  0.26927  0.300   0.7645  
factor(agem)2 -0.10850   0.89718  0.29427 -0.369   0.7123  
factor(agem)3  0.49968   1.64819  0.29021  1.722   0.0851 .
factor(agem)4  0.27275   1.31357  0.28741  0.949   0.3426  
factor(agem)5  0.50070   1.64988  0.29790  1.681   0.0928 .
factor(agem)6  0.29128   1.33814  0.41155  0.708   0.4791  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

              exp(coef) exp(-coef) lower .95 upper .95
factor(bcg)1     0.7315     1.3670    0.5189     1.031
factor(agem)1    1.0840     0.9225    0.6395     1.838
factor(agem)2    0.8972     1.1146    0.5040     1.597
factor(agem)3    1.6482     0.6067    0.9332     2.911
factor(agem)4    1.3136     0.7613    0.7478     2.307
factor(agem)5    1.6499     0.6061    0.9202     2.958
factor(agem)6    1.3381     0.7473    0.5973     2.998

Concordance= 0.572  (se = 0.022 )
Likelihood ratio test= 8.79  on 7 df,   p=0.3
Wald test            = 8.68  on 7 df,   p=0.3
Score (logrank) test = 8.79  on 7 df,   p=0.3

Code show/hide

proc import 
datafile="data/bissau.csv"
    out=bissau
    dbms = csv
    replace;
run;
data bissau; 
    set bissau; 
    agem = int(age/30.4);
run; 

/* 1.: Fit Cox model */
proc phreg data=bissau;
  class bcg(ref="0") agem;
  model fuptime*dead(0)=bcg agem / rl;
run;

                  Parameter    Standard
Parameter    DF    Estimate       Error

bcg       1   1    -0.34720     0.14605
agem      0   1     0.01053     0.35339
agem      1   1     0.12553     0.34494
agem      2   1    -0.24631     0.35903
agem      3   1     0.20946     0.34502
agem      4   1     0.34300     0.34265
agem      5   1     0.34118     0.34699


/* 2.: Create a 12.5 pct sub-cohort */

data bissau; 
    set bissau;
    seed=260452;
    s=ranbin(seed,1,0.125);
run;
proc freq; 
    tables s*dead; 
run;

/* 3.: Fit Cox model to case-cohort data */
data casecoho; set bissau;
    epsilon=0.001;
    if dead=1 and s=0 then do;
    d=1; start=fuptime-epsilon; stop=fuptime; w=1; output; end;
    if dead=0 and s=1 then do;
    d=0; start=0; stop=fuptime; w=1/0.125; output; end;
    if dead=1 and s=1 then do;
    d=0; start=0; stop=fuptime-epsilon; w=1/0.125; output;
    d=1; start=fuptime-epsilon; stop=fuptime; w=1; output; end;
run;

proc phreg data=casecoho covsandwich(aggregate);
    class bcg(ref="0") agem;
    model stop*d(0)=bcg agem/rl entry=start;
    weight w; id id;
run;

                 Parameter   Standard
 Parameter   DF   Estimate      Error

 bcg       1  1   -0.38946    0.16599
 agem      0  1   -0.00371    0.41199
 agem      1  1    0.18793    0.39633
 agem      2  1   -0.33429    0.40462
 agem      3  1    0.09770    0.39586
 agem      4  1    0.14055    0.38892
 agem      5  1    0.09873    0.39251


/* NCC MACRO */

data source; set bissau;
    study_id=id;
    age_entry=0;
    age_dlo=fuptime;
    censor=dead;
run;

%macro caseset;
    %let sampling = 1;
    %let ratio = 3;
    /* Enumerate Cases */
    data cases;
    set source ;
    if censor = 1;
    run;
    data cases;
    set cases end = eof;
    if eof then call symput ('ncases', put(_n_,6.));
    run;
    /* Create Risk Set */
    %do iter = 1 %to &ncases;
    data temp_case;
    set cases;

    if _n_ = &iter ;
    call symput ('rs', put(_n_,6.));
    call symput ('age_rs', put(age_dlo,8.)); call symput ('case_id',
    put(study_id,8.));
    run;
    data temp_control;
    set source;
    if age_entry  <= &age_rs  <= age_dlo;
    /* Exclude Index Case */
    if study_id = &case_id then delete;
    number = ranuni(0);
    age_rs = &age_rs;
    censor = 0;
    run;
    /* Sample Controls */
    %if &sampling = 1 %then %do;
    proc sort data = temp_control;
    by number;
    data temp_control;
    set temp_control;
    by age_rs;
    retain m;
    if first.age_rs then m = 0;
    m=m+1;
    if m <= &ratio then output temp_control;
    run;
    %end; 
    /* End If Sampling = 1 */
    /* Combine Case with Controls */
    data rs_&iter;
    set temp_case
    temp_control;
    rs = &rs;
    age_rs = &age_rs;
    run;
    /* DM Output 'Clear'; Log 'Clear'; */

    %end; 
    /* End Loop Creating Risk Set */
    /* Append Risk Sets */
    %do j = 2 %to &ncases;
    proc append base = rs_1 data = rs_&j;
    run;
    %end;
    data final; set rs_1; run;
%mend ; 
/* End Macro */

/* Invoke Macro */
%caseset;

proc phreg data=final;
class bcg(ref="0") agem;
  model fuptime*censor(0)=bcg agem / rl;
  strata rs;
run;