#' Fit a zero-truncated recurrent events model
#'
#' \code{fitrem} fits a zero-truncated recurrent events model on the data
#' provided and calculates the point and interval estimate of the population
#' size
#'
#' This is a specialised function, as it requires the data to be in a
#' three-dimensional format (persons x time x variables). The default
#' optimization procedure is a Newton-Raphson method, while numerical
#' optimization via \code{optim} is also possible
#'
#'
#' @param data An array of at least three dimensions
#' @param bstart A vector of starting values
#' @param ... Further arguments to be passed from or to other methods
#' @return A coefficient matrix, point and interval estimate of the population
#' size, loglikelihood and number of iterations
#' @export
fitrem <- function(data, bstart, ...) {
if(missing(data)) stop("No data supplied")
events <- data[, , 1]
X <- data[, , -1, drop = F]
p <- dim(X)[3]
if (missing(bstart)) {
bfirst <- bstart <- rep(0, p)
} else if (length(bstart) < p) {
warning("Starting values not provided for all parameters, using 0 as start value for parameters not specified")
bfirst <- c(bstart, rep(0, p - length(bstart)))
} else if (length(bstart) > p) {
warning("Too many starting values provided, discarding excess values")
bfirst <- bstart[1:p]
} else bfirst <- bstart
time <- dim(X)[2]
n <- dim(X)[1]
b <- matrix(bfirst, p, 1)
iter <- 0
tol <- 1
while (tol > 1e-07) {
u <- exp(tensor::tensor(X, b, alongA = 3, alongB = 1))
dim(u) <- c(n, time)
U <- c(rowSums(u))
p0 <- exp(-U)
pCap <- 1 - p0
eU <- exp(U)
logl <- sum(log(u[events == 1])) - sum(U) - sum(log(1 - p0))
upCap <- (u * p0)/pCap
EupCap <- (events - u - upCap)
db <- tensor::tensor(X, EupCap, alongA = c(1, 2), alongB = c(1, 2))
XuX <- tensor::tensor(X, c(u) * X, alongA = c(1, 2), alongB = c(1, 2))
ddb1 <- tensor::tensor(X * c((u)/(eU - 1)), X, alongA = c(1, 2), alongB = c(1, 2))
uX <- apply(c(u) * X, c(1, 3), sum)
ddb2 <- colSums(tensor::tensor(X * c(u * eU)/(eU - 1)^2, uX, alongA = c(1), alongB = c(1)))
ddb <- -XuX - ddb1 + ddb2
infor <- solve(ddb)
bnew <- b - infor %*% db
tol <- sum(abs(bnew - b))
b <- bnew
iter <- iter + 1
print(logl)
}
se <- sqrt(diag(-infor))
tvalue <- b/se
pvalue <- 2 * pnorm(-abs(tvalue))
nb <- length(b)
coef <- matrix(c(b, se, tvalue, pvalue), nrow = nb, ncol = 4)
colnames(coef) <- c("B", "SE", "t-value", "P-value")
rownames(coef) <- dimnames(data)[[3]][-1]
Nhat <- sum(1/pCap)
if (p != 1) {
var1 <- colSums((uX * p0)/pCap^2)
} else {
var1 <- sum((uX * p0)/pCap^2)
}
var1 <- -(t(var1) %*% infor %*% var1)
var2 <- sum(p0/pCap^2)
Nse <- sqrt(var1 + var2)
CI <- cbind(Nhat - 1.96*Nse, Nhat + 1.96*Nse)
AIC <- -2 * logl + 2 * length(b)
Rij <- cumsum(colSums(u))
output <- list(coef = coef, logl = logl, AIC = AIC, iter = iter, Nhat = Nhat, CI = CI,
Nse = Nse, infor = infor)
class(output) <- c("trem", "list")
output
}
#' Print results of \code{fitrem}
#'
#' \code{print.trem} is a S3 method to print the results generated by \code{fitrem}
#' @param x object of class \code{trem}
#' @param ... Further arguments to be passed from or to other methods
#' @export
print.trem <- function(x, ...){
printCoefmat(x$coef, has.Pvalue = T)
cat("\n", "Nhat (95%-CI): ", x$Nhat, " (" ,x$CI[1]," - ", x$CI[2], ")", sep = "")
}
#' Print a summary of the results from \code{fitrem}
#'
#' \code{summary.trem} is a S3 method that prints more information than \code{print.trem}
#' @param object object of class \code{trem}
#' @param ... Further arguments to be passed from or to other methods
#' @export
summary.trem <- function(object, ...){
printCoefmat(object$coef, has.Pvalue = T)
cat("\n", "Nhat (95%-CI) ", object$Nhat, " (" ,object$CI[1]," - ", object$CI[2], ")", sep = "")
cat("\n\n", "Loglikelihood:", object$logl)
cat("\n\n", "AIC:", object$AIC)
cat("\n\n", "Estimation converged in", object$iter, "iterations")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.