Import data:
load("./6Clusters.RData")Remove Patients without Extra.Mets:
data.analysis <- data.with.clusts %>% filter(!is.na(Extra.Mets))
data.analysis.id <- data.analysis[!duplicated(data.analysis$Patient),]data.no.missing <- data.analysis %>% filter(is.missing == FALSE)
data.no.missing.id <- data.analysis[!duplicated(data.analysis$Patient),]
model.surv.no.missing <- coxph(Surv(Survival.after.begining.ZA.mo, Death.Status) ~ Age.Diagnosis +
Sex + Extra.Mets, data = data.analysis.id, x = TRUE)# Obtain optimal tuning parameters using the *nlme* routine
rational.nlme <- nlme(model = log(NTX) ~ 1/(time + a)^g + i,
fixed = a + g + i ~ 1,
random = i ~ 1 | Patient,
start = c(a = 1, g = 1, i = 1),
data = data.no.missing)
# Compose the fixed and random effects formulas
fixed.rational.omit <- as.formula(sprintf("log(NTX) ~ I(1/(time + %f)^%f) + 1",
coef(rational.nlme)[1,"a"],
coef(rational.nlme)[1,"g"]))
random.rational.omit <- as.formula(sprintf("~ I(1/(time + %f)^%f) + 1 | Patient",
coef(rational.nlme)[1,"a"],
coef(rational.nlme)[1,"g"]))
rational.random.no.missing <- lme(fixed = fixed.rational.omit,
random = random.rational.omit,
method = "ML", data = data.no.missing,
control = lmeControl(opt = "optim"))
model.rational.no.missing <- jointModel(rational.random.no.missing, model.surv.no.missing,
timeVar = "time", method = "piecewise-PH-aGH",
lng.in.kn = 5)
summary(model.rational.no.missing)##
## Call:
## jointModel(lmeObject = rational.random.no.missing, survObject = model.surv.no.missing,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 524 Number of Events: 101 (80.2%)
## Number of Groups: 126
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with piecewise-constant
## baseline risk function
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -1244.717 2521.434 2566.814
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.1345 (Intr)
## I(1/(time + 1.223934)^1.14205) 0.2768 -0.9279
## Residual 0.8373
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.1157 0.1206 34.1287 <0.0001
## I(1/(time + 1.223934)^1.14205) 0.9575 0.1411 6.7874 <0.0001
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0185 0.0076 2.4280 0.0152
## Sex 0.3160 0.2195 1.4400 0.1499
## Extra.Mets 0.6655 0.2278 2.9217 0.0035
## Assoct 0.2093 0.1110 1.8857 0.0593
## log(xi.1) -6.7610 0.6963 -9.7099
## log(xi.2) -6.3122 0.6810 -9.2686
## log(xi.3) -6.0880 0.6763 -9.0020
## log(xi.4) -5.8937 0.6674 -8.8314
## log(xi.5) -5.9122 0.6390 -9.2527
## log(xi.6) -5.9400 0.5893 -10.0806
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0exponential.random.nlme <- nlme(log(NTX) ~ exp(-time*a) + int,
fixed = a + int ~ 1,
random = int ~ 1 | Patient,
data = data.no.missing,
start = c( a = 1, int = 3))
fixed.exponential.omit <- as.formula(sprintf("log(NTX) ~ I(exp(-time*%f)) + 1",
coef(exponential.random.nlme)$a[1]))
random.exponential.omit <- as.formula(sprintf("~ I(exp(-time*%f)) + 1 | Patient",
coef(exponential.random.nlme)$a[1]))
exponential.random.no.missing <- lme(fixed = fixed.exponential.omit,
random = random.exponential.omit,
method = "ML", data = data.no.missing,
control = lmeControl(opt = "optim"))
model.exponential.no.missing <- jointModel(exponential.random.no.missing, model.surv.no.missing,
timeVar = "time", method = "piecewise-PH-aGH",
lng.in.kn = 5)
summary(model.exponential.no.missing)##
## Call:
## jointModel(lmeObject = exponential.random.no.missing, survObject = model.surv.no.missing,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 524 Number of Events: 101 (80.2%)
## Number of Groups: 126
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with piecewise-constant
## baseline risk function
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -1246.003 2524.005 2569.386
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.0454 (Intr)
## I(exp(-time * 0.554024)) 0.0314 -0.8846
## Residual 0.8424
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.1536 0.1106 37.5429 <0.0001
## I(exp(-time * 0.554024)) 0.7015 0.0992 7.0725 <0.0001
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0184 0.0076 2.4078 0.0160
## Sex 0.3134 0.2193 1.4291 0.1530
## Extra.Mets 0.6649 0.2276 2.9214 0.0035
## Assoct 0.2175 0.1174 1.8536 0.0638
## log(xi.1) -6.7823 0.7121 -9.5246
## log(xi.2) -6.3328 0.6957 -9.1027
## log(xi.3) -6.1121 0.6919 -8.8343
## log(xi.4) -5.9202 0.6828 -8.6702
## log(xi.5) -5.9421 0.6544 -9.0803
## log(xi.6) -5.9782 0.6049 -9.8823
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0spline.random.no.missing <- lme(fixed = log(NTX) ~ ns(time,3) + 1,
random = ~ ns(time, 3) + 1 | Patient,
method = "ML", data = data.no.missing,
control = lmeControl(opt = "optim"))
# Needs more EM iterations to converge
model.spline.no.missing <- jointModel(spline.random.no.missing, model.surv.no.missing,
timeVar = "time", method = "piecewise-PH-aGH",
lng.in.kn = 5, iter.EM = 200)
summary(model.spline.no.missing)##
## Call:
## jointModel(lmeObject = spline.random.no.missing, survObject = model.surv.no.missing,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5,
## iter.EM = 200)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 524 Number of Events: 101 (80.2%)
## Number of Groups: 126
##
## Joint Model Summary:
## Longitudinal Process: Linear mixed-effects model
## Event Process: Relative risk model with piecewise-constant
## baseline risk function
## Parameterization: Time-dependent
##
## log.Lik AIC BIC
## -903.1675 1856.335 1927.242
##
## Variance Components:
## StdDev Corr
## (Intercept) 0.8701 (Intr) n(,3)1 n(,3)2
## ns(time, 3)1 0.5313 -0.7352
## ns(time, 3)2 0.4377 0.8764 -0.3203
## ns(time, 3)3 0.3792 -0.6283 0.8625 -0.2390
## Residual 0.6749
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.8562 0.0565 85.9593 <0.0001
## ns(time, 3)1 -0.3694 0.1408 -2.6235 0.0087
## ns(time, 3)2 -1.3193 0.1548 -8.5243 <0.0001
## ns(time, 3)3 -0.3441 0.1003 -3.4293 0.0006
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0184 0.0075 2.4402 0.0147
## Sex 0.2655 0.2173 1.2223 0.2216
## Extra.Mets 0.6051 0.2300 2.6302 0.0085
## Assoct 0.1095 0.0544 2.0119 0.0442
## log(xi.1) -6.2420 0.5859 -10.6528
## log(xi.2) -5.8026 0.5796 -10.0116
## log(xi.3) -5.5604 0.5773 -9.6310
## log(xi.4) -5.3427 0.5799 -9.2132
## log(xi.5) -5.2863 0.5659 -9.3411
## log(xi.6) -5.1188 0.5826 -8.7854
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 3
##
## Optimization:
## Convergence: 0