fauvernierma / survPen

hazard and excess hazard modelling with multidimensional penalized splines
GNU General Public License v3.0
8 stars 2 forks source link

Confidence Intervals for Hazard Ratios #21

Closed modedec closed 5 years ago

modedec commented 5 years ago

How can the Hazard Ratio confidence intervals be calculated

1)survPen(...) -> model 2)predict(model,...) -> pred_1 and pred_2 3)HR=pred_1$haz/pred_2$haz

Could one just do error propagation? f(x,y)=x/y var(HR)=(1/y*dy)+(x/y^2 dx)^2+(1/y^2 dy dx)^2

fauvernierma commented 5 years ago

Hello @modedec and thank you for your question

The simplest way to do it is with the log hazard ratio :

log(HR) = log(haz1) - log(haz2) = (X1 - X2)%*%model$coef

where X1 and X2 are the model matrix associated with each log-hazard

Thus, the variance for the log hazard ratio is: Var(log(HR)) = (X1 - X2) %% model$Vp %% t(X1 - X2)

With survPen, you can use the option: type="lpmatrix" for the predict function to get X1 and X2, this would give:

conf.int=0.95 qt.norm<-qnorm(1-(1-conf.int)/2)

newdata1 <- .... newdata2 <- ...

myMat1<- predict(model,newdata1,type="lpmatrix") myMat2<- predict(model,newdata2,type="lpmatrix")

myMat<-myMat1-myMat2

log.haz.ratio<-as.vector(myMat%*%model$coef)

haz.ratio<-as.vector(exp(pred.haz.ratio))

confidence intervals for hazard

std<-sqrt(rowSums((myMat%%model$Vp)myMat)) haz.ratio.inf<-as.vector(exp(log.haz.ratio-qt.normstd)) haz.ratio.sup<-as.vector(exp(log.haz.ratio+qt.normstd))

fauvernierma commented 5 years ago

Don't hesitate if something is not clear and close the issue if this is ok for you

Thanks and have a great day

fauvernierma commented 5 years ago

Erratum, in my answer, one must read :

haz.ratio<-as.vector(exp(log.haz.ratio)) instead of haz.ratio<-as.vector(exp(pred.haz.ratio))

Also, there are * missing at the very end:

haz.ratio.inf<-as.vector(exp(log.haz.ratio-qt.normstd)) haz.ratio.sup<-as.vector(exp(log.haz.ratio+qt.normstd))

modedec commented 5 years ago

Thank you for taking the time to answer our questions. We gonna try to implement your answer in the next couple of days and will come back to you.

Thx for the great package, have a great day

modedec commented 5 years ago

We checked the result with the cox timesplitting approach. Your solution works fine:

  lpmatrix.risk <- predict(model,df.risk,type="lpmatrix")
  lpmatrix.ref <- predict(model,df.ref,type="lpmatrix")
  myMat<-lpmatrix.risk-lpmatrix.ref
  std<-sqrt(rowSums((myMat%*%model$Vp)*myMat))
  log.haz.ratio<-as.vector((myMat)%*%model$coef)
  haz.ratio.upper<-as.vector(exp(log.haz.ratio+1.96*std))
  haz.ratio.lower<-as.vector(exp(log.haz.ratio-1.96*std))
  haz.ratio<-as.vector(exp(log.haz.ratio))

Thank you for your support, we highly appreciate it.