Basic statistics exercises (survival analysis)

par(mar=c(10,10,5,5))
fw <- prodlim(Surv(ttp,status)~smokeWoman,data=ttp)
plot(fw,
     type="cuminc",
     ylab="Probability of getting pregnant",
     legend.title="Woman smoking",
     atrisk.title="Woman smoking",
     legend.x="topleft",
     atrisk.at=seq(0,200,50),
     atrisk.pos=-20,
     xlab="Days")

logRankWoman <- survdiff(Surv(ttp,status)~smokeWoman,data=ttp)
publish(logRankWoman,org=TRUE)

X11cairo 
       2
 Log-rank test      N Observed Expected squared(O-E)/E squared(O-E)/V 
 smokeWoman=No 300.00   194.00   176.46           1.74           5.64 
 smokeWoman=Yes 123.00    62.00    79.54           3.87           5.64 

 Chisq= 5.6  on 1 degrees of freedom, p= 0.0176

The chance of getting pregnant was significantly higher when the woman was not smoking compared to when the woman was smoking (Log-rank p-value: 0.0176).

par(mar=c(10,10,5,5))
fm <- prodlim(Surv(ttp,status)~smokeMan,data=ttp)
plot(fm,
     type="cuminc",
     ylab="Probability of getting pregnant",
     legend.x="topleft",
     legend.title="Man smoking",
     atrisk.title="Man smoking",
     atrisk.at=seq(0,200,50),
     atrisk.pos=-20,
     xlab="Days")

logRankMan <- survdiff(Surv(ttp,status)~smokeMan,data=ttp)
publish(logRankMan,org=TRUE)

X11cairo 
       2
 Log-rank test      N Observed Expected squared(O-E)/E squared(O-E)/V 
 smokeMan=No 290.00   183.00   172.31           0.66           2.04 
 smokeMan=Yes 133.00    73.00    83.69           1.37           2.04 

 Chisq= 2  on 1 degrees of freedom, p= 0.1536

The chance of getting pregnant was higher when the man was not smoking compared to when the man was smoking (see Kaplan-Meier plot), however, this difference was not statistically significant (Log-rank p-value: 0.1536).

coxWoman <- coxph(Surv(ttp,status)~smokeWoman,data=ttp)
publish(coxWoman,org=TRUE)

Smoking of the woman decreased the hazard and hence increased the time to pregnancy significantly (hazard ratio=0.71, 95%-CI=[0.53;0.94], p=0.018).

coxMan <- coxph(Surv(ttp,status)~smokeMan,data=ttp)
publish(coxMan,org=TRUE)

Smoking of the man did not significantly increase the time to pregnancy (hazard ratio=0.82, 95%-CI=[0.63;1.08], p=0.15).

coxInter <- coxph(Surv(ttp,status)~smokeMan*smokeWoman,data=ttp,x=TRUE)
summary(coxInter)

Call:
coxph(formula = Surv(ttp, status) ~ smokeMan * smokeWoman, data = ttp, 
    x = TRUE)

  n= 423, number of events= 256 

                             coef exp(coef) se(coef)      z Pr(>|z|)  
smokeManYes               -0.1178    0.8889   0.1869 -0.630   0.5285  
smokeWomanYes             -0.3695    0.6911   0.2192 -1.686   0.0918 .
smokeManYes:smokeWomanYes  0.1210    1.1287   0.3209  0.377   0.7061  
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                          exp(coef) exp(-coef) lower .95 upper .95
smokeManYes                  0.8889      1.125    0.6162     1.282
smokeWomanYes                0.6911      1.447    0.4497     1.062
smokeManYes:smokeWomanYes    1.1287      0.886    0.6018     2.117

Concordance= 0.536  (se = 0.017 )
Rsquare= 0.015   (max possible= 0.999 )
Likelihood ratio test= 6.31  on 3 df,   p=0.1
Wald test            = 6.01  on 3 df,   p=0.1
Score (logrank) test = 6.07  on 3 df,   p=0.1

publish(coxInter,org=TRUE)

There is no significant effect modification. The hazard ratio for non-smoking couples compared to smoking couples is:

exp(coef_1 + coef_2 + coef_3) = exp(-0.1178-0.3695+0.121)=
 0.8889 * 0.6911 * 1.1287 = 0.693

This is similar to the hazard ratio for smoking woman compared to non-smoking woman when the man is not smoking:

exp(coef_2) = exp(0.3695)= 0.6911

The p-value is 0.7061.

We can obtain the comparison between both smoking and none smoking by including the variable coupleSmoking (which has 4 levels):

publish(coxph(Surv(ttp,status)~coupleSmoking,data=ttp),org=TRUE)

coxDrink <- coxph(Surv(ttp,status)~alcWoman,data=ttp)
publish(coxDrink,org=TRUE)

The effect of alcohol is not statistically significant (p=0.132) but this model assumes that the relation between a unit change of alcohol drinking behaviour (1 more beverage per week) and the hazard ratio is linear. This assumption may be violated.

ttp[,drinkWoman:=factor(cut(alcWoman,c(-1,0,10,100)),labels=c("no","moderate","high"))]
coxDrink2 <- coxph(Surv(ttp,status)~drinkWoman,data=ttp)
publish(coxDrink2,org=TRUE)

The probability of pregnancy is significantly lower for women drinking moderately as compared to women not drinking (hazard ratio: 0.73, 95%-CI: [0.55;0.95], p-value: 0.021) and for women drinking excessively as compared to women not drinking (hazard ratio: 0.53, 95%-CI: [0.32;0.89], p-value: 0.016).

The comparison of excessively to moderate is obtained by using relevel:

ttp$drinkWoman2 <- relevel(ttp$drinkWoman,"moderate")
coxDrink3 <- coxph(Surv(ttp,status)~drinkWoman2,data=ttp)
publish(coxDrink3,org=TRUE)

The probability of pregnancy is estimated as being lower, but not statistically significantly lower, for women drinking excessively as compared to women drinking moderately (hazard ratio: 0.73, 95%-CI: [0.45;0.95], p-value: 0.21).

coxDrink4 <- coxph(Surv(ttp,status)~drinkWoman+smokeWoman,data=ttp)
publish(coxDrink4,org=TRUE)

After adjusting for smoking the results are similar as above.

coxDrink5 <- coxph(Surv(ttp,status)~drinkWoman+spermConc,data=ttp)
publish(coxDrink5,org=TRUE)

After adjusting for sperm concentration the results become more significant and the hazard ratios for group comparisons are more extreme. The reasons are

1. sperm concentration is an important factor which when omitted leads to a systematic underestimation of effects. Like in logistic regression, this occurs even if the variable of interest (i.e. the drinking behavior of the woman) is independent of the factor (i.e. the sperm concentration).
2. sperm concentration is missing for 113 of the men. The default behavior of the software is to remove all the data of all couples where the sperm concentration is missing. This changes the data set quite a lot and there may be a selection bias.