My dataset is very similar to the dataset 'Melanoma' included in the RiskRegression package : 3307 patients, 502 events of interest (fracture), 264 deaths (competing risk). The time is the years after bone examination (DXA) and status is coded in this way O=censored,1=fracture,2=death). I am trying to fit a Fine-Gray model with interaction, but when I introduce an interaction term under the form of var1 * var2) I receive an error message :
« Error in design[pos, , drop = FALSE] : subscript out of bounds » .
Here is my code :
fgr<-FGR(Hist(time,status)~age+htot_bmd+tot_bmd+amof+PR+atcdtfam+AlcFR+PR+BMI3C+malchronFR+malchronFR*BMI3C+atcdtfam*PR,data=df2,cause=1)
I tried the code provided in the paper of Zhongheng et al. "Model validation for competing risks data" with the data set 'Melanoma' introducing an interaction but the same error message appears.
Is it possible to introduce an interaction with FGR and how to do it ?
Thanks
You can do with your data with the following code:
> library(riskRegression)
> library(survival)
> library(prodlim)
> library(cmprsk)
> library(readxl)
> df2 <- read_xlsx("/Users/zhang/Downloads/df2.xlsx")
New names:
* `` -> ...1
> df2
# A tibble: 300 x 14
...1 neck_bmd htot_bmd tot_bmd age AlcFR PR atcdtfam malchronFR amof BMI3C time event
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 1 0.960 0.953 1.04 79.1 0 0 0 0 2 3 9.00 Cen
2 2 0.612 0.620 0.988 79.2 0 0 0 0 0 3 4.76 MOF
3 3 0.880 0.990 0.827 78.6 0 0 0 1 1 2 9.14 Cen
4 4 0.869 0.905 0.866 79.0 0 0 0 0 0 2 9.11 Cen
5 5 0.863 0.991 1.17 79.0 1 0 1 0 0 2 10.2 Cen
6 6 0.722 0.902 0.842 78.8 0 0 0 0 0 2 9.09 Cen
7 7 0.853 0.929 1.33 76.9 0 0 0 0 0 3 10.1 Cen
8 8 0.830 0.912 0.947 77.0 0 0 0 1 0 2 8.13 Cen
9 9 0.872 0.968 1.22 77.2 1 0 0 0 0 2 8.12 Cen
10 10 0.639 0.776 0.822 76.7 0 0 0 1 0 2 8.12 Cen
# … with 290 more rows, and 1 more variable: status <dbl>
> modMatrix <- model.matrix(~age+htot_bmd+tot_bmd+amof+atcdtfam+
+ AlcFR+PR+BMI3C*malchronFR+neck_bmd,df2)[,-1]
> dtInteraction <- cbind(data.frame( modMatrix),
+ status=df2$status, time=df2$time)
> fgr.Interaction<- FGR(as.formula(paste("Hist(time,status)~",
+ paste(names(dtInteraction[1:11]),collapse = "+"))),
+ data = dtInteraction,cause = 1)
> score.cv<-riskRegression::Score(list("Fine-Gray"= fgr.Interaction),
+ formula = Hist(time,status)~1,
+ data=dtInteraction,times = sort(unique(dtInteraction$time))[25:200],
+ cens.method="jackknife",
+ se.fit=1L,plots="calibration")
> plotCalibration(score.cv,times = df2$time[11],
+ cens.method="local")