CecileProust-Lima / lcmm

R package lcmm
https://CecileProust-Lima.github.io/lcmm/
48 stars 13 forks source link

Estimating hazard ratios from JLCMM output #238

Open law852 opened 4 months ago

law852 commented 4 months ago

I'm looking for some guidance on estimating hazard ratios from a JLCMM object with class-specific Weibull hazard curves (e.g. jlcmm(fixed = y ~ x, mixture = ~time, random = ~1, subject = "pid", ng = 4, survival = Surv(time, event2) ~ b1 + b2, hazard = "Weibull", hazardtype = "Specific", link = "3-equi-splines", data = data).

Thank you!

VivianePhilipps commented 4 months ago

Hi,

With Weibull specific hazards, the ratio between the hazard of two classes will depend on time. You get the hazard of each latent class in the $predSurv output, for a grid of time. You can then compute the ratio based of these values. For example for the ratio between the (baseline) hazard of class 1 over class 2:

mj2 <- Jointlcmm(normMMSE ~ age65 + I(age65^2) + CEP, mixture =~ age65 + I(age65^2), 
              random =~ age65 + I(age65^2), survival = Surv(age_init, agedem, dem) ~ CEP + male,
              hazard = "Weibull", subject = 'ID', data = paquidS, ng=2,
              B=c(-0.615, 0.11, 4.695, 0.102, 5.324, -0.81, 0.173, 59.65, 63.876, 14.762, 5.538, 
                            -14.01, -4.234, 14.38, 227.631, -228.341, 435.744, 58.487, -129.354, 41.967, 10.019))

hr_mj2 <- data.frame(time = mj2$predSurv[,1], HR = mj2$predSurv[,2] / mj2$predSurv[,3])

plot(hr_mj2[,1], hr_mj2[,2], type="l")

If you want it for a specific time, you'll need the following risk function:

baselineriskWeibull <- function(t, w1, w2, logscale = FALSE)
{
  if(!logscale)
   res <- as.numeric(w1^2 * w2^2 * (w1^2 * t)^(w2^2 - 1))
  else
   res <-  as.numeric(w1^2 * w2^2 * t^(w2^2 - 1))
  return(res)
}

and then compute the HR with at time 70 for example with: baselineriskWeibull(70, mj2$best[2], mj2$best[3]) / baselineriskWeibull(70, mj2$best[4], mj2$best[5])

Best,

Viviane