The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.