jointModel <- function (lmeObject, survObject, timeVar, parameterization = c("value", "slope", "both"),
method = c("weibull-PH-aGH", "weibull-PH-GH", "weibull-AFT-aGH", "weibull-AFT-GH",
"piecewise-PH-aGH", "piecewise-PH-GH", "Cox-PH-aGH", "Cox-PH-GH", "spline-PH-aGH",
"spline-PH-GH", "ch-Laplace"), interFact = NULL, derivForm = NULL, lag = 0,
scaleWB = NULL, CompRisk = FALSE, init = NULL, control = list(), ...) {
cl <- match.call()
if (!inherits(lmeObject, "lme"))
stop("\n'lmeObject' must inherit from class lme.")
if (length(lmeObject$group) > 1)
stop("\nnested random-effects are not allowed in lme().")
if (!is.null(lmeObject$modelStruct$corStruct))
warning("correlation structure in 'lmeObject' is ignored.\n")
if (!is.null(lmeObject$modelStruct$varStruct))
warning("variance structure in 'lmeObject' is ignored.\n")
if (!inherits(survObject, "coxph") && !inherits(survObject, "survreg"))
stop("\n'survObject' must inherit from class coxph or class survreg.")
if (!is.matrix(survObject$x))
stop("\nuse argument 'x = TRUE' in ",
if (inherits(survObject, "coxph")) "'coxph()'." else "'survreg()'.")
if (length(timeVar) != 1 || !is.character(timeVar))
stop("\n'timeVar' must be a character string.")
method. <- match.arg(method)
method <- switch(method., "weibull-AFT-GH" = , "weibull-AFT-aGH" = "weibull-AFT-GH",
"weibull-PH-GH" = , "weibull-PH-aGH" = "weibull-PH-GH",
"piecewise-PH-GH" = , "piecewise-PH-aGH" = "piecewise-PH-GH",
"Cox-PH-GH" = , "Cox-PH-aGH" = "Cox-PH-GH", "spline-PH-GH" =,
"spline-PH-aGH" = "spline-PH-GH", "ch-Laplace" = "ch-Laplace")
parameterization <- match.arg(parameterization)
if (method == "Cox-PH-GH" && !inherits(survObject, "coxph"))
stop("\nfor 'method = Cox-PH-GH', 'survObject' must inherit from class coxph.")
if (parameterization %in% c("slope", "both") && method %in% c("Cox-PH-GH", "ch-Laplace"))
stop("\nthe slope parameterization is not currently available for methods 'Cox-PH-GH' & 'ch-Laplace'.")
if (parameterization %in% c("slope", "both") && is.null(derivForm)) {
stop("\nwhen parameterization is 'slope' or 'both' you need to specify the 'derivForm' argument.")
}
if (parameterization %in% c("slope", "both") && !is.list(derivForm)) {
stop("\nthe 'derivForm' argument must be a list with components 'fixed' (a formula),\n\t'indFixed'",
"(a numeric vector), 'random' (a formula) and 'indRandom' (a numeric vector).")
}
if (!is.null(interFact) && !is.list(interFact)) {
stop("\nthe 'interFact' argument must be a list -- check the help file for more info.")
}
if (!is.null(interFact) && method %in% c("Cox-PH-GH", "ch-Laplace")) {
stop("\nincluding interaction terms is not currently available for methods 'Cox-PH-GH' & 'ch-Laplace'.")
}
if (CompRisk && (method != "spline-PH-GH" || is.null(survObject$strata))) {
stop("\nto fit a competing risks joint model you must choose as method 'spline-PH-GH'",
" and include a strata() in the specification of the coxph().")
}
# survival process
formT <- formula(survObject)
if (inherits(survObject, "coxph")) {
W <- survObject$x
keepW <- suppressWarnings(!is.na(survObject$coefficients))
W <- W[, keepW, drop = FALSE]
if (CompRisk) {
nRisks <- length(unique(survObject$strata))
} else {
nRisks <- 1
}
surv <- survObject$y
if (attr(surv, "type") == "right") {
LongFormat <- FALSE
Time <- survObject$y[, 1]
d <- survObject$y[, 2]
} else if (attr(surv, "type") == "counting") {
LongFormat <- TRUE
if (is.null(survObject$model))
stop("\nplease refit the Cox model including in the ",
"call to coxph() the argument 'model = TRUE'.")
Time <- survObject$y[, 2]
d <- survObject$y[, 3]
}
idT <- if (!is.null(survObject$model$`(cluster)`)) {
as.vector(unclass(survObject$model$`(cluster)`))
} else {
if (!CompRisk) seq_along(Time)
else rep(seq_len(length(Time)/nRisks), each = nRisks)
}
idT <- match(idT, unique(idT))
} else {
W <- survObject$x[, -1, drop = FALSE]
Time <- exp(survObject$y[, 1])
d <- survObject$y[, 2]
idT <- seq_along(Time)
LongFormat <- FALSE
nRisks <- 1
}
nT <- length(unique(idT))
if (LongFormat && is.null(survObject$model$`(cluster)`))
stop("\nuse argument 'model = TRUE' and cluster() in coxph().")
if (!length(W))
W <- NULL
if (sum(d) < 5)
warning("\nmore than 5 events are required.")
WintF.vl <- WintF.sl <- as.matrix(rep(1, length(Time)))
if (!is.null(interFact)) {
if (!is.null(interFact$value))
WintF.vl <- if (is.null(survObject$model)) {
model.matrix(interFact$value, data = interFact$data)
} else {
model.matrix(interFact$value, data = survObject$model)
}
if (!is.null(interFact$slope))
WintF.sl <- if (is.null(survObject$model)) {
model.matrix(interFact$slope, data = interFact$data)
} else {
model.matrix(interFact$slope, data = survObject$model)
}
}
# longitudinal process
id <- lmeObject$groups[[1]]
id <- match(id, unique(id))
b <- data.matrix(ranef(lmeObject))
dimnames(b) <- NULL
nY <- nrow(b)
if (nY != nT)
stop("sample sizes in the longitudinal and event processes differ; ",
"maybe you forgot the cluster() argument.\n")
TermsX <- lmeObject$terms
data <- lmeObject$data[all.vars(TermsX)]
data <- data[complete.cases(data), ]
formYx <- formula(lmeObject)
mfX <- model.frame(TermsX, data = data)
X <- model.matrix(formYx, mfX)
formYz <- formula(lmeObject$modelStruct$reStruct[[1]])
mfZ <- model.frame(terms(formYz), data = data)
TermsZ <- attr(mfZ, "terms")
Z <- model.matrix(formYz, mfZ)
y.long <- model.response(mfX, "numeric")
data.id <- data[!duplicated(id), ]
data.id <- data.id[idT, ]
if (!timeVar %in% names(data))
stop("\n'timeVar' does not correspond to one of the columns in the model.frame of 'lmeObject'.")
# check if there are any longitudinal measurements after the event times
max.timeY <- tapply(data[[timeVar]], factor(id, unique(id)), max)
max.timeT <- tapply(Time, factor(idT, unique(idT)), max)
if (!isTRUE(all(max.timeT >= max.timeY))) {
idnams <- factor(lmeObject$groups[[1]])
#stop("\nit seems that there are longitudinal measurements taken after the event times for some subjects ",
# "(i.e., check subject(s): ", paste(levels(idnams)[(max.timeT < max.timeY)], collapse = ", "), ").")
}
# extra design matrices for the longitudinal part
data.id[[timeVar]] <- pmax(Time - lag, 0)
if (parameterization %in% c("value", "both")) {
mfX.id <- model.frame(TermsX, data = data.id)
mfZ.id <- model.frame(TermsZ, data = data.id)
Xtime <- model.matrix(formYx, mfX.id)
Ztime <- model.matrix(formYz, mfZ.id)
long <- c(X %*% fixef(lmeObject)) + rowSums(Z * b[id, ])
}
if (parameterization %in% c("slope", "both")) {
mfX.deriv <- model.frame(terms(derivForm$fixed), data = data)
TermsX.deriv <- attr(mfX.deriv, "terms")
mfZ.deriv <- model.frame(terms(derivForm$random), data = data)
TermsZ.deriv <- attr(mfZ.deriv, "terms")
mfX.deriv.id <- model.frame(TermsX.deriv, data = data.id)
mfZ.deriv.id <- model.frame(TermsZ.deriv, data = data.id)
Xtime.deriv <- model.matrix(derivForm$fixed, mfX.deriv.id)
Ztime.deriv <- model.matrix(derivForm$random, mfZ.deriv.id)
Xderiv <- model.matrix(derivForm$fixed, mfX.deriv)
Zderiv <- model.matrix(derivForm$random, mfZ.deriv)
long.deriv <- as.vector(c(Xderiv %*% fixef(lmeObject)[derivForm$indFixed]) +
if (length(derivForm$indRandom) > 1 || derivForm$indRandom)
rowSums(Zderiv * b[id, derivForm$indRandom, drop = FALSE])
else
rep(0, nrow(Zderiv)))
}
if (parameterization == "value")
long.deriv <- NULL
if (parameterization == "slope")
long <- NULL
# response vectors and design matrices
y <- list(y = y.long, logT = log(Time), d = d, lag = lag)
x <- list(X = X, Z = Z, W = W, WintF.vl = WintF.vl,
WintF.sl = WintF.sl, idT = idT, nRisks = nRisks)
x <- switch(parameterization,
"value" = c(x, list(Xtime = Xtime, Ztime = Ztime)),
"slope" = c(x, list(Xtime.deriv = Xtime.deriv, Ztime.deriv = Ztime.deriv)),
"both" = c(x, list(Xtime = Xtime, Ztime = Ztime, Xtime.deriv = Xtime.deriv,
Ztime.deriv = Ztime.deriv)))
# control values
ind.noadapt <- method. %in% c("weibull-AFT-GH", "weibull-PH-GH", "piecewise-PH-GH",
"Cox-PH-GH", "spline-PH-GH")
con <- list(only.EM = FALSE, iter.EM = if (method == "spline-PH-GH") 120 else 50,
iter.qN = 350, optimizer = "optim", tol1 = 1e-03, tol2 = 1e-04,
tol3 = if (!CompRisk) sqrt(.Machine$double.eps) else 1e-09, numeriDeriv = "fd", eps.Hes = 1e-06,
parscale = NULL, step.max = 0.1, backtrackSteps = 2,
knots = NULL, ObsTimes.knots = TRUE,
lng.in.kn = if (method == "piecewise-PH-GH") 6 else 5, ord = 4,
equal.strata.knots = TRUE, typeGH = if (ind.noadapt) "simple" else "adaptive",
GHk = if (ncol(Z) < 3 && nrow(Z) < 2000) 15 else 9,
GKk = if (method == "piecewise-PH-GH" || length(Time) > nRisks*nT) 7 else 15, verbose = FALSE)
if (method == "Cox-PH-GH") {
con$only.EM <- TRUE
con$iter.EM <- 200
con$GHk <- if (ncol(Z) == 1) 15 else if (ncol(Z) == 2) 11 else 9
}
control <- c(control, list(...))
namC <- names(con)
con[(namc <- names(control))] <- control
if (con$typeGH != "simple" && !"GHk" %in% namc) {
con$GHk <- if (ncol(Z) <= 3 && nrow(Z) < 2000) 5 else 3
}
if (length(noNms <- namc[!namc %in% namC]) > 0)
warning("unknown names in 'control': ", paste(noNms, collapse = ", "))
if (method == "Cox-PH-GH" && !con$only.EM)
stop("with method 'Cox-PH-GH' only the EM algorithm is used.\n")
if (method == "Cox-PH-GH" && any(!is.na(match(c("iter.qN", "optimizer"), namc))))
warning("method 'Cox-PH-GH' uses only the EM algorithm.\n")
# extra design matrices for 'method = "weibull-AFT-GH"' and 'method = "weibull-PH-GH"'
# extra design matrices for 'method = "spline-PH-GH"' and 'method = "spline-PH-Laplace"'
if (method %in% c("weibull-AFT-GH", "weibull-PH-GH", "spline-PH-GH", "spline-PH-Laplace")) {
wk <- gaussKronrod(con$GKk)$wk
sk <- gaussKronrod(con$GKk)$sk
if (LongFormat) {
Time0 <- survObject$y[, 1]
P <- (Time - Time0) / 2
P1 <- (Time + Time0) / 2
st <- outer(P, sk) + P1
} else {
P <- as.vector(Time)/2
st <- outer(P, sk + 1)
}
dimnames(st) <- names(P) <- NULL
id.GK <- rep(seq_along(Time), each = con$GKk)
data.id2 <- data.id[id.GK, , drop = FALSE]
data.id2[[timeVar]] <- pmax(c(t(st)) - lag, 0)
if (parameterization %in% c("value", "both")) {
mfX <- model.frame(TermsX, data = data.id2)
mfZ <- model.frame(TermsZ, data = data.id2)
Xs <- model.matrix(formYx, mfX)
Zs <- model.matrix(formYz, mfZ)
}
if (parameterization %in% c("slope", "both")) {
mfX.deriv <- model.frame(TermsX.deriv, data = data.id2)
mfZ.deriv <- model.frame(TermsZ.deriv, data = data.id2)
Xs.deriv <- model.matrix(derivForm$fixed, mfX.deriv)
Zs.deriv <- model.matrix(derivForm$random, mfZ.deriv)
}
Ws.intF.vl <- WintF.vl[id.GK, , drop = FALSE]
Ws.intF.sl <- WintF.sl[id.GK, , drop = FALSE]
x <- c(x, list(P = P, st = c(t(st)), wk = wk, Ws.intF.vl = Ws.intF.vl,
Ws.intF.sl = Ws.intF.sl))
x <- switch(parameterization,
"value" = c(x, list(Xs = Xs, Zs = Zs)),
"slope" = c(x, list(Xs.deriv = Xs.deriv, Zs.deriv = Zs.deriv)),
"both" = c(x, list(Xs.deriv = Xs.deriv, Zs.deriv = Zs.deriv, Xs = Xs, Zs = Zs)))
if (method == "spline-PH-GH" || method == "spline-PH-Laplace") {
strt <- if (is.null(survObject$strata)) gl(1, length(Time)) else survObject$strata
nstrt <- length(levels(strt))
split.Time <- split(Time, strt)
ind.t <- if (LongFormat) {
unlist(tapply(idT, idT,
function (x) c(rep(FALSE, length(x) - 1), TRUE)))
} else {
rep(TRUE, length(Time))
}
kn <- if (con$equal.strata.knots) {
kk <- if (is.null(con$knots)) {
pp <- seq(0, 1, length.out = con$lng.in.kn + 2)
pp <- tail(head(pp, -1), -1)
quantile(Time[ind.t], pp, names = FALSE)
} else {
con$knots
}
kk <- kk[kk < max(Time)]
rr <- rep(list(sort(c(rep(range(Time, st), con$ord), kk))), nstrt)
names(rr) <- names(split.Time)
rr
} else {
spt <- if (length(Time) > nT & !CompRisk)
mapply(function (x, y){
x[unlist(tapply(y, y,
function (z) c(rep(FALSE, length(z) - 1), TRUE)))]
}, split.Time, split(idT, strt), SIMPLIFY = FALSE)
else split.Time
lapply(spt, function (t) {
kk <- if (is.null(con$knots)) {
pp <- seq(0, 1, length.out = con$lng.in.kn + 2)
pp <- tail(head(pp, -1), -1)
quantile(t, pp, names = FALSE)
} else {
con$knots
}
kk <- kk[kk < max(t)]
sort(c(rep(range(Time, st), con$ord), kk))
})
}
con$knots <- kn
W2 <- mapply(function (k, t) splineDesign(k, t, ord = con$ord), kn,
split.Time, SIMPLIFY = FALSE)
if (any(sapply(W2, colSums) == 0))
stop("\nsome of the knots of the B-splines basis are set outside the range",
"\n of the observed event times for one of the strata; refit the model",
"\n setting the control argument 'equal.strata.knots' to FALSE.")
W2 <- mapply(function (w2, ind) {
out <- matrix(0, length(Time), ncol(w2))
out[strt == ind, ] <- w2
out
}, W2, levels(strt), SIMPLIFY = FALSE)
W2 <- do.call(cbind, W2)
strt.s <- rep(strt, each = con$GKk)
split.Time <- split(c(t(st)), strt.s)
W2s <- mapply(function (k, t) splineDesign(k, t, ord = con$ord),
kn, split.Time, SIMPLIFY = FALSE)
W2s <- mapply(function (w2s, ind) {
out <- matrix(0, length(Time) * con$GKk, ncol(w2s))
out[strt.s == ind, ] <- w2s
out
}, W2s, levels(strt), SIMPLIFY = FALSE)
W2s <- do.call(cbind, W2s)
y <- c(y, list(strata = strt))
x <- c(x, list(W2 = W2, W2s = W2s))
}
}
# extra design matrices for 'method = "piecewise-PH-GH"'
if (method == "piecewise-PH-GH") {
wk <- gaussKronrod(con$GKk)$wk
sk <- gaussKronrod(con$GKk)$sk
nk <- length(sk)
if (is.null(con$knots) || !is.numeric(con$knots)) {
Q <- con$lng.in.kn + 1
qs <- if (con$ObsTimes.knots) {
unique(quantile(Time, seq(0, 1, len = Q + 1),
names = FALSE)[-c(1, Q + 1)])
} else {
unique(quantile(Time[d == 1], seq(0, 1, len = Q - 1),
names = FALSE))
}
qs <- qs + 1e-06
if (max(qs) > max(Time))
qs[which.max(qs)] <- max(Time) - 1e-06
con$knots <- qs
qs <- c(0, qs, max(Time) + 1)
Q <- length(qs) - 1
} else {
qs <- c(0, sort(con$knots) + 1e-06, max(Time) + 1)
Q <- length(qs) - 1
}
ind <- findInterval(Time, qs, rightmost.closed = TRUE)
D <- matrix(0, length(ind), Q)
D[cbind(seq_along(ind), ind)] <- 1
D <- D * d
Tiq <- outer(Time, qs, pmin)
Lo <- Tiq[, 1:Q]
Up <- Tiq[, 2:(Q+1)]
T <- Up - Lo
P <- T / 2
P[P < con$tol3] <- as.numeric(NA)
P1 <- (Up + Lo) / 2
st <- matrix(0, nY, nk*Q)
skQ <- rep(sk, Q)
for (i in seq_len(nY)) {
st[i, ] <- rep(P[i, ], each = nk) * skQ + rep(P1[i, ], each = nk)
}
y <- c(y, list(ind.D = ind))
id.GK <- rep(seq_len(nY), rowSums(!is.na(st)))
P <- c(t(P))
data.id2 <- data.id[rep(seq_len(nY), each = nk*Q), ]
data.id2[[timeVar]] <- pmax(c(t(st)) - lag, 0)
data.id2 <- data.id2[!is.na(data.id2[[timeVar]]), ]
if (parameterization %in% c("value", "both")) {
mfX <- model.frame(TermsX, data = data.id2)
mfZ <- model.frame(TermsZ, data = data.id2)
Xs <- model.matrix(formYx, mfX)
Zs <- model.matrix(formYz, mfZ)
}
if (parameterization %in% c("slope", "both")) {
mfX.deriv <- model.frame(TermsX.deriv, data = data.id2)
mfZ.deriv <- model.frame(TermsZ.deriv, data = data.id2)
Xs.deriv <- model.matrix(derivForm$fixed, mfX.deriv)
Zs.deriv <- model.matrix(derivForm$random, mfZ.deriv)
}
Ws.intF.vl <- WintF.vl[id.GK, , drop = FALSE]
Ws.intF.sl <- WintF.sl[id.GK, , drop = FALSE]
x <- c(x, list(P = P[!is.na(P)], st = st[!is.na(st)], wk = wk,
id.GK = id.GK, Q = Q, Ws.intF.vl = Ws.intF.vl, Ws.intF.sl = Ws.intF.sl))
x <- switch(parameterization,
"value" = c(x, list(Xs = Xs, Zs = Zs)),
"slope" = c(x, list(Xs.deriv = Xs.deriv, Zs.deriv = Zs.deriv)),
"both" = c(x, list(Xs.deriv = Xs.deriv,
Zs.deriv = Zs.deriv, Xs = Xs, Zs = Zs))
)
}
# extra design matrices for 'method = "Cox-PH-GH"' with event times prior to observed time for the ith subject
if (method == "Cox-PH-GH") {
unqT <- sort(unique(Time[d == 1]))
times <- lapply(Time, function (t) unqT[t >= unqT])
ind.len <- sapply(times, length)
indT <- rep(1:nrow(data.id), ind.len)
data.id2 <- data.id[indT, ]
data.id2[timeVar] <- pmax(unlist(times, use.names = FALSE) - lag, 0)
if (parameterization %in% c("value", "both")) {
mfX <- model.frame(TermsX, data = data.id2)
mfZ <- model.frame(TermsZ, data = data.id2)
Xtime2 <- model.matrix(formYx, mfX)
Ztime2 <- model.matrix(formYz, mfZ)
}
if (parameterization %in% c("slope", "both")) {
mfX.deriv <- model.frame(TermsX.deriv, data = data.id2)
mfZ.deriv <- model.frame(TermsZ.deriv, data = data.id2)
Xtime2.deriv <- model.matrix(derivForm$fixed, mfX.deriv)
Ztime2.deriv <- model.matrix(derivForm$random, mfZ.deriv)
}
x <- c(x, list(indT = indT))
x <- switch(parameterization,
"value" = c(x, list(Xtime2 = Xtime2, Ztime2 = Ztime2)),
"slope" = c(x, list(Xtime2.deriv = Xtime2.deriv,
Ztime2.deriv = Ztime2.deriv)),
"both" = c(x, list(Xtime2.deriv = Xtime2.deriv,
Ztime2.deriv = Ztime2.deriv, Xtime2 = Xtime2, Ztime2 = Ztime2))
)
}
# initial values
VC <- lapply(pdMatrix(lmeObject$modelStruct$reStruct), "*",
lmeObject$sigma^2)[[1]]
if (con$typeGH != "simple") {
Vs <- vector("list", nY)
inv.VC <- solve(VC)
for (i in 1:nY) {
Z.i <- Z[id == i, , drop = FALSE]
Vs[[i]] <- solve(crossprod(Z.i) / lmeObject$sigma^2 + inv.VC)
}
con$inv.chol.VCs <- lapply(Vs, function (x) solve(chol(solve(x))))
con$det.inv.chol.VCs <- sapply(con$inv.chol.VCs, det)
}
con$inv.chol.VC <- solve(chol(solve(VC)))
con$det.inv.chol.VC <- det(con$inv.chol.VC)
con$ranef <- b
if (all(VC[upper.tri(VC)] == 0))
VC <- diag(VC)
init.surv <- initial.surv(Time, d, W, WintF.vl, WintF.sl, id,
times = data[[timeVar]], method, parameterization, long = long,
long.deriv = long.deriv,
extra = list(W2 = x$W2, control = con, ii = idT, strata = survObject$strata),
LongFormat = CompRisk | length(Time) > nT)
if (method == "Cox-PH-GH" && length(init.surv$lambda0) < length(unqT))
init.surv$lambda0 <- basehaz(survObject)$hazard
initial.values <- c(list(betas = fixef(lmeObject), sigma = lmeObject$sigma, D = VC), init.surv)
if (!is.null(init)) {
nams1 <- names(init)
nams2 <- names(initial.values)
if (!is.list(init) || length(noNms <- nams1[!nams1 %in% nams2])) {
warning("unknown names in 'init': ", paste(noNms, collapse = ", "))
} else {
initial.values[nams1] <- init
}
}
# remove objects
rmObjs <- c(names(x), "y.long", "mfX", "mfZ", "data.id2")
rm(list = rmObjs); gc()
# joint model fit
out <- switch(method,
"Cox-PH-GH" = phGH.fit(x, y, id, initial.values, parameterization, derivForm, con),
"weibull-AFT-GH" = weibullAFTGH.fit(x, y, id, initial.values, scaleWB, parameterization, derivForm, con),
"weibull-PH-GH" = weibullPHGH.fit(x, y, id, initial.values, scaleWB, parameterization, derivForm, con),
"piecewise-PH-GH" = piecewisePHGH.fit(x, y, id, initial.values, parameterization, derivForm, con),
"spline-PH-GH" = splinePHGH.fit(x, y, id, initial.values, parameterization, derivForm, con),
"ch-Laplace" = chLaplace.fit(x, y, id, initial.values, b, parameterization, derivForm, con))
# check for problems with the Hessian at convergence
H <- out$Hessian
if (any(is.na(H) | !is.finite(H))) {
warning("infinite or missing values in Hessian at convergence.\n")
} else {
ev <- eigen(H, symmetric = TRUE, only.values = TRUE)$values
if (!all(ev >= -1e-06 * abs(ev[1])))
warning("Hessian matrix at convergence is not positive definite.\n")
}
out$coefficients <- out$coefficients[!sapply(out$coefficients, is.null)]
out$x <- x
out$y <- y
out$times <- data[[timeVar]]
out$data <- data
out$data.id <- data.id
out$method <- method
out$termsYx <- TermsX
out$termsYz <- TermsZ
if (parameterization %in% c("slope", "both")) {
out$termsYx.deriv <- TermsX.deriv
out$termsYz.deriv <- TermsZ.deriv
}
out$termsT <- survObject$terms
out$formYx <- formYx
out$formYz <- formYz
out$formT <- formT
out$timeVar <- timeVar
out$control <- con
out$parameterization <- parameterization
out$derivForm <- derivForm
out$interFact <- interFact
out$CompRisk <- CompRisk
out$LongFormat <- LongFormat
out$assignY <- attr(lmeObject$fixDF, "assign")[-1]
out$assignT <- survObject$assign
out$call <- cl
class(out) <- "jointModel"
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.