Nothing
frailty.fit<- function (formula, data, dist.frail = "gamma", dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL)
{
if (!inherits(formula, "formula")) {
if (inherits(formula, "data.frame"))
warning("You gave a data.frame instead of a formula.")
stop("formula is not an object of type formula")
}
if (!inherits(data, "data.frame")) {
if (inherits(data, "formula"))
warning("You gave a formula instead of a data.frame.")
stop("data is not an object of type data.frame")
}
if (missing(formula))
stop("Missing formula")
if (missing(data))
stop("Missing data")
cluster <- function(x) x
mf <- model.frame(formula, data)
pos_cluster <- grep("cluster", names(mf))
if (length(pos_cluster) != 1)
stop("misspecified or non-specified cluster")
cluster <- mf[[pos_cluster]]
Y <- mf[[1]]
if (!inherits(Y, "Surv"))
stop("left hand side not a survival object")
X1 <- model.matrix(formula, data)
pos_cluster_X1 <- grep("cluster", colnames(X1))
x <- X1[, -c(1, pos_cluster_X1), drop = FALSE]
t <- Y[, 1]
delta <- Y[, 2]
if (is.null(dist))
stop("dist must be specified")
if (is.null(prec))
stop("prec must be specified")
if (is.null(max.iter))
stop("max.iter must be specified")
t <- c(t)
delta <- c(delta)
cluster <- c(cluster)
ind = cluster
max.iter <- round(max.iter)
if (!any(dist.frail == c("gamma", "GA", "IG", "WL", "BS",
"TN", "MIG","MBS")))
stop("frailty distribution is not recognized")
if (dist.frail == "gamma")
dist.frail = "GA"
if (!any(dist == c("np", "weibull", "pe", "exponential")))
stop("distribution for time-to-event is not recognized")
if (length(t) != length(delta))
stop("t and delta don't have the same length")
if (length(t) != length(cluster))
stop("t and cluster don't have the same length")
if (prec > 1)
stop("prec is too high")
if (max.iter <= 0)
stop("max.iter at least 1")
H.base <- function(t, lambda, rho = 1, dist, part = NULL) {
if (dist == "weibull") {
Lambda0 <- lambda * t^rho
}
if (dist == "pe") {
Lambda0 <- -ppexp(t, rate = lambda, t = part, lower.tail = FALSE,
log.p = TRUE)
}
Lambda0
}
h.base <- function(t, lambda, rho = 1, dist, log = TRUE,
part = NULL) {
if (dist == "weibull") {
log.lambda0 <- log(lambda) + log(rho) + (rho - 1) *
log(t)
}
if (dist == "pe") {
log.lambda0 <- dpexp(t, rate = lambda, t = part,
log = TRUE) - ppexp(t, rate = lambda, t = part,
lower.tail = FALSE, log.p = TRUE)
}
log.lambda0
}
frailtyGA <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
prof.GA <- function(theta, z, logz) {
ll = -(1/theta) * log(theta) - lgamma(1/theta) +
(1/theta - 1) * logz - z/theta
-sum(ll)
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
tau.GA <- function(theta) {
theta/(theta + 2)
}
if (dist == "weibull") {
observed.llike.0.ga.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
rho = eta[1]
lambda = eta[2]
theta = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
a = r + 1/theta
b = fail.cluster(Lambda0, ind) + 1/theta
P1 = -lgamma(1/theta) - (1/theta) * log(theta) +
lgamma(a) - a * log(b)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ga.dist <- function(eta, t, delta,
x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
theta = eta[ncol(x) + 3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
a = r + 1/theta
b = fail.cluster(Lambda0 * exp(x %*% beta), ind) +
1/theta
P1 = -lgamma(1/theta) - (1/theta) * log(theta) +
lgamma(a) - a * log(b)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ga.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.ga.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
logz.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
a.aux = r + 1/theta.last
b.aux = fail.cluster(Lambda0.last, cluster) +
1/theta.last
z.new = a.aux/b.aux
logz.new = digamma(a.aux) - log(b.aux)
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.ga.dist.Q1.0, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
theta.new = optim(theta.last, prof.GA, z = z.new,
logz = logz.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(rho.last, lambda.last, theta.last) -
c(rho.new, lambda.new, theta.new)))
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ga.dist, x0 = c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.ga.dist(c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster, dist = dist)
para = c(rho.new, lambda.new, theta.new)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
a.aux = r + 1/theta.last
b.aux = fail.cluster(Lambda0.last * exp(x %*%
beta.last), cluster) + 1/theta.last
z.new = a.aux/b.aux
logz.new = digamma(a.aux) - log(b.aux)
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.ga.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
theta.new = optim(theta.last, prof.GA, z = z.new,
logz = logz.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(beta.last, rho.last, lambda.last,
theta.last) - c(beta.new, rho.new, lambda.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.ga.dist, x0 = c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.ga.dist(c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
para = c(beta.last, rho.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "GA"
object.out$tau <- tau.GA(theta.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.ga.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
lambda = eta[1:length(part)]
theta = eta[length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
a = r + 1/theta
b = fail.cluster(Lambda0, ind) + 1/theta
P1 = -lgamma(1/theta) - (1/theta) * log(theta) +
lgamma(a) - a * log(b)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ga.dist <- function(eta, t, delta,
x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
theta = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
a = r + 1/theta
b = fail.cluster(Lambda0 * exp(x %*% beta), ind) +
1/theta
P1 = -lgamma(1/theta) - (1/theta) * log(theta) +
lgamma(a) - a * log(b)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ga.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.ga.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
logz.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
a.aux = r + 1/theta.last
b.aux = fail.cluster(Lambda0.last, cluster) +
1/theta.last
z.new = a.aux/b.aux
logz.new = digamma(a.aux) - log(b.aux)
aux.1 = optim(log(lambda.last), observed.llike.ga.dist.Q1.0,
method = "BFGS", theta = theta.last, z = z.new,
t = t, delta = delta, x = x, ind = ind, dist = dist,
part = part)
lambda.new = exp(aux.1$par)
theta.new = optim(theta.last, prof.GA, z = z.new,
logz = logz.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(lambda.last, theta.last) -
c(lambda.new, theta.new)))
theta.last = theta.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ga.dist, x0 = c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.ga.dist(c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
para = c(lambda.new, theta.new)
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
a.aux = r + 1/theta.last
b.aux = fail.cluster(Lambda0.last * exp(x %*%
beta.last), cluster) + 1/theta.last
z.new = a.aux/b.aux
logz.new = digamma(a.aux) - log(b.aux)
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.ga.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part = part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
theta.new = optim(theta.last, prof.GA, z = z.new,
logz = logz.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(beta.last, lambda.last, theta.last) -
c(beta.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.ga.dist, x0 = c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.ga.dist(c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
para = c(beta.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "GA"
object.out$tau <- tau.GA(theta.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0.ga <- function(eta, t, delta, ind,
cox.aux) {
theta = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
a = r + 1/theta
b = fail.cluster(Lambda0, ind) + 1/theta
P1 = -lgamma(1/theta) - (1/theta) * log(theta) +
lgamma(a) - a * log(b)
-sum(P1)
}
observed.llike.ga <- function(eta, t, delta, x, ind,
cox.aux) {
theta = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
a = r + 1/theta
b = fail.cluster(Lambda0 * exp(x %*% beta), ind) +
1/theta
P1 = -lgamma(1/theta) - (1/theta) * log(theta) +
lgamma(a) - a * log(b)
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
logz.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(z.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
a.aux = r + 1/theta.last
b.aux = fail.cluster(Lambda0.new, cluster) +
1/theta.last
z.new = a.aux/b.aux
logz.new = digamma(a.aux) - log(b.aux)
theta.new = optim(theta.last, prof.GA, z = z.new,
logz = logz.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(theta.last - theta.new))
theta.last = theta.new
z.last = z.new
logz.last = logz.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ga, x0 = c(theta.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(theta.new)
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(z.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
a.aux = r + 1/theta.last
b.aux = fail.cluster(Lambda0.new * exp(x %*%
beta.new), cluster) + 1/theta.last
z.new = a.aux/b.aux
logz.new = digamma(a.aux) - log(b.aux)
theta.new = optim(theta.last, prof.GA, z = z.new,
logz = logz.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(beta.last, theta.last) - c(beta.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.ga, x0 = c(beta.last,
theta.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(beta.new, theta.new)
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "GA"
object.out$tau <- tau.GA(theta.last)
}
object.out
}
frailtyWL <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
prof.WL <- function(theta, z, kappa) {
phi = 4/(theta * (theta + 4))
alpha = sqrt(phi * (phi + 1))
ll = (phi + 1) * log(alpha) - log(alpha + phi) -
lgamma(phi) + (phi - 1) * kappa - alpha * z
-sum(ll)
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
tau.WL <- function(theta) {
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
aux.int <- function(x, theta) {
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
x * (1 + a * x)^(-2 * b - 4) * (1 + theta * x/2) *
(1 + theta * (x + 1)/(theta + 2))
}
4 * a * (1 + b) * integrate(aux.int, lower = 0, upper = Inf,
theta = theta)$value - 1
}
if (dist == "weibull") {
observed.llike.0.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
rho = eta[1]
lambda = eta[2]
theta = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
a.i = 1/(fail.cluster(Lambda0, ind) + 1/a)
b.i = r + b
P1 = log(theta) - log(2) - (b + 1) * log(a) -
lgamma(b) + lgamma(b.i) + b.i * log(a.i) +
log1p(a.i * b.i)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.dist <- function(eta, t, delta, x,
ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
theta = eta[ncol(x) + 3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
a.i = 1/(fail.cluster(Lambda0 * exp(x %*% beta),
ind) + 1/a)
b.i = r + b
P1 = log(theta) - log(2) - (b + 1) * log(a) -
lgamma(b) + lgamma(b.i) + b.i * log(a.i) +
log1p(a.i * b.i)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.wl.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.wl.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = lambda * t^rho
log.lambda0 = log(lambda) + log(rho) + (rho -
1) * log(t)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
kappa.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
phi.aux = 4/(theta.last * (theta.last + 4)) +
r
alpha.aux = 2 * (theta.last + 2)/(theta.last *
(theta.last + 4)) + fail.cluster(Lambda0.last,
cluster)
z.new = phi.aux * (alpha.aux + phi.aux + 1)/(alpha.aux *
(alpha.aux + phi.aux))
kappa.new = -alpha.aux/(phi.aux * (alpha.aux +
phi.aux)) + digamma(phi.aux + 1) - log(alpha.aux)
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.wl.dist.Q1.0, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
theta.new = optim(theta.last, prof.WL, z = z.new,
kappa = kappa.new, method = "Brent", lower = 1e-07,
upper = 20)$par
dif = max(abs(c(rho.last, lambda.last, theta.last) -
c(rho.new, lambda.new, theta.new)))
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
kappa.last = kappa.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.dist, x0 = c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.dist(c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster, dist = dist)
para = c(rho.new, lambda.new, theta.new)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
phi.aux = 4/(theta.last * (theta.last + 4)) +
r
alpha.aux = 2 * (theta.last + 2)/(theta.last *
(theta.last + 4)) + fail.cluster(Lambda0.last *
exp(x %*% beta.last), cluster)
z.new = phi.aux * (alpha.aux + phi.aux + 1)/(alpha.aux *
(alpha.aux + phi.aux))
kappa.new = -alpha.aux/(phi.aux * (alpha.aux +
phi.aux)) + digamma(phi.aux + 1) - log(alpha.aux)
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.wl.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
theta.new = optim(theta.last, prof.WL, z = z.new,
kappa = kappa.new, method = "Brent", lower = 1e-07,
upper = 20)$par
dif = max(abs(c(beta.last, rho.last, lambda.last,
theta.last) - c(beta.new, rho.new, lambda.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
kappa.last = kappa.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.dist, x0 = c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.dist(c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
para = c(beta.last, rho.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "WL"
object.out$tau <- tau.WL(theta.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
lambda = eta[1:length(part)]
theta = eta[length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
a.i = 1/(fail.cluster(Lambda0, ind) + 1/a)
b.i = r + b
P1 = log(theta) - log(2) - (b + 1) * log(a) -
lgamma(b) + lgamma(b.i) + b.i * log(a.i) +
log1p(a.i * b.i)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.dist <- function(eta, t, delta, x,
ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
theta = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
a.i = 1/(fail.cluster(Lambda0 * exp(x %*% beta),
ind) + 1/a)
b.i = r + b
P1 = log(theta) - log(2) - (b + 1) * log(a) -
lgamma(b) + lgamma(b.i) + b.i * log(a.i) +
log1p(a.i * b.i)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.wl.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.wl.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
kappa.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
phi.aux = 4/(theta.last * (theta.last + 4)) +
r
alpha.aux = 2 * (theta.last + 2)/(theta.last *
(theta.last + 4)) + fail.cluster(Lambda0.last,
cluster)
z.new = phi.aux * (alpha.aux + phi.aux + 1)/(alpha.aux *
(alpha.aux + phi.aux))
kappa.new = -alpha.aux/(phi.aux * (alpha.aux +
phi.aux)) + digamma(phi.aux + 1) - log(alpha.aux)
aux.1 = optim(log(lambda.last), observed.llike.wl.dist.Q1.0,
method = "BFGS", theta = theta.last, z = z.new,
t = t, delta = delta, x = x, ind = ind, dist = dist,
part = part)
lambda.new = exp(aux.1$par[1:length(part)])
theta.new = optim(theta.last, prof.WL, z = z.new,
kappa = kappa.new, method = "Brent", lower = 1e-07,
upper = 20)$par
dif = max(abs(c(lambda.last, theta.last) -
c(lambda.new, theta.new)))
theta.last = theta.new
lambda.last = lambda.new
kappa.last = kappa.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.dist, x0 = c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.dist(c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
para = c(lambda.new, theta.new)
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(exp(-coef(cox.aux)[1]/cox.aux$scale),
length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
phi.aux = 4/(theta.last * (theta.last + 4)) +
r
alpha.aux = 2 * (theta.last + 2)/(theta.last *
(theta.last + 4)) + fail.cluster(Lambda0.last *
exp(x %*% beta.last), cluster)
z.new = phi.aux * (alpha.aux + phi.aux + 1)/(alpha.aux *
(alpha.aux + phi.aux))
kappa.new = -alpha.aux/(phi.aux * (alpha.aux +
phi.aux)) + digamma(phi.aux + 1) - log(alpha.aux)
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.wl.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part = part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
theta.new = optim(theta.last, prof.WL, z = z.new,
kappa = kappa.new, method = "Brent", lower = 1e-07,
upper = 20)$par
dif = max(abs(c(beta.last, lambda.last, theta.last) -
c(beta.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
lambda.last = lambda.new
kappa.last = kappa.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.dist, x0 = c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.dist(c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
para = c(beta.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "WL"
object.out$tau <- tau.WL(theta.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0 <- function(eta, t, delta, ind,
cox.aux) {
theta = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
a.i = 1/(fail.cluster(Lambda0, ind) + 1/a)
b.i = r + b
P1 = log(theta) - log(2) - (b + 1) * log(a) -
lgamma(b) + lgamma(b.i) + b.i * log(a.i) +
log1p(a.i * b.i)
-sum(P1)
}
observed.llike <- function(eta, t, delta, x, ind,
cox.aux) {
theta = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
a = theta * (theta + 4)/(2 * (theta + 2))
b = 4/(theta * (theta + 4))
a.i = 1/(fail.cluster(Lambda0 * exp(x %*% beta),
ind) + 1/a)
b.i = r + b
P1 = log(theta) - log(2) - (b + 1) * log(a) -
lgamma(b) + lgamma(b.i) + b.i * log(a.i) +
log1p(a.i * b.i)
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
kappa.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(z.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
phi.aux = 4/(theta.last * (theta.last + 4)) +
r
alpha.aux = 2 * (theta.last + 2)/(theta.last *
(theta.last + 4)) + fail.cluster(Lambda0.new,
cluster)
z.new = phi.aux * (alpha.aux + phi.aux + 1)/(alpha.aux *
(alpha.aux + phi.aux))
kappa.new = -alpha.aux/(phi.aux * (alpha.aux +
phi.aux)) + digamma(phi.aux + 1) - log(alpha.aux)
theta.new = optim(theta.last, prof.WL, z = z.new,
kappa = kappa.new, method = "Brent", lower = 1e-07,
upper = 20)$par
dif = max(abs(theta.last - theta.new))
theta.last = theta.new
kappa.last = kappa.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0, x0 = c(theta.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(theta.new)
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(z.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
phi.aux = 4/(theta.last * (theta.last + 4)) +
r
alpha.aux = 2 * (theta.last + 2)/(theta.last *
(theta.last + 4)) + fail.cluster(Lambda0.new *
exp(x %*% beta.new), cluster)
z.new = phi.aux * (alpha.aux + phi.aux + 1)/(alpha.aux *
(alpha.aux + phi.aux))
kappa.new = -alpha.aux/(phi.aux * (alpha.aux +
phi.aux)) + digamma(phi.aux + 1) - log(alpha.aux)
theta.new = optim(theta.last, prof.WL, z = z.new,
kappa = kappa.new, method = "Brent", lower = 1e-07,
upper = 20)$par
dif = max(abs(c(beta.last, theta.last) - c(beta.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
kappa.last = kappa.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike, x0 = c(beta.last,
theta.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(beta.new, theta.new)
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "WL"
object.out$tau <- tau.WL(theta.last)
}
object.out
}
frailtyIG <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
prof.IG <- function(theta, z, z1) {
ll = -(1/2) * log(theta) + 1/theta - (1/2) * (z/theta +
z1/theta)
-sum(ll)
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
tau.IG <- function(theta) {
0.5 - 1/theta + (2/theta^2) * exp(2/theta) * expint_E1(2/theta)
}
if (dist == "weibull") {
observed.llike.0.ig.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
rho = eta[1]
lambda = eta[2]
theta = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0, ind) + 1/theta
b = 1/theta
P1 = 0.5 * log(2) + 1/theta - 0.5 * log(pi *
theta) + log(besselK(sqrt(a * b), nu = p)) -
(p/2) * (log(a) - log(b))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ig.dist <- function(eta, t, delta,
x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
theta = eta[ncol(x) + 3]
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0 * exp(x %*% beta),
ind) + 1/theta
b = 1/theta
P1 = 0.5 * log(2) + 1/theta - 0.5 * log(pi *
theta) + log(besselK(sqrt(a * b), nu = p)) -
(p/2) * (log(a) - log(b))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ig.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.ig.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
z1.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
p.aux = r - 1/2
a.aux = 2 * fail.cluster(Lambda0.last, ind) +
1/theta.last
b.aux = 1/theta.last
z.new = sqrt(b.aux/a.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux)
z1.new = sqrt(a.aux/b.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux) - 2 * p.aux/b.aux
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.ig.dist.Q1.0, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
theta.new = optim(theta.last, prof.IG, z = z.new,
z1 = z1.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(rho.last, lambda.last, theta.last) -
c(rho.new, lambda.new, theta.new)))
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ig.dist, x0 = c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.ig.dist(c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster, dist = dist)
para = c(rho.new, lambda.new, theta.new)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
p.aux = r - 1/2
a.aux = 2 * fail.cluster(Lambda0.last * exp(x %*%
beta.last), ind) + 1/theta.last
b.aux = 1/theta.last
z.new = sqrt(b.aux/a.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux)
z1.new = sqrt(a.aux/b.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux) - 2 * p.aux/b.aux
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.ig.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
theta.new = optim(theta.last, prof.IG, z = z.new,
z1 = z1.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(beta.last, rho.last, lambda.last,
theta.last) - c(beta.new, rho.new, lambda.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.ig.dist, x0 = c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.ig.dist(c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
para = c(beta.last, rho.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "IG"
object.out$tau <- tau.IG(theta.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.ig.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
lambda = eta[1:length(part)]
theta = eta[length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0, ind) + 1/theta
b = 1/theta
P1 = 0.5 * log(2) + 1/theta - 0.5 * log(pi *
theta) + log(besselK(sqrt(a * b), nu = p)) -
(p/2) * (log(a) - log(b))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ig.dist <- function(eta, t, delta,
x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
theta = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0 * exp(x %*% beta),
ind) + 1/theta
b = 1/theta
P1 = 0.5 * log(2) + 1/theta - 0.5 * log(pi *
theta) + log(besselK(sqrt(a * b), nu = p)) -
(p/2) * (log(a) - log(b))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ig.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.ig.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
z1.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
p.aux = r - 1/2
a.aux = 2 * fail.cluster(Lambda0.last, ind) +
1/theta.last
b.aux = 1/theta.last
z.new = sqrt(b.aux/a.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux)
z1.new = sqrt(a.aux/b.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux) - 2 * p.aux/b.aux
aux.1 = optim(log(lambda.last), observed.llike.ig.dist.Q1.0,
method = "BFGS", theta = theta.last, z = z.new,
t = t, delta = delta, x = x, ind = ind, dist = dist,
part = part)
lambda.new = exp(aux.1$par[1:length(part)])
theta.new = optim(theta.last, prof.IG, z = z.new,
z1 = z1.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(lambda.last, theta.last) -
c(lambda.new, theta.new)))
theta.last = theta.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ig.dist, x0 = c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.ig.dist(c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
para = c(lambda.new, theta.new)
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
p.aux = r - 1/2
a.aux = 2 * fail.cluster(Lambda0.last * exp(x %*%
beta.last), ind) + 1/theta.last
b.aux = 1/theta.last
z.new = sqrt(b.aux/a.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux)
z1.new = sqrt(a.aux/b.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux) - 2 * p.aux/b.aux
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.ig.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part = part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
theta.new = optim(theta.last, prof.IG, z = z.new,
z1 = z1.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(beta.last, lambda.last, theta.last) -
c(beta.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.ig.dist, x0 = c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.ig.dist(c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
para = c(beta.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "IG"
object.out$tau <- tau.IG(theta.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0.ig <- function(eta, t, delta, ind,
cox.aux) {
theta = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0, ind) + 1/theta
b = 1/theta
P1 = 0.5 * log(2) + 1/theta - 0.5 * log(pi *
theta) + log(besselK(sqrt(a * b), nu = p)) -
(p/2) * (log(a) - log(b))
-sum(P1)
}
observed.llike.ig <- function(eta, t, delta, x, ind,
cox.aux) {
theta = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0 * exp(x %*% beta),
ind) + 1/theta
b = 1/theta
P1 = 0.5 * log(2) + 1/theta - 0.5 * log(pi *
theta) + log(besselK(sqrt(a * b), nu = p)) -
(p/2) * (log(a) - log(b))
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
logz.last = rep(0, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(z.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
p.aux = r - 1/2
a.aux = 2 * fail.cluster(Lambda0.new, ind) +
1/theta.last
b.aux = 1/theta.last
z.new = sqrt(b.aux/a.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux)
z1.new = sqrt(a.aux/b.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux) - 2 * p.aux/b.aux
theta.new = optim(theta.last, prof.IG, z = z.new,
z1 = z1.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(theta.last - theta.new))
theta.last = theta.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ig, x0 = c(theta.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(theta.new)
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(z.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
p.aux = r - 1/2
a.aux = 2 * fail.cluster(Lambda0.new * exp(x %*%
beta.new), ind) + 1/theta.last
b.aux = 1/theta.last
z.new = sqrt(b.aux/a.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux)
z1.new = sqrt(a.aux/b.aux) * besselK(sqrt(a.aux *
b.aux), nu = p.aux + 1)/besselK(sqrt(a.aux *
b.aux), nu = p.aux) - 2 * p.aux/b.aux
theta.new = optim(theta.last, prof.IG, z = z.new,
z1 = z1.new, method = "Brent", lower = 1e-04,
upper = 20)$par
dif = max(abs(c(beta.last, theta.last) - c(beta.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.ig, x0 = c(beta.last,
theta.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(beta.new, theta.new)
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "IG"
object.out$tau <- tau.IG(theta.last)
}
object.out
}
frailtyBS <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
prof.BS <- function(theta, z, z1, y) {
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
ll = (1/2) * (phi + log1p(phi)) + (1 - y) * (log(phi) -
log1p(phi)) - (1/4) * (z * (1 + phi) + z1 * phi^2/(1 +
phi))
-sum(ll)
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
tau.BS <- function(theta) {
aux.int <- function(x, theta) {
phi <- (1 + sqrt(3 * theta + 1) - theta)/theta
ll <- -log(2) + log1p((1 + 4 * x/(phi + 1))^(-1/2)) +
(phi/2) * (1 - (1 + 4 * x/(phi + 1))^(1/2))
L.z <- exp(ll)
p1 <- exp((phi/2) * (1 - (1 + 4 * x/(phi + 1))^(1/2)))
p2 <- (phi * (phi^2 + phi * (4 * x + 3) + 2) *
sqrt(1 + 4 * x/(phi + 1)) + phi^3 + phi^2 *
(4 * x - 5) - 2 * phi * (12 * x + 5) - 4) *
sqrt(1 + 4 * x/(phi + 1)) + 4 * (2 * (phi^2 +
phi * (4 * x + 3) + 2) * sqrt(1 + 4 * x/(phi +
1)) + phi * (phi + 4 * x + 1))
L2.z <- p1 * p2/(2 * (phi + 4 * x + 1)^3)
x * L.z * L2.z
}
4 * integrate(aux.int, lower = 0, upper = Inf, theta = theta)$value -
1
}
if (dist == "weibull") {
observed.llike.0.bs.dist <- function(eta, t, delta,
ind, part = NULL) {
rho = eta[1]
lambda = eta[2]
theta = eta[3]
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0, ind) + (phi + 1)/2
b = phi^2/(2 * (phi + 1))
P1 = phi/2 + (1/2) * log1p(phi) - log(2) - 0.5 *
log(pi) - (p/2) * (log(a) - log(b)) + log(sqrt(b/a) *
besselK(sqrt(a * b), nu = p + 1) + (phi/(1 +
phi)) * besselK(sqrt(a * b), nu = p))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.bs.dist <- function(eta, t, delta,
x, ind, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
theta = eta[ncol(x) + 3]
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0 * exp(x %*% beta),
ind) + (phi + 1)/2
b = phi^2/(2 * (phi + 1))
P1 = phi/2 + (1/2) * log1p(phi) - log(2) - 0.5 *
log(pi) - (p/2) * (log(a) - log(b)) + log(sqrt(b/a) *
besselK(sqrt(a * b), nu = p + 1) + (phi/(1 +
phi)) * besselK(sqrt(a * b), nu = p))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.bs.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.bs.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
z1.last = rep(1, m)
y.last = rep(0.5, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
phi.last = (1 + sqrt(3 * theta.last + 1) -
theta.last)/theta.last
p0.aux = r - 1/2
p1.aux = r + 1 - 1/2
a.aux = 2 * fail.cluster(Lambda0.last, cluster) +
(phi.last + 1)/2
b.aux = phi.last^2/(2 * (phi.last + 1))
xi0.aux = (phi.last/(1 + phi.last)) * besselK(sqrt(a.aux *
b.aux), nu = p0.aux)/(a.aux/b.aux)^(p0.aux/2)
xi1.aux = besselK(sqrt(a.aux * b.aux), nu = p1.aux)/(a.aux/b.aux)^(p1.aux/2)
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
z.new = sqrt(b.aux/a.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux))
z1.new = sqrt(a.aux/b.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux)) -
(2/b.aux) * ((1 - nu.aux) * p0.aux + nu.aux *
p1.aux)
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.bs.dist.Q1.0, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
theta.new = optim(theta.last, prof.BS, z = z.new,
z1 = z1.new, y = y.new, method = "Brent",
lower = 1e-04, upper = 4.9999)$par
dif = max(abs(c(rho.last, lambda.last, theta.last) -
c(rho.new, lambda.new, theta.new)))
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.bs.dist, x0 = c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.bs.dist(c(rho.last,
lambda.last, theta.last), t = t, delta = delta,
ind = cluster)
para = c(rho.new, lambda.new, theta.new)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
phi.last = (1 + sqrt(3 * theta.last + 1) -
theta.last)/theta.last
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
p0.aux = r - 1/2
p1.aux = r + 1 - 1/2
a.aux = 2 * fail.cluster(Lambda0.last * exp(x %*%
beta.last), cluster) + (phi.last + 1)/2
b.aux = phi.last^2/(2 * (phi.last + 1))
xi0.aux = (phi.last/(1 + phi.last)) * besselK(sqrt(a.aux *
b.aux), nu = p0.aux)/(a.aux/b.aux)^(p0.aux/2)
xi1.aux = besselK(sqrt(a.aux * b.aux), nu = p1.aux)/(a.aux/b.aux)^(p1.aux/2)
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
z.new = sqrt(b.aux/a.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux))
z1.new = sqrt(a.aux/b.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux)) -
(2/b.aux) * ((1 - nu.aux) * p0.aux + nu.aux *
p1.aux)
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.bs.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
theta.new = optim(theta.last, prof.BS, z = z.new,
z1 = z1.new, y = y.new, method = "Brent",
lower = 1e-04, upper = 4.9999)$par
dif = max(abs(c(beta.last, rho.last, lambda.last,
theta.last) - c(beta.new, rho.new, lambda.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.bs.dist, x0 = c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.bs.dist(c(beta.last,
rho.last, lambda.last, theta.last), t = t,
x = x, delta = delta, ind = cluster)
para = c(beta.last, rho.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "BS"
object.out$tau <- tau.BS(theta.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.bs.dist <- function(eta, t, delta,
ind, part = NULL) {
lambda = eta[1:length(part)]
theta = eta[length(part) + 1]
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0, ind) + (phi + 1)/2
b = phi^2/(2 * (phi + 1))
P1 = phi/2 + (1/2) * log1p(phi) - log(2) - 0.5 *
log(pi) - (p/2) * (log(a) - log(b)) + log(sqrt(b/a) *
besselK(sqrt(a * b), nu = p + 1) + (phi/(1 +
phi)) * besselK(sqrt(a * b), nu = p))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.bs.dist <- function(eta, t, delta,
x, ind, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
theta = eta[ncol(x) + length(part) + 1]
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0 * exp(x %*% beta),
ind) + (phi + 1)/2
b = phi^2/(2 * (phi + 1))
P1 = phi/2 + (1/2) * log1p(phi) - log(2) - 0.5 *
log(pi) - (p/2) * (log(a) - log(b)) + log(sqrt(b/a) *
besselK(sqrt(a * b), nu = p + 1) + (phi/(1 +
phi)) * besselK(sqrt(a * b), nu = p))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.bs.dist.Q1.0 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.bs.dist.Q1 <- function(eta, theta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
z1.last = rep(1, m)
y.last = rep(0.5, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
phi.last = (1 + sqrt(3 * theta.last + 1) -
theta.last)/theta.last
p0.aux = r - 1/2
p1.aux = r + 1 - 1/2
a.aux = 2 * fail.cluster(Lambda0.last, cluster) +
(phi.last + 1)/2
b.aux = phi.last^2/(2 * (phi.last + 1))
xi0.aux = (phi.last/(1 + phi.last)) * besselK(sqrt(a.aux *
b.aux), nu = p0.aux)/(a.aux/b.aux)^(p0.aux/2)
xi1.aux = besselK(sqrt(a.aux * b.aux), nu = p1.aux)/(a.aux/b.aux)^(p1.aux/2)
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
z.new = sqrt(b.aux/a.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux))
z1.new = sqrt(a.aux/b.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux)) -
(2/b.aux) * ((1 - nu.aux) * p0.aux + nu.aux *
p1.aux)
aux.1 = optim(log(lambda.last), observed.llike.bs.dist.Q1.0,
method = "BFGS", theta = theta.last, z = z.new,
t = t, delta = delta, x = x, ind = ind, dist = dist,
part = part)
lambda.new = exp(aux.1$par[1:length(part)])
theta.new = optim(theta.last, prof.BS, z = z.new,
z1 = z1.new, y = y.new, method = "Brent",
lower = 1e-04, upper = 4.9999)$par
dif = max(abs(c(lambda.last, theta.last) -
c(lambda.new, theta.new)))
theta.last = theta.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.bs.dist, x0 = c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.bs.dist(c(lambda.last,
theta.last), t = t, delta = delta, ind = cluster,
part = part)
para = c(lambda.new, theta.new)
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
while (i <= max.iter & dif > prec) {
phi.last = (1 + sqrt(3 * theta.last + 1) -
theta.last)/theta.last
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
p0.aux = r - 1/2
p1.aux = r + 1 - 1/2
a.aux = 2 * fail.cluster(Lambda0.last * exp(x %*%
beta.last), cluster) + (phi.last + 1)/2
b.aux = phi.last^2/(2 * (phi.last + 1))
xi0.aux = (phi.last/(1 + phi.last)) * besselK(sqrt(a.aux *
b.aux), nu = p0.aux)/(a.aux/b.aux)^(p0.aux/2)
xi1.aux = besselK(sqrt(a.aux * b.aux), nu = p1.aux)/(a.aux/b.aux)^(p1.aux/2)
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
z.new = sqrt(b.aux/a.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux))
z1.new = sqrt(a.aux/b.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux)) -
(2/b.aux) * ((1 - nu.aux) * p0.aux + nu.aux *
p1.aux)
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.bs.dist.Q1, method = "BFGS",
theta = theta.last, z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part = part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
theta.new = optim(theta.last, prof.BS, z = z.new,
z1 = z1.new, y = y.new, method = "Brent",
lower = 1e-04, upper = 4.9999)$par
dif = max(abs(c(beta.last, lambda.last, theta.last) -
c(beta.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
lambda.last = lambda.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.bs.dist, x0 = c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.bs.dist(c(beta.last,
lambda.last, theta.last), t = t, x = x, delta = delta,
ind = cluster, part = part)
para = c(beta.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "BS"
object.out$tau <- tau.BS(theta.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0.bs <- function(eta, t, delta, ind,
cox.aux) {
theta = eta
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0, ind) + (phi + 1)/2
b = phi^2/(2 * (phi + 1))
P1 = phi/2 + (1/2) * log1p(phi) - log(2) - 0.5 *
log(pi) - (p/2) * (log(a) - log(b)) + log(sqrt(b/a) *
besselK(sqrt(a * b), nu = p + 1) + (phi/(1 +
phi)) * besselK(sqrt(a * b), nu = p))
-sum(P1)
}
observed.llike.bs <- function(eta, t, delta, x, ind,
cox.aux) {
theta = eta[length(eta)]
phi = (1 + sqrt(3 * theta + 1) - theta)/theta
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
p = r - 1/2
a = 2 * fail.cluster(Lambda0 * exp(x %*% beta),
ind) + (phi + 1)/2
b = phi^2/(2 * (phi + 1))
P1 = phi/2 + (1/2) * log1p(phi) - log(2) - 0.5 *
log(pi) - (p/2) * (log(a) - log(b)) + log(sqrt(b/a) *
besselK(sqrt(a * b), nu = p + 1) + (phi/(1 +
phi)) * besselK(sqrt(a * b), nu = p))
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
m = length(t)
theta.last = 0.5
z.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
phi.last = (1 + sqrt(3 * theta.last + 1) -
theta.last)/theta.last
cox.aux = coxph(Surv(t, delta) ~ offset(log(z.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
p0.aux = r - 1/2
p1.aux = r + 1 - 1/2
a.aux = 2 * fail.cluster(Lambda0.new, cluster) +
(phi.last + 1)/2
b.aux = phi.last^2/(2 * (phi.last + 1))
xi0.aux = (phi.last/(1 + phi.last)) * besselK(sqrt(a.aux *
b.aux), nu = p0.aux)/(a.aux/b.aux)^(p0.aux/2)
xi1.aux = besselK(sqrt(a.aux * b.aux), nu = p1.aux)/(a.aux/b.aux)^(p1.aux/2)
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
z.new = sqrt(b.aux/a.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux))
z1.new = sqrt(a.aux/b.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux)) -
(2/b.aux) * ((1 - nu.aux) * p0.aux + nu.aux *
p1.aux)
theta.new = optim(theta.last, prof.BS, z = z.new,
z1 = z1.new, y = y.new, method = "Brent",
lower = 1e-04, upper = 4.9999)$par
dif = max(abs(theta.last - theta.new))
theta.last = theta.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.0.bs, x0 = c(theta.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(theta.new)
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
phi.last = (1 + sqrt(3 * theta.last + 1) -
theta.last)/theta.last
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(z.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
p0.aux = r - 1/2
p1.aux = r + 1 - 1/2
a.aux = 2 * fail.cluster(Lambda0.new * exp(x %*%
beta.new), cluster) + (phi.last + 1)/2
b.aux = phi.last^2/(2 * (phi.last + 1))
xi0.aux = (phi.last/(1 + phi.last)) * besselK(sqrt(a.aux *
b.aux), nu = p0.aux)/(a.aux/b.aux)^(p0.aux/2)
xi1.aux = besselK(sqrt(a.aux * b.aux), nu = p1.aux)/(a.aux/b.aux)^(p1.aux/2)
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
z.new = sqrt(b.aux/a.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux))
z1.new = sqrt(a.aux/b.aux) * ((1 - nu.aux) *
besselK(sqrt(a.aux * b.aux), nu = p0.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p0.aux) +
nu.aux * besselK(sqrt(a.aux * b.aux), nu = p1.aux +
1)/besselK(sqrt(a.aux * b.aux), nu = p1.aux)) -
(2/b.aux) * ((1 - nu.aux) * p0.aux + nu.aux *
p1.aux)
theta.new = optim(theta.last, prof.BS, z = z.new,
z1 = z1.new, y = y.new, method = "Brent",
lower = 1e-04, upper = 4.9999)$par
dif = max(abs(c(beta.last, theta.last) - c(beta.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
z.last = z.new
i = i + 1
}
aux.se = hessian(observed.llike.bs, x0 = c(beta.last,
theta.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
para = c(beta.new, theta.new)
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = z.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "BS"
object.out$tau <- tau.BS(theta.last)
}
object.out
}
frailtyTN <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
ratio <- function(x) exp(dnorm(x, log = TRUE) - pnorm(x,
log.p = TRUE))
prof.TN <- function(nu, z, z2) {
ll = log1p(nu + ratio(nu) - 1) - 0.5 * (z2 * (nu +
ratio(nu))^2 - 2 * z * (nu + ratio(nu)) * nu +
nu^2) - pnorm(nu, log.p = TRUE)
-sum(ll)
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
vari.tn <- function(nu) {
(nu + ratio(nu))^(-2) - ratio(nu) * (nu + ratio(nu))^(-1)
}
C.k.aux <- function(x, k, ri, bi, nu) {
exp((ri + k) * log(x) - x * bi - x^2 * 0.5 * (nu +
ratio(nu))^2)
}
C.k <- function(k, ri, bi, nu) {
integrate(C.k.aux, lower = 0, upper = 10, k = k,
ri = ri, bi = bi, nu = nu, abs.tol = 1e-15)$value
}
dg <- function(x) {
-2 * (x + ratio(x))^(-3) * (1 - ratio(x) * (x + ratio(x))) +
ratio(x) + ratio(x) * (x + ratio(x))^(-2) * (1 -
ratio(x) * (x + ratio(x)))
}
tau.TN.aux <- function(x, theta) {
eq.vari.tn = function(nu, theta) {
vari.tn(nu) - theta
}
nu = uniroot(eq.vari.tn, lower = -20, upper = 20,
theta = theta, extendInt = "yes")$root
b = nu + exp(dnorm(nu, log = TRUE) - pnorm(nu, log.p = TRUE))
L = exp(pnorm(nu - x/b, log.p = TRUE) - pnorm(nu,
log.p = TRUE) + (x/b) * (-nu + x/(2 * b)))
L2 = L/b^2 * ((nu - x/b) * (exp(dnorm(nu - x/b, log = TRUE) -
pnorm(nu - x/b, log.p = TRUE)) + nu - x/b) +
1)
pmax(0, x * L^2/b^2 * ((nu - x/b) * (exp(dnorm(nu -
x/b, log = TRUE) - pnorm(nu - x/b, log.p = TRUE)) +
nu - x/b) + 1))
}
L.deriv <- function(s, nu, n) {
gamma = nu + ratio(nu)
alpha = nu - s/gamma
L0 = exp(pnorm(alpha, log.p = TRUE) - pnorm(nu, log.p = TRUE) +
(s/gamma) * (0.5 * s/gamma - nu))
L1 = -(L0/gamma) * (alpha + ratio(alpha))
L2 = (L0/gamma^2) * (alpha * (ratio(alpha) + alpha) +
1)
if (n == 0)
L = L0
if (n == 1)
L = L1
if (n == 2)
L = L2
if (n > 2) {
L = c(L1, L2)
for (i in 3:n) {
L[i] = (i - 1) * L[i - 2]/gamma^2 - alpha *
L[i - 1]/gamma
}
L = L[n]
}
L
}
tau.TN = function(theta) 4 * integrate(tau.TN.aux, lower = 0,
upper = 1000, theta = theta, stop.on.error = FALSE)$value -
1
if (dist == "weibull") {
observed.llike.0.tn.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
rho = eta[1]
lambda = eta[2]
nu = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0, ind)
P1 = log((-1)^r * mapply(L.deriv, s = b, nu = nu,
n = r))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.tn.dist <- function(eta, t, delta,
x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
nu = eta[ncol(x) + 3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0 * exp(x %*% beta), ind)
P1 = log((-1)^r * mapply(L.deriv, s = b, nu = nu,
n = r))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.tn.dist.Q1.0 <- function(eta, z, t,
delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.tn.dist.Q1 <- function(eta, z, t,
delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
nu.last = 2
theta.last = vari.tn(nu.last)
z.last = rep(1, m)
z2.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
b.aux = fail.cluster(Lambda0.last, cluster) -
nu.last * (nu.last + ratio(nu.last))
C0 = mapply(C.k, k = 0, ri = r, bi = b.aux,
nu = nu.last)
C1 = mapply(C.k, k = 1, ri = r, bi = b.aux,
nu = nu.last)
C2 = mapply(C.k, k = 2, ri = r, bi = b.aux,
nu = nu.last)
z.new = C1/C0
z2.new = C2/C0
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.tn.dist.Q1.0, method = "BFGS",
z = z.new, t = t, delta = delta, x = x, ind = ind,
dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
nu.new = optim(nu.last, prof.TN, z = z.new,
z2 = z2.new, method = "Brent", lower = max(-40,
nu.last - 3), upper = min(40, nu.last +
3))$par
theta.new = vari.tn(nu.new)
dif = max(abs(c(rho.last, lambda.last, theta.last) -
c(rho.new, lambda.new, theta.new)))
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
nu.last = nu.new
z.last = z.new
z2.last = z2.new
i = i + 1
}
aux.se = hessian(observed.llike.0.tn.dist, x0 = c(rho.last,
lambda.last, nu.last), t = t, delta = delta,
ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(nu.last))
llike.obs = -observed.llike.0.tn.dist(c(rho.last,
lambda.last, nu.last), t = t, delta = delta,
ind = cluster, dist = dist)
para = c(rho.new, lambda.new, theta.new)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
b.aux = fail.cluster(Lambda0.last * exp(x %*%
beta.last), cluster) - nu.last * (nu.last +
ratio(nu.last))
C0 = mapply(C.k, k = 0, ri = r, bi = b.aux,
nu = nu.last)
C1 = mapply(C.k, k = 1, ri = r, bi = b.aux,
nu = nu.last)
C2 = mapply(C.k, k = 2, ri = r, bi = b.aux,
nu = nu.last)
z.new = C1/C0
z2.new = C2/C0
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.tn.dist.Q1, method = "BFGS",
z = z.new, t = t, delta = delta, x = x, ind = ind,
dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
nu.new = optim(nu.last, prof.TN, z = z.new,
z2 = z2.new, method = "Brent", lower = -40,
upper = 40)$par
theta.new = vari.tn(nu.new)
dif = max(abs(c(beta.last, rho.last, lambda.last,
theta.last) - c(beta.new, rho.new, lambda.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
nu.last = nu.new
z.last = z.new
z2.last = z2.new
i = i + 1
}
aux.se = hessian(observed.llike.tn.dist, x0 = c(beta.last,
rho.last, lambda.last, nu.last), t = t, x = x,
delta = delta, ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(nu.last))
llike.obs = -observed.llike.tn.dist(c(beta.last,
rho.last, lambda.last, nu.last), t = t, x = x,
delta = delta, ind = cluster, dist = dist)
para = c(beta.new, rho.new, lambda.new, theta.new)
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new, z2 = z2.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "TN"
object.out$tau <- tau.TN(theta.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.tn.dist <- function(eta, t, delta,
ind, dist, part = NULL) {
lambda = eta[1:length(part)]
nu = eta[length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0, ind)
P1 = log((-1)^r * mapply(L.deriv, s = b, nu = nu,
n = r))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.tn.dist <- function(eta, t, delta,
x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
nu = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0 * exp(x %*% beta), ind)
P1 = log((-1)^r * mapply(L.deriv, s = b, nu = nu,
n = r))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.tn.dist.Q1.0 <- function(eta, z, t,
delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.tn.dist.Q1 <- function(eta, z, t,
delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
m = length(t)
nu.last = 2
theta.last = vari.tn(nu.last)
z.last = rep(1, m)
z2.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
b.aux = fail.cluster(Lambda0.last, cluster) -
nu.last * (nu.last + ratio(nu.last))
C0 = mapply(C.k, k = 0, ri = r, bi = b.aux,
nu = nu.last)
C1 = mapply(C.k, k = 1, ri = r, bi = b.aux,
nu = nu.last)
C2 = mapply(C.k, k = 2, ri = r, bi = b.aux,
nu = nu.last)
z.new = C1/C0
z2.new = C2/C0
aux.1 = optim(log(lambda.last), observed.llike.tn.dist.Q1.0,
method = "BFGS", z = z.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part = part)
lambda.new = exp(aux.1$par)
nu.new = optim(nu.last, prof.TN, z = z.new,
z2 = z2.new, method = "Brent", lower = -40,
upper = 40)$par
theta.new = vari.tn(nu.new)
dif = max(abs(c(lambda.last, theta.last) -
c(lambda.new, theta.new)))
theta.last = theta.new
lambda.last = lambda.new
nu.last = nu.new
z.last = z.new
z2.last = z2.new
i = i + 1
}
aux.se = hessian(observed.llike.0.tn.dist, x0 = c(lambda.last,
nu.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.tn.dist(c(lambda.last,
nu.last), t = t, delta = delta, ind = cluster,
dist = dist, part = part)
se = abs(c(lambda.last, theta.last))
se[length(se)] = se[length(se)] * abs(dg(nu.last))
llike.obs = lambda.last
para = c(lambda.new, theta.new)
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
b.aux = fail.cluster(Lambda0.last * exp(x %*%
beta.last), cluster) - nu.last * (nu.last +
ratio(nu.last))
C0 = mapply(C.k, k = 0, ri = r, bi = b.aux,
nu = nu.last)
C1 = mapply(C.k, k = 1, ri = r, bi = b.aux,
nu = nu.last)
C2 = mapply(C.k, k = 2, ri = r, bi = b.aux,
nu = nu.last)
z.new = C1/C0
z2.new = C2/C0
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.tn.dist.Q1, method = "BFGS",
z = z.new, t = t, delta = delta, x = x, ind = ind,
dist = dist, part = part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
nu.new = optim(nu.last, prof.TN, z = z.new,
z2 = z2.new, method = "Brent", lower = -40,
upper = 40)$par
theta.new = vari.tn(nu.new)
dif = max(abs(c(beta.last, lambda.last, theta.last) -
c(beta.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
lambda.last = lambda.new
nu.last = nu.new
z.last = z.new
z2.last = z2.new
i = i + 1
}
aux.se = hessian(observed.llike.tn.dist, x0 = c(beta.last,
lambda.last, nu.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(nu.last))
llike.obs = -observed.llike.tn.dist(c(beta.last,
lambda.last, nu.last), t = t, x = x, delta = delta,
ind = cluster, dist = dist, part = part)
para = c(beta.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = z.new, z2 = z2.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "TN"
object.out$tau <- tau.TN(theta.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0.tn <- function(eta, t, delta, ind,
cox.aux) {
nu = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0, ind)
P1 = log((-1)^r * mapply(L.deriv, s = b, nu = nu,
n = r))
-sum(P1)
}
observed.llike.tn <- function(eta, t, delta, x, ind,
cox.aux) {
nu = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0 * exp(x %*% beta), ind)
P1 = log((-1)^r * mapply(L.deriv, s = b, nu = nu,
n = r))
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
m = length(t)
nu.last = 2
theta.last = vari.tn(nu.last)
z.last = rep(1, m)
z2.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(z.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
b.aux = fail.cluster(Lambda0.new, cluster) -
nu.last * (nu.last + ratio(nu.last))
C0 = mapply(C.k, k = 0, ri = r, bi = b.aux,
nu = nu.last)
C1 = mapply(C.k, k = 1, ri = r, bi = b.aux,
nu = nu.last)
C2 = mapply(C.k, k = 2, ri = r, bi = b.aux,
nu = nu.last)
z.new = C1/C0
z2.new = C2/C0
nu.new = optim(nu.last, prof.TN, z = z.new,
z2 = z2.new, method = "Brent", lower = -40,
upper = 40)$par
theta.new = vari.tn(nu.new)
dif = max(abs(theta.last - theta.new))
theta.last = theta.new
nu.last = nu.new
z.last = z.new
z2.last = z2.new
i = i + 1
}
aux.se = hessian(observed.llike.0.tn, x0 = c(nu.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(nu.last))
para = c(theta.new)
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(z.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
b.aux = fail.cluster(Lambda0.new * exp(x %*%
beta.last), cluster) - nu.last * (nu.last +
ratio(nu.last))
C0 = mapply(C.k, k = 0, ri = r, bi = b.aux,
nu = nu.last)
C1 = mapply(C.k, k = 1, ri = r, bi = b.aux,
nu = nu.last)
C2 = mapply(C.k, k = 2, ri = r, bi = b.aux,
nu = nu.last)
z.new = C1/C0
z2.new = C2/C0
nu.new = optim(nu.last, prof.TN, z = z.new,
z2 = z2.new, method = "Brent", lower = -40,
upper = 40)$par
theta.new = vari.tn(nu.new)
dif = max(abs(c(beta.last, theta.last) - c(beta.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
nu.last = nu.new
z.last = z.new
z2.last = z2.new
i = i + 1
}
aux.se = hessian(observed.llike.tn, x0 = c(beta.last,
nu.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(nu.last))
para = c(beta.new, theta.new)
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = z.new, z2 = z2.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "TN"
object.out$tau <- tau.TN(theta.last)
}
object.out
}
frailtyMIG <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
Var_U <- function(m) {(m-1)*(m+1)*(m^2-m-1)/m^2}
prof.MIG <- function(m, u1, y, a.0, a.1, b, p, nu.aux) {
kappa = (1-nu.aux)*m^(-2)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*(m+1)^(-2)*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
kappa1 = (m+1)^(-1)*nu.aux+m^(-1)*(1-nu.aux)
ll = log(m)-0.5*log1p(-m)-0.5*(m^2/(1-m))*(u1+kappa-2*kappa1)+
y*log1p(-m)+(1-y)*log(m)
-sum(ll)
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
tau.MIG <- function(m) {
aux.int <- function(x, m) {
a1<-(m/(1-m))*(1-sqrt(1+2*(1-m)*x))
a2<-(m^2/(1-m^2))*(1-sqrt(1+2*(1+m)^2*(1-m)*x/m^2))
L.z <- m*exp(a1)+(1-m)*exp(a2)
p1<-m^2*(exp(a1)*(1-m)*(1+2*(1-m)*x)^(-3/2)+(1+2*(1-m)*x)^(-1)*m*exp(a1))
p2<-(1-m^2)*(exp(a2)*(1+m)^2*(1-m)*m^(-2)*(1+2*(1+m)^2*(1-m)*x/m^2)^(-3/2)+(1+2*(1+m)^2*(1-m)*x/m^2)^(-1)*exp(a2)*(1+m))
L2.z <- p1+p2
x * L.z * L2.z
}
4 * integrate(aux.int, lower = 0, upper = Inf, m=m)$value - 1
}
if (dist == "weibull") {
observed.llike.0.mig.dist <- function(eta, t, delta,
ind, part = NULL) {
rho = eta[1]
lambda = eta[2]
m = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
a.0=2*A+1/(1-m)
a.1=2*A+m^2/((m+1)^2*(1-m))
b=m^2/(1-m)
p=r-1/2
aux1=m*exp(m/(1-m))*2*besselK(sqrt(a.0*b),nu=p)*(a.0/b)^(-p/2)
aux2=(1-m)*exp(m^2/(1-m^2))*2*besselK(sqrt(a.1*b),nu=p)*(a.1/b)^(-p/2)
P1 = log(m)-0.5*(log(2*pi)+log1p(-m))+log(aux1+aux2)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mig.dist <- function(eta, t, delta,
x, ind, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
m = eta[ncol(x) + 3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta), cluster)
a.0=2*A+1/(1-m)
a.1=2*A+m^2/((m+1)^2*(1-m))
b=m^2/(1-m)
p=r-1/2
aux1=m*exp(m/(1-m))*2*besselK(sqrt(a.0*b),nu=p)*(a.0/b)^(-p/2)
aux2=(1-m)*exp(m^2/(1-m^2))*2*besselK(sqrt(a.1*b),nu=p)*(a.1/b)^(-p/2)
P1 = log(m)-0.5*(log(2*pi)+log1p(-m))+log(aux1+aux2)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mig.dist.Q1.0 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.mig.dist.Q1 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
mm = length(t)
m.last = 0.75
u.last = rep(1, mm)
u1.last = rep(1, mm)
y.last = rep(0.5, mm)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
a.0=2*A+1/(1-m.last)
a.1=2*A+m.last^2/((m.last+1)^2*(1-m.last))
b=m.last^2/(1-m.last)
p=r-1/2
xi0.aux = m.last*besselK(sqrt(a.0*b), nu=p)*a.0^(-p/2)*exp(m.last/(1-m.last))
xi1.aux = (1-m.last)*besselK(sqrt(a.1*b), nu=p)*a.1^(-p/2)*exp(m.last^2/((m.last+1)*(1-m.last)))
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
u.new = (1-nu.aux)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
u1.new = (1-nu.aux)*sqrt(a.0)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(a.1)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.1*b), nu=p))-2*p/b
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.mig.dist.Q1.0, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
m.new = optimize(prof.MIG, c(0.0001,0.9999),
u1 = u1.new, y = y.new, a.0=a.0, a.1=a.1, b=b, p=p, nu.aux=nu.aux)$minimum
dif = max(abs(c(rho.last, lambda.last, m.last) -
c(rho.new, lambda.new, m.new)))
m.last = m.new
rho.last = rho.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
u1.last = u1.new
i = i + 1
}
aux.se = hessian(observed.llike.0.mig.dist, x0 = c(rho.last,
lambda.last, m.last), t = t, delta = delta,
ind = cluster)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.mig.dist(c(rho.last,
lambda.last, m.last), t = t, delta = delta,
ind = cluster)
para = c(rho.new, lambda.new, Var_U(m.new))
se[length(se)] = se[length(se)]*abs((2*m.new^4-m.new^3-m.new-2)/m.new^3)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull", data=data)
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta.last), cluster)
a.0=2*A+1/(1-m.last)
a.1=2*A+m.last^2/((m.last+1)^2*(1-m.last))
b=m.last^2/(1-m.last)
p=r-1/2
xi0.aux = m.last*besselK(sqrt(a.0*b), nu=p)*a.0^(-p/2)*exp(m.last/(1-m.last))
xi1.aux = (1-m.last)*besselK(sqrt(a.1*b), nu=p)*a.1^(-p/2)*exp(m.last^2/((m.last+1)*(1-m.last)))
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
u.new = (1-nu.aux)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
u1.new = (1-nu.aux)*sqrt(a.0)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(a.1)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.1*b), nu=p))-2*p/b
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.mig.dist.Q1, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
m.new = optimize(prof.MIG, c(0.0001,0.9999),
u1 = u1.new, y = y.new, a.0=a.0, a.1=a.1, b=b, p=p, nu.aux=nu.aux)$minimum
dif = max(abs(c(beta.last, rho.last, lambda.last,
m.last) - c(beta.new, rho.new, lambda.new,
m.new)))
beta.last = beta.new
m.last = m.new
rho.last = rho.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
u1.last = u1.new
i = i + 1
}
aux.se = hessian(observed.llike.mig.dist, x0 = c(beta.last,
rho.last, lambda.last, m.last), t = t,
x = x, delta = delta, ind = cluster)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs((2*m.new^4-m.new^3-m.new-2)/m.new^3)
llike.obs = -observed.llike.mig.dist(c(beta.last,
rho.last, lambda.last, m.last), t = t,
x = x, delta = delta, ind = cluster)
para = c(beta.new, rho.new, lambda.new, Var_U(m.new))
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z = u.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "MIG"
object.out$tau <- tau.MIG(m.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.mig.dist <- function(eta, t, delta,
ind, part = NULL) {
lambda = eta[1:length(part)]
m = eta[length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
a.0=2*A+1/(1-m)
a.1=2*A+m^2/((m+1)^2*(1-m))
b=m^2/(1-m)
p=r-1/2
aux1=m*exp(m/(1-m))*2*besselK(sqrt(a.0*b),nu=p)*(a.0/b)^(-p/2)
aux2=(1-m)*exp(m^2/(1-m^2))*2*besselK(sqrt(a.1*b),nu=p)*(a.1/b)^(-p/2)
P1 = log(m)-0.5*(log(2*pi)+log1p(-m))+log(aux1+aux2)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mig.dist <- function(eta, t, delta,
x, ind, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
m = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta), cluster)
a.0=2*A+1/(1-m)
a.1=2*A+m^2/((m+1)^2*(1-m))
b=m^2/(1-m)
p=r-1/2
aux1=m*exp(m/(1-m))*2*besselK(sqrt(a.0*b),nu=p)*(a.0/b)^(-p/2)
aux2=(1-m)*exp(m^2/(1-m^2))*2*besselK(sqrt(a.1*b),nu=p)*(a.1/b)^(-p/2)
P1 = log(m)-0.5*(log(2*pi)+log1p(-m))+log(aux1+aux2)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mig.dist.Q1.0 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.mig.dist.Q1 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
mm = length(t)
m.last = 0.75
u.last = rep(1, mm)
u1.last = rep(1, mm)
y.last = rep(0.5, mm)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
a.0=2*A+1/(1-m.last)
a.1=2*A+m.last^2/((m.last+1)^2*(1-m.last))
b=m.last^2/(1-m.last)
p=r-1/2
xi0.aux = m.last*besselK(sqrt(a.0*b), nu=p)*a.0^(-p/2)*exp(m.last/(1-m.last))
xi1.aux = (1-m.last)*besselK(sqrt(a.1*b), nu=p)*a.1^(-p/2)*exp(m.last^2/((m.last+1)*(1-m.last)))
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
u.new = (1-nu.aux)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
u1.new = (1-nu.aux)*sqrt(a.0)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(a.1)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.1*b), nu=p))-2*p/b
aux.1 = optim(c(log(lambda.last)),
observed.llike.mig.dist.Q1.0, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part=part)
lambda.new = exp(aux.1$par)
m.new = optimize(prof.MIG, c(0.0001,0.9999),
u1 = u1.new, y = y.new, a.0=a.0, a.1=a.1, b=b, p=p, nu.aux=nu.aux)$minimum
dif = max(abs(c(lambda.last, m.last) -
c(lambda.new, m.new)))
m.last = m.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
u1.last = u1.new
i = i + 1
}
aux.se = hessian(observed.llike.0.mig.dist, x0 = c(lambda.last,
m.last), t = t, delta = delta, ind = cluster,
part = part)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs((2*m.new^4-m.new^3-m.new-2)/m.new^3)
llike.obs = -observed.llike.0.mig.dist(c(lambda.last,
m.last), t = t, delta = delta, ind = cluster,
part = part)
para = c(lambda.new, Var_U(m.new))
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta.last), cluster)
a.0=2*A+1/(1-m.last)
a.1=2*A+m.last^2/((m.last+1)^2*(1-m.last))
b=m.last^2/(1-m.last)
p=r-1/2
xi0.aux = m.last*besselK(sqrt(a.0*b), nu=p)*a.0^(-p/2)*exp(m.last/(1-m.last))
xi1.aux = (1-m.last)*besselK(sqrt(a.1*b), nu=p)*a.1^(-p/2)*exp(m.last^2/((m.last+1)*(1-m.last)))
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
u.new = (1-nu.aux)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
u1.new = (1-nu.aux)*sqrt(a.0)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(a.1)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.1*b), nu=p))-2*p/b
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.mig.dist.Q1, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part=part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
m.new = optimize(prof.MIG, c(0.0001,0.9999),
u1 = u1.new, y = y.new, a.0=a.0, a.1=a.1, b=b, p=p, nu.aux=nu.aux)$minimum
dif = max(abs(c(beta.last, lambda.last,
m.last) - c(beta.new, lambda.new,
m.new)))
beta.last = beta.new
m.last = m.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
u1.last = u1.new
i = i + 1
}
aux.se = hessian(observed.llike.mig.dist, x0 = c(beta.last,
lambda.last, m.last), t = t, x = x, delta = delta,
ind = cluster, part = part)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs((2*m.new^4-m.new^3-m.new-2)/m.new^3)
llike.obs = -observed.llike.mig.dist(c(beta.last,
lambda.last, m.last), t = t, x = x, delta = delta,
ind = cluster, part = part)
para = c(beta.last, lambda.last, Var_U(m.last))
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = u.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "MIG"
object.out$tau <- tau.MIG(m.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0.mig <- function(eta, t, delta, ind,
cox.aux) {
m = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0, cluster)
a.0=2*A+1/(1-m)
a.1=2*A+m^2/((m+1)^2*(1-m))
b=m^2/(1-m)
p=r-1/2
aux1=m*exp(m/(1-m))*2*besselK(sqrt(a.0*b),nu=p)*(a.0/b)^(-p/2)
aux2=(1-m)*exp(m^2/(1-m^2))*2*besselK(sqrt(a.1*b),nu=p)*(a.1/b)^(-p/2)
P1 = log(m)-0.5*(log(2*pi)+log1p(-m))+log(aux1+aux2)
-sum(P1)
}
observed.llike.mig <- function(eta, t, delta, x, ind,
cox.aux) {
m = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0*exp(x %*% beta), cluster)
a.0=2*A+1/(1-m)
a.1=2*A+m^2/((m+1)^2*(1-m))
b=m^2/(1-m)
p=r-1/2
aux1=m*exp(m/(1-m))*2*besselK(sqrt(a.0*b),nu=p)*(a.0/b)^(-p/2)
aux2=(1-m)*exp(m^2/(1-m^2))*2*besselK(sqrt(a.1*b),nu=p)*(a.1/b)^(-p/2)
P1 = log(m)-0.5*(log(2*pi)+log1p(-m))+log(aux1+aux2)
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
mm = length(t)
m.last = 0.75
u.last = rep(1, mm)
u1.last = rep(1, mm)
y.last = rep(0.5, mm)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(u.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.new, cluster)
a.0=2*A+1/(1-m.last)
a.1=2*A+m.last^2/((m.last+1)^2*(1-m.last))
b=m.last^2/(1-m.last)
p=r-1/2
xi0.aux = m.last*besselK(sqrt(a.0*b), nu=p)*a.0^(-p/2)*exp(m.last/(1-m.last))
xi1.aux = (1-m.last)*besselK(sqrt(a.1*b), nu=p)*a.1^(-p/2)*exp(m.last^2/((m.last+1)*(1-m.last)))
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
u.new = (1-nu.aux)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
u1.new = (1-nu.aux)*sqrt(a.0)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(a.1)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.1*b), nu=p))-2*p/b
m.new = optimize(prof.MIG, c(0.0001,0.9999),
u1 = u1.new, y = y.new, a.0=a.0, a.1=a.1, b=b, p=p, nu.aux=nu.aux)$minimum
dif = max(abs(m.last - m.new))
m.last = m.new
y.last = y.new
u.last = u.new
u1.last = u1.new
i = i + 1
}
aux.se = hessian(observed.llike.0.mig, x0 = c(m.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs((2*m.new^4-m.new^3-m.new-2)/m.new^3)
para = c(Var_U(m.new))
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(u.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.new*exp(x %*% beta.last), cluster)
a.0=2*A+1/(1-m.last)
a.1=2*A+m.last^2/((m.last+1)^2*(1-m.last))
b=m.last^2/(1-m.last)
p=r-1/2
xi0.aux = m.last*besselK(sqrt(a.0*b), nu=p)*a.0^(-p/2)*exp(m.last/(1-m.last))
xi1.aux = (1-m.last)*besselK(sqrt(a.1*b), nu=p)*a.1^(-p/2)*exp(m.last^2/((m.last+1)*(1-m.last)))
nu.aux = xi1.aux/(xi0.aux + xi1.aux)
y.new = nu.aux
u.new = (1-nu.aux)*sqrt(b)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(a.0)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(b)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(a.1)*besselK(sqrt(a.1*b), nu=p))
u1.new = (1-nu.aux)*sqrt(a.0)*besselK(sqrt(a.0*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.0*b), nu=p))+nu.aux*sqrt(a.1)*besselK(sqrt(a.1*b), nu=p+1)/(sqrt(b)*besselK(sqrt(a.1*b), nu=p))-2*p/b
m.new = optimize(prof.MIG, c(0.0001,0.9999),
u1 = u1.new, y = y.new, a.0=a.0, a.1=a.1, b=b, p=p, nu.aux=nu.aux)$minimum
dif = max(abs(c(beta.last, m.last) - c(beta.new, m.new)))
beta.last = beta.new
m.last = m.new
y.last = y.new
u.last = u.new
u1.last = u1.new
i = i + 1
}
aux.se = hessian(observed.llike.mig, x0 = c(beta.last,
m.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs((2*m.new^4-m.new^3-m.new-2)/m.new^3)
para = c(beta.new, Var_U(m.new))
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = u.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "MIG"
object.out$tau <- tau.MIG(m.last)
}
object.out
}
frailtyMBS <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part = NULL) {
Var_U <- function(m) {(2 * m^2 / (1 - m) + 5) / ((m^2 / (1 - m) + 1)^2) * (m - m^2 + 1) + m * (1 - m)}
d.Var_U <- function(m) {((2 * m * (2 - m) / (1 - m)^2) * ((m^2 / (1 - m)) + 1)^2 -
((2 * m^2 / (1 - m)) + 5) * (2 * ((m^2 / (1 - m)) + 1) * (m * (2 - m) / (1 - m)^2))) /
(((m^2 / (1 - m)) + 1)^4) * (m - m^2 + 1) +
(((2 * m^2 / (1 - m)) + 5) / (((m^2 / (1 - m)) + 1)^2) * (1 - 2 * m) + (1 - m))-m}
prof.MBS <- function(m, u, y, v, a.0, a.1, nu.0, nu.1, b.0, b.1, p.0, p.1) {
gamma=m^2/(2*m^2+1-m)
##kappa es E((V_i-1/2)*log(m+Y_i))
kappa=E.g(m, type=4, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
##kappa1 es E(Ui/(m+Yi))
kappa1=E.g(m, type=2, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
##kappa2 es E((m+Yi)/Ui)
kappa2=E.g(m, type=3, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
ll = y*log1p(-m)+(1-y)*log(m)+kappa+(1/2-v)*log1p((1-m)/m^2)-log(besselK(0.5*m^2/(1-m), nu=1/2))-(m^2/(1-m)+1)*kappa1/4-m^4*(1-m)^(-1)*(m^2+1-m)^(-1)*kappa2/4
-sum(ll)
}
a<-function(y, a.0, a.1){ifelse(y==0, a.0, a.1)}
b<-function(y, b.0, b.1){ifelse(y==0, b.0, b.1)}
p<-function(y, p.0, p.1){ifelse(y==0, p.0, p.1)}
nu<-function(y, nu.0, nu.1){ifelse(y==0, nu.0, nu.1)}
g<-function(y, v, m, a.0, a.1, b.0, b.1, p.0, p.1, type=1)
{
##type=1 for Esp(U|...), type=2 for E(U/(m+Y)|...), type=3 for E((m+Y)/U|...) and type=4 for E((V-1/2)*log(m+Y)|...)
if(type==1){gg=sqrt(b(y,b.0, b.1)/a(y,a.0, a.1))*besselK(sqrt(a(y,a.0, a.1)*b(y,b.0, b.1)),nu=p(v,p.0, p.1)+1)/besselK(sqrt(a(y,a.0, a.1)*b(y,b.0, b.1)),nu=p(v,p.0, p.1))}
if(type==2){gg=(m+y)^(-1)*sqrt(b(y,b.0, b.1)/a(y,a.0, a.1))*besselK(sqrt(a(y,a.0, a.1)*b(y,b.0, b.1)),nu=p(v,p.0, p.1)+1)/besselK(sqrt(a(y,a.0, a.1)*b(y,b.0, b.1)),nu=p(v,p.0, p.1))}
if(type==3){gg=(m+y)*sqrt(a(y,a.0, a.1)/b(y,b.0, b.1))*besselK(sqrt(a(y,a.0, a.1)*b(y,b.0, b.1)),nu=p(v,p.0, p.1)+1)/besselK(sqrt(a(y,a.0, a.1)*b(y,b.0, b.1)),nu=p(v,p.0, p.1))-2*(m+y)*p(v,p.0, p.1)/b(y,b.0, b.1)}
if(type==4){gg=(v-1/2)*log(m+y)}
gg
}
E.g<-function(m, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
{
zeros=rep(0, length(a.0))
ones=rep(1, length(a.0))
g(zeros,zeros,m,a.0, a.1, b.0, b.1, p.0, p.1,type)*(1-gamma)+g(zeros,ones,m,a.0, a.1, b.0, b.1, p.0, p.1,type)*gamma+(g(ones,zeros,m,a.0, a.1, b.0, b.1, p.0, p.1,type)-g(zeros,zeros,m,a.0, a.1, b.0, b.1, p.0, p.1,type))*nu(zeros,nu.0,nu.1)*(1-gamma)+(g(ones,ones,m,a.0, a.1, b.0, b.1, p.0, p.1,type)-g(zeros,ones,m,a.0, a.1, b.0, b.1, p.0, p.1,type))*nu(ones,nu.0,nu.1)*gamma
}
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
###
##Hay que modificar este tau
###
tau.MBS <- function(m) {
L_MBS_deriv <- function(t, m, n) {
a <- (m^2 / (1 - m) + 1) / (2 * m)
b <- (m * (m^2 / (1 - m))^2) / (2 * (m^2 / (1 - m) + 1))
c <- (m^2 / (1 - m) + 1) / (2 * (m + 1))
d <- ((m + 1) * (m^2 / (1 - m))^2) / (2 * (m^2 / (1 - m) + 1))
# Define the modified Bessel functions
K <- function(nu, x) {
return(besselK(x, nu))
}
sqrt_term_m <- sqrt(m^2 / (1 - m) + 1)
sqrt_term_m1 <- sqrt((m^2 / (1 - m)) + 1 + 4 * m * t)
sqrt_term_m1_m1 <- sqrt((m^2 / (1 - m)) + 1 + 4 * (m + 1) * t)
L_m <- exp((m^2 / (1 - m) * (sqrt_term_m - sqrt_term_m1)) / (2 * sqrt_term_m))
L_m1 <- exp((m^2 / (1 - m) * (sqrt_term_m - sqrt_term_m1_m1)) / (2 * sqrt_term_m))
term1 <- (1 / 2) * m * L_m * (1 + sqrt_term_m / sqrt_term_m1)
term2 <- (1 / 2) * (1 - m) * L_m1 * (1 + sqrt_term_m / sqrt_term_m1_m1)
L0 <- term1 + term2
term1_deriv <- -(m / 2) *(( (a / b)^(1/4) / (sqrt((a + 2 * t) / b))^(3/2) * K(3/2, sqrt(b * (a + 2 * t))) / K(1/2, sqrt(a * b))) +
((a / b)^(-1/4) / (sqrt((a + 2 * t) / b))^(1/2) * K(1/2, sqrt(b * (a + 2 * t))) / K(-1/2, sqrt(a * b))))
term2_deriv <- -(1 - m) / 2 * (( (c / d)^(1/4) / (sqrt((c + 2 * t) / d))^(3/2) * K(3/2, sqrt(d * (c + 2 * t))) / K(1/2, sqrt(c * d))) +
((c / d)^(-1/4) / (sqrt((c + 2 * t) / d))^(1/2) * K(1/2, sqrt(d * (c + 2 * t))) / K(-1/2, sqrt(c * d))))
L1 <- term1_deriv + term2_deriv
term1_deriv2 <-(m / 2)*(( (a / b)^(1/4) / (sqrt((a + 2 * t) / b))^(5/2) * K(5/2, sqrt(b * (a + 2 * t))) / K(1/2, sqrt(a * b))) +
((a / b)^(-1/4) / (sqrt((a + 2 * t) / b))^(3/2) * K(3/2, sqrt(b * (a + 2 * t))) / K(-1/2, sqrt(a * b))))
term2_deriv2 <- (1 - m) / 2 * (( (c / d)^(1/4) / (sqrt((c + 2 * t) / d))^(5/2) * K(5/2, sqrt(d * (c + 2 * t))) / K(1/2, sqrt(c * d))) +
((c / d)^(-1/4) / (sqrt((c + 2 * t) / d))^(3/2) * K(3/2, sqrt(d * (c + 2 * t))) / K(-1/2, sqrt(c * d))))
L2 <- term1_deriv2 + term2_deriv2
if (n == 0) L = L0
if (n == 1) L = L1
if (n == 2) L = L2
L
}
aux.int <- function(x, m) {
x*L_MBS_deriv(x, m, 0)*L_MBS_deriv(x, m, 2)
}
4 * integrate(aux.int, lower = 0, upper = Inf, m=m)$value - 1
}
if (dist == "weibull") {
observed.llike.0.mbs.dist <- function(eta, t, delta,
ind, part = NULL) {
rho = eta[1]
lambda = eta[2]
m = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
theta=m^2/(1-m)
p=r-1/2
a.m=2*A+(theta+1)/(2*m)
b.m=m*theta^2/(2*(theta+1))
a.m1=2*A+(theta+1)/(2*(m+1))
b.m1=(m+1)*theta^2/(2*(theta+1))
aux1=m*(sqrt(pi*m))^(-1)*(b.m/a.m)^(p/2)*(m*theta*(theta+1)^(-1)*besselK(sqrt(a.m*b.m),nu=p)+sqrt(b.m/a.m)*besselK(sqrt(a.m*b.m),nu=p+1))
aux2=(1-m)*(sqrt(pi*(m+1)))^(-1)*(b.m1/a.m1)^(p/2)*((m+1)*theta*(theta+1)^(-1)*besselK(sqrt(a.m1*b.m1),nu=p)+sqrt(b.m1/a.m1)*besselK(sqrt(a.m1*b.m1),nu=p+1))
P1 = (theta/2)-log(2)+(1/2)*log1p(theta)+log(aux1+aux2)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mbs.dist <- function(eta, t, delta,
x, ind, part = NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
m = eta[ncol(x) + 3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta), cluster)
theta=m^2/(1-m)
p=r-1/2
a.m=2*A+(theta+1)/(2*m)
b.m=m*theta^2/(2*(theta+1))
a.m1=2*A+(theta+1)/(2*(m+1))
b.m1=(m+1)*theta^2/(2*(theta+1))
aux1=m*(sqrt(pi*m))^(-1)*(b.m/a.m)^(p/2)*(m*theta*(theta+1)^(-1)*besselK(sqrt(a.m*b.m),nu=p)+sqrt(b.m/a.m)*besselK(sqrt(a.m*b.m),nu=p+1))
aux2=(1-m)*(sqrt(pi*(m+1)))^(-1)*(b.m1/a.m1)^(p/2)*((m+1)*theta*(theta+1)^(-1)*besselK(sqrt(a.m1*b.m1),nu=p)+sqrt(b.m1/a.m1)*besselK(sqrt(a.m1*b.m1),nu=p+1))
P1 = (theta/2)-log(2)+(1/2)*log1p(theta)+log(aux1+aux2)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mbs.dist.Q1.0 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.mbs.dist.Q1 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
mm = length(t)
m.last = 0.5
u.last = rep(1, mm)
v.last = rep(0.5, mm)
y.last = rep(0.5, mm)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
gamma=m.last^2/(2*m.last^2+1-m.last)
a.0=2*A+(m.last^2/(1-m.last)+1)/(2*m.last)
a.1=2*A+(m.last^2/(1-m.last)+1)/(2*(m.last+1))
b.0=(m.last)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
b.1=(m.last+1)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
p.0=r+1/2
p.1=r+1/2-1
nu.0=besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(0-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.0)*(a.0/b.0)^(-p.0/2)*m.last^(0+1/2)+besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(0-1/2)*(1-m.last))
nu.1=besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.1)*(a.0/b.0)^(-p.1/2)*m.last^(1+1/2)+besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last))
u.new=E.g(m.last, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
y.new=nu.0*(1-gamma)+nu.1*gamma
v.new=gamma
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.mbs.dist.Q1.0, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
m.new = optimize(prof.MBS, c(0.0001,0.9999),
u = u.new, y = y.new, v = v.new, a.0=a.0, a.1=a.1, b.0=b.0, b.1=b.1, p.0=p.0, p.1=p.1, nu.0=nu.0, nu.1=nu.1)$minimum
dif = max(abs(c(rho.last, lambda.last, m.last) -
c(rho.new, lambda.new, m.new)))
m.last = m.new
rho.last = rho.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
v.last = v.new
i = i + 1
}
m.last=m.new=max(m.new, 0.001)
aux.se = hessian(observed.llike.0.mbs.dist, x0 = c(rho.last,
lambda.last, m.last), t = t, delta = delta,
ind = cluster)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.mbs.dist(c(rho.last,
lambda.last, m.last), t = t, delta = delta,
ind = cluster)
para = c(rho.new, lambda.new, Var_U(m.new))
se[length(se)] = se[length(se)]*abs(d.Var_U(m.new))
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull", data=data)
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last,
dist)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta.last), cluster)
gamma=m.last^2/(2*m.last^2+1-m.last)
a.0=2*A+(m.last^2/(1-m.last)+1)/(2*m.last)
a.1=2*A+(m.last^2/(1-m.last)+1)/(2*(m.last+1))
b.0=(m.last)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
b.1=(m.last+1)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
p.0=r+1/2
p.1=r+1/2-1
nu.0=besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.0)*(a.0/b.0)^(-p.0/2)*m.last^(1/2)+besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last))
nu.1=besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.1)*(a.0/b.0)^(-p.1/2)*m.last^(1+1/2)+besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last))
u.new=E.g(m.last, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
y.new=nu.0*(1-gamma)+nu.1*gamma
v.new=gamma
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.mbs.dist.Q1, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
m.new = optimize(prof.MBS, c(0.0001,0.9999),
u = u.new, y = y.new, v = v.new, a.0=a.0, a.1=a.1, b.0=b.0, b.1=b.1, p.0=p.0, p.1=p.1, nu.0=nu.0, nu.1=nu.1)$minimum
dif = max(abs(c(beta.last, rho.last, lambda.last,
m.last) - c(beta.new, rho.new, lambda.new,
m.new)))
beta.last = beta.new
m.last = m.new
rho.last = rho.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
v.last = v.new
i = i + 1
}
aux.se = hessian(observed.llike.mbs.dist, x0 = c(beta.last,
rho.last, lambda.last, m.last), t = t,
x = x, delta = delta, ind = cluster)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs(d.Var_U(m.new))
llike.obs = -observed.llike.mbs.dist(c(beta.last,
rho.last, lambda.last, m.last), t = t,
x = x, delta = delta, ind = cluster)
para = c(beta.new, rho.new, lambda.new, Var_U(m.new))
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
m.last=m.new=max(m.new, 0.001)
object.out <- list(coefficients = para, se = se,
z = u.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "MBS"
object.out$tau <- tau.MBS(m.last)
object.out$logLik <- llike.obs
}
if (dist == "pe" | dist == "exponential") {
dist.aux = "pe"
if (dist == "exponential") {
dist = "pe"
part = 0
dist.aux = "exponential"
}
observed.llike.0.mbs.dist <- function(eta, t, delta,
ind, part = NULL) {
lambda = eta[1:length(part)]
m = eta[length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
theta=m^2/(1-m)
p=r-1/2
a.m=2*A+(theta+1)/(2*m)
b.m=m*theta^2/(2*(theta+1))
a.m1=2*A+(theta+1)/(2*(m+1))
b.m1=(m+1)*theta^2/(2*(theta+1))
aux1=m*(sqrt(pi*m))^(-1)*(b.m/a.m)^(p/2)*(m*theta*(theta+1)^(-1)*besselK(sqrt(a.m*b.m),nu=p)+sqrt(b.m/a.m)*besselK(sqrt(a.m*b.m),nu=p+1))
aux2=(1-m)*(sqrt(pi*(m+1)))^(-1)*(b.m1/a.m1)^(p/2)*((m+1)*theta*(theta+1)^(-1)*besselK(sqrt(a.m1*b.m1),nu=p)+sqrt(b.m1/a.m1)*besselK(sqrt(a.m1*b.m1),nu=p+1))
P1 = (theta/2)-log(2)+(1/2)*log1p(theta)+log(aux1+aux2)
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mbs.dist <- function(eta, t, delta,
x, ind, part = NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x) + 1:length(part)]
m = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta), cluster)
theta=m^2/(1-m)
p=r-1/2
a.m=2*A+(theta+1)/(2*m)
b.m=m*theta^2/(2*(theta+1))
a.m1=2*A+(theta+1)/(2*(m+1))
b.m1=(m+1)*theta^2/(2*(theta+1))
aux1=m*(sqrt(pi*m))^(-1)*(b.m/a.m)^(p/2)*(m*theta*(theta+1)^(-1)*besselK(sqrt(a.m*b.m),nu=p)+sqrt(b.m/a.m)*besselK(sqrt(a.m*b.m),nu=p+1))
aux2=(1-m)*(sqrt(pi*(m+1)))^(-1)*(b.m1/a.m1)^(p/2)*((m+1)*theta*(theta+1)^(-1)*besselK(sqrt(a.m1*b.m1),nu=p)+sqrt(b.m1/a.m1)*besselK(sqrt(a.m1*b.m1),nu=p+1))
P1 = (theta/2)-log(2)+(1/2)*log1p(theta)+log(aux1+aux2)
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.mbs.dist.Q1.0 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
lambda = exp(eta[1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * log.lambda0 - z[ind] * Lambda0
-sum(P1)
}
observed.llike.mbs.dist.Q1 <- function(eta,
z, t, delta, x, ind, dist, part = NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x) + 1:length(part)])
Lambda0 = H.base(t, lambda = lambda, dist = dist,
part = part)
log.lambda0 = h.base(t, lambda = lambda, dist = dist,
part = part)
P1 = delta * x %*% beta + delta * log.lambda0 -
z[ind] * Lambda0 * exp(x %*% beta)
-sum(P1)
}
mm = length(t)
m.last = 0.5
u.last = rep(1, mm)
v.last = rep(0.5, mm)
y.last = rep(0.5, mm)
dif = 10
i = 1
lambda.last = rep(1/mean(t[which(delta==1)]),
length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last, cluster)
gamma=m.last^2/(2*m.last^2+1-m.last)
a.0=2*A+(m.last^2/(1-m.last)+1)/(2*m.last)
a.1=2*A+(m.last^2/(1-m.last)+1)/(2*(m.last+1))
b.0=(m.last)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
b.1=(m.last+1)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
p.0=r+1/2
p.1=r+1/2-1
nu.0=besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.0)*(a.0/b.0)^(-p.0/2)*m.last^(1/2)+besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last))
nu.1=besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.1)*(a.0/b.0)^(-p.1/2)*m.last^(1+1/2)+besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last))
u.new=E.g(m.last, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
y.new=nu.0*(1-gamma)+nu.1*gamma
v.new=gamma
aux.1 = optim(c(log(lambda.last)),
observed.llike.mbs.dist.Q1.0, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part=part)
lambda.new = exp(aux.1$par)
m.new = optimize(prof.MBS, c(0.0001,0.9999),
u = u.new, y = y.new, v = v.new, a.0=a.0, a.1=a.1, b.0=b.0, b.1=b.1, p.0=p.0, p.1=p.1, nu.0=nu.0, nu.1=nu.1)$minimum
dif = max(abs(c(lambda.last, m.last) -
c(lambda.new, m.new)))
m.last = m.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
v.last = v.new
i = i + 1
}
m.last=m.new=max(m.new, 0.001)
aux.se = hessian(observed.llike.0.mbs.dist, x0 = c(lambda.last,
m.last), t = t, delta = delta, ind = cluster,
part = part)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs(d.Var_U(m.new))
llike.obs = -observed.llike.0.mbs.dist(c(lambda.last,
m.last), t = t, delta = delta, ind = cluster,
part = part)
para = c(lambda.new, Var_U(m.new))
names(para) = names(se) = c(paste("lambda", 1:length(part),
sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(sum(delta)/sum(t),
length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda = lambda.last,
dist = dist, part = part)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.last*exp(x %*% beta.last), cluster)
gamma=m.last^2/(2*m.last^2+1-m.last)
a.0=2*A+(m.last^2/(1-m.last)+1)/(2*m.last)
a.1=2*A+(m.last^2/(1-m.last)+1)/(2*(m.last+1))
b.0=(m.last)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
b.1=(m.last+1)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
p.0=r+1/2
p.1=r+1/2-1
nu.0=besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.0)*(a.0/b.0)^(-p.0/2)*m.last^(1/2)+besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last))
nu.1=besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.1)*(a.0/b.0)^(-p.1/2)*m.last^(1+1/2)+besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last))
u.new=E.g(m.last, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
y.new=nu.0*(1-gamma)+nu.1*gamma
v.new=gamma
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.mbs.dist.Q1, method = "BFGS",
z = u.new, t = t, delta = delta,
x = x, ind = ind, dist = dist, part=part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
m.new = optimize(prof.MBS, c(0.0001,0.9999),
u = u.new, y = y.new, v = v.new, a.0=a.0, a.1=a.1, b.0=b.0, b.1=b.1, p.0=p.0, p.1=p.1, nu.0=nu.0, nu.1=nu.1)$minimum
dif = max(abs(c(beta.last, lambda.last,
m.last) - c(beta.new, lambda.new,
m.new)))
beta.last = beta.new
m.last = m.new
lambda.last = lambda.new
y.last = y.new
u.last = u.new
v.last = v.new
i = i + 1
}
m.last=m.new=max(m.new, 0.001)
aux.se = hessian(observed.llike.mbs.dist, x0 = c(beta.last,
lambda.last, m.last), t = t, x = x, delta = delta,
ind = cluster, part = part)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs(d.Var_U(m.new))
llike.obs = -observed.llike.mbs.dist(c(beta.last,
lambda.last, m.last), t = t, x = x, delta = delta,
ind = cluster, part = part)
para = c(beta.last, lambda.last, Var_U(m.last))
names(para) = names(se) = c(colnames(x), paste("lambda",
1:length(part), sep = ""), "theta")
if (dist.aux == "exponential")
names(para) = names(se) = c(colnames(x), "lambda",
"theta")
}
object.out <- list(coefficients = para, se = se,
z = u.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "MBS"
object.out$tau <- tau.MBS(m.last)
object.out$logLik <- llike.obs
if (dist.aux == "pe")
object.out$part <- part
}
if (dist == "np") {
observed.llike.0.mbs <- function(eta, t, delta, ind,
cox.aux) {
m = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0, cluster)
theta=m^2/(1-m)
p=r-1/2
a.m=2*A+(theta+1)/(2*m)
b.m=m*theta^2/(2*(theta+1))
a.m1=2*A+(theta+1)/(2*(m+1))
b.m1=(m+1)*theta^2/(2*(theta+1))
aux1=m*(sqrt(pi*m))^(-1)*(b.m/a.m)^(p/2)*(m*theta*(theta+1)^(-1)*besselK(sqrt(a.m*b.m),nu=p)+sqrt(b.m/a.m)*besselK(sqrt(a.m*b.m),nu=p+1))
aux2=(1-m)*(sqrt(pi*(m+1)))^(-1)*(b.m1/a.m1)^(p/2)*((m+1)*theta*(theta+1)^(-1)*besselK(sqrt(a.m1*b.m1),nu=p)+sqrt(b.m1/a.m1)*besselK(sqrt(a.m1*b.m1),nu=p+1))
P1 = (theta/2)-log(2)+(1/2)*log1p(theta)+log(aux1+aux2)
-sum(P1)
}
observed.llike.mbs <- function(eta, t, delta, x, ind,
cox.aux) {
m = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0*exp(x %*% beta), cluster)
theta=m^2/(1-m)
p=r-1/2
a.m=2*A+(theta+1)/(2*m)
b.m=m*theta^2/(2*(theta+1))
a.m1=2*A+(theta+1)/(2*(m+1))
b.m1=(m+1)*theta^2/(2*(theta+1))
aux1=m*(sqrt(pi*m))^(-1)*(b.m/a.m)^(p/2)*(m*theta*(theta+1)^(-1)*besselK(sqrt(a.m*b.m),nu=p)+sqrt(b.m/a.m)*besselK(sqrt(a.m*b.m),nu=p+1))
aux2=(1-m)*(sqrt(pi*(m+1)))^(-1)*(b.m1/a.m1)^(p/2)*((m+1)*theta*(theta+1)^(-1)*besselK(sqrt(a.m1*b.m1),nu=p)+sqrt(b.m1/a.m1)*besselK(sqrt(a.m1*b.m1),nu=p+1))
P1 = (theta/2)-log(2)+(1/2)*log1p(theta)+log(aux1+aux2)
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
mm = length(t)
m.last = 0.5
u.last = rep(1, mm)
v.last = rep(0.5, mm)
y.last = rep(0.5, mm)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(u.last[cluster])))
Lambda0.new = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.new, cluster)
gamma=m.last^2/(2*m.last^2+1-m.last)
a.0=2*A+(m.last^2/(1-m.last)+1)/(2*m.last)
a.1=2*A+(m.last^2/(1-m.last)+1)/(2*(m.last+1))
b.0=(m.last)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
b.1=(m.last+1)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
p.0=r+1/2
p.1=r+1/2-1
nu.0=besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.0)*(a.0/b.0)^(-p.0/2)*m.last^(1/2)+besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last))
nu.1=besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.1)*(a.0/b.0)^(-p.1/2)*m.last^(1+1/2)+besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last))
u.new=E.g(m.last, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
y.new=nu.0*(1-gamma)+nu.1*gamma
v.new=gamma
m.new = optimize(prof.MBS, c(0.0001,0.9999),
u = u.new, y = y.new, v = v.new, a.0=a.0, a.1=a.1, b.0=b.0, b.1=b.1, p.0=p.0, p.1=p.1, nu.0=nu.0, nu.1=nu.1)$minimum
dif = max(abs(m.last - m.new))
m.last = m.new
y.last = y.new
u.last = u.new
v.last = v.new
i = i + 1
}
m.last=m.new=max(m.new, 0.001)
aux.se = hessian(observed.llike.0.mbs, x0 = c(m.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs(d.Var_U(m.new))
para = c(Var_U(m.new))
names(para) = names(se) = c("theta")
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(u.last[cluster])))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
A = fail.cluster(Lambda0.new*exp(x %*% beta.last), cluster)
gamma=m.last^2/(2*m.last^2+1-m.last)
a.0=2*A+(m.last^2/(1-m.last)+1)/(2*m.last)
a.1=2*A+(m.last^2/(1-m.last)+1)/(2*(m.last+1))
b.0=(m.last)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
b.1=(m.last+1)*(m.last^2/(1-m.last))^2/(2*(m.last^2/(1-m.last)+1))
p.0=r+1/2
p.1=r+1/2-1
nu.0=besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.0)*(a.0/b.0)^(-p.0/2)*m.last^(1/2)+besselK(sqrt(a.1*b.1), nu=p.0)*(a.1/b.1)^(-p.0/2)*(m.last+1)^(-1/2)*(1-m.last))
nu.1=besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last)/(besselK(sqrt(a.0*b.0), nu=p.1)*(a.0/b.0)^(-p.1/2)*m.last^(1+1/2)+besselK(sqrt(a.1*b.1), nu=p.1)*(a.1/b.1)^(-p.1/2)*(m.last+1)^(1-1/2)*(1-m.last))
u.new=E.g(m.last, type=1, a.0, a.1, b.0, b.1, p.0, p.1, nu.0, nu.1, gamma)
y.new=nu.0*(1-gamma)+nu.1*gamma
v.new=gamma
m.new = optimize(prof.MBS, c(0.0001,0.9999),
u = u.new, y = y.new, v = v.new, a.0=a.0, a.1=a.1, b.0=b.0, b.1=b.1, p.0=p.0, p.1=p.1, nu.0=nu.0, nu.1=nu.1)$minimum
dif = max(abs(c(beta.last, m.last) - c(beta.new, m.new)))
beta.last = beta.new
m.last = m.new
y.last = y.new
u.last = u.new
v.last = v.new
i = i + 1
}
m.last=m.new=max(m.new, 0.001)
aux.se = hessian(observed.llike.mbs, x0 = c(beta.last,
m.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)]*abs(d.Var_U(m.new))
para = c(beta.new, Var_U(m.new))
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z = u.new)
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "MBS"
object.out$tau <- tau.MBS(m.last)
}
object.out
}
#--------------------------------#
# ADD GE FRAILTY MODEL #
#--------------------------------#
frailtyGE <- function(formula, data, dist = "np", prec = 1e-04,
max.iter = 1000, part=NULL) {
fail.cluster <- function(delta, indice) {
sum.fail <- function(ind, delta) {
sum(delta[which(indice == ind)])
}
unlist(lapply(1:max(indice), sum.fail, delta = delta))
}
vari.ge <- function(alpha){ (trigamma(1)-trigamma(alpha+1))/(digamma(alpha+1)-digamma(1))^2 }
dg <- function(x) {-psigamma(x+1,2)/(digamma(x+1)-digamma(1))^2-2*trigamma(x+1)*(trigamma(1)-trigamma(x+1))/(digamma(x+1)-digamma(1))^3} #derivative of theta in relacion to alpha
tau.GE.aux <- function(x, theta) {
eq.vari.ge=function(alpha, theta){vari.ge(alpha)-theta}
alpha=uniroot(eq.vari.ge, lower=0.00001, upper=40, theta=theta, extendInt="yes")$root
y = digamma(alpha+1)-digamma(1)
L = exp(lgamma(alpha+1) + lgamma(x/y+1) - lgamma(alpha+x/y+1))
aux = trigamma(x/y+1) - trigamma(alpha+x/y+1) + (digamma(x/y+1)-digamma(alpha+x/y+1))^2
x*L^2*aux/y^2
}
L.deriv<-function(alpha, s, n){
y = digamma(alpha+1)-digamma(1)
C <- function(x,i) {R = c(); for (j in 1:length(x)) R[j] = factorial(i-1)/(factorial(x[j]-1)*factorial(i-x[j])); R}
g <- function(x) {RR=c(); for(j in 1:length(x)) RR[j] = psigamma(x[j]-1,s/y+1)-psigamma(x[j]-1,alpha+s/y+1); RR}
L_i = c(); L0 = exp(lgamma(alpha+1)+lgamma(s/y+1)-lgamma(alpha+s/y+1))
if(n==0){return(L0)}
L_i[1] = L0*g(1)/y
if (n>1)
{
for (i in 2:n) L_i[i] = sum(c(C(1:(i-1),i)*L_i[(i-1):1]*g(1:(i-1))/y^(1:(i-1)), L0*g(i)/y^(i)))
}
L_i[n]
}
C.k.aux<-function(x, k, ri, bi, alpha){
exp((ri+k)*log(x)-x*bi+(alpha-1)*log(1-exp(-x)))}
C.k<-function(k, ri, bi, alpha){integrate(C.k.aux, lower=0, upper=Inf, k=k,
ri=ri, bi=bi, alpha=alpha, abs.tol=1e-15)$value}
C.k2.aux<-function(x, ri, bi, alpha){
log(1-exp(-x))*exp(ri*log(x)-x*bi+(alpha-1)*log(1-exp(-x)))}
C.k2<-function(ri, bi, alpha){integrate(C.k2.aux, lower=0, upper=Inf,
ri=ri, bi=bi, alpha=alpha, abs.tol=1e-15)$value}
tau.GE=function(theta) 4*integrate(tau.GE.aux,lower=0,upper=Inf,theta=theta, stop.on.error = FALSE)$value-1
#----------------------------------------------------
# BASELINE WEIBULL
#----------------------------------------------------
if (dist == "weibull") {
observed.llike.0.ge.dist <- function(eta, t, delta,
ind, dist, part=NULL) {
rho = eta[1]
lambda = eta[2]
alpha = eta[3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0, ind)
P1 = log((-1)^r*mapply(L.deriv, s=b, alpha=alpha, n=r))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ge.dist <- function(eta, t, delta,
x, ind, dist, part=NULL) {
beta = eta[1:ncol(x)]
rho = eta[ncol(x) + 1]
lambda = eta[ncol(x) + 2]
alpha = eta[ncol(x) + 3]
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0 * exp(x %*% beta), ind)
P1 = log((-1)^r*mapply(L.deriv, s=b, alpha=alpha, n=r))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ge.dist.Q1.0 <- function(eta, alpha,
w1, t, delta, x, ind, dist, part=NULL) {
rho = exp(eta[1])
lambda = exp(eta[2])
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * log.lambda0 - w1[ind] * Lambda0 / gamma
-sum(P1)
}
observed.llike.ge.dist.Q2.0 <- function(eta, rho=1, lambda,
w1, w2, t, delta, x, ind, dist, part=NULL) {
alpha = eta
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda, rho, dist)
P1 = - w1[ind] * Lambda0 / gamma - delta * log(gamma)
P2 = log(alpha) + alpha * w2
-(sum(P1)+sum(P2))
}
observed.llike.ge.dist.Q1 <- function(eta, alpha,
w1, t, delta, x, ind, dist, part=NULL) {
beta = eta[1:ncol(x)]
rho = exp(eta[ncol(x) + 1])
lambda = exp(eta[ncol(x) + 2])
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda, rho, dist)
log.lambda0 = h.base(t, lambda, rho, dist)
P1 = delta * x %*% beta + delta * log.lambda0 -
w1[ind] * Lambda0 * exp(x %*% beta) / gamma
-sum(P1)
}
observed.llike.ge.dist.Q2 <- function(eta, beta, rho=1, lambda,
w1, w2, t, delta, x, ind, dist, part=NULL) {
alpha = eta
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda, rho, dist)
P1 = - w1[ind] * Lambda0 * exp(x %*% beta) / gamma - delta*log(gamma)
P2 = log(alpha) + alpha * w2
-(sum(P1)+sum(P2))
}
m = length(t)
alpha.last = 2
theta.last = vari.ge(alpha.last)
w1.last = rep(1, m)
w2.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
rho.last = 1
lambda.last = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last, dist)
H0=fail.cluster(Lambda0.last, cluster)
bi = 1+H0/(digamma(alpha.last+1)-digamma(1))
C0=mapply(C.k, k=0, ri=r, bi=bi, alpha=alpha.last)
C1=mapply(C.k, k=1, ri=r, bi=bi, alpha=alpha.last)
C2=mapply(C.k2, ri=r, bi=bi, alpha=alpha.last)
w1.new = C1/C0
w2.new = C2/C0
aux.1 = optim(c(log(rho.last), log(lambda.last)),
observed.llike.ge.dist.Q1.0, method = "BFGS",
alpha = alpha.last, w1 = w1.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
rho.new = exp(aux.1$par[1])
lambda.new = exp(aux.1$par[2])
alpha.new = optim(alpha.last, observed.llike.ge.dist.Q2.0, w1 = w1.new, w2 = w2.new,
rho=rho.new, lambda=lambda.new, t=t, delta=delta, ind=ind, dist=dist,
method = "Brent", lower = 0, upper = min(10,alpha.last+3))$par
theta.new = vari.ge(alpha.new)
dif = max(abs(c(rho.last, lambda.last, theta.last) -
c(rho.new, lambda.new, theta.new)))
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
alpha.last = alpha.new
w1.last = w1.new
w2.last = w2.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ge.dist, x0 = c(rho.last,
lambda.last, alpha.last), t = t, delta = delta,
ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(alpha.last))
llike.obs = -observed.llike.0.ge.dist(c(rho.last,
lambda.last, alpha.last), t = t, delta = delta,
ind = cluster, dist = dist)
para = c(rho.new, lambda.new, theta.new)
names(para) = names(se) = c("rho", "lambda",
"theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = exp(-coef(cox.aux)[1]/cox.aux$scale)
rho.last = 1/cox.aux$scale
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda.last, rho.last, dist)
H0 = fail.cluster(Lambda0.last*exp(x%*%beta.last), cluster)
bi = 1+H0/(digamma(alpha.last+1)-digamma(1))
C0=mapply(C.k, k=0, ri=r, bi=bi, alpha=alpha.last)
C1=mapply(C.k, k=1, ri=r, bi=bi, alpha=alpha.last)
C2=mapply(C.k2, ri=r, bi=bi, alpha=alpha.last)
w1.new = C1/C0
w2.new = C2/C0
aux.1 = optim(c(beta.last, log(rho.last), log(lambda.last)),
observed.llike.ge.dist.Q1, method = "BFGS",
alpha = alpha.last, w1 = w1.new, t = t, delta = delta,
x = x, ind = ind, dist = dist)
beta.new = aux.1$par[1:ncol(x)]
rho.new = exp(aux.1$par[ncol(x) + 1])
lambda.new = exp(aux.1$par[ncol(x) + 2])
alpha.new = optim(alpha.last, observed.llike.ge.dist.Q2, w1 = w1.new, w2 = w2.new,
rho=rho.new, lambda=lambda.new, beta=beta.new, t=t, delta=delta, x=x, ind=ind, dist=dist,
method = "Brent", lower = 0, upper = min(10,alpha.last+3))$par
theta.new = vari.ge(alpha.new)
dif = max(abs(c(beta.last, rho.last, lambda.last, theta.last) -
c(beta.new, rho.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
rho.last = rho.new
lambda.last = lambda.new
alpha.last = alpha.new
w1.last = w1.new
w2.last = w2.new
i = i + 1
}
aux.se = hessian(observed.llike.ge.dist, x0 = c(beta.last,
rho.last, lambda.last, alpha.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(alpha.last))
llike.obs = -observed.llike.ge.dist(c(beta.last,
rho.last, lambda.last, alpha.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist)
para = c(beta.new, rho.new, lambda.new, theta.new)
names(para) = names(se) = c(colnames(x), "rho",
"lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z1 = w1.new/(digamma(alpha.new+1)-digamma(1)))
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "GE"
object.out$tau <- tau.GE(theta.last)
object.out$logLik <- llike.obs
}
#----------------------------------------------------
# BASELINE PE O EXP
#----------------------------------------------------
if (dist == "pe" | dist=="exponential") {
dist.aux="pe"
if(dist=="exponential"){dist="pe"; part=0; dist.aux="exponential"}
observed.llike.0.ge.dist <- function(eta, t, delta,
ind, dist, part=NULL) {
lambda = eta[1:length(part)]
alpha = eta[length(part)+1]
Lambda0 = H.base(t, lambda=lambda, dist=dist, part=part)
log.lambda0 = h.base(t, lambda=lambda, dist=dist, part=part)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0, ind)
P1 = log((-1)^r*mapply(L.deriv, s=b, alpha=alpha, n=r))
P2 = delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ge.dist <- function(eta, t, delta,
x, ind, dist, part=NULL) {
beta = eta[1:ncol(x)]
lambda = eta[ncol(x)+1:length(part)]
alpha = eta[ncol(x) + length(part) + 1]
Lambda0 = H.base(t, lambda=lambda, dist=dist, part=part)
log.lambda0 = h.base(t, lambda=lambda, dist=dist, part=part)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0 * exp(x %*% beta), ind)
P1 = log((-1)^r*mapply(L.deriv, s=b, alpha=alpha, n=r))
P2 = delta * x %*% beta + delta * log.lambda0
-(sum(P1) + sum(P2))
}
observed.llike.ge.dist.Q1.0 <- function(eta, alpha,
w1, t, delta, x, ind, dist, part=NULL) {
lambda = exp(eta[1:length(part)])
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda=lambda, dist=dist, part=part)
log.lambda0 = h.base(t, lambda=lambda, dist=dist, part=part)
P1 = delta * log.lambda0 - w1[ind] * Lambda0 / gamma
-sum(P1)
}
observed.llike.ge.dist.Q2.0 <- function(eta, rho=1, lambda,
w1, w2, t, delta, x, ind, dist, part=NULL) {
alpha = eta
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda, rho, dist)
P1 = - w1[ind] * Lambda0 / gamma - delta * log(gamma)
P2 = log(alpha) + alpha * w2
-(sum(P1)+sum(P2))
}
observed.llike.ge.dist.Q1 <- function(eta, alpha,
w1, t, delta, x, ind, dist, part=NULL) {
beta = eta[1:ncol(x)]
lambda = exp(eta[ncol(x)+1:length(part)])
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda=lambda, dist=dist, part=part)
log.lambda0 = h.base(t, lambda=lambda, dist=dist, part=part)
P1 = delta * x %*% beta + delta * log.lambda0 -
w1[ind] * Lambda0 * exp(x %*% beta) / gamma
-sum(P1)
}
observed.llike.ge.dist.Q2 <- function(eta, beta, rho=1, lambda,
w1, w2, t, delta, x, ind, dist, part=NULL) {
alpha = eta
gamma = digamma(alpha+1)-digamma(1)
Lambda0 = H.base(t, lambda, rho, dist)
P1 = - w1[ind] * Lambda0 * exp(x %*% beta) / gamma - delta*log(gamma)
P2 = log(alpha) + alpha * w2
-(sum(P1)+sum(P2))
}
m = length(t)
alpha.last=2
theta.last = vari.ge(alpha.last)
w1.last = rep(1, m)
w2.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
lambda.last = rep(1/mean(t), length(part))
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda=lambda.last, dist=dist, part=part)
H0=fail.cluster(Lambda0.last, cluster)
bi = 1+H0/(digamma(alpha.last+1)-digamma(1))
C0=mapply(C.k, k=0, ri=r, bi=bi, alpha=alpha.last)
C1=mapply(C.k, k=1, ri=r, bi=bi, alpha=alpha.last)
C2=mapply(C.k2, ri=r, bi=bi, alpha=alpha.last)
w1.new = C1/C0
w2.new = C2/C0
aux.1 = optim(log(lambda.last),
observed.llike.ge.dist.Q1.0, method = "BFGS",
w1 = w1.new, t = t, delta = delta, alpha = alpha.last,
x = x, ind = ind, dist = dist, part = part)
lambda.new = exp(aux.1$par)
alpha.new = optim(alpha.last, observed.llike.ge.dist.Q2.0, w1 = w1.new, w2 = w2.new,
method = "Brent", lower = 0, upper = min(10,alpha.last+3))$par
theta.new = vari.ge(alpha.new)
dif = max(abs(c(lambda.last, theta.last) -
c(lambda.new, theta.new)))
theta.last = theta.new
lambda.last = lambda.new
alpha.last = alpha.new
w1.last = w1.new
w2.last = w2.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ge.dist, x0 = c(lambda.last,
alpha.last), t = t, delta = delta,
ind = cluster, dist = dist, part=part)
se = sqrt(diag(solve(aux.se)))
llike.obs = -observed.llike.0.ge.dist(c(lambda.last, alpha.last), t = t, delta = delta,
ind = cluster, dist = dist, part=part)
se[length(se)] = se[length(se)] * abs(dg(alpha.last))
para = c(lambda.new, theta.new)
names(para) = names(se) = c(paste("lambda",1:length(part),sep=""), "theta")
if(dist.aux=="exponential") names(para) = names(se) = c("lambda", "theta")
}
if (ncol(x) > 0) {
cox.aux = survreg(Surv(t, delta) ~ x, dist = "weibull")
beta.last = -coef(cox.aux)[-1]/cox.aux$scale
lambda.last = rep(exp(-coef(cox.aux)[1]/cox.aux$scale), length(part))
while (i <= max.iter & dif > prec) {
Lambda0.last = H.base(t, lambda=lambda.last, dist=dist, part=part)
H0=fail.cluster(Lambda0.last*exp(x%*%beta.last), cluster)
bi = 1+H0/(digamma(alpha.last+1)-digamma(1))
C0=mapply(C.k, k=0, ri=r, bi=bi, alpha=alpha.last)
C1=mapply(C.k, k=1, ri=r, bi=bi, alpha=alpha.last)
C2=mapply(C.k2, ri=r, bi=bi, alpha=alpha.last)
w1.new = C1/C0
w2.new = C2/C0
aux.1 = optim(c(beta.last, log(lambda.last)),
observed.llike.ge.dist.Q1, method = "BFGS",
w1 = w1.new, t = t, delta = delta, alpha = alpha.last,
x = x, ind = ind, dist = dist, part=part)
beta.new = aux.1$par[1:ncol(x)]
lambda.new = exp(aux.1$par[ncol(x) + 1:length(part)])
alpha.new = optim(alpha.last, observed.llike.ge.dist.Q2, w1 = w1.new, w2 = w2.new,
lambda=lambda.new, beta=beta.new, t=t, delta=delta, x=x, ind=ind, dist=dist, part=part,
method = "Brent", lower = 0, upper = min(10,alpha.last+3))$par
theta.new = vari.ge(alpha.new)
dif = max(abs(c(beta.last, lambda.last, theta.last) -
c(beta.new, lambda.new, theta.new)))
beta.last = beta.new
theta.last = theta.new
lambda.last = lambda.new
alpha.last = alpha.new
w1.last = w1.new
w2.last = w2.new
i = i + 1
}
aux.se = hessian(observed.llike.ge.dist, x0 = c(beta.last,
lambda.last, alpha.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist, part=part)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(alpha.last))
llike.obs = -observed.llike.ge.dist(c(beta.last,
lambda.last, alpha.last), t = t,
x = x, delta = delta, ind = cluster, dist = dist, part=part)
para = c(beta.last, lambda.last, theta.last)
names(para) = names(se) = c(colnames(x), paste("lambda",1:length(part),sep=""), "theta")
if(dist.aux=="exponential") names(para) = names(se) = c(colnames(x), "lambda", "theta")
}
object.out <- list(coefficients = para, se = se,
z1 = w1.new/(digamma(alpha.new+1)-digamma(1)))
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$x <- x
object.out$dist <- dist.aux
object.out$dist.frail <- "GE"
object.out$tau <- tau.GE(theta.last)
object.out$logLik <- llike.obs
if(dist.aux=="pe") object.out$part <- part
}
#----------------------------------------------------
# BASELINE NP
#----------------------------------------------------
if (dist == "np") {
observed.llike.0.ge <- function(eta, t, delta, ind,
cox.aux) {
alpha = eta
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0, ind)
P1 = log((-1)^r*mapply(L.deriv, s=b, alpha=alpha, n=r))
-sum(P1)
}
observed.llike.ge <- function(eta, t, delta, x, ind,
cox.aux) {
alpha = eta[length(eta)]
beta = eta[-length(eta)]
Lambda0 = cumhazard.basal(t, cox.aux)
r = fail.cluster(delta, ind)
b = fail.cluster(Lambda0 * exp(x %*% beta), ind)
P1 = log((-1)^r*mapply(L.deriv, s=b, alpha=alpha, n=r))
P2 = delta * x %*% beta
-(sum(P1) + sum(P2))
}
cumhazard.basal <- function(t, coxph.object) {
ind.min <- function(t0, time) {
min(which(time >= t0))
}
bb = basehaz(coxph.object)
tt = bb$time
bb$hazard[unlist(lapply(t, ind.min, time = tt))]
}
observed.llike.ge.Q2.0 <- function(eta, Lambda0,
w1, w2, t, delta, x, ind, dist, part=NULL) {
alpha = eta
gamma = digamma(alpha+1)-digamma(1)
P1 = - w1[ind] * Lambda0 / gamma - delta * log(gamma)
P2 = log(alpha) + alpha * w2
-(sum(P1)+sum(P2))
}
observed.llike.ge.Q2 <- function(eta, beta, Lambda0,
w1, w2, t, delta, x, ind, dist, part=NULL) {
alpha = eta
gamma = digamma(alpha+1)-digamma(1)
P1 = - w1[ind] * Lambda0 * exp(x %*% beta) / gamma - delta*log(gamma)
P2 = log(alpha) + alpha * w2
-(sum(P1)+sum(P2))
}
m = length(t)
alpha.last=2
theta.last = vari.ge(alpha.last)
w1.last = rep(1, m)
w2.last = rep(1, m)
r = fail.cluster(delta, cluster)
dif = 10
i = 1
if (ncol(x) == 0) {
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ offset(log(w1.last[cluster]/(digamma(alpha.last+1)-digamma(1)))))
Lambda0.new = cumhazard.basal(t, cox.aux)
H0 = fail.cluster(Lambda0.new, cluster)
bi = 1+H0/(digamma(alpha.last+1)-digamma(1))
C0=mapply(C.k, k=0, ri=r, bi=bi, alpha=alpha.last)
C1=mapply(C.k, k=1, ri=r, bi=bi, alpha=alpha.last)
C2=mapply(C.k2, ri=r, bi=bi, alpha=alpha.last)
w1.new = C1/C0
w2.new = C2/C0
alpha.new = optim(alpha.last, observed.llike.ge.Q2.0, w1 = w1.new, w2 = w2.new,
Lambda0=Lambda0.new, t=t, delta=delta, ind=ind, dist=dist,
method = "Brent", lower = 0, upper = min(10,alpha.last+3))$par
theta.new=vari.ge(alpha.new)
dif = abs(theta.last - theta.new)
theta.last = theta.new
alpha.last = alpha.new
w1.last = w1.new
w2.last = w2.new
i = i + 1
}
aux.se = hessian(observed.llike.0.ge, x0 = c(alpha.last),
t = t, delta = delta, ind = cluster, cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(alpha.last))
para = c(theta.new)
names(para) = names(se) = "theta"
}
if (ncol(x) > 0) {
cox.aux = coxph(Surv(t, delta) ~ x)
beta.last = coef(cox.aux)
while (i <= max.iter & dif > prec) {
cox.aux = coxph(Surv(t, delta) ~ x + offset(log(w1.last[cluster]/(digamma(alpha.last+1)-digamma(1)))))
beta.new = coef(cox.aux)
Lambda0.new = cumhazard.basal(t, cox.aux)
H0 = fail.cluster(Lambda0.new*exp(x%*%beta.last), cluster)
bi = 1+H0/(digamma(alpha.last+1)-digamma(1))
C0=mapply(C.k, k=0, ri=r, bi=bi, alpha=alpha.last)
C1=mapply(C.k, k=1, ri=r, bi=bi, alpha=alpha.last)
C2=mapply(C.k2, ri=r, bi=bi, alpha=alpha.last)
w1.new = C1/C0
w2.new = C2/C0
alpha.new = optim(alpha.last, observed.llike.ge.Q2, w1 = w1.new, w2 = w2.new,
beta=beta.new, Lambda0=Lambda0.new, t=t, delta=delta, ind=ind, dist=dist, x=x,
method = "Brent", lower = 0, upper = min(10,alpha.last+3))$par
theta.new=vari.ge(alpha.new)
dif = max(abs(c(beta.last, theta.last) - c(beta.new,
theta.new)))
beta.last = beta.new
theta.last = theta.new
alpha.last=alpha.new
w1.last = w1.new
w2.last = w2.new
i = i + 1
}
aux.se = hessian(observed.llike.ge, x0 = c(beta.last,
alpha.last), t = t, delta = delta, x = x, ind = cluster,
cox.aux = cox.aux)
se = sqrt(diag(solve(aux.se)))
se[length(se)] = se[length(se)] * abs(dg(alpha.last))
para = c(beta.new, theta.new)
names(para) = names(se) = c(colnames(x), "theta")
}
bb = basehaz(cox.aux)
Lambda0 = cbind(bb$time, bb$hazard)
colnames(Lambda0) = c("time", "hazard")
object.out <- list(coefficients = para, se = se,
z1 = w1.new/(digamma(alpha.new+1)-digamma(1)))
class(object.out) <- "extrafrail"
object.out$t <- t
object.out$delta <- delta
object.out$id <- cluster
object.out$Lambda0 <- Lambda0
object.out$x <- x
object.out$dist <- dist
object.out$dist.frail <- "GE"
object.out$tau <- tau.GE(theta.last)
}
object.out
}
if (dist.frail == "GA")
val <- frailtyGA(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "IG")
val <- frailtyIG(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "WL")
val <- frailtyWL(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "BS")
val <- frailtyBS(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "TN")
val <- frailtyTN(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "MIG")
val <- frailtyMIG(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "MBS")
val <- frailtyMBS(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
if (dist.frail == "GE")
val <- frailtyGE(formula, data, dist = dist, prec = prec,
max.iter = max.iter, part = part)
val
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.