R/EMalgRandprop.R
In mesudell/joineRmeta: Joint Modelling for Meta-Analytic (Multi-Study) Data
Defines functions
EMalgRandprop
Documented in
EMalgRandprop
#' EM algorithm function used in jointmeta1
#'
#' Function to run EM algorithm during one stage model fit. Used when the
#' jointmeta1 function is called.
#'
#' @param data the original \code{jointdata} as supplied to the
#' \code{jointmeta1}] call
#' @param longdat the longitudinal data with factors and interaction terms
#' expanded, ordered by increasing survival time
#' @param survdat the survival data with factors and interaction terms expanded,
#' ordered by increasing survival time
#' @param id.name character string specifying the id variable in the dataset
#' @param time.long the name of the variable holding the longitudinal time
#' covariate
#' @param paraests a list of the estimates present from the separate
#' longitudinal and survival fits. Same structure as \code{sepests} if
#' requested in a \code{\link{jointmeta1.object}}
#' @param studies the names of the studies present in the supplied data
#' @param p1 the number of fixed effects included in the longitudinal sub-model
#' @param p2 the number of fixed effects included in the survival sub-model
#' @param q the number of individual level random effects
#' @param r the number of study level random effects, set to \code{NULL} if no
#' study level random effects included in the model
#' @inheritParams jointmeta1
#' @inheritParams longst
#'
#' @return This function returns a list of the estimates of parameters and other
#' information from the run of the EM algorithm. The list has the following
#' components: \describe{
#'
#' \item{\code{beta1}}{a data frame containing the estimates of the fixed
#' effect parameters from the longitudinal sub-model.}
#'
#' \item{\code{beta2}}{a data frame containing the estimates of the fixed
#' effect parameters from the survival sub-model.}
#'
#' \item{\code{sigma.e}}{the estimate of the variance of the measurement
#' errors.}
#'
#' \item{\code{haz}}{the estimated baseline hazard. If \code{strat = TRUE} in
#' the function call to \code{jointmeta1} then this is a list of length equal
#' to the number of studies in the supplied dataset, each element of the list
#' being the baseline hazard for the corresponding study. Otherwise there is a
#' common baseline across all studies in the dataset and this is one vector.}
#'
#' \item{\code{D}}{the estimated covariance matrix for the individual level
#' random effects}
#'
#' \item{\code{A}}{the estimated covariance matrix for the study level random
#' effects. This is only present in the output if study level random effects
#' were specified in the function call to \code{jointmeta1}.}
#'
#' \item{\code{random2}}{a list of matrices containing the conditional modes
#' of the individual level random effects given the supplied data and the
#' estimated parameters of the joint model. The list is of length equal to the
#' number of studies in the dataset, and each element of the list has number
#' of rows equal to the number of individuals in the study, and number of
#' columns equal to the number of specified individual level random effects.}
#'
#' \item{\code{random3}}{a matrix containing the conditional modes of the
#' study level random effects given the supplied data and the estimated
#' parameters of the joint model. The matrix has number of rows equal to the
#' number of studies, and number of columns equal to the number of specified
#' study level random effects.}
#'
#' \item{\code{long.rand.ind.form}}{a character string giving the formulation
#' of the individual level random effects.}
#'
#' \item{\code{long.rand.stud.form}}{a character string giving the formulation
#' of the study level random effects if included in the model.}
#'
#' \item{\code{conv}}{a logical value indicating whether convergence of the EM
#' algorithm was achieved or not.}
#'
#' \item{\code{iters}}{the number of iterations completed by the EM algorithm}
#'
#' \item{\code{n.bystudy}}{the number of individuals present in each study in
#' the data supplied to the function.} }
#'
#' @seealso \code{\link{jointmeta1}}, \code{\link{tojointdata}},
#' \code{\link[joineR]{jointdata}},\code{\link[lme4]{lmer}},
#' \code{\link[survival]{coxph}}
#'
#' @keywords internal
#'
#' @import survival
#' @importFrom stats as.formula model.frame model.matrix
#' @importFrom statmod gauss.quad.prob
#'
EMalgRandprop <- function(data, longdat, survdat, long.rand.ind, long.rand.stud = NULL,
id.name, study.name, gpt, max.it, tol, time.long, surv.formula, long.formula,
long.formula.orig, paraests, studies, p1, p2, strat, print.detail,
bootrun = FALSE, q, r = NULL) {
numstudies <- length(studies)
ids.surv <- survdat[, 1]
ids.bystudy <- lapply(1:numstudies, function(u) {
ids.surv[which(survdat[, 4] == studies[u])]
})
names(ids.bystudy) <- studies
id.long <- longdat[, 1]
id.long.bystudy <- lapply(1:numstudies, function(u) {
id.long[which(id.long %in% ids.bystudy[[u]])]
})
Y <- longdat[, 2]
Y.bystudy <- lapply(1:numstudies, function(u) {
Y[id.long %in% ids.bystudy[[u]]]
})
X1 <- as.matrix(longdat[, 5:ncol(longdat)])
X1.bystudy <- lapply(1:numstudies, function(u) {
X1[id.long %in% ids.bystudy[[u]], ]
})
n <- nrow(survdat)
n.bystudy <- sapply(ids.bystudy, length)
names(n.bystudy) <- studies
s <- survdat[, 2]
s.bystudy <- lapply(1:numstudies, function(u) {
s[ids.surv %in% ids.bystudy[[u]]]
})
names(s.bystudy) <- studies
cen <- survdat[, 3]
cen.bystudy <- lapply(1:numstudies, function(u) {
cen[ids.surv %in% ids.bystudy[[u]]]
})
names(cen.bystudy) <- studies
X2 <- rep(0, numstudies)
if (p2 > 0) {
X2 <- as.matrix(survdat[, 5:dim(survdat)[2]])
X2.bystudy <- lapply(1:numstudies, function(u) {
as.matrix(X2[ids.surv %in% ids.bystudy[[u]], ])
})
} else {
beta2x <- matrix(0, n, 1)
beta2x.bystudy <- lapply(1:numstudies, function(u) {
as.matrix(rep(0, n.bystudy[[u]]))
})
}
beta1 <- paraests$beta1[, 1]
names(beta1) <- rownames(paraests$beta1)
D <- as.matrix(paraests$D)
sigma.e <- paraests$sigma.e
if (is.null(long.rand.stud) == FALSE) {
beta2 <- c(paraests$beta2, 0, 0)
} else {
beta2 <- c(paraests$beta2, 0)
}
haz <- paraests$haz
sf <- paraests$sf
rs <- paraests$rs
nev <- paraests$nev
if (strat) {
rs <- rs[match(names(ids.bystudy), names(rs))]
sf <- sf[match(names(ids.bystudy), names(sf))]
haz <- haz[match(names(ids.bystudy), names(haz))]
nev <- nev[match(names(ids.bystudy), names(nev))]
}
nn <- lapply(1:numstudies, function(u) {
nn <- diff(match(unique(id.long.bystudy[[u]]), id.long.bystudy[[u]]))
nn <- c(nn, length(id.long.bystudy[[u]]) - sum(nn))
})
N <- sum(do.call(sum, nn))
g <- gauss.quad.prob(gpt, "normal", sigma = sqrt(0.5))
ab <- g$nodes
w2 <- g$weights * sqrt(pi)
gmat2 <- matrix(0, gpt^q, q)
gmat2[, 1] <- rep(ab, each = gpt^(q - 1))
if (q > 1) {
gmat2[, 2] <- rep(ab, gpt)
w2 <- as.vector(w2 %x% w2)
if (q > 2) {
gmat2[, 3] <- rep(ab, each = gpt)
w2 <- as.vector(w2 %x% g$weights * sqrt(pi))
}
}
EU.2 <- lapply(1:numstudies, function(u) {
matrix(0, n.bystudy[[u]], q)
})
EUU.2 <- lapply(1:numstudies, function(u) {
matrix(0, n.bystudy[[u]], sum(1:q))
})
EexpU.2 <- lapply(1:numstudies, function(u) {
if (strat) {
matrix(0, n.bystudy[[u]], length(haz[[u]]))
} else {
matrix(0, n.bystudy[[u]], length(haz))
}
})
if ("(Intercept)" %in% long.rand.ind) {
long.rand.ind2 <- long.rand.ind
long.rand.ind2[which(long.rand.ind2 == "(Intercept)")] <- "1"
long.rand.ind.form <- paste(long.rand.ind2, collapse = "+")
}
if ("noint" %in% long.rand.ind) {
long.rand.ind2 <- long.rand.ind[-which(long.rand.ind == "noint")]
long.rand.ind.form <- paste("-1", paste(long.rand.ind2, collapse = "+"),
sep = "+")
}
Z2.form <- as.formula(paste("~", long.rand.ind.form, sep = ""))
Z2.frame <- model.frame(Z2.form, data = longdat)
tZ2 <- model.matrix(Z2.form, Z2.frame)
Z2 <- t(tZ2)
tZ2.bystudy <- lapply(1:numstudies, function(u) {
tZ2[which(id.long %in% ids.bystudy[[u]]), ]
})
Z2.bystudy <- lapply(1:numstudies, function(u) {
Z2[, which(id.long %in% ids.bystudy[[u]])]
})
s1.2 <- rep(1:(q - 1), (q - 1):1)
s2.2 <- sequence((q - 1):1) + rep(1:(q - 1), (q - 1):1)
Z2.form.surv <- as.formula(gsub(time.long, names(survdat)[2], Z2.form))
survdat$survdatorder <- 1:nrow(survdat)
tempsurv <- survdat[, which(!(names(survdat) %in% names(longdat)))]
tempsurv <- cbind(survdat[, 1], tempsurv)
names(tempsurv)[1] <- names(survdat)[1]
tempdat <- merge(tempsurv, longdat[match(unique(longdat[, which(names(longdat) %in%
id.name)]), longdat[, which(names(longdat) %in% id.name)]), c(1,
5:ncol(longdat))], by = id.name)
tempdat <- tempdat[order(tempdat$survdatorder), ]
tempdat <- tempdat[, -which(names(tempdat) %in% "survdatorder")]
survdat <- survdat[, -which(names(survdat) %in% "survdatorder")]
Z2.frame.surv <- model.frame(Z2.form.surv, data = tempdat)
tZ2.surv <- model.matrix(Z2.form.surv, Z2.frame.surv)
tZ2.surv.bystudy <- lapply(1:numstudies, function(u) {
tZ2.surv[which(ids.surv %in% ids.bystudy[[u]]), ]
})
Z2.event.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[[u]], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
if (q > 1) {
temp <- matrix(rep(tZ2.bystudy[[u]][v, ], rstemp), ncol = q,
byrow = TRUE)
} else {
temp <- matrix(rep(tZ2.bystudy[[u]][v], rstemp), ncol = q,
byrow = TRUE)
}
colnames(temp) <- colnames(tZ2)
if (time.long %in% colnames(temp)) {
if (strat) {
temp[, colnames(temp) == time.long] <- sf[[u]][1:rstemp]
} else {
temp[, colnames(temp) == time.long] <- sf[1:rstemp]
}
}
t(temp)
} else {
matrix(rep(0, q))
}
})
})
Z2.event.sqr.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
temp <- (Z2.event.bystudy[[u]][[v]])^2
if (q > 1) {
if (rstemp == 1) {
temp <- rbind(temp, matrix(Z2.event.bystudy[[u]][[v]][s1.2,
] * Z2.event.bystudy[[u]][[v]][s2.2, ]))
} else {
temp <- rbind(temp, Z2.event.bystudy[[u]][[v]][s1.2,
] * Z2.event.bystudy[[u]][[v]][s2.2, ])
}
}
temp
} else {
matrix(rep(0, sum(1:q)))
}
})
})
Z2b2event <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = q)
} else {
matrix(0, nrow = length(haz), ncol = q)
}
})
})
Z2b2eventsqr <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = sum(1:q))
} else {
matrix(0, nrow = length(haz), ncol = sum(1:q))
}
})
})
newu.2.all <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
idtemp <- ids.bystudy[[u]][v]
B2temp <- paraests$randstart.ind.cov[[which(names(paraests$randstart.ind.cov) %in%
idtemp)]]
cm.2 <- matrix(as.numeric(as.character(paraests$randstart.ind[which(rownames(paraests$randstart.ind) %in%
idtemp), ])), gpt^q, q, TRUE)
gmat2 %*% solve(chol(solve(B2temp))) + cm.2
})
})
randstart.ind <- paraests$randstart.ind
randstart.ind.bystudy <- lapply(1:numstudies, function(u) {
out <- as.matrix(randstart.ind[which(rownames(randstart.ind) %in%
ids.bystudy[[u]]), ])
rownames(out) <- rownames(randstart.ind)[which(rownames(randstart.ind) %in%
ids.bystudy[[u]])]
out
})
randstart.ind.bystudy <- lapply(1:numstudies, function(u) {
randstart.ind.bystudy[[u]][match(ids.bystudy[[u]], rownames(randstart.ind.bystudy[[u]])),
]
})
if (is.null(long.rand.stud) == FALSE) {
A <- as.matrix(paraests$A)
w3 <- g$weights * sqrt(pi)
gmat3 <- matrix(0, gpt^r, r)
gmat3[, 1] <- rep(ab, each = gpt^(r - 1))
if (r > 1) {
gmat3[, 2] <- rep(ab, gpt)
w3 <- as.vector(w3 %x% w3)
if (r > 2) {
gmat3[, 3] <- rep(ab, each = gpt)
w3 <- as.vector(w3 %x% g$weights * sqrt(pi))
}
}
EU.3 <- matrix(0, numstudies, r)
EUU.3 <- matrix(0, numstudies, sum(1:r))
EexpU.3 <- lapply(1:numstudies, function(u) {
if (strat) {
matrix(0, n.bystudy[[u]], length(haz[[u]]))
} else {
matrix(0, n.bystudy[[u]], length(haz))
}
})
if (study.name %in% long.rand.stud) {
long.rand.stud2 <- long.rand.stud
long.rand.stud2[which(long.rand.stud2 == study.name)] <- "1"
long.rand.stud.form <- paste(long.rand.stud2, collapse = "+")
} else {
long.rand.stud2 <- long.rand.stud
long.rand.stud.form <- paste("-1", paste(long.rand.stud2, collapse = "+"),
sep = "+")
}
Z3.form <- as.formula(paste("~", long.rand.stud.form, sep = ""))
Z3.frame <- model.frame(Z3.form, data = longdat)
tZ3 <- model.matrix(Z3.form, Z3.frame)
Z3 <- t(tZ3)
tZ3.bystudy <- lapply(1:numstudies, function(u) {
tZ3[which(id.long %in% ids.bystudy[[u]]), ]
})
Z3.bystudy <- lapply(1:numstudies, function(u) {
Z3[, which(id.long %in% ids.bystudy[[u]])]
})
s1.3 <- rep(1:(r - 1), (r - 1):1)
s2.3 <- sequence((r - 1):1) + rep(1:(r - 1), (r - 1):1)
Z3.form.surv <- as.formula(gsub(time.long, names(survdat)[2], Z3.form))
Z3.frame.surv <- model.frame(Z3.form.surv, data = tempdat)
tempdat <- NULL
tZ3.surv <- model.matrix(Z3.form.surv, Z3.frame.surv)
Z3.surv <- t(tZ3.surv)
tZ3.surv.bystudy <- lapply(1:numstudies, function(u) {
tZ3.surv[which(ids.surv %in% ids.bystudy[[u]]), ]
})
Z3.surv.bystudy <- lapply(1:numstudies, function(u) {
Z3.surv[, which(ids.surv %in% ids.bystudy[[u]])]
})
Z3.event.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[[u]], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
if (r > 1) {
temp <- matrix(rep(tZ3.bystudy[[u]][v, ], rstemp),
ncol = r, byrow = TRUE)
} else {
temp <- matrix(rep(tZ3.bystudy[[u]][v], rstemp), ncol = r,
byrow = TRUE)
}
colnames(temp) <- colnames(tZ3)
if (time.long %in% colnames(temp)) {
if (strat) {
temp[, colnames(temp) == time.long] <- sf[[u]][1:rstemp]
} else {
temp[, colnames(temp) == time.long] <- sf[1:rstemp]
}
}
t(temp)
} else {
matrix(rep(0, r))
}
})
})
Z3.event.sqr.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
temp <- (Z3.event.bystudy[[u]][[v]])^2
if (r > 1) {
if (rstemp == 1) {
temp <- rbind(temp, matrix(Z3.event.bystudy[[u]][[v]][s1.3] *
Z3.event.bystudy[[u]][[v]][s2.3, ]))
} else {
temp <- rbind(temp, Z3.event.bystudy[[u]][[v]][s1.3,
] * Z3.event.bystudy[[u]][[v]][s2.3, ])
}
}
temp
} else {
matrix(rep(0, sum(1:r)))
}
})
})
randselector1 <- rep(1:r, each = q)
randselector2 <- rep(1:q, times = r)
Z3Z2.event.sqr.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
Z3.event.bystudy[[u]][[v]][randselector1, ] * Z2.event.bystudy[[u]][[v]][randselector2,
]
} else {
matrix(rep(0, length(randselector2)))
}
})
})
Z3b3event <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = r)
} else {
matrix(0, nrow = length(haz), ncol = r)
}
})
})
Z3b3eventsqr <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = sum(1:r))
} else {
matrix(0, nrow = length(haz), ncol = sum(1:r))
}
})
})
newu.3.all <- lapply(1:numstudies, function(u) {
studtemp <- studies[u]
B3temp <- paraests$randstart.stud.cov[[which(names(paraests$randstart.stud.cov) %in%
studtemp)]]
cm.3 <- matrix(as.numeric(as.character(paraests$randstart.stud[which(rownames(paraests$randstart.stud) %in%
studtemp), ])), gpt^r, r, TRUE)
gmat3 %*% solve(chol(solve(B3temp))) + cm.3
})
names(newu.3.all) <- studies
randstart.stud <- paraests$randstart.stud
survnames <- NULL
if (p2 > 0) {
survnames <- paste("T.", names(beta2)[1:(length(beta2) - 2)],
sep = "")
}
paranames <- c(paste("Y.", names(beta1), sep = ""), survnames,
"gamma_ind_0", "gamma_stud_0", "sigma.e", paste("D", paste(rep(1:q,
times = q), rep(1:q, each = q), sep = ""), sep = ""), paste("A",
paste(rep(1:r, times = r), rep(1:r, each = r), sep = ""),
sep = ""))
} else {
r <- NULL
survnames <- NULL
if (p2 > 0) {
survnames <- paste("T.", names(beta2)[1:(length(beta2) - 1)],
sep = "")
}
paranames <- c(paste("Y.", names(beta1), sep = ""), survnames,
"gamma_ind_0", "sigma.e", paste("D", paste(rep(1:q, times = q),
rep(1:q, each = q), sep = ""), sep = ""))
}
conv <- FALSE
for (it in 1:max.it) {
if (p2 > 0) {
beta2x <- lapply(1:numstudies, function(u) {
X2.bystudy[[u]] %*% beta2[1:p2]
})
} else {
beta2x <- lapply(1:numstudies, function(u) {
matrix(rep(0, n.bystudy[u]))
})
}
ebeta2x <- lapply(1:numstudies, function(u) {
exp(beta2x[[u]])
})
for (k in 1:numstudies) {
if (is.null(long.rand.stud) == FALSE) {
if (it == 1) {
b2est <- randstart.ind.bystudy[[k]]
} else {
b2est <- EU.2[[k]]
}
if (it == 1) {
b3est <- as.numeric(randstart.stud[k, ])
} else {
b3est <- as.numeric(EU.3[k, ])
}
newu.3 <- newu.3.all[[k]]
newu.3sqr <- newu.3^2
if (r > 1) {
newu.3sqr <- cbind(newu.3sqr, newu.3[, s1.3] * newu.3[,
s2.3])
}
egDUs.3 <- t(apply(exp(newu.3 %*% (Z3.surv.bystudy[[k]] *
beta2[(p2 + 2)])), 1, function(u) {
u^cen.bystudy[[k]]
}))
egDUsf.3 <- lapply(1:n.bystudy[k], function(v) {
exp(newu.3 %*% (Z3.event.bystudy[[k]][[v]] * beta2[(p2 +
2)]))
})
if (q > 1) {
egDUsf.2 <- lapply(1:n.bystudy[k], function(v) {
exp(as.numeric(b2est[v, ]) %*% (Z2.event.bystudy[[k]][[v]] *
beta2[(p2 + 1)]))
})
} else {
egDUsf.2 <- lapply(1:n.bystudy[k], function(v) {
exp(as.numeric(b2est[v]) %*% (Z2.event.bystudy[[k]][[v]] *
beta2[(p2 + 1)]))
})
}
ess.3 <- lapply(1:n.bystudy[k], function(v) {
if (strat) {
exp(-(ebeta2x[[k]][v, ] * egDUsf.3[[v]]) %*% ((haz[[k]][1:rs[[k]][which(names(rs[[k]]) %in%
ids.bystudy[[k]][[v]])]] * as.numeric(egDUsf.2[[v]]))))
} else {
exp(-(ebeta2x[[k]][v, ] * egDUsf.3[[v]]) %*% ((haz[1:rs[which(names(rs) %in%
ids.bystudy[[k]][[v]])]] * as.numeric(egDUsf.2[[v]]))))
}
})
f.3 <- rowSums(do.call(cbind, lapply(1:n.bystudy[k], function(v) {
rowSums(egDUs.3[, v] * ess.3[[v]] * w3)
})))
den.3 <- sum(f.3)
EU.3[k, 1:r] <- f.3 %*% newu.3/den.3
EUU.3[k, 1:sum(1:r)] <- f.3 %*% newu.3sqr/den.3
}
for (i in 1:n.bystudy[k]) {
if (strat) {
rstemp <- rs[[k]][which(names(rs[[k]]) %in% ids.bystudy[[k]][[i]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[k]][[i]])]
}
newu.2 <- newu.2.all[[k]][[i]]
newu.2sqr <- newu.2^2
if (q > 1) {
newu.2sqr <- cbind(newu.2sqr, newu.2[, s1.2] * newu.2[,
s2.2])
}
egDUs.2.ind <- 1
egDUs.3.const <- 1
egDUs.3.ind <- 1
egDUs.2.const <- 1
if (cen.bystudy[[k]][i] == 1) {
if (q > 1) {
egDUs.2.ind <- exp(newu.2 %*% (tZ2.surv.bystudy[[k]][i,
] * beta2[(p2 + 1)]))
} else {
egDUs.2.ind <- exp(newu.2 %*% (tZ2.surv.bystudy[[k]][i] *
beta2[(p2 + 1)]))
}
}
egDUsf.2.ind <- exp(newu.2 %*% (Z2.event.bystudy[[k]][[i]] *
beta2[(p2 + 1)]))
if (is.null(long.rand.stud) == FALSE) {
if (cen.bystudy[[k]][i] == 1) {
if (r > 1) {
egDUs.3.const <- exp(b3est %*% (tZ3.surv.bystudy[[k]][i,
] * beta2[(p2 + 2)]))
egDUs.3.ind <- exp(newu.3 %*% (tZ3.surv.bystudy[[k]][i,
] * beta2[(p2 + 2)]))
} else {
egDUs.3.const <- exp(b3est %*% (tZ3.surv.bystudy[[k]][i] *
beta2[(p2 + 2)]))
egDUs.3.ind <- exp(newu.3 %*% (tZ3.surv.bystudy[[k]][i] *
beta2[(p2 + 2)]))
}
if (q > 1) {
egDUs.2.const <- exp(b2est[i, ] %*% (tZ2.surv.bystudy[[k]][i,
] * beta2[(p2 + 1)]))
} else {
egDUs.2.const <- exp(b2est[i] %*% (tZ2.surv.bystudy[[k]][i] *
beta2[(p2 + 1)]))
}
}
if (q > 1) {
egDUsf.2.const <- exp(b2est[i, ] %*% (Z2.event.bystudy[[k]][[i]] *
beta2[(p2 + 1)]))
} else {
egDUsf.2.const <- exp(b2est[i] %*% (Z2.event.bystudy[[k]][[i]] *
beta2[(p2 + 1)]))
}
egDUsf.3.const <- exp(b3est %*% (Z3.event.bystudy[[k]][[i]] *
beta2[(p2 + 2)]))
egDUsf.3.ind <- exp(newu.3 %*% (Z3.event.bystudy[[k]][[i]] *
beta2[(p2 + 2)]))
if (strat) {
ess.3.ind <- exp(-(ebeta2x[[k]][i, ] * egDUsf.3.ind) %*%
t(haz[[k]][1:rstemp] * egDUsf.2.const))
} else {
ess.3.ind <- exp(-(ebeta2x[[k]][i, ] * egDUsf.3.ind) %*%
t(haz[1:rstemp] * egDUsf.2.const))
}
f.3.ind <- egDUs.3.ind * ess.3.ind * w3
den.3.ind <- sum(f.3.ind)
C.3 <- egDUsf.3.ind[, 1:rstemp]
EexpU.3[[k]][i, 1:rstemp] <- f.3.ind[, 1] %*% C.3/den.3.ind
Z3b3event[[k]][[i]][1:ncol(Z3.event.bystudy[[k]][[i]]),
] <- (t(Z3.event.bystudy[[k]][[i]] * as.numeric(f.3.ind[,
1] %*% newu.3)))/den.3.ind
Z3b3eventsqr[[k]][[i]][1:ncol(Z3.event.sqr.bystudy[[k]][[i]]),
] <- (as.numeric(f.3.ind[, 1] %*% newu.3sqr) * Z3.event.sqr.bystudy[[k]][[i]])/den.3.ind
}
if (strat) {
if (is.null(long.rand.stud) == FALSE) {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
t(haz[[k]][1:rstemp] * egDUsf.3.const))
} else {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
(haz[[k]][1:rstemp]))
}
} else {
if (is.null(long.rand.stud) == FALSE) {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
t(haz[1:rstemp] * egDUsf.3.const))
} else {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
haz[1:rstemp])
}
}
f.2 <- egDUs.2.ind * ess.2 * w2
den.2 <- sum(f.2)
EU.2[[k]][i, 1:q] <- f.2[, 1] %*% newu.2/den.2
EUU.2[[k]][i, 1:sum(1:q)] <- f.2[, 1] %*% newu.2sqr/den.2
C.2 <- egDUsf.2.ind[, 1:rstemp]
EexpU.2[[k]][i, 1:rstemp] <- f.2[, 1] %*% C.2/den.2
Z2b2event[[k]][[i]][1:ncol(Z2.event.bystudy[[k]][[i]]),
] <- (t(Z2.event.bystudy[[k]][[i]] * as.numeric(f.2[,
1] %*% newu.2)))/den.2
Z2b2eventsqr[[k]][[i]][1:ncol(Z2.event.sqr.bystudy[[k]][[i]]),
] <- (as.numeric(f.2[, 1] %*% newu.2sqr) * Z2.event.sqr.bystudy[[k]][[i]])/den.2
}
}
parac <- data.frame(c(beta1, beta2, sigma.e, D))
if (is.null(long.rand.stud) == FALSE) {
parac <- data.frame(c(beta1, beta2, sigma.e, D, A))
}
EexpU <- lapply(1:numstudies, function(u) {
if (is.null(long.rand.stud) == FALSE) {
EexpU.2[[u]] * EexpU.3[[u]]
} else {
EexpU.2[[u]]
}
})
if (strat == FALSE) {
haz <- nev/colSums(do.call(rbind, EexpU) * do.call(rbind, ebeta2x)[,
1])
} else {
haz <- lapply(1:numstudies, function(u) {
nev[[u]]/colSums(EexpU[[u]] * ebeta2x[[u]][, 1])
})
}
EUmat.2 <- lapply(1:numstudies, function(u) {
apply(EU.2[[u]], 2, rep, nn[[u]])
})
EUUmat.2 <- lapply(1:numstudies, function(u) {
apply(EUU.2[[u]], 2, rep, nn[[u]])
})
Ut.2 <- lapply(1:numstudies, function(u) {
rowSums(EUmat.2[[u]] * tZ2.bystudy[[u]])
})
if (is.null(long.rand.stud) == FALSE) {
EUmat.3 <- lapply(1:numstudies, function(u) {
matrix(rep(EU.3[u, ], sum(nn[[u]])), ncol = r)
})
EUUmat.3 <- lapply(1:numstudies, function(u) {
matrix(rep(EUU.3[u, ], sum(nn[[u]])), ncol = sum(1:r))
})
Ut.3 <- lapply(1:numstudies, function(u) {
rowSums(EUmat.3[[u]] * tZ3.bystudy[[u]])
})
beta1 <- as.numeric(solve(crossprod(do.call(rbind, X1.bystudy)),
crossprod(do.call(rbind, X1.bystudy), do.call(c, Y.bystudy) -
(do.call(c, Ut.2) + do.call(c, Ut.3)))))
if (q == 1) {
if (r == 1) {
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (do.call(c, Z2.bystudy) *
do.call(c, lapply(1:numstudies, function(u) {
rep(EU.2[[u]], nn[[u]])
}))) + (do.call(c, Z3.bystudy) * do.call(c, lapply(1:numstudies,
function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]])]
})))))^2)/sum(N)
} else {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z3.bystudy)
out2 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (do.call(c, Z2.bystudy) *
do.call(c, lapply(1:numstudies, function(u) {
rep(EU.2[[u]], nn[[u]])
}))) + out2))^2)/sum(N)
}
} else {
if (r == 1) {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.2[[u]][rep(1:nrow(EU.2[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z2.bystudy)
out1 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + out1 + (do.call(c, Z3.bystudy) *
do.call(c, lapply(1:numstudies, function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]])]
})))))^2)/sum(N)
} else {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.2[[u]][rep(1:nrow(EU.2[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z2.bystudy)
out1 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z3.bystudy)
out2 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + out1 + out2))^2)/sum(N)
}
}
diag(A) <- colMeans(EUU.3)[1:r]
if (r > 1) {
A[lower.tri(A)] <- colMeans(EUU.3)[-(1:r)]
A[upper.tri(A)] <- t(A)[upper.tri(A)]
}
} else {
beta1 <- as.numeric(solve(crossprod(do.call(rbind, X1.bystudy)),
crossprod(do.call(rbind, X1.bystudy), do.call(c, Y.bystudy) -
(do.call(c, Ut.2)))))
if (q == 1) {
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (do.call(c, Z2.bystudy) * do.call(c,
lapply(1:numstudies, function(u) {
rep(EU.2[[u]], nn[[u]])
})))))^2)/sum(N)
} else {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.2[[u]][rep(1:nrow(EU.2[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z2.bystudy)
out <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (out)))^2)/sum(N)
}
}
diag(D) <- colMeans(do.call(rbind, EUU.2))[1:q]
if (q > 1) {
D[lower.tri(D)] <- colMeans(do.call(rbind, EUU.2))[-(1:q)]
D[upper.tri(D)] <- t(D)[upper.tri(D)]
}
if (is.null(long.rand.stud) == FALSE) {
EexpU <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
EexpU.2[[u]][v, ] * EexpU.3[[u]][v, ]
})
})
fd <- vector("numeric", p2 + q + r)
sd <- matrix(0, p2 + q + r, p2 + q + r)
} else {
EexpU <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
EexpU.2[[u]][v, ]
})
})
fd <- vector("numeric", p2 + q)
sd <- matrix(0, p2 + q, p2 + q)
}
if (q > 1) {
fd[(p2 + 1):(p2 + q)] <- colSums(do.call(c, cen.bystudy) *
(do.call(rbind, EU.2) * do.call(rbind, tZ2.surv.bystudy))) -
colSums(do.call(rbind, ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
fd[(p2 + 1):(p2 + q)] <- sum(do.call(c, cen.bystudy) * (do.call(rbind,
EU.2) * do.call(c, tZ2.surv.bystudy))) - colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
})))
}
if (is.null(long.rand.stud) == FALSE) {
if (r > 1) {
fd[(p2 + q + 1):(p2 + q + r)] <- colSums(do.call(c, cen.bystudy) *
(EU.3[rep(1:nrow(EU.3), n.bystudy), ]) * do.call(rbind,
tZ3.surv.bystudy)) - colSums(do.call(rbind, ebeta2x)[,
1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
fd[(p2 + q + 1):(p2 + q + r)] <- sum(do.call(c, cen.bystudy) *
(EU.3[rep(1:nrow(EU.3), n.bystudy), ]) * do.call(c, tZ3.surv.bystudy)) -
colSums(do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]] * EexpU[[u]][[v]] *
haztemp)
}))
})))
}
}
if (q > 1) {
if (q == 2) {
sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)][upper.tri(sd[(p2 +
1):(p2 + q), (p2 + 1):(p2 + q)])] <- -sum(do.call(rbind,
ebeta2x)[, 1] * 0.5 * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z2b2eventsqr[[u]][[v]][, (q + 1):sum(1:q)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)][upper.tri(sd[(p2 +
1):(p2 + q), (p2 + 1):(p2 + q)])] <- -colSums(do.call(rbind,
ebeta2x)[, 1] * 0.5 * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z2b2eventsqr[[u]][[v]][, (q + 1):sum(1:q)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
}
}
if (is.null(long.rand.stud) == FALSE) {
if (r > 1) {
if (r == 2) {
sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 + q +
r)][upper.tri(sd[(p2 + q + 1):(p2 + q + r), (p2 + q +
1):(p2 + q + r)])] <- -sum(do.call(rbind, ebeta2x)[,
1] * 0.5 * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z3b3eventsqr[[u]][[v]][, (r + 1):sum(1:r)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 + q +
r)][upper.tri(sd[(p2 + q + 1):(p2 + q + r), (p2 + q +
1):(p2 + q + r)])] <- -colSums(do.call(rbind, ebeta2x)[,
1] * 0.5 * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z3b3eventsqr[[u]][[v]][, (r + 1):sum(1:r)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
}
}
if (q == 1 && r == 1) {
sd[(p2 + 1):(p2 + q), ((p2 + q + 1):(p2 + q + r))] <- -colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z3b3event[[u]][[v]] * Z2b2event[[u]][[v]] * EexpU.2[[u]][v,
] * EexpU.3[[u]][v, ] * haztemp)
}))
})))
} else {
sd[(p2 + 1):(p2 + q), ((p2 + q + 1):(p2 + q + r))] <- -colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]][, randselector1] * Z2b2event[[u]][[v]][,
randselector2] * EexpU[[u]][[v]] * haztemp)
}))
})))
}
}
if (p2 > 0) {
fd[1:p2] <- c(colSums((do.call(c, cen.bystudy) * do.call(rbind,
X2.bystudy)) - (do.call(rbind, X2.bystudy) * as.numeric(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(EexpU[[u]][[v]] * haztemp)
}))
}))))))
sd[(1:p2), (p2 + 1):(p2 + q)] <- -(t(do.call(rbind, X2.bystudy)) %*%
(do.call(rbind, ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))
if (is.null(long.rand.stud) == FALSE) {
sd[(1:p2), (p2 + q + 1):(p2 + q + r)] <- -(t(do.call(rbind,
X2.bystudy)) %*% (do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))
}
sd <- sd + t(sd)
for (i in 1:p2) {
for (j in 1:p2) {
sd[i, j] <- -(sum(do.call(rbind, X2.bystudy)[, i] * do.call(rbind,
X2.bystudy)[, j] * do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(EexpU[[u]][[v]] * haztemp)
}))
}))))
}
}
}
if (q == 1) {
sd[(p2 + 1), (p2 + 1)] <- -(colSums(do.call(rbind, ebeta2x)[,
1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2eventsqr[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))[1:q]
} else {
diag(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)]) <- -(colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2eventsqr[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))[1:q]
}
if (is.null(long.rand.stud) == FALSE) {
if (r == 1) {
sd[(p2 + q + 1), (p2 + q + 1)] <- -(colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3eventsqr[[u]][[v]] * EexpU[[u]][[v]] *
haztemp)
}))
}))))[1:r]
} else {
diag(sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 + q +
r)]) <- -(colSums(do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3eventsqr[[u]][[v]] * EexpU[[u]][[v]] *
haztemp)
}))
}))))[1:r]
}
}
if (is.null(long.rand.stud) == FALSE) {
if (r == 1 && q == 1) {
fd <- fd
sd <- sd
} else {
if (p2 > 0) {
fd <- c(fd[1:p2], sum(fd[(p2 + 1):(p2 + q)]), sum(fd[(p2 +
q + 1):(p2 + q + r)]))
sd.p2 <- sd[1:p2, 1:p2]
sd.q <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
sd.r <- sum(sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 +
q + r)])
if (p2 > 1) {
if (q > 1) {
sd.p2q <- rowSums(sd[1:p2, (p2 + 1):(p2 + q)])
} else {
sd.p2q <- sd[1:p2, (p2 + 1):(p2 + q)]
}
if (r > 1) {
sd.p2r <- rowSums(sd[1:p2, (p2 + q + 1):(p2 + q +
r)])
} else {
sd.p2r <- sd[1:p2, (p2 + q + 1):(p2 + q + r)]
}
} else {
sd.p2q <- sum(sd[1:p2, (p2 + 1):(p2 + q)])
sd.p2r <- sum(sd[1:p2, (p2 + q + 1):(p2 + q + r)])
}
sd.qr <- sum(sd[(p2 + 1):(p2 + q), (p2 + q + 1):(p2 +
q + r)])
sd <- rbind(cbind(sd.p2, sd.p2q, sd.p2r), c(sd.p2q, sd.q,
sd.qr), c(sd.p2r, sd.qr, sd.r))
rownames(sd) <- colnames(sd) <- NULL
} else {
fd <- c(sum(fd[(p2 + 1):(p2 + q)]), sum(fd[(p2 + q +
1):(p2 + q + r)]))
sd.q <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
sd.r <- sum(sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 +
q + r)])
sd.qr <- sum(sd[(p2 + 1):(p2 + q), (p2 + q + 1):(p2 +
q + r)])
sd <- rbind(c(sd.q, sd.qr), c(sd.qr, sd.r))
rownames(sd) <- colnames(sd) <- NULL
}
}
} else {
if (q == 1) {
fd <- fd
sd <- sd
} else {
if (p2 > 0) {
fd <- c(fd[1:p2], sum(fd[(p2 + 1):(p2 + q)]))
sd.p2 <- sd[1:p2, 1:p2]
sd.q <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
if (p2 > 1) {
sd.p2q <- rowSums(sd[1:p2, (p2 + 1):(p2 + q)])
} else {
sd.p2q <- sum(sd[1:p2, (p2 + 1):(p2 + q)])
}
sd <- rbind(cbind(sd.p2, sd.p2q), c(sd.p2q, sd.q))
rownames(sd) <- colnames(sd) <- NULL
} else {
fd <- c(sum(fd[(p2 + 1):(p2 + q)]))
sd <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
rownames(sd) <- colnames(sd) <- NULL
}
}
}
beta2 <- beta2 - solve(sd, fd)
para <- data.frame(c(beta1, beta2, sigma.e, D))
if (is.null(long.rand.stud) == FALSE) {
para <- data.frame(c(beta1, beta2, sigma.e, D, A))
}
dd <- abs(parac - para)
if (print.detail) {
print(paste("Iteration: ", it, sep = ""))
print("Current parameter estimates:")
detail <- data.frame(para)
rownames(detail) <- paranames
colnames(detail) <- NULL
print(detail)
}
if (max(dd) < tol) {
conv <- TRUE
break
}
}
if (conv != TRUE && bootrun == FALSE) {
print("Not converged")
}
if (strat) {
names(haz) <- studies
}
results <- list(beta1 = data.frame(beta1), beta2 = data.frame(beta2),
sigma.e = sigma.e, D = D, haz = haz, random2 = EU.2, conv = conv,
iters = it, long.rand.ind.form = long.rand.ind.form, n.bystudy = n.bystudy)
if (is.null(long.rand.stud) == FALSE) {
results$A <- A
results$random3 <- EU.3
results$long.rand.stud.form <- long.rand.stud.form
}
results
}
mesudell/joineRmeta documentation built on Jan. 24, 2020, 6:06 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions EMalgRandprop
Documented in EMalgRandprop
#' EM algorithm function used in jointmeta1
#'
#' Function to run EM algorithm during one stage model fit. Used when the
#' jointmeta1 function is called.
#'
#' @param data the original \code{jointdata} as supplied to the
#' \code{jointmeta1}] call
#' @param longdat the longitudinal data with factors and interaction terms
#' expanded, ordered by increasing survival time
#' @param survdat the survival data with factors and interaction terms expanded,
#' ordered by increasing survival time
#' @param id.name character string specifying the id variable in the dataset
#' @param time.long the name of the variable holding the longitudinal time
#' covariate
#' @param paraests a list of the estimates present from the separate
#' longitudinal and survival fits. Same structure as \code{sepests} if
#' requested in a \code{\link{jointmeta1.object}}
#' @param studies the names of the studies present in the supplied data
#' @param p1 the number of fixed effects included in the longitudinal sub-model
#' @param p2 the number of fixed effects included in the survival sub-model
#' @param q the number of individual level random effects
#' @param r the number of study level random effects, set to \code{NULL} if no
#' study level random effects included in the model
#' @inheritParams jointmeta1
#' @inheritParams longst
#'
#' @return This function returns a list of the estimates of parameters and other
#' information from the run of the EM algorithm. The list has the following
#' components: \describe{
#'
#' \item{\code{beta1}}{a data frame containing the estimates of the fixed
#' effect parameters from the longitudinal sub-model.}
#'
#' \item{\code{beta2}}{a data frame containing the estimates of the fixed
#' effect parameters from the survival sub-model.}
#'
#' \item{\code{sigma.e}}{the estimate of the variance of the measurement
#' errors.}
#'
#' \item{\code{haz}}{the estimated baseline hazard. If \code{strat = TRUE} in
#' the function call to \code{jointmeta1} then this is a list of length equal
#' to the number of studies in the supplied dataset, each element of the list
#' being the baseline hazard for the corresponding study. Otherwise there is a
#' common baseline across all studies in the dataset and this is one vector.}
#'
#' \item{\code{D}}{the estimated covariance matrix for the individual level
#' random effects}
#'
#' \item{\code{A}}{the estimated covariance matrix for the study level random
#' effects. This is only present in the output if study level random effects
#' were specified in the function call to \code{jointmeta1}.}
#'
#' \item{\code{random2}}{a list of matrices containing the conditional modes
#' of the individual level random effects given the supplied data and the
#' estimated parameters of the joint model. The list is of length equal to the
#' number of studies in the dataset, and each element of the list has number
#' of rows equal to the number of individuals in the study, and number of
#' columns equal to the number of specified individual level random effects.}
#'
#' \item{\code{random3}}{a matrix containing the conditional modes of the
#' study level random effects given the supplied data and the estimated
#' parameters of the joint model. The matrix has number of rows equal to the
#' number of studies, and number of columns equal to the number of specified
#' study level random effects.}
#'
#' \item{\code{long.rand.ind.form}}{a character string giving the formulation
#' of the individual level random effects.}
#'
#' \item{\code{long.rand.stud.form}}{a character string giving the formulation
#' of the study level random effects if included in the model.}
#'
#' \item{\code{conv}}{a logical value indicating whether convergence of the EM
#' algorithm was achieved or not.}
#'
#' \item{\code{iters}}{the number of iterations completed by the EM algorithm}
#'
#' \item{\code{n.bystudy}}{the number of individuals present in each study in
#' the data supplied to the function.} }
#'
#' @seealso \code{\link{jointmeta1}}, \code{\link{tojointdata}},
#' \code{\link[joineR]{jointdata}},\code{\link[lme4]{lmer}},
#' \code{\link[survival]{coxph}}
#'
#' @keywords internal
#'
#' @import survival
#' @importFrom stats as.formula model.frame model.matrix
#' @importFrom statmod gauss.quad.prob
#'
EMalgRandprop <- function(data, longdat, survdat, long.rand.ind, long.rand.stud = NULL,
id.name, study.name, gpt, max.it, tol, time.long, surv.formula, long.formula,
long.formula.orig, paraests, studies, p1, p2, strat, print.detail,
bootrun = FALSE, q, r = NULL) {
numstudies <- length(studies)
ids.surv <- survdat[, 1]
ids.bystudy <- lapply(1:numstudies, function(u) {
ids.surv[which(survdat[, 4] == studies[u])]
})
names(ids.bystudy) <- studies
id.long <- longdat[, 1]
id.long.bystudy <- lapply(1:numstudies, function(u) {
id.long[which(id.long %in% ids.bystudy[[u]])]
})
Y <- longdat[, 2]
Y.bystudy <- lapply(1:numstudies, function(u) {
Y[id.long %in% ids.bystudy[[u]]]
})
X1 <- as.matrix(longdat[, 5:ncol(longdat)])
X1.bystudy <- lapply(1:numstudies, function(u) {
X1[id.long %in% ids.bystudy[[u]], ]
})
n <- nrow(survdat)
n.bystudy <- sapply(ids.bystudy, length)
names(n.bystudy) <- studies
s <- survdat[, 2]
s.bystudy <- lapply(1:numstudies, function(u) {
s[ids.surv %in% ids.bystudy[[u]]]
})
names(s.bystudy) <- studies
cen <- survdat[, 3]
cen.bystudy <- lapply(1:numstudies, function(u) {
cen[ids.surv %in% ids.bystudy[[u]]]
})
names(cen.bystudy) <- studies
X2 <- rep(0, numstudies)
if (p2 > 0) {
X2 <- as.matrix(survdat[, 5:dim(survdat)[2]])
X2.bystudy <- lapply(1:numstudies, function(u) {
as.matrix(X2[ids.surv %in% ids.bystudy[[u]], ])
})
} else {
beta2x <- matrix(0, n, 1)
beta2x.bystudy <- lapply(1:numstudies, function(u) {
as.matrix(rep(0, n.bystudy[[u]]))
})
}
beta1 <- paraests$beta1[, 1]
names(beta1) <- rownames(paraests$beta1)
D <- as.matrix(paraests$D)
sigma.e <- paraests$sigma.e
if (is.null(long.rand.stud) == FALSE) {
beta2 <- c(paraests$beta2, 0, 0)
} else {
beta2 <- c(paraests$beta2, 0)
}
haz <- paraests$haz
sf <- paraests$sf
rs <- paraests$rs
nev <- paraests$nev
if (strat) {
rs <- rs[match(names(ids.bystudy), names(rs))]
sf <- sf[match(names(ids.bystudy), names(sf))]
haz <- haz[match(names(ids.bystudy), names(haz))]
nev <- nev[match(names(ids.bystudy), names(nev))]
}
nn <- lapply(1:numstudies, function(u) {
nn <- diff(match(unique(id.long.bystudy[[u]]), id.long.bystudy[[u]]))
nn <- c(nn, length(id.long.bystudy[[u]]) - sum(nn))
})
N <- sum(do.call(sum, nn))
g <- gauss.quad.prob(gpt, "normal", sigma = sqrt(0.5))
ab <- g$nodes
w2 <- g$weights * sqrt(pi)
gmat2 <- matrix(0, gpt^q, q)
gmat2[, 1] <- rep(ab, each = gpt^(q - 1))
if (q > 1) {
gmat2[, 2] <- rep(ab, gpt)
w2 <- as.vector(w2 %x% w2)
if (q > 2) {
gmat2[, 3] <- rep(ab, each = gpt)
w2 <- as.vector(w2 %x% g$weights * sqrt(pi))
}
}
EU.2 <- lapply(1:numstudies, function(u) {
matrix(0, n.bystudy[[u]], q)
})
EUU.2 <- lapply(1:numstudies, function(u) {
matrix(0, n.bystudy[[u]], sum(1:q))
})
EexpU.2 <- lapply(1:numstudies, function(u) {
if (strat) {
matrix(0, n.bystudy[[u]], length(haz[[u]]))
} else {
matrix(0, n.bystudy[[u]], length(haz))
}
})
if ("(Intercept)" %in% long.rand.ind) {
long.rand.ind2 <- long.rand.ind
long.rand.ind2[which(long.rand.ind2 == "(Intercept)")] <- "1"
long.rand.ind.form <- paste(long.rand.ind2, collapse = "+")
}
if ("noint" %in% long.rand.ind) {
long.rand.ind2 <- long.rand.ind[-which(long.rand.ind == "noint")]
long.rand.ind.form <- paste("-1", paste(long.rand.ind2, collapse = "+"),
sep = "+")
}
Z2.form <- as.formula(paste("~", long.rand.ind.form, sep = ""))
Z2.frame <- model.frame(Z2.form, data = longdat)
tZ2 <- model.matrix(Z2.form, Z2.frame)
Z2 <- t(tZ2)
tZ2.bystudy <- lapply(1:numstudies, function(u) {
tZ2[which(id.long %in% ids.bystudy[[u]]), ]
})
Z2.bystudy <- lapply(1:numstudies, function(u) {
Z2[, which(id.long %in% ids.bystudy[[u]])]
})
s1.2 <- rep(1:(q - 1), (q - 1):1)
s2.2 <- sequence((q - 1):1) + rep(1:(q - 1), (q - 1):1)
Z2.form.surv <- as.formula(gsub(time.long, names(survdat)[2], Z2.form))
survdat$survdatorder <- 1:nrow(survdat)
tempsurv <- survdat[, which(!(names(survdat) %in% names(longdat)))]
tempsurv <- cbind(survdat[, 1], tempsurv)
names(tempsurv)[1] <- names(survdat)[1]
tempdat <- merge(tempsurv, longdat[match(unique(longdat[, which(names(longdat) %in%
id.name)]), longdat[, which(names(longdat) %in% id.name)]), c(1,
5:ncol(longdat))], by = id.name)
tempdat <- tempdat[order(tempdat$survdatorder), ]
tempdat <- tempdat[, -which(names(tempdat) %in% "survdatorder")]
survdat <- survdat[, -which(names(survdat) %in% "survdatorder")]
Z2.frame.surv <- model.frame(Z2.form.surv, data = tempdat)
tZ2.surv <- model.matrix(Z2.form.surv, Z2.frame.surv)
tZ2.surv.bystudy <- lapply(1:numstudies, function(u) {
tZ2.surv[which(ids.surv %in% ids.bystudy[[u]]), ]
})
Z2.event.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[[u]], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
if (q > 1) {
temp <- matrix(rep(tZ2.bystudy[[u]][v, ], rstemp), ncol = q,
byrow = TRUE)
} else {
temp <- matrix(rep(tZ2.bystudy[[u]][v], rstemp), ncol = q,
byrow = TRUE)
}
colnames(temp) <- colnames(tZ2)
if (time.long %in% colnames(temp)) {
if (strat) {
temp[, colnames(temp) == time.long] <- sf[[u]][1:rstemp]
} else {
temp[, colnames(temp) == time.long] <- sf[1:rstemp]
}
}
t(temp)
} else {
matrix(rep(0, q))
}
})
})
Z2.event.sqr.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
temp <- (Z2.event.bystudy[[u]][[v]])^2
if (q > 1) {
if (rstemp == 1) {
temp <- rbind(temp, matrix(Z2.event.bystudy[[u]][[v]][s1.2,
] * Z2.event.bystudy[[u]][[v]][s2.2, ]))
} else {
temp <- rbind(temp, Z2.event.bystudy[[u]][[v]][s1.2,
] * Z2.event.bystudy[[u]][[v]][s2.2, ])
}
}
temp
} else {
matrix(rep(0, sum(1:q)))
}
})
})
Z2b2event <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = q)
} else {
matrix(0, nrow = length(haz), ncol = q)
}
})
})
Z2b2eventsqr <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = sum(1:q))
} else {
matrix(0, nrow = length(haz), ncol = sum(1:q))
}
})
})
newu.2.all <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
idtemp <- ids.bystudy[[u]][v]
B2temp <- paraests$randstart.ind.cov[[which(names(paraests$randstart.ind.cov) %in%
idtemp)]]
cm.2 <- matrix(as.numeric(as.character(paraests$randstart.ind[which(rownames(paraests$randstart.ind) %in%
idtemp), ])), gpt^q, q, TRUE)
gmat2 %*% solve(chol(solve(B2temp))) + cm.2
})
})
randstart.ind <- paraests$randstart.ind
randstart.ind.bystudy <- lapply(1:numstudies, function(u) {
out <- as.matrix(randstart.ind[which(rownames(randstart.ind) %in%
ids.bystudy[[u]]), ])
rownames(out) <- rownames(randstart.ind)[which(rownames(randstart.ind) %in%
ids.bystudy[[u]])]
out
})
randstart.ind.bystudy <- lapply(1:numstudies, function(u) {
randstart.ind.bystudy[[u]][match(ids.bystudy[[u]], rownames(randstart.ind.bystudy[[u]])),
]
})
if (is.null(long.rand.stud) == FALSE) {
A <- as.matrix(paraests$A)
w3 <- g$weights * sqrt(pi)
gmat3 <- matrix(0, gpt^r, r)
gmat3[, 1] <- rep(ab, each = gpt^(r - 1))
if (r > 1) {
gmat3[, 2] <- rep(ab, gpt)
w3 <- as.vector(w3 %x% w3)
if (r > 2) {
gmat3[, 3] <- rep(ab, each = gpt)
w3 <- as.vector(w3 %x% g$weights * sqrt(pi))
}
}
EU.3 <- matrix(0, numstudies, r)
EUU.3 <- matrix(0, numstudies, sum(1:r))
EexpU.3 <- lapply(1:numstudies, function(u) {
if (strat) {
matrix(0, n.bystudy[[u]], length(haz[[u]]))
} else {
matrix(0, n.bystudy[[u]], length(haz))
}
})
if (study.name %in% long.rand.stud) {
long.rand.stud2 <- long.rand.stud
long.rand.stud2[which(long.rand.stud2 == study.name)] <- "1"
long.rand.stud.form <- paste(long.rand.stud2, collapse = "+")
} else {
long.rand.stud2 <- long.rand.stud
long.rand.stud.form <- paste("-1", paste(long.rand.stud2, collapse = "+"),
sep = "+")
}
Z3.form <- as.formula(paste("~", long.rand.stud.form, sep = ""))
Z3.frame <- model.frame(Z3.form, data = longdat)
tZ3 <- model.matrix(Z3.form, Z3.frame)
Z3 <- t(tZ3)
tZ3.bystudy <- lapply(1:numstudies, function(u) {
tZ3[which(id.long %in% ids.bystudy[[u]]), ]
})
Z3.bystudy <- lapply(1:numstudies, function(u) {
Z3[, which(id.long %in% ids.bystudy[[u]])]
})
s1.3 <- rep(1:(r - 1), (r - 1):1)
s2.3 <- sequence((r - 1):1) + rep(1:(r - 1), (r - 1):1)
Z3.form.surv <- as.formula(gsub(time.long, names(survdat)[2], Z3.form))
Z3.frame.surv <- model.frame(Z3.form.surv, data = tempdat)
tempdat <- NULL
tZ3.surv <- model.matrix(Z3.form.surv, Z3.frame.surv)
Z3.surv <- t(tZ3.surv)
tZ3.surv.bystudy <- lapply(1:numstudies, function(u) {
tZ3.surv[which(ids.surv %in% ids.bystudy[[u]]), ]
})
Z3.surv.bystudy <- lapply(1:numstudies, function(u) {
Z3.surv[, which(ids.surv %in% ids.bystudy[[u]])]
})
Z3.event.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[[u]], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
if (r > 1) {
temp <- matrix(rep(tZ3.bystudy[[u]][v, ], rstemp),
ncol = r, byrow = TRUE)
} else {
temp <- matrix(rep(tZ3.bystudy[[u]][v], rstemp), ncol = r,
byrow = TRUE)
}
colnames(temp) <- colnames(tZ3)
if (time.long %in% colnames(temp)) {
if (strat) {
temp[, colnames(temp) == time.long] <- sf[[u]][1:rstemp]
} else {
temp[, colnames(temp) == time.long] <- sf[1:rstemp]
}
}
t(temp)
} else {
matrix(rep(0, r))
}
})
})
Z3.event.sqr.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
temp <- (Z3.event.bystudy[[u]][[v]])^2
if (r > 1) {
if (rstemp == 1) {
temp <- rbind(temp, matrix(Z3.event.bystudy[[u]][[v]][s1.3] *
Z3.event.bystudy[[u]][[v]][s2.3, ]))
} else {
temp <- rbind(temp, Z3.event.bystudy[[u]][[v]][s1.3,
] * Z3.event.bystudy[[u]][[v]][s2.3, ])
}
}
temp
} else {
matrix(rep(0, sum(1:r)))
}
})
})
randselector1 <- rep(1:r, each = q)
randselector2 <- rep(1:q, times = r)
Z3Z2.event.sqr.bystudy <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
rstemp <- rs[[u]][which(names(rs[[u]]) %in% ids.bystudy[[u]][[v]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[u]][[v]])]
}
if (rstemp > 0) {
Z3.event.bystudy[[u]][[v]][randselector1, ] * Z2.event.bystudy[[u]][[v]][randselector2,
]
} else {
matrix(rep(0, length(randselector2)))
}
})
})
Z3b3event <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = r)
} else {
matrix(0, nrow = length(haz), ncol = r)
}
})
})
Z3b3eventsqr <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
if (strat) {
matrix(0, nrow = length(haz[[u]]), ncol = sum(1:r))
} else {
matrix(0, nrow = length(haz), ncol = sum(1:r))
}
})
})
newu.3.all <- lapply(1:numstudies, function(u) {
studtemp <- studies[u]
B3temp <- paraests$randstart.stud.cov[[which(names(paraests$randstart.stud.cov) %in%
studtemp)]]
cm.3 <- matrix(as.numeric(as.character(paraests$randstart.stud[which(rownames(paraests$randstart.stud) %in%
studtemp), ])), gpt^r, r, TRUE)
gmat3 %*% solve(chol(solve(B3temp))) + cm.3
})
names(newu.3.all) <- studies
randstart.stud <- paraests$randstart.stud
survnames <- NULL
if (p2 > 0) {
survnames <- paste("T.", names(beta2)[1:(length(beta2) - 2)],
sep = "")
}
paranames <- c(paste("Y.", names(beta1), sep = ""), survnames,
"gamma_ind_0", "gamma_stud_0", "sigma.e", paste("D", paste(rep(1:q,
times = q), rep(1:q, each = q), sep = ""), sep = ""), paste("A",
paste(rep(1:r, times = r), rep(1:r, each = r), sep = ""),
sep = ""))
} else {
r <- NULL
survnames <- NULL
if (p2 > 0) {
survnames <- paste("T.", names(beta2)[1:(length(beta2) - 1)],
sep = "")
}
paranames <- c(paste("Y.", names(beta1), sep = ""), survnames,
"gamma_ind_0", "sigma.e", paste("D", paste(rep(1:q, times = q),
rep(1:q, each = q), sep = ""), sep = ""))
}
conv <- FALSE
for (it in 1:max.it) {
if (p2 > 0) {
beta2x <- lapply(1:numstudies, function(u) {
X2.bystudy[[u]] %*% beta2[1:p2]
})
} else {
beta2x <- lapply(1:numstudies, function(u) {
matrix(rep(0, n.bystudy[u]))
})
}
ebeta2x <- lapply(1:numstudies, function(u) {
exp(beta2x[[u]])
})
for (k in 1:numstudies) {
if (is.null(long.rand.stud) == FALSE) {
if (it == 1) {
b2est <- randstart.ind.bystudy[[k]]
} else {
b2est <- EU.2[[k]]
}
if (it == 1) {
b3est <- as.numeric(randstart.stud[k, ])
} else {
b3est <- as.numeric(EU.3[k, ])
}
newu.3 <- newu.3.all[[k]]
newu.3sqr <- newu.3^2
if (r > 1) {
newu.3sqr <- cbind(newu.3sqr, newu.3[, s1.3] * newu.3[,
s2.3])
}
egDUs.3 <- t(apply(exp(newu.3 %*% (Z3.surv.bystudy[[k]] *
beta2[(p2 + 2)])), 1, function(u) {
u^cen.bystudy[[k]]
}))
egDUsf.3 <- lapply(1:n.bystudy[k], function(v) {
exp(newu.3 %*% (Z3.event.bystudy[[k]][[v]] * beta2[(p2 +
2)]))
})
if (q > 1) {
egDUsf.2 <- lapply(1:n.bystudy[k], function(v) {
exp(as.numeric(b2est[v, ]) %*% (Z2.event.bystudy[[k]][[v]] *
beta2[(p2 + 1)]))
})
} else {
egDUsf.2 <- lapply(1:n.bystudy[k], function(v) {
exp(as.numeric(b2est[v]) %*% (Z2.event.bystudy[[k]][[v]] *
beta2[(p2 + 1)]))
})
}
ess.3 <- lapply(1:n.bystudy[k], function(v) {
if (strat) {
exp(-(ebeta2x[[k]][v, ] * egDUsf.3[[v]]) %*% ((haz[[k]][1:rs[[k]][which(names(rs[[k]]) %in%
ids.bystudy[[k]][[v]])]] * as.numeric(egDUsf.2[[v]]))))
} else {
exp(-(ebeta2x[[k]][v, ] * egDUsf.3[[v]]) %*% ((haz[1:rs[which(names(rs) %in%
ids.bystudy[[k]][[v]])]] * as.numeric(egDUsf.2[[v]]))))
}
})
f.3 <- rowSums(do.call(cbind, lapply(1:n.bystudy[k], function(v) {
rowSums(egDUs.3[, v] * ess.3[[v]] * w3)
})))
den.3 <- sum(f.3)
EU.3[k, 1:r] <- f.3 %*% newu.3/den.3
EUU.3[k, 1:sum(1:r)] <- f.3 %*% newu.3sqr/den.3
}
for (i in 1:n.bystudy[k]) {
if (strat) {
rstemp <- rs[[k]][which(names(rs[[k]]) %in% ids.bystudy[[k]][[i]])]
} else {
rstemp <- rs[which(names(rs) %in% ids.bystudy[[k]][[i]])]
}
newu.2 <- newu.2.all[[k]][[i]]
newu.2sqr <- newu.2^2
if (q > 1) {
newu.2sqr <- cbind(newu.2sqr, newu.2[, s1.2] * newu.2[,
s2.2])
}
egDUs.2.ind <- 1
egDUs.3.const <- 1
egDUs.3.ind <- 1
egDUs.2.const <- 1
if (cen.bystudy[[k]][i] == 1) {
if (q > 1) {
egDUs.2.ind <- exp(newu.2 %*% (tZ2.surv.bystudy[[k]][i,
] * beta2[(p2 + 1)]))
} else {
egDUs.2.ind <- exp(newu.2 %*% (tZ2.surv.bystudy[[k]][i] *
beta2[(p2 + 1)]))
}
}
egDUsf.2.ind <- exp(newu.2 %*% (Z2.event.bystudy[[k]][[i]] *
beta2[(p2 + 1)]))
if (is.null(long.rand.stud) == FALSE) {
if (cen.bystudy[[k]][i] == 1) {
if (r > 1) {
egDUs.3.const <- exp(b3est %*% (tZ3.surv.bystudy[[k]][i,
] * beta2[(p2 + 2)]))
egDUs.3.ind <- exp(newu.3 %*% (tZ3.surv.bystudy[[k]][i,
] * beta2[(p2 + 2)]))
} else {
egDUs.3.const <- exp(b3est %*% (tZ3.surv.bystudy[[k]][i] *
beta2[(p2 + 2)]))
egDUs.3.ind <- exp(newu.3 %*% (tZ3.surv.bystudy[[k]][i] *
beta2[(p2 + 2)]))
}
if (q > 1) {
egDUs.2.const <- exp(b2est[i, ] %*% (tZ2.surv.bystudy[[k]][i,
] * beta2[(p2 + 1)]))
} else {
egDUs.2.const <- exp(b2est[i] %*% (tZ2.surv.bystudy[[k]][i] *
beta2[(p2 + 1)]))
}
}
if (q > 1) {
egDUsf.2.const <- exp(b2est[i, ] %*% (Z2.event.bystudy[[k]][[i]] *
beta2[(p2 + 1)]))
} else {
egDUsf.2.const <- exp(b2est[i] %*% (Z2.event.bystudy[[k]][[i]] *
beta2[(p2 + 1)]))
}
egDUsf.3.const <- exp(b3est %*% (Z3.event.bystudy[[k]][[i]] *
beta2[(p2 + 2)]))
egDUsf.3.ind <- exp(newu.3 %*% (Z3.event.bystudy[[k]][[i]] *
beta2[(p2 + 2)]))
if (strat) {
ess.3.ind <- exp(-(ebeta2x[[k]][i, ] * egDUsf.3.ind) %*%
t(haz[[k]][1:rstemp] * egDUsf.2.const))
} else {
ess.3.ind <- exp(-(ebeta2x[[k]][i, ] * egDUsf.3.ind) %*%
t(haz[1:rstemp] * egDUsf.2.const))
}
f.3.ind <- egDUs.3.ind * ess.3.ind * w3
den.3.ind <- sum(f.3.ind)
C.3 <- egDUsf.3.ind[, 1:rstemp]
EexpU.3[[k]][i, 1:rstemp] <- f.3.ind[, 1] %*% C.3/den.3.ind
Z3b3event[[k]][[i]][1:ncol(Z3.event.bystudy[[k]][[i]]),
] <- (t(Z3.event.bystudy[[k]][[i]] * as.numeric(f.3.ind[,
1] %*% newu.3)))/den.3.ind
Z3b3eventsqr[[k]][[i]][1:ncol(Z3.event.sqr.bystudy[[k]][[i]]),
] <- (as.numeric(f.3.ind[, 1] %*% newu.3sqr) * Z3.event.sqr.bystudy[[k]][[i]])/den.3.ind
}
if (strat) {
if (is.null(long.rand.stud) == FALSE) {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
t(haz[[k]][1:rstemp] * egDUsf.3.const))
} else {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
(haz[[k]][1:rstemp]))
}
} else {
if (is.null(long.rand.stud) == FALSE) {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
t(haz[1:rstemp] * egDUsf.3.const))
} else {
ess.2 <- exp(-(ebeta2x[[k]][i, ] * egDUsf.2.ind) %*%
haz[1:rstemp])
}
}
f.2 <- egDUs.2.ind * ess.2 * w2
den.2 <- sum(f.2)
EU.2[[k]][i, 1:q] <- f.2[, 1] %*% newu.2/den.2
EUU.2[[k]][i, 1:sum(1:q)] <- f.2[, 1] %*% newu.2sqr/den.2
C.2 <- egDUsf.2.ind[, 1:rstemp]
EexpU.2[[k]][i, 1:rstemp] <- f.2[, 1] %*% C.2/den.2
Z2b2event[[k]][[i]][1:ncol(Z2.event.bystudy[[k]][[i]]),
] <- (t(Z2.event.bystudy[[k]][[i]] * as.numeric(f.2[,
1] %*% newu.2)))/den.2
Z2b2eventsqr[[k]][[i]][1:ncol(Z2.event.sqr.bystudy[[k]][[i]]),
] <- (as.numeric(f.2[, 1] %*% newu.2sqr) * Z2.event.sqr.bystudy[[k]][[i]])/den.2
}
}
parac <- data.frame(c(beta1, beta2, sigma.e, D))
if (is.null(long.rand.stud) == FALSE) {
parac <- data.frame(c(beta1, beta2, sigma.e, D, A))
}
EexpU <- lapply(1:numstudies, function(u) {
if (is.null(long.rand.stud) == FALSE) {
EexpU.2[[u]] * EexpU.3[[u]]
} else {
EexpU.2[[u]]
}
})
if (strat == FALSE) {
haz <- nev/colSums(do.call(rbind, EexpU) * do.call(rbind, ebeta2x)[,
1])
} else {
haz <- lapply(1:numstudies, function(u) {
nev[[u]]/colSums(EexpU[[u]] * ebeta2x[[u]][, 1])
})
}
EUmat.2 <- lapply(1:numstudies, function(u) {
apply(EU.2[[u]], 2, rep, nn[[u]])
})
EUUmat.2 <- lapply(1:numstudies, function(u) {
apply(EUU.2[[u]], 2, rep, nn[[u]])
})
Ut.2 <- lapply(1:numstudies, function(u) {
rowSums(EUmat.2[[u]] * tZ2.bystudy[[u]])
})
if (is.null(long.rand.stud) == FALSE) {
EUmat.3 <- lapply(1:numstudies, function(u) {
matrix(rep(EU.3[u, ], sum(nn[[u]])), ncol = r)
})
EUUmat.3 <- lapply(1:numstudies, function(u) {
matrix(rep(EUU.3[u, ], sum(nn[[u]])), ncol = sum(1:r))
})
Ut.3 <- lapply(1:numstudies, function(u) {
rowSums(EUmat.3[[u]] * tZ3.bystudy[[u]])
})
beta1 <- as.numeric(solve(crossprod(do.call(rbind, X1.bystudy)),
crossprod(do.call(rbind, X1.bystudy), do.call(c, Y.bystudy) -
(do.call(c, Ut.2) + do.call(c, Ut.3)))))
if (q == 1) {
if (r == 1) {
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (do.call(c, Z2.bystudy) *
do.call(c, lapply(1:numstudies, function(u) {
rep(EU.2[[u]], nn[[u]])
}))) + (do.call(c, Z3.bystudy) * do.call(c, lapply(1:numstudies,
function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]])]
})))))^2)/sum(N)
} else {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z3.bystudy)
out2 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (do.call(c, Z2.bystudy) *
do.call(c, lapply(1:numstudies, function(u) {
rep(EU.2[[u]], nn[[u]])
}))) + out2))^2)/sum(N)
}
} else {
if (r == 1) {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.2[[u]][rep(1:nrow(EU.2[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z2.bystudy)
out1 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + out1 + (do.call(c, Z3.bystudy) *
do.call(c, lapply(1:numstudies, function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]])]
})))))^2)/sum(N)
} else {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.2[[u]][rep(1:nrow(EU.2[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z2.bystudy)
out1 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.3[rep(rep(u, n.bystudy[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z3.bystudy)
out2 <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + out1 + out2))^2)/sum(N)
}
}
diag(A) <- colMeans(EUU.3)[1:r]
if (r > 1) {
A[lower.tri(A)] <- colMeans(EUU.3)[-(1:r)]
A[upper.tri(A)] <- t(A)[upper.tri(A)]
}
} else {
beta1 <- as.numeric(solve(crossprod(do.call(rbind, X1.bystudy)),
crossprod(do.call(rbind, X1.bystudy), do.call(c, Y.bystudy) -
(do.call(c, Ut.2)))))
if (q == 1) {
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (do.call(c, Z2.bystudy) * do.call(c,
lapply(1:numstudies, function(u) {
rep(EU.2[[u]], nn[[u]])
})))))^2)/sum(N)
} else {
temp1 <- do.call(rbind, lapply(1:numstudies, function(u) {
EU.2[[u]][rep(1:nrow(EU.2[[u]]), nn[[u]]), ]
}))
temp2 <- do.call(cbind, Z2.bystudy)
out <- unlist(lapply(1:nrow(temp1), function(u) {
sum(temp1[u, ] * temp2[, u])
}))
sigma.e <- colSums((do.call(c, Y.bystudy) - ((do.call(rbind,
X1.bystudy) %*% beta1) + (out)))^2)/sum(N)
}
}
diag(D) <- colMeans(do.call(rbind, EUU.2))[1:q]
if (q > 1) {
D[lower.tri(D)] <- colMeans(do.call(rbind, EUU.2))[-(1:q)]
D[upper.tri(D)] <- t(D)[upper.tri(D)]
}
if (is.null(long.rand.stud) == FALSE) {
EexpU <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
EexpU.2[[u]][v, ] * EexpU.3[[u]][v, ]
})
})
fd <- vector("numeric", p2 + q + r)
sd <- matrix(0, p2 + q + r, p2 + q + r)
} else {
EexpU <- lapply(1:numstudies, function(u) {
lapply(1:n.bystudy[u], function(v) {
EexpU.2[[u]][v, ]
})
})
fd <- vector("numeric", p2 + q)
sd <- matrix(0, p2 + q, p2 + q)
}
if (q > 1) {
fd[(p2 + 1):(p2 + q)] <- colSums(do.call(c, cen.bystudy) *
(do.call(rbind, EU.2) * do.call(rbind, tZ2.surv.bystudy))) -
colSums(do.call(rbind, ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
fd[(p2 + 1):(p2 + q)] <- sum(do.call(c, cen.bystudy) * (do.call(rbind,
EU.2) * do.call(c, tZ2.surv.bystudy))) - colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
})))
}
if (is.null(long.rand.stud) == FALSE) {
if (r > 1) {
fd[(p2 + q + 1):(p2 + q + r)] <- colSums(do.call(c, cen.bystudy) *
(EU.3[rep(1:nrow(EU.3), n.bystudy), ]) * do.call(rbind,
tZ3.surv.bystudy)) - colSums(do.call(rbind, ebeta2x)[,
1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
fd[(p2 + q + 1):(p2 + q + r)] <- sum(do.call(c, cen.bystudy) *
(EU.3[rep(1:nrow(EU.3), n.bystudy), ]) * do.call(c, tZ3.surv.bystudy)) -
colSums(do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]] * EexpU[[u]][[v]] *
haztemp)
}))
})))
}
}
if (q > 1) {
if (q == 2) {
sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)][upper.tri(sd[(p2 +
1):(p2 + q), (p2 + 1):(p2 + q)])] <- -sum(do.call(rbind,
ebeta2x)[, 1] * 0.5 * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z2b2eventsqr[[u]][[v]][, (q + 1):sum(1:q)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)][upper.tri(sd[(p2 +
1):(p2 + q), (p2 + 1):(p2 + q)])] <- -colSums(do.call(rbind,
ebeta2x)[, 1] * 0.5 * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z2b2eventsqr[[u]][[v]][, (q + 1):sum(1:q)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
}
}
if (is.null(long.rand.stud) == FALSE) {
if (r > 1) {
if (r == 2) {
sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 + q +
r)][upper.tri(sd[(p2 + q + 1):(p2 + q + r), (p2 + q +
1):(p2 + q + r)])] <- -sum(do.call(rbind, ebeta2x)[,
1] * 0.5 * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z3b3eventsqr[[u]][[v]][, (r + 1):sum(1:r)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
} else {
sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 + q +
r)][upper.tri(sd[(p2 + q + 1):(p2 + q + r), (p2 + q +
1):(p2 + q + r)])] <- -colSums(do.call(rbind, ebeta2x)[,
1] * 0.5 * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z3b3eventsqr[[u]][[v]][, (r + 1):sum(1:r)] *
EexpU[[u]][[v]] * haztemp)
}))
})))
}
}
if (q == 1 && r == 1) {
sd[(p2 + 1):(p2 + q), ((p2 + q + 1):(p2 + q + r))] <- -colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(Z3b3event[[u]][[v]] * Z2b2event[[u]][[v]] * EexpU.2[[u]][v,
] * EexpU.3[[u]][v, ] * haztemp)
}))
})))
} else {
sd[(p2 + 1):(p2 + q), ((p2 + q + 1):(p2 + q + r))] <- -colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]][, randselector1] * Z2b2event[[u]][[v]][,
randselector2] * EexpU[[u]][[v]] * haztemp)
}))
})))
}
}
if (p2 > 0) {
fd[1:p2] <- c(colSums((do.call(c, cen.bystudy) * do.call(rbind,
X2.bystudy)) - (do.call(rbind, X2.bystudy) * as.numeric(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(EexpU[[u]][[v]] * haztemp)
}))
}))))))
sd[(1:p2), (p2 + 1):(p2 + q)] <- -(t(do.call(rbind, X2.bystudy)) %*%
(do.call(rbind, ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies,
function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))
if (is.null(long.rand.stud) == FALSE) {
sd[(1:p2), (p2 + q + 1):(p2 + q + r)] <- -(t(do.call(rbind,
X2.bystudy)) %*% (do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3event[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))
}
sd <- sd + t(sd)
for (i in 1:p2) {
for (j in 1:p2) {
sd[i, j] <- -(sum(do.call(rbind, X2.bystudy)[, i] * do.call(rbind,
X2.bystudy)[, j] * do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
sum(EexpU[[u]][[v]] * haztemp)
}))
}))))
}
}
}
if (q == 1) {
sd[(p2 + 1), (p2 + 1)] <- -(colSums(do.call(rbind, ebeta2x)[,
1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2eventsqr[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))[1:q]
} else {
diag(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)]) <- -(colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z2b2eventsqr[[u]][[v]] * EexpU[[u]][[v]] * haztemp)
}))
}))))[1:q]
}
if (is.null(long.rand.stud) == FALSE) {
if (r == 1) {
sd[(p2 + q + 1), (p2 + q + 1)] <- -(colSums(do.call(rbind,
ebeta2x)[, 1] * do.call(rbind, lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3eventsqr[[u]][[v]] * EexpU[[u]][[v]] *
haztemp)
}))
}))))[1:r]
} else {
diag(sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 + q +
r)]) <- -(colSums(do.call(rbind, ebeta2x)[, 1] * do.call(rbind,
lapply(1:numstudies, function(u) {
if (strat) {
haztemp <- haz[[u]]
} else {
haztemp <- haz
}
do.call(rbind, lapply(1:n.bystudy[u], function(v) {
colSums(Z3b3eventsqr[[u]][[v]] * EexpU[[u]][[v]] *
haztemp)
}))
}))))[1:r]
}
}
if (is.null(long.rand.stud) == FALSE) {
if (r == 1 && q == 1) {
fd <- fd
sd <- sd
} else {
if (p2 > 0) {
fd <- c(fd[1:p2], sum(fd[(p2 + 1):(p2 + q)]), sum(fd[(p2 +
q + 1):(p2 + q + r)]))
sd.p2 <- sd[1:p2, 1:p2]
sd.q <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
sd.r <- sum(sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 +
q + r)])
if (p2 > 1) {
if (q > 1) {
sd.p2q <- rowSums(sd[1:p2, (p2 + 1):(p2 + q)])
} else {
sd.p2q <- sd[1:p2, (p2 + 1):(p2 + q)]
}
if (r > 1) {
sd.p2r <- rowSums(sd[1:p2, (p2 + q + 1):(p2 + q +
r)])
} else {
sd.p2r <- sd[1:p2, (p2 + q + 1):(p2 + q + r)]
}
} else {
sd.p2q <- sum(sd[1:p2, (p2 + 1):(p2 + q)])
sd.p2r <- sum(sd[1:p2, (p2 + q + 1):(p2 + q + r)])
}
sd.qr <- sum(sd[(p2 + 1):(p2 + q), (p2 + q + 1):(p2 +
q + r)])
sd <- rbind(cbind(sd.p2, sd.p2q, sd.p2r), c(sd.p2q, sd.q,
sd.qr), c(sd.p2r, sd.qr, sd.r))
rownames(sd) <- colnames(sd) <- NULL
} else {
fd <- c(sum(fd[(p2 + 1):(p2 + q)]), sum(fd[(p2 + q +
1):(p2 + q + r)]))
sd.q <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
sd.r <- sum(sd[(p2 + q + 1):(p2 + q + r), (p2 + q + 1):(p2 +
q + r)])
sd.qr <- sum(sd[(p2 + 1):(p2 + q), (p2 + q + 1):(p2 +
q + r)])
sd <- rbind(c(sd.q, sd.qr), c(sd.qr, sd.r))
rownames(sd) <- colnames(sd) <- NULL
}
}
} else {
if (q == 1) {
fd <- fd
sd <- sd
} else {
if (p2 > 0) {
fd <- c(fd[1:p2], sum(fd[(p2 + 1):(p2 + q)]))
sd.p2 <- sd[1:p2, 1:p2]
sd.q <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
if (p2 > 1) {
sd.p2q <- rowSums(sd[1:p2, (p2 + 1):(p2 + q)])
} else {
sd.p2q <- sum(sd[1:p2, (p2 + 1):(p2 + q)])
}
sd <- rbind(cbind(sd.p2, sd.p2q), c(sd.p2q, sd.q))
rownames(sd) <- colnames(sd) <- NULL
} else {
fd <- c(sum(fd[(p2 + 1):(p2 + q)]))
sd <- sum(sd[(p2 + 1):(p2 + q), (p2 + 1):(p2 + q)])
rownames(sd) <- colnames(sd) <- NULL
}
}
}
beta2 <- beta2 - solve(sd, fd)
para <- data.frame(c(beta1, beta2, sigma.e, D))
if (is.null(long.rand.stud) == FALSE) {
para <- data.frame(c(beta1, beta2, sigma.e, D, A))
}
dd <- abs(parac - para)
if (print.detail) {
print(paste("Iteration: ", it, sep = ""))
print("Current parameter estimates:")
detail <- data.frame(para)
rownames(detail) <- paranames
colnames(detail) <- NULL
print(detail)
}
if (max(dd) < tol) {
conv <- TRUE
break
}
}
if (conv != TRUE && bootrun == FALSE) {
print("Not converged")
}
if (strat) {
names(haz) <- studies
}
results <- list(beta1 = data.frame(beta1), beta2 = data.frame(beta2),
sigma.e = sigma.e, D = D, haz = haz, random2 = EU.2, conv = conv,
iters = it, long.rand.ind.form = long.rand.ind.form, n.bystudy = n.bystudy)
if (is.null(long.rand.stud) == FALSE) {
results$A <- A
results$random3 <- EU.3
results$long.rand.stud.form <- long.rand.stud.form
}
results
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.