.packageName <- "event"
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
# autointensity(times, window=NULL, maxlag=total/10,
#	ylab="Conditional probability", xlab="Lag",
#	main="Autointensity function", xlim=c(0,max(times)),
#	ylim=c(0,if(plotse)max(se1)else max(z$density)), lty=1,
#	plot=TRUE, plotse=TRUE, add=FALSE, ...)
#
#  DESCRIPTION
#
#    Plots the autointensity function of a point process

autointensity <- function(times, window=NULL, maxlag=max(times),
	ylab="Conditional probability", xlab="Lag",
	main="Autointensity function", xlim=c(0,max(times)),
	ylim=c(0,if(plotse)max(se1)else max(z$density)), lty=1,
	plot=TRUE, plotse=TRUE, add=FALSE, ...){
if(any(times<=0))stop("Times must be strictly positive")
n <- length(times)
total <- sum(times)
if(is.null(window))window <- total/n/2
#
# calculate histogram and standard errors
#
tmp <- cumsum(times)
times <- tmp[2:n]-tmp[1]
for(i in 3:n)times <- c(times,tmp[i:n]-tmp[i-1])
times <- times[times<maxlag]
z <- hist(times,plot=FALSE,breaks=seq(0,max(times+window),by=window))
z$density <- z$counts*total/(total-z$mids)/window/n
se1 <- ifelse(z$density>0,(sqrt(z$density)+2/sqrt(window*n))^2,0)
se2 <- sqrt(z$density)-2/sqrt(window*n)
se2 <- ifelse(z$density>0,ifelse(se2>0,se2,0)^2,0)
#
# plot if required
#
if(plot){
	if(add)plot(z$mids,z$density,type="l",lty=lty)
	else plot(z$mids,z$density,type="l",main=main,ylab=ylab,xlab=xlab,
		xlim=xlim,ylim=ylim,lty=lty,...)
	lines(xlim,rep(n/total,2))
	if(plotse){
		lines(z$mids,se1,lty=3)
		lines(z$mids,se2,lty=3)}}
invisible(list(mids=z$mids,density=z$density,se=rbind(se1,se2)))}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#	coxre(response, censor, nest=NULL, cov=NULL, stratified=FALSE,
#	  cumul=FALSE, estimate=1, iter=10, print.level=0, ndigit=10,
#	  gradtol=0.00001, steptol=0.00001, iterlim=100, fscale=1,
#	  typsiz=abs(estimate), stepmax=estimate)
#
#  DESCRIPTION
#
#    Cox model with frailites.
# glim code for coxre translated from Clayton, D. (1987) The analysis of
# event history data: a review of progress and outstanding problems.
# Statistics in Medicine 7: 819-841

coxre <- function(response, censor, nest=NULL, cov=NULL, stratified=FALSE,
	cumul=FALSE,estimate=1, iter=10, print.level=0, ndigit=10,
	gradtol=0.00001, steptol=0.00001, iterlim=100, fscale=1,
	typsiz=abs(estimate), stepmax=estimate){
#
# set up likelihood function to iterative call glm: SLOW!!
#
like <- function(g){
	g <- exp(g)
	l <- 0
	for(i in 1:iter){
		fexp <- NULL
		for(ij in split(fit,nnest))fexp <- c(fexp,sum(ij))
		fexp <- fexp/rand
		ffex <- NULL
		for(ij in split(fit,interval))ffex <- c(ffex,sum(ij))
		ffex <- ffex/fixd
		tmp <- l
		l <- sum(counts*log(fit))-sum((fobs+1/g)*log(1+g*fexp)+fobs*log(rand))
		rand <- (1+g*fobs)/(1+g*fexp)
		fixd <- ffob/ffex
		oset <- log(ntimes*rand[nnest]*fixd[interval])
		lmod <- if(!fcov){if(stratified)counts~1 else counts~resp}
			else {if(stratified)counts~cov else counts~cov+resp}
		z <- glm(lmod,family=poisson,offset=oset)
		fit <- z$fit
		coef <- z$coef
		rm(z)
		if((l-tmp)^2<0.00002)break}
	l <- l+sum(lgamma(1/g+fobs)-lgamma(1/g)+fobs*log(g))
	names(fixd) <- 1:length(fixd)
	names(rand) <- 1:length(rand)
	list(like=-l,
		fixed=fixd,
		random=rand,
		inthaz=ffex,
		coefficients=coef)}
#
# likelihood for nlm
#
likel <- function(g){
	g <- exp(g)
	l <- 0
	for(i in 1:iter){
		fexp <- NULL
		for(ij in split(fit,nnest))fexp <- c(fexp,sum(ij))
		fexp <- fexp/rand
		ffex <- NULL
		for(ij in split(fit,interval))ffex <- c(ffex,sum(ij))
		ffex <- ffex/fixd
		tmp <- l
		l <- sum(counts*log(fit))-sum((fobs+1/g)*log(1+g*fexp)+fobs*log(rand))
		rand <- (1+g*fobs)/(1+g*fexp)
		fixd <- ffob/ffex
		oset <- log(ntimes*rand[nnest]*fixd[interval])
		lmod <- if(!fcov){if(stratified)counts~1 else counts~resp}
			else {if(stratified)counts~cov else counts~cov+resp}
		z <- glm(lmod,family=poisson,offset=oset)
		fit <- z$fit
		rm(z)
		if((l-tmp)^2<0.00002)break}
	l <- l+sum(lgamma(1/g+fobs)-lgamma(1/g)+fobs*log(g))
	-l}
call <- sys.call()
if(estimate<=0)stop("estimate must be positive")
#
# from response and covariate can tell what model to fit
#
fcov <- !missing(cov)
if(is.vector(response,mode="numeric")){
	nind1 <- length(response)
	nc <- 1
	# check censoring and
	# transform response and censor to 1 column matrices
	if(is.null(censor))censor <- matrix(1,nrow=nind1,ncol=1)
	else if(is.vector(censor,mode="numeric")&&
		length(censor)==length(response))
		censor <- matrix(censor,ncol=1)
	else stop("censor must be a vector the same length as response")
	response <- matrix(response,ncol=1)
	stratified <- TRUE}
else if(is.matrix(response)){
	nind1 <- dim(response)[1]
	# check censoring
	if(is.null(censor))censor <- matrix(1,ncol=dim(response)[2],nrow=nind1)
	else if(!is.matrix(censor)||dim(censor)[2]!=dim(response)[2])
		stop("response and censor must have the same number of columns")
	# if not stratified model, create nesting indicator and
	# transform response and censor to 1 column matrices
	if(!stratified){
		resp1 <- gl(dim(response)[2],dim(response)[1],length(response))
		nest <- as.integer(rep(nest,dim(response)[2]))
		if(fcov)cov <- rep(cov,dim(response)[2])
		response <- matrix(response,ncol=1)
		censor <- matrix(censor,ncol=1)
		nc <- 1}
	else nc <- dim(response)[2]}
else stop("response must be a vector or matrix")
#
# if cumulative times, transform to times between events
#
if(cumul)response <- cbind(response[,1],response[,2:dim(response)[2]]-response[,1:(dim(response)[2]-1)])
#
# check data all agree in size
#
nind <- dim(response)[1]
if(dim(censor)[1]!=nind)
	stop("response and censor must have the same number of rows")
if(!is.null(nest)&&length(nest)!=nind)
	stop(paste("nest must have length",nind1))
if(!is.null(cov)&&length(cov)!=nind)stop(paste("cov must have length",nind1))
#
# find unique times taking into account censoring
#
tot <- 0
ntimes <- interval <- nnest <- counts <- event <- ncov <- resp <- NULL
for(k in 1:dim(response)[2]){
	tt <- rep(0,sum(censor[,k]))
	tt[cumsum(censor[,k])*censor[,k]] <- response[(1:nind)*censor[,k],k]
	o <- order(response[,k])
	cc <- censor[o,k]
	i <- nest[o]
	tt <- unique(sort(tt))
	tt1 <- sort(response[,k])
	tt1 <- c(cc[1],(cc[2:nind]*(tt1[2:nind]>tt1[1:(nind-1)])))
	tt2 <- c(tt[1],tt[2:length(tt)]-tt[1:(length(tt)-1)])
	tt <- rep(0,sum(tt2>0))
	tt[cumsum(tt2>0)*(tt2>0)] <- tt2[(1:length(tt2))*(tt2>0)]
	int <- tt[1]
	nint <- 1
	id <- i[1]
	nevent <- cc[1]
	nz <- 0
	if(fcov){
		tcov <- cov[o]
		ncov1 <- tcov[1]}
	if(!stratified){
		tresp <- resp1[o]
		nresp <- tresp[1]}
	for(j in 2:nind){
		nz <- nz+(tt1[j]==0)
		ev <- rep(0,j-nz)
		ev[j-nz] <- cc[j]
		nevent <- c(nevent,ev)
		id <- c(id,rep(i[j],j-nz))
		nint <- c(nint,1:(j-nz))
		int <- c(int,tt[1:(j-nz)])
		if(fcov)ncov1 <- c(ncov1,rep(tcov[j],j-nz))
		if(!stratified)nresp <- c(nresp,rep(tresp[j],j-nz))}
	ntimes <- c(ntimes,int)
	interval <- c(interval,nint+tot)
	nnest <- c(nnest,id)
	counts <- c(counts,nevent)
	event <- c(event,rep(k,length(nint)))
	if(fcov)ncov <- c(ncov,ncov1)
	if(!stratified)resp <- c(resp,nresp)
	tot <- max(interval)}
if(!stratified)resp <- as.factor(resp)
cov <- ncov
#
# set up formula to fit
#
if(!fcov){
	if(stratified)lmod <- counts~1
	else lmod <- counts~resp}
else {
	if(stratified)lmod <- counts~cov
	else lmod <- counts~cov+resp}
fobs <- capply(counts,nnest)
ffob <- capply(counts,interval)
#
# obtain initial estimates of random and fixed effects, and fitted values
# then call nlm
#
rand <- rep(1,max(nnest))
fixd <- rep(1,max(interval))
z <- glm(lmod,family=poisson)
fit <- z$fit
rm(z)
z1 <- nlm(likel,log(estimate), hessian=TRUE, print.level=print.level,
	typsiz=typsiz, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
z2 <- like(z1$est)
z3 <- list(maxlike=z2$like,
	aic=z2$like+length(z2$coef)+length(z2$fixed)+1,
	df=length(response)-length(z2$coef)-length(z2$fixed)-1,
	iterations=z1$iter,
	code=z1$code,
	call=call,
	gamma=exp(z1$est),
	fixed=z2$fixed,
	random=z2$random,
	inthaz=z2$inthaz,
	coefficients=z2$coef,
	stratified=stratified)
class(z3) <- "llrf"
z3}

### print method
###
print.llrf <- function(z){
	if(z$stratified)cat("Stratified ")
	cat("Cox proportional hazards model with gamma frailty\n")
	cat("\nCall:",deparse(z$call),sep="\n")
	cat("\n")
	if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
	cat("-Log likelihood   ",z$maxlike,"\n")
	cat("Degrees of freedom",z$df,"\n")
	cat("AIC               ",z$aic,"\n")
	cat("Iterations        ",z$iterations,"\n\n")
	tmp <- trigamma(1/z$gamma)
	cat("gamma =      ",z$gamma,"\n")
	cat("correlation =",tmp/(tmp+trigamma(1)),"\n")
	cat("\nRegression coefficients:\n")
	print(z$coef)
	cat("\nFixed effects:\n")
	print(z$fixed)
	cat("\nRandom effects:\n")
	print(z$random)}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#	cprocess(times=NULL, events=NULL, number=FALSE, lty=NULL,
#		xlim=NULL, ylim=NULL, xlab="Time", ylab="Counts", ...)
#
#  DESCRIPTION
#
#    Plots counting process data

cprocess <- function(times=NULL, events=NULL, number=FALSE, lty=NULL,
	xlim=NULL, ylim=NULL, xlab="Time", ylab="Counts", ...){
if(is.null(times)&&is.null(events))
	stop("Either times or events must be supplied")
#
# if events are supplied, create the times
#
if(!is.null(events)){
	if(is.matrix(events)){
		if(is.null(times))tim <- list()
		ev <- list()
		for(i in 1:dim(events)[1]){
			if(is.null(times))
				tim <- c(tim,list(rep(1,dim(events)[2])))
			ev <- c(ev,list(events[i,]))}
		times <- tim
		rm(tim)
		events <- ev
		rm(ev)}
	else if(is.vector(events,mode="numeric")){
		if(is.null(times))times <- list(rep(1,length(events)))
		events <- list(events)}
	else if(inherits(events,"response")){
		tmp <- list()
		for(i in 1:length(nobs(events)))
			tmp <- c(tmp,list(c(events$times[covind(events)==i][1],diff(events$times[covind(events)==i]))))
		times <- tmp
		tmp <- list()
		for(i in 1:length(nobs(events)))
			tmp <- c(tmp,list(events$y[covind(times)==i]))
		events <- tmp
		rm(tmp)}
	else if(inherits(events,"repeated")){
		tmp <- list()
		for(i in 1:length(nobs(events)))
			tmp <- c(tmp,list(c(events$response$times[covind(events)==i][1],diff(events$response$times[covind(events)==i]))))
		times <- tmp
		tmp <- list()
		for(i in 1:length(nobs(events)))
			tmp <- c(tmp,list(events$response$y[covind(events)==i]))
		events <- tmp
		rm(tmp)}
	else if(is.list(events)){
		if(is.null(times))for(i in events)
			times <- c(times,list(rep(1,length(i))))}
	else stop("events must be a matrix, vector, or list of vectors")}
#
# if times are supplied, transform them to a list
#
if(!is.list(times)){
	if(!is.vector(times,mode="numeric"))
		stop("times must be a vector or a list of vectors")
	times <- list(times)
	if(!is.null(events)){
		if(length(times)!=length(events))
			stop("numbers of individuals for events and times differ")
		for(i in 1:length(times))
		if(length(events[[i]])!=length(times[[i]]))
			stop(paste("individual ",i,"does not have the same number of events as times"))}}
else if(inherits(times,"response")){
	tmp <- list()
	for(i in 1:length(nobs(times)))
		tmp <- c(tmp,list(times$y[covind(times)==i]))
	times <- tmp
	rm(tmp)}
k <- xl <- yl <- 0
for(i in times){
	k <- k+1
	if(is.matrix(i))i <- i[,1]
	else if(!is.vector(i,mode="numeric"))
		stop("times must be a list of vectors or matrices")
	if(any(i<0))stop("negative times")
	jx <- length(i)
	if(is.null(events))jy <- jx
	else jy <- sum(events[[k]])
	tmp <- cumsum(i)
	yl <- max(jy,yl)
	xl <- max(tmp[jx],xl)}
if(missing(xlim))xlim <- c(0,xl)
if(missing(ylim))ylim <- c(0,yl)
#
# plot the event histories
#
z <- list()
k <- j <- 0
for(i in times){
	j <- j+1
	if(is.null(events))s <- 1:length(i)
	else s <- cumsum(events[[j]])
	time <- cumsum(i)
	# double the values to create steps
	n <- 2*length(s)-1
	count <- rep(0,n)
	count[seq(1,n,by=2)] <- s
	count[seq(2,n-1,by=2)] <- s[1:(length(s)-1)]
	count <- c(0,0,count)
	tmp <- rep(0,n)
	tmp[seq(1,n,by=2)] <- time
	tmp[seq(2,n-1,by=2)] <- time[1:(length(s)-1)]
	tmp <- c(0,tmp,tmp[length(tmp)])
	if(missing(lty))llty <- k%%4+1
	else llty <- lty
	if(k)lines(tmp,count,lty=llty)
	else plot(tmp,count,type="l",lty=llty,xlim=xlim,ylim=ylim,
		xlab=xlab,ylab=ylab,...)
	if(number)text(tmp[length(tmp)],count[length(count)]+1,paste(k+1))
	k <- k+1
	z <- c(z,list(cbind(tmp,count)))}
invisible(z)}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#	pp(y, censor=1)
#	ident(y, id)
#	tpast(y)
#	ttime(y, id)
#	bp(y, id, censor=1)
#	tccov(y, x, id)
#	tvcov(y, x, tx)
#	vdm(y, x, id=NULL, tx=NULL, factor=FALSE, time=FALSE)
#	ehr(point, lambda=NULL, linear=NULL, plambda=NULL, wt=1,
#		envir=parent.frame(), print.level=0,
#		typsiz=rep(1,length(plambda)), ndigit=10,
#		gradtol=0.00001, stepmax=max(10*sqrt(plambda%*%plambda),10),
#		steptol=0.0004, iterlim=100, fscale=1)
#
#  DESCRIPTION
#
#    Functions for setting up and fitting counting process models

### point process created from times (y) between events
###
pp <- function(y, censor=1){
#
# y must contain integers
#
if(!is.vector(y,mode="numeric"))stop("y must be a numeric vector")
if(min(y)<=0)stop("All times must be positive")
if(any(round(y)!=y))stop("Times must be integers")
if(any(censor!=0&&censor!=1))stop("Censor indicator must be zeros and ones")
if(length(censor)!=1&&length(censor)!=length(y))
	stop("Time and censor vectors must be the same length")
#
# create vector of zeros with ones at events
#
point <- rep(0, sum(y))
point[cumsum(y)] <- censor
point}

### individual identification vector
###
ident <- function(y, id){
if(min(y)<=0)stop("All times must be positive")
if(!is.vector(id,mode="numeric"))stop("id must be a numeric vector")
if(length(y)!=length(id))stop("Time and id vectors must be the same length")
if(length(unique(id))!=max(id)-min(id)+1||min(id)!=1)
	stop("id must be consecutive numbers")
rep(id, y)}

### time past since previous event
###
tpast <- function(y){
if(min(y)<=0)stop("All times must be positive")
unlist(lapply(as.list(y), seq))}

### total time elapsed for each individual
###
ttime <- function(y, id){
if(length(idd <- ident(y,id))==1)return(idd)
z <- capply(rep(1,length(idd)),idd,cumsum)
names(z) <- NULL
z}

### number of previous events for each individual, for birth processes
### one should be added, if process starts at an event
###
bp <- function(y, id, censor=1){
bp1 <- function(i) c(0,cumsum(i)[1:(length(i)-1)])
if(length(point <- pp(y, censor=censor))==1)return(point)
if(length(idd <- ident(y, id))==1)return(idd)
z <- capply(point, idd, bp1)
names(z) <- NULL
z}

### time-constant covariate - id must be numbered consecutively
### x has one value for each distinct id
###
tccov <- function(y, x, id){
if(length(y)!=length(id))stop("Time and id must be the same length")
if(length(x)!=length(unique(id)))
	stop("There must be one covariate value per individual")
if(length(idd <- ident(y, id))==1)return(idd)
x[idd]}

### time-varying covariate - tx gives the times at which x changes
### may also be used to create weight vector
###
tvcov <- function(y, x, tx){
if(min(y)<=0||min(tx)<0)stop("All times must be positive")
if(length(x)!=length(tx))
	stop("Covariate and time vectors must be the same length")
if(sum(y)!=sum(tx))stop("Total response time must equal total covariate time")
rep(x, tx)}

### design matrix
###
vdm <- function(y, x, id=NULL, tx=NULL, factor=FALSE, time=FALSE) {
if(time){
	if(length(xx <- tvcov(y, x, tx))==1)return(xx)}
else if(length(xx <- tccov(y, x, id))==1)return(xx)
if(factor)xx <- factor(xx)
wr(~xx)$design}

### fit an intensity function to event histories, where point is
### produced by point <- pp(y) and lambda is the log intensity function
###
ehr <- function(point, lambda=NULL, linear=NULL, plambda=NULL, delta=1,
	envir=parent.frame(), print.level=0, typsiz=rep(1,length(plambda)),
	ndigit=10, gradtol=0.00001, iterlim=100, fscale=1,
	stepmax=max(10*sqrt(plambda%*%plambda),10), steptol=0.0004){
call <- sys.call()
if(any(point<0))stop("Response vector must be non-negative integers")
n <- length(point)
npl <- length(plambda)
#
# set up deltas if necessary
#
dt <- any(delta>1)
if(dt){
	if(length(point)!=length(delta))stop("point and delta must be the same length")
	delta <- log(delta)}
#
# find linear part of each regression and save model for printing
#
if(inherits(linear,"formula")&&is.null(lambda)){
	lambda <- linear
	linear <- NULL}
if(inherits(linear,"formula")){
	linmodel <- if(!is.null(attr(finterp(linear,.envir=envir),"parameters")))
			attr(finterp(linear,.envir=envir),"model")}
else linmodel <- NULL
#
# check if linear contains W&R formula
#
if(inherits(linear,"formula")){
	tmp <- attributes(finterp(linear,.envir=envir))
	lf1 <- length(tmp$parameters)
	if(!is.character(tmp$model))stop("linear must be a W&R formula")
	else if(length(tmp$model)==1)stop("linear must contain covariates")
	rm(tmp)}
else lf1 <- 0
#
# transform formula to function and check number of parameters
#
if(inherits(lambda,"formula")){
	linarg <- if(lf1>0) "linear" else NULL
	lambda4 <- finterp(lambda,.envir=envir,.args=linarg)
	npt1 <- length(attr(lambda4,"parameters"))
	if(is.character(attr(lambda4,"model"))){
	# W&R formula
		if(length(attr(lambda4,"model"))==1){
		# intercept model
			lambda3 <- if(dt)
				function(p) lambda(p,p[lf1]*rep(1,n))+delta
			else function(p) lambda(p,p[lf1]*rep(1,n))
			attributes(lambda3) <- attributes(lambda4)}
		else {
			lambda3 <- lambda4
			rm(lambda4)}}
	else {
	# formula with unknowns
		if(npl!=npt1&&lf1==0){
			cat("\nParameters are ")
			cat(attr(lambda4,"parameters"),"\n")
			stop(paste("plambda should have",npt1,"estimates"))}
		if(dt){
			lambda3 <- if(lf1>0)
				function(p,linear) lambda4(p,linear)+delta
				else function(p) lambda4(p)+delta
			attributes(lambda3) <- attributes(lambda4)}
		else lambda3 <- lambda4
		if(is.list(plambda)){
			if(!is.null(names(plambda))){
				o <- match(attr(lambda4,"parameters"),names(plambda))
				plambda <- unlist(plambda)[o]
				if(sum(!is.na(o))!=length(plambda))stop("invalid estimates for lambda - probably wrong names")}
			else plambda <- unlist(plambda)}}}
else if(!is.function(lambda)){
	lambda3 <- if(dt)function(p) p[1]*rep(1,n)+delta
		else function(p) p[1]*rep(1,n)
	npt1 <- 1}
else {
	lambda4 <- fnenvir(lambda,.envir=envir)
	npt1 <- length(attr(lambda4,"parameters"))-(lf1>0)
	if(dt)lambda3 <- if(lf1>0)function(p,linear) lambda(p,linear)+delta
		else function(p) lambda(p)+delta
	else lambda3 <- lambda
	attributes(lambda3) <- attributes(lambda4)
	rm(lambda4)}
#
# if linear part, modify location function appropriately
#
if(lf1>0){
	if(is.character(attr(lambda3,"model")))
		stop("lambda cannot be a W&R formula if linear is supplied")
	dm1 <- wr(linear,data=envir)$design
	if(is.null(lambda2))lambda2 <- lambda3
	lambda1 <- function(p)lambda3(p,dm1%*%p[(npt1+1):(npt1+lf1)])}
else {
	if(lf1==0&&length(lambda3(plambda))==1){
		lambda1 <- function(p) lambda3(p)*rep(1,n)
		attributes(lambda1) <- attributes(lambda3)}
	else lambda1 <- lambda3}
#
# give appropriate attributes to lambda1 for printing
#
if(is.null(attr(lambda1,"parameters"))){
	attributes(lambda1) <- if(is.function(lambda)){
		if(!inherits(lambda,"formulafn"))
			attributes(fnenvir(lambda,.envir=envir))
		else attributes(lambda)}
		else attributes(fnenvir(lambda1))}
#
# check that correct number of estimates was supplied
#
nlp <- npt1+lf1
if(nlp!=npl)stop(paste("plambda should have",nlp,"initial estimates"))
#
# set up likelihood function, check that it returns appropriate values
# and call nlm
#
fn <- function(p) {
	l <- lambda1(p)
	sum(exp(l)-point*l)}
if(fscale==1)fscale <- fn(plambda)
if(is.na(fn(plambda)))
	stop("Likelihood returns NAs: probably invalid initial values")
z0 <- nlm(fn, p=plambda, hessian=TRUE, print.level=print.level, typsiz=typsiz,
	ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
if(any(point>1))z0$minimum <- z0$minimum+sum(lgamma(point+1))
#
# calculate se's
#
if(length(plambda)==1)cov <- 1/z0$hessian
else {
	a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
		else qr(z0$hessian)$rank
	if(a==length(plambda))cov <- solve(z0$hessian)
	else cov <- matrix(NA,ncol=length(plambda),nrow=length(plambda))}
se <- sqrt(diag(cov))
z1 <- list(
	call=call,
	intensity=lambda1,
	linear=linear,
	linmodel=linmodel,
	maxlike=z0$minimum,
	aic=z0$minimum+length(plambda),
	coefficients=z0$estimate,
	se=se,
	cov=cov,
	corr=cov/(se%o%se),
	gradient=z0$gradient,
	iterations=z0$iter,
	error=z0$error,
	code=z0$code)
class(z1) <- "intensity"
return(z1)}

### standard method

deviance.intensity <- function(z) 2*z$maxlike

### print method
###
print.intensity <- function(z) {
np <- length(z$coefficients)
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
cat("Log intensity function:\n")
if(!is.null(attr(z$intensity,"formula")))
	cat(deparse(attr(z$intensity,"formula")),sep="\n")
else if(!is.null(attr(z$intensity,"model"))){
	t <- deparse(attr(z$intensity,"model"))
	t[1] <- sub("expression\\(","",t[1])
	t[length(t)] <- sub("\\)$","",t[length(t)])
	cat(t,sep="\n")}
if(inherits(z$linear,"formulafn"))
	cat("Linear part: ",deparse(attr(z$linear,"formula")),sep="\n")
cat("\n-Log likelihood   ",z$maxlike,"\n")
cat("AIC               ",z$aic,"\n")
cat("Iterations        ",z$iterations,"\n\n")
cat("Coefficients:\n")
cname <- if(is.character(attr(z$intensity,"model")))attr(z$intensity,"model")
	else if(length(grep("linear",attr(z$intensity,"parameters")))>0)
	attr(z$intensity,"parameters")[grep("\\[",attr(z$intensity,"parameters"))]
	else attr(z$intensity,"parameters")
if(!is.null(z$linmodel))cname <- c(cname,z$linmodel)
coef.table <- cbind(z$coefficients, z$se)
dimnames(coef.table) <- list(cname, c("estimate", "se"))
print.default(coef.table, digits=4, print.gap=2)
if(np>1){
	cat("\nCorrelations:\n")
	dimnames(z$corr) <- list(seq(1,np),seq(1,np))
	print.default(z$corr, digits=4)}
invisible(z)}

# examples of linear log intensity functions
#exponential <- ~1
#Weibull <- ~log(time(y))
#extreme.value <- ~time(y)
#birth1 <- ~bp(y,id)
#birth2 <- ~log(1+bp(y,id))

# examples of nonlinear log intensity functions
#negative.binomial <- function(p) p[1]+log(p[2]+bp(y,id))
#gen.negative.binomial <- function(p) p[1]+p[3]*log(p[2]+bp(y,id))
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#	hboxcox(y,m,s,f)
#	hburr(y,m,s,f)
#	hcauchy(y,m,s)
#	hexp(y,m)
#	hgextval(y,s,m,f)
#	hgamma(y,s,m)
#	hggamma(y,s,m,f)
#	hhjorth(y,m,s,f)
#	hinvgauss(y,m,s)
#	hlaplace(y,m,s)
#	hlnorm(y,m,s)
#	hlogis(y,m,s)
#	hglogis(y,m,s,f)
#	hnorm(y,m,s)
#	hpareto(y,m,s)
#	hstudent(y,m,s,f)
#	hweibull(y,s,m)
#	hgweibull(y,s,m,f)
#	hskewlaplace(y,s,m,f)
#
#  DESCRIPTION
#
#    Functions for various log hazards or intensities

# for log distributions, subtract log(y) from the intensity function

### Box-Cox intensity
###
# f=1 gives truncated normal
hboxcox <- function(y,m,s,f) {
	y1 <- y^f/f
	-(y1-m)^2/s/2+(f-1)*log(y)-log(2*pi*s)/2-log(1-pnorm(y1,m,sqrt(s))+(f<0)*(1-pnorm(0,m,sqrt(s))))}

### Burr intensity
###
hburr <- function(y,m,s,f) {
	y1 <- y/m
	y2 <- y1^s
	log(f*s/m)+(s-1)*log(y1)-log(1+y2)}

### Cauchy intensity
###
hcauchy <- function(y,m,s) log(dcauchy(y,m,s))-log(1-pcauchy(y,m,s))

### exponential intensity
###
hexp <- function(y,rate)
	if(length(rate)==1)rep(log(rate),length(y)) else log(rate)

### generalized extreme value intensity
###
# f=1 gives truncated extreme value
hgextval <- function(y,s,m,f) {
	y1 <- y^f/f
	ey <-exp(y1)
	log(s)+s*(y1-log(m))-(ey/m)^s+(f-1)*log(y)-log(1-pweibull(ey,s,m)-(f<0)*exp(-m^-s))}

### gamma intensity
###
hgamma <- function(y,shape,rate=1,scale=1/rate)
	dgamma(y,shape,scale=scale,log=TRUE)-
		pgamma(y,shape,scale=scale,lower=FALSE,log=TRUE)

### generalized gamma intensity
###
hggamma <- function(y,s,m,f) {
	t <- m/s
	u <- t^f
	y1 <- y^f
	v <- s*f
	-v*log(t)-y1/u+log(f)+(v-1)*log(y)-lgamma(s)-
		log(1-pgamma(y1,s,scale=u))}

### Hjorth intensity
###
hhjorth <- function(y,m,s,f) log(y/m^2+f/(1+s*y))

### inverse Gaussian intensity
###
hinvgauss <- function(y,m,s) {
	t <- y/m
	v <- sqrt(y*s)
	-((t-1)^2/(y*s)+log(2*s*pi*y^3))/2-log(1-pnorm((t-1)/v)
		-exp(2/(m*s))*pnorm(-(t+1)/v))}
### Laplace intensity
###
hlaplace <- function(y,m=0,s=1){
	plp <- function(u){
		t <- exp(-abs(u))/2
		ifelse(u<0,t,1-t)}
	-abs(y-m)/s-log(2*s)-log(1-plp((y-m)/s))}

### log normal intensity
###
hlnorm <- function(y,m,s) dlnorm(y,m,s,TRUE)-plnorm(y,m,s,FALSE,TRUE)

### logistic intensity
###
hlogis <- function(y,m,s) dlogis(y,m,s,TRUE)-plogis(y,m,s,FALSE,TRUE)

### generalized logistic intensity
###
# f=1 gives hlogis
hglogis <- function(y,m,s,f) {
	y1 <- (y-m)/s
	ey <- exp(-y1)
	-log(s/f)-y1-(f+1)*log(1+ey)-log(1-(1+ey)^-f)}

### normal intensity
###
hnorm <- function(y,m,s) dnorm(y,m,s,TRUE)-pnorm(y,m,s,FALSE,TRUE)

### Pareto intensity
###
hpareto <- function(y,m,s) (s+1)/(m*s+y)

### Student t intensity
###
hstudent <- function(y,m,s,f){
	pst <- function(u,f){
		t <- 0.5*pbeta(f/(f+u^2),f/2,0.5)
		ifelse(u<0,t,1-t)}
	t <- (f+1)/2
	u <- (y-m)/s
	lgamma(t)-lgamma(f/2)-log(f)/2-(t)*log(1+u^2/f)
		-log(pi)/2-log(1-pst(u,f))}
### Weibull intensity
###
hweibull <- function(y,s,m) log(s)+(s-1)*log(y)-s*log(m)

### generalized Weibull intensity
###
# Mudholkar, Srivastava, & Freimer (1995) Technometrics 37: 436-445
hgweibull <- function(y,s,m,f) {
	y1 <- y/m
	y2 <- y1^s
	y3 <- exp(-y2)
	log(s*f/m)+(s-1)*log(y1)+(f-1)*log(1-y3)-y2-log(1-(1-y3)^f)}

### skew Laplace intensity
###
hskewlaplace <- function(y,m=0,s=1,f=1){
	plp <- function(u)
		ifelse(u>0,1-exp(-f*abs(u))/(1+f^2),f^2*exp(-abs(u)/f)/(1+f^2))
	log(f)+ifelse(y>m,-f*(y-m),(y-m)/f)/s-log((1+f^2)*s)-
		log(1-plp((y-m)/s))}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  You should have received a copy of the GNU General Public License
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#  kalsurv(response=NULL, intensity="exponential", distribution="Pareto",
#	depend="independence", update="Markov", mu=NULL, shape=NULL,
#	renewal=TRUE, density=FALSE, censor=NULL, delta=NULL, ccov=NULL,
#	tvcov=NULL,preg=NULL, ptvc=NULL, pbirth=NULL, pintercept=NULL,
#	pshape=NULL, pinitial=1, pdepend=NULL, pfamily=NULL,
#	envir=parent.frame(), print.level=0, ndigit=10, gradtol=0.00001,
#	steptol=0.00001, iterlim=100, fscale=1, typsiz=abs(p),
#	stepmax=10*sqrt(p%*%p))
#
#  DESCRIPTION
#
#    Function to fit various distributions inserted in a Pareto, gamma, or
# Weibull distribution with serial dependence or gamma frailties using
# Kalman-type update for event histories.

kalsurv <- function(response=NULL, intensity="exponential",
	distribution="Pareto", depend="independence", update="Markov",
	mu=NULL, shape=NULL, renewal=TRUE, density=FALSE, censor=NULL,
	delta=NULL, ccov=NULL, tvcov=NULL, preg=NULL, ptvc=NULL, pbirth=NULL,
	pintercept=NULL, pshape=NULL, pinitial=1, pdepend=NULL,
	pfamily=NULL, envir=parent.frame(), print.level=0, ndigit=10,
	gradtol=0.00001, steptol=0.00001, iterlim=100, fscale=1,
	typsiz=abs(p), stepmax=10*sqrt(p%*%p)){
#
# Pareto (beta) likelihood function for serial dependence
#
ksurvb <- function(p){
	if(rf)b <- mu1(p)
	if(sf)v <- sh1(p[nps1:np])
	z <- .C("ksurvb",
		p=as.double(p),
		y=as.double(resp$response$y),
		x=as.double(resp$ccov$ccov),
		cens=as.integer(resp$response$censor),
		nind=as.integer(nind),
		nobs=as.integer(nobs),
		nbs=as.integer(n),
		nccov=as.integer(nccov),
		model=as.integer(mdl),
		density=as.integer(density),
		pfamily=as.integer(!is.null(pfamily)),
		dep=as.integer(dep),
		birth=as.integer(birth),
		tvc=as.integer(tvc),
		tvcov=as.double(resp$tvcov$tvcov),
		fit=as.integer(0),
		pred=double(n),
		rpred=double(n),
		renewal=as.integer(renewal),
		rf=as.integer(rf),
		bb=as.double(b),
		sf=as.integer(sf),
		vv=as.double(v),
		like=double(1),
		DUP=FALSE,
		PACKAGE="event")
	z$like}
#
# gamma and Weibull likelihood function for serial dependence
#
ksurvg <- function(p){
	if(rf)b <- mu1(p)
	if(sf)v <- sh1(p[nps1:np])
	z <- .C("ksurvg",
		p=as.double(p),
		y=as.double(resp$response$y),
		x=as.double(resp$ccov$ccov),
		cens=as.integer(resp$response$censor),
		nind=as.integer(nind),
		nobs=as.integer(nobs),
		nbs=as.integer(n),
		nccov=as.integer(nccov),
		model=as.integer(mdl),
		distribution=as.integer(dst),
		density=as.integer(density),
		dep=as.integer(dep),
		birth=as.integer(birth),
		tvc=as.integer(tvc),
		tvcov=as.double(resp$tvcov$tvcov),
		fit=as.integer(0),
		pred=double(n),
		rpred=double(n),
		renewal=as.integer(renewal),
		rf=as.integer(rf),
		bb=as.double(b),
		sf=as.integer(sf),
		vv=as.double(v),
		like=double(1),
		DUP=FALSE,
		PACKAGE="event")
	z$like}
#
# Pareto (beta) likelihood function for frailty
#
frailb <- function(p){
	if(rf)b <- mu1(p)
	if(sf)v <- sh1(p[nps1:np])
	z <- .C("frailb",
		p=as.double(p),
		y=as.double(resp$response$y),
		x=as.double(resp$ccov$ccov),
		cens=as.integer(resp$response$censor),
		nind=as.integer(nind),
		nobs=as.integer(nobs),
		nbs=as.integer(n),
		nccov=as.integer(nccov),
		model=as.integer(mdl),
		density=as.integer(density),
		dep=as.integer(dep),
		birth=as.integer(birth),
		tvc=as.integer(tvc),
		tvcov=as.double(resp$tvcov$tvcov),
		fit=as.integer(0),
		pred=double(n),
		rpred=double(n),
		rf=as.integer(rf),
		bb=as.double(b),
		sf=as.integer(sf),
		vv=as.double(v),
		frser=as.integer(frser),
		like=double(1),
		DUP=FALSE,
		PACKAGE="event")
	z$like}
call <- sys.call()
#
# check model
#
tmp <- c("exponential","Weibull","gamma","gen logistic",
	"log normal","log logistic","log Cauchy","log Laplace")
mdl <- match(intensity <- match.arg(intensity,tmp),tmp)
tmp <- c("Pareto","gamma","Weibull")
dst <- match(distribution <- match.arg(distribution,tmp),tmp)
if(distribution=="Pareto"){
	if(depend=="frailty")surv <- frailb
	else surv <- ksurvb}
else {
	if(depend=="frailty")stop("not yet available") # surv <- frailg 
	else surv <- ksurvg}
depend <- match.arg(depend,c("independence","serial","frailty"))
tmp <- c("elapsed Markov","serial","Markov","event","cumulated","count","kalman","time")
dep <- match(update <- match.arg(update,tmp),tmp)
rf <- !is.null(mu)
sf <- !is.null(shape)
#
# check if a data object is being supplied
#
type <- "unknown"
respenv <- exists(deparse(substitute(response)),envir=parent.frame())&&
	inherits(response,"repeated")&&!inherits(envir,"repeated")
if(respenv){
	if(dim(response$response$y)[2]>1)
		stop("kalsurv only handles univariate responses")
	if(!is.null(response$NAs)&&any(response$NAs))
		stop("kalsurv does not handle data with NAs")
	type <- response$response$type}
envname <- if(respenv)deparse(substitute(response))
	else if(!is.null(class(envir)))deparse(substitute(envir))
	else NULL
#
# if envir, remove extra (multivariate) responses
#
if(!respenv&&inherits(envir,"repeated")){
	if(!is.null(envir$NAs)&&any(envir$NAs))
		stop("kalsurv does not handle data with NAs")
	cn <- deparse(substitute(response))
	if(length(grep("\"",cn))>0)cn <- response
	if(length(cn)>1)stop("only one response variable allowed")
	response <- envir
	col <- match(cn,colnames(response$response$y))
	if(is.na(col))stop(paste("response variable",cn,"not found"))
	type <- response$response$type[col]
	if(dim(response$response$y)[2]>1){
		response$response$y <- response$response$y[,col,drop=FALSE]
		if(!is.null(response$response$delta)){
			response$response$delta <- response$response$delta[,col,drop=FALSE]
			if(all(response$response$delta==1)||all(is.na(response$response$delta)))response$response$delta <- NULL}
		if(!is.null(response$response$censor)){
			response$response$censor <- response$response$censor[,col,drop=FALSE]
			if(all(response$response$censor==1)||all(is.na(response$response$censor)))response$response$censor <- NULL}}}
#
# find the covariates
#
if(respenv||inherits(envir,"repeated")){
	if(rf)resp <- rmna(response$response)
	else {
		resp <- response
		# time-constant covariates
		if(is.null(ccov))resp$ccov <- NULL
		else if(inherits(ccov,"formula")){
			if(any(is.na(match(rownames(attr(terms(ccov,data=response),"factors")),colnames(response$ccov$ccov)))))
				stop("ccov covariate(s) not found")
			tmp <- wr(ccov,data=response,expand=FALSE)$design
			resp$ccov$ccov <- tmp[,-1,drop=FALSE]
			rm(tmp)}
		else stop("ccov must be a W&R formula")
		# time-varying covariates
		if(is.null(tvcov))resp$tvcov <- NULL
		else if(inherits(tvcov,"formula")){
			if(any(is.na(match(rownames(attr(terms(tvcov,data=response),"factors")),colnames(response$tvcov$tvcov)))))
				stop("tvcov covariate(s) not found")
			tmp <- wr(tvcov,data=response)$design
			resp$tvcov$tvcov <- tmp[,-1,drop=FALSE]
			rm(tmp)}
		else stop("tvcov must be a W&R formula")}
	nccov <- if(rf||is.null(resp$ccov$ccov)) 0
		 else  dim(resp$ccov$ccov)[2]
	ttvc <- if(rf||is.null(resp$tvcov$tvcov)) 0
		 else  dim(resp$tvcov$tvcov)[2]}
else {
	# make response object
	if(inherits(response,"response")){
		resp <- response
		type <- response$type}
	else {
		if(is.null(censor)){
			if(is.matrix(response)||is.data.frame(response))
				censor <- rep(1,dim(response)[1])
			else if(is.list(response))
				censor <- rep(1,length(response))}
		resp <- restovec(response,censor=censor,delta=delta)}
	# time-constant covariates
	if(is.null(ccov))nccov <- 0
	else {
		if(!inherits(ccov,"tccov")){
			ccname <- deparse(substitute(ccov))
			if((is.matrix(ccov)&&is.null(colnames(ccov)))){
				ccname <- deparse(substitute(ccov))
				if(dim(ccov)[2]>1){
					tmp <- NULL
					for(i in 1:dim(ccov)[2])
						tmp <- c(tmp,paste(ccname,i,sep=""))
					ccname <- tmp}}
			ccov <- tcctomat(ccov,names=ccname)}
		nccov <- if(rf) 0 else dim(ccov$ccov)[2]}
	# time-varying covariates
	if(is.null(tvcov))ttvc <- 0
	else {
		if(!inherits(tvcov,"tvcov")){
			tvcname <- deparse(substitute(tvcov))
			if(is.list(tvcov)&&dim(tvcov[[1]])[2]>1){
				if(is.null(colnames(tvcov[[1]]))){
					tvcname <- deparse(substitute(tvcov))
					tmp <- NULL
					for(i in 1:dim(tvcov[[1]])[2])tmp <- c(tmp,paste(tvcname,i,sep=""))
					tvcname <- tmp}
				else tvcname <- colnames(tvcov[[1]])}
			tvcov <- tvctomat(tvcov, names=tvcname)}
		ttvc <- if(rf) 0 else dim(tvcov$tvcov)[2]}
	resp <- rmna(response=resp, tvcov=tvcov, ccov=ccov)
	if(!is.null(ccov))rm(ccov)
	if(!is.null(tvcov))rm(tvcov)}
if(any(resp$response$y<0))stop("All times must be non-negative")
n <- dim(resp$response$y)[1]
nobs <- nobs(resp)
nind <- length(nobs)
if((inherits(envir,"repeated")&&(length(nobs)!=length(nobs(envir))||
	any(nobs!=nobs(envir))))||(inherits(envir,"tvcov")&&
	(length(nobs)!=length(envir$tvcov$nobs)||any(nobs!=envir$tvcov$nobs))))
	stop("response and envir objects are incompatible")
if(type!="unknown"&&type!="duration"&&type!="continuous")
	stop("duration data required")
#
# if a data object was supplied, modify formulae or functions to read from it
#
mu3 <- sh3 <- NULL
if(respenv||inherits(envir,"repeated")||inherits(envir,"tccov")||inherits(envir,"tvcov")){
	if(is.null(envname))envname <- deparse(substitute(envir))
	if(inherits(mu,"formula")){
		mu3 <- if(respenv)finterp(mu,.envir=response,.name=envname)
			else finterp(mu,.envir=envir,.name=envname)}
	else if(is.function(mu)){
		if(is.null(attr(mu,"model"))){
		        tmp <- parse(text=deparse(mu)[-1])
		        mu <- if(respenv)fnenvir(mu,.envir=response,.name=envname)
		        	else fnenvir(mu,.envir=envir,.name=envname)
		        mu3 <- mu
		        attr(mu3,"model") <- tmp}
		else mu3 <- mu}
	if(inherits(shape,"formula")){
		sh3 <- if(respenv)finterp(shape,.envir=response,.name=envname)
			else finterp(shape,.envir=envir,.name=envname)}
	else if(is.function(shape)){
		if(is.null(attr(shape,"model"))){
		        tmp <- parse(text=deparse(shape)[-1])
		        shape <- if(respenv)fnenvir(shape,.envir=response,.name=envname)
		        	else fnenvir(shape,.envir=envir,.name=envname)
		        sh3 <- shape
		        attr(sh3,"model") <- tmp}
		else sh3 <- shape}}
#
# transform location formula to function and check number of parameters
#
npreg <- length(preg)
mu1 <- sh1 <- v <- b <- NULL
if(inherits(mu,"formula")){
	pr <- if(npreg>0)preg else ptvc
	npr <- length(pr)
	mu2 <- if(respenv)finterp(mu,.envir=response,.name=envname,.expand=is.null(preg))
		else finterp(mu,.envir=envir,.name=envname,.expand=is.null(preg))
	npt1 <- length(attr(mu2,"parameters"))
	if(is.character(attr(mu2,"model"))){
	# W&R formula
		if(length(attr(mu2,"model"))==1){
			mu1 <- function(p) p[1]*rep(1,n)
			attributes(mu1) <- attributes(mu2)
			mu2 <- NULL}}
	else {
	# formula with unknowns
		if(npr!=npt1&&length(ptvc)!=npt1){
			cat("\nParameters are ")
			cat(attr(mu2,"parameters"),"\n")
			stop(paste("preg or ptvc should have",npt1,"estimates"))}
		if(is.list(pr)){
			if(!is.null(names(pr))){
				o <- match(attr(mu2,"parameters"),names(pr))
				pr <- unlist(pr)[o]
				if(sum(!is.na(o))!=length(pr))
					stop("invalid estimates for mu - probably wrong names")}
			else pr <- unlist(pr)
			if(npreg>0)preg <- pr else ptvc <- pr}}
	if(!is.null(mu2)){
		if(inherits(envir,"tccov")){
			cv <- covind(response)
			mu1 <- function(p) mu2(p)[cv]
			attributes(mu1) <- attributes(mu2)}
		else {
			mu1 <- mu2
			rm(mu2)}}}
else if(is.function(mu))mu1 <- mu
#
# give appropriate attributes to mu1 for printing
#
if(!is.null(mu1)&&is.null(attr(mu1,"parameters"))){
	attributes(mu1) <- if(is.function(mu)){
		if(!inherits(mu,"formulafn")){
			if(respenv)attributes(fnenvir(mu,.envir=response))
			else attributes(fnenvir(mu,.envir=envir))}
		else attributes(mu)}
		else {
			if(respenv)attributes(fnenvir(mu1,.envir=response))
			else attributes(fnenvir(mu1,.envir=envir))}}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(mu1))length(attr(mu1,"parameters"))
	else if(is.null(mu1))NULL
	else npt1
if(!is.null(nlp)){
	if(is.null(ptvc)&&nlp!=npreg)
		stop(paste("preg should have",nlp,"initial estimates"))
	else if(!is.null(ptvc)&&length(ptvc)!=nlp)
		stop(paste("ptvc should have",nlp,"initial estimates"))}
#
# transform shape formula to function and check number of parameters
#
nps <- length(pshape)
if(inherits(shape,"formula")){
	sh2 <- if(respenv)finterp(shape,.envir=response,.name=envname)
		else finterp(shape,.envir=envir,.name=envname)
	npt2 <- length(attr(sh2,"parameters"))
	if(is.character(attr(sh2,"model"))){
	# W&R formula
		if(length(attr(sh2,"model"))==1){
			sh1 <- function(p) p[npl1]*rep(1,n)
			attributes(sh1) <- attributes(sh2)
			sh2 <- NULL}}
	else {
	# formula with unknowns
		if(nps!=npt2){
			cat("\nParameters are ")
			cat(attr(sh2,"parameters"),"\n")
			stop(paste("pshape should have",npt2,"estimates"))}
		if(is.list(pshape)){
			if(!is.null(names(pshape))){
				o <- match(attr(sh2,"parameters"),names(pshape))
				pshape <- unlist(pshape)[o]
				if(sum(!is.na(o))!=length(pshape))
					stop("invalid estimates for shape - probably wrong names")}
			else pshape <- unlist(pshape)}}
	if(!is.null(sh2)){
		if(inherits(envir,"tccov")){
			cv <- covind(response)
			sh1 <- function(p) sh2(p)[cv]
			attributes(sh1) <- attributes(sh2)}
		else {
			sh1 <- sh2
			rm(sh2)}}}
else if(is.function(shape))sh1 <- shape
#
# give appropriate attributes to sh1 for printing
#
if(!is.null(sh1)&&is.null(attr(sh1,"parameters")))
	attributes(sh1) <- if(is.function(shape)){
		if(!inherits(shape,"formulafn")){
			if(respenv)attributes(fnenvir(shape,.envir=response))
			else attributes(fnenvir(shape,.envir=envir))}
		else attributes(shape)}
		else {
			if(respenv)attributes(fnenvir(sh1,.envir=response))
			else attributes(fnenvir(sh1,.envir=envir))}
#
# check that correct number of estimates was supplied
#
nlp <- if(is.function(shape))length(attr(sh1,"parameters"))
	else if(is.null(shape))NULL
	else npt2
if(!is.null(nlp)&&nlp!=nps)
	stop(paste("pshape should have",nlp,"initial estimates"))
#
# check that appropriate functions supplied
#
if(rf&&!is.function(mu1))stop("mu must be a formula or function")
if(sf&&!is.function(sh1))stop("shape must be a formula or function")
birth <- !is.null(pbirth)
tvc <- length(ptvc)
if(!rf&&(ttvc>0&&tvc!=ttvc||ttvc==0&&tvc>0))
	stop(paste(ttvc,"initial estimates of coefficients for time-varying covariates must be supplied"))
if(rf&&birth)stop("Birth models cannot be fitted with a mean function")
if(intensity=="exponential"){
	sf <- FALSE
	pshape <- NULL}
else {
	if(is.null(pshape))
		stop("Initial value of the shape parameter must be supplied")
	if(!sf){
		if(pshape<=0)stop("Shape must be positive")
		else pshape <- log(pshape)}}
if(intensity=="gen logistic"){
	if(is.null(pintercept))
		stop("Initial value of the intercept parameter must be supplied")}
else pintercept <- NULL
np <- 1+npreg+(depend=="serial")+birth+nps+(!is.null(pintercept))+
	(!is.null(pfamily))
if(pinitial<=0)stop("Estimate of initial parameter must be > 0")
else pinitial <- log(pinitial)
#
# check what dependence is required
#
frser <- FALSE
if(depend=="independence"){
	pdepend <- NULL
	dep <- 0}
else if(depend=="serial"){
	if(update=="time")stop("time update can only be used with frailty")
	if(is.null(pdepend))
		stop("An estimate of the dependence parameter must be supplied")
	else if(pdepend<=0|pdepend>=1)
		stop("Dependence parameter must be between 0 and 1")
	else pdepend <- log(pdepend/(1-pdepend))}
else if(depend=="frailty"){
	frser <- !is.null(pdepend)
	if(frser){
		if(pdepend<=0)stop("Dependence parameter must be positive")
		np <- np+1
		pdepend <- log(pdepend)}
	if(update=="time")dep <- 1
	else {
		dep <- 0
		update <- "no"}}
if(is.null(resp$response$censor))
	resp$response$censor <- rep(1,n)
#
# check number of parameter estimates
#
if(!is.null(pfamily)){
	if(depend=="frailty")
		stop("pfamily cannot be used with frailty dependence")
	if(distribution!="Pareto")
		stop("pfamily can only be used with Pareto distribution")}
if(rf&&npreg>0)nccov <- npreg-1
if(!rf&&nccov+1!=npreg)
	stop(paste(nccov+1,"regression estimates must be supplied"))
if(ttvc>0)np <- np+tvc
#
# set up for location function and make sure it gives correct number
# of values
#
if(rf){
	if(tvc>0&&nccov>0)
		stop("With a mean function, initial estimates must be supplied either in preg or in ptvc")
	if(tvc>0){
		if(length(mu1(ptvc))!=n)
			stop("The mu function must provide an estimate for each observation")
		tvc <- tvc-1
		np <- np+tvc+1}
	else {
		lp <- length(mu1(preg))
		if(lp==1){
			if(nccov==0)mu1 <- function(p) rep(p[1],nind)
			else stop("Number of estimates does not correspond to mu function")}
		else if(lp!=nind)
			stop("The mu function must provide an estimate for each individual")}}
if(sf&&length(sh1(pshape))!=n)
	stop("The shape function must provide an estimate for each observation")
#
# check that the likelihood function gives appropriate values and call nlm
#
nps1 <- np-nps+1
p <- c(preg,pbirth,ptvc,pinitial,pdepend,pfamily,pintercept,pshape)
if(fscale==1)fscale <- surv(p)
if(is.na(surv(p)))
	stop("Likelihood returns NAs: probably invalid initial values")
z0 <- nlm(surv, p=p, hessian=TRUE, print.level=print.level,
	typsiz=typsiz, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
p <- z0$estimate
like <- z0$minimum
if(!is.null(resp$response$delta))like <- like-
	if(length(resp$response$delta)==1)n*log(resp$response$delta)
	else sum(log(resp$response$delta))
#
# calculate se's
#
a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
	else qr(z0$hessian)$rank
if(a==np)cov <- solve(z0$hessian)
else cov <- matrix(NA,ncol=np,nrow=np)
se <- sqrt(diag(cov))
corr <- cov/(se%o%se)
dimnames(corr) <- list(1:np,1:np)
#
# calculate recursive fitted values
#
if(mdl==4)z <- list()
else {
	z <- if(depend=="frailty"){
		if(rf)b <- mu1(p)
		if(sf)v <- sh1(p[nps1:np])
		z <- .C("frailb",
			p=as.double(p),
			y=as.double(resp$response$y),
			x=as.double(resp$ccov$ccov),
			cens=as.integer(resp$response$censor),
			nind=as.integer(nind),
			nobs=as.integer(nobs),
			nbs=as.integer(n),
			nccov=as.integer(nccov),
			model=as.integer(mdl),
			density=as.integer(density),
			dep=as.integer(dep),
			birth=as.integer(birth),
			tvc=as.integer(tvc),
			tvcov=as.double(resp$tvcov$tvcov),
			fit=as.integer(1),
			pred=double(n),
			rpred=double(n),
			rf=as.integer(rf),
			bb=as.double(b),
			sf=as.integer(sf),
			vv=as.double(v),
			frser=as.integer(frser),
			like=double(1),
			DUP=FALSE,
			PACKAGE="event")}
	else if(distribution=="Pareto"){
		if(rf)b <- mu1(p)
		if(sf)v <- sh1(p[nps1:np])
		z <- .C("ksurvb",
			p=as.double(p),
			y=as.double(resp$response$y),
			x=as.double(resp$ccov$ccov),
			cens=as.integer(resp$response$censor),
			nind=as.integer(nind),
			nobs=as.integer(nobs),
			nbs=as.integer(n),
			nccov=as.integer(nccov),
			model=as.integer(mdl),
			density=as.integer(density),
			pfamily=as.integer(!is.null(pfamily)),
			dep=as.integer(dep),
			birth=as.integer(birth),
			tvc=as.integer(tvc),
			tvcov=as.double(resp$tvcov$tvcov),
			fit=as.integer(1),
			pred=double(n),
			rpred=double(n),
			renewal=as.integer(renewal),
			rf=as.integer(rf),
			bb=as.double(b),
			sf=as.integer(sf),
			vv=as.double(v),
			like=double(1),
			DUP=FALSE,
			PACKAGE="event")}
	else {
		if(rf)b <- mu1(p)
		if(sf)v <- sh1(p[nps1:np])
		z <- .C("ksurvg",
			p=as.double(p),
			y=as.double(resp$response$y),
			x=as.double(resp$ccov$ccov),
			cens=as.integer(resp$response$censor),
			nind=as.integer(nind),
			nobs=as.integer(nobs),
			nbs=as.integer(n),
			nccov=as.integer(nccov),
			model=as.integer(mdl),
			distribution=as.integer(dst),
			density=as.integer(density),
			dep=as.integer(dep),
			birth=as.integer(birth),
			tvc=as.integer(tvc),
			tvcov=as.double(resp$tvcov$tvcov),
			fit=as.integer(1),
			pred=double(n),
			rpred=double(n),
			renewal=as.integer(renewal),
			rf=as.integer(rf),
			bb=as.double(b),
			sf=as.integer(sf),
			vv=as.double(v),
			like=double(1),
			DUP=FALSE,
			PACKAGE="event")}
	for(i in 1:n)if(resp$response$y[i]==0){
		z$pred[i] <- z$pred[i-1]
		z$rpred[i] <- z$rpred[i-1]}}
#
# return appropriate attributes on functions
#
if(!is.null(mu3))mu1 <- mu3
if(!is.null(sh3))sh1 <- sh3
z <- list(
	call=call,
	intensity=intensity,
	distribution=distribution,
	pfamily=!is.null(pfamily),
	frser=frser,
	mu=mu1,
	npr=1+nccov+tvc+birth,
	shape=sh1,
	nps=nps,
	density=density,
	depend=depend,
	update=update,
	birth=birth,
	renewal=renewal,
	response=resp$response,
	pred=z$pred,
	rpred=z$rpred,
	ccov=resp$ccov,
	tvcov=resp$tvcov,
	maxlike=like,
	aic=like+np,
	df=n-np,
	npt=np,
	npv=npreg,
	coefficients=p,
	se=se,
	cov=cov,
	corr=corr,
	grad=z0$gradient,
	iterations=z0$iterations,
	code=z0$code)
class(z) <- if(mdl==4)"kalsurv" else c("kalsurv","recursive")
return(z)}

### standard methods
###

deviance.kalsurv <- function(z) 2*z$maxlike

fitted.kalsurv <- function(z, recursive=TRUE)
if(recursive) z$rpred else z$pred

residuals.kalsurv <- function(z, recursive=TRUE)
if(recursive) z$response$y-z$rpred else z$response$y-z$pred

### print method
###
print.kalsurv <- function(z,digits=max(3,.Options$digits-3),correlation=TRUE){
tvc <- !is.null(z$tvcov)
expm <- z$intensity!="exponential"&&!is.function(z$shape)
glm <- z$intensity=="gen logistic"
npt <- if(is.function(z$shape)) z$npt-z$nps else z$npt
deppar <- (z$depend!="independence")&&(z$depend!="frailty")
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
if(z$code>2)cat("Warning: no convergence - error",z$code,"\n\n")
cat("Number of subjects    ",length(nobs(z)),"\n")
cat("Number of observations",length(z$response$y),"\n")
cat(z$distribution,"distribution ")
if(z$renewal){
	if(!z$birth)cat("with renewal process\n")
	else cat("with birth process\n")}
else {
	cat("with zero origin\n")
	if(z$birth)cat("and birth process\n")}
if(z$density)cat(z$intensity," density",sep="")
else cat(z$intensity," intensity",sep="")
if(z$depend=="independence")cat(" with independence\n")
else if(z$depend=="frailty")
	cat(" with",z$depend,"dependence",if(!z$frser)"and",
		z$update,"weight",if(z$frser)"and AR","\n")
else cat(" with ",z$update," update\n",sep="")
cat("\n-Log likelihood   ",z$maxlike,"\n")
cat("Degrees of freedom",z$df,"\n")
cat("AIC               ",z$aic,"\n")
cat("Iterations        ",z$iterations,"\n\n")
cat("Location parameters\n")
if(!is.null(attr(z$mu,"formula")))
	cat(deparse(attr(z$mu,"formula")),sep="\n")
else if(!is.null(attr(z$mu,"model"))){
	t <- deparse(attr(z$mu,"model"))
	t[1] <- sub("expression\\(","",t[1])
	t[length(t)] <- sub("\\)$","",t[length(t)])
	cat(t,sep="\n")}
coef.table <- cbind(z$coef[1:z$npr],z$se[1:z$npr])
if(inherits(z$mu,"formulafn"))
	cname <- if(is.character(attr(z$mu,"model")))attr(z$mu,"model")
		else attr(z$mu,"parameters")
else {
	cname <- "(Intercept)"
	if(z$npv>0)cname <- c(cname,colnames(z$ccov$ccov))
	if(z$birth)cname <- c(cname,"birth")
	if(tvc)cname <- c(cname,colnames(z$tvcov$tvcov))}
dimnames(coef.table) <- list(cname, c("estimate","se"))
print.default(coef.table, digits=digits, print.gap=2)
if(is.function(z$shape))cat("\nDependence parameters\n")
else cat("\nNonlinear parameters\n")
coef <- exp(z$coef[(npt-deppar-z$frser-z$pfamily-glm-expm):npt])
cname <- "initial"
if(deppar){
	cname <- c(cname,"depend")
	coef[2] <- coef[2]/(1+coef[2])}
if(z$frser)cname <- c(cname,"AR")
if(z$pfamily){
	cname <- c(cname,"family")
	coef[length(coef)-glm-expm] <- z$coef[npt-glm-expm]}
if(glm){
	cname <- c(cname,"intercept")
	coef[length(coef)-expm] <- z$coef[npt-expm]
	if(expm){
		cname <- c(cname,"asymptote")
		coef[length(coef)] <- 1/coef[length(coef)]}}
else if(expm)cname <- c(cname,"shape")
coef.table <- cbind(z$coef[(npt-deppar-z$frser-z$pfamily-glm-expm):npt],
	z$se[(npt-deppar-z$frser-z$pfamily-glm-expm):npt],coef)
dimnames(coef.table) <- list(cname, c("estimate","se","parameter"))
print.default(coef.table, digits=digits, print.gap=2)
if(z$depend=="frailty"){
	tmp <- trigamma(exp(-z$coef[npt-deppar-expm-z$frser]))
	cat("Correlation =",tmp/(tmp+trigamma(1)),"\n")}
if(is.function(z$sh1)||is.function(z$shape)){
	cat("\nShape parameters\n")
	if(!is.null(attr(z$sh1,"formula")))
		cat(deparse(attr(z$sh1,"formula")),sep="\n")
	else if(!is.null(attr(z$sh1,"model"))){
		t <- deparse(attr(z$sh1,"model"))
		t[1] <- sub("expression\\(","",t[1])
		t[length(t)] <- sub("\\)$","",t[length(t)])
		cat(t,sep="\n")}
	cname <- if(is.character(attr(z$sh1,"model")))
			attr(z$sh1,"model")
		else attr(z$sh1,"parameters")
	coef.table <- cbind(z$coef[(z$npt-z$nps+1):z$npt],
		z$se[(z$npt-z$nps+1):z$npt])
	dimnames(coef.table) <- list(cname, c("estimate","se"))
	print.default(coef.table, digits=digits, print.gap=2)}
if(correlation){
	cat("\nCorrelation matrix\n")
	print.default(z$corr, digits=digits)}}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#  km(times, censor=1, group=1, freq=1, cdf=FALSE)
#  print.km(z)
#  plot.km(z, surv, times=NULL, group=1, cdf=FALSE, plot=TRUE, add=FALSE,
#	xlim, ylim=c(0,1), main=NULL, xlab="Time",
#	ylab=NULL, lty=NULL, ...)
#  plot.surv(times, survivor, group=1, cdf=FALSE, add=FALSE, xlim=NULL,
#	ylim=c(0,1), main=NULL, xlab="Time", ylab=NULL, lty=NULL, ...)
#  plot.intensity(z, ...)
#  plot.intensity.default(times, censor=1, group=1, cl=1, mix=1, ylim=c(0,1),
#	ylab="p", xlab="Time", main="Empirical Hazard Function(s)")
#  plot.intensity.km(z, add=FALSE, xlab="Time", ylab="Hazard", type="l",
#	lty=NULL, ...)
#  plot.dist.km(z)

#
#  DESCRIPTION
#
#    Functions to compute and plot Kaplan-Meier estimates

### function to compute Kaplan-Meier estimates
###
km <- function(times, censor=1, group=1, freq=1, cdf=FALSE){
#
# list supplied, make vectors
#
cens <- gp <- tt <- NULL
j <- 0
if(is.list(times)){
	for(i in times){
		j <- j+1
		tt <- c(tt,i)
		cens <- c(cens,rep(1,length(i)))
		if(!missing(censor))cens[length(cens)] <- censor[j]
		if(!missing(group))gp <- c(gp,rep(group[j],length(i)))}
	times <- tt
	censor <- cens
	if(!missing(group))group <- gp
	rm(tt,cens,gp)}
#
# check values supplied are appropriate
#
group <- as.numeric(group)
if(!is.vector(times,mode="numeric"))stop("times must be a vector")
else if(any(times<0))stop("negative times")
if(length(group)==1)group <- rep(1,length(times))
if(length(censor)==1)censor <- rep(1,length(times))
if(length(group)!=length(times)||length(times)!=length(censor))
	stop("All vectors must be the same length")
if(any(censor!=0&censor!=1))
	stop("Censor vector must be zeros and ones")
#
# replicate data if grouped
#
if(length(freq)>1&any(freq>1)){
	times <- rep(times,freq)
	censor <- rep(censor,freq)
	group <- rep(group,freq)}
#
# calculate unique times of events
#
ss <- cens <- gg <- ff <- cg <- NULL
for(i in unique(group)){
	n <- sum(group==i)
	o <- order(times[group==i])
	cens2 <- censor[group==i][o]
	times2 <- times[group==i][o]
	f2 <- as.vector(table(list(1-cens2,times2)))
	if(!any(cens2==0)||!any(cens2==1)){
		ff2 <- rep(0,2*length(f2))
		ff2[seq(1,2*length(f2)-1,by=2)] <- f2
		f2 <- ff2}
	l2 <- length(f2)
	g2 <- rep(i,l2)
	c2 <- rep(c(1,0),l2/2)
	s2 <- vector(mode="numeric",length=l2)
	s2[seq(1,l2-1,by=2)] <- unique(times2)
	s2[seq(2,l2,by=2)] <- unique(times2)
	ss <- c(ss,s2)
	cens <- c(cens,c2)
	gg <- c(gg,g2)
	ff <- c(ff,f2)
	cg <- c(cg,f2[l2]==0)}
#
# calculate survivor curve
#
j <- (ff!=0)|(cens==1)
k <- seq(1:length(ff))
n <- sum(j)
t <- cc <- f <- g <- tt <- vector(mode="numeric",length=n)
t <- ss[k*j]
cc <- cens[k*j]
f <- ff[k*j]
g <- gg[k*j]
v <- r <- s <- rep(0,n)
for(i in unique(group)){
	tt[n:1] <- cumsum(f[n:1]*(g[n:1]==i))
	r <- r+tt*(g==i)
	tmp <- cumsum(log(ifelse(tt==f|tt==0,1,(tt-f)/tt))*(g==i)*cc)
	s <- s+exp(ifelse(is.na(tmp),0,tmp))*(g==i)*(tt!=f)
	tmp <- cumsum(f/tt/(tt-f)*(g==i)*cc)
	v <- v+ifelse(is.na(tmp),0,tmp)*(g==i)}
v <- s*s*v
#
# set up matrix to return
#
m <- NULL
for(i in unique(g))m <- c(m,sum(g==i))
j <- seq(1,length(g))
j[cumsum(m)] <- cg*j[cumsum(m)]
z <- cbind(t[j],g[j],r[j],s[j],v[j])
dp <- rep(TRUE,dim(z)[1])
for(i in 2:dim(z)[1])if(all(z[i,]==z[i-1,],na.rm=TRUE))dp[i] <- FALSE
z <- z[dp,]
colnames(z) <- c("Time","Group","At risk","S(t)","Var(S)")
rownames(z) <- paste(rep("",dim(z)[1]))
class(z) <- "km"
attr(z,"cdf") <- cdf
z}

### print method
###
print.km <- function(z){
attr(z,"class") <- attr(z,"cdf") <- NULL
print.default(z)}

### function to plot survivor curves output from km
###
plot.km <- function(z, add=FALSE, xlim=NULL, ylim=c(0,1), main=NULL, xlab="Time",
	ylab=NULL, lty=NULL, ...){
surv <- z[,4]
times <- z[,1]
group <- z[,2]
cdf <- attr(z,"cdf")
kms <- ttt <- NULL
#
# plot one curve per group
#
k <- ln <- lt <- 0
for(i in unique(group)){
	if(is.null(lty))lt <- lt%%4+1
	else lt <- lty[k <- k+1]
	if(length(group)>1&&length(unique(group))>1){
		j <- (cumsum(group==i)+ln)*(group==i)
		s <- surv[j]
		if(!missing(times))t <- times[j]
		ln <- ln+sum(group==i)}
	else {
		s <- surv
		t <- times}
	n <- 2*length(s)-1
	# prepare survivor function
	km <- rep(0,n)
	km[seq(1,n,by=2)] <- s
	km[seq(2,n-1,by=2)] <- s[1:(length(s)-1)]
	km <- c(1,1,km)
	if(cdf){
		km <- 1-km
		if(is.null(ylab))ylab <- "Failure probability"
		if(missing(main))main <- "Kaplan-Meier cumulative probability curve"}
	else {
		if(is.null(ylab))ylab <- "Survival probability"
		if(missing(main))main <- "Kaplan-Meier survivor curve"}
	# prepare times
	tt <- rep(0,n)
	tt[seq(1,n,by=2)] <- t
	tt[seq(2,n-1,by=2)] <- t[1:(length(s)-1)]
	tt <- c(0,tt,tt[length(tt)])
	# plot
	if(add) lines(tt, km, lty=lt,...)
	else {
		if(is.null(xlim))
			xlim <- c(min(times)-1,max(times+1))
		plot(tt, km, type="l", xlim=xlim, ylim=ylim,
			xlab=xlab, ylab=ylab, main=main,
			lty=lt, ...)}
	add <- TRUE
	ttt <- c(ttt,tt)
	kms <- c(kms,km)}
invisible(cbind(ttt,kms))}

### function to plot survivor curves
###
plot.surv <- function(times=NULL, survivor, group=1, cdf=FALSE, add=FALSE,
	xlim=NULL, ylim=c(0,1), main=NULL, xlab="Time", ylab=NULL, lty=NULL,
	...){
if(missing(times)||missing(survivor))
	stop("times and survivor curve must be supplied")
#
# plot one curve per group
#
kms <- ttt <- NULL
k <- ln <- lt <- 0
for(i in unique(group)){
	if(is.null(lty))lt <- lt%%4+1
	else lt <- lty[k <- k+1]
	if(length(group)>1&&length(unique(group))>1){
		j <- (cumsum(group==i)+ln)*(group==i)
		s <- survivor[j]
		if(!missing(times))t <- times[j]
		ln <- ln+sum(group==i)}
	else {
		s <- survivor
		t <- times}
	n <- 2*length(s)-1
	# prepare survivor function
	km <- rep(0,n)
	km[seq(1,n,by=2)] <- s
	km[seq(2,n-1,by=2)] <- s[1:(length(s)-1)]
	km <- c(1,1,km)
	if(cdf){
		km <- 1-km
		if(is.null(ylab))ylab <- "Failure probability"
		if(missing(main))main <- "Kaplan-Meier cumulative probability curve"}
	else {
		if(is.null(ylab))ylab <- "Survival probability"
		if(missing(main))main <- "Kaplan-Meier survivor curve"}
	# prepare times
	if(is.null(times)){
		tt <- 1 + (ceiling(1:(n+1)/2) - 1)%%length(s)
		tt <- c(0,tt)}
	else {
		tt <- rep(0,n)
		tt[seq(1,n,by=2)] <- t
		tt[seq(2,n-1,by=2)] <- t[1:(length(s)-1)]
		tt <- c(0,tt,tt[length(tt)])}
	# plot
	if(add)lines(tt, km, lty=lt,...)
	else {
		if(is.null(xlim))
			xlim <- c(min(times)-1,max(times+1))
		plot(tt, km, type="l", xlim=xlim, ylim=ylim,
			xlab=xlab, ylab=ylab, main=main,
			lty=lt, ...)}
	add <- TRUE
	ttt <- c(ttt,tt)
	kms <- c(kms,km)}
invisible(cbind(ttt,kms))}

### functions to plot intensity function
###
plot.intensity <- function(z, ...) UseMethod("plot.intensity")

plot.intensity.default <- function(times, censor=1, group=1, colour=TRUE,
	mix=1, ylim=c(0,1), ylab="p", xlab="Time",
	main="Empirical Hazard Function(s)"){
censor2 <- censor
if(length(group)==1)group <- rep(1,length(times))
if(length(censor)==1)censor <- rep(1,length(times))
group2 <- as.numeric(group)
group <- as.factor(group)
#
# set up times
#
index <- order(times)
tim.gr <- cbind(times[index],censor[index],group2[index])
listim <- vector(mode="list",length=nlevels(group))
for(i in 1:nlevels(group))
	listim[[i]] <- tim.gr[(tim.gr[,3]==i),1:2]
tc <- vector(mode="list",length=nlevels(group))
for(i in 1:nlevels(group))
	tc[[i]] <- listim[[i]][(listim[[i]][,2]==0),1]
tnc <- vector(mode="list",length=nlevels(group))
for(i in 1:nlevels(group))
	tnc[[i]] <- listim[[i]][(listim[[i]][,2]==1),1]
breaks <- seq(0,floor(max(times+1)),1)
tccat <- vector(mode="list",length=nlevels(group))
tnccat <- vector(mode="list",length=nlevels(group))
for(i in 1:nlevels(group)){
	tccat[[i]] <- cut(as.numeric(tc[[i]]),breaks,right=FALSE)
	tnccat[[i]] <- cut(as.numeric(tnc[[i]]),breaks,right=FALSE)}
#
# set up frequencies
#
tncfreq <- matrix(ncol=nlevels(group),nrow=(floor(max(times))+1))
tfreq <- matrix(ncol=nlevels(group),nrow=(floor(max(times))+1))
for(i in 1:nlevels(group))tncfreq[,i] <- table(tnccat[[i]])
if(length(censor2)!=1){
	tcfreq <- matrix(ncol=nlevels(group),nrow=(floor(max(times))+1))
	for(i in 1:nlevels(group)){
		tcfreq[,i] <- table(tccat[[i]])
		tfreq[,i] <- tncfreq[,i]+tcfreq[,i]}}
if(length(censor2)==1)tfreq <- tncfreq
cumfreq <- matrix(ncol=nlevels(group),nrow=(floor(max(times))+1))
cumfreq2 <- matrix(ncol=nlevels(group),nrow=(floor(max(times))+1))
risk <- matrix(ncol=nlevels(group),nrow=(floor(max(times))+1))
for(i in 1:nlevels(group)){
	cumfreq[,i] <- cumsum(tfreq[,i])
	cumfreq2[,i] <- c(0,cumfreq[1:floor(max(times)),i])
	risk[,i] <- length(listim[[i]][,1])-cumfreq2[,i]}
#
# calcuate intensities
#
risk <- risk*mix
haz <- tncfreq/risk
time <- 0:floor(max(times))
haz2 <- vector(mode="list",length=nlevels(group))
hm <- vector(mode="list",length=nlevels(group))
for(i in 1:nlevels(group)){
	haz2[[i]] <- cbind(haz[,i],time)
	hm[[i]] <- haz2[[i]][(haz2[[i]][,1]!=0),]}
#
# plot curves
#
if(colour){
	plot(hm[[1]][,2],hm[[1]][,1],ylim=ylim,col=1,xlab=xlab,
		ylab=ylab,main=main,type="l")
	if(nlevels(group)!=1){
		for(i in 2:nlevels(group)){
			lines(hm[[i]][,2],hm[[i]][,1],col=i)}}}
else {
	plot(hm[[1]][,2],hm[[1]][,1],ylim=ylim,col=gray(0),
		xlab=xlab,ylab=ylab,main=main,type="l")
	if(nlevels(group)!=1){
		for(i in 2:nlevels(group)){
			lines(hm[[i]][,2],hm[[i]][,1],
				col=gray(i/(2*nlevels(group))))}}}}

### function to plot intensity function from km
###
plot.intensity.km <- function(z, add=FALSE, xlab="Time", ylab="Hazard",
	type="l", lty=NULL, ...){
hazt <- NULL
group <- unique(z[,2])
if(!is.null(lty)&&length(lty)!=length(group))
	stop("lty must have one value per group")
#
# plot one curve per group
#
k <- lt <- 0
for(i in 1:length(group)){
	if(is.null(lty))lt <- lt%%4+1
	else lt <- lty[k <- k+1]
	# calculate intensities and times
	z1t <- z[z[,2]==group[i],1]
	z4t <- z[z[,2]==group[i],4]
	nt <- length(z1t)
	haz <- 2*(z4t[1:(nt-1)]-z4t[2:nt])*(z1t[2:nt]-z1t[1:(nt-1)])/
		(z4t[1:(nt-1)]+z4t[2:nt])
	# plot
	if(add)	lines(z1t[1:(nt-1)],haz,lty=lt)
	else plot(z1t[1:(nt-1)],haz,type=type,lty=lt,xlab=xlab,
		ylab=ylab,...)
	add <- TRUE
	hazt <- c(hazt,haz)}
invisible(cbind(z[1:(dim(z)[1]-length(group)),1],hazt))}

### function to plot transforms of survivor function to determine what
### distribution it might have
###
plot.dist <- function(z, ...) UseMethod("plot.dist")

plot.dist.km <- function(z){
#
# set up parameters
#
oldpar <- par(mfrow=c(3,3),mar=c(5,4,4,2),font.main=1)
group <- unique(z[,2])
mn <- min(z[,4],na.rm=TRUE)
if(mn<=0)mn <- 0.01
mx <- max(z[,4],na.rm=TRUE)
if(mx>=1)mx <- 0.999
#
# exponential distribution
#
plot(z[z[,2]==group[1],1],log(z[z[,2]==group[1],4]),
	main="log[S(t)] vs t",ylab="",type="l",ylim=c(log(mn),log(mx)),
	xlab="Linear through origin if Exponential Distribution")
if(length(group)>1)for(i in 2:length(group))
	lines(z[z[,2]==group[i],1],log(z[z[,2]==group[i],4]),
	lty=(i-1)%%4+1)
#
# Pareto distribution
#
plot(log(z[z[,2]==group[1],1]),log(z[z[,2]==group[1],4]),
	main="log S(t) vs log(t)",type="l",ylim=c(log(mn),log(mx)),
	xlab="Linear if Pareto Distribution",ylab="")
if(length(group)>1)for(i in 2:length(group))
	lines(log(z[z[,2]==group[i],1]),log(z[z[,2]==group[i],4]),
	lty=(i-1)%%4+1)
#
# extreme value distribution
#
plot(z[z[,2]==group[1],1],log(-log(z[z[,2]==group[1],4])),
	main="log{-Log[S(t)]} vs t",ylab="",type="l",
	ylim=c(log(-log(mx)),log(-log(mn))),
	xlab="Linear if Extreme Value Distribution")
if(length(group)>1)for(i in 2:length(group))
	lines(z[z[,2]==group[i],1],log(-log(z[z[,2]==group[i],4])),
	lty=(i-1)%%4+1)
#
# Weibull distribution
#
plot(log(z[z[,2]==group[1],1]),log(-log(z[z[,2]==group[1],4])),
	main="log{-Log[S(t)]} vs log(t)",type="l",
	ylim=c(log(-log(mx)),log(-log(mn))),
	ylab="",xlab="Linear if Weibull Distribution")
if(length(group)>1)for(i in 2:length(group))
	lines(log(z[z[,2]==group[i],1]),log(-log(z[z[,2]==group[i],4]))
	,lty=(i-1)%%4+1)
#
# normal distribution
#
plot(z[z[,2]==group[1],1],qnorm(1-z[z[,2]==group[1],4]),
	main="qnorm[1-S(t)] vs t",ylab="",type="l",
	ylim=c(qnorm(1-mx),qnorm(1-mn)),
	xlab="Linear if Normal Distribution")
if(length(group)>1)for(i in 2:length(group))
	lines(z[z[,2]==group[i],1],qnorm(1-z[z[,2]==group[i],4]),
	lty=(i-1)%%4+1)
#
# log normal distribution
#
plot(log(z[z[,2]==group[1],1]),qnorm(1-z[z[,2]==group[1],4]),
	main="qnorm[1-S(t)] vs log(t)",ylab="",type="l",
	ylim=c(qnorm(1-mx),qnorm(1-mn)),
	xlab="Linear if Log Normal Distribution")
if(length(group)>1)for(i in 2:length(group))
	lines(log(z[z[,2]==group[i],1]),qnorm(1-z[z[,2]==group[i],4]),
	lty=(i-1)%%4+1)
#
# gamma distribution
#
plot(sqrt(z[z[,2]==group[1],1]),qnorm(1-z[z[,2]==group[1],4]),
	main="qnorm[1-S(t)] vs sqrt(t)",type="l",
	ylim=c(qnorm(1-mx),qnorm(1-mn)),
	ylab="",xlab="Linear if Gamma Distribution")
if(length(group)>1)for(i in 2:length(group))
	lines(sqrt(z[z[,2]==group[i],1]),qnorm(1-z[z[,2]==group[i],4]),
	lty=(i-1)%%4+1)
#
# logistic distribution
#
plot(z[z[,2]==group[1],1],log((1-z[z[,2]==group[1],4])/
	z[z[,2]==group[1],4]),main="log{[1-S(t)]/S(t)} vs t",
	ylim=c(log((1-mx)/mx),log((1-mn)/mn)),
	ylab="",xlab="Linear if Logistic Distribution",type="l")
if(length(group)>1)for(i in 2:length(group))
	lines(z[z[,2]==group[i],1],log((1-z[z[,2]==group[i],4])/
	z[z[,2]==group[i],4]),
	lty=(i-1)%%4+1)
#
# log logistic distribution
#
plot(log(z[z[,2]==group[1],1]),log((1-z[z[,2]==group[1],4])/
	z[z[,2]==group[1],4]),ylab="",
	ylim=c(log((1-mx)/mx),log((1-mn)/mn)),
	main="log{[1-S(t)]/S(t)} vs log(t)",
	xlab="Linear if Log Logistic Distribution",type="l")
if(length(group)>1)for(i in 2:length(group))
	lines(log(z[z[,2]==group[i],1]),log((1-z[z[,2]==group[i],4])/
	z[z[,2]==group[i],4]),lty=(i-1)%%4+1)
par(oldpar)}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#  pbirth(frequencies, p, intensity="negative binomial",
#	type="spectral decomposition", print.level=0, ndigit=10,
#	gradtol=0.00001, steptol=0.00001, fscale=1, iterlim=100,
#	typsiz=abs(p), stepmax=10*sqrt(p%*%p))
#
#  DESCRIPTION
#
#    Function to fit overdispersed count data as a birth process

pbirth <- function(frequencies, p, intensity="negative binomial",
	type="spectral decomposition", print.level=0, ndigit=10,
	gradtol=0.00001, steptol=0.00001, fscale=1, iterlim=100,
	typsiz=abs(p), stepmax=10*sqrt(p%*%p)){
call <- sys.call()
#
# find distribution and method
#
intensity <- match.arg(intensity,c("binomial","binomial exponential",
	"binomial logistic","binomial total","Poisson","Poisson exponential",
	"negative binomial","gen negative binomial"))
type <- match.arg(type,c("spectral decomposition","series approximation"))
#
# set up frequencies
#
if(is.vector(frequencies))frequencies <- matrix(frequencies,nrow=1)
n <- dim(frequencies)[2]-1
nr <- dim(frequencies)[1]
np <- length(p)
n1 <- NULL
for(i in 1:nr)n1 <- c(n1,length(frequencies[i,!is.na(frequencies[i,])])-1)
#
# create appropriate intensity function
#
lambda <- switch(intensity,
	"binomial"=function(p,nn,n) (n-nn)*exp(p[1]),
	"binomial exponential"=function(p,nn,n) (n-nn)*exp(p[1]+p[2]*nn),
	"binomial logistic"=function(p,nn,n)
		(n-nn)*exp(p[1])/(1+exp(p[2]+p[3]*nn)),
	"binomial total"=function(p,nn,n)
		(n-nn)*exp(p[1]+p[4]*n)/(1+exp(p[2]+p[3]*nn)),
	"Poisson"=function(p,nn,n) rep(exp(p[1]),length(nn)),
	"Poisson exponential"=function(p,nn,n) exp(p[1]+p[2]*nn),
	"negative binomial"=function(p,nn,n) exp(p[1])*(exp(p[2])+nn),
	"gen negative binomial"=
		function(p,nn,n) exp(p[1])*(exp(p[2])+nn)^(1-exp(p[3])))
#
# create probability functions
#
x <- matrix(0,n+1,n+1)
i1 <- matrix(1:(n+1),n+1,2)
i2 <- cbind(1:n,2:(n+1))
prob <- function(p,n) {
	x[i1[1:(n+1),]] <- -lambda(p,0:n,n)
	x[i2[1:n,,drop=FALSE]] <- lambda(p,0:(n-1),n)
	pr <- mexp(x[1:(n+1),1:(n+1)],type=type)[1,]}
#	z <- eigen(x[1:(n+1),1:(n+1)],sym=FALSE)
#	pr <- (z$vectors%*%diag(exp(z$values))%*%solve(z$vectors))[1,]}
prob1 <- function(p){
	pr <- NULL
	for(i in 1:nr)pr <- c(pr,prob(p,n1[i]))}
#
# call nlm to optimize
#
freq <- as.vector(t(frequencies))
freq <- freq[!is.na(freq)]
like <- function(p) -sum(freq*log(prob1(p)),na.rm=TRUE)
z0 <- nlm(like,p=p, hessian=TRUE, print.level=print.level,
	typsiz=typsiz, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
#
# calculate se's
#
if(np==1)cov <- 1/z0$hessian
else {
	a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
		else qr(z0$hessian)$rank
	if(a==np)cov <- solve(z0$hessian)
	else cov <- matrix(NA,ncol=np,nrow=np)}
se <- sqrt(diag(cov))
nn <- sum(frequencies)
#
# back transform coefficients
#
p <- z0$estimate
pr <- prob(p,n)
pt <- vector(mode="double",np)
pr1 <- exp(-exp(p[1]))
if(intensity=="binomial")pt <- 1-pr1
else if(intensity=="Poisson")pt <- sum(pr*0:n)
else if(intensity=="negative binomial"||intensity=="gen negative binomial"){
	pr2 <- exp(p[2])
	pt[1] <- pr2/pr1-pr2
	pt[2] <- pr2
	if(intensity=="gen negative binomial")pt[3] <- 1-exp(p[3])}
else pt <- exp(-exp(p))
#
# create values from intensity function
#
an <- switch(intensity,
	"binomial"=rep(exp(p[1]),n+1),
	"binomial exponential"=exp(p[1]+p[2]*0:n),
	"binomial logistic"=exp(p[1])/(1+exp(p[2]+p[3]*0:n)),
	"binomial total"=exp(p[1]+p[4]*n)/(1+exp(p[2]+p[3]*0:n)),
	"Poisson"=rep(exp(p[1]),n+1),
	"Poisson exponential"=exp(p[1]+p[2]*0:n),
	"negative binomial"=exp(p[1])*(exp(p[2])+0:n),
	"gen negative binomial"=
		exp(p[1])*(exp(p[2])+0:n)^(1-exp(p[3])))
fitted.values <- nn*pr
residuals <- (frequencies-fitted.values)/sqrt(fitted.values)
z1 <- list(
	call=call,
	intensity=intensity,
	lambda=lambda,
	an=an,
	frequencies=frequencies,
	maxlike=z0$minimum,
	aic=z0$minimum+np,
	fitted.values=fitted.values,
	prob=pr,
	residuals=residuals,
	initial.values=p,
	coefficients=p,
	pt=pt,
	se=se,
	cov=cov,
	corr=cov/(se%o%se),
	gradient=z0$gradient,
	iterations=z0$iterations,
	error=z0$error,
	code=z0$code)
class(z1) <- "pbirth"
return(z1)}

### standard method
###

deviance.pbirth <- function(z) 2*z$maxlike

### print method
###
print.pbirth <- function(z){
np <- length(z$coefficients)
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
t <- deparse(z$lambda)
cat(z$intensity,"intensity function:",t[2:length(t)],sep="\n")
cat("-Log likelihood   ",z$maxlike,"\n")
cat("AIC               ",z$aic,"\n")
cat("Iterations        ",z$iterations,"\n\n")
cat("Coefficients:\n")
coef.table <- cbind(z$coefficients,z$se,z$pt)
dn <- paste("p",1:np,sep="")
dimnames(coef.table) <- list(dn,c("estimate","se","parameter"))
print.default(coef.table,digits=4,print.gap=2)
if(np>1){
	cat("\nCorrelations:\n")
	dimnames(z$corr) <- list(seq(1,np),seq(1,np))
	print.default(z$corr,digits=4)}
invisible(z)}
#
#  event : A Library of Special Functions for Event Histories
#  Copyright (C) 1998, 1999, 2000, 2001 J.K. Lindsey
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public Licence as published by
#  the Free Software Foundation; either version 2 of the Licence, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public Licence for more details.
#
#  You should have received a copy of the GNU General Public Licence
#  along with this program; if not, write to the Free Software
#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
#  SYNOPSIS
#
#  survkit(times, censor=NULL, ccov=NULL, tvcov=NULL,
#	strata=NULL, id=NULL, model="Weibull", baseline=FALSE,
#	residuals=FALSE, survival=NULL, svalues=NULL, valrho=NULL,
#	constraints=NULL, impose=NULL, dist=NULL, random=NULL,
#	estimate=NULL, moments=FALSE, rule=NULL, pedigree=NULL,
#	integrate=NULL, jointmode=FALSE, within=NULL, converge=1.e-8,
#	iterlim=100)
#
#  DESCRIPTION
#
#    Function to fit  Weibull and Cox models with frailties

.First.lib <- function(lib, pkg)
	library.dynam("event", pkg, lib)
require(rmutil)

survkit <- function(times, censor=NULL, ccov=NULL, tvcov=NULL,
	strata=NULL, id=NULL, model="Weibull", baseline=FALSE,
	residuals=FALSE, survival=NULL, svalues=NULL, valrho=NULL,
	constraints=NULL, impose=NULL, dist=NULL, random=NULL,
	estimate=NULL, moments=FALSE, rule=NULL, pedigree=NULL,
	integrate=NULL, jointmode=FALSE, within=NULL, converge=1.e-8,
	iterlim=100){
#
# fixed Fortran parameters
#
#-- MXSTRA= max. number of strata, i.e. levels of the strata variable
#--    defined in the STRATA statement.
mxstra <- 25
#-- NDIMAX = max. total number of levels of effects in the model
#--    (size of the vector of solutions)
ndimax <- 100
#-- MXEF = max. number of covariates. These are discrete
#--    (CLASS) covariates and continuous covariates.
mxef <- 30
#-- NSTIMAX = maximum number of distinct times or quantiles that
#--    can be defined in the SURVIVAL statement options SPECIFIC
#--    and QUANTILE.
nstimax <- 20
#-- NTIMMAX = largest possible value of the (dependent) time variable
ntimmax <- 6000
#-- EPS_BS = convergence criterion
eps.bf <- converge
#--    no_log = 1 -> rho is used in the maximisation routine
#--    no_log = 0 -> log(rho) is used in the maximisation routine
no.log <- 0
#-- ITER_BF = maximum number of iterations.
iter.bf <- iterlim
call <- sys.call()
model <- match.arg(model,c("Weibull","Cox","Kaplan.Meier"))
#
# time, censor, id
#
if(max(times)>ntimmax)stop(paste("maximum time must be <",ntimmax))
if(any(trunc(times)!=times))warning("some times are not integers")
if(any(times<=0))stop("times must be positive integers")
nrr <- length(times)
if(is.null(censor))censor <- rep(1,nrr)
else if(length(censor)!=nrr)stop("censor must have same length as times")
if(!is.null(id)){
	if(length(id)!=nrr)stop("id must have same length as times")
	if(is.factor(id))id <- as.numeric(as.numeric(id))
	nunique <- length(unique(id))}
else nunique <- nrr
info <- cbind(times,censor)
squant <- xtabl <- itabl <- fact <- cont <- xlabelc <- xtabl2 <-
	itabl2 <- info2 <- xlabelf <- xlabelrg <- ind1 <- ind2 <-
	xlabelrn <- istime <- NULL
integam <- nwithin <- ninttdep <- nrule <- nsurv <- nrr2 <-
	nstime <- ncons <- nrelmat <- nquant <- npest <-
	xinclu <- nall.tim <- nccov <- ndcov <- ncoef <- 0
#
# covariates
# model : assume all variables in model
#
if(!is.null(ccov)){
	if(!inherits(ccov,"formula"))stop("ccov must be a model formula")
	mf <- model.frame(terms(ccov),parent.frame(),
		na.action=na.fail)
	if(dim(mf)[1]!=nunique&&dim(mf)[1]!=nrr)
		stop("time-constant covariates must have one value per individual or one per time")
	for(i in 1:dim(mf)[2]){
		if(nunique!=nrr&&dim(mf)[1]<nrr)mf[,i] <- mf[,i][id]
		if(is.factor(mf[,i])){
			xlabelf <- c(xlabelf,colnames(mf)[i])
			ndcov <- ndcov+1
			fact <- c(fact,length(levels(mf[,i])))
			itabl <- cbind(itabl,as.numeric(mf[,i]))}
		else if(is.numeric(mf[,i])){
			xlabelc <- c(xlabelc,colnames(mf)[i])
			nccov <- nccov+1
			cont <- c(cont,1)
			xtabl <- cbind(xtabl,as.numeric(mf[,i]))}}}
includ1 <- rep(1,nccov)
includ2 <- rep(1,ndcov)
ind1 <- cont
ind2 <- fact
if(model!="Weibull"){
	iplus1 <- rep(0,nccov)
	iplus2 <- rep(0,ndcov)}
if(!is.null(tvcov)){
	if(!inherits(tvcov,"formula"))stop("tvcov must be a model formula")
	mf <- model.frame(terms(tvcov),parent.frame(),
		na.action=na.fail)
	if(dim(mf)[1]!=nrr)
		stop("time-constant covariates must have one value per time")
	for(i in 1:dim(mf)[2]){
		if(is.factor(mf[,i])){
			xlabelf <- c(xlabelf,colnames(mf)[i])
			fact <- c(fact,length(levels(mf[,i])))
			ndcov <- ndcov+1
			includ2 <- c(includ2,1)
			itabl <- cbind(itabl,as.numeric(mf[,i]))
			if(model=="Cox"){
				iplus2 <- c(iplus2,1,-1)
				ndcov <- ndcov+1
				ind2 <- c(ind2,0,length(levels(mf[,i])))
				includ2 <- c(includ2,0)
				itabl <- cbind(itabl,if(censor!=-1)0
					else c(as.numeric(mf[,i])[2:nrr],0))}
			else ind2 <- c(ind2,length(levels(mf[,i])))}
		else if(is.numeric(mf[,i])){
			xlabelc <- c(xlabelc,colnames(mf)[i])
			cont <- c(cont,1)
			nccov <- nccov+1
			includ1 <- c(includ1,1)
			xtabl <- cbind(xtabl,as.numeric(mf[,i]))
			if(model=="Cox"){
				iplus1 <- c(iplus1,1,-1)
				nccov <- nccov+1
				ind1 <- c(ind1,0,1)
				includ1 <- c(includ1,0)
				xtabl <- cbind(xtabl,if(censor!=-1)0
					else c(as.numeric(mf[,i])[2:nrr],0))}
			else ind1 <- c(ind1,1)}}}

ind <- c(0,cumsum(c(ind1,ind2)))
includ <- c(includ1,includ2)
xnor <- xgam <- rep(0,nccov+ndcov+1)
inor <- igam <- matrix(0,ncol=3,nrow=nccov+ndcov+1)
if(nccov+ndcov>0)bound <- matrix(0,nrow=nccov+ndcov,ncol=3)
else bound <- NULL
#
# constraints
#
icons <- rep(0,ndimax+1)
if(!is.null(constraints)&&is.null(impose)){
	tmp <- c("none","find","largest","first")
	ncons <- match(constraints <- match.arg(constraints,tmp),tmp)-3
	if(is.na(ncons))stop("unknown constraint")
	if(ncons==1){
		icons[1] <- length(fact)
		if(icons[1]==1) icons[2] <- length(ind1)+1
		else icons[2:(length(fact)+1)] <- cumsum(c(0,fact[2:length(fact)]))+length(ind1)+1}}
#
# impose is a list of a vector of variable names and a vector of their
#	constraint levels
#
if(!is.null(impose)){
	ncons <- 1
	icons[1] <- length(fact)
	if(icons[1]==1) icons[2] <- length(ind1)+1
	else icons[2:(length(fact)+1)] <- cumsum(c(0,fact[2:length(fact)]))+length(ind1)+1
	for(i in 1:length(impose[[1]])){
		tmp <- match(impose[[1]][i],xlabelf)
		if(is.na(tmp))stop("imposed constraint variable not found")
		icons[1] <- icons[1]+1
		icons[tmp+1] <- ind[tmp+length(ind1)]+impose[[2]][i]}}
#
# baseline, kaplan, and residuals  for Cox
nbase <- as.integer(baseline)
nkaplan <- as.integer(model=="Kaplan.Meier")
if(residuals){
	if(model=="Cox")nsurv <- 2
	else stop("residuals only available for Cox model")}
#
# random requires dist, estimate, and maybe rule
#
if(!is.null(random)){
	if(!is.factor(random))stop("random must be a factor variable")
	if(length(random)!=nrr)
		stop("random covariates must have one value per time")
	if(is.null(dist))
		stop("a dist must be provided for the random covariate")
	if(is.null(estimate))
		stop("estimates must be provided for the random covariate")
	ndcov <- ndcov+1
	itabl <- cbind(itabl,as.numeric(random))
	includ <- c(includ,1)
	if(is.na(match(dist,c("loggamma","normal","multivariate"))))
		stop("unknown dist in random")
	xlabelf <- c(xlabelf,paste(deparse(substitute(random))))
	fact <- c(fact,length(levels(random)))
	if(dist=="loggamma"){
		igam[1,1] <- 1
		colgam <- ndcov
		xlabelrg <- c(xlabelrg,paste(deparse(substitute(random))))
		igam <- rbind(igam,c(1,0,0))
		inor <- rbind(inor,c(0,0,0))}
	else if(dist=="normal"){
		inor[1,1] <- 1
		colnor <- ndcov
		xlabelrn <- c(xlabelrn,paste(deparse(substitute(random))))
		igam <- rbind(igam,c(0,0,0))
		inor <- rbind(inor,c(1,0,0))}
	else if(dist=="multivariate"){
		if(is.null(pedigree)||!is.matrix(pedigree))
			stop("a pedigree matrix must be supplied with multivariate random effects")
		if(dim(pedigree)[2]!=4)
			stop("the pedigree matrix must have 4 columns: id, sex, father, mother")
		pedigree <- t(rbind(pedigree,matrix(0,ncol=4,nrow=ndimax-dim(pedigree)[1])))
		if(is.null(rule))rule <- "usual"
		tmp <- c("usual","mgs","sire.dam")
		nrule <- match(rule <- match(dist,tmp),tmp)
		nrelmat <- nccov+ndcov}
	if(length(estimate)==3){
		if(moments)npest <- 100
		else npest <- 1
		if(estimate[1]<=0)
			stop("lower estimate bound negative")
		else if(estimate[1]>=estimate[2])
			stop("incorrect estimate bounds")
		bound <- rbind(bound,estimate)
		if(dist=="loggamma"){
			xgam <- c(xgam,estimate[2])
			xnor <- c(xnor,0)}
		else {
			xgam <- c(xgam,0)
			xnor <- c(xnor,estimate[1])}}
	else {
		if(dist=="loggamma"){
			xgam <- c(xgam,estimate)
			xnor <- c(xnor,0)}
		else if(dist=="normal"){
			xnor <- c(xnor,estimate)
			xgam <- c(xgam,0)}
		bound <- rbind(bound,rep(0,3))}}
#
# integrate
#
if(!is.null(integrate)){
	if(model=="Cox")stop("integrate only possible with Weibull model")
	if(!is.factor(integrate))stop("integrate must be a factor variable")
	if(length(integrate)!=nrr)
		stop("integrate must be the same length as other variables")
	integam <- nccov+ndcov+1
	itabl <- cbind(itabl,as.numeric(integrate))
	includ <- c(includ,0)
	xlabelrg <- c(xlabelrg,paste(deparse(substitute(integrate))))
	if(!is.null(xlabelf)&&any(xlabelrg[length(xlabelrg)]==xlabelf))
		stop("the integrate variable must not be used elsewhere in the model")
	igam[1,1] <- 1
	igam <- rbind(igam,c(1,0,0))
	inor <- rbind(inor,c(0,0,0))
	if(is.null(estimate)||length(estimate)!=1)
		stop("one estimate must be supplied with integrate")
	xgam <- c(xgam,estimate)
	xnor <- c(xnor,0)
	bound <- rbind(bound,rep(0,3))
	if(!is.null(within)){
		if(!is.factor(within))stop("within must be a factor variable")
		if(nunique!=nrr&&length(within)<nrr)
			within <- within[id]
		nwithin <- 1}}
else jointmode <- FALSE
#
# strata
#
if(!is.null(strata)){
	if(!is.factor(strata))stop("strata must be a factor variable")
	else if(length(strata)!=nrr)
		stop("strata variable must have the same length as other variables")
	nstrata <- length(levels(strata))
	if(nstrata>mxstra)stop(paste("maximum",mxstra,"strata allowed"))
	info <- cbind(info,as.numeric(strata))}
else {
	nstrata <- 1
	if(nwithin==1)info <- cbind(info,as.numeric(within))
	else info <- cbind(info,as.numeric(rep(1,nrr)))}

if(model=="Weibull"&&integam>0)info <- cbind(info,itabl[,integam-nccov])
else if(is.null(id))info <- cbind(info,1:nrr)
else info <- cbind(info,id)
nind <- length(unique(info[,4]))
#
# rho_fixed
#
if(model=="Cox")icrho <- 0
else if(!is.null(valrho)){
	icrho <- 0
	if(is.vector(valrho,mode="numeric")&&length(valrho)!=nstrata)
		stop("valrho must have a value for each	stratum")
	valrho <- c(valrho,rep(0,mxstra-nstrata))}
else {
	icrho <- 1
	valrho <- if(no.log==0)rep(1,mxstra)
		else rep(0,mxstra)}
#
# survival
#
if(!is.null(survival)){
	if(nsurv==2)stop("survival and residuals cannot both be chosen")
	nsurv <- 1
	survival <- match.arg(survival,c("equal","quantiles","specific","all.times"))
	if(survival!="all.times"){
		if(is.null(svalues))
			stop("svalues must be supplied with survival")
		if(!is.vector(svalues,mode="numeric"))
			stop("svalues must be a vector")
		if(length(svalues)>nstimax)stop("too many svalues")
		if(survival!="equal"&&any(diff(svalues)<=0))
			stop("svalues must be strictly increasing")}
	else if(model!="Cox")
		stop("all.times survival only valid with Cox model")
	if(survival=="quantiles"){
		if(any(svalues<0|svalues>100))
			stop("quantiles must be between 0 and 100")
		squant <- 0.01*svalues
		nquant <- length(squant)}
	else if(survival=="equal"){
		if(length(svalues)!=2)stop("2 svalues must be given")
		else if(svalues[1]>=svalues[2])
			stop("time interval must be smaller than time limit")
		istime <- seq(svalues[1],svalues[2],by=svalues[1])
		if(length(istime)>nstimax)
			stop("time interval too small")
		else nstime <- length(istime)}
	else if(survival=="all.times")nall.tim <- 1
	else {
		if(any(svalues<=0))stop("svalues must be positive")
		istime <- svalues
		if(length(istime)>nstimax)
			stop("time interval too small")
		else nstime <- length(istime)}}

if(nsurv>0){
	info2 <- info
	itabl2 <- itabl
	xtabl2 <- xtabl
	nrr2 <- nrr}
if(model=="Weibull"){
	if(nstrata>1||integam>0){
		if(!is.null(integrate)&&nwithin==0)
			o <- order(itabl[,integam-nccov],info[,4],times)
			else o <- order(info[,3],info[,4],times)}
	else o <- order(info[,4],times)}
else {
	if(nstrata>1)o <- order(info[,3],-times,censor)
	else o <- order(-times,censor)}
info <- info[o,]
itabl <- itabl[o,,drop=FALSE]
xtabl<- xtabl[o,]
#
# call Fortran
#
iconst <- c(ncons,nsurv,nstrata,no.log,nrelmat,integam,nwithin,
	nccov+sum(includ),nall.tim,nbase,icrho,ninttdep,jointmode,
	npest,model=="Cox",nrule,nkaplan,iter.bf)
rconst <- c(xinclu,eps.bf)
if(model=="Cox"&&length(iplus2)==0)iplus2 <- 0
if(is.null(istime))istime <- 0
z0 <- if(model=="Weibull")
	.Fortran("weibull",
	iconst=as.integer(iconst),
	rconst=as.double(rconst),
	istime=as.integer(istime),
	squant=as.double(squant),
	valrho=as.double(valrho),
	includ=as.integer(includ),
	info=as.integer(info),
	itabl=as.integer(itabl),
	xtabl=as.double(xtabl),
	nrr=as.integer(nrr),
	info2=as.integer(info2),
	itabl2=as.integer(itabl2),
	xtabl2=as.double(xtabl2),
	ipedig=as.integer(pedigree),
	nrr2=as.integer(nrr2),
	nccov=as.integer(nccov),
	ndcov=as.integer(ndcov),
	nstime=as.integer(nstime),
	nquant=as.integer(nquant),
	igam=as.integer(igam),
	inor=as.integer(inor),
	xgam=as.double(xgam),
	xnor=as.double(xnor),
	bound=as.double(bound),
	icons=as.integer(icons),
	beta=double(ndimax),
	stder=double(ndimax),
	grad=double(ndimax),
	hess=double(ndimax*ndimax),
	xmom=double(4),
	surv=double(3*nstimax*nrr),
	like=double(2),
	df=integer(2),
	DUP=FALSE,
	PACKAGE="event")
else
	.Fortran("cox",
	iconst=as.integer(iconst),
	rconst=as.double(rconst),
	istime=as.integer(istime),
	squant=as.double(squant),
	includ=as.integer(includ),
	info=as.integer(info),
	itabl=as.integer(itabl),
	xtabl=as.double(xtabl),
	nrr=as.integer(nrr),
	info2=as.integer(info2),
	itabl2=as.integer(itabl2),
	xtabl2=as.double(xtabl2),
	iplus1=as.integer(iplus1),
	iplus2=as.integer(iplus2),
	ipedig=as.integer(pedigree),
	nrr2=as.integer(nrr2),
	nccov=as.integer(nccov),
	ndcov=as.integer(ndcov),
	nstime=as.integer(nstime),
	nquant=as.integer(nquant),
	igam=as.integer(igam),
	inor=as.integer(inor),
	xgam=as.double(xgam),
	xnor=as.double(xnor),
	bound=as.double(bound),
	icons=as.integer(icons),
	beta=double(ndimax),
	stder=double(ndimax),
	grad=double(ndimax),
	hess=double(ndimax*ndimax),
	xmom=double(4),
	surv=double(3*nstimax*nrr),
	like=double(2),
	df=integer(2),
	km=double(nrr*nstrata*4),
	resid=double(nrr*nstrata*3),
	ut=as.integer(nrr*nstrata),
	DUP=FALSE,
	PACKAGE="event")
if(z0$iconst[1]>0)switch(as.character(z0$iconst[1]),
	"1"=stop("too many error constraints"),
	"2"=stop("negative value on the diagonal of the Hessian"),
	"3"=stop("algorithm diverged"),
	"4"=stop("survival curve changes too many times"),
	"5"=stop("matrix is not semi-positive definite"),
	"6"=stop("too many parameters"))
code <- z0$iconst[3]>=iterlim+2*(z0$iconst[1]==11)
if(nkaplan>0){
	km <- matrix(z0$km,ncol=4)
	rownames(km) <- km[,1]
	km <- km[km[,1]>0,2:4]
	colnames(km) <- c("survival function","log variance","cumulative hazard")
	km <- list(km=km,nstrata=nstrata)
	class(km) <- c("survivalkit",model)
	return(km)}
else if(model=="Cox"){
	if(baseline){
		base <- matrix(z0$km,ncol=4)
		rownames(base) <- base[,1]
		base <- base[base[,1]>0,2:3]
		colnames(base) <- c("cumulative hazard","survival function")}
	else base <- NULL
	if(residuals){
		o <- order(info2[,1])
		info2 <- info2[o,]
		resid <- matrix(z0$resid,ncol=3)
		o <- resid[,1]>0
		resid <- resid[resid[,1]>0,1:3]
		resid <- cbind(info2[info2[,2]>-1&o,1:3],resid)
		colnames(resid) <- c("time","censor","stratum","generalized","martingale","deviance")}
	else resid <- NULL}
else resid <- base <- km <- NULL
if(nsurv==1){
	surv <- matrix(z0$surv,ncol=3)
	surv <- surv[surv[,1]>0,]
	tmp <- surv[,1]
	surv <- surv[,2:3]
	colnames(surv) <- c("time","S(t)")
	rownames(surv) <- tmp}
else surv <- NULL
ndim <- z0$iconst[2]-nstrata*(model=="Weibull")-jointmode
ndimf <- z0$iconst[2]
if(!is.null(xlabelc)||!is.null(xlabelf)){
	xlabel <- c(xlabelc,xlabelf)
	xlabel <- xlabel[rep(1:length(xlabel),c(cont,fact))]}
else xlabel <- NULL
cov <- matrix(z0$hess,ncol=ndimax)[1:ndimf,1:ndimf,drop=FALSE]
if(model=="Weibull"){
	rho <- z0$beta[(ndim+1):(ndimf-jointmode)]
	rhose <- z0$stder[(ndim+1):(ndimf-jointmode)]
	tmp <- (ndim-nstrata+1):ndim
	if(ndim>nstrata)tmp <- c(tmp,1:(ndim-nstrata))
	beta <- z0$beta[tmp]
	se <- z0$stder[tmp]
	tmp <- c(tmp,(ndim+1):ndimf)
	cov <- cov[tmp,tmp,drop=FALSE]}
else {
	rho <- rhose <- NULL
	beta <- z0$beta[1:ndim]
	se <- z0$stder[1:ndim]}
if(jointmode){
	beta <- c(beta,z0$beta[ndimf])
	se <- c(se,z0$stder[ndimf])}
if(npest==100)moments <- z0$xmom
else moments <- NULL
if(!is.null(xlabelf)){
	if(model=="Weibull")
		xlabel[(length(xlabelc)+1):(ndim-nstrata)] <- paste(xlabel[(length(xlabelc)+1):(ndim-nstrata)],sequence(fact),sep="")
	else xlabel[(length(xlabelc)+1):ndim] <- paste(xlabel[(length(xlabelc)+1):ndim],sequence(fact),sep="")}
if(jointmode)xlabel <- c(xlabel,integrate)
ndim <- ndim+jointmode
tmp <- se[1:ndim]>0
beta <- beta[tmp]
se <- se[tmp]
if(!is.null(xlabel)){
	if(model=="Weibull")xlabel <- xlabel[tmp[(nstrata+1):ndim]]
	else xlabel <- xlabel[tmp]}
ndim <- sum(tmp)
if(model=="Weibull")tmp <- c(tmp,rep(TRUE,nstrata))
if(ndim>0)cov <- cov[tmp,tmp,drop=FALSE]
if(model=="Weibull"){
	df <- c(z0$df[1],ndim+nstrata)
	tmp <- "(Intercept)"
	if(nstrata>1)tmp <- paste(tmp,1:nstrata,sep="")}
else {
	df <- c(z0$df[1],ndim)
	tmp <- NULL}
blabel <- c(tmp,xlabel)
betag <- betan <- seg <- sen <- blabelg <- blabeln <- NULL
ndimg <- ndimn <- 0
if(integam>0&&jointmode){
	tmp <- ndim
	betag <- beta[ndim]
	seg <- se[ndim]
	blabelg <- blabel[ndim]
	ndimg <- 1}
else if(igam[1,1]>0&&integam==0){
	tmp <- (ndim-length(unique(itabl[,colgam]))+1):ndim
	betag <- beta[tmp]
	seg <- se[tmp]
	blabelg <- blabel[tmp]
	ndimg <- length(tmp)}
else if(inor[1,1]>0){
	tmp <- (ndim-length(unique(itabl[,colnor]))+1):ndim
	betan <- beta[tmp]
	sen <- se[tmp]
	blabeln <- blabel[tmp]
	ndimn <- length(tmp)}
if((igam[1,1]>0&&integam==0)||inor[1,1]>0||(integam>0&&jointmode)){
	tmp <- 1:(tmp[1]-1)
	beta <- beta[tmp]
	se <- se[tmp]
	blabel <- blabel[tmp]
	ndim <- ndim-ndimg-ndimn}
z <- list(
	ccov=ccov,
	tvcov=tvcov,
	beta=beta,
	betag=betag,
	betan=betan,
	rho=rho,
	iterations=z0$iconst[3],
	igam=igam,
	xgam=z0$xgam,
	inor=inor,
	xnor=z0$xnor,
	integam=integam,
	moments=moments,
	ndim=ndim,
	ndimg=ndimg,
	ndimn=ndimn,
	ncov=nccov+ndcov+nstrata-inor[1,1]-igam[1,1]*(integam==0),
	nind=nind,
	nstrata=nstrata,
	jointmode=jointmode,
	blabel=blabel,
	blabelg=blabelg,
	blabeln=blabeln,
	xlabelrg=xlabelrg,
	xlabelrn=xlabelrn,
	se=se,
	seg=seg,
	sen=sen,
	rhose=rhose,
	like=z0$like,
	nrr=nrr,
	nevents=sum(censor==1),
	ncensor=sum(censor==0),
	df=df,
	grad=z0$grad,
	cov=cov,
	baseline=base,
	residuals=resid,
	survival=surv,
	call=call,
	code=code)
class(z) <- c("survivalkit",model)
z}

### standard methods
###
residuals.survivalkit <- function(z){
if(inherits(z,"Cox")){
	if(!is.null(z$residuals))print(z$residuals)
	else print("residuals not available: use residuals=TRUE")}
else print("residuals only available for Cox model")}

print.survivalkit <- function(z){
if(inherits(z,"Kaplan.Meier")){
	cat("\nKaplan-Meier estimates with 95% interval and Nelson estimates\n\n")
	u <- log(z$km[,1])
	r <- sqrt(z$km[,2])/abs(u)
	tab <- cbind(z$km[,1],exp(exp(1.96*r)*u),exp(exp(-1.96*r)*u),z$km[,3])
	tmp <- c(colnames(z$km)[1],"lower","upper",colnames(z$km)[3])
	if(z$nstrata>1){
		tmp <- c(tmp,"stratum")
		strata <- cumsum(z$km[,1]==0)+1
		strata <- c(1,strata[1:(length(strata)-1)])
		tab <- cbind(tab,strata)}
	colnames(tab) <- tmp
	print(tab)
	return(tab)}
else if(inherits(z,"Weibull")){
	mdl <- 1
	cat("\nWeibull model\n\n")}
else {
	mdl <- 2
	cat("\nCox model\n\n")}
cat("Call:\n",deparse(z$call),"\n\n",sep="")
if(z$code%%2>0)warning("algorithm did not converge")
if(z$code>1)warning("no convergence for Gauss-Hermite integration")
if(z$code>0)cat("\n")
if(z$ndimg==0&&z$ndimn==0){
	cat("-Log likelihood (null)",z$like[1],"\n")
	cat("AIC (null)            ",z$like[1]+z$df[1],"\n")}
cat("-Log likelihood       ",z$like[2],"\n")
cat("AIC                   ",z$like[2]+z$df[2],"\n")
cat("Number of parameters  ",z$df[2],"\n")
if(z$nind<z$nevents+z$ncensor)cat("Number of individuals ",z$nind,"\n")
cat("Number of events      ",z$nevents,"\n")
if(z$ncensor>0)cat("Number censored       ",z$ncensor,"\n")
if(z$nstrata>1)cat("Number of strata      ",z$nstrata,"\n")
cat("Number of iterations  ",z$iter,"\n")
if(z$igam[1,1]>0&&!z$jointmode){
	cat("\nLog gamma random variance (mode)\n")
	tmp <- NULL
	for(i in 1:z$igam[1,1])tmp <- c(tmp,z$xgam[z$ncov-z$jointmode+i])
	names(tmp) <- z$xlabelrg
	print.default(tmp, digits=4, print.gap=2)}
if(z$inor[1,1]>0){
	cat("\nNormal random variance (mode)\n")
	tmp <- NULL
	for(i in 1:z$inor[1,1])tmp <- c(tmp,z$xnor[z$ncov+i])
	names(tmp) <- z$xlabelrn
	print.default(tmp, digits=4, print.gap=2)}
if(!is.null(z$moments)){
	tmp <- matrix(z$moments[2:4],nrow=1)
	rownames(tmp) <- "moments"
	colnames(tmp) <- c("mean","sd","skew")
	print.default(tmp, digits=4, print.gap=2)}
if(z$ndim>0){
	cat("\nRegression coefficients\n")
	coef.table <- cbind(z$beta,z$se)
	if(z$ndim==1)coef.table <- matrix(coef.table,ncol=2)
	colnames(coef.table) <- c("estimate","se")
	rownames(coef.table) <- z$blabel
	print.default(coef.table, digits=4, print.gap=2)}
if(z$ndimg>0){
	if(!z$jointmode)cat("\nLog gamma random effects\n")
	else cat("\nLog gamma random variance\n")
	coef.table <- cbind(z$betag,z$seg)
	if(z$ndimg==1)coef.table <- matrix(coef.table,ncol=2)
	colnames(coef.table) <- c("estimate","se")
	rownames(coef.table) <- z$blabelg
	print.default(coef.table, digits=4, print.gap=2)}
else if(z$ndimn>0){
	cat("\nNormal random effects\n")
	coef.table <- cbind(z$betan,z$sen)
	colnames(coef.table) <- c("estimate","se")
	rownames(coef.table) <- z$blabeln
	print.default(coef.table, digits=4, print.gap=2)}
if(inherits(z,"Weibull")){
	cat("\nWeibull power parameter\n")
	coef.table <- cbind(z$rho,z$rhose)
	colnames(coef.table) <- c("estimate","se")
	rownames(coef.table) <- paste(1:z$nstrata)
	print.default(coef.table, digits=4, print.gap=2)}
if((inherits(z,"Weibull")&&z$ndim+z$nstrata>1)||
	(any(class(z)=="Cox")&&z$ndim>1)){
	cat("\nCorrelations:\n")
	tmp <- c(z$se,z$seg,z$sen,z$rhose)
	corr <- z$cov/(tmp%o%tmp)
	if(inherits(z,"Cox"))tmp <- seq(1,z$ndim+z$ndimg+z$ndimn)
	else tmp <- seq(1,z$ndim+z$ndimg+z$ndimn+z$nstrata)
	dimnames(corr) <- list(tmp,tmp)
	print.default(corr, digits=4)}}

### special methods
###
baseline <- function(z, ...) UseMethod("baseline")

baseline.survivalkit <- function(z){
if(inherits(z,"Cox")){
	if(!is.null(z$baseline))print(z$baseline)
	else print("baseline not available: use baseline=TRUE")}
else print("baseline only available for Cox model")}

survival <- function(z, ...) UseMethod("survival")

survival.survivalkit <- function(z){
if(!is.null(z$survival))print(z$survival)
else print("survival not available: use survival option")}
