Author: Andreas Karlsson, 2015-03-06
These data were used to study a possible effect of marital bereavement (loss of husband or wife) on all–cause mortality in the elderly. The dataset was extracted from a larger follow-up study of an elderly population and concerns subjects whose husbands or wives were alive at entry to the study. Thus all subjects enter as not bereaved but may become bereaved at some point during follow–up. The variable dosp records the date of death of each subject’s spouse and takes the value 1/1/2000 where this has not yet happened.
You may have to install the required packages the first time you use them. You can install a package by install.packages("package_of_interest") for each package you require.
require(foreign) # for reading data set from Stata
require(survival) # for Surv and survfit
require(dplyr) # for data manipulation
###########################################
### A help function to calculate ###
### and print incidence (hazard) ratios
### from a fitted poisson regression
### from glm
###########################################
IRR <- function(fit){
summfit <- summary(fit )$coefficients
IRfit <- exp( cbind( summfit[, 1:2], summfit[, 1] - 1.96*summfit[, 2], summfit[, 1] +
1.96*summfit[, 2] ) )
colnames(IRfit) <- c("IRR", "Std. err", "CI_lower", "CI_upper")
print(IRfit)
}(a) Load the bereavement data and explore it.
brv_raw <- read.dta("http://biostat3.net/download/brv.dta")head(brv_raw)## id couple dob doe dox dosp fail group
## 1 14464 93 1901-01-31 1981-03-11 1990-11-10 2000-01-01 1 1
## 2 2000 18 1896-06-27 1981-03-19 1981-10-02 1982-03-12 1 2
## 3 24486 187 1905-06-03 1981-04-13 1991-01-01 1982-03-03 0 2
## 4 24646 195 1895-11-20 1981-04-29 1985-08-17 2000-01-01 1 1
## 5 16810 119 1898-12-26 1981-03-18 1989-08-31 2000-01-01 1 2
## 6 6344 49 1904-10-31 1981-05-22 1989-02-19 1986-10-05 1 1
## disab health sex
## 1 0 2 1
## 2 3 1 1
## 3 0 2 1
## 4 0 2 1
## 5 0 2 1
## 6 3 2 1##Look at the five first couples
brv_raw %>% select(couple, id, sex, doe, dosp, dox, fail) %>% filter(couple<=5) %>% arrange(couple, id)## couple id sex doe dosp dox fail
## 1 1 48 1 1981-01-15 2000-01-01 1983-04-18 1
## 2 2 49 2 1981-04-27 1983-01-24 1984-11-20 1
## 3 2 263 1 1981-04-27 1984-11-20 1983-01-24 1
## 4 3 60 1 1981-01-20 1981-12-31 1981-08-03 1
## 5 3 63 2 1981-01-20 1981-08-03 1981-12-31 1
## 6 4 156 1 1981-01-20 1988-11-23 1991-01-01 0
## 7 4 220 2 1981-01-20 2000-01-01 1988-11-23 1
## 8 5 361 1 1981-05-07 1984-02-04 1984-02-04 1(b) Calculate the mortality rate per 1000 years for men and for women, and find the rate ratio comparing women (coded 2) with men (coded 1).
brv <- brv_raw %>%
mutate(age_entry = as.numeric(doe - dob) / 365.24, # Calc age at entry
att_age = as.numeric(dox - dob) / 365.24, # Calc attained age
t_at_risk = att_age - age_entry) # Calc time at risk
## crude rates
brv %>%
select(fail, t_at_risk, sex) %>%
group_by(sex) %>%
summarise(D = sum(fail), Y = sum(t_at_risk)/1000, Rate = D/Y) %>%
mutate(CI_low = Rate + qnorm(0.025) * Rate / sqrt(D),
CI_high = Rate + qnorm(0.975) * Rate / sqrt(D))## Source: local data frame [2 x 6]
##
## sex D Y Rate CI_low CI_high
## 1 1 181 1.340521 135.02210 115.3517 154.6925
## 2 2 97 1.095187 88.56937 70.9437 106.1950poisson22b <- glm( fail ~ sex + offset( log( t_at_risk) ), family=poisson, data=brv)
summary( poisson22b )##
## Call:
## glm(formula = fail ~ sex + offset(log(t_at_risk)), family = poisson,
## data = brv)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6389 -1.3080 0.2636 0.8742 3.1495
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.5807 0.1800 -8.780 < 2e-16 ***
## sex -0.4217 0.1258 -3.351 0.000806 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 526.28 on 398 degrees of freedom
## Residual deviance: 514.64 on 397 degrees of freedom
## AIC: 1074.6
##
## Number of Fisher Scoring iterations: 6IRR(poisson22b)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.2058383 1.197246 0.1446393 0.2929314
## sex 0.6559620 1.134094 0.5125861 0.8394417brv %>%
select(sex, age_entry) %>%
group_by(sex) %>%
summarise(meanAgeAtEntry = mean(age_entry))## Source: local data frame [2 x 2]
##
## sex meanAgeAtEntry
## 1 1 79.07153
## 2 2 78.65995(c) Breaking records into pre and post bereavement. In these data a subject changes exposure status from not bereaved to bereaved when his or her spouse dies. The first stage of the analysis therefore is to partition each follow–up into a record describing the period of follow-up pre–bereavement and (for subjects who were bereaved during the study) the period post–bereavement.
## Creating times relativ to spouse death (year=0)
brv2 <- mutate(brv_raw,
y_before_sp_dth = as.numeric(doe -dosp) / 365.24,
y_after_sp_dth = as.numeric(dox - dosp) / 365.24)
## Splitting at spouse death (year=0)
brvSplit <- survSplit(brv2, cut = 0, end="y_after_sp_dth", start="y_before_sp_dth", id="id",event="fail")
## Calculating risk times
brvSplit <- mutate(brvSplit,
t_sp_at_risk = y_after_sp_dth - y_before_sp_dth,
brv = ifelse(y_after_sp_dth > 0, 1, 0))
## Look at the five first couples
brvSplit %>% select(couple, id, sex, doe, dosp, dox, fail, y_before_sp_dth, y_after_sp_dth, t_sp_at_risk) %>% filter(couple<=5) %>% arrange(couple, id)## couple id sex doe dosp dox fail y_before_sp_dth
## 1 1 197 1 1981-01-15 2000-01-01 1983-04-18 1 -18.9601358
## 2 2 177 1 1981-04-27 1984-11-20 1983-01-24 1 -3.5675172
## 3 2 270 2 1981-04-27 1983-01-24 1984-11-20 0 -1.7440587
## 4 2 270 2 1981-04-27 1983-01-24 1984-11-20 1 0.0000000
## 5 3 168 1 1981-01-20 1981-12-31 1981-08-03 1 -0.9445844
## 6 3 384 2 1981-01-20 1981-08-03 1981-12-31 0 -0.5338955
## 7 3 384 2 1981-01-20 1981-08-03 1981-12-31 1 0.0000000
## 8 4 12 1 1981-01-20 1988-11-23 1991-01-01 0 -7.8414193
## 9 4 12 1 1981-01-20 1988-11-23 1991-01-01 0 0.0000000
## 10 4 300 2 1981-01-20 2000-01-01 1988-11-23 1 -18.9464462
## 11 5 212 1 1981-05-07 1984-02-04 1984-02-04 1 -2.7461395
## y_after_sp_dth t_sp_at_risk
## 1 -16.7068229 2.2533129
## 2 -1.8234585 1.7440587
## 3 0.0000000 1.7440587
## 4 1.8234585 1.8234585
## 5 -0.4106889 0.5338955
## 6 0.0000000 0.5338955
## 7 0.4106889 0.4106889
## 8 0.0000000 7.8414193
## 9 2.1054649 2.1054649
## 10 -11.1050268 7.8414193
## 11 0.0000000 2.7461395(d) Now find the (crude) effect of bereavement.
poisson22d <- glm( fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson, data=brvSplit)
summary(poisson22d)##
## Call:
## glm(formula = fail ~ brv + offset(log(t_sp_at_risk)), family = poisson,
## data = brvSplit)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5442 -1.0377 -0.0657 0.8531 3.4313
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.20379 0.07125 -30.932 <2e-16 ***
## brv 0.11970 0.13199 0.907 0.364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 704.42 on 554 degrees of freedom
## Residual deviance: 703.61 on 553 degrees of freedom
## AIC: 1263.6
##
## Number of Fisher Scoring iterations: 6IRR(poisson22d)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.1103837 1.073846 0.09599711 0.1269263
## brv 1.1271537 1.141098 0.87022190 1.4599443(e) Since there is a strong possibility that the effect of bereavement is not the same for men as for women, use streg to estimate the effect of bereavement separately for men and women. Do this both by fitting separate models for males and females (e.g. streg brv if sex==1) as well as by using a single model with an interaction term (you may need to create dummy variables). Confirm that the estimates are identical for these two approaches.
## Poisson regression for sex==1
poisson22e1 <- glm( fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson, data=filter(brvSplit, sex==1))
summary(poisson22e1)##
## Call:
## glm(formula = fail ~ brv + offset(log(t_sp_at_risk)), family = poisson,
## data = filter(brvSplit, sex == 1))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6325 -1.0016 0.1241 0.8594 3.1501
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.00434 0.08248 -24.301 <2e-16 ***
## brv 0.01080 0.19030 0.057 0.955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 341.95 on 294 degrees of freedom
## Residual deviance: 341.95 on 293 degrees of freedom
## AIC: 707.95
##
## Number of Fisher Scoring iterations: 6IRR(poisson22e1)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.1347495 1.085975 0.1146358 0.1583925
## brv 1.0108630 1.209614 0.6961532 1.4678436## Poisson regression for sex==2
poisson22e2 <- glm( fail ~ brv + offset( log(t_sp_at_risk) ), family=poisson, data=filter(brvSplit, sex==2))
summary(poisson22e2)##
## Call:
## glm(formula = fail ~ brv + offset(log(t_sp_at_risk)), family = poisson,
## data = filter(brvSplit, sex == 2))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4851 -0.9557 -0.5097 0.8569 3.4490
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6301 0.1414 -18.598 <2e-16 ***
## brv 0.4853 0.2032 2.389 0.0169 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 350.83 on 259 degrees of freedom
## Residual deviance: 345.21 on 258 degrees of freedom
## AIC: 543.21
##
## Number of Fisher Scoring iterations: 6IRR(poisson22e2)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.07206987 1.151909 0.05462283 0.09508965
## brv 1.62461324 1.225271 1.09097543 2.41927370## Poisson regression, interaction with sex
brvSplit2 <- mutate(brvSplit,
sex = as.factor(sex),
brv = as.factor(brv))
poisson22e3 <- glm( fail ~ sex + brv:sex + offset( log(t_sp_at_risk) ), family=poisson, data=brvSplit2)
summary(poisson22e3)##
## Call:
## glm(formula = fail ~ sex + brv:sex + offset(log(t_sp_at_risk)),
## family = poisson, data = brvSplit2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6325 -0.9659 -0.2355 0.8583 3.4490
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.00434 0.08248 -24.301 < 2e-16 ***
## sex2 -0.62578 0.16371 -3.822 0.000132 ***
## sex1:brv1 0.01080 0.19030 0.057 0.954724
## sex2:brv1 0.48527 0.20316 2.389 0.016913 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 704.42 on 554 degrees of freedom
## Residual deviance: 687.16 on 551 degrees of freedom
## AIC: 1251.2
##
## Number of Fisher Scoring iterations: 6IRR(poisson22e3)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.1347495 1.085975 0.1146358 0.1583925
## sex2 0.5348431 1.177878 0.3880340 0.7371962
## sex1:brv1 1.0108630 1.209614 0.6961532 1.4678436
## sex2:brv1 1.6246132 1.225271 1.0909754 2.4192737(f) Controlling for age. There is strong confounding by age. Expand the data by 5 year age–bands, and check that the rate is increasing with age. Find the effect of bereavement controlled for age. If you wish to study the distribution of age then it is useful to know that age at entry and exit. Translate your time scale to age.
## Translate time scale from years from spouse death to ages
brvSplit3 <- brvSplit2 %>%
mutate(age_sp_dth = as.numeric(dosp - dob) / 365.24, # Age at spouse death
age_start = age_sp_dth + y_before_sp_dth, # Age at start of timeband
age_end = age_sp_dth + y_after_sp_dth) # Age at end of timeband
age_cat <- seq(70,100,5) # Split at these ages
brvSplit4 <- survSplit(brvSplit3, cut=age_cat, start="age_start", end="age_end", event="fail", zero = 0)
brvSplit4 <- mutate(brvSplit4,
t_at_risk = age_end- age_start, # Creating new time at risk
age = cut(age_end, age_cat)) # Creating age band category
## Calculate crude rates
brvSplit4 %>%
select(fail, age, t_at_risk) %>%
group_by(age) %>%
summarise(D = sum(fail), Y = sum(t_at_risk), Rate = D/Y) %>%
mutate(CI_low = Rate + qnorm(0.025) * Rate / sqrt(D),
CI_high = Rate + qnorm(0.975) * Rate / sqrt(D))## Source: local data frame [5 x 6]
##
## age D Y Rate CI_low CI_high
## 1 (75,80] 45 703.612419 0.06395566 0.04526947 0.08264186
## 2 (80,85] 123 1184.684043 0.10382515 0.08547676 0.12217355
## 3 (85,90] 95 490.021356 0.19386910 0.15488434 0.23285386
## 4 (90,95] 12 55.090352 0.21782399 0.09458073 0.34106724
## 5 (95,100] 3 2.299858 1.30442857 -0.17164419 2.78050133poisson22f1 <- glm( fail ~ brv + age + offset( log(t_at_risk) ), family=poisson, data=brvSplit4)
summary(poisson22f1)##
## Call:
## glm(formula = fail ~ brv + age + offset(log(t_at_risk)), family = poisson,
## data = brvSplit4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8226 -0.7876 -0.5159 0.4976 3.3740
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7359 0.1495 -18.299 < 2e-16 ***
## brv1 -0.1515 0.1371 -1.105 0.26928
## age(80,85] 0.5106 0.1757 2.907 0.00365 **
## age(85,90] 1.1627 0.1869 6.220 4.98e-10 ***
## age(90,95] 1.2847 0.3290 3.905 9.42e-05 ***
## age(95,100] 3.0431 0.5967 5.100 3.40e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1131.0 on 1035 degrees of freedom
## Residual deviance: 1074.4 on 1030 degrees of freedom
## AIC: 1642.4
##
## Number of Fisher Scoring iterations: 6IRR(poisson22f1)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.06483296 1.161264 0.04836499 0.08690817
## brv1 0.85941225 1.146995 0.65684160 1.12445589
## age(80,85] 1.66632969 1.192030 1.18096760 2.35116919
## age(85,90] 3.19848093 1.205546 2.21729440 4.61385743
## age(90,95] 3.61371336 1.389570 1.89630372 6.88651512
## age(95,100] 20.97060543 1.816161 6.51131081 67.53882670poisson22f2 <- glm( fail ~ brv +age + sex + offset( log(t_at_risk) ), family=poisson, data=brvSplit4)
summary(poisson22f2)##
## Call:
## glm(formula = fail ~ brv + age + sex + offset(log(t_at_risk)),
## family = poisson, data = brvSplit4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7766 -0.7973 -0.4988 0.4702 3.3896
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.57737 0.15364 -16.776 < 2e-16 ***
## brv1 -0.02676 0.14019 -0.191 0.848603
## age(80,85] 0.51641 0.17567 2.940 0.003286 **
## age(85,90] 1.15434 0.18626 6.197 5.74e-10 ***
## age(90,95] 1.29672 0.32900 3.941 8.10e-05 ***
## age(95,100] 3.32531 0.60227 5.521 3.37e-08 ***
## sex2 -0.49188 0.13054 -3.768 0.000165 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1131.0 on 1035 degrees of freedom
## Residual deviance: 1059.6 on 1029 degrees of freedom
## AIC: 1629.6
##
## Number of Fisher Scoring iterations: 6IRR(poisson22f2)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.07597334 1.166067 0.05621898 0.1026690
## brv1 0.97359228 1.150493 0.73968081 1.2814743
## age(80,85] 1.67599651 1.192045 1.18778920 2.3648677
## age(85,90] 3.17193748 1.204741 2.20177467 4.5695808
## age(90,95] 3.65729031 1.389577 1.91915306 6.9696225
## age(95,100] 27.80767040 1.826262 8.54084990 90.5374221
## sex2 0.61147403 1.139443 0.47343510 0.7897608(g) Now estimate the effect of bereavement (controlled for age) separately for each sex.
poisson22g <- glm( fail ~ age + sex + brv:sex + offset( log(t_at_risk) ), family=poisson, data=brvSplit4)
summary(poisson22g)##
## Call:
## glm(formula = fail ~ age + sex + brv:sex + offset(log(t_at_risk)),
## family = poisson, data = brvSplit4)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7246 -0.7874 -0.4931 0.4611 3.3591
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.5398 0.1554 -16.343 < 2e-16 ***
## age(80,85] 0.5176 0.1757 2.946 0.003221 **
## age(85,90] 1.1410 0.1866 6.113 9.76e-10 ***
## age(90,95] 1.2962 0.3291 3.939 8.19e-05 ***
## age(95,100] 3.3586 0.6031 5.568 2.57e-08 ***
## sex2 -0.6221 0.1656 -3.756 0.000172 ***
## sex1:brv1 -0.1940 0.1925 -1.008 0.313628
## sex2:brv1 0.1823 0.2085 0.874 0.381981
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1131.0 on 1035 degrees of freedom
## Residual deviance: 1057.8 on 1028 degrees of freedom
## AIC: 1629.8
##
## Number of Fisher Scoring iterations: 6IRR(poisson22g)## IRR Std. err CI_lower CI_upper
## (Intercept) 0.07888104 1.168132 0.05816856 0.1069688
## age(80,85] 1.67794302 1.192075 1.18911036 2.3677304
## age(85,90] 3.12991468 1.205198 2.17099247 4.5123906
## age(90,95] 3.65549666 1.389698 1.91788439 6.9673938
## age(95,100] 28.74862893 1.827859 8.81474206 93.7615259
## sex2 0.53681350 1.180123 0.38801249 0.7426790
## sex1:brv1 0.82368705 1.212270 0.56481703 1.2012038
## sex2:brv1 1.19991663 1.231788 0.79744516 1.8055159(h) We have assumed that any effect of bereavement is both immediate and permanent. This is not realistic and we might wish to improve the analysis by further subdividing the post–bereavement follow–up. How might you do this? (you are not expected to actually do it)
(i) Analysis using Cox regression. We can also model these data using Cox regression. Provided we use the attained age as the time scale and split the data to obtain separate observations for the bereaved and non-bereaved person-time the following command will estimate the effect of bereavement adjusted for attained age.
summary(coxph(Surv(age_start, age_end, fail) ~ brv,
method=c("breslow"), data = brvSplit4))## Call:
## coxph(formula = Surv(age_start, age_end, fail) ~ brv, data = brvSplit4,
## method = c("breslow"))
##
## n= 1036, number of events= 278
##
## coef exp(coef) se(coef) z Pr(>|z|)
## brv1 -0.2062 0.8136 0.1390 -1.483 0.138
##
## exp(coef) exp(-coef) lower .95 upper .95
## brv1 0.8136 1.229 0.6196 1.069
##
## Concordance= 0.511 (se = 0.014 )
## Rsquare= 0.002 (max possible= 0.93 )
## Likelihood ratio test= 2.25 on 1 df, p=0.1337
## Wald test = 2.2 on 1 df, p=0.138
## Score (logrank) test = 2.21 on 1 df, p=0.1374summary(coxph(Surv(age_start, age_end, fail) ~ brv + sex,
method=c("breslow"), data = brvSplit4))## Call:
## coxph(formula = Surv(age_start, age_end, fail) ~ brv + sex, data = brvSplit4,
## method = c("breslow"))
##
## n= 1036, number of events= 278
##
## coef exp(coef) se(coef) z Pr(>|z|)
## brv1 -0.07782 0.92513 0.14244 -0.546 0.584814
## sex2 -0.47273 0.62330 0.13074 -3.616 0.000299 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## brv1 0.9251 1.081 0.6998 1.2231
## sex2 0.6233 1.604 0.4824 0.8053
##
## Concordance= 0.56 (se = 0.018 )
## Rsquare= 0.015 (max possible= 0.93 )
## Likelihood ratio test= 15.82 on 2 df, p=0.0003662
## Wald test = 15.19 on 2 df, p=0.0005036
## Score (logrank) test = 15.49 on 2 df, p=0.0004338That is, we do not have to split the data by attained age (although we can fit the model to data split by attained age and the results will be the same).
(j) Use the Cox model to estimate the effect of bereavement separately for males and females and compare the estimates to those obtained using Poisson regression.
summary(coxph(Surv(age_start, age_end, fail) ~ sex + sex:brv,
method=c("breslow"), data = brvSplit4))## Call:
## coxph(formula = Surv(age_start, age_end, fail) ~ sex + sex:brv,
## data = brvSplit4, method = c("breslow"))
##
## n= 1036, number of events= 278
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex2 -0.58154 0.55904 0.16556 -3.513 0.000444 ***
## sex1:brv1 -0.21638 0.80543 0.19302 -1.121 0.262274
## sex2:brv1 0.09877 1.10381 0.21190 0.466 0.641143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex2 0.5590 1.789 0.4041 0.7733
## sex1:brv1 0.8054 1.242 0.5517 1.1758
## sex2:brv1 1.1038 0.906 0.7287 1.6721
##
## Concordance= 0.568 (se = 0.018 )
## Rsquare= 0.016 (max possible= 0.93 )
## Likelihood ratio test= 17.09 on 3 df, p=0.0006759
## Wald test = 16.52 on 3 df, p=0.0008852
## Score (logrank) test = 16.9 on 3 df, p=0.0007426