I am a new user of the 'simsurv' package in R, and I am attempting to replicate a simulation study described in NICE DSU TSD 21 (link: https://www.sheffield.ac.uk/nice-dsu/tsds/flexible-methods-survival-analysis). To do this, I'm attempting to re-produce the same results, and one of these results relies on the specification of a mixture-Weibull model.
When I run this code:
# Load relevant packages
library(survival)
library(simsurv)
# Set sample size for the truth (i.e., a large number)
ss_truth <- 100000
# Sample age and frailty
age <- rnorm(ss_truth, 60, 6)
z <- rnorm(ss_truth, 0, 1)
t_lim_truth <- 80
# Other parameters
dist <- "weib"
lambda1 <- 0.6
lambda2 <- 0.7
gamma1 <- 1.8
gamma2 <- 0.5
pmix <- 0.3
frailty_beta <- 0.5
# Set seed for reproducibility
set.seed(1)
# Adjust for mixture versus non-mixture models
pmixopt <- pmix
mixopt <- TRUE
# Adjust length of lambda for mixture versus non-mixture models
if(is.na(lambda2)==TRUE){
lambda_in <- lambda1
} else {
lambda_in <- c(lambda1, lambda2)
}
# Adjust length of gamma for mixture versus non-mixture models
if(is.na(gamma2)==TRUE){
gamma_in <- gamma1
} else {
gamma_in <- c(gamma1, gamma2)
}
# Generate "truth"
td <- simsurv(dist = dist, lambda = lambda_in, gamma = gamma_in,
pmix = pmixopt, maxt = t_lim_truth, x = data.frame(z, frailty_beta), mixture = mixopt)
I obtain this error:
Error in stats::uniroot(rootfn_surv, survival = survival, x = x_i, betas = betas_i, :
f() values at end points not of opposite sign
Please could someone let me know what I am doing incorrectly to generate this error?
I have tried manually replacing the inputs in the line above with the hard-coded entries, which also didn't seem to help:
td <- simsurv(dist = "weib", lambda = c(0.6, 0.7), gamma = c(1.8, 0.5),
pmix = 0.3, maxt = 80, x = data.frame(z, frailty_beta), mixture = TRUE)
When I try different values for gamma, I have managed to get the code to work (e.g., gamma values of 1.2 and 0.8). Therefore, it might be something to do with these specific values of gamma. However, these values are exactly the same as the values used in the DSU study, so I'm unsure why it isn't working in R, but seemed to work ok in Stata. Please see attached image referencing the original gamma values.
Following an exchange with the package developer, Sam Brilleman, a solution has been identified! In brief, the interval argument needs to be defined such that it is large enough to find a solution within. Below is the line I ran which now works correctly (important part right at the end).
td <- simsurv(dist = dist, lambda = lambda_in, gamma = gamma_in, pmix = pmixopt, maxt = t_lim_truth, x = data.frame(z, frailty_beta), mixture = mixopt, interval = c(1E-20, 1E20))