-1

Below is some simulated data. I would like some help in a) Plotting observed vs. predicted risk by 3 quantile groups Q1-Q3 and then run H-L test to report Chi squared statistic and p-value for goodness of fit for the subdistribution model. b) Aligning the x-axis labels to display horizontally in one line. (X-axis labels are hardcoded in the package function.) c) Getting the baseline subhazard Cumulative Incidence Function CIF at time=30.

library(simstudy)
library(data.table)
library(survival)
library(riskRegression)
library(cmprsk)
library(prodlim)

#Generate time to event data
d1 <- defData(varname = "x1", formula = .5, dist = "binary")
d1 <- defData(d1, "x2", .5, dist = "binary")
dS <- defSurv(varname = "event_1", formula = "-12 - 0.1*x1 - 0.2*x2", shape = 0.3)
dS <- defSurv(dS, "event_2", "-12 - 0.3*x1 - 0.2*x2", shape = 0.3)
dS <- defSurv(dS, "event_3", "-12 - 0.4*x1 - 0.3*x2", shape = 0.3)
dS <- defSurv(dS, "censor", "-13", shape = 0.3)
set.seed(2140)
dtCov <- genData(3001, d1)
dtSurv <- genSurv(dtCov, dS)
head(dtSurv)
f <- "(time==censor)*0 + (time==event_1)*1 + (time==event_2)*2 + (time==event_3)*3"
cdef <- defDataAdd(varname = "time", 
                   formula = "pmin(censor, event_1, event_2, event_3)", dist = "nonrandom")
cdef <- defDataAdd(cdef, varname = "event", 
                   formula = f, 
                   dist = "nonrandom")
dtSurv_final <- addCompRisk(dtSurv, 
                            events = c("event_1", "event_2", "event_3", "censor"), 
                            timeName = "time", censorName = "censor")
head(dtSurv_final)

#Fit subdistribution model
crr<-riskRegression::FGR(Hist(time, event)~ x1+x2 ,data=dtSurv_final, cause=1)
crr.object<-riskRegression::Score(list(model1=crr),
          Hist(time,event)~1,data=dtSurv_final, cause=1, times=30, plots="cal")
#Store model output 
cal<-crr.object[["Calibration"]]$plotframe %>% as.data.frame()

#Need Hosmer lemeshow Chi square statistic and p-value to determine goodness of fit by comparing the risk groups in the plot below:
riskRegression::plotCalibration(crr.object,
                                bars=TRUE,
                                q=3,
                                cens.method="local",
                                ylim=c(0,.50),
                                names=paste0("Q",c(1:3)),
                                names.cex=1.1)

I tried getting the observed and predicted risks from cal but have been unable to do so.

Mahim
  • 11
  • 3
  • `Error in Hist(time, event) : could not find function "Hist"`. (Need better list of needed packages to be loaded.) Also, what was the reason you did not use one of the existing H-L test functions? (I found two such functions in the packages that I have already installed, but asking for package recommendations is off-topic here.) – IRTFM Aug 12 '23 at 22:09
  • @IRTFM , thanks for catching that! The package for Hist() is "prodlim". Could you please guide which packages provide the H-L test functions that can be used for this model output? Thanks! – Mahim Aug 14 '23 at 18:04
  • Learn to search with Google using the stem "CRAN ...".. I had no trouble with the glaringly obvious search terms using the `??` operator, but since you probably have fewer installed packages than I do, Google will probably present you with even more options.. – IRTFM Aug 14 '23 at 21:30

0 Answers0