First, we read the data:
library(haven)
worktrips <- read_dta("worktrips.dta")
attach(worktrips)dil = sd >= 0
d2l = sd >= flex
sgl = nperson == 1
cp = cpnew
fl = flex > 0
sdlx = pmax(sd - flex, rep(0, length(sd)))
correct = alt == choice
timsgl = tim * sgl
timcp = tim * cp
sdesgl = sde * sgl
sdecp = sde * cp
sdesde = sde * sde
sdlwc = sdl * wc
sdlfl = sdl * fl
sdlsdl = sdl * sdl
dilwc = dil * wc
dilfl = dil * fllibrary(survival)model1 <- clogit(correct ~ r15 + r10 + tim + sde + sdl + dil + strata(id))
model2 <- clogit(correct ~ r15 + r10 + tim + timsgl + sde + sdesgl + sdl + sdlwc + dil + dilwc + strata(id))
model3 <- clogit(correct ~ r15 + r10 + tim + timcp + sde + sdecp + sdl + sdlfl + dil + dilfl + strata(id))
model4 <- clogit(correct ~ r15 + r10 + tim + timcp + sde + sdecp + sdl + sdlx + d2l + strata(id))
model5 <- clogit(correct ~ r15 + r10 + tim + timsgl + timcp + sde + sdesgl + sdecp + sdl + sdlwc + sdlx + dilwc + d2l + strata(id))
model4a <- clogit(correct ~ r15 + r10 + tim + timcp + sde + sdecp + sdesde + sdl + sdlx + sdlsdl + d2l + strata(id))
model4b <- clogit(correct ~ r15 + r10 + tim + timcp + sdl + sdlx + d2l + strata(id))
model4c <- clogit(correct ~ r15 + r10 + tim + timcp + sde + sdecp + d2l + strata(id))
model4d <- clogit(correct ~ r15 + r10 + tim + sde + sdl + sdlx + d2l + strata(id)) summary(model1)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## sde + sdl + dil + strata(id), method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.085933 2.962203 0.116568 9.316 < 2e-16 ***
## r10 0.377249 1.458268 0.118455 3.185 0.00145 **
## tim -0.106414 0.899052 0.037926 -2.806 0.00502 **
## sde -0.065307 0.936779 0.006994 -9.337 < 2e-16 ***
## sdl -0.253824 0.775828 0.029837 -8.507 < 2e-16 ***
## dilTRUE -0.584162 0.557573 0.206587 -2.828 0.00469 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 2.9622 0.3376 2.3572 3.7225
## r10 1.4583 0.6857 1.1561 1.8394
## tim 0.8991 1.1123 0.8346 0.9684
## sde 0.9368 1.0675 0.9240 0.9497
## sdl 0.7758 1.2889 0.7318 0.8226
## dilTRUE 0.5576 1.7935 0.3719 0.8359
##
## Concordance= 0.8 (se = 0.011 )
## Likelihood ratio test= 570.6 on 6 df, p=<2e-16
## Wald test = 430.2 on 6 df, p=<2e-16
## Score (logrank) test = 646.9 on 6 df, p=<2e-16summary(model2)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timsgl + sde + sdesgl + sdl + sdlwc + dil + dilwc + strata(id),
## method = "exact")
##
## n= 6264, number of events= 522
## (60 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.06999 2.91536 0.11902 8.990 < 2e-16 ***
## r10 0.37865 1.46031 0.12036 3.146 0.001656 **
## tim -0.03821 0.96252 0.04174 -0.915 0.359970
## timsgl -0.27226 0.76166 0.10663 -2.553 0.010675 *
## sde -0.06326 0.93870 0.00721 -8.774 < 2e-16 ***
## sdesgl -0.02150 0.97872 0.01150 -1.870 0.061544 .
## sdl -0.53033 0.58841 0.18351 -2.890 0.003853 **
## sdlwc 0.29006 1.33651 0.18425 1.574 0.115423
## dilTRUE -1.10842 0.33008 0.25520 -4.343 1.4e-05 ***
## dilwc 0.78223 2.18635 0.20289 3.855 0.000116 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 2.9154 0.3430 2.3088 3.6813
## r10 1.4603 0.6848 1.1534 1.8488
## tim 0.9625 1.0389 0.8869 1.0446
## timsgl 0.7617 1.3129 0.6180 0.9387
## sde 0.9387 1.0653 0.9255 0.9521
## sdesgl 0.9787 1.0217 0.9569 1.0010
## sdl 0.5884 1.6995 0.4107 0.8431
## sdlwc 1.3365 0.7482 0.9314 1.9178
## dilTRUE 0.3301 3.0296 0.2002 0.5443
## dilwc 2.1863 0.4574 1.4690 3.2540
##
## Concordance= 0.806 (se = 0.011 )
## Likelihood ratio test= 613.9 on 10 df, p=<2e-16
## Wald test = 430.4 on 10 df, p=<2e-16
## Score (logrank) test = 696.9 on 10 df, p=<2e-16summary(model3)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timcp + sde + sdecp + sdl + sdlfl + dil + dilfl + strata(id),
## method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.068173 2.910057 0.117042 9.126 < 2e-16 ***
## r10 0.364180 1.439333 0.119045 3.059 0.00222 **
## tim -0.137788 0.871283 0.053080 -2.596 0.00944 **
## timcp 0.106028 1.111853 0.075614 1.402 0.16085
## sde -0.072054 0.930481 0.007778 -9.264 < 2e-16 ***
## sdecp 0.022723 1.022984 0.008805 2.581 0.00986 **
## sdl -0.389379 0.677477 0.077290 -5.038 4.71e-07 ***
## sdlfl 0.181943 1.199546 0.080423 2.262 0.02368 *
## dilTRUE -1.005277 0.365943 0.225256 -4.463 8.09e-06 ***
## dilfl 1.185794 3.273284 0.196298 6.041 1.53e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 2.9101 0.3436 2.3135 3.6604
## r10 1.4393 0.6948 1.1398 1.8176
## tim 0.8713 1.1477 0.7852 0.9668
## timcp 1.1119 0.8994 0.9587 1.2895
## sde 0.9305 1.0747 0.9164 0.9448
## sdecp 1.0230 0.9775 1.0055 1.0408
## sdl 0.6775 1.4761 0.5822 0.7883
## sdlfl 1.1995 0.8336 1.0246 1.4043
## dilTRUE 0.3659 2.7327 0.2353 0.5691
## dilfl 3.2733 0.3055 2.2279 4.8092
##
## Concordance= 0.803 (se = 0.011 )
## Likelihood ratio test= 630.7 on 10 df, p=<2e-16
## Wald test = 443.2 on 10 df, p=<2e-16
## Score (logrank) test = 733.6 on 10 df, p=<2e-16summary(model4)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timcp + sde + sdecp + sdl + sdlx + d2l + strata(id), method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.105662 3.021225 0.100690 10.981 < 2e-16 ***
## r10 0.397988 1.488826 0.101581 3.918 8.93e-05 ***
## tim -0.140818 0.868647 0.052754 -2.669 0.00760 **
## timcp 0.105238 1.110975 0.075675 1.391 0.16433
## sde -0.074685 0.928036 0.006095 -12.253 < 2e-16 ***
## sdecp 0.022521 1.022776 0.008826 2.552 0.01072 *
## sdl -0.174396 0.839964 0.029097 -5.994 2.05e-09 ***
## sdlx -0.216150 0.805615 0.080719 -2.678 0.00741 **
## d2lTRUE -1.057030 0.347486 0.169984 -6.218 5.02e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 3.0212 0.3310 2.4801 3.6804
## r10 1.4888 0.6717 1.2201 1.8168
## tim 0.8686 1.1512 0.7833 0.9633
## timcp 1.1110 0.9001 0.9578 1.2886
## sde 0.9280 1.0775 0.9170 0.9392
## sdecp 1.0228 0.9777 1.0052 1.0406
## sdl 0.8400 1.1905 0.7934 0.8893
## sdlx 0.8056 1.2413 0.6877 0.9437
## d2lTRUE 0.3475 2.8778 0.2490 0.4849
##
## Concordance= 0.802 (se = 0.011 )
## Likelihood ratio test= 629.4 on 9 df, p=<2e-16
## Wald test = 426.2 on 9 df, p=<2e-16
## Score (logrank) test = 638.4 on 9 df, p=<2e-16summary(model5)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timsgl + timcp + sde + sdesgl + sdecp + sdl + sdlwc + sdlx +
## dilwc + d2l + strata(id), method = "exact")
##
## n= 6264, number of events= 522
## (60 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.002270 2.724459 0.106913 9.375 < 2e-16 ***
## r10 0.308630 1.361559 0.107990 2.858 0.00426 **
## tim -0.060062 0.941706 0.059780 -1.005 0.31503
## timsgl -0.223491 0.799722 0.113686 -1.966 0.04932 *
## timcp 0.068932 1.071363 0.079236 0.870 0.38433
## sde -0.065149 0.936927 0.006881 -9.468 < 2e-16 ***
## sdesgl -0.019092 0.981089 0.011812 -1.616 0.10602
## sdecp 0.020814 1.021032 0.008924 2.332 0.01968 *
## sdl -0.466322 0.627306 0.185719 -2.511 0.01204 *
## sdlwc 0.277764 1.320175 0.186192 1.492 0.13575
## sdlx -0.176544 0.838162 0.081957 -2.154 0.03123 *
## dilwc 0.504850 1.656737 0.174387 2.895 0.00379 **
## d2lTRUE -1.166821 0.311355 0.175058 -6.665 2.64e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 2.7245 0.3670 2.2094 3.3596
## r10 1.3616 0.7345 1.1018 1.6825
## tim 0.9417 1.0619 0.8376 1.0588
## timsgl 0.7997 1.2504 0.6400 0.9993
## timcp 1.0714 0.9334 0.9173 1.2514
## sde 0.9369 1.0673 0.9244 0.9496
## sdesgl 0.9811 1.0193 0.9586 1.0041
## sdecp 1.0210 0.9794 1.0033 1.0390
## sdl 0.6273 1.5941 0.4359 0.9027
## sdlwc 1.3202 0.7575 0.9165 1.9016
## sdlx 0.8382 1.1931 0.7138 0.9842
## dilwc 1.6567 0.6036 1.1771 2.3318
## d2lTRUE 0.3114 3.2118 0.2209 0.4388
##
## Concordance= 0.808 (se = 0.01 )
## Likelihood ratio test= 658.5 on 13 df, p=<2e-16
## Wald test = 445.1 on 13 df, p=<2e-16
## Score (logrank) test = 735.5 on 13 df, p=<2e-16summary(model4a)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timcp + sde + sdecp + sdesde + sdl + sdlx + sdlsdl + d2l +
## strata(id), method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.1429535 3.1360169 0.1048154 10.904 < 2e-16 ***
## r10 0.3924253 1.4805672 0.1086180 3.613 0.000303 ***
## tim -0.1361214 0.8727367 0.0530121 -2.568 0.010236 *
## timcp 0.1054547 1.1112158 0.0753353 1.400 0.161572
## sde -0.0843917 0.9190712 0.0165520 -5.099 3.42e-07 ***
## sdecp 0.0223924 1.0226450 0.0087450 2.561 0.010450 *
## sdesde 0.0002991 1.0002991 0.0004577 0.653 0.513495
## sdl -0.0442571 0.9567080 0.0907775 -0.488 0.625880
## sdlx -0.2428262 0.7844079 0.0887560 -2.736 0.006221 **
## sdlsdl -0.0115881 0.9884787 0.0078007 -1.486 0.137403
## d2lTRUE -1.1307766 0.3227825 0.1849109 -6.115 9.64e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 3.1360 0.3189 2.5536 3.8512
## r10 1.4806 0.6754 1.1967 1.8318
## tim 0.8727 1.1458 0.7866 0.9683
## timcp 1.1112 0.8999 0.9587 1.2880
## sde 0.9191 1.0881 0.8897 0.9494
## sdecp 1.0226 0.9779 1.0053 1.0403
## sdesde 1.0003 0.9997 0.9994 1.0012
## sdl 0.9567 1.0453 0.8008 1.1430
## sdlx 0.7844 1.2748 0.6592 0.9335
## sdlsdl 0.9885 1.0117 0.9735 1.0037
## d2lTRUE 0.3228 3.0981 0.2247 0.4638
##
## Concordance= 0.803 (se = 0.011 )
## Likelihood ratio test= 632.4 on 11 df, p=<2e-16
## Wald test = 427.1 on 11 df, p=<2e-16
## Score (logrank) test = 700.3 on 11 df, p=<2e-16summary(model4b)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timcp + sdl + sdlx + d2l + strata(id), method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.14388 3.13892 0.09675 11.823 < 2e-16 ***
## r10 0.20987 1.23352 0.09823 2.137 0.03263 *
## tim -0.02549 0.97484 0.04108 -0.620 0.53500
## timcp 0.03065 1.03113 0.06107 0.502 0.61571
## sdl -0.11949 0.88738 0.02689 -4.443 8.87e-06 ***
## sdlx -0.28813 0.74966 0.08239 -3.497 0.00047 ***
## d2lTRUE 0.30695 1.35927 0.13560 2.264 0.02360 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 3.1389 0.3186 2.5967 3.7943
## r10 1.2335 0.8107 1.0175 1.4954
## tim 0.9748 1.0258 0.8994 1.0566
## timcp 1.0311 0.9698 0.9148 1.1622
## sdl 0.8874 1.1269 0.8418 0.9354
## sdlx 0.7497 1.3339 0.6379 0.8810
## d2lTRUE 1.3593 0.7357 1.0420 1.7731
##
## Concordance= 0.751 (se = 0.013 )
## Likelihood ratio test= 414.5 on 7 df, p=<2e-16
## Wald test = 270.8 on 7 df, p=<2e-16
## Score (logrank) test = 422.3 on 7 df, p=<2e-16summary(model4c)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## timcp + sde + sdecp + d2l + strata(id), method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.388717 4.009702 0.095049 14.611 < 2e-16 ***
## r10 0.870564 2.388258 0.094356 9.226 < 2e-16 ***
## tim -0.073009 0.929592 0.046858 -1.558 0.11921
## timcp 0.084036 1.087668 0.068525 1.226 0.22006
## sde -0.072420 0.930140 0.006148 -11.779 < 2e-16 ***
## sdecp 0.022926 1.023191 0.008651 2.650 0.00804 **
## d2lTRUE -1.986469 0.137179 0.157241 -12.633 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 4.0097 0.2494 3.3282 4.8308
## r10 2.3883 0.4187 1.9850 2.8734
## tim 0.9296 1.0757 0.8480 1.0190
## timcp 1.0877 0.9194 0.9510 1.2440
## sde 0.9301 1.0751 0.9190 0.9414
## sdecp 1.0232 0.9773 1.0060 1.0407
## d2lTRUE 0.1372 7.2897 0.1008 0.1867
##
## Concordance= 0.769 (se = 0.012 )
## Likelihood ratio test= 464.2 on 7 df, p=<2e-16
## Wald test = 367.8 on 7 df, p=<2e-16
## Score (logrank) test = 433.4 on 7 df, p=<2e-16summary(model4d)## Call:
## coxph(formula = Surv(rep(1, 6324L), correct) ~ r15 + r10 + tim +
## sde + sdl + sdlx + d2l + strata(id), method = "exact")
##
## n= 6324, number of events= 527
##
## coef exp(coef) se(coef) z Pr(>|z|)
## r15 1.101108 3.007496 0.100566 10.949 < 2e-16 ***
## r10 0.395952 1.485798 0.101460 3.903 9.52e-05 ***
## tim -0.092756 0.911416 0.038785 -2.392 0.0168 *
## sde -0.066761 0.935419 0.005064 -13.184 < 2e-16 ***
## sdl -0.174547 0.839837 0.029081 -6.002 1.95e-09 ***
## sdlx -0.216847 0.805053 0.080821 -2.683 0.0073 **
## d2lTRUE -1.029510 0.357182 0.168840 -6.098 1.08e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## r15 3.0075 0.3325 2.4695 3.6628
## r10 1.4858 0.6730 1.2179 1.8127
## tim 0.9114 1.0972 0.8447 0.9834
## sde 0.9354 1.0690 0.9262 0.9447
## sdl 0.8398 1.1907 0.7933 0.8891
## sdlx 0.8051 1.2422 0.6871 0.9432
## d2lTRUE 0.3572 2.7997 0.2566 0.4973
##
## Concordance= 0.805 (se = 0.011 )
## Likelihood ratio test= 622.1 on 7 df, p=<2e-16
## Wald test = 425.6 on 7 df, p=<2e-16
## Score (logrank) test = 635 on 7 df, p=<2e-16