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),]
model.surv <- coxph(Surv(Survival.after.begining.ZA.mo, Death.Status) ~ Age.Diagnosis + Sex + Extra.Mets,
data = data.analysis.id, x = TRUE)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.analysis)
fixed.rational.ocs <- as.formula(sprintf("log(NTX) ~ I(1/(time + %f)^%f) + 1",
coef(rational.nlme)[1,"a"],
coef(rational.nlme)[1,"g"]))
random.rational.ocs <- as.formula(sprintf("~ I(1/(time + %f)^%f) + 1 | Patient",
coef(rational.nlme)[1,"a"],
coef(rational.nlme)[1,"g"]))
rational.random <- lme(fixed = fixed.rational.ocs,
random = random.rational.ocs,
method = "ML", data = data.analysis,
control = lmeControl(opt = "optim"))
model.rational.random <- jointModel(rational.random, model.surv, timeVar = "time",
method = "piecewise-PH-aGH", lng.in.kn = 5)
summary(model.rational.random)##
## Call:
## jointModel(lmeObject = rational.random, survObject = model.surv,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 698 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
## -1395.13 2822.26 2867.641
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.2227 (Intr)
## I(1/(time + 1.351057)^2.049934) 0.5752 -0.9750
## Residual 0.7125
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.3451 0.1145 37.9506 <0.0001
## I(1/(time + 1.351057)^2.049934) 0.9768 0.1493 6.5429 <0.0001
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0182 0.0076 2.3998 0.0164
## Sex 0.3048 0.2176 1.4011 0.1612
## Extra.Mets 0.6586 0.2264 2.9091 0.0036
## Assoct 0.1955 0.0966 2.0233 0.0430
## log(xi.1) -6.6908 0.6691 -10.0002
## log(xi.2) -6.2567 0.6596 -9.4860
## log(xi.3) -6.0353 0.6561 -9.1981
## log(xi.4) -5.8409 0.6488 -9.0025
## log(xi.5) -5.8599 0.6212 -9.4328
## log(xi.6) -5.8933 0.5729 -10.2876
##
## 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.analysis,
start = c( a = 1, int = 3))
fixed.exponential.ocs <- as.formula(sprintf("log(NTX) ~ I(exp(-time*%f)) + 1",
coef(exponential.random.nlme)$a[1]))
random.exponential.ocs <- as.formula(sprintf("~ I(exp(-time*%f)) + 1 | Patient",
coef(exponential.random.nlme)$a[1]))
exponential.random <- lme(fixed = fixed.exponential.ocs,
random = random.exponential.ocs,
method = "ML", data = data.analysis,
control = lmeControl(opt = "optim"))
model.exponential.random <- jointModel(exponential.random, model.surv, timeVar = "time",
method = "piecewise-PH-aGH", lng.in.kn = 5)
summary(model.exponential.random)##
## Call:
## jointModel(lmeObject = exponential.random, survObject = model.surv,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 698 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
## -1395.08 2822.16 2867.541
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.2154 (Intr)
## I(exp(-time * 1.092705)) 0.3075 -0.9495
## Residual 0.7118
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.3580 0.1136 38.3766 <0.0001
## I(exp(-time * 1.092705)) 0.5137 0.0783 6.5571 <0.0001
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0181 0.0076 2.3956 0.0166
## Sex 0.3046 0.2175 1.4001 0.1615
## Extra.Mets 0.6575 0.2263 2.9050 0.0037
## Assoct 0.1962 0.0970 2.0233 0.0430
## log(xi.1) -6.6874 0.6698 -9.9836
## log(xi.2) -6.2591 0.6608 -9.4713
## log(xi.3) -6.0383 0.6575 -9.1835
## log(xi.4) -5.8437 0.6501 -8.9884
## log(xi.5) -5.8630 0.6225 -9.4178
## log(xi.6) -5.8975 0.5742 -10.2709
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0spline.random <- lme(fixed = log(NTX) ~ ns(time, 3) + 1,
random = ~ ns(time, 3) + 1 | Patient,
method = "ML", data = data.analysis,
control = lmeControl(opt = "optim"))
# Needs more EM iterations to avoid convergence problems
model.spline.random <- jointModel(spline.random, model.surv, timeVar = "time",
method = "piecewise-PH-aGH",
lng.in.kn = 5, iter.EM = 500)
summary(model.spline.random)##
## Call:
## jointModel(lmeObject = spline.random, survObject = model.surv,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5,
## iter.EM = 500)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 698 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
## -1280.885 2611.771 2682.678
##
## Variance Components:
## StdDev Corr
## (Intercept) 0.9675 (Intr) n(,3)1 n(,3)2
## ns(time, 3)1 0.4994 -0.4523
## ns(time, 3)2 0.8039 0.5191 0.5185
## ns(time, 3)3 0.2122 -0.4303 0.9445 0.4443
## Residual 0.6231
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.8073 0.0528 91.0829 <0.0001
## ns(time, 3)1 -0.0960 0.1065 -0.9021 0.3670
## ns(time, 3)2 -1.0021 0.1352 -7.4126 <0.0001
## ns(time, 3)3 -0.1198 0.0694 -1.7255 0.0844
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0170 0.0075 2.2688 0.0233
## Sex 0.3153 0.2188 1.4408 0.1496
## Extra.Mets 0.6681 0.2264 2.9509 0.0032
## Assoct 0.1966 0.0611 3.2200 0.0013
## log(xi.1) -6.6290 0.6395 -10.3661
## log(xi.2) -6.1809 0.6308 -9.7982
## log(xi.3) -5.9290 0.6295 -9.4193
## log(xi.4) -5.6756 0.6298 -9.0120
## log(xi.5) -5.5331 0.5997 -9.2264
## log(xi.6) -5.0740 0.5996 -8.4625
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 3
##
## Optimization:
## Convergence: 0