search

chjackson / flexsurv

The flexsurv R package for flexible parametric survival and multi-state modelling
http://chjackson.github.io/flexsurv/
53 stars 28 forks source link

path detection in fmixmsm() #137

Closed kkmann closed 2 years ago

kkmann commented 2 years ago

Hi,

I get the follwoing unexpected behaviour with fmixmxm().

Here is my code

tbl_data <- data.frame(
  to = factor(1:3, levels = 1:3, labels = c("response", "progression", "death")),
  dt = rep(1, 3),
  status = rep(1, 3)
)
fit1 <- flexsurvmix(
  Surv(dt, status) ~ 1,
  event = tbl_data$to,
  data = tbl_data,
  dists = c(
    "response" = "exponential",
    "progression" = "exponential",
    "death" = "exponential"
  )
)

tbl_data <- data.frame(
  to = factor(1:2, levels = 1:2, labels = c("progression", "death")),
  dt = rep(1, 2),
  status = rep(1, 2)
)
fit2 <- flexsurvmix(
  Surv(dt, status) ~ 1,
  event = tbl_data$to,
  data = tbl_data,
  dists = c(
    "progression" = "exponential",
    "death" = "exponential"
  )
)

tbl_data <- data.frame(
  to = factor(1:1, levels = 1:1, labels = c("death")),
  dt = rep(1, 1),
  status = rep(1, 1)
)
fit3 <- flexsurvmix(
  Surv(dt, status) ~ 1,
  event = tbl_data$to,
  data = tbl_data,
  dists = c(
    "death" = "exponential"
  )
)

fit <- fmixmsm(
  "stable" = fit1,
  "response" = fit2,
  "progression" = fit3
)

Here is the output from inspecting fit:

> fit
$stable
Call:
flexsurvmix(formula = Surv(dt, status) ~ 1, data = tbl_data, 
    event = tbl_data$to, dists = c(response = "exponential", 
        progression = "exponential", death = "exponential"))

Estimates: 
     component         dist  terms    est      est.t    se
1     response               prob1  0.333         NA  1.73
2  progression               prob2  0.333  -7.41e-09  1.41
3        death               prob3  0.333  -7.41e-09  1.41
4     response  exponential   rate  1.000   0.00e+00  1.00
5  progression  exponential   rate  1.000   0.00e+00  1.00
6        death  exponential   rate  1.000   0.00e+00  1.00

Log-likelihood = -6.295837, df = 5
AIC = 22.59167

$response
Call:
flexsurvmix(formula = Surv(dt, status) ~ 1, data = tbl_data, 
    event = tbl_data$to, dists = c(progression = "exponential", 
        death = "exponential"))

Estimates: 
     component         dist  terms  est      est.t    se
1  progression               prob1  0.5         NA  1.41
2        death               prob2  0.5  -1.11e-13  1.41
3  progression  exponential   rate  1.0   0.00e+00  1.00
4        death  exponential   rate  1.0   0.00e+00  1.00

Log-likelihood = -3.386294, df = 3
AIC = 12.77259

$progression
Call:
flexsurvmix(formula = Surv(dt, status) ~ 1, data = tbl_data, 
    event = tbl_data$to, dists = c(death = "exponential"))

Estimates: 
   component         dist  terms  est  est.t  se
1      death               prob1    1     NA  NA
2      death  exponential   rate    1      0   1

Log-likelihood = -1, df = 1
AIC = 4

attr(,"pathway_str")
list()

I am confused that the pathways attribute is empty since there are many clear pathways here (s->d; s->r->p->d; s->p->d). If I look at the reduced model with the first two components only, all seems to be fine:

fmixmsm(
+     "stable" = fit1,
+     "response" = fit2
+ )
$stable
Call:
flexsurvmix(formula = Surv(dt, status) ~ 1, data = tbl_data, 
    event = tbl_data$to, dists = c(response = "exponential", 
        progression = "exponential", death = "exponential"))

Estimates: 
     component         dist  terms    est      est.t    se
1     response               prob1  0.333         NA  1.73
2  progression               prob2  0.333  -7.41e-09  1.41
3        death               prob3  0.333  -7.41e-09  1.41
4     response  exponential   rate  1.000   0.00e+00  1.00
5  progression  exponential   rate  1.000   0.00e+00  1.00
6        death  exponential   rate  1.000   0.00e+00  1.00

Log-likelihood = -6.295837, df = 5
AIC = 22.59167

$response
Call:
flexsurvmix(formula = Surv(dt, status) ~ 1, data = tbl_data, 
    event = tbl_data$to, dists = c(progression = "exponential", 
        death = "exponential"))

Estimates: 
     component         dist  terms  est      est.t    se
1  progression               prob1  0.5         NA  1.41
2        death               prob2  0.5  -1.11e-13  1.41
3  progression  exponential   rate  1.0   0.00e+00  1.00
4        death  exponential   rate  1.0   0.00e+00  1.00

Log-likelihood = -3.386294, df = 3
AIC = 12.77259

attr(,"pathways")
attr(,"pathways")[[1]]
[1] "stable"      "response"    "progression"

attr(,"pathways")[[2]]
[1] "stable"   "response" "death"   

attr(,"pathways")[[3]]
[1] "stable"      "progression"

attr(,"pathways")[[4]]
[1] "stable" "death" 

attr(,"pathway_str")
[1] "stable-response-progression" "stable-response-death"      
[3] "stable-progression"          "stable-death"

Is this expected behaviour (if so why?)?

Thanks,

Kevin

kkmann commented 2 years ago

commenting out

https://github.com/chjackson/flexsurv-dev/blob/9dbfef739e37b4b8a04ff36000843a727f173d2d/R/fmixmsm.R#L45

resolves the problem but some form of protection against infinite recursion is probably required

chjackson commented 2 years ago

Thanks for the report. I've reworked this pathway detection function now, and it should be working in https://github.com/chjackson/flexsurv-dev/commit/39465114cbf1b9f68b9548e3c10bf4559250462a

kkmann commented 2 years ago

Hi Chris,

thanks for the quick reply - I checked it and the pathway detection works for me too now. I then ran 

> ppath_fmixmsm(fit, B = 20)
  final                           pathway       val      lower     upper
1 death stable-response-progression-death 0.1666667 0.02346216 0.4203543
2 death             stable-response-death 0.1666667 0.02854172 0.4040497
3 death                      stable-death 0.3333333 0.04255448 0.719964

should the path probabilities not add to one? Maybe my data example is too simple? I thought in the mixture model, artificial censoring events are not required for non-observed transitions, right?

chjackson commented 2 years ago

With your example above, I get another row with stable-progression-death

  final                           pathway       val       lower     upper
1 death stable-response-progression-death 0.1666667 0.009222188 0.5191553
2 death             stable-response-death 0.1666667 0.002622412 0.4903672
3 death          stable-progression-death 0.3333333 0.092747042 0.7017576
4 death                      stable-death 0.3333333 0.067210043 0.8036272
kkmann commented 2 years ago

Sorry, only now got back to checking. I can confirm it's working. Thanks again for the quick fix. Really impressive to see what flexsurv can do now!