Getting familiar with dMrs
Paul Little, Reuben Adatorwovor
Last Updated: 01/19/2025
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: gsubfnLoading required package: protoLoading required package: RSQLiteLoading required package: survivalLoading required package: date2 Application
The code below will load relsurv’s working dataset
rdata and import Slovenia’s latest ratetable from HMD.
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"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[1] 1040 6rdata$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 [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"[1] 872 73 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 1Prep 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.87rDAT = 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.00051Perform matching to calculate log density and log cdf.
.872 out of 872 done 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 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.8746426Prep 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.9957323Run 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 doneSun Jan 19 01:01:52 2025: End GRIDSun Jan 19 01:01:52 2025: NR optimization w/ 1 profile point(s) ...#---# Copula=Independent, Dist=weibull, Prof. point=1 out of 1iPARS = (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.000117290238296447FALSEc("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 doneSun 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 1iPARS = (-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 doneSun 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 1iPARS = (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.00190746698854164FALSEc("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 doneSun 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 2iPARS = (-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 2iPARS = (-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 doneSun 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 1iPARS = (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.00113166168428236FALSEc("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 doneSun 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 2iPARS = (-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 2iPARS = (-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
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[1] 1 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 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[1] -1524.81 log_alpha1 log_lambda1 unc_kappa1 log_theta
-1.043190e-04 -5.361471e-05 0.000000e+00 0.000000e+00 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 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 05 Net-survival
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
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