DEM 7223 - Event History Analysis - Models of Frailty
Corey Sparks, PhD
November 6, 2017
This example will illustrate how to fit the extended Cox Proportional hazards model with Gaussian frailty to continuous duration data (i.e. person-level data) and a discrete-time (longitudinal) data set. In this example, I will use the event of a child dying before age 5. The data for this example come from the model.data Demographic and Health Survey for 2012 birth history recode file. This file contains information for all births to women in the survey.
The longitudinal data example uses data from the ECLS-K. Specifically, we will examine the transition into poverty between kindergarten and third grade.
#Load required libraries
library(foreign)
library(survival)
library(car)
library(survey)
library(coxme)
library(knitr)
library(lme4)
opts_chunk$set(tidy.opts=list(width.cutoff=60),tidy=TRUE)# load the data
model.dat <- haven::read_dta("~/Google Drive/zzbr62dt/ZZBR62FL.DTA")# We form a subset of variables
sub <- data.frame(CASEID = model.dat$caseid, v008 = model.dat$v008,
bord = model.dat$bidx, csex = model.dat$b4, b2 = model.dat$b2,
b3 = model.dat$b3, b5 = model.dat$b5, b7 = model.dat$b7,
ibint = model.dat$b11, rural = model.dat$v025, educ = model.dat$v106,
age = model.dat$v012, partneredu = model.dat$v701, partnerage = model.dat$v730,
hhses = model.dat$v190, weight = model.dat$v005/1e+06, psu = model.dat$v021,
strata = model.dat$v022, region = model.dat$v023)
sub$death.age <- ifelse(sub$b5 == 1, ((((sub$v008)) + 1900) -
(((sub$b3)) + 1900)), sub$b7)
# censoring indicator for death by age 5, in months (<=60
# months)
sub$d.event <- ifelse(is.na(sub$b7) == T | sub$b7 > 60, 0, 1)
sub$d.eventfac <- factor(sub$d.event)
levels(sub$d.eventfac) <- c("Alive at Age 5", "Dead by Age 5")
table(sub$d.eventfac)##
## Alive at Age 5 Dead by Age 5
## 19224 4442# recodes
sub$male <- ifelse(sub$csex == 1, 1, 0)
sub$educ.high <- ifelse(sub$educ %in% c(2, 3), 1, 0)
sub$age2 <- sub$age^2
sub$partnerhiedu <- ifelse(sub$partneredu < 3, 0, ifelse(sub$partneredu %in%
c(8, 9), NA, 1))Fit the ordinary Cox model
Here I fit the ordinary Cox model without frailty, just for comparison sake.
# using coxph in survival library
fit.cox2 <- coxph(Surv(death.age, d.event) ~ bord + male + educ.high +
I(age/5) + I(hhses > 3), data = sub, weights = weight)
summary(fit.cox2)## Call:
## coxph(formula = Surv(death.age, d.event) ~ bord + male + educ.high +
## I(age/5) + I(hhses > 3), data = sub, weights = weight)
##
## n= 23666, number of events= 4442
##
## coef exp(coef) se(coef) z Pr(>|z|)
## bord 0.165766 1.180297 0.007043 23.536 < 2e-16 ***
## male 0.075614 1.078547 0.030602 2.471 0.0135 *
## educ.high -0.043896 0.957053 0.053125 -0.826 0.4086
## I(age/5) -0.057113 0.944488 0.010958 -5.212 1.87e-07 ***
## I(hhses > 3)TRUE -0.184026 0.831914 0.035336 -5.208 1.91e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## bord 1.1803 0.8472 1.1641 1.1967
## male 1.0785 0.9272 1.0158 1.1452
## educ.high 0.9571 1.0449 0.8624 1.0621
## I(age/5) 0.9445 1.0588 0.9244 0.9650
## I(hhses > 3)TRUE 0.8319 1.2020 0.7762 0.8916
##
## Concordance= 0.607 (se = 0.005 )
## Rsquare= 0.026 (max possible= 0.972 )
## Likelihood ratio test= 623.2 on 5 df, p=0
## Wald test = 679.2 on 5 df, p=0
## Score (logrank) test = 690.1 on 5 df, p=0plot(survfit(fit.cox2), ylim = c(0.8, 1), xlim = c(0, 60), ylab = "S(t)",
xlab = "Age in Months")
title(main = "Survival Function for Child Mortality")
Fit the Cox model with frailty at the regional level
The coxme() function in the coxme library Link will fit the Cox model with shared frailty, assuming a Gaussian frailty term. The model would look like:
\(h_j (t) = h_{0j} e^{(x'\beta + u_j)}\)
where
\(u_j \sim N(0, \sigma^2)\)
is a Normally distributed random effect, identical for each person in the jth group. This term raises or lowers the average hazard function the same way for each person within each group, but allows the overall risk for people in different groups to be different. This would be considered to be a random intercept model, if we were considering an ordinary linear or genearlized linear model.
fit.cox.f <- coxme(Surv(death.age, d.event) ~ bord + male + educ.high +
I(age/5) + I(hhses > 3) + (1 | region), data = sub, weights = weight)
summary(fit.cox.f)## Cox mixed-effects model fit by maximum likelihood
## Data: sub
## events, n = 4442, 23666
## Iterations= 10 53
## NULL Integrated Fitted
## Log-likelihood -42496.25 -42168.52 -42160.02
##
## Chisq df p AIC BIC
## Integrated loglik 655.46 6.00 0 643.46 605.06
## Penalized loglik 672.45 10.29 0 651.87 586.03
##
## Model: Surv(death.age, d.event) ~ bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 | region)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## bord 0.16407385 1.1783013 0.00704951 23.27 0.0e+00
## male 0.07622015 1.0792001 0.03060558 2.49 1.3e-02
## educ.high -0.05391452 0.9475131 0.05376093 -1.00 3.2e-01
## I(age/5) -0.05596300 0.9455741 0.01096637 -5.10 3.3e-07
## I(hhses > 3)TRUE -0.19942922 0.8191982 0.04246931 -4.70 2.7e-06
##
## Random effects
## Group Variable Std Dev Variance
## region Intercept 0.10319525 0.01064926This gives us the variance in child mortality by region, which is honestly pretty substantial. We can use a likelihood ratio test to see if the frailty model is significantly better than the ordinary Cox model:
anova(fit.cox.f, fit.cox2)## Analysis of Deviance Table
## Cox model: response is Surv(death.age, d.event)
## Model 1: ~bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 | region)
## Model 2: ~bord + male + educ.high + I(age/5) + I(hhses > 3)
## loglik Chisq Df P(>|Chi|)
## 1 -42169
## 2 -42185 32.242 1 1.361e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(fit.cox.f)## [1] 84340.63AIC(fit.cox2)## [1] 84379.28Which it is, and this is supported by the AIC difference between the two models of 38.66 points.
So, what are the frailties in this case? We can get those from the frail portion of the model structure:
fit.cox.f$frail## $region
## 1 2 3 4 5
## -0.060602442 -0.160882759 0.001678816 0.013684448 -0.048005072
## 6 7 8
## 0.074802258 0.133004772 0.046319979hist(fit.cox.f$frail$region)
Which shows the region region.7 has the highest frailty, which means that the average level of childhood mortality is highest in that region, while the lowest frailty is in region.2 .
Random slopes
If we were interested in whether a predictor variable had heterogeneous effects across the various groups within our data, we could include that in our model as well, and we would have effectively a random slope model:
\(h_j (t) = h_{0j} e^{(x'\beta + u_j+\gamma_j 'x)}\)
where \(\gamma_j\) is a group-specific effect of a particular predictor variable, and these two random effects will be distributed as:
\[ \left[\begin{array} {rrr} u_j \\ \gamma_j \end{array}\right] \sim \text{MVN}(0, \Sigma) \]
# See if higher birth order children face equal disadvantage
# in all regions
fit.cox.f2 <- coxme(Surv(death.age, d.event) ~ bord + male +
educ.high + I(age/5) + I(hhses > 3) + (1 + bord | region),
data = sub, weights = weight)
summary(fit.cox.f2)## Cox mixed-effects model fit by maximum likelihood
## Data: sub
## events, n = 4442, 23666
## Iterations= 10 53
## NULL Integrated Fitted
## Log-likelihood -42496.25 -42168.24 -42156.91
##
## Chisq df p AIC BIC
## Integrated loglik 656.01 8.00 0 640.01 588.82
## Penalized loglik 678.68 12.38 0 653.92 574.73
##
## Model: Surv(death.age, d.event) ~ bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 + bord | region)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## bord 0.16344882 1.1775651 0.008861874 18.44 0.0e+00
## male 0.07631351 1.0793009 0.030611447 2.49 1.3e-02
## educ.high -0.05409460 0.9473425 0.053852298 -1.00 3.2e-01
## I(age/5) -0.05580822 0.9457205 0.010967384 -5.09 3.6e-07
## I(hhses > 3)TRUE -0.19967351 0.8189981 0.042493875 -4.70 2.6e-06
##
## Random effects
## Group Variable Std Dev Variance Corr
## region Intercept 0.1162372113 0.0135110893 -0.4586696457
## bord 0.0129927062 0.0001688104anova(fit.cox.f, fit.cox.f2)## Analysis of Deviance Table
## Cox model: response is Surv(death.age, d.event)
## Model 1: ~bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 | region)
## Model 2: ~bord + male + educ.high + I(age/5) + I(hhses > 3) + (1 + bord | region)
## loglik Chisq Df P(>|Chi|)
## 1 -42169
## 2 -42168 0.559 2 0.7562And it looks like there is not significant variation in this effect, because the model with the additional term does not fit the data significantly better than the model with only the random “intercept”.
Using Longitudinal Data
As in the other examples, I illustrate fitting these models to data that are longitudinal, instead of person-duration. In this example, we will examine how to fit the Cox model to a longitudinally collected data set.
First we load our data First we load our data
load("~/Google Drive/classes/dem7223/class17/data/eclsk.Rdata")
names(eclsk) <- tolower(names(eclsk))
# get out only the variables I'm going to use for this
# example
myvars <- c("childid", "gender", "race", "r1_kage", "r3age",
"r4age", "r5age", "r6age", "r7age", "c1r4mtsc", "c4r4mtsc",
"c5r4mtsc", "c6r4mtsc", "c7r4mtsc", "wkpov_r", "w1povrty",
"w3povrty", "w5povrty", "w8povrty", "wkmomed", "w1momed",
"w3momed", "w5momed", "w8momed", "s2_id", "c1_7fp0", "c17fpstr",
"c17fppsu")
eclsk <- eclsk[, myvars]
eclsk$age1 <- ifelse(eclsk$r1_kage == -9, NA, eclsk$r1_kage/12)
eclsk$age2 <- ifelse(eclsk$r3age == -9, NA, eclsk$r3age/12)
# for the later waves, the NCES group the ages into ranges of
# months, so 1= <105 months, 2=105 to 108 months. So, I fix
# the age at the midpoint of the interval they give, and make
# it into years by dividing by 12
eclsk$age3 <- recode(eclsk$r5age, recodes = "1=105; 2=107; 3=109; 4=112; 5=115; 6=117; -9=NA")/12
eclsk$age4 <- recode(eclsk$r7age, recodes = "1=155; 2=166; 3=172; 4=178; 5=192; -9=NA")/12
eclsk$pov1 <- ifelse(eclsk$wkpov_r == 1, 1, 0)
eclsk$pov2 <- ifelse(eclsk$w3povrty == 1, 1, 0)
eclsk$pov3 <- ifelse(eclsk$w5povrty == 1, 1, 0)
eclsk$pov4 <- ifelse(eclsk$w8povrty == 1, 1, 0)
# mother's education-time varying > hs ==1
eclsk$momedu1 <- recode(eclsk$wkmomed, recodes = "1:2='primarysch';3='highschool';4:5='jcollege'; 6:9='college'; else =NA",
as.factor.result = T)
eclsk$momedu1 <- relevel(eclsk$momedu1, ref = "highschool")
eclsk$momedu2 <- recode(eclsk$w3momed, recodes = "1:2='primarysch';3='highschool';4:5='jcollege'; 6:9='college'; else =NA",
as.factor.result = T)
eclsk$momedu2 <- relevel(eclsk$momedu2, ref = "highschool")
eclsk$momedu3 <- recode(eclsk$w5momed, recodes = "1:2='primarysch';3='highschool';4:5='jcollege'; 6:9='college'; else =NA",
as.factor.result = T)
eclsk$momedu3 <- relevel(eclsk$momedu3, ref = "highschool")
# Recode race with white, non Hispanic as reference using
# dummy vars
eclsk$hisp <- recode(eclsk$race, recodes = "3:4=1;-9=NA; else=0")
eclsk$black <- recode(eclsk$race, recodes = "2=1;-9=NA; else=0")
eclsk$asian <- recode(eclsk$race, recodes = "5=1;-9=NA; else=0")
eclsk$nahn <- recode(eclsk$race, recodes = "6:7=1;-9=NA; else=0")
eclsk$other <- recode(eclsk$race, recodes = "8=1;-9=NA; else=0")
eclsk$race_gr <- recode(eclsk$race, recodes = "3:4='hisp'; 2='nh black'; 5='nh asian'; 6:7='nahn'; 8='other'; 1='nh white'; else=NA",
as.factor.result = T)
eclsk$race_gr <- relevel(eclsk$race_gr, ref = "nh white")
eclsk$male <- recode(eclsk$gender, recodes = "1=1; 2=0; -9=NA")Now, I need to form the transition variable, this is my event variable, and in this case it will be 1 if a child enters poverty between the first wave of the data and the third grade wave, and 0 otherwise. NOTE I need to remove any children who are already in poverty age wave 1, because they are not at risk of experiencing this particular transition.
eclsk <- subset(eclsk, is.na(pov1) == F & is.na(pov2) == F &
is.na(pov3) == F & is.na(age1) == F & is.na(age2) == F &
is.na(age3) == F & pov1 != 1 & is.na(eclsk$c1_7fp0) == F &
is.na(eclsk$s2_id) == F)# make an id that is the combination of state and strata
eclsk$sampleid <- paste(eclsk$c17fpstr, eclsk$c17fppsu)
# within each sampling unit, sum the weights
wts <- tapply(eclsk$c1_7fp0, eclsk$sampleid, sum)
# make a data frame from this
wts <- data.frame(id = names(unlist(wts)), wt = unlist(wts))
# get the unique sampling location ids'
t1 <- as.data.frame(table(eclsk$sampleid))
# put all of this into a data set
wts2 <- data.frame(ids = wts$id, sumwt = wts$wt, jn = t1$Freq)
# merge all of this back to the original data file
eclsk2 <- merge(eclsk, wts2, by.x = "sampleid", by.y = "ids",
all.x = T)
# In the new data set, multiply the original weight by the
# fraction of the sampling unit total population each person
# represents
eclsk2$swts <- eclsk2$c1_7fp0 * (eclsk2$jn/eclsk2$sumwt)Now we do the entire data set. To analyze data longitudinally, we need to reshape the data from the current “wide” format (repeated measures in columns) to a “long” format (repeated observations in rows). The reshape() function allows us to do this easily. It allows us to specify our repeated measures, time varying covariates as well as time-constant covariates.
e.long <- reshape(eclsk2, idvar = "childid", varying = list(age = c("age1",
"age2", "age3"), age2 = c("age2", "age3", "age4"), momedu = c("momedu1",
"momedu2", "momedu3")), times = 1:3, direction = "long",
v.names = c("age_enter", "age_exit", "momedu"))
e.long <- e.long[order(e.long$childid, e.long$time), ]
head(e.long)## sampleid childid gender race r1_kage r3age r4age r5age r6age
## 0015006C.1 77 1 0015006C 1 1 76.8 87.83 95.53 6 NA
## 0015006C.2 77 1 0015006C 1 1 76.8 87.83 95.53 6 NA
## 0015006C.3 77 1 0015006C 1 1 76.8 87.83 95.53 6 NA
## 0015014C.1 77 1 0015014C 2 1 64.2 75.23 82.93 2 2
## 0015014C.2 77 1 0015014C 2 1 64.2 75.23 82.93 2 2
## 0015014C.3 77 1 0015014C 2 1 64.2 75.23 82.93 2 2
## r7age c1r4mtsc c4r4mtsc c5r4mtsc c6r4mtsc c7r4mtsc wkpov_r
## 0015006C.1 4 56.187 51.941 52.684 NA 50.263 2
## 0015006C.2 4 56.187 51.941 52.684 NA 50.263 2
## 0015006C.3 4 56.187 51.941 52.684 NA 50.263 2
## 0015014C.1 2 57.739 55.178 51.977 59.293 60.001 2
## 0015014C.2 2 57.739 55.178 51.977 59.293 60.001 2
## 0015014C.3 2 57.739 55.178 51.977 59.293 60.001 2
## w1povrty w3povrty w5povrty w8povrty wkmomed w1momed w3momed
## 0015006C.1 2 2 2 2 5 5 6
## 0015006C.2 2 2 2 2 5 5 6
## 0015006C.3 2 2 2 2 5 5 6
## 0015014C.1 2 2 2 2 5 6 6
## 0015014C.2 2 2 2 2 5 6 6
## 0015014C.3 2 2 2 2 5 6 6
## w5momed w8momed s2_id c1_7fp0 c17fpstr c17fppsu pov1 pov2 pov3
## 0015006C.1 8 8 0015 934.11 77 1 0 0 0
## 0015006C.2 8 8 0015 934.11 77 1 0 0 0
## 0015006C.3 8 8 0015 934.11 77 1 0 0 0
## 0015014C.1 6 6 0015 129.55 77 1 0 0 0
## 0015014C.2 6 6 0015 129.55 77 1 0 0 0
## 0015014C.3 6 6 0015 129.55 77 1 0 0 0
## pov4 hisp black asian nahn other race_gr male sumwt jn
## 0015006C.1 0 0 0 0 0 0 nh white 1 13743.12 58
## 0015006C.2 0 0 0 0 0 0 nh white 1 13743.12 58
## 0015006C.3 0 0 0 0 0 0 nh white 1 13743.12 58
## 0015014C.1 0 0 0 0 0 0 nh white 0 13743.12 58
## 0015014C.2 0 0 0 0 0 0 nh white 0 13743.12 58
## 0015014C.3 0 0 0 0 0 0 nh white 0 13743.12 58
## swts time age_enter age_exit momedu
## 0015006C.1 3.942218 1 6.400000 7.319167 jcollege
## 0015006C.2 3.942218 2 7.319167 9.750000 college
## 0015006C.3 3.942218 3 9.750000 14.833333 college
## 0015014C.1 0.546739 1 5.350000 6.269167 jcollege
## 0015014C.2 0.546739 2 6.269167 8.916667 college
## 0015014C.3 0.546739 3 8.916667 13.833333 collegeSince I have 3 risk periods this time, the previous way of coding the transitions has become burdensome. So, instead I hard code it. To do this, I initialize the povtran variable as missing, then hard code the specific transitions I’m interested in. This lets me identify when values between waves stay the same, versus when they change. Also, since I’m initilizing the variable as missing, I dont’ have to recode the missing values as I’ve done in previous examples.
# find which kids failed in the first time period and remove
# them from the second risk period risk set
e.long$povtran <- NA
e.long$povtran[e.long$pov1 == 0 & e.long$pov2 == 1 & e.long$time ==
1] <- 1
e.long$povtran[e.long$pov2 == 0 & e.long$pov3 == 1 & e.long$time ==
2] <- 1
e.long$povtran[e.long$pov3 == 0 & e.long$pov4 == 1 & e.long$time ==
3] <- 1
e.long$povtran[e.long$pov1 == 0 & e.long$pov2 == 0 & e.long$time ==
1] <- 0
e.long$povtran[e.long$pov2 == 0 & e.long$pov3 == 0 & e.long$time ==
2] <- 0
e.long$povtran[e.long$pov3 == 0 & e.long$pov4 == 0 & e.long$time ==
3] <- 0
# find which kids failed in earlier time periods and remove
# them from the second & third period risk set
failed1 <- which(is.na(e.long$povtran) == T)
e.long <- e.long[-failed1, ]
# round age to the nearest year
e.long$age_enter_r <- round(e.long$age_enter, 0)
e.long$age_exit_r <- round(e.long$age_exit, 0)
head(e.long[, c("childid", "time", "age_enter", "age_exit", "pov1",
"pov2", "pov3", "time", "povtran")], n = 10)## childid time age_enter age_exit pov1 pov2 pov3 time.1 povtran
## 0015006C.1 0015006C 1 6.400000 7.319167 0 0 0 1 0
## 0015006C.2 0015006C 2 7.319167 9.750000 0 0 0 2 0
## 0015006C.3 0015006C 3 9.750000 14.833333 0 0 0 3 0
## 0015014C.1 0015014C 1 5.350000 6.269167 0 0 0 1 0
## 0015014C.2 0015014C 2 6.269167 8.916667 0 0 0 2 0
## 0015014C.3 0015014C 3 8.916667 13.833333 0 0 0 3 0
## 0015019C.1 0015019C 1 5.685833 6.605833 0 0 0 1 0
## 0015019C.2 0015019C 2 6.605833 9.333333 0 0 0 2 0
## 0015019C.3 0015019C 3 9.333333 14.333333 0 0 0 3 0
## 0015020C.1 0015020C 1 6.583333 7.802500 0 0 0 1 0Fit basic Cox model
Now we fit the Cox model using full survey design. In the ECLS-K, I use the longitudinal weight for waves 1-7, as well as the associated psu and strata id’s for the longitudinal data from these waves from the parents of the child, since no data from the child themselves are used in the outcome.
# Fit the model
fitl1 <- coxph(Surv(time = time, event = povtran) ~ momedu +
race_gr, data = e.long, weights = swts)
summary(fitl1)## Call:
## coxph(formula = Surv(time = time, event = povtran) ~ momedu +
## race_gr, data = e.long, weights = swts)
##
## n= 4757, number of events= 166
## (1084 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## momeducollege -2.77946 0.06207 0.34086 -8.154 3.33e-16 ***
## momedujcollege -0.88433 0.41299 0.16636 -5.316 1.06e-07 ***
## momeduprimarysch 0.64833 1.91234 0.22085 2.936 0.003328 **
## race_grhisp 0.64218 1.90062 0.18921 3.394 0.000689 ***
## race_grnahn 1.01589 2.76183 0.31121 3.264 0.001097 **
## race_grnh asian 0.32817 1.38842 0.41285 0.795 0.426678
## race_grnh black 0.87163 2.39079 0.21311 4.090 4.31e-05 ***
## race_grother 0.30120 1.35148 0.42328 0.712 0.476729
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## momeducollege 0.06207 16.1103 0.03182 0.1211
## momedujcollege 0.41299 2.4214 0.29808 0.5722
## momeduprimarysch 1.91234 0.5229 1.24046 2.9481
## race_grhisp 1.90062 0.5261 1.31173 2.7539
## race_grnahn 2.76183 0.3621 1.50071 5.0827
## race_grnh asian 1.38842 0.7202 0.61817 3.1184
## race_grnh black 2.39079 0.4183 1.57450 3.6303
## race_grother 1.35148 0.7399 0.58953 3.0982
##
## Concordance= 0.79 (se = 0.022 )
## Rsquare= 0.044 (max possible= 0.479 )
## Likelihood ratio test= 214.4 on 8 df, p=0
## Wald test = 161.8 on 8 df, p=0
## Score (logrank) test = 255.5 on 8 df, p=0Cox model with additive frailty
Now we fit the Cox model and doing fraily by the school the child attends. I use the weights calculated above, standardized to the within cluster sample size.
# Fit the model
fitl2 <- coxme(Surv(time = time, event = povtran) ~ momedu +
race_gr + (1 | s2_id), e.long, weights = swts)
summary(fitl2)## Cox mixed-effects model fit by maximum likelihood
## Data: e.long
## events, n = 166, 4757 (1084 observations deleted due to missingness)
## Iterations= 21 154
## NULL Integrated Fitted
## Log-likelihood -1549.173 -1385.129 -1239.421
##
## Chisq df p AIC BIC
## Integrated loglik 328.09 9.0 0 310.09 282.08
## Penalized loglik 619.50 126.2 0 367.10 -25.64
##
## Model: Surv(time = time, event = povtran) ~ momedu + race_gr + (1 | s2_id)
## Fixed coefficients
## coef exp(coef) se(coef) z p
## momeducollege -2.7392476 0.06461895 0.3702210 -7.40 1.4e-13
## momedujcollege -0.8606695 0.42287885 0.1888406 -4.56 5.2e-06
## momeduprimarysch 1.0056573 2.73370341 0.2725945 3.69 2.2e-04
## race_grhisp 0.5413800 1.71837655 0.2370643 2.28 2.2e-02
## race_grnahn 0.8579181 2.35824592 0.4578572 1.87 6.1e-02
## race_grnh asian 0.7495192 2.11598242 0.4839161 1.55 1.2e-01
## race_grnh black 1.1923449 3.29479809 0.2895393 4.12 3.8e-05
## race_grother -0.1063573 0.89910338 0.5511786 -0.19 8.5e-01
##
## Random effects
## Group Variable Std Dev Variance
## s2_id Intercept 1.396481 1.950160anova(fitl1, fitl2)## Analysis of Deviance Table
## Cox model: response is Surv(time = time, event = povtran)
## Model 1: ~momedu + race_gr
## Model 2: ~momedu + race_gr + (1 | s2_id)
## loglik Chisq Df P(>|Chi|)
## 1 -1442.0
## 2 -1385.1 113.71 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(fitl1)## [1] 2899.966AIC(fitl2)## [1] 2731.246hist(fitl2$frail$s2_id)
So, we see that the frailty model has a large variance, and that the AIC is lower than the ordinary Cox model, suggesting better model fit. The model likelihood ratio test also confirms that the frailty model fits better.
Discrete time frailty model
Basic Discrete time model
First, we fit the basic discrete time model for our outcome. Here I’m going to use the new standardized weights I created above in a regular glm() model:
fit1 <- glm(povtran ~ as.factor(time) + momedu + race_gr - 1,
data = e.long, weights = swts, family = binomial(link = "cloglog"))## Warning in eval(family$initialize): non-integer #successes in a binomial
## glm!arm::display(fit1, detail = T)## glm(formula = povtran ~ as.factor(time) + momedu + race_gr -
## 1, family = binomial(link = "cloglog"), data = e.long, weights = swts)
## coef.est coef.se z value Pr(>|z|)
## as.factor(time)1 -2.47 0.15 -16.07 0.00
## as.factor(time)2 -2.81 0.18 -15.55 0.00
## as.factor(time)3 -3.00 0.19 -15.47 0.00
## momeducollege -2.64 0.34 -7.72 0.00
## momedujcollege -0.75 0.17 -4.55 0.00
## momeduprimarysch 0.61 0.22 2.78 0.01
## race_grhisp 0.65 0.19 3.41 0.00
## race_grnahn 0.99 0.31 3.15 0.00
## race_grnh asian 0.35 0.42 0.85 0.39
## race_grnh black 0.79 0.21 3.74 0.00
## race_grother 0.24 0.42 0.56 0.57
## ---
## n = 4757, k = 11
## residual deviance = 1394.0, null deviance = 9233.2 (difference = 7839.2)