I am fitting a continuous-time Markov model to a panel dataset using the R package MSM. Because I am interested in sex-differences in transition rates, I fit the model with covariate sex ("M" or "F") by running
model_object <- msm(
formula = state ~ nr_years,
subject = id_var,
qmatrix = M, # matrix encoding allowed transitions between states
data = panel_data,
covariates = ~ sex,
control = list(fnscale = 40000, maxit = 1e6) # got these from the help pages
)
After fitting the model I obtain the transition rate matrix using
qmatrix.msm(model_object, covariates = list(sex = "M"))
qmatrix.msm(model_object, covariates = list(sex = "F"))
These lines the exact same transition rate matrix. This is a bit unexpected to be, because when I use the hazard.msm function to extract hazard ratios, there are some differences between sexes. Some results are even statistically significant.
Am I using the function wrong?
This is probably because your sex variable is coded as a character instead of a factor.
Just converting sex into factor before to run the model should do the trick here:
panel_data$sex <- factor(panel_data$sex)
Here is a reproducible example to illustrate this tricky msm package behaviour:
Let's first create a sex variable on cav dataset as character (sex.chr) and as factor (sex.fct)
library(msm)
library(dplyr)
cav2 <-
cav %>%
mutate(
sex.chr = ifelse(sex == 1, 'M', 'F'),
sex.fct = factor(sex.chr)
)
we have to define the initial qmatrix too:
twoway4.q <- rbind(c(-0.5, 0.25, 0, 0.25), c(0.166, -0.498, 0.166, 0.166),
c(0, 0.25, -0.5, 0.25), c(0, 0, 0, 0))
If you use the character version of sex as covariate you will get tthe same results for both covariates (the covariates are basicaly ignored)
cav.msm.chr <- msm( state ~ years, subject = PTNUM, data = cav2,
qmatrix = twoway4.q, deathexact = 4, covariates = ~ sex.chr,
control = list ( trace = 2, REPORT = 1 ) )
qmatrix.msm(cav.msm.chr, covariates = list(sex.chr = "M"))
State 1 State 2 State 3
State 1 -0.17747 (-0.19940,-0.15795) 0.13612 ( 0.11789, 0.15718) 0
State 2 0.21992 ( 0.16100, 0.30040) -0.60494 (-0.70972,-0.51563) 0.33409 ( 0.26412, 0.42261)
State 3 0 0.12913 ( 0.07728, 0.21578) -0.41050 (-0.52552,-0.32065)
State 4 0 0 0
State 4
State 1 0.04135 ( 0.03245, 0.05268)
State 2 0.05092 ( 0.01808, 0.14343)
State 3 0.28137 ( 0.21509, 0.36808)
State 4 0
qmatrix.msm(cav.msm.chr, covariates = list(sex.chr = "F"))
State 1 State 2 State 3
State 1 -0.17747 (-0.19940,-0.15795) 0.13612 ( 0.11789, 0.15718) 0
State 2 0.21992 ( 0.16100, 0.30040) -0.60494 (-0.70972,-0.51563) 0.33409 ( 0.26412, 0.42261)
State 3 0 0.12913 ( 0.07728, 0.21578) -0.41050 (-0.52552,-0.32065)
State 4 0 0 0
State 4
State 1 0.04135 ( 0.03245, 0.05268)
State 2 0.05092 ( 0.01808, 0.14343)
State 3 0.28137 ( 0.21509, 0.36808)
State 4 0
But using the factorial version of sex will give you the expected results
cav.msm.fct <- msm( state ~ years, subject = PTNUM, data = cav2,
qmatrix = twoway4.q, deathexact = 4, covariates = ~ sex.fct,
control = list ( trace = 2, REPORT = 1 )
qmatrix.msm(cav.msm.fct, covariates = list(sex.fct = "M"))
State 1 State 2
State 1 -1.234e-01 (-1.750e-01,-8.693e-02) 7.667e-02 ( 4.630e-02, 1.270e-01)
State 2 2.838e-01 ( 1.139e-01, 7.075e-01) -6.435e-01 (-1.115e+00,-3.714e-01)
State 3 0 1.416e-01 ( 1.852e-02, 1.083e+00)
State 4 0 0
State 3 State 4
State 1 0 4.668e-02 ( 2.727e-02, 7.989e-02)
State 2 3.597e-01 ( 1.804e-01, 7.170e-01) 1.938e-05 ( 3.702e-66, 1.014e+56)
State 3 -8.207e-01 (-1.587e+00,-4.244e-01) 6.791e-01 ( 3.487e-01, 1.323e+00)
State 4 0
qmatrix.msm(cav.msm.fct, covariates = list(sex.fct = "F"))
State 1 State 2 State 3
State 1 -0.17747 (-0.19940,-0.15795) 0.13612 ( 0.11789, 0.15718) 0
State 2 0.21992 ( 0.16100, 0.30040) -0.60494 (-0.70972,-0.51563) 0.33409 ( 0.26412, 0.42261)
State 3 0 0.12913 ( 0.07728, 0.21578) -0.41050 #=(-0.52552,-0.32065)
State 4 0 0 0
State 4
State 1 0.04135 ( 0.03245, 0.05268)
State 2 0.05092 ( 0.01808, 0.14343)
State 3 0.28137 ( 0.21509, 0.36808)
State 4 0
```