Provide support for handling categorical variables in both sub-models.
Provide the anova() function to compare two fitted joint models.
Add the simulate argument in the survfitjmcs() function to obtain the conditional probabilities using the Gauss-Hermite quadrature rule for numerical integration.
Adjust the label position of y axis for clarity purposes when include.y = TRUE in the plot.survfitjmcs() function.