EM algorithm R code on Cox PH model with frailty

Question

Let say I have a 'kidney catheter' data set. Data are about the recurrence times to infection, at the point of insertion of the catheter, for kidney patients using portable dialysis equipment. Catheters may be removed for reasons other than infection, in which case the observation is censored. Each patient has exactly 2 observations.

If I want to fit a Cox PH model with random effect (here Gamma frailty) using EM algorithm. By using built in survival package R code coxph(), I can easily do it like this

library(survival)
data(kidney)
fit<-coxph(Surv(time, status) ~ age + sex + 
           frailty(id, dist='gamma', method='em'), kidney)

But if I want to write a step by step function for EM algorithm, how can I proceed?

Answer 1

coxph() actually implements a penalised log-likelihood approach which turns out to return the same estimates as the EM algorithm in the case of gamma frailties when method="em"; see Therneau and Grambsch (2000, Section 9.6). (method actually refers to the method used to select a solution for theta, the heterogeneity parameter, not to the estimation procedure). Both algorithms are clearly detailed in Duchateau and Janssen (2008, Chapter 5). Implementing the EM algorithm might require quite a lot of work, but might be possible following these lines. By the way, there is a SAS macro called gamfrail written by Klein that already does the job, even though it is not very user-friendly. It can be downloaded here together with a guide.