Import data:
load("./6Clusters.RData")data.analysis <- data.with.clusts
# Perform LOCF using the *na.locf* procedure from the **zoo** package
data.locf <- data.analysis %>% mutate(NTX = ifelse(is.missing, NA, NTX)) %>%
group_by(Patient) %>% mutate(NTX = zoo::na.locf(NTX, na.rm = FALSE)) %>% filter(!is.na(NTX)) %>% filter(!is.na(Extra.Mets))
data.locf.id <- data.locf[!duplicated(data.locf$Patient),]
library(survival)
model.surv.no.missing <- coxph(Surv(Survival.after.begining.ZA.mo, Death.Status) ~ Age.Diagnosis +
Sex + Extra.Mets, data = data.locf.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.locf)
fixed.rational.locf <- as.formula(sprintf("log(NTX) ~ I(1/(time + %f)^%f) + 1",
coef(rational.nlme)[1,"a"],
coef(rational.nlme)[1,"g"]))
random.rational.locf <- as.formula(sprintf("~ I(1/(time + %f)^%f) + 1 | Patient",
coef(rational.nlme)[1,"a"],
coef(rational.nlme)[1,"g"]))
rational.random.locf <- lme(fixed = fixed.rational.locf,
random = random.rational.locf,
method = "ML", data = data.locf,
control = lmeControl(opt = "optim"))
summary(rational.random.locf)## Linear mixed-effects model fit by maximum likelihood
## Data: data.locf
## AIC BIC logLik
## 1900.544 1927.816 -944.2718
##
## Random effects:
## Formula: ~I(1/(time + 1.370338)^1.468602) + 1 | Patient
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.2870805 (Intr)
## I(1/(time + 1.370338)^1.468602) 1.0385953 -0.669
## Residual 0.7139692
##
## Fixed effects: list(fixed.rational.locf)
## Value Std.Error DF t-value p-value
## (Intercept) 4.261362 0.1214442 569 35.08905 0
## I(1/(time + 1.370338)^1.468602) 0.978851 0.1572206 569 6.22597 0
## Correlation:
## (Intr)
## I(1/(time + 1.370338)^1.468602) -0.555
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -9.16561090 -0.42174586 0.05402667 0.37222123 4.73461615
##
## Number of Observations: 696
## Number of Groups: 126model.rational.locf <- jointModel(rational.random.locf, model.surv.no.missing,
timeVar = "time", method = "piecewise-PH-aGH",
lng.in.kn = 5)
summary(model.rational.locf)##
## Call:
## jointModel(lmeObject = rational.random.locf, survObject = model.surv.no.missing,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 696 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
## -1415.341 2862.681 2908.062
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.2891 (Intr)
## I(1/(time + 1.370338)^1.468602) 1.0385 -0.6712
## Residual 0.7140
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.2650 0.1215 35.1032 <0.0001
## I(1/(time + 1.370338)^1.468602) 0.9683 0.1572 6.1602 <0.0001
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0186 0.0076 2.4479 0.0144
## Sex 0.3006 0.2176 1.3811 0.1672
## Extra.Mets 0.6687 0.2266 2.9504 0.0032
## Assoct 0.1674 0.0941 1.7776 0.0755
## log(xi.1) -6.5896 0.6587 -10.0035
## log(xi.2) -6.1567 0.6499 -9.4735
## log(xi.3) -5.9373 0.6471 -9.1755
## log(xi.4) -5.7468 0.6400 -8.9798
## log(xi.5) -5.7666 0.6115 -9.4307
## log(xi.6) -5.8020 0.5625 -10.3151
##
## 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.locf,
start = c( a = 1, int = 3))
fixed.exponential.locf <- as.formula(sprintf("log(NTX) ~ I(exp(-time * %f)) + 1",
coef(exponential.random.nlme)$a[1]))
random.exponential.locf <- as.formula(sprintf("~ I(exp(-time * %f)) + 1 | Patient",
coef(exponential.random.nlme)$a[1]))
exponential.random.locf <- lme(fixed = fixed.exponential.locf,
random = random.exponential.locf,
method = "ML", data = data.locf,
control = lmeControl(opt = "optim"))
model.exponential.locf <- jointModel(exponential.random.locf, model.surv.no.missing,
timeVar = "time", method = "piecewise-PH-aGH",
lng.in.kn = 5)
summary(model.exponential.locf)##
## Call:
## jointModel(lmeObject = exponential.random.locf, survObject = model.surv.no.missing,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 696 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
## -1415.475 2862.949 2908.33
##
## Variance Components:
## StdDev Corr
## (Intercept) 1.2592 (Intr)
## I(exp(-time * 0.861081)) 0.6035 -0.6527
## Residual 0.7146
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.3055 0.1178 36.5377 <0.0001
## I(exp(-time * 0.861081)) 0.5693 0.0919 6.1976 <0.0001
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0186 0.0076 2.4406 0.0147
## Sex 0.3005 0.2175 1.3812 0.1672
## Extra.Mets 0.6673 0.2264 2.9474 0.0032
## Assoct 0.1692 0.0952 1.7784 0.0753
## log(xi.1) -6.5924 0.6616 -9.9651
## log(xi.2) -6.1635 0.6536 -9.4301
## log(xi.3) -5.9455 0.6511 -9.1313
## log(xi.4) -5.7555 0.6439 -8.9377
## log(xi.5) -5.7762 0.6155 -9.3851
## log(xi.6) -5.8129 0.5665 -10.2611
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 5
##
## Optimization:
## Convergence: 0spline.random.locf <- lme(fixed = log(NTX) ~ ns(time,3) + 1,
random = ~ ns(time, 3) + 1 | Patient,
method = "ML", data = data.locf,
control = lmeControl(opt = "optim"))
model.spline.locf <- jointModel(spline.random.locf, model.surv.no.missing,
timeVar = "time", method = "piecewise-PH-aGH",
lng.in.kn = 5)
summary(model.spline.locf)##
## Call:
## jointModel(lmeObject = spline.random.locf, survObject = model.surv.no.missing,
## timeVar = "time", method = "piecewise-PH-aGH", lng.in.kn = 5)
##
## Data Descriptives:
## Longitudinal Process Event Process
## Number of Observations: 696 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
## -1376.519 2803.037 2873.945
##
## Variance Components:
## StdDev Corr
## (Intercept) 0.9582 (Intr) n(,3)1 n(,3)2
## ns(time, 3)1 0.6542 -0.3155
## ns(time, 3)2 0.9078 0.5543 0.6110
## ns(time, 3)3 0.5545 -0.1383 0.9777 0.7357
## Residual 0.6508
##
## Coefficients:
## Longitudinal Process
## Value Std.Err z-value p-value
## (Intercept) 4.8851 0.0634 77.1063 <0.0001
## ns(time, 3)1 -0.2576 0.1191 -2.1629 0.0305
## ns(time, 3)2 -1.0792 0.1454 -7.4211 <0.0001
## ns(time, 3)3 -0.2157 0.0806 -2.6762 0.0074
##
## Event Process
## Value Std.Err z-value p-value
## Age.Diagnosis 0.0175 0.0077 2.2547 0.0242
## Sex 0.2930 0.2195 1.3347 0.1820
## Extra.Mets 0.6322 0.2265 2.7908 0.0053
## Assoct 0.1353 0.0619 2.1852 0.0289
## log(xi.1) -6.3419 0.5850 -10.8412
## log(xi.2) -5.9094 0.5798 -10.1919
## log(xi.3) -5.6709 0.5787 -9.7988
## log(xi.4) -5.4442 0.5802 -9.3836
## log(xi.5) -5.3798 0.5643 -9.5329
## log(xi.6) -5.1128 0.5761 -8.8747
##
## Integration:
## method: (pseudo) adaptive Gauss-Hermite
## quadrature points: 3
##
## Optimization:
## Convergence: 0