I am trying to fit a generalized gamma parametric function on a survival curve via flexsurv package (restricted mean t = 365):
library(dplyr)
library(ggplot2)
library(survival)
library(survminer)
library(survRM2)
library(flexsurv)
DD = c(1:365)
CNSR = rbinom(365, 1, 0.5)
data = as.data.frame(cbind(DD,CNSR))
data$SurvObj <- with(data, Surv(DD, CNSR == 0))
res.ggamma <- tryCatch(flexsurvreg(SurvObj ~ 1, data=data, dist="gengamma"), error = function(e) NA)
I then try to find the restricted mean and the lower and upper confidence intervals:
mean.predicted <- round(rmst_gengamma(365, mu = res.ggamma$res[1, 1], sigma = res.ggamma$res[2, 1], Q = res.ggamma$res[3, 1], start = 0)*12/365, 1)
mean.l95 <- round(rmst_gengamma(365, mu = res.ggamma$res[1, 2], sigma = res.ggamma$res[2, 2], Q = res.ggamma$res[3, 1], start = 0)*12/365, 1)
mean.u95 <- round(rmst_gengamma(365, mu = res.ggamma$res[1, 3], sigma = res.ggamma$res[2, 3], Q = res.ggamma$res[3, 1], start = 0)*12/365, 1)
summary <- paste(mean.predicted, " (", mean.l95, " - ", mean.u95, ")", sep = "")
For some of the summary outputs with certain data sets, I get something normal looking such as the below:
7.5 (6.3 - 8.2)
But some of the outputs on confidence intervals look really wacky:
7.7 (9.9 - 2.8) [this one I get from the randomly generated dataset in the code above]
I am wondering what's going on here? Does it have to do with the shape parameter? It doesn't happen with any other parametric fitting function (exp, gompertz, etc). When I look at the survival curve raw data, it looks fine (all the lower bound and upper bound points contain the mean number).