# R/hessian.R In joineRML: Joint Modelling of Multivariate Longitudinal Data and Time-to-Event Outcomes

#### Defines functions hessian

```#' @keywords internal
hessian <- function(theta, l, t, z, m) {

# MLE parameter estimates from EM algorithm
D <- theta\$D
beta <- theta\$beta
sigma2 <- theta\$sigma2
haz <- theta\$haz
gamma <- theta\$gamma

# Multivariate longitudinal data
yi <- l\$yi
Xi <- l\$Xi
Zi <- l\$Zi
Zit <- l\$Zit
nik <- l\$nik
yik <- l\$yik
Xik.list <- l\$Xik.list
Zik.list <- l\$Zik.list
p <- l\$p    # vector of fixed effect dims
r <- l\$r    # vector of random effect dims
K <- l\$K    # number of longitudinal markers
n <- l\$n    # number of subjects
nk <- l\$nk  # vector of number of observations per outcome

# Time-to-event data
V <- t\$V
survdat2 <- t\$survdat2
survdat2.list <- t\$survdat2.list
q <- t\$q
nev <- t\$nev
nev.uniq <- t\$nev.uniq

# Covariate data for W(u, b)
Zi.fail <- z\$Zi.fail
Zit.fail <- z\$Zit.fail
IW.fail <- z\$IW.fail

# MCEM calculations
Sigmai.inv <- m\$Sigmai.inv
Eb <- m\$Eb
EbbT <- m\$EbbT
bi.y <- m\$bi.y
expvstargam <- m\$expvstargam
pb.yt <- m\$pb.yt
haz.hat <- m\$haz.hat

#*****************************************************
# Complete data score components
#*****************************************************

# beta

sbeta <- mapply(function(x, sinv, y, z, b) {
(t(x) %*% sinv) %*% (y - x %*% beta - z %*% b)
},
x = Xi, sinv = Sigmai.inv, y = yi, z = Zi, b = Eb,
SIMPLIFY = TRUE)

rownames(sbeta) <- names(beta)

#-----------------------------------------------------

# D

Dinv <- solve(D)

D.inds <- which(lower.tri(D, diag = TRUE), arr.ind = TRUE)
dimnames(D.inds) <- NULL

delta.D <- lapply(1:nrow(D.inds), function(x, ind) {
mat <- matrix(0, nrow = nrow(D), ncol = ncol(D))
ii <- ind[x, , drop = FALSE]
mat[ii[1], ii[2]] <- mat[ii[2], ii[1]] <- 1
mat
}, ind = D.inds[, 2:1, drop = FALSE])

term1 <- sapply(delta.D, function(d) {
-0.5 * sum(diag(Dinv %*% d))
})

sDi <- function(i) {
mapply(function(b, pb) {
out <- 0.5 * crossprod(b, b * pb) %*% (Dinv %*% delta.D[[i]] %*% Dinv) / nrow(b)
term1[i] + sum(diag(out))
},
b = bi.y, pb = pb.yt,
SIMPLIFY = TRUE)
}

sD <- sapply(1:nrow(D.inds), sDi)
sD <- t(sD)
rownames(sD) <- paste0("D", D.inds[, 1], ",", D.inds[, 2])

#-----------------------------------------------------

# gamma

sgamma <- gammaUpdate_approx(bi.y, Zit.fail, expvstargam, pb.yt, haz.hat,
V, survdat2.list, K, q, nev.uniq)\$scorei

rownames(sgamma) <- names(gamma)

#-----------------------------------------------------

# sigma2

beta.inds <- cumsum(c(0, p))
b.inds <- cumsum(c(0, r))

ssigma2 <- matrix(nrow = K, ncol = n)
for (k in 1:K) {
beta.k <- beta[(beta.inds[k] + 1):(beta.inds[k + 1])]
ssigma2[k, ] <- mapply(function(y, x, z, b, b2, nik) {
b.k <- b[(b.inds[k] + 1):(b.inds[k + 1])]
bbT.k <- b2[(b.inds[k] + 1):(b.inds[k + 1]), (b.inds[k] + 1):(b.inds[k + 1])]
residFixed <- (y - x %*% beta.k)
resids <- t(residFixed) %*% (residFixed - 2*(z %*% b.k)) + sum(diag(crossprod(z) %*% bbT.k))
(-0.5 * nik[k] / sigma2[k]) + (0.5 * resids / sigma2[k]^2)
},
y = yik[[k]], x = Xik.list[[k]], z = Zik.list[[k]], b = Eb, b2 = EbbT, nik = nik)
}

rownames(ssigma2) <- paste0("sigma2_", 1:K)

#-----------------------------------------------------

si <- rbind(sD, sbeta, ssigma2, sgamma)

H <- matrix(0, nrow(si), nrow(si))
for (j in 1:ncol(si)) {
H <- H + tcrossprod(si[, j])
}
# Although RHS term = 0 in theory, in practice with MC integration
# not all terms are vanishingly small, so we add it in
H <- H - (rowSums(si) %*% t(rowSums(si))) / ncol(si)
rownames(H) <- colnames(H) <- rownames(si)

return(H)

}
```

## Try the joineRML package in your browser

Any scripts or data that you put into this service are public.

joineRML documentation built on Jan. 22, 2023, 1:18 a.m.