"recurrent"<-
function(MAT, tau, Xc = NULL, eps = 0.0001, times, mtime, cutpts = NULL)
{
#****************************************************************************
# Function name: recurrent
# ------------------------
# Written by: Melanie Wigg, Oct 1998  (University of Waterloo)
# Revised by: Melanie Wigg, Jan 1999
#
# This function is based on the methods discussed in the paper "Some
# Simple Robust Methods for the Analysis of Recurrent Events" by Lawless
# and Nadeau (Technometrics, May 1995, Vol. 37, No. 2).  
#
# See README file for details on how to use this function and interpret
# its results.
#****************************************************************************
# input variables:
#	MAT    = matrix of event times for each subject
#	tau    = vector of censoring times for each subject
#	Xc     = covariate values for each subject
#	eps    = parameter for the criterion for terminating Newton's Method
#	times  = vector of times defining the CMF covariance matrix 
#                   (only used if no covariates)
#	mtime  = single value "t" in the covariance matrix of (B, Mo(t))
#                   (only used if covariates given)
#	cutpts = vector of cut points defining intervals for piecewise
#		    constant intensity functions (all cut points, including
#		    lower and upper bounds, should be given; thus, a vector
#		    of length r+1 will result in r intervals)
#****************************************************************************
# output variables:  (if no covariates)
#	dft    = vector of distinct failure times
#	CMF    = estimated cumulative mean function at values in dft
#	times  = "times" vector defined by user in function call
#	Vr     = robust covariance matrix of CMF at specified "times"
#	Vp     = covariance matrix of CMF under the Poisson process assumption
#
# output variables:  (if covariates given)
#	B      = estimated beta vector
#	dft    = vector of distinct failure times
#	m0     = estimated baseline mean function at values in dft
#	M0     = estimated baseline cumulative mean function at values in dft
#	mtime  = "mtime" value defined by user in function call
#	Vr     = robust covariance matrix of (B, M0(t))
#	Vp     = covariance matrix under the Poisson process assumption
#	count  = number of iterations of Newton's Method until convergence
#****************************************************************************
# temporarily used variables:
#	k      = number of subjects
#	m      = number of columns of MAT
#	r      = number of distinct failure times
#	Event  = k x r matrix; Event(i,s) is the number of events occurring
#                at time s for subject i
#	idft   = indicates the column of each distinct failure time in the
#                Event matrix; e.g. if dft[2] = 16, then idft[17] = 2
#		 (idft is offset by 1 to allow for failures at time 0)
#	Delt   = k x r matrix; Delt(i,s) indicates whether subject i is 
#                observed (1) or not (0) at time s
#	nd     = vector; number of events at each distinct failure time
#	ds     = vector; number of subjects observed at each distinct
#                failure time
#	mf     = mean function, defined at each distinct failure time
#	dtimes = vector; dtimes[j] = max(i; dft[i] <= times[j])
#	maxt   = largest value in dtimes
#	mfd    = estimated variance of estimated CMF at each distinct
#                failure time
#	p      = number of covariates
#	check  = value compared to "eps" at the end of each iteration
#                of Newton's Method
#	XB     = matrix; XB[i,s] is the value of the linear predictor
#                for subject i at time s
#	R      = R(t;B) defined at each distinct failure time (see 3.4 in paper)
#	m0     = mean function at each distinct failure time (see 3.4)	
#	dg     = array with dim(r,p,k); dg[j,l,i] is the derivative
#                of g_i(j) with respect to B_l (see 3.3)
#	dR     = matrix; dR[j,l] is the derivative of R(j;B) w.r.t. B_l
#	W      = array with dim(r,p,k); W[j,l,i] is the value of
#                W_i(B,j) for B_l  (see 3.7)
#	U      = value of estimating equation for B at each B_l (see 3.3, 3.5)	
#	A1     = matrix, p x p  (see 3.9)
#	B.old  = old value of B used to compute new value
#	tx     = integer; tx = max(i; dft[i] <= mtime)
#	A2, A, C, B1i, B11, B12, B21, B22, BB  = building blocks for
#                calculating cov(B, Mo(t))  (see appendix A.2)
#****************************************************************************
#----------------------------------------------------------------------------
# create Delt and Event matrices 
#----------------------------------------------------------------------------
	k <- nrow(MAT)	
	if(k != length(tau)) {
		stop("How many subjects?  nrow(MAT) must equal length(tau)")
	}
	m <- ncol(MAT)
	dft <- sort(unique(as.vector(MAT)))
	if(sum(dft < 0) > 0) {
		stop("All values in MAT must be greater than or equal to 0")
	}
	if(sum(tau < 0) > 0) {
		stop("All values in tau must be greater than or equal to 0")
	}
	if(is.null(cutpts)) {
		if(sum(trunc(dft) != dft) > 0) {
			stop("All values in MAT must be integers when piecewise constant intensity model is not used")
		}
		r <- length(dft)
		Event <- matrix(0,k,r)
		idft <- numeric(max(dft)+1)
		for (i in 1:r) {
			idft[dft[i]+1] <- i
		}
		for (i in 1:k) {
			for (j in 1:m) {
				s <- MAT[i,j]
				u <- idft[s+1]
				Event[i,u] <- Event[i,u] + 1
			}
		}
		Delt <- t(apply(as.matrix(tau),1,makeDelt,dft=dft))
	}
	else {
		if(is.null(Xc)) {
			stop("Piecewise constant intensity model only allowed if covariates are specified")
		}
		cutpts <- sort(cutpts)
		if(cutpts[1] < 0) {
			stop("All cutpts values must be greater than or equal to 0")
		}
		r <- length(cutpts) - 1
		if(sum(dft < cutpts[1] | dft > cutpts[r+1]) > 0) {
			stop("Some event times do not fall within the intervals defined by cutpts")
		}
		if(sum(tau < cutpts[1] | tau > cutpts[r+1]) > 0) {
			stop("Some censoring times do not fall within the intervals defined by cutpts")
		}
		if(length(unique(cut(dft,breaks=cutpts,include.lowest=T))) != r) {
			stop("Some intervals defined by cutpts do not contain any events")
		}
		Delt <- t(apply(as.matrix(tau),1,makeDelt2,cutpts=cutpts))
		Event <- t(apply(MAT,1,makeEvent,cutpts=cutpts,r=r))
		Delt.v <- as.vector(Delt)
		Event <- as.vector(Event)
		Event <- ifelse(Delt.v == 0, 0, Event/Delt.v)
		Event <- matrix(Event,k,r)
		dft <- cutpts[-1]
		if(length(mtime) > 1) {
			stop("mtime must be a single time value")
		}
		if(sum(dft == mtime) == 0) {
			stop("mtime must be the upper bound of one of the intervals defined by cutpts")
		}
	}
	nd <- apply(Event*Delt, 2, sum)
#
#----------------------------------------------------------------------------
# if no covariates are specified, compute CMF and the robust covariance
# matrix of CMF at specified "times"
#----------------------------------------------------------------------------
	if(is.null(Xc)) {
		if(sum(times < min(dft)) > 0) {
			stop(paste("Some elements in times are less than",min(dft),"(the first distinct failure time)"))
		}
		ds <- apply(Delt, 2, sum)
		mf <- nd/ds
		CMF <- cumsum(mf)
		times <- sort(times)
		dtimes <- numeric(length(times))
		for (i in 1:length(times)) {
			dtimes[i] <- sum(dft <= times[i])
		}
		maxt <- max(dtimes)
		V <- apply(as.matrix(1:k), 1, Vhat, Delt = Delt[,1:maxt], Event = Event[,1:maxt], ds = ds[1:maxt], mf = mf[1:maxt], dtimes = dtimes, k=k, maxt=maxt) 
		V <- matrix(V,length(times),k)
		Vr <- V %*% t(V)
		Vp <- matrix(0,length(times),length(times))
		mfd <- cumsum(mf/ds)	
		for(i in 1:length(times)) {
			for(j in i:length(times)) {
				Vp[i,j] <- Vp[j,i] <- mfd[dtimes[i]]
			}
		}	
		return(dft, CMF, times, Vr, Vp)
	}
	else {
#
#----------------------------------------------------------------------------
# estimates based on regression model
#----------------------------------------------------------------------------
		if(is.vector(Xc)) {
			Xc <- as.matrix(Xc)
		}
		if(is.matrix(Xc)) {
			p <- ncol(Xc)
			Xc.new <- array(0,c(r,p,nrow(Xc)))
			Xc.new <- aperm(aperm(Xc.new,c(3,2,1)) + as.vector(Xc),c(3,2,1))
			Xc <- Xc.new
		}
		if(length(dim(Xc)) != 3) {
			stop("Xc must be a vector, matrix, or 3-dimensional array")
		}
		if(k != dim(Xc)[3]) {
			stop("MAT and Xc do not agree on number of subjects")
		}
		if(r != dim(Xc)[1]) {
			stop("MAT and Xc do not agree on number of distinct failure times (or intervals)")
		}
		if(length(mtime) > 1) {
			stop("mtime must be a single time value")
		}
		if(mtime < min(dft)) {
			stop(paste("mtime is less than",min(dft),"(the first distinct failure time)"))
		}
		p <- dim(Xc)[2]
		B <- rep(0,p)      # initial values of beta
		count <- 0         # count number of iterations
		check <- 1000
		while(check > eps) {
			fit <- getbeta(count, Xc, B, k, r, p, Delt, nd, Event)
			B <- fit$B
			count <- fit$count
			check <- fit$check
		}
#
#----------------------------------------------------------------------------
# calculate estimate of Mo(t) based on the estimated beta
#----------------------------------------------------------------------------
		XB.t <- aperm(aperm(Xc,c(2,1,3)) * B, c(3,2,1))
		XB <- matrix(0,k,r)
		for(i in 1:p) {
			XB <- XB + XB.t[,,i]
		}
		R <- diag(t(Delt) %*% exp(XB)) 
		m0 <- nd/R
		M0 <- cumsum(m0)
#
#----------------------------------------------------------------------------
# calculate estimated covariance matrix of (beta, Mo(t))
#----------------------------------------------------------------------------
		dg <- array(0,c(r,p,k))
		for (i in 1:p) {
			dg[,i,] <- t(exp(XB)) * Xc[,i,]
		}
		dR <- matrix(apply(aperm(aperm(dg,c(3,1,2)) * as.vector(Delt),c(2,3,1)),c(1,2),sum),r,p)
		W <- array(as.vector(apply(as.matrix(1:k),1,makeW,Xc=Xc,dR=dR,R=R)),c(r,p,k))
		A1 <- matrix(apply(as.matrix(1:k),1,makeA1,dg=dg,Delt=Delt,m0=m0,W=W),p^2,k)
		A1 <- matrix(apply(A1,1,sum),p,p)
		
		tx <- sum(dft <= mtime)	
		A2 <- matrix(m0[1:tx],1,tx) %*% matrix(dR[1:tx,]/R[1:tx],tx,p)
		A <- cbind(rbind(A1,A2),c(rep(0,p),1))

		C <- (1/R[1:tx]) %*% (t(Delt[,1:tx]) * (t(Event[,1:tx]) - (t(exp(XB[,1:tx])) * m0[1:tx])))

		B1i <- t(matrix(apply(as.matrix(1:k),1,makeB1i,Delt=Delt,Event=Event,m0=m0,XB=XB,W=W),p,k))
		B11 <- t(B1i) %*% B1i
		B21 <- C %*% B1i
		B12 <- t(B21)
		B22 <- C %*% t(C)
		BB <- cbind(rbind(B11,B21),rbind(B12,B22))
		Vr <- solve(A) %*% BB %*% t(solve(A))
		temp <- sum(m0[1:tx]/R[1:tx]) + (A2 %*% solve(A1) %*% t(A2))
		Vp <- cbind(rbind(solve(A1),A2),rbind(t(A2),temp))
		return(B,dft,m0,M0,mtime,Vr,Vp,count)
	}
}
"Vhat"<-
function(i, Delt, Event, ds, mf, dtimes, k, maxt)
{
	Delt <- matrix(Delt,k,maxt)
	Event <- matrix(Event,k,maxt)
	temp <- (Delt[i,  ]/ds) * (Event[i,  ] - mf)
	temp2 <- cumsum(temp)
	temp2[dtimes]
}
"makeDelt"<-
function(taui, dft)
{
	as.numeric(dft <= taui)
}
"makeDelt2"<-
function(taui, cutpts)
{
	Deltrow <- as.numeric(cutpts <= taui)
	last <- sum(Deltrow)
	if(last != length(cutpts)) {
		Deltrow[last+1] <- (taui - cutpts[last]) / (cutpts[last+1] - cutpts[last])
	}
	Deltrow[-1]
}
"makeEvent"<-
function(events, cutpts, r)
{
# events are sorted to remove any NA values
	events <- sort(events)
	group <- cut(events, breaks=cutpts, include.lowest=T)
	Eventrow <- numeric(r)
	for(i in 1:length(group)) {
		Eventrow[group[i]] <- Eventrow[group[i]] + 1
	}
	Eventrow
}
"makeW"<-
function(i, Xc, dR, R)
{
        Xc[,  , i] - (dR/R)
}
"makeU"<-
function(i, Delt, Event, W)
{
        (Delt[i,  ] * Event[i,  ]) %*% W[,  , i]
}
"makeA1"<-
function(i, dg, Delt, m0, W)
{
	t(dg[,,i] * (Delt[i,] * m0)) %*% W[,,i]
}
"makeB1i"<-
function(i, Delt, Event, m0, XB, W)
{
	(Delt[i,] * (Event[i,] - (m0 * exp(XB[i,])))) %*% W[,,i]
}
"getbeta"<-
function(count, Xc, B, k, r, p, Delt, nd, Event)
{
	count <- count + 1
	XB.t <- aperm(aperm(Xc,c(2,1,3)) * B, c(3,2,1))
	XB <- matrix(0,k,r)
	for(i in 1:p) {
		XB <- XB + XB.t[,,i]
	}
	R <- diag(t(Delt) %*% exp(XB)) 
	m0 <- nd/R
	dg <- array(0,c(r,p,k))
	for (i in 1:p) {
		dg[,i,] <- t(exp(XB)) * Xc[,i,]
	}
	dR <- matrix(apply(aperm(aperm(dg,c(3,1,2)) * as.vector(Delt),c(2,3,1)),c(1,2),sum),r,p)
	W <- array(as.vector(apply(as.matrix(1:k),1,makeW,Xc=Xc,dR=dR,R=R)),c(r,p,k))
	U <- apply(matrix(apply(as.matrix(1:k),1,makeU,Delt=Delt,Event=Event,W=W),p,k),1,sum)

	A1 <- matrix(apply(as.matrix(1:k),1,makeA1,dg=dg,Delt=Delt,m0=m0,W=W),p^2,k)
	A1 <- matrix(apply(A1,1,sum),p,p)
#
#----------------------------------------------------------------------------
# get an updated beta via Newton's Method
#----------------------------------------------------------------------------
	B.old <- B
	B <- B.old + (solve(A1) %*% U)
	B <- as.vector(B)
#
#----------------------------------------------------------------------------
# check if iterations should continue
#----------------------------------------------------------------------------
	check <- sum(abs((B-B.old)/B))
	return(B, count, check)
}