Frailty survival — R cross-check (parfm, coxme)¶

Companion to frailty_python.ipynb / frailty_pymc.ipynb. survreg has no frailty term, so we use:

  • parfm — parametric Weibull with log-normal frailty: the exact match to our model.
  • coxme — Cox model with a Gaussian random intercept (1|id): the semiparametric frailty.
In [1]:
ok <- suppressWarnings(suppressMessages(require(parfm)))
if (!ok) { install.packages('parfm', repos='https://cloud.r-project.org'); library(parfm) }
d <- read.csv('kidney.csv'); d$disease <- relevel(factor(d$disease), ref='Other')
d$female <- as.integer(d$sex==2); d$age_c <- d$age - mean(d$age)
fit <- parfm(Surv(time, status) ~ age_c + female + disease, cluster='id', data=d,
             dist='weibull', frailty='lognormal')
print(fit)
cat(sprintf('\nWeibull shape (rho) = %.3f   frailty sd = sqrt(sigma2) = %.3f\n',
            fit['rho','ESTIMATE'], sqrt(fit['sigma2','ESTIMATE'])))
cat(sprintf('female HR = %.2f   PKD HR = %.2f\n', exp(fit['female','ESTIMATE']), exp(fit['diseasePKD','ESTIMATE'])))

Frailty distribution: lognormal 
Baseline hazard distribution: Weibull 
Loglikelihood: -329.885 

           ESTIMATE SE    p-val    
sigma2      0.315   0.325          
rho         1.162   0.153          
lambda      0.015   0.010          
age_c       0.003   0.014 0.825    
female     -1.877   0.473 <.001 ***
diseaseAN   0.606   0.489 0.216    
diseaseGN   0.121   0.492 0.806    
diseasePKD -1.142   0.763 0.134    
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kendall's Tau: 0.13 

Weibull shape (rho) = 1.162   frailty sd = sqrt(sigma2) = 0.562
female HR = 0.15   PKD HR = 0.32
In [2]:
ok2 <- suppressWarnings(suppressMessages(require(coxme)))
if (!ok2) { install.packages('coxme', repos='https://cloud.r-project.org'); library(coxme) }
cm <- coxme(Surv(time, status) ~ age_c + female + disease + (1|id), data=d)
print(cm)
cat(sprintf('\ncoxme: female HR = %.2f   frailty sd = %.3f\n', exp(fixef(cm)['female']), sqrt(as.numeric(VarCorr(cm)$id))))

Cox mixed-effects model fit by maximum likelihood
  Data: d
  events, n = 58, 76
  Iterations= 5 28 
                    NULL Integrated    Fitted
Log-likelihood -187.9028  -179.0217 -173.7424

                  Chisq   df          p   AIC   BIC
Integrated loglik 17.76 6.00 0.00685470  5.76 -6.60
 Penalized loglik 28.32 9.11 0.00090634 10.10 -8.66

Model:  Surv(time, status) ~ age_c + female + disease + (1 | id) 
Fixed coefficients
                   coef exp(coef)   se(coef)     z       p
age_c       0.003755551 1.0037626 0.01206409  0.31 0.76000
female     -1.552913501 0.2116305 0.39038849 -3.98 0.00007
diseaseAN   0.358074826 1.4305727 0.43752739  0.82 0.41000
diseaseGN   0.119156978 1.1265467 0.44126751  0.27 0.79000
diseasePKD -1.329715978 0.2645524 0.69630860 -1.91 0.05600

Random effects
 Group Variable  Std Dev   Variance 
 id    Intercept 0.3428171 0.1175235

coxme: female HR = 0.21   frailty sd = 0.343

Results¶

  • parfm (parametric Weibull, log-normal frailty) is the direct match and confirms the Bayesian fits: shape ≈ 1.16, female HR ≈ 0.15, PKD HR ≈ 0.32, frailty sd ≈ 0.56 — right alongside from-scratch (0.65) and PyMC (0.76).
  • coxme (semiparametric Cox frailty) gives the same covariate story (female HR ≈ 0.21) with a smaller frailty sd ≈ 0.34 — the Cox baseline absorbs some structure the Weibull attributes to the frailty.
  • The frailty is robustly non-zero across parfm, coxme, and both Bayesian engines (sd ≈ 0.3–0.8). Only a Cox-frailty boundary ML fit (coxph(... + frailty(id))) collapses it to ~0 — the usual variance-component-at-the-boundary artefact with few small clusters.
  • Engines agree: from-scratch ≈ PyMC ≈ parfm on the Weibull frailty model; coxme the semiparametric cross-check.