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

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 1982-06-24 62-70
2 1097    1  63 female 1982-08-31 62-70
3 3764    1  60   male 1982-08-07 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

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)

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)

####
Sat Mar 22 07:18:44 2025: Try copula = Independent, dist = Weibull, ...

Sat Mar 22 07:18:44 2025: Assume independence ...

Sat Mar 22 07:18:44 2025: Grid search for initial parameters ...

Sat Mar 22 07:18:44 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

Sat Mar 22 07:18:45 2025: End GRID

Sat Mar 22 07:18:45 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.0142913; diff.PARS=0.0385319; nGRAD=0.0378659; meth=0; reach=0; PARS = (-0.364617, 4.11379, 0, -inf)
iter=10; LL=-1524.81; diff.LL=2.27374e-13; diff.PARS=6.62323e-12; nGRAD=3.05054e-07; meth=0; reach=5; PARS = (-0.364951, 4.11519, 0, -inf)
No more update
####
Num Iter = 14
Params = (-0.364951, 4.11519, 0, -inf)
LL = -1524.81
GRAD = (-1.18234e-06, -2.27374e-07, 0, 0)
Convergence Indicators: 
   NormGrad = 1.20401e-06
   NormIHessGrad = 7.21501e-09
Var = (0.00686627, 0.0566411, 0, 0)

Sat Mar 22 07:18:45 2025: Get covariance ...

02.683359413789010-Inf-1589.91334418337-0.36495050417084.115191847211570-Inf-1524.811.20400748481859e-06FALSE

c("alpha1", "lambda1", "kappa1", "theta")c(0.69423101156719, 61.2639664565374, 1, 0)c(0.0575260437845105, 14.5804554567041, 0, 0)c(0.581482037576475, 32.6867988832069, 1, 0)c(0.806979985557905, 89.841134029868, 1, 0)c(0.590161564242663, 38.4259191554175, 1, 0)c(0.816652128201682, 97.6755707731349, 1, 0)

####
Sat Mar 22 07:18:45 2025: Try copula = Independent, dist = Exp-Weibull, ...

Sat Mar 22 07:18:45 2025: Assume independence ...

Sat Mar 22 07:18:45 2025: Grid search for initial parameters ...

Sat Mar 22 07:18:45 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

Sat Mar 22 07:18:45 2025: End GRID

Sat Mar 22 07:18:45 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.316201
Var = (-0.0103416, -0.0377196, -0.0164523, 0)

Sat Mar 22 07:18:45 2025: No converged solutions! ...

####
Sat Mar 22 07:18:45 2025: Try copula = Clayton, dist = Weibull, ...

Sat Mar 22 07:18:45 2025: Assume dependence: ...

Sat Mar 22 07:18:45 2025: theta = MLE ...

Sat Mar 22 07:18:45 2025: Grid search for initial parameters ...

Sat Mar 22 07:18:45 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

Sat Mar 22 07:18:46 2025: End GRID

Sat Mar 22 07:18:46 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.17; diff.LL=0.233356; diff.PARS=0.95717; nGRAD=0.804172; meth=0; reach=0; PARS = (-0.362994, 4.10715, 0, -4.2269)
iter=10; LL=-1524.83; diff.LL=0.0014019; diff.PARS=0.0677835; nGRAD=0.028504; meth=0; reach=0; PARS = (-0.364653, 4.11301, 0, -6.03426)
iter=15; LL=-1524.82; diff.LL=0.000189503; diff.PARS=0.0138843; nGRAD=0.0182257; meth=0; reach=0; PARS = (-0.364839, 4.11386, 0, -6.41294)
iter=20; LL=-1524.82; diff.LL=0.000834955; diff.PARS=0.0892889; nGRAD=0.018853; meth=0; reach=0; PARS = (-0.364881, 4.11448, 0, -6.82646)
iter=25; LL=-1524.82; diff.LL=0.000290977; diff.PARS=0.0423673; nGRAD=0.0114016; meth=0; reach=0; PARS = (-0.364747, 4.11413, 0, -7.12106)
iter=30; LL=-1524.82; diff.LL=0.000290615; diff.PARS=0.0503968; nGRAD=0.0130991; meth=0; reach=0; PARS = (-0.364938, 4.11469, 0, -7.28006)
iter=35; LL=-1524.81; diff.LL=8.64693e-06; diff.PARS=0.00223018; nGRAD=0.0159636; meth=0; reach=1; PARS = (-0.364882, 4.11483, 0, -7.65433)
iter=40; LL=-1524.81; diff.LL=7.65586e-06; diff.PARS=0.00209933; nGRAD=0.00381238; meth=0; reach=1; PARS = (-0.364887, 4.11473, 0, -7.71737)
iter=45; LL=-1524.81; diff.LL=2.65642e-05; diff.PARS=0.00730097; nGRAD=0.00813958; meth=0; reach=0; PARS = (-0.364769, 4.11439, 0, -7.81518)
iter=50; LL=-1524.81; diff.LL=3.98524e-05; diff.PARS=0.0125517; nGRAD=0.00394193; meth=0; reach=0; PARS = (-0.364919, 4.11487, 0, -7.85972)
iter=55; LL=-1524.81; diff.LL=1.56537e-05; diff.PARS=0.00164498; nGRAD=0.0247852; meth=0; reach=1; PARS = (-0.364504, 4.11357, 0, -7.97989)
iter=60; LL=-1524.81; diff.LL=2.8224e-05; diff.PARS=0.0102841; nGRAD=0.00529869; meth=0; reach=0; PARS = (-0.364883, 4.11483, 0, -8.05981)
iter=65; LL=-1524.81; diff.LL=7.84174e-05; diff.PARS=0.0318346; nGRAD=0.00927572; meth=0; reach=0; PARS = (-0.364982, 4.11519, 0, -8.10808)
iter=70; LL=-1524.81; diff.LL=5.86835e-06; diff.PARS=0.00244917; nGRAD=0.00271871; meth=0; reach=2; PARS = (-0.364885, 4.11481, 0, -8.1352)
iter=75; LL=-1524.81; diff.LL=9.10569e-05; diff.PARS=0.0455383; nGRAD=0.0223242; meth=0; reach=0; PARS = (-0.365212, 4.11604, 0, -8.21097)
iter=80; LL=-1524.81; diff.LL=1.45111e-05; diff.PARS=0.00943657; nGRAD=0.00342286; meth=0; reach=0; PARS = (-0.36492, 4.11496, 0, -8.58583)
iter=85; LL=-1524.81; diff.LL=9.04297e-06; diff.PARS=0.0060557; nGRAD=0.00451498; meth=0; reach=0; PARS = (-0.364941, 4.11509, 0, -8.59767)
iter=90; LL=-1524.81; diff.LL=5.38822e-05; diff.PARS=0.0447334; nGRAD=0.0349745; meth=0; reach=0; PARS = (-0.365163, 4.11609, 0, -8.64807)
iter=95; LL=-1524.81; diff.LL=1.54866e-06; diff.PARS=0.00105039; nGRAD=0.00798328; meth=0; reach=2; PARS = (-0.364876, 4.11478, 0, -8.68612)
iter=100; LL=-1524.81; diff.LL=1.71485e-06; diff.PARS=0.00124459; nGRAD=0.00138718; meth=0; reach=4; PARS = (-0.364916, 4.11498, 0, -8.69718)
iter=105; LL=-1524.81; diff.LL=5.44619e-06; diff.PARS=0.00408193; nGRAD=0.00133929; meth=0; reach=2; PARS = (-0.36492, 4.115, 0, -8.7217)
iter=110; LL=-1524.81; diff.LL=6.38264e-07; diff.PARS=0.000441659; nGRAD=0.00225744; meth=0; reach=7; PARS = (-0.364947, 4.1151, 0, -8.72638)
iter=115; LL=-1524.81; diff.LL=7.00424e-06; diff.PARS=0.00541221; nGRAD=0.00351999; meth=0; reach=0; PARS = (-0.36497, 4.11516, 0, -8.7364)
iter=120; LL=-1524.81; diff.LL=1.56103e-06; diff.PARS=0.00120649; nGRAD=0.00214321; meth=0; reach=5; PARS = (-0.364927, 4.11504, 0, -8.74966)
iter=125; LL=-1524.81; diff.LL=5.61715e-06; diff.PARS=0.00451461; nGRAD=0.00430436; meth=0; reach=2; PARS = (-0.364854, 4.11478, 0, -8.77794)
iter=130; LL=-1524.81; diff.LL=3.39617e-06; diff.PARS=0.00272035; nGRAD=0.00253897; meth=0; reach=4; PARS = (-0.364926, 4.11504, 0, -8.80434)
iter=135; LL=-1524.81; diff.LL=8.51199e-05; diff.PARS=0.08292; nGRAD=0.0206622; meth=0; reach=0; PARS = (-0.364582, 4.11389, 0, -8.89487)
iter=140; LL=-1524.81; diff.LL=1.47454e-08; diff.PARS=1.38588e-05; nGRAD=0.00317088; meth=0; reach=2; PARS = (-0.364917, 4.11503, 0, -8.94657)
iter=145; LL=-1524.81; diff.LL=3.08442e-06; diff.PARS=0.00293849; nGRAD=0.00168014; meth=0; reach=1; PARS = (-0.364919, 4.11503, 0, -8.96757)
iter=150; LL=-1524.81; diff.LL=1.37945e-06; diff.PARS=0.00132082; nGRAD=0.00108888; meth=0; reach=6; PARS = (-0.364924, 4.11503, 0, -8.96986)
iter=155; LL=-1524.81; diff.LL=4.06321e-08; diff.PARS=3.92815e-05; nGRAD=0.00289312; meth=0; reach=11; PARS = (-0.364923, 4.11505, 0, -8.97614)
iter=160; LL=-1524.81; diff.LL=2.51836e-06; diff.PARS=0.00242641; nGRAD=0.00136607; meth=0; reach=2; PARS = (-0.364913, 4.115, 0, -8.98795)
iter=165; LL=-1524.81; diff.LL=1.1717e-06; diff.PARS=0.00062223; nGRAD=0.00583775; meth=0; reach=2; PARS = (-0.365014, 4.11532, 0, -9.00423)
iter=170; LL=-1524.81; diff.LL=7.11486e-07; diff.PARS=0.000722478; nGRAD=0.00151391; meth=0; reach=7; PARS = (-0.364938, 4.11509, 0, -9.01232)
iter=175; LL=-1524.81; diff.LL=6.73501e-08; diff.PARS=6.83239e-05; nGRAD=0.00125183; meth=0; reach=12; PARS = (-0.364925, 4.11504, 0, -9.02043)
Optimization criteria met
####
Num Iter = 178
Params = (-0.364926, 4.11505, 0, -9.02546)
LL = -1524.81
GRAD = (-0.000517457, -8.52197e-05, 0, -0.000980572)
Convergence Indicators: 
   NormGrad = 0.001112
   NormIHessGrad = 0.00134964
Var = (0.00679209, 0.0553299, 0, 1.37823)

Sat Mar 22 07:18:49 2025: Get covariance ...

02.6833594137890100.419649743100121-1567.33915264465-0.3649260520085194.11504644276680-9.02546235522217-1524.8110.0011120004635448FALSE

c("alpha1", "lambda1", "kappa1", "theta")c(0.69424798722409, 61.2550590511141, 1, 0.000120307167658827)c(0.0572158610097444, 14.4086059003788, 0, 0.000141238508120364)c(0.582106960300542, 33.0147104189403, 1, -0.000156515221487254)c(0.806389014147639, 89.4954076832878, 1, 0.000397129556804908)c(0.590695379361215, 38.6295530198729, 1, 1.20502427627685e-05)c(0.815954017256611, 97.1324275335304, 1, 0.00120112224085714)

####
Sat Mar 22 07:18:49 2025: Try copula = Clayton, dist = Exp-Weibull, ...

Sat Mar 22 07:18:49 2025: Assume dependence: ...

Sat Mar 22 07:18:49 2025: theta = MLE ...

Sat Mar 22 07:18:49 2025: Grid search for initial parameters ...

Sat Mar 22 07:18:49 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

Sat Mar 22 07:18:57 2025: End GRID

Sat Mar 22 07:18:57 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)

Sat Mar 22 07:18:58 2025: No converged solutions! ...

####
Sat Mar 22 07:18:58 2025: Try copula = Gumbel, dist = Weibull, ...

Sat Mar 22 07:18:58 2025: Assume dependence: ...

Sat Mar 22 07:18:58 2025: theta = MLE ...

Sat Mar 22 07:18:58 2025: Grid search for initial parameters ...

Sat Mar 22 07:18:58 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

Sat Mar 22 07:18:58 2025: End GRID

Sat Mar 22 07:18:58 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.41; diff.LL=0.444019; diff.PARS=1; nGRAD=0.375437; meth=0; reach=0; PARS = (-0.357865, 4.08924, 0, -4.40813)
iter=10; LL=-1524.84; diff.LL=0.00141428; diff.PARS=0.0558812; nGRAD=0.0350472; meth=0; reach=0; PARS = (-0.364474, 4.11273, 0, -6.24722)
iter=15; LL=-1524.82; diff.LL=0.000708864; diff.PARS=0.116743; nGRAD=0.0288871; meth=0; reach=0; PARS = (-0.364709, 4.11445, 0, -7.70475)
iter=20; LL=-1524.82; diff.LL=8.62891e-06; diff.PARS=0.00157192; nGRAD=0.0157446; meth=0; reach=2; PARS = (-0.364788, 4.11462, 0, -7.74245)
iter=25; LL=-1524.81; diff.LL=7.97181e-05; diff.PARS=0.01777; nGRAD=0.00655421; meth=0; reach=0; PARS = (-0.364869, 4.11484, 0, -7.95171)
iter=30; LL=-1524.81; diff.LL=0.00103469; diff.PARS=0.379006; nGRAD=0.0259302; meth=0; reach=0; PARS = (-0.364473, 4.11352, 0, -8.61682)
iter=35; LL=-1524.81; diff.LL=2.54646e-05; diff.PARS=0.0119697; nGRAD=0.00366808; meth=0; reach=0; PARS = (-0.364941, 4.11513, 0, -8.7115)
iter=40; LL=-1524.81; diff.LL=2.11058e-05; diff.PARS=0.0102827; nGRAD=0.00270159; meth=0; reach=0; PARS = (-0.364906, 4.115, 0, -8.73225)
iter=45; LL=-1524.81; diff.LL=2.99556e-05; diff.PARS=0.0154; nGRAD=0.00213657; meth=0; reach=0; PARS = (-0.364916, 4.11503, 0, -8.78563)
iter=50; LL=-1524.81; diff.LL=3.29644e-06; diff.PARS=0.00176719; nGRAD=0.00257685; meth=0; reach=2; PARS = (-0.36493, 4.11509, 0, -8.82149)
iter=55; LL=-1524.81; diff.LL=3.10109e-05; diff.PARS=0.0169431; nGRAD=0.00268026; meth=0; reach=0; PARS = (-0.36493, 4.11509, 0, -8.8465)
iter=60; LL=-1524.81; diff.LL=1.86735e-05; diff.PARS=0.0129221; nGRAD=0.00269816; meth=0; reach=0; PARS = (-0.36494, 4.11514, 0, -9.1464)
iter=65; LL=-1524.81; diff.LL=1.11309e-05; diff.PARS=0.00832736; nGRAD=0.00142632; meth=0; reach=0; PARS = (-0.364923, 4.11506, 0, -9.20601)
iter=70; LL=-1524.81; diff.LL=9.88828e-06; diff.PARS=0.00798686; nGRAD=0.00126887; meth=0; reach=0; PARS = (-0.36493, 4.11509, 0, -9.23933)
iter=75; LL=-1524.81; diff.LL=3.00459e-05; diff.PARS=0.0259241; nGRAD=0.00240972; meth=0; reach=0; PARS = (-0.364907, 4.11503, 0, -9.30926)
iter=80; LL=-1524.81; diff.LL=2.00248e-06; diff.PARS=0.00180796; nGRAD=0.00208774; meth=0; reach=1; PARS = (-0.364932, 4.11511, 0, -9.33943)
iter=85; LL=-1524.81; diff.LL=3.54459e-06; diff.PARS=0.00329927; nGRAD=0.00137856; meth=0; reach=2; PARS = (-0.364929, 4.1151, 0, -9.37302)
iter=90; LL=-1524.81; diff.LL=2.61673e-07; diff.PARS=0.000248396; nGRAD=0.00128663; meth=0; reach=4; PARS = (-0.364937, 4.11513, 0, -9.39591)
iter=95; LL=-1524.81; diff.LL=1.95186e-05; diff.PARS=0.0192507; nGRAD=0.00246229; meth=0; reach=0; PARS = (-0.364969, 4.11523, 0, -9.43039)
iter=100; LL=-1524.81; diff.LL=1.13904e-05; diff.PARS=0.0118894; nGRAD=0.00202989; meth=0; reach=0; PARS = (-0.364907, 4.11503, 0, -9.48858)
iter=105; LL=-1524.81; diff.LL=9.55886e-06; diff.PARS=0.0103023; nGRAD=0.00181012; meth=0; reach=0; PARS = (-0.364942, 4.11515, 0, -9.52129)
iter=110; LL=-1524.81; diff.LL=1.09776e-06; diff.PARS=0.00119046; nGRAD=0.00121236; meth=0; reach=4; PARS = (-0.364945, 4.11515, 0, -9.53792)
iter=115; LL=-1524.81; diff.LL=1.19919e-06; diff.PARS=0.00135018; nGRAD=0.00176916; meth=0; reach=3; PARS = (-0.364957, 4.1152, 0, -9.5655)
iter=120; LL=-1524.81; diff.LL=3.05603e-06; diff.PARS=0.00346618; nGRAD=0.000947819; meth=0; reach=8; PARS = (-0.364936, 4.11512, 0, -9.57133)
iter=125; LL=-1524.81; diff.LL=1.70661e-06; diff.PARS=0.00197633; nGRAD=0.000978901; meth=0; reach=1; PARS = (-0.364938, 4.11513, 0, -9.59249)
iter=130; LL=-1524.81; diff.LL=1.16761e-05; diff.PARS=0.0146799; nGRAD=0.000831749; meth=0; reach=0; PARS = (-0.364938, 4.11513, 0, -9.67963)
iter=135; LL=-1524.81; diff.LL=1.41904e-06; diff.PARS=0.00180022; nGRAD=0.00123578; meth=0; reach=5; PARS = (-0.36494, 4.11515, 0, -9.68919)
iter=140; LL=-1524.81; diff.LL=3.38537e-09; diff.PARS=4.33578e-06; nGRAD=0.00175262; meth=0; reach=1; PARS = (-0.36495, 4.11518, 0, -9.70446)
iter=145; LL=-1524.81; diff.LL=8.48947e-06; diff.PARS=0.0117674; nGRAD=0.00174313; meth=0; reach=0; PARS = (-0.364961, 4.11522, 0, -9.78748)
iter=150; LL=-1524.81; diff.LL=4.73052e-06; diff.PARS=0.0068919; nGRAD=0.000917812; meth=0; reach=0; PARS = (-0.364943, 4.11515, 0, -9.82226)
iter=155; LL=-1524.81; diff.LL=4.56234e-05; diff.PARS=0.0783402; nGRAD=0.0116916; meth=0; reach=0; PARS = (-0.364804, 4.11461, 0, -9.97164)
iter=160; LL=-1524.81; diff.LL=5.17397e-06; diff.PARS=0.00904655; nGRAD=0.00163183; meth=0; reach=0; PARS = (-0.36494, 4.11516, 0, -10.0061)
iter=165; LL=-1524.81; diff.LL=4.34918e-06; diff.PARS=0.0077611; nGRAD=0.000939268; meth=0; reach=0; PARS = (-0.364931, 4.11511, 0, -10.0274)
iter=170; LL=-1524.81; diff.LL=3.39295e-06; diff.PARS=0.00606478; nGRAD=0.000590271; meth=0; reach=0; PARS = (-0.364937, 4.11513, 0, -10.0507)
iter=175; LL=-1524.81; diff.LL=4.00366e-06; diff.PARS=0.0074547; nGRAD=0.000863739; meth=0; reach=0; PARS = (-0.364951, 4.11518, 0, -10.0656)
iter=180; LL=-1524.81; diff.LL=3.36712e-07; diff.PARS=0.000652645; nGRAD=0.000814971; meth=0; reach=1; PARS = (-0.364941, 4.11515, 0, -10.1008)
iter=185; LL=-1524.81; diff.LL=1.70003e-08; diff.PARS=3.50468e-05; nGRAD=0.00055913; meth=0; reach=1; PARS = (-0.364939, 4.11514, 0, -10.1672)
iter=190; LL=-1524.81; diff.LL=8.4592e-07; diff.PARS=0.00177835; nGRAD=0.000675394; meth=0; reach=2; PARS = (-0.364946, 4.11517, 0, -10.181)
iter=195; LL=-1524.81; diff.LL=7.00924e-07; diff.PARS=0.00148166; nGRAD=0.00072548; meth=0; reach=7; PARS = (-0.364947, 4.11517, 0, -10.185)
iter=200; LL=-1524.81; diff.LL=9.0321e-08; diff.PARS=0.000191298; nGRAD=0.000855151; meth=0; reach=12; PARS = (-0.364945, 4.11517, 0, -10.1919)
####
Num Iter = 201
Params = (-0.364945, 4.11517, 0, -10.1919)
LL = -1524.81
GRAD = (-0.000545515, -0.000458522, 0, -0.00047271)
Convergence Indicators: 
   NormGrad = 0.000855151
   NormIHessGrad = 0.00746616
Var = (0.00683391, 0.0559384, 0, 15.8525)

Sat Mar 22 07:19:03 2025: Get covariance ...

02.6833594137890100-1565.92618810742-0.36494488243964.115166390088510-10.191920509471-1524.810.000855151069362377FALSE

c("alpha1", "lambda1", "kappa1", "theta")c(0.694234914358298, 61.2624068720554, 1, 1.0000374718524)c(0.0573906572640301, 14.4893503561323, 0, 0.000149195039020658)c(0.581751293071717, 32.8638020146535, 1, 0.999745054949248)c(0.806718535644879, 89.6610117294574, 1, 1.00032988875555)c(0.590391038853449, 38.5366446286816, 1, 1.00000001529772)c(0.816343888366012, 97.3899656267417, 1, 1.0917875151197)

####
Sat Mar 22 07:19:03 2025: Try copula = Gumbel, dist = Exp-Weibull, ...

Sat Mar 22 07:19:03 2025: Assume dependence: ...

Sat Mar 22 07:19:03 2025: theta = MLE ...

Sat Mar 22 07:19:03 2025: Grid search for initial parameters ...

Sat Mar 22 07:19: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

Sat Mar 22 07:19:13 2025: End GRID

Sat Mar 22 07:19:13 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)

Sat Mar 22 07:19:13 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.6942310 61.26397      1 0.0000000000 -1524.810  2
2   2     Clayton weibull 0.6942480 61.25506      1 0.0001203072 -1524.811  3
3   3      Gumbel weibull 0.6942349 61.26241      1 1.0000374719 -1524.810  3
        BIC       POST
1 -3063.162 0.99771395
2 -3069.934 0.00114188
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.3649505 0.08286297 -0.5273589 -0.2025421
2 log_lambda1  4.1151918 0.23799398  3.6487322  4.5816515
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.694231  0.05752604  0.581482  0.80698  0.5901616  0.8166521
2 lambda1 61.263966 14.58045546 32.686799 89.84113 38.4259192 97.6755708
3  kappa1  1.000000  0.00000000  1.000000  1.00000  1.0000000  1.0000000
4   theta  0.000000  0.00000000  0.000000  0.00000  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.182343e-06 -2.273737e-07  0.000000e+00  0.000000e+00 

# Hessian
res[[idx]]$RES$HESS

            log_alpha1 log_lambda1 unc_kappa1 log_theta
log_alpha1   -390.5097   -107.6660          0         0
log_lambda1  -107.6660    -47.3392          0         0
unc_kappa1      0.0000      0.0000          0         0
log_theta       0.0000      0.0000          0         0

# Covariance matrix
res[[idx]]$RES$COVAR

              log_alpha1 log_lambda1 unc_kappa1 log_theta
log_alpha1   0.006866272 -0.01561632          0         0
log_lambda1 -0.015616316  0.05664114          0         0
unc_kappa1   0.000000000  0.00000000          0         0
log_theta    0.000000000  0.00000000          0         0

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.03488153  -0.0355 0.9651227 0.006074835 0.9532161 0.9770294
3  1.0 -2.885480 0.05582800  -0.0574 0.9442163 0.007933657 0.9286663 0.9597663
         AA log_low_surv log_high_surv low_surv2 high_surv2
1 0.0000000   0.00000000    0.00000000 1.0000000  1.0000000
2 0.3474782  -0.05025002   -0.02507959 0.9509916  0.9752323
3 0.2866894  -0.07645740   -0.04309276 0.9263924  0.9578225

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

Session information

R version 4.4.3 (2025-02-28)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

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

other attached packages:
[1] dMrs_1.0.0        data.table_1.17.0 ggplot2_3.5.1     relsurv_2.3-2    
[5] survival_3.8-3    sqldf_0.4-11      RSQLite_2.3.9     gsubfn_0.7       
[9] proto_1.0.0      

loaded via a namespace (and not attached):
 [1] gtable_0.3.6        pspline_1.0-21      xfun_0.51          
 [4] bslib_0.9.0         caTools_1.18.3      lattice_0.22-6     
 [7] numDeriv_2016.8-1.1 vctrs_0.6.5         tools_4.4.3        
[10] bitops_1.0-9        generics_0.1.3      stats4_4.4.3       
[13] tibble_3.2.1        chron_2.3-62        blob_1.2.4         
[16] pkgconfig_2.0.3     Matrix_1.7-3        KernSmooth_2.23-26 
[19] lifecycle_1.0.4     farver_2.1.2        compiler_4.4.3     
[22] gplots_3.2.0        munsell_0.5.1       htmltools_0.5.8.1  
[25] sys_3.4.3           buildtools_1.0.0    sass_0.4.9         
[28] yaml_2.3.10         gmp_0.7-5           pillar_1.10.1      
[31] jquerylib_0.1.4     cachem_1.1.0        viridis_0.6.5      
[34] gtools_3.9.5        tidyselect_1.2.1    digest_0.6.37      
[37] mvtnorm_1.3-3       dplyr_1.1.4         labeling_0.4.3     
[40] maketools_1.3.2     splines_4.4.3       pcaPP_2.0-5        
[43] gsl_2.1-8           fastmap_1.2.0       grid_4.4.3         
[46] copula_1.1-6        colorspace_2.1-1    cli_3.6.4          
[49] magrittr_2.0.3      withr_3.0.2         Rmpfr_1.0-0        
[52] scales_1.3.0        bit64_4.6.0-1       rmarkdown_2.29     
[55] bit_4.6.0           gridExtra_2.3       memoise_2.0.1      
[58] evaluate_1.0.3      knitr_1.50          tcltk_4.4.3        
[61] ADGofTest_0.3       stabledist_0.7-2    viridisLite_0.4.2  
[64] rlang_1.1.5         Rcpp_1.0.14         glue_1.8.0         
[67] DBI_1.2.3           jsonlite_1.9.1      R6_2.6.1           

References