.packageName <- "growth"
#
#  growth : A Library of Normal Distribution Growth Curve Models
#  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
#
#     pergram(y)
#     corgram(y, wt=1, add=FALSE, lty=1, xlim=NULL, ylim=NULL, xlab=NULL,
#		ylab=NULL, main=NULL)
#
#  DESCRIPTION
#
#    Functions to compute and plot correlograms and periodograms

### periodogram functions
###
pergram <- function(y){
	ll <- length(y)
	len <- trunc(ll/2)
	len1 <- if(ll==2*len)len+1 else len+2
	fc <- fft(y)
	z <- cbind(2*pi*(1:len)/ll,
		(Re(fc[ll:len1])^2+Im(fc[ll:len1])^2)*ll/4/pi)
	class(z) <- "pergram"
	z}

#pergram <- function(y){
#	ll <- length(y)
#	len <- trunc(ll/2)
#	fc <- matrix(0,ncol=2,nrow=ll)
#	for(i in 1:ll){
#		a <- 2*pi*i*(0:(ll-1))/ll
#		fc[i,1] <- sum(y*sin(a))
#		fc[i,2] <- sum(y*cos(a))}
#	invisible(cbind(2*pi*(1:len)/ll,(fc[,1]^2+fc[,2]^2)*ll/4/pi))}

plot.pergram <- function(y, add=FALSE, lty=1, xlab="Frequency",
	ylab="Periodogram", main="Periodogram", ylim=c(0,max(po[,2]))){
	if(inherits(y,"pergram"))po <- y
	else po <- pergram(y)
	if(add)lines(po[,1],po[,2],lty=lty)
	else plot(po[,1],po[,2],xlim=c(0,3.6),xlab=xlab,ylab=ylab,type="l",
		lty=lty,ylim=ylim)
	invisible(po)}

plot.cum <- function(z, ...) 
	UseMethod("plot.cum")

plot.cum.pergram <- function(y, xlab="Frequency", ylab="Periodogram",
	main="Cumulative periodogram",
	ylim=c(0,max(cpo+1.358/(a+0.12+0.11/a)))){
	if(inherits(y,"pergram")){
		len <- 2*dim(y)[1]
		po <- y}
	else {
		len <- length(y)
		po <- pergram(y)}
	cpo <- cumsum(po[,2])
	cpo <- cpo/max(cpo)
	a <- sqrt(len)
	pa <- 2*pi*(1:length(cpo))/len
	x <- (1:dim(po)[1])/dim(po)[1]
	plot(po[,1],cpo,xlim=c(0,3.6),xlab=xlab,ylab=ylab,type="l",
		ylim=ylim,main=main)
	lines(po[,1],x+1.358/(a+0.12+0.11/a),lty=3)
	lines(po[,1],x-1.358/(a+0.12+0.11/a),lty=3)}
###
### correlogram function
###
corgram <- function(y, wt=1, maxlag=NULL, partial=FALSE, add=FALSE, lty=1,
	xlim=NULL, ylim=NULL, xlab=NULL, ylab=NULL, main=NULL, ...){
	if(any(wt!=0&&wt!=1))stop("weigts must be 0 or 1")
	len <- length(y)
	if(length(wt)==1)wt <- rep(1,len)
	if(any(is.na(y)))wt[is.na(y)] <- 0
	wt <- wt>0
	num <- sum(wt)
	limit <- 2/sqrt(num)
	if(is.null(maxlag))maxlag <- trunc(num/4)
	ave <- sum(y*wt,na.rm=TRUE)/num
	var <- sum((y-ave)^2*wt,na.rm=TRUE)/(num-1)
	co <- rep(1,maxlag+1)
	for(lag in 1:maxlag)
		co[lag+1] <- sum((y[1:(len-lag)]-ave)*(y[(lag+1):len]-ave)*
			wt[1:(len-lag)]*wt[(lag+1):len],na.rm=TRUE)/
			var/(sum(wt[1:(len-lag)]*wt[(lag+1):len],na.rm=TRUE)-1)
	if(partial){
		a <- tmp <- rep(0,maxlag)
		a[1] <- tmp[1] <- co[2]
		for(lag in 2:maxlag){
			a[lag] <- tmp[lag] <- (co[lag+1]-sum(a[1:(lag-1)]*
				co[lag:2]))/(1-sum(a[1:(lag-1)]*co[2:lag]))
			a[1:(lag-1)] <- a[1:(lag-1)]-a[lag]*a[(lag-1):1]}
		co <- tmp}
	start <- if(partial)1 else 0
	if(is.null(xlim))xlim <- c(0,maxlag)
	if(is.null(ylim))ylim <- c(if(min(co)<0)min(c(min(co),-limit))else 0,1)
	if(is.null(xlab))xlab <- "Lag"
	if(is.null(ylab))ylab <- expression(rho)
	if(is.null(main))main <- if(partial)"Partial Correlogram"
		else"Correlogram"
	if(add)lines(start:maxlag,co,lty=lty)
	else {
		plot(start:maxlag,co,type="l",lty=lty,xlim=xlim,ylim=ylim,
			xlab=xlab,ylab=ylab,main=main,...)
		abline(h=limit,lty=3)
		if(min(co)<0)abline(h=-limit,lty=3)
	abline(h=0)}
	invisible(cbind(start:maxlag,co))}
#
#  growth : A Library of Normal Distribution Growth Curve Models
#  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
#
#     carma(response=NULL, ccov=NULL, times=NULL, torder=0, interaction,
#	arma=c(0,0,0), parma=NULL, pre=NULL, position=NULL, iopt=TRUE,
#	resid=TRUE, transform="identity", delta=NULL, envir=parent.frame(),
#	print.level=0, typsiz=abs(p), ndigit=10, gradtol=0.00001,
#	steptol=0.00001, iterlim=100, fscale=1, stepmax=10*sqrt(p%*%p))
#
#  DESCRIPTION
#
#    Function to fit the multivariate normal distribution with ARMA
#  and random intercept using Kalman filtering in continuous time

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

carma <- function(response=NULL, ccov=NULL, times=NULL, torder=0, interaction,
	arma=c(0,0,0), parma=NULL, pre=NULL, position=NULL, iopt=TRUE,
	resid=TRUE, transform="identity", delta=NULL, envir=parent.frame(),
	print.level=0, typsiz=abs(p), ndigit=10, gradtol=0.00001,
	steptol=0.00001, iterlim=100, fscale=1, stepmax=10*sqrt(p%*%p)){
#
# likelihood function returning parameter values
#
kalman <- function(p){
	z <- .Fortran("kalman",
		np=as.integer(length(p)),
		par=as.double(p),
		like=double(1),
		xx=matrix(0,nrow=nlp+1,ncol=nlp+1),
		y=as.double(y),
		sse=double(1),
		nq=as.integer(nre),
		nlp=as.integer(nlp),
		ns=as.integer(ns),
		nt=as.integer(nt),
		model=as.integer(arma),
		t=as.double(times),
		nobs=as.integer(nobs(response)),
		nod=as.integer(nod),
		as.integer(position),
		cv=as.double(ccov),
		ncv=as.integer(ncv),
		nxcv=as.integer(interaction),
		nx=as.integer(torder),
		p=double(length(p)),
		x=double(nlp+1),
		state=double(maxre*(nlp+1)),
		innov=double(nlp+1),
		cstate=complex(maxar*(nlp+1)),
		exx=complex(nlp+1),
		DUP=FALSE,
		PACKAGE="growth")
	list(like=z$like/2,
		sse=z$sse,
		xx=z$xx)}
#
# likelihood function for nlm
#
kalmanl <- function(p){
	z <- .Fortran("kalman",
		np=as.integer(length(p)),
		par=as.double(p),
		like=double(1),
		xx=matrix(0,nrow=nlp+1,ncol=nlp+1),
		y=as.double(y),
		sse=double(1),
		nq=as.integer(nre),
		nlp=as.integer(nlp),
		ns=as.integer(ns),
		nt=as.integer(nt),
		model=as.integer(arma),
		t=as.double(times),
		nobs=as.integer(nobs(response)),
		nod=as.integer(nod),
		as.integer(position),
		cv=as.double(ccov),
		ncv=as.integer(ncv),
		nxcv=as.integer(interaction),
		nx=as.integer(torder),
		p=double(length(p)),
		x=double(nlp+1),
		state=double(maxre*(nlp+1)),
		innov=double(nlp+1),
		cstate=complex(maxar*(nlp+1)),
		exx=complex(nlp+1),
		DUP=FALSE,
		PACKAGE="growth")
	z$like/2}
#
# Fortran constants for dependence parameters
#
maxre <- 6
maxar <- 6
maxma <- 5
call <- sys.call()
#
# check transformation
#
if(!missing(transform))transform <- match.arg(transform,c("identity",
	"exp","square","sqrt","log"))
#
# if envir, remove extra (multivariate) responses
#
type <- "unknown"
respenv <- exists(deparse(substitute(response)),env=parent.frame())&&
	inherits(response,"repeated")&&!inherits(envir,"repeated")
if(!respenv&&inherits(envir,"repeated")){
	if(!is.null(envir$NAs)&&any(envir$NAs))
		stop("carma 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 response is not a data object, make one and handle covariates
#
if(inherits(response,"repeated")){
	if(dim(response$response$y)[2]>1)
		stop("carma only handles univariate responses")
	if(!is.null(response$NAs)&&any(response$NAs))
		stop("carma does not handle data with NAs")
	type <- response$response$type
	# only keep covariates used in model in data object
	if(is.null(ccov))response$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
		response$ccov$ccov <- tmp[,-1,drop=FALSE]
		rm(tmp)}
	else stop("ccov must be a W&R formula")
	ncv <- if(is.null(response$ccov$ccov)) 0
		 else  dim(response$ccov$ccov)[2]}
else {
	if(!inherits(response,"response"))
		response <- restovec(response,times,delta=delta)
	type <- response$type
	if(dim(response$y)[2]>1)stop("carma only handles univariate responses")
	if(!is.null(ccov)){
	# find covariates and put in data object
		if(!inherits(ccov,"tccov")){
			ccname <- deparse(substitute(ccov))
			if(ccname=="~1")stop("ccov must contain variables")
			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)}
		ncv <- dim(ccov$ccov)[2]}
	else ncv <- 0
	response <- rmna(response=response, ccov=ccov)
	if(!is.null(ccov))rm(ccov)}
if(type!="unknown"&&type!="continuous")stop("continuous data required")
if(!is.null(response$response$nest)&&arma[1]>0)
	stop("carma does not handle two levels of nesting with an AR")
if(!is.null(response$response$censor))stop("carma does not handle censoring")
y <- response$response$y
ns <- length(nobs(response))
if(!is.null(response$ccov))ccov <- t(response$ccov$ccov)
#
# check times
#
torder <- torder+1
if(!is.null(response$response$times)){
	if(torder>length(unique(response$response$times)))
		stop("torder is too large for the number of distinct times")
	ave <- mean(response$response$times)
	times <- response$response$times-ave}
else {
	if(any(arma>0)||!is.null(parma))stop("No times. ARMA cannot be fitted")
	else if(torder>1)stop("No times. Time trends cannot be fitted.")
	ave <- 0
	times <- 1:length(y)}
#
# check interactions of covariates with times
#
if(missing(interaction)) {
	if(ncv==0)interaction <- NULL
	else interaction <- rep(1,ncv)}
else if(length(interaction)!=ncv)
	stop(paste("The vector, interaction, must have length, ", ncv))
else interaction <- interaction+1
if(!is.null(interaction)&&any(interaction>torder))
	stop("Interactions cannot have higher order than torder")
#
# check arma
#
if(length(arma)!=3){
	if(length(arma)==1)arma <- c(arma,0,0)
	else stop("The model specification must have length 3")}
if(sum(arma)==0&&!missing(parma)){
	if(length(parma)==1)arma[1] <- 1
	else stop("If an ARMA is required, arma must be supplied")}
if(arma[1]>maxar)stop(paste("Maximum order of AR is ",maxar))
if(arma[2]>maxma)stop(paste("Maximum order of MA is ",maxma))
if(arma[2]>=arma[1]&&arma[2]!=0)
	stop(paste("Order of MA(",arma[2],") >= order of AR(",arma[1],")"))
if(arma[1]==0&&arma[3]!=0)
	stop("Cannot have observation errors without AR structure")
np <- arma[1]+arma[2]+arma[3]
p <- parma
if(np>0&&length(parma)!=np)
	stop(paste(np,"initial arma estimates must be supplied"))
#
# check random coefficients
#
if(!is.null(pre))pre[pre==0] <- 0.01
nre <- length(pre)
if(nre>maxre)stop(paste("Only ",maxre," random effects allowed"))
if(is.null(position)){
	if(nre==1)position <- c(1,1)
	else if(nre==2)position <- rbind(c(1,1),c(2,2))
	else if(nre>2)stop("If random effects are required, position must be supplied")}
if(!is.null(position)){
	if(!is.matrix(position))position <- t(as.matrix(position))
	if(dim(position)[1]!=nre)
		stop(paste("Random effects position must have",nre,"rows"))
	if(dim(position)[2]!=2)stop("Position matrix must have two columns")
	if((any(position[,1]<1)||any(position[,1]>nre)||any(position[,2]<1)||
		any(position[,2]>nre))||position[,1]>position[,2])
		stop("Position for covariance matrix out of range")
	if(max(position[,2])>torder)
		warning("number of random effects greater than order of polynomial in time")
	nod <- dim(position)[1]
	if(nod>nre*(nre+1)/2)
		stop(paste("Only ",nre*(nre+1)/2," elements allowed in the random effects covariance matrix"))
	else if(nod<nre)stop("Not enough initial estimates of random effects")
	p <- c(p,pre)
	np <- length(p)}
else {
	position <- NULL
	nod <- 0}
nt <- sum(nobs(response))
nlp <- torder
if(ncv>0) nlp <- nlp+sum(interaction)
#
# calculate transformation and Jacobian
#
if(transform=="identity")jacob <- 0
else if(transform=="exp"){
	jacob <- -sum(y)
	y <- exp(y)}
else if(any(y==0))stop("Zero response values: invalid transformation")
else if(transform=="square"){
	jacob <- -sum(log(abs(y)))-length(resp$response$y)*log(2)
	y  <- y^2}
else if(any(y<0))stop("Nonpositive response values: invalid transformation")
else if(transform=="sqrt"){
	jacob <- sum(log(y))/2+length(resp$response$y)*log(2)
	y <- sqrt(y)}
else if(transform=="log"){
	jacob <- sum(log(y))
	y <- log(y)}
#
# include delta in Jacobian
#
if(!is.null(response$response$delta)){
	if(length(response$response$delta)==1)jacob <- jacob-nt*log(response$response$delta)
	else jacob <- jacob -sum(log(response$response$delta))}
#
# estimate model
#
if(iopt&&np>0){
	if(fscale==1)fscale <- kalmanl(p)
	z0 <- nlm(kalmanl, 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}
z2 <- kalman(p)
like <- z2$like+jacob
sse <- z2$sse
mse <- sse/(nt-nlp-np)
aic <- like+np+nlp+1
#
# get beta coefficients
#
z1 <- .Fortran("back",
	xx=z2$xx,
	as.integer(nlp),
	DUP=FALSE,
	PACKAGE="growth")
beta <- z1$xx[1:nlp,nlp+1]
#
# calculate se's
#
z1 <- .Fortran("ttvert",
	xx=z1$xx,
	as.integer(nlp),
	DUP=FALSE,
	PACKAGE="growth")
vbeta <- z1$xx[1:nlp,1:nlp]*mse
if(nlp>1){
	betase <- sqrt(diag(vbeta))
	corr <- vbeta/(betase%o%betase)}
else {
     betase <- sqrt(vbeta)
     corr <- vbeta/betase^2}
if(np>0){
	if(np==1){
		nlcov <- 1/z0$hessian
		nlse <- sqrt(nlcov)}
	else {
		a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
			else qr(z0$hessian)$rank
		if(a==np)nlcov <- solve(z0$hessian)
		else nlcov <- matrix(NA,ncol=np,nrow=np)
		nlse <- sqrt(diag(nlcov))}
	if(length(nlse)>1)nlcorr <- nlcov/(nlse%o%nlse)
	else nlcorr <- nlcov/nlse^2
	coef <- z0$estimate
	grad <- z0$gradient
	code <- z0$code
	iterations <- z0$iterations}
else nlcov <- nlse <- nlcorr <- coef <- grad <- code <- iterations <- NULL
#
# calculate random coefficient covariance matrix
#
if(nod>0){
	a <- matrix(0,nrow=max(position[,2]),ncol=max(position[,2]))
	ii <- sum(arma)
	for(i in 1:nod){
		ii <- ii+1
		a[position[i,1],position[i,2]] <- p[ii]}}
else a <- NULL
#
# calculate recursive fitted values
#
if(resid){
	coln <- NULL
	if(ncv>0)for(i in 1:ns)for(j in 1:nobs(response)[i])
		coln <- rbind(coln,ccov[1:ncv,i])
	pred <- rep(beta[1],nt)
	cc <- matrix(rep(1,nt),ncol=1)
	if(torder>1)for(i in 1:(torder-1)) {
		cc <- cbind(cc,times^i)
		pred <- pred+beta[i+1]*times^i}
	if(ncv>0){
		jj <- torder
		for(i in 1:ncv)for(j in 1:interaction[i]){
			jj <- jj+1
			cc <- cbind(cc,cc[,j]*coln[,i])
			pred <- pred+beta[jj]*cc[,j]*coln[,i]}}
	z <- .Fortran("resid",
		np=as.integer(length(p)),
		par=as.double(coef),
		beta=as.double(beta),
		ave=as.double(ave),
		pred=double(nt),
		sdr=double(nt),
		res=double(nt),
		y=as.double(y),
		sse=as.double(sse),
		nq=as.integer(nre),
		nlp=as.integer(nlp),
		ns=as.integer(ns),
		nt=as.integer(nt),
		model=as.integer(arma),
		t=as.double(times),
		nobs=as.integer(nobs(response)),
		nod=as.integer(nod),
		as.integer(position),
		cv=as.double(ccov),
		nxcv=as.integer(interaction),
		nx=as.integer(torder),
		ncv=as.integer(ncv),
		p=double(length(p)),
		x=double(nlp+1),
		DUP=FALSE,
		PACKAGE="growth")
	rpred <- z$pred
	sdr <- z$sdr
	rres <- z$res}
else pred <- res <- rpred <- sdr <- rres <- NULL
z <- list(
	call=call,
	response=response$response,
	transform=transform,
	ccov=response$ccov,
	torder=torder-1,
	interaction=interaction-1,
	np=np,
	nlp=nlp,
	nre=nre,
	nod=nod,
	arma=arma,
	maxlike=like,
	aic=aic,
	df=nt-nlp-np-1,
	coefficients=coef,
	nlse=nlse,
	nlcov=nlcov,
	nlcorr=nlcorr,
	are=a,
	position=position,
	beta=beta,
	betase=betase,
	vbeta=vbeta,
	corr=corr,
	sse=sse,
	mse=mse,
	pred=pred,
	rpred=rpred,
	sdrpred=sdr,
	rresiduals=rres,
	gradient=grad,
	iterations=iterations,
	code=code)
class(z) <- "carma"
if(resid)class(z) <- c(class(z),"recursive")
return(z)}

### standard methods
###
coef.carma <- function(z) list(beta=z$beta,coef=z$coefficients)

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

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

residuals.carma <- function(z, recursive=TRUE){
if(recursive)return(z$rresiduals)
else {
	if(z$transform=="exp")z$response$y <- exp(z$response$y)
	else if(z$transform=="square")z$response$y  <- z$response$y^2
	else if(z$transform=="sqrt")z$response$y <- sqrt(z$response$y)
	else if(z$transform=="log")z$response$y <- log(z$response$y)
	return((z$response$y-z$pred)/sqrt(z$mse))}}

### print method
###
print.carma <- function(z,digits=max(3,.Options$digits-3),correlation=TRUE){
if(!is.null(z$ccov$ccov))nccov <- dim(z$ccov$ccov)[2]
else nccov <- 0
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
if(!is.null(z$code)&&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")
if(!is.null(z$response$times))cat("Mean time             ",mean(z$response$times),"\n")
cat("Transformation        ",z$trans,"\n\n")
cat("-Log likelihood   ",z$maxlike,"\n")
cat("Degrees of freedom",z$df,"\n")
cat("AIC               ",z$aic,"\n")
if(!is.null(z$iterations))cat("Iterations        ",z$iterations,"\n\n")
cat("Estimates of linear parameters\n")
if(inherits(z$ccov$linear,"formula"))
	cat("Formula: ",deparse(z$ccov$linear[[2]]),"\n")
coef.table <- cbind(z$beta, z$betase)
cname <- "(Intercept)"
if(z$torder>0)cname <- c(cname,paste("t^",1:z$torder,sep=""))
if(nccov>0)for(i in 1:nccov){
	cname <- c(cname,colnames(z$ccov$ccov)[i])
	if(!is.na(z$interaction)&&z$interaction[i]>0){
		cname <- c(cname,paste(colnames(z$ccov$ccov)[i],".t^",1:z$interaction[i],sep=""))}}
dimnames(coef.table) <- list(cname, c("estimate", "se"))
print.default(coef.table, digits=digits, print.gap=2)
if(z$nlp>1&&correlation){
	cat("\nCorrelation matrix of linear parameters\n")
	dimnames(z$corr) <- list(seq(1,z$nlp),seq(1,z$nlp))
	print.default(z$corr, digits=digits)}
cat("\n(REML) Variance ",z$mse,"\n")
if(z$np>0){
	cat("\nEstimates of nonlinear parameters\n")
	n1 <- z$arma[1]+z$arma[2]
	if(n1>1) {
		z0 <- .Fortran("roots",
			as.integer(n1),
			as.double(z$coef[1:n1]),
			r=complex(n1),
			DUP=FALSE,
			PACKAGE="growth")
		tmp <- if(any(Im(z0$r)!=0))z0$r else Re(z0$r)
		title <- "Roots"}
	if(z$arma[1]>0){
		cat("Autoregression\n")
		n2 <- z$arma[1]
		if(n2==1&&z$arma[2]==0){
			tmp <- exp(-exp(z$coef[1]))
			title <- "AR"}
		coef.table <- cbind(z$coef[1:n2],z$nlse[1:n2],
			tmp[1:n2])
		dimnames(coef.table) <- list(seq(1,z$arma[1]),
			c("estimate","se",title))
		print.default(coef.table, digits=digits, print.gap=2)}
	if(z$arma[2]>0){
		cat("Moving average\n")
		n1 <- z$arma[1]+1
		n2 <- z$arma[1]+z$arma[2]
		coef.table <- cbind(z$coef[n1:n2],z$nlse[n1:n2],
			tmp[n1:n2])
		dimnames(coef.table) <- list(seq(1,z$arma[2]),
			c("estimate","se",title))
		print.default(coef.table, digits=digits, print.gap=2)}
	if(z$arma[3]>0){
		cat("Measurement error\n")
		n1 <- z$arma[1]+z$arma[2]+1
		coef.table <- cbind(z$coef[n1],z$nlse[n1],exp(z$coef[n1]))
		dimnames(coef.table) <- list("1",c("estimate","se",""))
		print.default(coef.table, digits=digits, print.gap=2)}}
if(z$nod>0){
	cat("Estimated unscaled factored between subject covariance matrix\n")
	n1 <- z$arma[1]+z$arma[2]+z$arma[3]+1
	n2 <- n1+z$nre-1
	coef.table <- cbind(z$position[1:z$nre,,drop=FALSE],
		z$are[z$position[1:z$nre,,drop=FALSE]],z$nlse[n1:n2])
	dimnames(coef.table) <- list(rep("",z$nre),
		c("i1","i2","estimate","se"))
	print.default(coef.table, digits=digits, print.gap=2)
	cat("Estimated between subject covariance matrix\n")
	tmps <- max(z$position[,2])
	wrk <- matrix(0,ncol=tmps,nrow=tmps)
	for(j in 1:tmps)for(i in 1:j){
		for(l in 1:i)wrk[i,j] <- wrk[i,j]+z$are[l,i]*
			z$are[l,j]*z$mse
		wrk[j,i] <- wrk[i,j]}
	dimnames(wrk) <- list(1:tmps,1:tmps)
	print.default(wrk, digits=digits)}
if(z$np>1&&correlation){
	cat("\nCorrelation matrix of nonlinear parameters\n")
	dimnames(z$nlcorr) <- list(seq(1,z$np),seq(1,z$np))
	print.default(z$nlcorr, digits=digits)}}

### special marginal profiles with se's
###
mprofile.carma <- function(z, times=NULL, ccov, plotse=TRUE){
#
# if there are time-constant covariates, calculate how many and check
# for what values plot is to be made
#
if(!is.null(z$ccov$ccov)){
	nccov <- dim(z$ccov$ccov)[2]
	if(missing(ccov))stop("Covariate values must be supplied")
	else if(nccov!=length(ccov))
		stop("Incorrect number of covariate values")}
else nccov <- 0
#
# obtain plotting times
#
z$ptimes <- if(is.null(times))seq(min(z$response$times),
	max(z$response$times),length.out=25) else times
npts <- length(z$ptimes)
#
# calculate the polynomial
#
z$pred <- rep(z$beta[1],npts)
cc <- matrix(rep(1,npts),ncol=1)
if(z$torder>0)for(i in 1:z$torder) {
	cc <- cbind(cc,(z$ptimes-mean(z$response$times))^i)
	z$pred <- z$pred+z$beta[i+1]*(z$ptimes-mean(z$response$times))^i}
if(nccov>0){
	jj <- z$torder+1
	for(i in 1:nccov)for(j in 1:(z$interaction[i]+1)){
		jj <- jj+1
		cc <- cbind(cc,cc[,j]*ccov[i])
		z$pred <- z$pred+z$beta[jj]*cc[,j]*ccov[i]}}
se <- NULL
for(i in 1:npts) se <- c(se,sqrt(cc[i,]%*%z$vbeta%*%cc[i,]))
#
# response was transformed, transform fitted values
#
if(z$transform=="exp"){
	if(plotse){
		se1 <- log(z$pred+2*se)
		se2 <- log(z$pred-2*se)}
	z$pred <- log(z$pred)}
else if(z$transform=="square"){
	if(plotse){
		se1 <- sqrt(z$pred+2*se)
		se2 <- sqrt(z$pred-2*se)}
	z$pred  <- sqrt(z$pred)}
else if(z$transform=="sqrt"){
	if(plotse){
		se1 <- (z$pred+2*se)^2
		se2 <- (z$pred-2*se)^2}
	z$pred <- z$pred^2}
else if(z$transform=="log"){
	if(plotse){
		se1 <- exp(z$pred+2*se)
		se2 <- exp(z$pred-2*se)}
	z$pred <- exp(z$pred)}
else {
	if(plotse){
		se1 <- z$pred+2*se
		se2 <- z$pred-2*se}}
#
# standard errors, if requested
#
if(plotse)z$pse <- cbind(se1,se2)
class(z) <- "mprofile"
invisible(z)}
#
#  growth : A Library of Normal Distribution Growth Curve Models
#  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
#
#     elliptic(response=NULL, model="linear", distribution="normal",
#	times=NULL, dose=NULL, ccov=NULL, tvcov=NULL, nest=NULL, torder=0,
#	interaction=NULL, transform="identity", link="identity",
#	autocorr="exponential", pell=NULL, preg=NULL, pvar=var(y), varfn=NULL,
#	par=NULL, pre=NULL, delta=NULL, shfn=FALSE, common=FALSE,
#	envir=parent.frame(), print.level=0, ndigit=10,
#	gradtol=0.00001, steptol=0.00001, iterlim=100, fscale=1,
#	stepmax=10*sqrt(theta%*%theta), typsiz=abs(c(theta)))
#
#  DESCRIPTION
#
#    Function to fit the multivariate elliptical distribution with
# various autocorrelation functions, one or two levels of random
# effects, and nonlinear regression.

elliptic <- function(response=NULL, model="linear", distribution="normal",
	times=NULL, dose=NULL, ccov=NULL, tvcov=NULL, nest=NULL,
	torder=0, interaction=NULL, transform="identity",
	link="identity", autocorr="exponential", pell=NULL, preg=NULL,
	pvar=var(y), varfn=NULL, par=NULL, pre=NULL, delta=NULL,
	shfn=FALSE, common=FALSE, envir=parent.frame(),
	print.level=0, ndigit=10, gradtol=0.00001, steptol=0.00001,
	iterlim=100, fscale=1, stepmax=10*sqrt(theta%*%theta),
	typsiz=abs(c(theta))){
#
# likelihood function returning parameter values
#
plra <- function(theta){
	if(mdl==2)mu <- mu1(theta)
	if(cvar==1)varn <- sh1(theta)
	z <- .Fortran("plra",
		theta=as.double(theta),
		like=double(1),
		dist=as.integer(dst),
		as.double(rxl),
		x=as.double(times),
		as.double(y),
		tvcov=as.double(resp$tvcov$tvcov),
		ccov=as.double(resp$ccov$ccov),
		dose=as.double(dose),
		nobs=as.integer(nobs),
		nbs=as.integer(n),
		nest=as.integer(resp$response$nest),
		lnest=as.integer(lnest),
		dev=double(nm),
		nind=as.integer(nind),
		nld=as.integer(nld),
		nxrl=as.integer(nxrl),
		np=as.integer(np),
		npell=as.integer(npell),
		npv=as.integer(npv),
		npvl=as.integer(npvl),
		nccov=as.integer(nccov),
		npvar=as.integer(npvar),
		cvar=as.integer(cvar),
		npre=as.integer(npre),
		npar=as.integer(npar),
		link=as.integer(lnk),
		torder=as.integer(torder),
		inter=as.integer(interaction),
		model=as.integer(mdl),
		ar=as.integer(ar),
		tvc=as.integer(ntvc),
		beta=double(npv2),
		betacov=double(npv2*npv2),
		v=double(nld*nld),
		sigsq=double(nld),
		ey=double(nld),
		tb=double(npvl),
		as.double(mu),
		as.double(varn),
		DUP=FALSE,
		PACKAGE="growth")
	list(like=z$like,res=z$dev,beta=z$beta,betacov=z$betacov)}
#
# likelihood function for nlm
#
plral <- function(theta){
	if(mdl==2)mu <- mu1(theta)
	if(cvar==1)varn <- sh1(theta)
	z <- .Fortran("plra",
		theta=as.double(theta),
		like=double(1),
		dist=as.integer(dst),
		as.double(rxl),
		x=as.double(times),
		as.double(y),
		tvcov=as.double(resp$tvcov$tvcov),
		ccov=as.double(resp$ccov$ccov),
		dose=as.double(dose),
		nobs=as.integer(nobs),
		nbs=as.integer(n),
		nest=as.integer(resp$response$nest),
		lnest=as.integer(lnest),
		dev=double(nm),
		nind=as.integer(nind),
		nld=as.integer(nld),
		nxrl=as.integer(nxrl),
		np=as.integer(np),
		npell=as.integer(npell),
		npv=as.integer(npv),
		npvl=as.integer(npvl),
		nccov=as.integer(nccov),
		npvar=as.integer(npvar),
		cvar=as.integer(cvar),
		npre=as.integer(npre),
		npar=as.integer(npar),
		link=as.integer(lnk),
		torder=as.integer(torder),
		inter=as.integer(interaction),
		model=as.integer(mdl),
		ar=as.integer(ar),
		tvc=as.integer(ntvc),
		beta=double(npv2),
		betacov=double(npv2*npv2),
		v=double(nld*nld),
		sigsq=double(nld),
		ey=double(nld),
		tb=double(npvl),
		as.double(mu),
		as.double(varn),
		DUP=FALSE,
		PACKAGE="growth")
	z$like}
call <- sys.call()
#
# check type of model to be fitted
#
if(!is.function(model)&&!inherits(model,"formula")&&!is.null(model))
	model <- match.arg(model,c("linear","logistic","pkpd"))
tmp <- c("exponential","gaussian","cauchy","spherical","IOU")
ar <- match(autocorr <- match.arg(autocorr,tmp),tmp)
tmp <- c("normal","power exponential","Student t","Laplace")
dst <- match(distribution <- match.arg(distribution,tmp),tmp)
if(dst==1)pell <- NULL
npell <- !is.null(pell)
if(!npell){
	if(dst==2)
		stop("An estimate of the elliptic power parameter must be supplied")
	else if(dst==3)
		stop("An estimate of the degrees of freedom must be supplied")}
if(npell&&pell<=0){
	if(dst==2)stop("The elliptic power parameter must be positive.")
	else if(dst==3)
		stop("The degrees of freedom parameter must be positive.")}
#
# check transformation and link
#
transform <- match.arg(transform,c("identity","exp","square","sqrt","log"))
tmp <- c("identity","exp","square","sqrt","log")
lnk <- match(link <- match.arg(link,tmp),tmp)
#
# check if mean and variance can have common parameters
#
if(common){
	if(!is.function(model)&&!inherits(model,"formula"))
		stop("with common parameters, model must be a function or formula")
	if(!is.function(varfn)&&!inherits(varfn,"formula"))
		stop("with common parameters, varfn must be a function or formula")
	pvar <- NULL}
npar <- length(par)
ntvc <- if(!is.null(tvcov))1 else 0
#
# check if a data object is being supplied
#
respenv <- exists(deparse(substitute(response)),env=parent.frame())&&
	inherits(response,"repeated")&&!inherits(envir,"repeated")
if(respenv){
	if(dim(response$response$y)[2]>1)
		stop("elliptic only handles univariate responses")
	if(!is.null(response$NAs)&&any(response$NAs))
		stop("elliptic does not handle data with NAs")}
envname <- if(respenv)deparse(substitute(response))
	else if(!is.null(class(envir)))deparse(substitute(envir))
	else NULL
#
# if envir, remove extra (multivariate) responses
#
type <- "unknown"
if(!respenv&&inherits(envir,"repeated")){
	if(!is.null(envir$NAs)&&any(envir$NAs))
		stop("elliptic 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 response is not a data object, make one, and prepare covariates
#
if(inherits(response,"repeated")){
	if(dim(response$response$y)[2]>1)
		stop("elliptic only handles univariate responses")
	if(!is.null(response$NAs)&&any(response$NAs))
		stop("elliptic does not handle data with NAs")
	if(is.character(model)){
		resp <- response
		if(is.null(ccov))resp$ccov <- NULL
		else if(inherits(ccov,"formula")){
		# find time-constant covariates and put in data object
			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")
		if(is.null(tvcov))resp$tvcov <- NULL
		else if(inherits(tvcov,"formula")){
		# find time-varying covariates and put in data object
#			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")}
	else resp <- rmna(response$response)
	type <- resp$response$type
	nccov <- if(is.null(resp$ccov$ccov)) 0 else  dim(resp$ccov$ccov)[2]
	ntvc <- if(is.null(resp$tvcov$tvcov)) 0 else  dim(resp$tvcov$tvcov)[2]}
else {
	resp <- if(inherits(response,"response")) response
		else restovec(response,times,nest=nest,delta=delta)
	type <- resp$type
	if(dim(resp$y)[2]>1)stop("elliptic only handles univariate responses")
	if(is.null(ccov))nccov <- 0
	else {
	# set up time-constant covariates
		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 <- dim(ccov$ccov)[2]}
	if(is.null(tvcov))ntvc <- 0
	else {
	# set up time-varying covariates
		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)}
		ntvc <- 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((inherits(envir,"repeated")&&(length(nobs(resp))!=length(nobs(envir))||
	any(nobs(resp)!=nobs(envir))))||(inherits(envir,"tvcov")&&
	(length(nobs(resp))!=length(envir$tvcov$nobs)||
	any(nobs(resp)!=envir$tvcov$nobs))))
	stop("response and envir objects are incompatible")
if(type!="unknown"&&type!="continuous"&&type!="duration")
	stop("continuous or duration data required")
if(!is.null(resp$response$censor))stop("elliptic does not handle censoring")
full <- !is.null(resp$ccov$ccov)
y <- resp$response$y
n <- length(y)
times <- resp$response$times
nobs <- nobs(resp)
nind <- length(nobs)
nld <- max(nobs)
mu <- varn <- NULL
npre <- length(pre)
#
# 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(inherits(model,"formula")){
		mu3 <- if(respenv)finterp(model,.envir=response,.name=envname)
			else finterp(model,.envir=envir,.name=envname)}
	else if(is.function(model)){
		if(is.null(attr(model,"model"))){
		        tmp <- parse(text=deparse(model)[-1])
		        model <- if(respenv)fnenvir(model,.envir=response,.name=envname)
		        	else fnenvir(model,.envir=envir,.name=envname)
		        mu3 <- model
		        attr(mu3,"model") <- tmp}
		else mu3 <- model}
	if(inherits(varfn,"formula")){
		sh3 <- if(respenv)finterp(varfn,.envir=response,.name=envname)
			else finterp(varfn,.envir=envir,.name=envname)}
	else if(is.function(varfn)){
		if(is.null(attr(varfn,"model"))){
		        tmp <- parse(text=deparse(varfn)[-1])
		        varfn <- if(respenv)fnenvir(varfn,.envir=response,.name=envname)
		        	else fnenvir(varfn,.envir=envir,.name=envname)
		        sh3 <- varfn
		        attr(sh3,"model") <- tmp}
		else sh3 <- varfn}}
npr <- length(preg)
#
# transform model formula to function and check number of parameters
#
if(inherits(model,"formula")){
	mu2 <- if(respenv)finterp(model,.envir=response,.name=envname)
		else finterp(model,.envir=envir,.name=envname)
	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&&!common){
			cat("\nParameters are ")
			cat(attr(mu2,"parameters"),"\n")
			stop(paste("preg should have",npt1,"estimates"))}
		if(is.list(preg)){
			if(!is.null(names(preg))){
				o <- match(attr(mu2,"parameters"),names(preg))
				preg <- unlist(preg)[o]
				if(sum(!is.na(o))!=length(preg))stop("invalid estimates for model - probably wrong names")}
			else preg <- unlist(preg)}}
	if(!is.null(mu2)){
		if(inherits(envir,"tccov")){
		# fix length if only time-constant covariates
			cv <- covind(response)
			mu1 <- function(p) mu2(p)[cv]
			attributes(mu1) <- attributes(mu2)}
		else {
			mu1 <- mu2
			rm(mu2)}}}
else if(is.function(model))mu1 <- model
else mu1 <- NULL
#
# give appropriate attributes to mu1 for printing
#
if(!is.null(mu1)&&is.null(attr(mu1,"parameters"))){
	attributes(mu1) <- if(is.function(model)){
		if(!inherits(model,"formulafn")){
			if(respenv)attributes(fnenvir(model,.envir=response))
			else attributes(fnenvir(model,.envir=envir))}
		else attributes(model)}
		else {
			if(respenv)attributes(fnenvir(mu1,.envir=response))
			else attributes(fnenvir(mu1,.envir=envir))}}
#
# if possible, 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)&&!common&&nlp!=npr)
	stop(paste("preg should have",nlp,"initial estimates"))
if(is.function(mu1))mdl <- 2
else if(model=="linear"){
	if(dst==4)stop("linear model option not allowed with Laplace distribution")
	mdl <- 1
	torder <- torder+1}
else if(model=="logistic")mdl <- 3
else if(model=="pkpd")mdl <- 4
if(mdl==2&&length(mu1(preg))!=sum(nobs(resp)))
	stop("The mean function must provide an estimate for each observation")
npvar <- length(pvar)
n1 <- if(common){if(inherits(varfn,"formula"))
	length(attr(mu1,"parameters"))+1 else 1} else length(preg)+1
n2 <- length(c(preg,pvar))
#
# transform variance formula to function and check number of parameters
#
if(inherits(varfn,"formula")){
	old <- if(common)mu1 else NULL
	mufn <- if(shfn)"mu" else NULL
	sh3 <- if(respenv)finterp(varfn,.envir=response,.start=n1,.name=envname,.old=old,.args=mufn)
		else finterp(varfn,.envir=envir,.start=n1,.name=envname,.old=old,.args=mufn)
	tmp <- attributes(sh3)
	sh2 <- if(shfn)function(p) sh3(p,mu1(p)) else sh3
	attributes(sh2) <- tmp
	npt2 <- length(attr(sh2,"parameters"))
	if(is.character(attr(sh2,"model"))){
	# W&R formula
		if(length(attr(sh2,"model"))==1){
			sh1 <- function(p) p[n1]*rep(1,n)
			attributes(sh1) <- attributes(sh2)
			sh2 <- NULL}}
	else {
	# formula with unknowns
		if(npvar!=npt2&&!common){
			cat("\nParameters are ")
			cat(attr(sh2,"parameters"),"\n")
			stop(paste("pvar should have",npt2,"estimates"))}
		if(is.list(pvar)){
			if(!is.null(names(pvar))){
				o <- match(attr(sh2,"parameters"),names(pvar))
				pvar <- unlist(pvar)[o]
				if(sum(!is.na(o))!=length(pvar))stop("invalid estimates for varfn - probably wrong names")}
			else pvar <- unlist(pvar)}}
	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(varfn))sh1 <- if(shfn) function(p) varfn(p[n1:n2],mu1(p))
		else function(p) varfn(p[n1:n2])
else sh1 <- NULL
#
# give appropriate attributes to sh1 for printing
#
if(!is.null(sh1)&&is.null(attr(sh1,"parameters")))
	attributes(sh1) <- if(is.function(varfn)){
		if(!inherits(varfn,"formulafn")){
			if(respenv)attributes(fnenvir(varfn,.envir=response))
			else attributes(fnenvir(varfn,.envir=envir))}
		else attributes(varfn)}
		else {
			if(respenv)attributes(fnenvir(sh1,.envir=response))
			else attributes(fnenvir(sh1,.envir=envir))}
nlp <- if(is.function(varfn))length(attr(sh1,"parameters"))-shfn
	else if(is.null(varfn)||is.character(varfn))NULL
	else if(!inherits(varfn,"formula"))stop("varfn must be a function or formula")
	else npt2
#
# if possible, check that correct number of estimates was supplied
#
if(!is.null(nlp)&&!common&&nlp!=npvar)
	stop(paste("pvar should have",nlp,"initial estimates"))
if(common){
	nlp <- length(unique(c(attr(mu1,"parameters"),attr(sh1,"parameters"))))
	if(nlp!=npr)stop(paste("with a common parameter model, preg should contain",nlp,"estimates"))}
#
# check interactions of covariates with times
#
if(mdl==1&&!is.null(interaction)){
	if(length(interaction)!=nccov)
		stop(paste(nccov,"interactions with time must be specified"))
	else if(any(interaction>torder-1))
		stop(paste("Interactions can be at most of order ",torder-1))
	else interaction <- interaction+1}
else interaction <- rep(1,nccov)
if(mdl==4&&ntvc==0&&is.null(dose))stop("Doses required for PKPD model")
#
# set up nesting indicator
#
if(!is.null(resp$response$nest))lnest <- max(resp$response$nest)
else {
	lnest <- 0
	resp$response$nest <- rep(1,length(y))}
#
# check times
#
if(!is.null(times)){
	if(mdl==1){
		if(torder>length(unique(times)))
			stop("torder is too large for the number of distinct times")
		ave <- mean(times)
		times <- times-ave}
	else ave <- 0}
else {
	if(!is.null(par))stop("No times. AR cannot be fitted")
	if(torder>1)stop("No times. Time trends cannot be fitted.")
	ave <- 0}
#
# set up for full logistic model
#
if(full){
	rxl <- NULL
	for(i in 1:nind){
		tmp <- 1
		for(j in 1:dim(resp$ccov$ccov)[2])
			tmp <- tmp+resp$ccov$ccov[i,j]*2^(j-1)
		rxl <- c(rxl,tmp)}
	nxrl <- max(rxl)
	p <- preg}
else {
	nxrl <- 1
	rxl <- matrix(1,ncol=1,nrow=nind)
	p <- NULL}
nm <- sum(nobs(resp))
#
# calculate transformation and Jacobian
#
if(transform=="identity")jacob <- 0
else if(transform=="exp"){
	jacob <- -sum(y)
	y <- exp(y)}
else if(any(y==0))stop("Zero response values: invalid transformation")
else if(transform=="square"){
	jacob <- -sum(log(asb(y)))-length(resp$response$y)*log(2)
	y  <- y^2}
else if(any(y<0))stop("Nonpositive response values: invalid transformation")
else if(transform=="sqrt"){
	jacob <- sum(log(y))/2+length(resp$response$y)*log(2)
	y <- sqrt(y)}
else if(transform=="log"){
	jacob <- sum(log(y))
	y <- log(y)}
#
# include delta in Jacobian
#
if(!is.null(resp$response$delta)){
	if(length(resp$response$delta)==1)jacob <- jacob-nm*log(resp$response$delta)
	else jacob <- jacob -sum(log(resp$response$delta))}
#
# special cases for various models
#
if(mdl==1){
# linear model
	if(lnk==1&&is.null(varfn)){
		npvl <- torder+sum(interaction)+ntvc
		npv <- 0}
	else {
		npv <- torder+sum(interaction)+ntvc
		npvl <- 0
		if(is.null(preg)){
			if(lnk!=1)stop(paste("Initial regression estimates must be supplied for the linear model with ",link,"link"))
			else stop("Initial regression estimates must be supplied for the linear model when there is a variance function")}}}
else {
	npvl <- 0
	if(mdl==2)npv <- length(preg)
	else if(mdl==3)npv <- if(ntvc==0)4*nxrl else 4+nxrl-1
	else if(mdl==4){
		if(!is.function(varfn))npv <- 2+nxrl
		else npv <- 3}}
npv2 <- max(npv,npvl)
if(length(preg)!=npv)
	stop(paste(npv,"initial parameter estimates for the mean regression must be supplied"))
#
# set up variance function if there is one and check initial estimates
#
if(is.null(varfn))cvar <- 0
else if(is.function(sh1)){
	cvar <- 1
	if(length(sh1(if(common)preg else pvar))!=sum(nobs(resp)))
		stop("The variance function or formula must provide an estimate for each observation")}
else if(varfn=="identity"){
	if(any(y<0))warning("identity variance function not recommended with negative responses")
	cvar <- 2}
else if(varfn=="square")cvar <- 3
else stop("Unknown variance function: choices are identity and square")
if(!(mdl==4&&npvar==4)&&cvar==0){
	if(any(pvar<=0))stop("All variance parameters must be positive")
	pvar <- log(pvar)}
theta <- c(preg,pvar)
#
# check AR and random effect
#
if(npre>0){
	if(any(pre<=0))stop("All variance components must be positive")
	theta <- c(theta,log(pre))}
if(npar>0){
	if(par<=0||(par>=1&&ar!=5))
		stop("Estimate of autocorrelation must lie between 0 and 1")
	theta <- if(ar!=5)c(theta,log(par/(1-par)))
		else c(theta,log(par))}
if(npell)theta <- c(theta,if(dst==4)pell else log(pell))
np <- npv+npvl*(mdl==1)+npvar+npre+npar+npell
#
# estimate model
#
if(fscale==1)fscale <- plral(theta)
z0 <- nlm(plral, p=theta, hessian=TRUE, print.level=print.level,
	typsiz=typsiz, ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
	steptol=steptol, iterlim=iterlim, fscale=fscale)
p <- z0$estimate
z <- plra(p)
like <- z$like+nm*log(pi)/2+jacob
pred <- y-z$res
#
# get parameters
#
sigsq <- p[(npv+1):(npv+npvar)]
if(!(mdl==4&&npvar==4)&&!common)sigsq <- exp(sigsq)
if(npre>0)tausq <- exp(p[(npv+npvar+1):(npv+npvar+npre)])
else tausq <- 0
if(npar>0){
	rho <- exp(p[npv+npvar+npre+1])
	if(ar!=5)rho <- rho/(1+rho)}
else rho <- 0
#
# calculate se's
#
if(np-npvl==1){
	nlcov <- 1/z0$hessian
	nlse <- sqrt(nlcov)}
else {
	a <- if(any(is.na(z0$hessian))||any(abs(z0$hessian)==Inf))0
		else qr(z0$hessian)$rank
	if(a==np-npvl)nlcov <- solve(z0$hessian)
	else nlcov <- matrix(NA,ncol=np-npvl,nrow=np-npvl)
	nlse <- sqrt(diag(nlcov))}
if(length(nlse)>1)nlcorr <- nlcov/(nlse%o%nlse)
else nlcorr <- as.matrix(nlcov/nlse^2)
dimnames(nlcorr) <- list(seq(1,np-npvl),seq(1,np-npvl))
betase <- betacorr <- NULL
if(mdl==1){
	if(lnk==1&&cvar==0){
		betacov <- matrix(z$betacov,ncol=npvl)
		bt <- exp(p[length(p)])
		if(npell>0)betacov <- betacov/(npvl*gamma(npvl/(2*bt)))*
			2^(1/bt)*gamma((npvl+2)/(2*bt))
		if(npvl==1)betase <- sqrt(betacov)
		else if(npvl>1)betase <- sqrt(diag(betacov))
		if(npvl>1){
			betacorr <- betacov/(betase%o%betase)
			dimnames(betacorr) <- list(seq(1,npvl),seq(1,npvl))}}
	else {
		betase <- nlse[1:npv]
		betacov <- nlcov[1:npv,1:npv]
		if(npvl>1){
			betacorr <- nlcorr[1:npv,1:npv]
			dimnames(betacorr) <- list(seq(1,npv),seq(1,npv))}}}
else betacov <- NULL
#
# if possible, calculate recursive fitted values
#
if(dst==1&&!common&&length(sigsq)==1&&((npar>0&&
	autocorr=="exponential")||(npar==0&&length(pre)==1))){
	nt <- sum(nobs(resp))
	mse <- sigsq*length(y)/(length(y)-npv-npvl*(mdl==1))
	coef <- NULL
	if(npar>0)coef <- c(coef,log(-log(rho)))
	if(length(pre)==1)coef <-
	c(coef,sqrt(exp(p[npv+npvar+1])/mse))
	z1 <- .Fortran("resid2",
		np=as.integer(npar+length(pre)),
		par=as.double(coef),
		ave=as.double(ave),
		pred=as.double(pred),
		rpred=double(nt),
		sdr=double(nt),
		res=double(nt),
		y=as.double(y),
		mse=as.double(mse),
		ns=as.integer(nind),
		nt=as.integer(nt),
		model=as.integer(npar),
		t=as.double(times),
		nobs=as.integer(nobs(resp)),
		nod=as.integer(length(pre)),
		p=double(length(p)),
		DUP=FALSE,
		PACKAGE="growth")
	rpred <- z1$rpred
	ii <- covind(resp$response)
	kk <- unique(resp$response$nest)
	if(lnest)for(i in 1:nind)for(k in kk)
		z1$rpred[i==ii&k==resp$response$nest][1] <- z1$pred[i==ii&k==resp$response$nest][1]
	sdr <- z1$sdr
	rres <- z1$res}
else rpred <- sdr <- rres <- NULL
if(!is.null(mu3))mu1 <- mu3
if(!is.null(sh3))sh1 <- sh3
z <- list(
	call=call,
	model=model,
	distribution=distribution,
	mu1=mu1,
	varfn=varfn,
	common=common,
	sh1=sh1,
	shfn=shfn,
	autocorr=autocorr,
	response=resp$response,
	transform=transform,
	torder=torder-1,
	interaction=interaction-1,
	ccov=resp$ccov,
	tvcov=resp$tvcov,
	full=full,
	link=link,
	maxlike=like,
	aic=like+np,
	df=nm-np,
	np=np,
	npell=npell,
	npv=npv,
	npvl=npvl,
	npvar=npvar,
	npar=npar,
	npre=npre,
	coefficients=p,
	beta=z$beta,
	betacov=betacov,
	betacorr=betacorr,
	betase=betase,
	nlse=nlse,
	nlcov=nlcov,
	nlcorr=nlcorr,
	sigsq=sigsq,
	tausq=tausq,
	rho=rho,
	residuals=z$res,
	pred=pred,
	rpred=rpred,
	sdrpred=sdr,
	rresiduals=rres,
	grad=z0$gradient,
	iterations=z0$iterations,
	code=z0$code)
class(z) <- "elliptic"
if(!is.null(rpred))class(z) <- c(class(z),"recursive")
return(z)}

### standard methods
###

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

fitted.elliptic <- function(z, recursive=FALSE){
if(recursive){
	if(is.null(z$rpred))stop("recursive fitted values not available")
	z$rpred}
else z$pred}

residuals.elliptic <- function(z, recursive=FALSE){
if(recursive){
	if(is.null(z$rresiduals))stop("recursive residuals not available")
	return(z$rresiduals)}
else {
	if(z$transform=="exp")z$response$y <- exp(z$response$y)
	else if(z$transform=="square")z$response$y  <- z$response$y^2
	else if(z$transform=="sqrt")z$response$y <- sqrt(z$response$y)
	else if(z$transform=="log")z$response$y <- log(z$response$y)
	return(z$response$y-z$pred)}}

### print method
###
print.elliptic <- function(z,digits=max(3,.Options$digits-3),correlation=TRUE){
if(!is.null(z$ccov$ccov))nccov <- dim(z$ccov$ccov)[2]
else nccov <- 0
cat("\nMultivariate",z$distribution,"distribution\n")
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")
if(!inherits(z$mu1,"formulafn")){
	if(z$model=="linear"){
		if(z$torder>0){
			cat("\nPolynomial model\n")
			cat("Times centred at  ",mean(z$response$times),"\n\n")}
		else cat("\nLinear model\n\n")}
	else if(z$model=="logistic")cat("\nGeneralized logistic model\n\n")
	else if(z$model=="pkpd")cat("\nPKPD model\n\n")}
cat("Transformation:",z$trans,"\n")
cat("Link function: ",z$link,"\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")
ntvc <- !is.null(z$tvcov)
if(!inherits(z$mu1,"formulafn")){
	cat("Location parameters\n")
	if(inherits(z$ccov$linear,"formula"))
		cat("Linear part: ",deparse(z$ccov$linear),sep="\n")}
if(inherits(z$mu1,"formulafn")){
	if(z$common)cat("Location function\n")
	else cat("Location function parameters\n")
	if(!is.null(attr(z$mu1,"formula")))
		cat(deparse(attr(z$mu1,"formula")),sep="\n")
	else if(!is.null(attr(z$mu1,"model"))){
		t <- deparse(attr(z$mu1,"model"))
		t[1] <- sub("expression\\(","",t[1])
		t[length(t)] <- sub("\\)$","",t[length(t)])
		cat(t,sep="\n")}
	cname <- if(is.character(attr(z$mu1,"model")))attr(z$mu1,"model")
		else attr(z$mu1,"parameters")
	coef.table <- cbind(z$coef[1:z$npv],z$nlse[1:z$npv])
	if(!z$common){
		dimnames(coef.table) <- list(cname,c("estimate","se"))
		print.default(coef.table,digits=digits,print.gap=2)}}
else if(z$model=="linear"){
	if(z$npell==1&&z$link=="identity"&&!is.function(z$varfn))
		cat("(Approximate s.e.)\n")
	tord <- z$torder+1+nccov
	if(ntvc)tord <- tord+1
	coef.table <- cbind(z$beta,z$betase)
	cname <- "(Intercept)"
	if(z$torder>0)cname <- c(cname,paste("t^",1:z$torder,sep=""))
	if(nccov>0)for(i in 1:nccov){
		cname <- c(cname,colnames(z$ccov$ccov)[i])
		if(z$interaction[i]>0){
			cname <- c(cname,paste(colnames(z$ccov$ccov)[i],".t^",1:z$interaction[i],sep=""))}}
	if(ntvc)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(z$npvl>1&&z$link=="identity"&&is.null(z$varfn)&&correlation){
		cat("\nCorrelation matrix of linear parameters\n")
		print.default(z$betacorr,digits=digits)}}
else if(z$model=="logistic"){
	coef.table <- cbind(z$coef[1:4],z$nlse[1:4])
	if(ntvc)cname  <- c("kappa1","kappa3","kappa4","beta")
	else cname <- c("kappa1","kappa2","kappa3","kappa4")
	dimnames(coef.table) <- list(cname,c("estimate","se"))
	print.default(coef.table,digits=digits,print.gap=2)
	if(z$full){
		if(ntvc){
			coef.table <- cbind(z$coef[5:z$npv],z$nlse[5:z$npv])
			cname <- colnames(z$ccov$ccov)
			dimnames(coef.table) <- list(cname,c("estimate","se"))
			print.default(coef.table,digits=digits,print.gap=2)}
		else {
			for(i in 1:nccov){
				cat("   ",colnames(z$ccov$ccov)[i],"\n")
				coef.table <- cbind(z$coef[(i*4+1):((i+1)*4)],z$nlse[(i*4+1):((i+1)*4)])
				dimnames(coef.table) <- list(cname,c("estimate", "se"))
				print.default(coef.table,digits=digits,print.gap=2)}}}}
else if(z$model=="pkpd"){
	coef.table <- cbind(z$coef[1:3], z$nlse[1:3])
	cname <- c("log k_a","log k_e","log V")
	dimnames(coef.table) <- list(cname,c("estimate","se"))
	print.default(coef.table,digits=digits,print.gap=2)
	if(z$full){
		for(i in 1:nccov){
			cat("   ",colnames(z$ccov$ccov)[i],"\n")
			coef.table <- cbind(z$coef[i+3],z$nlse[i+3])
			dimnames(coef.table) <- list(cname[3],c("estimate","se"))
			print.default(coef.table,digits=digits,print.gap=2)}}}
if(inherits(z$sh1,"formulafn")){
	if(z$common){
		if(z$distribution=="normal")cat("\nVariance function\n")
		else cat("\nDispersion function\n")}
	else {
		if(z$distribution=="normal")
			cat("\nVariance function parameters\n")
		else cat("\nDispersion function parameters\n")
		cname <- NULL}
	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 <- c(cname,if(is.character(attr(z$sh1,"model")))
		attr(z$sh1,"model")
		else if(is.expression(attr(z$sh1,"model")))
		attr(z$sh1,"parameters")
		else attr(z$sh1,"parameters")[1:z$npvar])
	if(z$shfn)cname <- cname[cname!="mu"]
	if(!z$common)coef.table <- cbind(z$coef[(z$npv+1):(z$npv+z$npvar)], z$nlse[(z$npv+1):(z$npv+z$npvar)])
	else cat("\nCommon parameters\n")
	dimnames(coef.table) <- list(unique(cname),c("estimate","se"))
	print.default(coef.table,digits=digits,print.gap=2)}
else if(!is.function(z$model)&&!inherits(z$mu1,"formulafn")&&z$model=="pkpd"
	&&z$npvar==4){
	if(z$distribution=="normal")cat("\nVariance parameters\n")
	else cat("\nDispersion parameters\n")
	coef.table <- cbind(z$sigsq, z$nlse[(z$npv+1):(z$npv+4)])
	cname <- c("log k_a","log k_e","log V","power")
	dimnames(coef.table) <- list(cname,c("estimate","se"))
	print.default(coef.table,digits=digits,print.gap=2)}
else {
	if(z$distribution=="normal"){
		if(z$npvar==1)cat("\nVariance\n")
		else cat("\nVariance parameters\n")}
	else {
		if(z$npvar==1)cat("\nDispersion\n")
		else cat("\nDispersion parameters\n")}
	if(!is.null(z$varfn)){
		cat(z$varfn,"function of the location parameter\n")
		vname <- "factor"
		if(z$npvar==2)vname <- c("constant",vname)}
	else if(z$npvar>1)
		vname <- c("(Intercept)",paste("t^",1:(z$npvar-1),sep=""))
	else vname <- ""
	coef.table <- cbind(z$coef[(z$npv+1):(z$npv+z$npvar)],
		z$nlse[(z$npv+1):(z$npv+z$npvar)],z$sigsq)
	dimnames(coef.table) <- list(vname,c("estimate","se","sigsq"))
	print.default(coef.table,digits=digits,print.gap=2)}
if(z$npre>0){
	cat("\nVariance components\n")
	coef.table <- cbind(z$coef[(z$npv+z$npvar+1):(z$npv+z$npvar+z$npre)],
		z$nlse[(z$npv+z$npvar+1):(z$npv+z$npvar+z$npre)],z$tausq)
	if(z$npre==1)cname <- "tausq"
	else cname <- c("Level 1","Level 2")
	dimnames(coef.table) <- list(cname,c("estimate","se",""))
	print.default(coef.table,digits=digits,print.gap=2)}
if(z$rho!=0){
	cat("\n",z$autocorr," autocorrelation\n",sep="")
	coef.table <- cbind(z$coef[z$npv+z$npvar+z$npre+1],
		z$nlse[z$npv+z$npvar+z$npre+1],z$rho)
	dimnames(coef.table) <- list("rho",c("estimate","se",""))
	print.default(coef.table,digits=digits,print.gap=2)}
if(z$npell==1){
	coef.table <- cbind(z$coef[z$np-z$npvl],z$nlse[z$np-z$npvl],
		if(z$distribution=="Laplace")NULL
			else exp(z$coef[z$np-z$npvl]))
	if(z$distribution=="power exponential"){
		cat("\nPower exponential parameter\n")
		dimnames(coef.table) <- list("",c("estimate","se","power"))}
	else if(z$distribution=="Student t"){
		cat("\nDegrees of freedom parameter\n")
		dimnames(coef.table) <- list("",c("estimate","se","d.f."))}
	else if(z$distribution=="Laplace"){
		cat("\nAsymmetry parameter\n")
		dimnames(coef.table) <- list("",c("estimate","se"))}
	print.default(coef.table,digits=digits,print.gap=2)}
if(z$np-z$npvl>1&&correlation){
	cat("\nCorrelation matrix of nonlinear parameters\n")
	print.default(z$nlcorr,digits=digits)}}
#
#  growth : A Library of Normal Distribution Growth Curve Models
#  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
#
#     potthoff(response, x=NULL, ccov=NULL, times=NULL, torder=0,
#	orthogonal=TRUE)
#
#  DESCRIPTION
#
#    Function to fit the Potthoff and Roy (1964) growth curve model

potthoff <- function(response, x=NULL, ccov=NULL, times=NULL, torder=0,
	orthogonal=TRUE){
pcov <- function(z, s, x){
	p <- t(x)%x%z
	solve(p%*%(diag(n)%x%s)%*%t(p))}
call <- sys.call()
#
# check response
#
if(is.data.frame(response))response <- as.matrix(response)
if(!is.matrix(response))stop("response must be a matrix or dataframe")
n <- dim(response)[1]
r <- dim(response)[2]
#
# set up complete set of intersubject differences
#
if(inherits(x,"formula")){
	formula <- x
	mt <- terms(x)
	if(is.numeric(mt[[2]]))x <- matrix(1,nrow=n)
	else {
		mf <- model.frame(mt,parent.frame(),na.action=na.fail)
		x <- model.matrix(mt, mf)}}
if(is.null(x)||dim(x)[1]!=n)stop(paste("x must be a matrix with",n,"rows"))
#
# set up covariates to be fitted
#
if(inherits(ccov,"formula")){
	formula <- ccov
	mt <- terms(ccov)
	if(is.numeric(mt[[2]]))ccov <- matrix(1,nrow=n)
	else {
		mf <- model.frame(mt,parent.frame(),na.action=na.fail)
		ccov <- model.matrix(mt, mf)}}
if(is.null(ccov)) xx <- matrix(rep(1,n),ncol=1)
else xx <- x
#
# set up times and orthogonalize if necessary
#
if(is.null(times))times <- 1:r
if(orthogonal)z <- rbind(rep(1,length(times)),t(contr.poly(times)))
else {
	z <- rep(1,r)
	tt <- times-sum(times)/r
	z <- rbind(z,tt)
	for(i in 2:(r-1))z <- rbind(z,tt^i)}
#
# calculate covariance matrices
#
s <- t(response)%*%(diag(n)-xx%*%solve(t(xx)%*%xx)%*%t(xx))%*%response/n
ss <- solve(t(response)%*%(diag(n)-x%*%solve(t(x)%*%x)%*%t(x))%*%response/n)
b <- solve(t(xx)%*%xx)%*%t(xx)%*%response
if(!is.matrix(b))b <- matrix(b,nrow=1)
#
# model with polynomial in time
#
if(is.numeric(torder)){
	zz <- z[1:(torder+1),]
	if(!is.matrix(zz))zz <- matrix(zz,nrow=1)
	b1 <- b%*%ss%*%t(zz)%*%solve(zz%*%ss%*%t(zz))
	s1 <- response-xx%*%b1%*%zz
	s1 <- t(s1)%*%s1/n}
else {
	zz <- diag(r)
	b1 <- b
	s1 <- s}
res <- response-xx%*%b1%*%zz
if(!is.matrix(b1))b1 <- matrix(b1,nrow=1)
if(is.matrix(b1)&&dim(b1)[1]>1){
	if(is.null(colnames(ccov))){
		tn <- "(Intercept)"
		if(dim(ccov)[2]>1)
			tn <- c(tn,paste("ccov",1:(dim(ccov)[2]-1),sep=""))
		colnames(ccov) <- tn}
	rownames(b1) <- colnames(ccov)}
else {
	b1 <- matrix(b1,nrow=1)
	rownames(b1) <- "Mean"}
like <- n*(r*(1+log(2*pi))+log(abs(det(s1))))/2
aic <- like+length(b1)+r*(r+1)/2
#
# calculate se's
#
pc <- pcov(zz,solve(s1),xx)
if(is.matrix(pc)){
	d <- sqrt(diag(pc))
	se <- matrix(d,ncol=dim(b1)[2],byrow=TRUE)}
else se <- sqrt(pc)
corr <- pc
if(is.matrix(pc))corr <- pc/(d%o%d)
if(is.matrix(se)&&dim(se)[1]>1)rownames(se) <- colnames(ccov)
else {
	se <- matrix(se,nrow=1)
	rownames(se) <- "Mean"}
if(!is.numeric(torder)){
	if(is.null(colnames(response)))
		colnames(b1) <- paste("t",1:dim(b1)[2],sep="")
	else colnames(b1) <- colnames(response)
	if(is.null(colnames(response)))
		colnames(se) <- paste("t",1:dim(se)[2],sep="")
	else colnames(se) <- colnames(response)}
d <- sqrt(diag(s1))
c1 <- s1
for(i in 2:r)for(j in 1:(i-1))c1[i,j] <- s1[i,j]/d[i]/d[j]
z <- list(
	call=call,
	y=response,
	x=x,
	time=z,
	torder=torder,
	ns=n,
	nt=n*r,
	df=n*r-(length(b1)+r*(r+1)/2),
	beta=b1,
	ccov=s1,
	maxlike=like,
	aic=aic,
	pcov=pc,
	pcorr=corr,
	se=se,
	corr=c1,
	residuals=res)
class(z) <- "potthoff"
return(z)}

### standard methods
###
coef.potthoff <- function(z) z$beta

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

residuals.potthoff <- function(z) z$residuals

### print method
###
print.potthoff <- function(z, digits = max(3, .Options$digits - 3)) {
cat("\nCall:",deparse(z$call),sep="\n")
cat("\n")
cat("Number of subjects    ",z$ns,"\n")
cat("Number of observations",z$nt,"\n")
cat("-Log likelihood   ",z$maxlike,"\n")
cat("Degrees of freedom",z$df,"\n")
cat("AIC               ",z$aic,"\n\n")
if(is.null(colnames(z$beta))&&z$torder<=4){
	tn <- "(Intercept)"
	if(z$torder>0)tn <- c(tn,paste("t^",1:z$torder,sep=""))
	colnames(z$beta) <- tn
	colnames(z$se) <- tn}
cat("Estimates of linear parameters\n")
print(z$beta)
cat("Standard errors\n")
print(z$se)
nlp <- length(z$beta)
if(nlp>1){
	cat("\nCorrelation matrix of linear parameters\n")
	dimnames(z$pcorr) <- list(seq(1,nlp),seq(1,nlp))
	print.default(z$pcorr,digits=digits)}
cat("\nEstimated covariance/correlation matrix\n")
if(is.null(colnames(z$corr))){
	tn <- paste("t^",1:dim(z$corr)[2],sep="")
	dimnames(z$corr) <- list(tn,tn)}
print(z$corr)}
#
#  growth : A Library of Normal Distribution Growth Curve Models
#  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
#
#     rmaov(response, tvcov=NULL, ccov=NULL, analysis=TRUE)
#
#  DESCRIPTION
#
#    Function to calculate the ANOVA table for a balanced split-plot design
#    Copyright (C) 1999 Ralf Goertz (ralf.goertz@uni-jena.de)

rmaov <- function(response, tvcov=NULL, ccov=NULL, analysis=TRUE){
	if(!is.matrix(response)||!is.numeric(response))stop("response must be a matrix")
	if(is.null(tvcov)){
		nfac <- 1
		nrep <- nlev <- dim(response)[2]
		rfacname <- "repfac"}
	else if(is.vector(tvcov,mode="numeric")||is.factor(tvcov)){
		nfac <- 1
		nrep <- nlev <- length(unique(tvcov))
		rfacname <- paste(deparse(substitute(tvcov)))
		tvcov <- as.data.frame(tvcov)}
	else if(is.data.frame(tvcov)||is.matrix(tvcov)){
		nfac <- dim(tvcov)[2]
		if(!is.null(colnames(tvcov)))rfacname <- colnames(tvcov)
		else rfacname <- paste("repfac",1:dim(tvcov)[2],sep="")
		nlev <- NULL
		nrep <- 1
		for(i in 1:dim(tvcov)[2]){
			nlev <- c(nlev,length(unique(tvcov[,i])))
			nrep <- nrep*nlev[i]}}
	else stop("If specified, tvcov must be a vector, matrix, or dataframe containing the repeated factors")
	if(is.vector(ccov,mode="numeric")||is.factor(ccov)){
		tmp <- paste(deparse(substitute(ccov)))
		ccov <- list(ccov)
		names(ccov) <- tmp}
	else if(is.matrix(ccov))ccov <- data.frame(ccov)
	else if(is.list(ccov)&&!is.data.frame(ccov))for(i in 1:length(ccov))
		if(is.matrix(ccov[[i]]))ccov[[i]] <- as.vector(t(ccov[[i]]))
	else if(!is.null(ccov))print("ccov must be a vector, matrix, dataframe, or list of time-constant covariates")
	ncases <- dim(response)[2]*dim(response)[1]
	if(nrep!=dim(response)[2])
		stop("Number of columns in response must match product of levels")
	block <- gl(dim(response)[1],nrep)
	nc <- n <- dim(response)[2]
	response <- as.vector(t(response))
	for(i in 1:nfac){
		n <- n/nlev[i]
		com <- paste("fm",as.character(i),"<-gl(nlev[i],n,ncases",sep="")
		if(!is.null(levels(tvcov[,i])))
			com <- paste(com,",labels=levels(tvcov[,i])",sep="")
		com <- paste(com,")\n",sep="")
		eval(parse(text=com))}
	if(!is.null(ccov)){
		tccov <- ccov
		if(length(tccov[[1]])!=ncases)for(i in 1:length(tccov)){
			if(!is.factor(tccov[[i]]))
				tccov[[i]] <- as.factor(tccov[[i]])
			ccov[[i]] <- gl(nlevels(tccov[[i]]),1,ncases,labels=levels(tccov[[i]]))
			ccov[[i]][(rep(1:length(tccov[[i]]),rep(nc,length(tccov[[i]])))-1)*nc+rep(1:nc,length(tccov[[i]]))] <- tccov[[i]][rep(1:length(tccov[[i]]),rep(nc,length(tccov[[i]])))]}}
	com <- paste("b<-data.frame(response=response,block=block")
	for(i in 1:nfac)com <- paste(com,",",rfacname[i],"=fm",as.character(i),sep="")
	if(!is.null(ccov)){
		if(is.null(names(ccov)))
			names(ccov) <- paste("ccov",1:length(ccov),sep="")
		for(i in 1:length(ccov))
			com <- paste(com,",",names(ccov)[i],"=ccov$",names(ccov)[i],sep="")}
	com <- paste(com,")\n",sep="")
	eval(parse(text=com))
	if(!analysis)return(b)
	efacterm <- facterm <- paste(rfacname,collapse="*")
	if(!is.null(ccov))
		for(i in 1:length(ccov))
			facterm <- paste(facterm,names(ccov)[i],sep="*")
	com <- paste("res<-aov(response~",facterm,"+Error(block/(",efacterm,")),data=b)\n",sep="")
	eval(parse(text=com))
	res}


