8.4 Parametric frailty models

It is possible utilize the connection between Poisson regression and the piecewise constant proportional hazards model discussed in Chapter 6 to fit parametric frailty models. We look at the fertility data again. The standard analysis without frailty effects is

fit0 <- phreg(Surv(next.ivl, event) ~ parity + ses,
              dist = "pch", cuts = 1:13, data = fe)
summary(fit0)

## Covariate           Mean       Coef     Rel.Risk   S.E.    LR p
## parity              4.242    -0.131     0.878     0.005   0.0000 
## ses                                                       0.0000 
##           farmer    0.491     0         1 (reference)
##          unknown    0.186     0.094     1.099     0.029
##            upper    0.018     0.086     1.089     0.083
##            lower    0.305    -0.150     0.860     0.026
## 
## Events                    8458 
## Total time at risk         29806 
## Max. log. likelihood      -14087 
## LR test statistic         838.09 
## Degrees of freedom        4 
## Overall p-value           0

For testing the presence of a random effect over women, an alternative is to use glmmML in eha.

library(glmmML)
fe13 <- survSplit(fe, end = "next.ivl", event = "event", 
                  cut = 1:13, episode = "years", 
                  start = "start") #1
fe13$years <- as.factor(fe13$years) # 2
fe13$offs <- log(fe13$next.ivl - fe13$start) # 3
fit1 <- glmmML(event ~ parity + ses + years + offset(offs), 
               cluster = id, family = poisson, method = "ghq", 
               data = fe13, n.points = 9)
summary(fit)

## Cox mixed-effects model fit by maximum likelihood
##   Data: oldmort
##   events, n = 888, 3340 (3155 observations deleted due to missingness)
##   Iterations= 16 84 
##                     NULL Integrated    Fitted
## Log-likelihood -5702.662   -5693.43 -5646.252
## 
##                    Chisq    df          p   AIC     BIC
## Integrated loglik  18.46  4.00 1.0012e-03 10.46   -8.69
##  Penalized loglik 112.82 49.23 6.7155e-07 14.36 -221.39
## 
## Model:  Surv(enter, exit, event) ~ sex + civ + (1 | m.id) 
## Fixed coefficients
##                  coef exp(coef)   se(coef)     z       p
## sexfemale  -0.2026325 0.8165783 0.07190559 -2.82 0.00480
## civmarried -0.4653648 0.6279060 0.12144548 -3.83 0.00013
## civwidow   -0.3305040 0.7185615 0.11964284 -2.76 0.00570
## 
## Random effects
##  Group Variable  Std Dev    Variance  
##  m.id  Intercept 0.23719205 0.05626007

Note the use of the function survSplit (from the survival package) (# 1), the transformation to factor for the slices of time (years) created by survSplit, and the creation of an offset (# 3). This is described in detail in Chapter 6.

The conclusion is very much the same as with the nonparametric (coxme) approach: The clustering effect of mother is very strong and must be taken into account in the analysis of birth intervals. The nonparametric approach is easier to use and recommended, but see the next session.