Search code examples
rregressionsurvival-analysis

Fitting a Fine-Gray model with interaction in R (package RiskRegression - function FGR)


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


Solution

  • 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")
    

    enter image description here