Because this is such a long question I've broken it down into 2 parts; the first being just the basic question and the second providing details of what I've attempted so far.

Question - Short

How do you fit an individual frailty survival model in R? In particular I am trying to re-create the coefficient estimates and SE's in the table below that were found from fitting the a semi-parametric frailty model to this dataset link. The model takes the form:

h_i(t) = z_i h_0(t) exp(\beta'X_i)

where z_i is the unknown frailty parameter per each patient, X_i is a vector of explanatory variables, \beta is the corresponding vector of coefficients and h_0(t) is the baseline hazard function using the explanatory variables disease, gender, bmi & age ( I have included code below to clean up the factor reference levels).

Question - Long

I am attempting to follow and re-create the Modelling Survival Data in Medical Research text book example for fitting frailty mdoels. In particular I am focusing on the semi parametric model for which the textbook provides parameter and variance estimates for the normal cox model, lognormal frailty and Gamma frailty which are shown in the above table

I am able to recreate the no frailty model estimates using

library(dplyr)
library(survival)
dat <- read.table(
    "./Survival of patients registered for a lung transplant.dat",
    header = T
) %>% 
    as_data_frame %>% 
    mutate( disease = factor(disease, levels = c(3,1,2,4))) %>% 
    mutate( gender = factor(gender, levels = c(2,1)))

mod_cox <- coxph( Surv(time, status) ~ age + gender + bmi + disease   ,data = dat)
mod_cox

however I am really struggling to find a package that can reliably re-create the results of the second 2 columns. Searching online I found this table which attempts to summarise the available packages:

source

Below I have posted my current findings as well as the code I've used encase it helps someone identify if I have simply specified the functions incorrectly:

frailtyEM - Seems to work the best for gamma however doesn't offer log-normal models

frailtyEM::emfrail(
    Surv(time, status) ~ age + gender + bmi + disease + cluster(patient), 
    data = dat ,
    distribution = frailtyEM::emfrail_dist(dist = "gamma")
)

survival - Gives warnings on the gamma and from everything I've read it seems that its frailty functionality is classed as depreciated with the recommendation to use coxme instead.

coxph( 
    Surv(time, status) ~ age + gender + bmi + disease + frailty.gamma(patient), 
    data = dat 
)
coxph( 
    Surv(time, status) ~ age + gender + bmi + disease + frailty.gaussian(patient), 
    data = dat 
)

coxme - Seems to work but provides different estimates to those in the table and doesn't support gamma distribution

coxme::coxme(
    Surv(time, status) ~ age + gender + bmi + disease + (1|patient),
    data = dat
)

frailtySurv - I couldn't get to work properly and seemed to always fit the variance parameter with a flat value of 1 and provide coefficient estimates as if a no frailty model had been fitted. Additionally the documentation doesn't state what strings are support for the frailty argument so I couldn't work out how to get it to fit a log-normal

frailtySurv::fitfrail(
    Surv(time, status) ~ age + gender + bmi + disease + cluster(patient),
    dat = dat,
    frailty = "gamma"
)

frailtyHL - Produce warning messages saying "did not converge" however it still produced coeficiant estimates however they were different to that of the text books

mod_n <- frailtyHL::frailtyHL(
    Surv(time, status) ~ age + gender + bmi + disease + (1|patient),
    data = dat,
    RandDist = "Normal"
)
mod_g <- frailtyHL::frailtyHL(
    Surv(time, status) ~ age + gender + bmi + disease + (1|patient),
    data = dat,
    RandDist = "Gamma"
)

frailtypack - I simply don't understand the implementation (or at least its very different from what is taught in the text book). The function requires the specification of knots and a smoother which seem to greatly impact the resulting estimates.

parfm - Only fits parametric models; having said that everytime I tried to use it to fit a weibull proportional hazards model it just errored.

phmm - Have not yet tried

I fully appreciate given the large number of packages that I've gotten through unsuccessfully that it is highly likely that the problem is myself not properly understanding the implementation and miss using the packages. Any help or examples on how to successfully re-create the above estimates though would be greatly appreciated.

Regarding