1 Introduction

This dMrs package is designed to fit survival data to the corresponding manuscript.

req_packs = c("sqldf","relsurv","ggplot2","data.table","dMrs")
for(pack in req_packs){
  
  chk_pack = tryCatch(find.package(pack),
    warning = function(ww){NULL},
    error = function(ee){NULL})
  
  if( !is.null(chk_pack) ){
      library(pack,character.only = TRUE)
      next
  }
  
  stop(sprintf("Install %s",pack))
  
}
Loading required package: gsubfn
Loading required package: proto
Loading required package: RSQLite
Loading required package: survival
Loading required package: date

2 Application

The code below will load relsurv’s working dataset rdata and import Slovenia’s latest ratetable from HMD.

data(rdata)

rdata$sex = ifelse(rdata$sex == 1,"male","female")
rdata[1:3,]
  time cens age    sex year agegr
1 2657    1  68 female 8210 62-70
2 1097    1  63 female 8278 62-70
3 3764    1  60   male 8254 54-61
# Slovenia population death tables
female_fn = "./slovenia_female.txt"
male_fn = "./slovenia_male.txt"
slotab = transrate.hmd(female = female_fn,male = male_fn)
dimnames(slotab)
$age
  [1] "0"   "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11" 
 [13] "12"  "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22"  "23" 
 [25] "24"  "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33"  "34"  "35" 
 [37] "36"  "37"  "38"  "39"  "40"  "41"  "42"  "43"  "44"  "45"  "46"  "47" 
 [49] "48"  "49"  "50"  "51"  "52"  "53"  "54"  "55"  "56"  "57"  "58"  "59" 
 [61] "60"  "61"  "62"  "63"  "64"  "65"  "66"  "67"  "68"  "69"  "70"  "71" 
 [73] "72"  "73"  "74"  "75"  "76"  "77"  "78"  "79"  "80"  "81"  "82"  "83" 
 [85] "84"  "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95" 
 [97] "96"  "97"  "98"  "99"  "100" "101" "102" "103" "104" "105" "106" "107"
[109] "108" "109" "110"

$year
 [1] "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990" "1991" "1992"
[11] "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000" "2001" "2002"
[21] "2003" "2004" "2005" "2006" "2007" "2008" "2009" "2010" "2011" "2012"
[31] "2013" "2014" "2015" "2016" "2017" "2018" "2019"

$sex
[1] "male"   "female"
# Hazards calculated per day within age, year, sex strata
slotab[1:3,1:3,]
Rate table with dimension(s): age year sex 
, , sex = male

   year
age         1983         1984         1985
  0 3.890717e-05 4.168568e-05 4.315943e-05
  1 2.985959e-06 3.643855e-06 2.821509e-06
  2 9.036621e-07 9.310505e-07 9.584391e-07

, , sex = female

   year
age         1983         1984         1985
  0 3.452286e-05 3.297061e-05 2.701924e-05
  1 4.164802e-06 1.807623e-06 8.488862e-07
  2 2.492640e-06 7.941114e-07 1.013217e-06
# Subset rdata for years captured by slotab
dim(rdata)
[1] 1040    6
rdata$datediag_yr = as.Date(rdata$year)
rdata$datediag_yr = as.character(rdata$datediag_yr)
rdata$datediag_yr = sapply(rdata$datediag_yr,
  function(xx) strsplit(xx,"-")[[1]][1])
rdata$datediag_yr = rdata$datediag_yr

table(rdata$datediag_yr)

1982 1983 1984 1985 1986 
 168  237  243  220  172 
dimnames(slotab)[[2]]
 [1] "1983" "1984" "1985" "1986" "1987" "1988" "1989" "1990" "1991" "1992"
[11] "1993" "1994" "1995" "1996" "1997" "1998" "1999" "2000" "2001" "2002"
[21] "2003" "2004" "2005" "2006" "2007" "2008" "2009" "2010" "2011" "2012"
[31] "2013" "2014" "2015" "2016" "2017" "2018" "2019"
rdata = rdata[which(rdata$datediag_yr %in% dimnames(slotab)[[2]]),]
dim(rdata)
[1] 872   7

3 Relsurv Approach

test = rs.surv(Surv(time,cens) ~ 1,
  ratetable = slotab,data = rdata,
  rmap = list(age = age*365.241))

str(test)
List of 14
 $ n        : int 872
 $ time     : num [1:804] 9 15 19 23 26 29 30 33 37 40 ...
 $ n.risk   : num [1:804] 872 871 870 869 868 867 866 865 864 862 ...
 $ n.event  : num [1:804] 1 1 1 1 1 1 1 1 2 1 ...
 $ n.censor : num [1:804] 0 0 0 0 0 0 0 0 0 0 ...
 $ std.err  : num [1:804] 0.00115 0.00162 0.00199 0.0023 0.00258 ...
 $ surv     : num [1:804] 1 0.999 0.998 0.998 0.997 ...
 $ lower    : num [1:804] 0.997 0.996 0.995 0.993 0.992 ...
 $ upper    : num [1:804] 1 1 1 1 1 ...
 $ conf.type: chr "log"
 $ conf.int : num 0.95
 $ method   : int 1
 $ call     : language rs.surv(formula = Surv(time, cens) ~ 1, data = rdata, ratetable = slotab,      rmap = list(age = age * 365.241))
 $ type     : chr "right"
 - attr(*, "class")= chr [1:2] "survfit" "rs.surv"
COMP = data.frame(Time = test$time / 365.241,
  SurvEst = test$surv,
  SurvLow = test$lower,
  SurvHigh = test$upper)
plot(COMP$Time,COMP$SurvEst,
  xlab = "Time (yrs)",type = "l",
  ylab = "Net Survival",main = "relsurv method",
  ylim = c(min(COMP$SurvLow),1))
lines(COMP$Time,COMP$SurvLow,lty = 2)
lines(COMP$Time,COMP$SurvHigh,lty = 2)

4 dMrs Approach

Prep wDAT, working dataset’s initial fields.

wDAT = rdata[,c("datediag_yr","time","cens","age","sex")]

wDAT$delta = wDAT$cens

wDAT$datediag_yr = as.integer(wDAT$datediag_yr)

# time in years
wDAT$time = wDAT$time / 365.241

wDAT$age = as.integer(wDAT$age)

wDAT[1:5,]
    datediag_yr       time cens age  sex delta
123        1983 11.4225949    1  56 male     1
124        1983  6.3492324    1  76 male     1
125        1983  0.8295892    1  61 male     1
126        1983 12.8791675    1  52 male     1
127        1983 10.4478960    1  82 male     1

Prep rDAT, the reference data.frame.

mm = fread(file = male_fn,data.table = FALSE)
ff = fread(file = female_fn,data.table = FALSE)

rDAT = rbind(data.frame(sex = "male",mm,stringsAsFactors = FALSE),
  data.frame(sex = "female",ff,stringsAsFactors = FALSE))
rDAT[1:5,]
   sex Year Age      mx      qx   ax     lx   dx    Lx      Tx    ex
1 male 1983   0 0.01429 0.01411 0.12 100000 1411 98759 6677791 66.78
2 male 1983   1 0.00109 0.00109 0.50  98589  107 98535 6579032 66.73
3 male 1983   2 0.00033 0.00033 0.50  98482   32 98465 6480497 65.80
4 male 1983   3 0.00070 0.00070 0.50  98449   69 98415 6382031 64.83
5 male 1983   4 0.00051 0.00051 0.50  98380   50 98355 6283617 63.87
rDAT = rDAT[,c("Year","Age","sex","qx")]
# table(rDAT$Age)
rDAT$Age = ifelse(rDAT$Age == "110+",110,rDAT$Age)
rDAT$Age = as.integer(rDAT$Age)
rDAT[1:5,]
  Year Age  sex      qx
1 1983   0 male 0.01411
2 1983   1 male 0.00109
3 1983   2 male 0.00033
4 1983   3 male 0.00070
5 1983   4 male 0.00051

Perform matching to calculate log density and log cdf.

aa = refData_match(wDAT = wDAT,rDAT = rDAT)
.872 out of 872 done
head(aa)
  log_dens_t2 log_cdf_t2
1   -3.547065 -1.3748526
2   -2.641917 -0.6638762
3   -3.698327 -3.8746426
4   -3.812026 -1.4622169
5   -3.416823 -0.1247204
6   -3.290098 -0.8723565
wDAT = cbind(wDAT,aa)
wDAT[1:3,]
    datediag_yr       time cens age  sex delta log_dens_t2 log_cdf_t2
123        1983 11.4225949    1  56 male     1   -3.547065 -1.3748526
124        1983  6.3492324    1  76 male     1   -2.641917 -0.6638762
125        1983  0.8295892    1  61 male     1   -3.698327 -3.8746426

Prep dMrs inputs.

len1 = 10
len2 = 15
A_range = c(0.4,4)
L_range = quantile(wDAT$time,c(0.5,1))
K_range = c(0.1,2)
T_range = c(0.1,20)

# Less fine grid for alpha/lambda
A_ugrid = log(seq(A_range[1],A_range[2],length.out = len1))
L_ugrid = log(seq(L_range[1],L_range[2],length.out = len1))
# Finer grid for kappa/theta
K_ugrid = log(seq(K_range[1],K_range[2],length.out = len2))
T_ugrid = log(seq(T_range[1],T_range[2],length.out = len2))

param_grid = list(A = A_ugrid,
    L = L_ugrid,K = K_ugrid,T = T_ugrid)
param_grid
$A
 [1] -0.9162907 -0.2231436  0.1823216  0.4700036  0.6931472  0.8754687
 [7]  1.0296194  1.1631508  1.2809338  1.3862944

$L
 [1] 2.060987 2.152588 2.236497 2.313908 2.385754 2.452783 2.515599 2.574702
 [9] 2.630506 2.683359

$K
 [1] -2.30258509 -1.44513486 -0.99039870 -0.67896255 -0.44183275 -0.25029454
 [7] -0.08961216  0.04879016  0.17034537  0.27871340  0.37647757  0.46552935
[13]  0.54729530  0.62287798  0.69314718

$T
 [1] -2.3025851  0.4196497  1.0793809  1.4734545  1.7553918  1.9750726
 [7]  2.1550790  2.3075726  2.4398595  2.5566734  2.6612580  2.7559329
[13]  2.8424146  2.9220088  2.9957323

Run data fit with dMrs’s main function.

res = run_analyses(DATA = wDAT,
    param_grid = param_grid,
    vec_time = seq(0,100,0.5),
    ncores = 1,
    verb = TRUE,
    PLOT = TRUE)

####
Sun Jan 19 01:01:52 2025: Try copula = Independent, dist = Weibull, ...

Sun Jan 19 01:01:52 2025: Assume independence ...
Sun Jan 19 01:01:52 2025: Grid search for initial parameters ...
Sun Jan 19 01:01:52 2025: Start GRID
#ALPHA grid points = 11
#LAMBDA grid points = 10
#KAPPA grid points = 1
#THETA grid points = 1
Num grid points = 110
 110 out of 110 done
Sun Jan 19 01:01:52 2025: End GRID
Sun Jan 19 01:01:52 2025: NR optimization w/ 1 profile point(s) ...
#---# Copula=Independent, Dist=weibull, Prof. point=1 out of 1
iPARS = (0, 2.68336, 0, -inf)
iter=5; LL=-1524.82; diff.LL=0.0142658; diff.PARS=0.0383152; nGRAD=0.0455929; meth=0; reach=0; PARS = (-0.364581, 4.1136, 0, -inf)
iter=10; LL=-1524.81; diff.LL=1.13687e-12; diff.PARS=2.03047e-11; nGRAD=0.00011729; meth=0; reach=5; PARS = (-0.364951, 4.11519, 0, -inf)
No more update
####
Num Iter = 11
Params = (-0.364951, 4.11519, 0, -inf)
LL = -1524.81
GRAD = (-0.000104319, -5.36147e-05, 0, 0)
Convergence Indicators: 
   NormGrad = 0.00011729
   NormIHessGrad = 1.40873e-06
Var = (0.00686217, 0.0565211, 0, 0)
Sun Jan 19 01:01:52 2025: Get covariance ...
02.683359413789010-Inf-1589.91334418337-0.3649506242326234.115193240084370-Inf-1524.810.000117290238296447FALSE
c("alpha1", "lambda1", "kappa1", "theta")c(0.694230928216555, 61.2640517895095, 1, 0)c(0.0575088452067365, 14.5650189952957, 0, 0)c(0.581515662818862, 32.7171391245881, 1, 0)c(0.806946193614247, 89.8109644544309, 1, 0)c(0.590190138044431, 38.4449787938715, 1, 0)c(0.816612394252063, 97.6274187011909, 1, 0)

####
Sun Jan 19 01:01:52 2025: Try copula = Independent, dist = Exp-Weibull, ...

Sun Jan 19 01:01:52 2025: Assume independence ...
Sun Jan 19 01:01:52 2025: Grid search for initial parameters ...
Sun Jan 19 01:01:52 2025: Start GRID
#ALPHA grid points = 11
#LAMBDA grid points = 10
#KAPPA grid points = 16
#THETA grid points = 1
Num grid points = 1760
... 1760 out of 1760 done
Sun Jan 19 01:01:54 2025: End GRID

Sun Jan 19 01:01:54 2025: NR optimization w/ 1 profile point(s) ...
#---# Copula=Independent, Dist=expweibull, Prof. point=1 out of 1
iPARS = (-0.916291, 2.68336, 0.693147, -inf)
No more update
####
Num Iter = 1
Params = (-0.916291, 2.68336, 0.693147, -inf)
LL = -1531.8
GRAD = (31.6921, 31.5899, 83.6041, 0)
Convergence Indicators: 
   NormGrad = 94.8259
   NormIHessGrad = 0.316124
Var = (-0.0103344, -0.0376988, -0.0164437, 0)
Sun Jan 19 01:01:54 2025: No converged solutions! ...

####
Sun Jan 19 01:01:54 2025: Try copula = Clayton, dist = Weibull, ...

Sun Jan 19 01:01:54 2025: Assume dependence: ...
Sun Jan 19 01:01:54 2025: theta = MLE ...
Sun Jan 19 01:01:54 2025: Grid search for initial parameters ...
Sun Jan 19 01:01:54 2025: Start GRID
#ALPHA grid points = 11
#LAMBDA grid points = 10
#KAPPA grid points = 1
#THETA grid points = 16
Num grid points = 1760
... 1760 out of 1760 done
Sun Jan 19 01:01:56 2025: End GRID

Sun Jan 19 01:01:56 2025: NR optimization w/ 1 profile point(s) ...
#---# Copula=Clayton, Dist=weibull, Prof. point=1 out of 1
iPARS = (0, 2.68336, 0, 0.41965)
iter=5; LL=-1525.22; diff.LL=0.26815; diff.PARS=1; nGRAD=2.16244; meth=0; reach=0; PARS = (-0.343036, 4.04924, 0, -4.28443)
iter=10; LL=-1524.86; diff.LL=0.0108357; diff.PARS=0.25; nGRAD=0.0405005; meth=0; reach=0; PARS = (-0.364222, 4.11026, 0, -5.37775)
iter=15; LL=-1524.84; diff.LL=0.00228082; diff.PARS=0.0936096; nGRAD=0.0255952; meth=0; reach=0; PARS = (-0.364416, 4.1119, 0, -5.86946)
iter=20; LL=-1524.82; diff.LL=5.50145e-05; diff.PARS=0.00712546; nGRAD=0.00946472; meth=0; reach=0; PARS = (-0.364761, 4.11407, 0, -6.97022)
iter=25; LL=-1524.82; diff.LL=3.27969e-06; diff.PARS=0.000453249; nGRAD=0.0156254; meth=0; reach=1; PARS = (-0.364899, 4.11465, 0, -7.05907)
iter=30; LL=-1524.82; diff.LL=1.82724e-05; diff.PARS=0.00306553; nGRAD=0.0148505; meth=0; reach=1; PARS = (-0.365051, 4.11519, 0, -7.3019)
iter=35; LL=-1524.82; diff.LL=7.28853e-05; diff.PARS=0.0135864; nGRAD=0.00619153; meth=0; reach=0; PARS = (-0.364788, 4.11429, 0, -7.33468)
iter=40; LL=-1524.81; diff.LL=0.000155857; diff.PARS=0.0322258; nGRAD=0.00842428; meth=0; reach=0; PARS = (-0.364929, 4.11483, 0, -7.44394)
iter=45; LL=-1524.81; diff.LL=0.000854395; diff.PARS=0.260748; nGRAD=0.06494; meth=0; reach=0; PARS = (-0.363941, 4.11183, 0, -7.87016)
iter=50; LL=-1524.81; diff.LL=0.000113742; diff.PARS=0.043695; nGRAD=0.00962198; meth=0; reach=0; PARS = (-0.364749, 4.11439, 0, -8.28866)
iter=55; LL=-1524.81; diff.LL=1.35936e-05; diff.PARS=0.0067181; nGRAD=0.00339567; meth=0; reach=0; PARS = (-0.364894, 4.11489, 0, -8.31562)
iter=60; LL=-1524.81; diff.LL=1.09125e-05; diff.PARS=0.00559048; nGRAD=0.0037728; meth=0; reach=0; PARS = (-0.364957, 4.11509, 0, -8.33561)
iter=65; LL=-1524.81; diff.LL=1.86089e-05; diff.PARS=0.00955709; nGRAD=0.00364773; meth=0; reach=0; PARS = (-0.364886, 4.11487, 0, -8.35142)
iter=70; LL=-1524.81; diff.LL=6.09548e-06; diff.PARS=0.00309094; nGRAD=0.00227521; meth=0; reach=3; PARS = (-0.364929, 4.11499, 0, -8.36763)
iter=75; LL=-1524.81; diff.LL=8.72633e-06; diff.PARS=0.0046756; nGRAD=0.00186553; meth=0; reach=3; PARS = (-0.36491, 4.11493, 0, -8.38638)
iter=80; LL=-1524.81; diff.LL=2.01537e-08; diff.PARS=3.56576e-06; nGRAD=0.0160512; meth=0; reach=2; PARS = (-0.365172, 4.1158, 0, -8.51842)
iter=85; LL=-1524.81; diff.LL=5.25293e-06; diff.PARS=0.00323047; nGRAD=0.00266432; meth=0; reach=7; PARS = (-0.364914, 4.11498, 0, -8.52886)
iter=90; LL=-1524.81; diff.LL=2.67506e-06; diff.PARS=0.00166774; nGRAD=0.00194806; meth=0; reach=12; PARS = (-0.364923, 4.115, 0, -8.54045)
Optimization criteria met
####
Num Iter = 93
Params = (-0.36493, 4.11502, 0, -8.54285)
LL = -1524.81
GRAD = (-0.000122009, -0.00104369, 0, -0.00159193)
Convergence Indicators: 
   NormGrad = 0.00190747
   NormIHessGrad = 0.00329959
Var = (0.00695345, 0.057652, 0, 2.09317)
Sun Jan 19 01:02:02 2025: Get covariance ...
02.6833594137890100.419649743100121-1567.33915264465-0.3649297168259994.115017798610590-8.54285470023963-1524.8120.00190746698854164FALSE
c("alpha1", "lambda1", "kappa1", "theta")c(0.694245442936594, 61.2533044767636, 1, 0.000194932989852757)c(0.0578912931614559, 14.7074211844081, 0, 0.000282025206382938)c(0.58078059332169, 32.4272886498623, 1, -0.000357826257390278)c(0.807710292551497, 90.0793203036648, 1, 0.000747692237095792)c(0.589567574777206, 38.2603594210566, 1, 1.1438886049511e-05)c(0.81750889237822, 98.0640894674452, 1, 0.0033219030566844)

####
Sun Jan 19 01:02:03 2025: Try copula = Clayton, dist = Exp-Weibull, ...

Sun Jan 19 01:02:03 2025: Assume dependence: ...
Sun Jan 19 01:02:03 2025: theta = MLE ...
Sun Jan 19 01:02:03 2025: Grid search for initial parameters ...
Sun Jan 19 01:02:03 2025: Start GRID
#ALPHA grid points = 11
#LAMBDA grid points = 10
#KAPPA grid points = 16
#THETA grid points = 16
Num grid points = 28160
.......... 5000.......... 10000 out of 28160 done
.......... 15000.......... 20000 out of 28160 done
.......... 25000...... 28160 out of 28160 done
Sun Jan 19 01:02:35 2025: End GRID

Sun Jan 19 01:02:35 2025: NR optimization w/ 2 profile point(s) ...
#---# Copula=Clayton, Dist=expweibull, Prof. point=1 out of 2
iPARS = (-0.916291, 2.68336, 0.693147, -2.30259)
No more update
####
Num Iter = 1
Params = (-0.916291, 2.68336, 0.693147, -2.30259)
LL = -1530.67
GRAD = (30.6387, 27.1626, 74.4244, 0.846512)
Convergence Indicators: 
   NormGrad = 84.9485
   NormIHessGrad = 1.05686
Var = (-0.0124412, -0.0586044, -0.0204214, -1.45715)
#---# Copula=Clayton, Dist=expweibull, Prof. point=2 out of 2
iPARS = (-0.223144, 2.68336, 0.278713, 0)
No more update
####
Num Iter = 1
Params = (-0.223144, 2.68336, 0.278713, 0)
LL = -1550.64
GRAD = (-52.2677, 44.7771, 2.73584, -2.50489)
Convergence Indicators: 
   NormGrad = 68.925
   NormIHessGrad = 1.17749
Var = (-0.0105215, -0.00903437, -0.015051, 0.187011)
Sun Jan 19 01:02:37 2025: No converged solutions! ...

####
Sun Jan 19 01:02:37 2025: Try copula = Gumbel, dist = Weibull, ...

Sun Jan 19 01:02:37 2025: Assume dependence: ...
Sun Jan 19 01:02:37 2025: theta = MLE ...
Sun Jan 19 01:02:37 2025: Grid search for initial parameters ...
Sun Jan 19 01:02:37 2025: Start GRID
#ALPHA grid points = 11
#LAMBDA grid points = 10
#KAPPA grid points = 1
#THETA grid points = 16
Num grid points = 1760
... 1760 out of 1760 done
Sun Jan 19 01:02:39 2025: End GRID

Sun Jan 19 01:02:39 2025: NR optimization w/ 1 profile point(s) ...
#---# Copula=Gumbel, Dist=weibull, Prof. point=1 out of 1
iPARS = (0, 2.68336, 0, 0)
iter=5; LL=-1525.42; diff.LL=0.444766; diff.PARS=0.937506; nGRAD=0.235658; meth=0; reach=0; PARS = (-0.359236, 4.09115, 0, -4.34478)
iter=10; LL=-1524.83; diff.LL=0.00490524; diff.PARS=0.352634; nGRAD=0.0190056; meth=0; reach=0; PARS = (-0.364631, 4.11394, 0, -7.028)
iter=15; LL=-1524.82; diff.LL=0.000723369; diff.PARS=0.10248; nGRAD=0.00686415; meth=0; reach=0; PARS = (-0.3648, 4.11451, 0, -7.54197)
iter=20; LL=-1524.82; diff.LL=0.000872367; diff.PARS=0.161166; nGRAD=0.00844451; meth=0; reach=0; PARS = (-0.364801, 4.11461, 0, -7.83735)
iter=25; LL=-1524.81; diff.LL=1.40578e-05; diff.PARS=0.00288327; nGRAD=0.00630609; meth=0; reach=2; PARS = (-0.364806, 4.11459, 0, -7.85752)
iter=30; LL=-1524.81; diff.LL=2.33073e-05; diff.PARS=0.00511904; nGRAD=0.0100628; meth=0; reach=0; PARS = (-0.364996, 4.11528, 0, -7.98445)
iter=35; LL=-1524.81; diff.LL=3.24447e-06; diff.PARS=0.000856284; nGRAD=0.00464296; meth=0; reach=1; PARS = (-0.364872, 4.11485, 0, -8.11161)
iter=40; LL=-1524.81; diff.LL=7.35308e-06; diff.PARS=0.00214414; nGRAD=0.00472779; meth=0; reach=1; PARS = (-0.364918, 4.11502, 0, -8.2136)
iter=45; LL=-1524.81; diff.LL=0.000128156; diff.PARS=0.0424082; nGRAD=0.00358772; meth=0; reach=0; PARS = (-0.364903, 4.11497, 0, -8.36078)
iter=50; LL=-1524.81; diff.LL=1.6813e-06; diff.PARS=0.000685932; nGRAD=0.00455423; meth=0; reach=1; PARS = (-0.364952, 4.11516, 0, -8.56591)
iter=55; LL=-1524.81; diff.LL=0.000117103; diff.PARS=0.0524827; nGRAD=0.00827748; meth=0; reach=0; PARS = (-0.365043, 4.11545, 0, -8.65155)
iter=60; LL=-1524.81; diff.LL=1.73855e-05; diff.PARS=0.0103164; nGRAD=0.0025298; meth=0; reach=0; PARS = (-0.364925, 4.11508, 0, -8.92867)
iter=65; LL=-1524.81; diff.LL=2.85778e-05; diff.PARS=0.0177819; nGRAD=0.0026209; meth=0; reach=0; PARS = (-0.364913, 4.11504, 0, -8.9814)
iter=70; LL=-1524.81; diff.LL=1.60999e-05; diff.PARS=0.010547; nGRAD=0.00157054; meth=0; reach=0; PARS = (-0.364927, 4.11508, 0, -9.02594)
iter=75; LL=-1524.81; diff.LL=4.43144e-06; diff.PARS=0.00346161; nGRAD=0.00447832; meth=0; reach=1; PARS = (-0.364964, 4.11524, 0, -9.2094)
iter=80; LL=-1524.81; diff.LL=3.92854e-07; diff.PARS=0.000316689; nGRAD=0.00135826; meth=0; reach=2; PARS = (-0.364922, 4.11507, 0, -9.22849)
iter=85; LL=-1524.81; diff.LL=9.24128e-06; diff.PARS=0.00777069; nGRAD=0.00148635; meth=0; reach=0; PARS = (-0.364941, 4.11513, 0, -9.27101)
iter=90; LL=-1524.81; diff.LL=5.12218e-06; diff.PARS=0.00448773; nGRAD=0.00168266; meth=0; reach=4; PARS = (-0.364914, 4.11503, 0, -9.31396)
iter=95; LL=-1524.81; diff.LL=3.29298e-05; diff.PARS=0.030586; nGRAD=0.00284653; meth=0; reach=0; PARS = (-0.364887, 4.11494, 0, -9.37928)
iter=100; LL=-1524.81; diff.LL=2.89599e-05; diff.PARS=0.0284176; nGRAD=0.00159831; meth=0; reach=0; PARS = (-0.36493, 4.11509, 0, -9.43851)
iter=105; LL=-1524.81; diff.LL=1.28676e-06; diff.PARS=0.00129304; nGRAD=0.00129885; meth=0; reach=2; PARS = (-0.364946, 4.11515, 0, -9.45401)
iter=110; LL=-1524.81; diff.LL=1.3031e-05; diff.PARS=0.013437; nGRAD=0.00174012; meth=0; reach=0; PARS = (-0.364908, 4.11502, 0, -9.47997)
iter=115; LL=-1524.81; diff.LL=2.95512e-06; diff.PARS=0.00318446; nGRAD=0.00154572; meth=0; reach=2; PARS = (-0.364944, 4.11515, 0, -9.52054)
iter=120; LL=-1524.81; diff.LL=3.50795e-06; diff.PARS=0.0038511; nGRAD=0.000999176; meth=0; reach=2; PARS = (-0.364937, 4.11513, 0, -9.53946)
iter=125; LL=-1524.81; diff.LL=4.29556e-05; diff.PARS=0.0495436; nGRAD=0.00672694; meth=0; reach=0; PARS = (-0.36488, 4.11498, 0, -9.60671)
iter=130; LL=-1524.81; diff.LL=3.89602e-06; diff.PARS=0.00475033; nGRAD=0.000843891; meth=0; reach=3; PARS = (-0.364938, 4.11513, 0, -9.64472)
iter=135; LL=-1524.81; diff.LL=5.56522e-06; diff.PARS=0.00743135; nGRAD=0.000819621; meth=0; reach=0; PARS = (-0.364934, 4.11511, 0, -9.74889)
iter=140; LL=-1524.81; diff.LL=7.31939e-09; diff.PARS=1.03345e-05; nGRAD=0.000887789; meth=0; reach=1; PARS = (-0.364942, 4.11515, 0, -9.79488)
iter=145; LL=-1524.81; diff.LL=3.66894e-07; diff.PARS=0.000544256; nGRAD=0.000748335; meth=0; reach=1; PARS = (-0.364939, 4.11514, 0, -9.83907)
iter=150; LL=-1524.81; diff.LL=1.38975e-07; diff.PARS=0.000192945; nGRAD=0.00211247; meth=0; reach=1; PARS = (-0.364906, 4.11502, 0, -9.91011)
iter=155; LL=-1524.81; diff.LL=4.15586e-05; diff.PARS=0.0714183; nGRAD=0.00608527; meth=0; reach=0; PARS = (-0.365041, 4.11549, 0, -9.99336)
iter=160; LL=-1524.81; diff.LL=3.14753e-07; diff.PARS=0.000448372; nGRAD=0.00461267; meth=0; reach=1; PARS = (-0.365013, 4.1154, 0, -10.0395)
iter=165; LL=-1524.81; diff.LL=1.93038e-06; diff.PARS=0.00350327; nGRAD=0.000743859; meth=0; reach=1; PARS = (-0.364943, 4.11516, 0, -10.0545)
iter=170; LL=-1524.81; diff.LL=3.93386e-07; diff.PARS=0.000733345; nGRAD=0.00071975; meth=0; reach=3; PARS = (-0.364943, 4.11516, 0, -10.0735)
iter=175; LL=-1524.81; diff.LL=1.97834e-06; diff.PARS=0.00382138; nGRAD=0.000652825; meth=0; reach=1; PARS = (-0.364945, 4.11516, 0, -10.1006)
iter=180; LL=-1524.81; diff.LL=9.06025e-07; diff.PARS=0.00177848; nGRAD=0.000684617; meth=0; reach=3; PARS = (-0.364934, 4.11512, 0, -10.1168)
iter=185; LL=-1524.81; diff.LL=5.03537e-07; diff.PARS=0.000996948; nGRAD=0.000575152; meth=0; reach=8; PARS = (-0.364943, 4.11516, 0, -10.1243)
iter=190; LL=-1524.81; diff.LL=5.56865e-06; diff.PARS=0.0110993; nGRAD=0.000782013; meth=0; reach=0; PARS = (-0.364935, 4.11513, 0, -10.1399)
iter=195; LL=-1524.81; diff.LL=3.72539e-07; diff.PARS=0.000731757; nGRAD=0.00125031; meth=0; reach=1; PARS = (-0.364922, 4.11508, 0, -10.1607)
iter=200; LL=-1524.81; diff.LL=2.32239e-06; diff.PARS=0.00486576; nGRAD=0.00113166; meth=0; reach=2; PARS = (-0.364924, 4.11509, 0, -10.1788)
####
Num Iter = 201
Params = (-0.364924, 4.11509, 0, -10.1788)
LL = -1524.81
GRAD = (-0.000566615, 0.000854607, 0, -0.000478803)
Convergence Indicators: 
   NormGrad = 0.00113166
   NormIHessGrad = 0.00451256
Var = (0.00680882, 0.0559227, 0, 9.26726)
Sun Jan 19 01:02:54 2025: Get covariance ...
02.6833594137890100-1565.92618810742-0.364924411774364.115092274826450-10.178835752156-1524.810.00113166168428236FALSE
c("alpha1", "lambda1", "kappa1", "theta")c(0.694249125954288, 61.2578665609712, 1, 1.00003796538433)c(0.0572863756056385, 14.4862485931599, 0, 0.000115574879391569)c(0.581969892962402, 32.8653410472839, 1, 0.999811442783202)c(0.806528358946173, 89.6503920746586, 1, 1.00026448798545)c(0.59057892488115, 38.5362888861027, 1, 1.00000009730728)c(0.816117590016081, 97.3764294401225, 1, 1.01481256454601)

####
Sun Jan 19 01:02:54 2025: Try copula = Gumbel, dist = Exp-Weibull, ...

Sun Jan 19 01:02:54 2025: Assume dependence: ...
Sun Jan 19 01:02:54 2025: theta = MLE ...
Sun Jan 19 01:02:54 2025: Grid search for initial parameters ...
Sun Jan 19 01:02:54 2025: Start GRID
#ALPHA grid points = 11
#LAMBDA grid points = 10
#KAPPA grid points = 16
#THETA grid points = 16
Num grid points = 28160
.......... 5000.......... 10000 out of 28160 done
.......... 15000.......... 20000 out of 28160 done
.......... 25000...... 28160 out of 28160 done
Sun Jan 19 01:03:29 2025: End GRID

Sun Jan 19 01:03:29 2025: NR optimization w/ 2 profile point(s) ...
#---# Copula=Gumbel, Dist=expweibull, Prof. point=1 out of 2
iPARS = (-0.916291, 2.68336, 0.693147, -2.30259)
No more update
####
Num Iter = 1
Params = (-0.916291, 2.68336, 0.693147, -2.30259)
LL = -1531.56
GRAD = (32.9268, 27.5153, 76.4866, -0.0741572)
Convergence Indicators: 
   NormGrad = 87.701
   NormIHessGrad = 3.75999
Var = (-0.0122108, -0.0231552, -0.019186, 2.69044)
#---# Copula=Gumbel, Dist=expweibull, Prof. point=2 out of 2
iPARS = (-0.223144, 2.68336, 0.376478, 0)
No more update
####
Num Iter = 1
Params = (-0.223144, 2.68336, 0.376478, 0)
LL = -1549.28
GRAD = (-50.5325, 26.4988, -13.5553, -4.95901)
Convergence Indicators: 
   NormGrad = 58.8563
   NormIHessGrad = 1.07062
Var = (-0.0156799, -0.0157895, -0.0253651, 0.124703)
Sun Jan 19 01:03:31 2025: No converged solutions! ...

Check dMrs output

# See all solutions
OO = opt_sum(OPT = res)
OO
  IDX      COPULA    DIST    alpha1  lambda1 kappa1       theta        LL NP
1   1 Independent weibull 0.6942309 61.26405      1 0.000000000 -1524.810  2
2   2     Clayton weibull 0.6942454 61.25330      1 0.000194933 -1524.812  3
3   3      Gumbel weibull 0.6942491 61.25787      1 1.000037965 -1524.810  3
        BIC       POST
1 -3063.162 0.99771623
2 -3069.936 0.00113960
3 -3069.932 0.00114417
# Select best model
idx = which(OO$BIC == max(OO$BIC))
idx
[1] 1
# MLEs (unconstrained)
res[[idx]]$RES$out
         PARS        EST         SE      lowCI     highCI
1  log_alpha1 -0.3649506 0.08283821 -0.5273105 -0.2025907
2 log_lambda1  4.1151932 0.23774169  3.6492281  4.5811584
3  unc_kappa1  0.0000000 0.00000000  0.0000000  0.0000000
4   log_theta       -Inf 0.00000000       -Inf       -Inf
# MLEs (constrained)
res[[idx]]$RES$cout
     PARS        EST          SE      lowCI     highCI    lowCI_2   highCI_2
1  alpha1  0.6942309  0.05750885  0.5815157  0.8069462  0.5901901  0.8166124
2 lambda1 61.2640518 14.56501900 32.7171391 89.8109645 38.4449788 97.6274187
3  kappa1  1.0000000  0.00000000  1.0000000  1.0000000  1.0000000  1.0000000
4   theta  0.0000000  0.00000000  0.0000000  0.0000000  0.0000000  0.0000000
# Log-likelihood
res[[idx]]$RES$LL
[1] -1524.81
# Gradient
res[[idx]]$RES$GRAD
   log_alpha1   log_lambda1    unc_kappa1     log_theta 
-1.043190e-04 -5.361471e-05  0.000000e+00  0.000000e+00 
# Hessian
res[[idx]]$RES$HESS
            log_alpha1 log_lambda1 unc_kappa1 log_theta
log_alpha1   -390.6644  -107.78422          0         0
log_lambda1  -107.7842   -47.43015          0         0
unc_kappa1      0.0000     0.00000          0         0
log_theta       0.0000     0.00000          0         0
# Covariance matrix
res[[idx]]$RES$COVAR
              log_alpha1 log_lambda1 unc_kappa1 log_theta
log_alpha1   0.006862169 -0.01559416          0         0
log_lambda1 -0.015594163  0.05652111          0         0
unc_kappa1   0.000000000  0.00000000          0         0
log_theta    0.000000000  0.00000000          0         0

5 Net-survival

# Predicted survivals
res[[idx]]$PRED[1:3,]
  time    log_F1         F1 log_surv      surv          SE  low_surv high_surv
1  0.0      -Inf 0.00000000   0.0000 1.0000000 0.000000000 1.0000000 1.0000000
2  0.5 -3.355798 0.03488151  -0.0355 0.9651227 0.006074747 0.9532162 0.9770292
3  1.0 -2.885480 0.05582797  -0.0574 0.9442163 0.007933013 0.9286676 0.9597650
         AA log_low_surv log_high_surv low_surv2 high_surv2
1 0.0000000   0.00000000    0.00000000 1.0000000  1.0000000
2 0.3474733  -0.05024977   -0.02507972 0.9509919  0.9752322
3 0.2866663  -0.07645563   -0.04309375 0.9263940  0.9578216
plot_SURVs(run_ANA = res[idx],MULTIPLE = FALSE)

Compare dMrs vs relsurv

tmp_pred = res[[idx]]$PRED

out = sqldf("
select
    COMP.*,
    'Pohar-Perme' as Method
from
    COMP

union

select
    DMRS.time as Time,
    DMRS.surv as SurvEst,
    DMRS.low_surv2 as SurvLow,
    DMRS.high_surv2 as SurvHigh,
    'dMrs' as Method
from
    tmp_pred as DMRS
")

my_themes = theme(text = element_text(size = 28),
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5),
    panel.background = element_blank(),
    panel.grid.major = element_line(colour = "grey50",
        linewidth = 0.5,linetype = "dotted"),
    panel.border = element_rect(colour = "black",
        fill = NA,linewidth = 1),
    legend.key.width = unit(1.5, "cm"),
    legend.key.size = unit(0.5, "inch"),
    legend.text = element_text(size = 20))
        
ggplot(data = out,
    mapping = aes(x = Time,y = SurvEst,group = Method,fill = Method)) +
    geom_line(linewidth = 1.25,alpha = 1,
        aes(color = Method),show.legend = FALSE) +
    geom_ribbon(mapping = aes(ymin = SurvLow,
        ymax = SurvHigh),alpha = 0.5) +
    ylim(c(0.4,1)) + xlim(0,20) +
    xlab("Time (yrs)") + ylab("Net Survival") +
    my_themes

6 Session information

R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale:
[1] LC_COLLATE=C                          
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] dMrs_1.0.0        data.table_1.14.8 ggplot2_3.4.0     relsurv_2.2-9    
 [5] date_1.2-42       survival_3.4-0    sqldf_0.4-11      RSQLite_2.3.1    
 [9] gsubfn_0.7        proto_1.0.0      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.10         ADGofTest_0.3       mvtnorm_1.1-3      
 [4] lattice_0.20-45     gtools_3.9.4        digest_0.6.31      
 [7] utf8_1.2.3          gmp_0.7-2           R6_2.5.1           
[10] chron_2.3-61        stats4_4.2.2        pcaPP_2.0-3        
[13] evaluate_0.20       highr_0.10          pillar_1.9.0       
[16] gplots_3.1.3        rlang_1.1.1         jquerylib_0.1.4    
[19] blob_1.2.4          Matrix_1.5-1        rmarkdown_2.20     
[22] labeling_0.4.3      splines_4.2.2       bit_4.0.5          
[25] munsell_0.5.0       compiler_4.2.2      numDeriv_2016.8-1.1
[28] xfun_0.42           pkgconfig_2.0.3     gsl_2.1-8          
[31] htmltools_0.5.4     tcltk_4.2.2         tidyselect_1.2.0   
[34] gridExtra_2.3       tibble_3.2.1        stabledist_0.7-1   
[37] viridisLite_0.4.2   fansi_1.0.4         pspline_1.0-19     
[40] dplyr_1.1.2         withr_2.5.0         bitops_1.0-7       
[43] grid_4.2.2          jsonlite_1.8.5      gtable_0.3.4       
[46] lifecycle_1.0.3     DBI_1.1.3           magrittr_2.0.3     
[49] scales_1.2.1        KernSmooth_2.23-20  cli_3.6.1          
[52] cachem_1.0.8        farver_2.1.1        Rmpfr_0.9-3        
[55] viridis_0.6.2       bslib_0.4.2         generics_0.1.3     
[58] vctrs_0.6.2         tools_4.2.2         bit64_4.0.5        
[61] copula_1.1-2        glue_1.6.2          fastmap_1.1.1      
[64] yaml_2.3.7          colorspace_2.1-0    caTools_1.18.2     
[67] memoise_2.0.1       knitr_1.42          sass_0.4.5

7 References