Plotting an exposure response curve from qgcomp.cox.boot

joint_exposure_response_qgcomp_cox_boot.R

library(qgcomp)
library(survival)
library(ggplot2)
set.seed(50)
N=200
dat <- data.frame(time=(tmg <- pmin(.1,rweibull(N, 10, 0.1))), 
                  d=1.0*(tmg<0.1), x1=runif(N), x2=runif(N), z=runif(N))
expnms=paste0("x", 1:2)
f = survival::Surv(time, d)~x1 + x2 + I(x2*x1)
## Not run: 
# Note that MCsize should be set much higher generally - it is set to 2000 so that this can run in a short period of time
(obj <- qgcomp.cox.boot(f, expnms = expnms, data = dat, B=10, MCsize=2000, degree=2, q=4))

obj$psi
obj$covmat.psi


q = 4
qoff = 0.5/q
testX = seq(0, q-1, 0.1)
plotpercentile = seq(qoff, 1-qoff, length.out=length(testX))

testDesign = cbind(psi1=testX, psi2=testX^2)

lnhr = testDesign %*% obj$psi
lnhr_stderr = rep(0, length(lnhr))
for(j in 1:length(lnhr_stderr)){
  #Var(aX+bY)=a2Var(X)+b2Var(Y)+2abCov(X,Y)
  lnhr_stderr[j] = sqrt(sum((testDesign[j,] %*% t(testDesign[j,])) * obj$covmat.psi))
  # qgcomp internal function for doing the same thing
  #lnhr_stderr[j] = qgcomp::se_comb(c("psi1", "psi2"), obj$covmat.psi, testDesign[j,])
}
  
ggplot()+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  geom_ribbon(aes(x=plotpercentile, ymin=lnhr-1.96*lnhr_stderr, ymax=lnhr+1.96*lnhr_stderr), alpha=0.4) +
  geom_line(aes(x=plotpercentile, y=lnhr)) + 
  labs(x="Percentile of all exposures", y="ln(HR)") + 
  lims(x=c(0,1)) + 
  theme_classic()