Nothing
aftgomp <-
function (X, Y, strata, offset, init, control, center)
{
if (!is.matrix(X))
X <- matrix(X, ncol = 1)
nn <- NROW(X)
ncov <- NCOL(X)
if (NROW(Y) != nn)
stop("Y NROW error")
if (NCOL(Y) != 3)
stop("Y NCOL error", )
if (ncov) {
wts <- Y[, 2] - Y[, 1]
means <- apply(X, 2, weighted.mean, w = wts)
means <- apply(X, 2, mean)
for (i in 1:ncov) {
X[, i] <- X[, i] - means[i]
}
}
if (missing(strata) || is.null(strata)) {
strata <- rep(1, nn)
ns <- 1
}
else {
strata <- as.integer(factor(strata))
ns <- max(strata)
}
if (length(strata) != nn)
stop("Error in stratum variable")
if (missing(offset) || is.null(offset))
offset <- rep(0, nn)
if (missing(init) || is.null(init))
init <- rep(0, ncov)
if (length(init) != ncov)
stop("Error in init")
printlevel <- control$trace
iter <- control$maxiter
nstra <- c(0, cumsum(table(strata)))
bdim <- ncov + 2 * ns
dGomp <- function(beta) {
enter <- Y[, 1]
exit <- Y[, 2]
event <- Y[, 3]
if (ncov) {
bz <- offset + X %*% beta[1:ncov]
}
else {
bz <- offset
}
grad <- numeric(ncov + 2 * ns)
for (i in 1:ns) {
T0 <- enter[strata == i]
T <- exit[strata == i]
D <- event[strata == i]
z <- X[strata == i, , drop = FALSE]
ezb <- exp(bz[strata == i])
alpha <- beta[ncov + 2 * i]
gamma <- beta[ncov + 2 * i - 1]
scale <- exp(gamma) * ezb
shape <- exp(alpha)
emgamma <- exp(-gamma)
ezbH <- ##ezb *
(Hgompertz(T, scale = scale, shape = shape,
param = "canonical") - Hgompertz(T0, scale = scale,
shape = shape, param = "canonical"))
grad[ncov + 2 * i] <- sum(D) - sum(ezbH)
grad[ncov + 2 * i - 1] <- -sum(D * T) * emgamma -
sum(D) + sum(ezb * (T * hgompertz(T, scale = scale,
shape = shape, param = "canonical") - T0 * hgompertz(T0,
scale = scale, shape = shape, param = "canonical")))
if (ncov) {
for (j in 1:ncov) {
grad[j] <- grad[j] + sum(D * z[, j]) - sum(ezbH *
z[, j])
}
}
}
grad
}
Fmin <- function(beta) {
total <- 0
for (i in 1:ns) {
scale <- exp(beta[ncov + 2 * i - 1])
shape <- exp(beta[ncov + 2 * i])
if (ncov) {
bz <- offset[strata == i] + X[strata == i, ,
drop = FALSE] %*% beta[1:ncov]
}
else {
bz <- offset[strata == i]
}
ebz <- exp(bz)
H1 <- Hgompertz(Y[strata == i, 2], shape = shape,
scale = scale, param = "canonical")
H0 <- Hgompertz(Y[strata == i, 1], shape = shape,
scale = scale, param = "canonical")
h <- hgompertz(Y[strata == i, 2], shape = shape,
scale = scale, log = TRUE, param = "canonical")
ret1 <- sum(Y[strata == i, 3] * (h + bz))
ret2 <- sum(ebz * (H1 - H0))
total <- total + ret1 - ret2
}
return(total)
}
beta0 <- numeric(2 * ns)
ncov.save <- ncov
ncov <- 0
for (i in 1:ns) {
enter <- Y[strata == i, 1]
exit <- Y[strata == i, 2]
event <- Y[strata == i, 3]
score <- offset[strata == i]
beta0[(2 * i - 1):(2 * i)] <- gompstartCanonical(enter,
exit, event, score)
}
res0 <- optim(beta0, Fmin, gr = dGomp, method = "BFGS", control = list(fnscale = -1,
reltol = 1e-10), hessian = FALSE)
ncov <- ncov.save
beta <- numeric(bdim)
if (ncov)
beta[1:ncov] <- init
for (i in 1:ns) {
enter <- Y[strata == i, 1]
exit <- Y[strata == i, 2]
event <- Y[strata == i, 3]
score <- offset[strata == i] + X[strata == i, , drop = FALSE] %*%
init
beta[(ncov + 2 * i - 1):(ncov + 2 * i)] <- gompstartCanonical(enter,
exit, event, score)
}
res <- optim(beta, Fmin, gr = dGomp, method = "BFGS", control = list(fnscale = -1,
reltol = 1e-10), hessian = TRUE)
if (res$convergence != 0)
stop("[gompreg]: No convergence")
coefficients <- res$par
coef.names <- colnames(X)
if (ns > 1) {
for (i in 1:ns) {
coef.names <- c(coef.names, paste("log(scale)", as.character(i),
sep = ":"), paste("log(shape)", as.character(i),
sep = ":"))
}
}
else {
coef.names <- c(coef.names, "log(scale)", "log(shape)")
}
names(coefficients) <- coef.names
fit <- list(coefficients = coefficients, loglik = c(res0$value,
res$value))
fit$gradient <- dGomp(coefficients)
fit$pfixed <- FALSE
fit$var <- tryCatch(solve(-res$hessian), error = function(e) e)
fit$hessian <- res$hessian
if (is.matrix(fit$hessian)) {
colnames(fit$hessian) <- coef.names
rownames(fit$hessian) <- coef.names
}
if (ncov) {
dxy <- diag(2 * ns + ncov)
for (i in 1:ns) {
row <- ncov + 2 * i
shape.corr <- sum(means * fit$coefficients[1:ncov])
fit$coefficients[row] <- fit$coefficients[row] -
shape.corr
dxy[row, 1:ncov] <- means
}
if (is.numeric(fit$var))
fit$var <- dxy %*% fit$var %*% t(dxy)
fit$dxy <- dxy
}
if (is.matrix(fit$var)) {
colnames(fit$var) <- coef.names
rownames(fit$var) <- coef.names
}
fit$n.strata <- ns
fit$df <- ncov
fit$fail <- FALSE
fit$param <- "canonical"
fit
}
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.