Competing risk analysis with interaction term I am trying to do a Fine & Gray competing risk analysis for breast cancer death in patients with breast cancer. I want to assess the interaction betweeen race and chemotherapy on breast cancer death.
My database is df2 and has all categorical variables. survival is survival time, and cod_breast is the variable containing reason for death (failcode = 1, censoring = 2)
library(cmprsk)
library(tidyverse)

My overal model is:
cov1 <- model.matrix(~  chemotherapy + race +  age_di + stage + grade+
                                        radiation + surgery, data = df2)[, -1]

crr.cov1 <- crr(ftime=df2$survival, fstatus=df2$cod_breast, failcode=1, cencode=2, cov1=cov1)

My model with interaction is:
cov1int <- model.matrix(~  chemotherapy + race + chemotherapy:race + age_di + stage + grade+
                                        radiation + surgery, data = df2)[, -1]

crr.cov1int <- crr(ftime=df2$survival, fstatus=df2$cod_breast, failcode=1, cencode=2, cov1=cov1int)

My question is how do I assess the significance of the interaction in this sceneario? For Cox models, I usually do anova(crr.cov1,crr.cov1int,test="Chisq"). However, this doesnt work with Fine & Gray.

If you accept the pseudolikelihood ratio test (the default method of inference in the Fine Gray model see https://www.rdocumentation.org/packages/cmprsk/versions/2.2-11/topics/summary.crr), then testing the interaction effect is no different than testing the marginal effect. Your specific problem here is that the test of interest is larger than 1 degree of freedom, because chemo and race together (presumably) have more than 1 effect. OR at least I have to assume so, because you don't provide any data tabulations (please provide data tabulations when asking on CV). Luckily, it's not too hard to build this out from default. As an example using simulation:
set.seed(10)
ftime <- rexp(200)
fstatus <- sample(0:2,200,replace=TRUE)
X <- matrix(sample(letters[1:3], 400, replace=T), 200)
colnames(X) <- c('x', 'w')
Xf <- model.matrix(~ x*w, data=as.data.frame(X))[, -1]
z <- crr(ftime,fstatus,Xf)

gives
> z
convergence:  TRUE 
coefficients:
      xb       xc       wb       wc    xb:wb    xc:wb    xb:wc    xc:wc 
-0.14510 -0.83210  0.30300  0.17500  0.35890  0.45140 -0.05315  0.67440 
standard errors:
[1] 0.5927 0.6785 0.6103 0.5814 0.7481 0.8635 0.7623 0.8167
two-sided p-values:
   xb    xc    wb    wc xb:wb xc:wb xb:wc xc:wc 
 0.81  0.22  0.62  0.76  0.63  0.60  0.94  0.41 

to test the 'x' and 'w' interaction terms (4 of), fit the lesser model:
Xr <- model.matrix(~ x + w, data=as.data.frame(X))[, -1]
y <- crr(ftime,fstatus,Xr)
p <- pchisq(-2 * (y$loglik - z$loglik), df=4, lower.tail = F)

Each object of class crr has a loglik value, so the by-hand (p)LRT test statistic (and p-value) is:
Ts <- -2 * (y$loglik - z$loglik)
p <- pchisq(Ts, df=4, lower.tail = F)*2

which gives:
> p
[1] 0.2029832

