.packageName <- "ordinal"
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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
#
#  kalord(response,times=NULL,distribution="multinomial",
#         depend="independence",mu=NULL,ccov=NULL,tvcov=NULL,torder=0,
#         interaction=NULL,preg=NULL,ptvc=NULL,pinitial=1,pdepend=NULL,
#         envir=sys.frame(sys.parent()),optimize=T,print.level=0,
#         ndigit=10,gradtol=0.00001,steptol=0.00001,fscale=1,
#         iterlim=100,typsiz=abs(p),stepmax=10*sqrt(p%*%p))
#
#  DESCRIPTION
#
#    A function to fit various parameterization of the logistic
#  distribution inserted into a Pareto distribution with serial
#  dependence or gamma frailties using Kalman-type update for
#  ordinal longitudinal data.
#

kalord <- function(response,times=NULL,distribution="multinomial",depend="independence",mu=NULL,ccov=NULL,tvcov=NULL,torder=0,interaction=NULL,preg=NULL,
                   ptvc=NULL,pinitial=1,pdepend=NULL,envir=sys.frame(sys.parent()),optimize=T,print.level=0,ndigit=10,gradtol=0.00001,steptol=0.00001,
                   fscale=1,iterlim=100,typsiz=abs(p),stepmax=10*sqrt(p%*%p)) {
  series <- function(p) {
    if(rf)
      b <- mu1(p[-(1:(nc-1))])
    z <- .C("kord",
            p=as.double(p),
            y=as.integer(y),
            t=as.double(times),
            x=as.double(resp$ccov$ccov),
            nind=as.integer(nind),
            nobs=as.integer(nobs),
            nbs=as.integer(n),
            nc=as.integer(nc),
            nccov=as.integer(nccov),
            npv=as.integer(npv),
            init=as.integer(init),
            model=as.integer(mdl),
            dep=as.integer(dep),
            torder=as.integer(torder),
            inter=as.integer(interaction),
            tvc=as.integer(tvc),
            tvcov=as.double(resp$tvcov$tvcov),
            fit=as.integer(0),
            pred=double(n),
            rpred=double(n),
            cprob=double(n*(nc-1)),
            ppred=double(n),
            cpred=double(n*(nc-1)),
            creg=double(n*(nc-1)),
            rf=as.integer(rf),
            bbb=as.double(b),
            beta=double(nc-1),
            bt=double(nc-1),
            bt2=double(nc-1),
            bet=double(nc-1),
            cum=double(nc-1),
            like=double(1),
            DUP=F)
    z$like
  }
  call <- sys.call()
  tmp <- c("binary","multinomial","continuation-ratio","proportional-odds")
  mdl <- match(distribution <- match.arg(distribution,tmp),tmp)
  tmp <- c("independence","Markov","serial","frailty")
  dep <- match(depend <- match.arg(depend,tmp),tmp)-1
  rf <- !is.null(mu)
  respenv <- exists(paste(deparse(substitute(response))))&&inherits(response,"repeated")&&!inherits(envir,"repeated")
  if(respenv) {
    if(dim(response$response$y)[2]>1)
      stop("kalseries only handles univariate responses")
    if(!is.null(response$NAs)&&any(response$NAs))
      stop("kalseries does not handle data with NAs")
  }
  envname <- if(respenv) paste(deparse(substitute(response))) else NULL
  if(!respenv&&inherits(envir,"repeated")) {
    if(!is.null(envir$NAs)&&any(envir$NAs))
      stop("kalseries does not handle data with NAs")
    cn <- paste(deparse(substitute(response)))
    if(length(cn)>1)
      stop("only one response variable allowed")
    response <- envir
    if(dim(response$response$y)[2]>1) {
      col <- match(cn,colnames(response$response$y))
      if(is.na(col))
        stop("response variable not found")
      response$response$y <- response$response$y[,col,drop=F]
    }
  }
  if(respenv||inherits(envir,"repeated")) {
    if(!rf) {
      resp <- response
      if(is.null(ccov))
        resp$ccov <- NULL
      else
        if(inherits(ccov,"formula")) {
          tmp <- wr(ccov,data=response,expand=F)$design
          if(!is.matrix(tmp))
            stop("ccov covariate(s) not found")
          resp$ccov$ccov <- tmp[,-1,drop=F]
          rm(tmp)
        }
        else
          stop("ccov must be a W&R formula")
      if(is.null(tvcov))
        resp$tvcov <- NULL
      else
        if(inherits(tvcov,"formula")) {
          tmp <- wr(tvcov,data=response)$design
          if(!is.matrix(tmp))
            stop("tvcov covariate(s) not found")
          resp$tvcov$tvcov <- tmp[,-1,drop=F]
          rm(tmp)
        }
        else
          stop("tvcov must be a W&R formula")
    }
    else
      resp <- rmna(response$response)
    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 {
    if(inherits(response,"response"))
      resp <- response
    else
      resp <- restovec(response,times)
    if(is.null(ccov))
      nccov <- 0
    else {
      if(!inherits(ccov,"tccov")) {
        ccname <- paste(deparse(substitute(ccov)))
        if((is.matrix(ccov)&&is.null(colnames(ccov)))) {
          ccname <- paste(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]
    }
    if(is.null(tvcov))
      ttvc <- 0
    else {
      if(!inherits(tvcov,"tvcov")) {
        tvcname <- paste(deparse(substitute(tvcov)))
        if(is.list(tvcov)&&dim(tvcov[[1]])[2]>1) {
          if(is.null(colnames(tvcov[[1]]))) {
            tvcname <- paste(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)
  }
  n <- length(resp$response$y)
  nc <- length(unique(resp$response$y))
  if(min(resp$response$y)!=0||max(resp$response$y)!=nc-1)
    stop(paste("Ordinal values must be numbered from 0 to ",nc-1,".",sep=""))
  if(nc<2)
    stop("The response must at least have 2 categories.")
  else
    if(nc==2) {
      mdl <- 1
      distribution <- "binary"
    }
    else
      if(distribution=="binary") {
        mdl <- 2
        distribution <- "multinomial"
      }
  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")
  mu3 <- NULL
  if(respenv||inherits(envir,"repeated")||inherits(envir,"tccov")||inherits(envir,"tvcov")) {
    type <- if(respenv||inherits(envir,"repeated")) "repeated" else if(inherits(envir,"tccov")) "tccov" else "tvcov"
    if(is.null(envname))
      envname <- paste(deparse(substitute(envir)))
    if(inherits(mu,"formula")) {
      mu3 <- if(respenv) finterp(mu,.envir=response,.name=envname,.intercept=F) else finterp(mu,.envir=envir,.name=envname,.intercept=F)
      class(mu) <- c(class(mu),type)
    }
    else
      if(is.function(mu)) {
        tmp <- parse(text=paste(deparse(mu))[-1])
        class(mu) <- type
        mu <- if(respenv) fnenvir(mu,.envir=response,.name=envname) else fnenvir(mu,.envir=envir,.name=envname)
        mu3 <- mu
        attr(mu3,"model") <- tmp
      }
  }
  if(!is.null(preg)) {
    if(distribution=="proportional-odds")
      if(length(unique(preg[1:(nc-1)]))!=(nc-1)) {
        preg[1:(nc-1)] <- 1:(nc-1)
        warning("Changing initial preg estimates for the intercepts as they must be different.")
      }
  }
  else
    if(!is.null(ptvc))
      if(distribution=="proportional-odds")
        if(length(unique(ptvc[1:(nc-1)]))!=(nc-1)) {
          ptvc[1:(nc-1)] <- 1:(nc-1)
          warning("Changing initial preg estimates for the intercepts as they must be different.")
        }
  npreg <- length(preg)
  mu1 <- 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),.intercept=F) else finterp(mu,.envir=envir,.name=envname,.expand=is.null(preg),.intercept=F)
    npt1 <- length(attr(mu2,"parameters"))+nc-1
    if(is.character(attr(mu2,"model"))) {
      if(length(attr(mu2,"model"))==0) {
        mu2 <- NULL
        rf <- F
      }
    }
    else {
      if(npr!=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
  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))
    }
  }
  nlp <- if(is.function(mu1)) length(attr(mu1,"parameters"))+nc-1 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"))
  }
  if(rf&&!is.function(mu1))
    stop("mu must be a formula or function")
  tvc <- length(ptvc)
  init <- !is.null(pinitial)
  if(init)
    if(pinitial<=0)
      stop("Estimate of initial parameter must greater than 0")
    else
      pinitial <- log(pinitial)
  else
    if(depend!="independence")
      stop(paste("To introduce a ",depend," dependence, pinit must have an initial estimate.",sep=""))
  if(depend=="independence"||depend=="frailty")
    pdepend <- NULL
  else {
    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))
  }
  if(is.null(resp$response$times)) {
    if(depend=="serial"||depend=="Markov")
      stop("No times. Serial and Markov dependence cannot be fitted.")
    ave <- times <- 0
  }
  else {
    ave <- mean(resp$response$times)
    times <- resp$response$times-ave
  }
  if(!is.null(interaction)) {
    if(length(interaction)!=nccov)
      stop(paste(nccov,"interactions with time must be specified"))
    else
      if(any(interaction>torder))
        stop(paste("Interactions can be at most of order ",torder))
  }
  else
    interaction <- rep(0,nccov)
  npv <- torder+sum(interaction)
  if(rf&&npreg>0)
    nccov <- npreg-nc+1
  if(!rf&&nc-1+nccov+npv!=npreg)
    stop(paste(nc-1+nccov+npv,"regression estimates must be supplied"))
  y <- resp$response$y
  if(any(y<0))
    stop("All responses must be non-negative.")
  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) {
    if(tvc>0&&nccov>0)
      stop("With a mean function or formula, initial estimates must be supplied either in preg or in ptvc")
    ow <- options("warn")
    options(warn=-1)
    if(tvc>0) {
      if(length(mu1(ptvc))!=n)
        stop("The mu function must provide an estimate for each observation")
      tvc <- tvc-nc+1
    }
    else
      if(length(mu1(preg))==1) {
        if(nccov==0)
          mu1 <- function(p) rep(p[1],nind)
        else
          stop("Number of estimates does not correspond to mu function")
      }
      else
        if(length(mu1(preg))==n) { # &&nind!=n
          if(inherits(mu,"formula"))
            mu1 <- if(respenv) finterp(mu,.envir=response,.name=envname,.intercept=F) else finterp(mu,.envir=envir,.name=envname,.intercept=F)
          ptvc <- preg
          preg <- NULL
          npreg <- nccov <- 0
          tvc <- length(ptvc)-nc+1
        }
        else
          if(length(mu1(preg))!=nind)
            stop("The mu function or formula must provide an estimate for each individual")
    options(ow)
  }
  np <- nc-1+nccov+npv+tvc+init+(depend=="serial"||depend=="Markov")
  p <- c(preg,ptvc,pinitial,pdepend)
  serie <- series
  if(fscale==1)
    fscale <- serie(p)
  if(is.na(serie(p)))
    stop("Likelihood returns NAs: probably invalid initial values")
  if(optimize) {
    z0 <- nlm(serie,p=p,hessian=T,print.level=print.level,typsiz=typsiz,ndigit=ndigit,gradtol=gradtol,stepmax=stepmax,steptol=steptol,iterlim=iterlim,fscale=fscale)
    p <- z0$estimate
    like <- z0$minimum
    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)
    gradient <- z0$gradient
    iterations <- z0$iterations
    code <- z0$code
  }
  z <- {
    if(rf)
      b <- mu1(p[-(1:(nc-1))])
    .C("kord",
       p=as.double(p),
       y=as.integer(y),
       t=as.double(times),
       x=as.double(resp$ccov$ccov),
       nind=as.integer(nind),
       nobs=as.integer(nobs),
       nbs=as.integer(n),
       nc=as.integer(nc),
       nccov=as.integer(nccov),
       npv=as.integer(npv),
       init=as.integer(init),
       model=as.integer(mdl),
       dep=as.integer(dep),
       torder=as.integer(torder),
       inter=as.integer(interaction),
       tvc=as.integer(tvc),
       tvcov=as.double(resp$tvcov$tvcov),
       fit=as.integer(1),
       pred=double(n),
       rpred=double(n),
       cprob=double(n*(nc-1)),
       ppred=double(n),
       cpred=double(n*(nc-1)),
       creg=double(n*(nc-1)),
       rf=as.integer(rf),
       bbb=as.double(b),
       beta=double(nc-1),
       bt=double(nc-1),
       bt2=double(nc-1),
       bet=double(nc-1),
       cum=double(nc-1),
       like=double(1),
       DUP=F)
  }
  if(!optimize) {
    like <- z$like
    se <- rep(0,np)
    cov <- corr <- NULL
    gradient <- iterations <- code <- NULL
  }
  if(rf&&tvc>0) {
    nccov <- tvc
    tvc <- 0
  }
  if(!is.null(mu3))
    mu1 <- mu3
  z <- list(call=call,
            distribution=distribution,
            mdl=mdl,
            mu=mu1,
            npr=nc-1+nccov+tvc+torder+sum(interaction),
            depend=depend,
            torder=torder,
            interaction=interaction,
            response=resp$response,
            pred=z$pred,   # Fitted probabilities
            rpred=z$rpred, # Fitted recursive means
            cprob=z$cprob, # Fitted recursive cumulative probabilities
            ppred=z$ppred, # Fitted recursive probabilities
            cpred=z$cpred, # Fitted cumulative probabilities
            ccov=resp$ccov,
            creg=z$creg,
            tvcov=resp$tvcov,
            maxlike=like,
            aic=like+np,
            df=n-np,
            npt=np,
            npv=npv,
            init=init,
            coefficients=p,
            se=se,
            cov=cov,
            corr=corr,
            optimize=optimize,
            grad=gradient,
            iterations=iterations,
            code=code)
  class(z) <- c("kalordinal","recursive","repeated")
  return(z)
}
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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
#
#  lcr(response=NULL,frequencies=NULL,mu=NULL,ccov=NULL,tvcov=NULL,
#      times=F,pcoef=NULL,plink=NULL,distribution=NULL,direction=NULL,
#      transformation=NULL,link=NULL,envir=parent.frame(),
#      steptol=1e-04,iterlim=100)
#
#  DESCRIPTION
#
#    A function to fit various linear categorical regression models.
#

lcr <- function(response=NULL,frequencies=NULL,mu=NULL,ccov=NULL,tvcov=NULL,times=F,pcoef=NULL,plink=NULL,distribution=NULL,direction=NULL,transformation=NULL,
                link=NULL,envir=parent.frame(),steptol=1e-04,iterlim=100) {
  respenv <- exists(deparse(substitute(response)),envir=parent.frame())&&inherits(response,"repeated")&&!inherits(envir,"repeated")
  if(respenv) {
    if(dim(response$response$y)[2]>1)
      stop("Only univariate responses are handled.")
    if(!is.null(response$NAs)&&any(response$NAs))
      stop("Data with NAs cannot be handled.")
    if(response$response$type!="ordinal")
      stop("The response must be of type \'ordinal\'.")
  }
  if(!respenv&&inherits(envir,"repeated")) {
    if(!is.null(envir$NAs)&&any(envir$NAs))
      stop("Data with NAs cannot be handled.")
    cn <- deparse(substitute(response))
    if(length(grep("\"",cn))>0)
      cn <- response
    if(length(cn)>1)
      stop("Only one response variable is allowed.")
    response <- envir
    if(dim(response$response$y)[2]>1) {
      col <- match(cn,colnames(response$response$y))
      if(is.na(col))
        stop("The response variable was not found.")
      if(response$response$type[col]!="ordinal")
        stop("The response must be of type \'ordinal\'.")
      response$response$y <- response$response$y[,col,drop=F]
      rm(col)
    }
    else
      if(response$response$type!="ordinal")
        stop("The response must be of type \'ordinal\'.")
    rm(cn)
  }
  if(respenv||inherits(envir,"repeated")) {
    rm(frequencies,times,ccov,tvcov)
    if(is.null(mu)) {
      response$ccov <- response$tvcov <- NULL
      ncv <- 0
    }
    else
      if(inherits(mu,"formula")) {
        ccname <- rownames(attr(terms(mu,data=response),"factors"))
        if(!is.null(response$response$times))
          ccnames <- c(colnames(response$ccov$ccov),"times",colnames(response$tvcov$tvcov))
        else
          ccnames <- c(colnames(response$ccov$ccov),colnames(response$tvcov$tvcov))
        if(any(is.na(match(ccname,ccnames))))
          stop(paste("Covariate(s) \'",paste(ccname[is.na(match(ccname,ccnames))],collapse=", "),"\' was(were) not found.",sep=""))
        rm(ccname,ccnames)
        tmp <- wr(mu,data=response,expand=T)$design
        response$ccov <- NULL
        response$tvcov$tvcov <- tmp[,-1,drop=F]
        rm(tmp)
        ncv <- dim(response$tvcov$tvcov)[2]
      }
      else
        stop("Covariate(s) must be specified as a W&R formula.")
    rm(mu)
  }
  else {
    mu <- NULL
    if(!inherits(response,"response"))
      response <- restovec(response,times=times,weights=frequencies,type="ordinal")
    rm(frequencies,times)
    if(dim(response$y)[2]>1)
      stop("Multivariate responses are not handled.")
    if(response$type!="ordinal")
      stop("The response must be of type \'ordinal\'.")
    if(!is.null(ccov)) {
      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
            rm(tmp)
          }
        }
        ccov <- tcctomat(ccov,names=ccname)
        rm(ccname)
      }
      ncv <- dim(ccov$ccov)[2]
    }
    else
      ncv <- 0
    if(!is.null(tvcov)) {
      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]])
          rm(tvcname)
        }
        tvcov <- tvctomat(tvcov,names=tvcname)
      }
      ncv <- ncv+dim(tvcov$tvcov)[2]
    }
    response <- rmna(response=response,tvcov=tvcov,ccov=ccov)
    if(!is.null(ccov)) {
      response$ccov$ccov <- apply(response$ccov$ccov,2,function(x) rep(x,nobs(response)))
      if(!is.null(tvcov))
        response$tvcov$tvcov <- cbind(response$ccov$ccov,response$tvcov$tvcov)
      else
        response$tvcov$tvcov <- response$ccov$ccov
      response$ccov <- NULL
    }
    rm(ccov,tvcov)
    if(!is.null(response$response$times)) {
      if(is.null(response$tvcov$tvcov))
        response$tvcov$tvcov <- matrix(response$response$times,ncol=1,dimnames=list(NULL,"times"))
      else
        response$tvcov$tvcov <- cbind(response$tvcov$tvcov,matrix(response$response$times,ncol=1,dimnames=list(NULL,"times")))
      ncv <- ncv+1
    }
  }
  n <- dim(response$response$y)[1]
  nobs <- nobs(response)
  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(response$response$type!="ordinal")
    stop("The response must be of type \'ordinal\'.")
  if(is.null(weights(response))) {
    response$response$wt <- rep(1,n)
    individual <- T
  }
  else
    if(length(unique(weights(response)))==1&&unique(weights(response))==1)
      individual <- T
    else
      individual <- F
  distributions <- c("binary","multinomial","simplified multinomial","continuation-ratio","proportional odds","adjacent categories")
  directions <- c("upwards","downwards")
  transformations <- c("identity","sqrt","log","logit","cloglog","square","exp","loglog")
  if(is.null(transformation))
    transformation <- 1
  else
    transformation <- match(match.arg(transformation,transformations),transformations)
  links <- c("logit","probit","loglog","cloglog","uniform","log-normal","exponential","Pareto","Cauchy","Laplace","Levy","simplex",
             "gamma","Weibull","inverse Gauss","t","chi-square","gen logistic","gen extreme value","Box-Cox","power exponential",
             "Burr","Hjorth","beta","stable","gen gamma","gen Weibull","gen inverse Gauss","F","nc t","nc chi-square","Tukey",
             "nc beta","nc F")
  levels <- sort(unique(response$response$y))
  nc <- length(levels)
  if(nc==2) {
    if(is.null(distribution))
      distribution <- "binary"
    else {
      distribution <- match.arg(distribution,distributions)
      if(distribution!="binary") {
        distribution <- "binary"
        warning("The distribution has been changed to \'binary\' as the response has only 2 levels.")
      }
    }
    if(is.null(pcoef))
      pcoef <- rep(0,(1+ncv))
    else
      if(length(pcoef)!=(1+ncv))
        stop(paste("The initial estimates must be a vector of length ",(1+ncv),".",sep=""))
    if(is.null(link))
      link <- 1
    else
      link <- match(match.arg(link,links),links)
    if(!is.null(direction)) {
      direction <- NULL
      warning("The direction is only required by the \'continuation-ratio\' distribution.")
    }
  }
  else
    if(is.null(distribution)) {
      distribution <- "multinomial"
      if(is.null(pcoef))
        pcoef <- rep(0,(1+ncv)*(nc-1))
      else
        if(length(pcoef)!=(1+ncv)*(nc-1))
          stop(paste("The initial estimates must be a vector of length ",(1+ncv)*(nc-1),".",sep=""))
      if(!is.null(link)&&link!="logit")
        warning("The link argument is only allowed with the \'proportional odds\' and \'binary\' distributions.")
      link <- NULL
      if(!is.null(direction)) {
        direction <- NULL
        warning("The direction is only required by the \'continuation-ratio\' distribution.")
      }
    }
    else {
      distribution <- match.arg(distribution,distributions)
      if(distribution=="binary")
        stop("The option \'distribution\' can only be set to binary if the \'response\' vector has 2 levels.")
      if(distribution=="proportional odds") {
        if(is.null(pcoef))
          pcoef <- c(1:(nc-1),rep(0,ncv))
        else
          if(length(pcoef)!=(nc-1+ncv))
            stop(paste("The initial estimates must be a vector of length ",(nc-1+ncv),".",sep=""))
          else
            if(length(unique(pcoef[1:(nc-1)]))!=(nc-1)) {
              pcoef[1:(nc-1)] <- 1:(nc-1)
              warning("Changed initial estimates for the intercepts as they must be different.")
            }
        if(is.null(link))
          link <- 1
        else
          link <- match(match.arg(link,links),links)
        if(!is.null(direction)) {
          direction <- NULL
          warning("The direction is only required by the \'continuation-ratio\' distribution.")
        }
      }
      else {
        if(!is.null(link)) {
          link <- NULL
          warning("The link argument is only allowed with the \'proportional odds\' and \'binary\' distributions.")
        }
        if(distribution=="multinomial") {
          if(is.null(pcoef))
            pcoef <- rep(rep(0,1+ncv),nc-1)
          else
            if(length(pcoef)!=(1+ncv)*(nc-1))
              stop(paste("The initial estimates must be a vector of length ",(1+ncv)*(nc-1),".",sep=""))
          if(!is.null(direction)) {
            direction <- NULL
            warning("The direction is only required by the \'continuation-ratio\' distribution.")
          }
        }
        else {
          if(is.null(pcoef))
            pcoef <- rep(0,(nc-1+ncv))
          else
            if(length(pcoef)!=(nc-1+ncv))
              stop(paste("The initial estimates must be a vector of length ",(nc-1+ncv),".",sep=""))
          if(distribution=="continuation-ratio")
            if(is.null(direction))
              direction <- "upwards"
            else
              direction <- match.arg(direction,directions)
          else
            if(!is.null(direction)) {
              direction <- NULL
              warning("The direction is only required by the \'continuation-ratio\' distribution.")
            }
        }
      }
    }
  if(distribution=="binary"||distribution=="proportional odds") {
    if(is.null(plink))
      switch(links[link],
             "gamma"=plink <- 2,
             "Weibull"=plink <- 1,
             "inverse Gauss"=plink <- 1,
             "t"=plink <- 1,
             "chi-square"=plink <- 3,
             "gen logistic"=plink <- 4,
             "gen extreme value"=plink <- -2,
             "Box-Cox"=plink <- 1,
             "power exponential"=plink <- 1,
             "Burr"=plink <- 3,
             "Hjorth"=plink <- 1,
             "beta"=plink <- c(2,2),
             "stable"=plink <- c(0.5,1),
             "gen gamma"=plink <- c(2,0.5),
             "gen Weibull"=plink <- c(1,3),
             "gen inverse Gauss"=plink <- c(1,1),
             "F"=plink <- c(10,10),
             "nc t"=plink <- c(1,0),
             "nc chi-square"=plink <- c(3,0),
             "Tukey"=plink <- c(1,2,2),
             "nc beta"=plink <- c(2,2,0),
             "nc F"=plink <- c(10,10,0)
             )
    else
      if(link<=12) {
        plink <- NULL
        warning(paste("No additional parameters are needed for a ",links[link]," link.",sep=""))
      }
      else
        switch(links[link],
               "gamma"=if(length(plink)!=1) {
                 plink <- 2
                 warning("One additional parameter is needed for a \'gamma\' link.")
               }
               else {
                 if(plink<=0) {
                   plink <- 2
                   warning("The additional parameter of the \'gamma\' link (shape) was ignored as it must be greater than 0.")
                 }
               },
               "Weibull"=if(length(plink)!=1) {
                 plink <- 1
                 warning("One additional parameter is needed for a \'Weibull\' link.")
               }
               else {
                 if(plink<=0) {
                   plink <- 1
                   warning("The additional parameter of the \'Weibull\' link (shape) was ignored as it must be greater than 0.")
                 }
               },
               "inverse Gauss"=if(length(plink)!=1) {
                 plink <- 1
                 warning("One additional parameter is needed for a \'inverse Gauss\' link.")
               }
               else {
                 if(plink<0) {
                   plink <- 1
                   warning("The additional parameter of the \'inverse Gauss\' link (shape) was ignored as it must be greater than or equal to 0.")
                 }
               },
               "t"=if(length(plink)!=1) {
                 plink <- 1
                 warning("One additional parameter is needed for a \'t\' link.")
               }
               else {
                 if(plink<=0) {
                   plink <- 1
                   warning("The additional parameter of the \'t\' link (df) was ignored as it must be greater than 0.")
                 }
               },
               "chi-square"=if(length(plink)!=1) {
                 plink <- 3
                 warning("One additional parameter is needed for a \'chi-square\' link.")
               }
               else {
                 if(plink<0) {
                   plink <- 3
                   warning("The first additional parameter of the \'chi-square\' link (df) was ignored as it must be greater than or equal to 0.")
                 }
               },
               "gen logistic"=if(length(plink)!=1) {
                 plink <- 4
                 warning("One additional parameters are needed for a \'gen logistic\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 4
                   warning("The additional parameter of the \'gen logistic\' link (family) was ignored as it must be greater than 0.")
                 }
               },
               "gen extreme value"=if(length(plink)!=1) {
                 plink <- -2
                 warning("One additional parameter is needed for a \'gen extreme value\' link.")
               },
               "Box-Cox"=if(length(plink)!=1) {
                 plink <- 1
                 warning("One additional paramete is needed for a \'Box-Cox\' link.")
               },
               "power exponential"=if(length(plink)!=1) {
                 plink <- 1
                 warning("One additional parameter is needed for a \'power exponential\' link.")
               }
               else {
                 if(plink<=0) {
                   plink <- 1
                   warning("The additional parameter of the \'power exponential\' link (family) was ignored as it must be greater than 0.")
                 }
               },
               "Burr"=if(length(plink)!=1) {
                 plink <- 3
                 warning("One additional parameter is needed for a \'Burr\' link.")
               }
               else {
                 if(plink<=0) {
                   plink <- 3
                   warning("The additional parameter of the \'Burr\' link (family) was ignored as it must be greater than 0.")
                 }
               },
               "Hjorth"=if(length(plink)!=1) {
                 plink <- 1
                 warning("One additional parameter is needed for a \'Hjorth\' link.")
               },
               "beta"=if(length(plink)!=2) {
                 plink <- c(2,2)
                 warning("Three additional parameters are needed for a \'beta\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 2
                   warning("The first additional parameter of the \'beta\' link (shape1) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]<=0) {
                   plink[2] <- 2
                   warning("The second additional parameter of the \'beta\' link (shape2) was ignored as it must be greater than 0.")
                 }
               },
               "stable"=if(length(plink)!=2) {
                 plink <- c(0.5,1)
                 warning("Two additional parameters are needed for a \'stable\' link.")
               }
               else {
                 if(-1>plink[1]||plink[1]>1) {
                   plink[1] <- 0.5
                   warning("The first additional parameter of the \'stable\' link (skew) was ignored as it must take a value in (-1,1).")
                 }
                 if(plink[2]<=0||plink[2]>2) {
                   plink[2] <- 1
                   warning("The second additional parameter of the \'stable\' link (tail) was ignored as it must take a value in (0,2).")
                 }
               },
               "gen gamma"=if(length(plink)!=2) {
                 plink <- c(2,0.5)
                 warning("Two additional parameters are needed for a \'gen gamma\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 2
                   warning("The first additional parameter of the \'gen gamma\' link (shape) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]<=0) {
                   plink[2] <- 0.5
                   warning("The second additional parameter of the \'gen gamma\' link (family) was ignored as it must be greater than 0.")
                 }
               },
               "gen Weibull"=if(length(plink)!=2) {
                 plink <- c(1,3)
                 warning("Two additional parameters are needed for a \'gen Weibull\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 1
                   warning("The first additional parameter of the \'gen Weibull\' link (shape) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]<=0) {
                   plink[2] <- 3
                   warning("The second additional parameter of the \'gen Weibull\' link (family) was ignored as it must be greater than 0.")
                 }
               },
               "gen inverse Gauss"=if(length(plink)!=2) {
                 plink <- c(1,1)
                 warning("Two additional parameters are needed for a \'gen inverse Gauss\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 1
                   warning("The first additional parameter of the \'gen inverse Gauss\' link (shape) was ignored as it must be greater than 0.")
                 }
               },
               "F"=if(length(plink)!=2) {
                 plink <- c(10,10)
                 warning("Two additional parameters are needed for a \'F\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 10
                   warning("The first additional parameter of the \'F\' link (df1) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]<=0) {
                   plink[2] <- 10
                   warning("The second additional parameter of the \'F\' link (df2) was ignored as it must be greater than 0.")
                 }
               },
               "nc t"=if(length(plink)!=2) {
                 plink <- c(1,0)
                 warning("Two additional parameters are needed for a \'nc t\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 1
                   warning("The first additional parameter of the \'nc t\' link (df) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]>37.62) {
                   plink[2] <- 0
                   warning("The second additional parameter of the \'nc t\' link (ncp) was ignored as it must be smaller than 37.62.")
                 }
               },
               "nc chi-square"=if(length(plink)!=2) {
                 plink <- c(3,0)
                 warning("Two additional parameters are needed for a \'nc chi-square\' link.")
               }
               else {
                 if(plink[1]<0) {
                   plink[1] <- 3
                   warning("The first additional parameter of the \'nc chi-square\' link (df) was ignored as it must be greater than or equal to 0.")
                 }
                 if(plink[2]<0) {
                   plink[2] <- 0
                   warning("The second additional parameter of the \'nc chi-square\' link (ncp) was ignored as it must be greater than or equal to 0.")
                 }
               },
               "Tukey"=if(length(plink)!=3) {
                 plink <- c(1,2,2)
                 warning("Three additional parameters are needed for a \'Tukey\' link.")
               }
               else {
                 if(plink[1]<1) {
                   plink[1] <- 1
                   warning("The first additional parameter of the \'Tukey\' link (nranges) was ignored as it must be greater than or equal to 1.")
                 }
                 if(plink[2]<2) {
                   plink[2] <- 2
                   warning("The second additional parameter of the \'Tukey\' link (nmeans) was ignored as it must be greater than or equal to 2.")
                 }
                 if(plink[3]<2) {
                   plink[3] <- 2
                   warning("The third additional parameter of the \'Tukey\' link (df) was ignored as it must be greater than or equal to 2.")
                 }
               },
               "nc beta"=if(length(plink)!=3) {
                 plink <- c(2,2,0)
                 warning("Three additional parameters are needed for a \'nc beta\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 2
                   warning("The first additional parameter of the \'nc beta\' link (shape1) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]<=0) {
                   plink[2] <- 2
                   warning("The second additional parameter of the \'nc beta\' link (shape2) was ignored as it must be greater than 0.")
                 }
                 if(plink[3]<0) {
                   plink[3] <- 0
                   warning("The third additional parameter of the \'nc beta\' link (ncp) was ignored as it must be greater than or equal to 0.")
                 }
               },
               "nc F"=if(length(plink)!=3) {
                 plink <- c(10,10,0)
                 warning("Three additional parameters are needed for a \'nc F\' link.")
               }
               else {
                 if(plink[1]<=0) {
                   plink[1] <- 10
                   warning("The first additional parameter of the \'nc F\' link (df1) was ignored as it must be greater than 0.")
                 }
                 if(plink[2]<=0) {
                   plink[2] <- 10
                   warning("The second additional parameter of the \'nc F\' link (df2) was ignored as it must be greater than 0.")
                 }
                 if(plink[3]<0) {
                   plink[3] <- 0
                   warning("The third additional parameter of the \'nc F\' link (ncp) was ignored as it must be greater than or equal to 0.")
                 }
               }
               )
  }
  switch(distribution,
         "binary" = z <- .C("lbr",
           as.integer(response$response$y),
           as.integer(response$response$wt),
           as.double(response$tvcov$tvcov),
           as.integer(ncv),
           beta=double(1+ncv),
           as.double(pcoef),
           as.integer(n),
           double(ncv),
           iter=integer(1),
           as.integer(iterlim),
           as.double(steptol),
           double((1+ncv)*(1+ncv)),
           double(1+ncv),
           as.integer(transformation),
           as.integer(link),
           as.double(plink),
           rank=integer(1),
           double(1+ncv),
           as.integer(1:(1+ncv)),
           double(2*(1+ncv)),
           double(1+ncv),
           info=integer(1),
           as.double(diag(1,1+ncv)),
           hessian=double((1+ncv)*(1+ncv)),
           loglik=double(1),
           fit=double(n),
           pred=integer(n),
           cpred=double(n),
           DUP=F,
           PACKAGE="ordinal")[c(5,9,17,22,24:28)],
         "multinomial" = z <- .C("lmr",
           as.integer(response$response$y),
           as.integer(response$response$wt),
           as.double(response$tvcov$tvcov),
           as.integer(nc),
           as.integer(ncv),
           beta=double((nc-1)*(1+ncv)),
           as.double(pcoef),
           as.integer(n),
           integer(nc-1),
           double((nc-1)*(1+ncv)),
           iter=integer(1),
           as.integer(iterlim),
           as.double(steptol),
           double(((nc-1)*(1+ncv))*((nc-1)*(1+ncv))),
           double((nc-1)*(1+ncv)),
           double(nc-1),
           as.integer(transformation),
           double(nc-1),
           rank=integer(1),
           double((nc-1)*(1+ncv)),
           as.integer(1:((nc-1)*(1+ncv))),
           double(2*((nc-1)*(1+ncv))),
           double((nc-1)*(1+ncv)),
           info=integer(1),
           as.double(diag(1,(nc-1)*(1+ncv))),
           hessian=double(((nc-1)*(1+ncv))*((nc-1)*(1+ncv))),
           loglik=double(1),
           fit=double(n),
           pred=integer(n),
           cpred=double(n*(nc-1)),
           DUP=F,
           PACKAGE="ordinal")[c(6,11,19,24,26:30)],
         "simplified multinomial" = z <- .C("lsmr",
           as.integer(response$response$y),
           as.integer(response$response$wt),
           as.double(response$tvcov$tvcov),
           as.integer(nc),
           as.integer(ncv),
           beta=double(nc-1+ncv),
           as.double(pcoef),
           as.integer(n),
           integer(nc-1),
           double(ncv),
           iter=integer(1),
           as.integer(iterlim),
           as.double(steptol),
           double((nc-1+ncv)*(nc-1+ncv)),
           double(nc-1+ncv),
           double(nc-1),
           as.integer(transformation),
           double(nc-1),
           rank=integer(1),
           double(nc-1+ncv),
           as.integer(1:(nc-1+ncv)),
           double(2*(nc-1+ncv)),
           double(nc-1+ncv),
           info=integer(1),
           as.double(diag(1,nc-1+ncv)),
           hessian=double((nc-1+ncv)*(nc-1+ncv)),
           loglik=double(1),
           fit=double(n),
           pred=integer(n),
           cpred=double(n*(nc-1)),
           DUP=F,
           PACKAGE="ordinal")[c(6,11,19,24,26:30)],
         "continuation-ratio" = {
           if(direction=="upwards")
             z <- .C("lucrr",
                     as.integer(response$response$y),
                     as.integer(response$response$wt),
                     as.double(response$tvcov$tvcov),
                     as.integer(nc),
                     as.integer(ncv),
                     beta=double(nc-1+ncv),
                     as.double(pcoef),
                     as.integer(n),
                     integer(nc-1),
                     double(ncv),
                     iter=integer(1),
                     as.integer(iterlim),
                     as.double(steptol),
                     double((nc-1+ncv)*(nc-1+ncv)),
                     double(nc-1+ncv),
                     double(nc-1),
                     as.integer(transformation),
                     double(nc-1),
                     double(nc-1),
                     rank=integer(1),
                     double(nc-1+ncv),
                     as.integer(1:(nc-1+ncv)),
                     double(2*(nc-1+ncv)),
                     double(nc-1+ncv),
                     info=integer(1),
                     as.double(diag(1,nc-1+ncv)),
                     hessian=double((nc-1+ncv)*(nc-1+ncv)),
                     loglik=double(1),
                     fit=double(n),
                     pred=integer(n),
                     cpred=double(n*(nc-1)),
                     DUP=F,
                     PACKAGE="ordinal")[c(6,11,20,25,27:31)]
           else
             z <- .C("ldcrr",
                     as.integer(response$response$y),
                     as.integer(response$response$wt),
                     as.double(response$tvcov$tvcov),
                     as.integer(nc),
                     as.integer(ncv),
                     beta=double(nc-1+ncv),
                     as.double(pcoef),
                     as.integer(n),
                     integer(nc-1),
                     double(ncv),
                     iter=integer(1),
                     as.integer(iterlim),
                     as.double(steptol),
                     double((nc-1+ncv)*(nc-1+ncv)),
                     double(nc-1+ncv),
                     double(nc-1),
                     as.integer(transformation),
                     double(nc-1),
                     double(nc-1),
                     rank=integer(1),
                     double(nc-1+ncv),
                     as.integer(1:(nc-1+ncv)),
                     double(2*(nc-1+ncv)),
                     double(nc-1+ncv),
                     info=integer(1),
                     as.double(diag(1,nc-1+ncv)),
                     hessian=double((nc-1+ncv)*(nc-1+ncv)),
                     loglik=double(1),
                     fit=double(n),
                     pred=integer(n),
                     cpred=double(n*(nc-1)),
                     DUP=F,
                     PACKAGE="ordinal")[c(6,11,20,25,27:31)]
         },
         "proportional odds" = z <- .C("lpor",
           as.integer(response$response$y),
           as.integer(response$response$wt),
           as.double(response$tvcov$tvcov),
           as.integer(nc),
           as.integer(ncv),
           beta=double(nc-1+ncv),
           as.double(pcoef),
           as.integer(n),
           integer(nc-1),
           double(ncv),
           iter=integer(1),
           as.integer(iterlim),
           as.double(steptol),
           double((nc-1+ncv)*(nc-1+ncv)),
           double(nc-1+ncv),
           double(nc-1),
           as.integer(transformation),
           as.integer(link),
           double(nc-1),
           as.double(plink),
           double(nc-1),
           rank=integer(1),
           double(nc-1+ncv),
           as.integer(1:(nc-1+ncv)),
           double(2*(nc-1+ncv)),
           double(nc-1+ncv),
           info=integer(1),
           as.double(diag(1,nc-1+ncv)),
           hessian=double((nc-1+ncv)*(nc-1+ncv)),
           loglik=double(1),
           fit=double(n),
           pred=integer(n),
           cpred=double(n*(nc-1)),
           DUP=F,
           PACKAGE="ordinal")[c(6,11,22,27,29:33)],
         "adjacent categories" = z <- .C("lacr",
           as.integer(response$response$y),
           as.integer(response$response$wt),
           as.double(response$tvcov$tvcov),
           as.integer(nc),
           as.integer(ncv),
           beta=double(nc-1+ncv),
           as.double(pcoef),
           as.integer(n),
           integer(nc-1),
           double(ncv),
           iter=integer(1),
           as.integer(iterlim),
           as.double(steptol),
           double((nc-1+ncv)*(nc-1+ncv)),
           double(nc-1+ncv),
           double(nc-1),
           as.integer(transformation),
           double(nc-1),
           rank=integer(1),
           double(nc-1+ncv),
           as.integer(1:(nc-1+ncv)),
           double(2*(nc-1+ncv)),
           double(nc-1+ncv),
           info=integer(1),
           as.double(diag(1,nc-1+ncv)),
           hessian=double((nc-1+ncv)*(nc-1+ncv)),
           loglik=double(1),
           fit=double(n),
           pred=integer(n),
           cpred=double(n*(nc-1)),
           DUP=F,
           PACKAGE="ordinal")[c(6,11,19,24,26:30)]
         )
  covar <- matrix(-z$hessian,ncol=length(z$beta))
  se <- sqrt(diag(covar))
  corr <- covar/(se%o%se)
  res <- list(call=sys.call(),
              response=response$response,
              tvcov=response$tvcov,
              frequencies=response$response$wt,
              individual=individual,
              distribution=distribution,
              direction=direction,
              transformation=transformations[transformation],
              link=link,
              n=n, # Number of non-empty cells (or number of loops in c-code)
              nc=nc,
              ncv=ncv,
              events=if(individual) length(nobs(response)) else sum(response$response$wt), # Total number of individuals or events
              likelihood=z$loglik,
              coefficients=z$beta,
              se=se,
              covariance=covar,
              correlations=corr,
              fitted=z$fit,
              pred=z$pred,
              cpred=z$cpred,
              iterations=z$iter,
              iterlim=iterlim,
              info=z$info,
              rank=z$rank)
  class(res) <- c("lcr","repeated")
  return(res)
}
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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.
#

.First.lib <- function(lib,pkg)
  library.dynam("ordinal",pkg,lib)
require(rmutil)
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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
#
#  type(z)
#  type.default(z)
#  type.response(z)
#  extractAIC.lcr(z)
#  deviance.lcr(z)
#  fitted.lcr(z,cumulative=F)
#  names.lcr(z,data=F)
#  deviance.kalordinal(z)
#  fitted.kalordinal(z,recursive=T,cumulative=F,means=T)
#
#  DESCRIPTION
#
#    Standard methods for linear categorical regression models.
#

type <- function(z) UseMethod("type")

type.default <- function(z) {
  if(is.null(z$response)||is.null(z$response$type))
    return(NULL)
  z$response$type
}

type.response <- function(z) {
  if(is.null(z$type))
    return(NULL)
  z$type
}

extractAIC.lcr <- function(z) {
  if(z$distribution=="multinomial")
    np <- (z$nc-1)*(z$ncv+1)
  else
    np <- z$nc-1+z$ncv
  -z$likelihood+np
}

deviance.lcr <- function(z) -2*z$likelihood

fitted.lcr <- function(z,cumulative=F) {
  if(cumulative)
    print(moprofile(z,"cumulative"))
  else
    print(moprofile(a,"probability"))
}

names.lcr <- function(z,data=F) {
  if(data)
    names.repeated(z)
  else
    names.default(z)
}

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

fitted.kalordinal <- function(z,recursive=T,cumulative=F,means=T) {
  if(recursive)
    if(cumulative)
      z$cprob
    else
      if(means)
        z$rpred
      else
        z$ppred
  else
    if(cumulative)
      z$cpred
    else
      z$pred
}
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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
#
#  logit(x)
#  invlogit(x)
#
#  DESCRIPTION
#
#    Miscellaneous tools for categorical data.
#

logit <- function(x) exp(x)/(1+exp(x))

invlogit <- function(x) log(x/(1-x))
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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
#
#  plot.ordinal(z,ccov=NULL,main=NULL,xlab=NULL,ylab=NULL,xlim=NULL,
#               ylim=NULL,lty=NULL,add=F,axes=F,bty="n",at=NULL,
#               touch=F,...)
#  moprofile(z,...)
#  moprofile.default(z,curve.type="probability")
#  moprofile.lcr(z,curve.type="probability")
#  moprofile.kalordinal(z,curve.type="probability")
#  plot.moprofile(z,nind=1,observed=T,main=NULL,xlab=NULL,ylab=NULL,
#                 xlim=NULL,ylim=NULL,lty=NULL,pch=NULL,add=F,axes=F,
#                 bty="n",at=NULL,touch=F,...)
#  ioprofile(z,...)
#  ioprofile.default(z,curve.type="mean")
#  ioprofile.kalordinal(z,curve.type="mean")
#  plot.ioprofile(z,nind=1,observed=T,main=NULL,xlab=NULL,ylab=NULL,
#                 xlim=NULL,ylim=NULL,lty=NULL,pch=NULL,add=F,axes=F,
#                 bty="n",at=NULL,touch=F,...)
#  poprofile(mu,...)
#  poprofile.default(mu,pintercept,preg,pinitial=NULL,
#                    depend="independence",times=NULL,
#                    distribution="multinomial",
#                    curve.type="probability")
#  plot.poprofile(z,main="Predicted profile",xlab=NULL,ylab=NULL,
#                 xlim=NULL,ylim=NULL,lty=NULL,add=F,axes=F,bty="n",
#                 at=NULL,touch=F,...)
#  hist.repeated(z,ccov=NULL,breaks=NULL,nclass=NULL,col=NULL,
#                border=par("fg"),levels=NULL,side=F,main=NULL,
#                xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,add=F,
#                axes=F,bty="n",at=NULL,touch=F,...)
#  hist.ordinal(z,ccov=NULL,breaks=NULL,nclass=NULL,col=NULL,
#               border=par("fg"),levels=F,side=F,main=NULL,
#               xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,add=F,
#               axes=F,bty="n",at=NULL,touch=F,...)
#
#  DESCRIPTION
#
#    Utility functions for ploting ordinal repeated measurements (the
#  response, marginal profiles, individual profiles, predicted
#  profiles, and respectively the equivalent histograms).
#

plot.ordinal <- function(z,ccov=NULL,main=NULL,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,lty=NULL,add=F,axes=F,bty="n",at=NULL,touch=F,...) {
  if(inherits(z,"response")) {
    if(!is.null(ccov))
      warning("No covariates available in response objects.")
    ccov <- NULL
    z <- rmna(z)
  }
  if(!inherits(z,"repeated"))
    stop("The object z must either be a response or a repeated object.")
  nc <- length(unique(z$response$y))-1
  if(is.null(z$response$times))
    stop("No times availlable.")
  if(is.null(xlim))
    xlim <- range(z$response$times)
  if(is.null(ylim))
    ylim <- c(0,1)
  if(axes)
    xaxt <- yaxt <- "s"
  else
    xaxt <- yaxt <- "n"
  if(touch) {
    xaxs <- yaxs <- "i"
    xaxt <- yaxt <- "n"
  }
  else
    xaxs <- yaxs <- "r"
  if(!is.null(lty))
    if(length(lty)!=1&&length(lty)!=nc) {
      warning(paste("Dropping contents of option lty as it must have lentgh ",nc,".",sep=""))
      lty <- NULL
    }
    else
      if(length(lty)==1)
        lty <- rep(lty,nc)
  if(!is.null(ccov)) {
    if(length(ccov)!=1) {
      warning("Only using first covariate.")
      ccov <- ccov[1]
    }
    cov.name <- ccov
    ccov <- covariates(z,names=ccov)
    if(is.data.frame(ccov))
      ccov <- codes(ccov)
    np <- unique(ccov)
    ns <- length(np)
    if(ns>25)
      stop("The covariate chosen has more than 25 levels (too many categories for plotting).")
    if(length(ccov)==length(nobs(z)))
      ccov <- rep(ccov,nobs(z))
    if(is.null(main))
      main <- paste(cov.name," ",np,".",sep="")
    else
      if(length(main)!=1&&length(main)!=ns) {
        warning(paste("Dropping contents of option main as it must have lentgh ",ns,".",sep=""))
        main <- paste(cov.name," ",np,".",sep="")
      }
      else
        if(length(main)==1)
          main <- rep(main,ns)
    if(is.null(xlab))
      xlab <- rep("Time",ns)
    else
      if(length(xlab)!=1&&length(xlab)!=ns) {
        warning(paste("Dropping contents of option xlab as it must have lentgh ",ns,".",sep=""))
        xlab <- rep("Time",ns)
      }
      else
        if(length(xlab)==1)
          xlab <- rep(xlab,ns)
    if(is.null(ylab))
      ylab <- rep("Cumulative probabilities",ns)
    else
      if(length(ylab)!=1&&length(ylab)!=ns) {
        warning(paste("Dropping contents of option ylab as it must have lentgh ",ns,".",sep=""))
        ylab <- rep("Cumulative probabilities",ns)
      }
      else
        if(length(ylab)==1)
          ylab <- rep(ylab,ns)
  }
  else {
    ccov <- np <- ns <- 1
    if(is.null(main))
      main <- "Cumulative time profile."
    else
      if(length(main)!=1) {
        warning("Dropping contents of option main as it must have lentgh 1.")
        main <- "Cumulative time profile."
      }
    if(is.null(xlab))
      xlab <- "Time"
    else
      if(length(xlab)!=1) {
        warning("Dropping contents of option xlab as it must have lentgh 1.")
        xlab <- "Time"
      }
    if(is.null(ylab))
      ylab <- "Cumulative probabilities"
    else
      if(length(ylab)!=1) {
        warning("Dropping contents of option ylab as it must have lentgh 1.")
        ylab <- "Cumulative probabilities"
      }
  }
  if(ns!=1&&!add)
    if(ns==2)
      oldpar <- par(mfrow=c(1,2))
    else {
      ng <- ceiling(sqrt(ns))
      ng <- c(ceiling(ns/ng),ng)
      oldpar <- par(mfrow=ng)
    }
  i <- 0
  for(k in np) {
    i <- i+1
    times <- sort(unique(times(z)[ccov==k]))
    ntimes <- length(times)
    if(is.null(weights(z)))
      tab <- table(response(z)[ccov==k],times(z)[ccov==k])
    else
      tab <- tapply(weights(z)[ccov==k],list(as.vector(response(z)[ccov==k]),times(z)[ccov==k]),sum)
    lt <- if(is.null(lty)) 1 else lty[1]
    plot(times,apply(matrix(tab[1,],ncol=ntimes),2,sum)/apply(tab,2,sum),main=main[i],xlim=xlim,ylim=ylim,ylab=ylab[i],xlab=xlab[i],type="l",lty=lt,
         xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
    if(!axes||touch) {
      if(!axes)
        if(is.null(at))
          axis(1,at=times,...)
        else
          axis(1,at=at,...)
      else
        axis(1,at=pretty(xlim),...)
      if(touch)
        axis(2,at=pretty(ylim),...)
      else
        axis(2,...)
    }
    if(nc>1)
      for(j in 2:nc) {
        lt <- if(is.null(lty)) lt%%6+1 else lty[j]
        lines(times,apply(matrix(tab[1:j,],ncol=ntimes),2,sum)/apply(tab,2,sum),lty=lt,...)
      }
  }
  if(ns!=1&&!add)
    par(oldpar)
}

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

moprofile.default <- function(z,curve.type="probability") {
  if(inherits(z$response,"response"))
    stop("Object z must contain the response.")
  tmp <- c("probability","cumulative")
  curve.type <- match.arg(curve.type,tmp)
  if(curve.type=="probability")
    if(is.null(z$pred))
      stop("Highest probabilities not available (pred).")
    else
      if(is.null(z$cpred))
        stop("Cumulative probabilities not available (cpred).")
  z$curve.type <- curve.type
  class(z) <- c("moprofile",class(z))
  invisible(z)
}

moprofile.lcr <- function(z,curve.type="probability") {
  z$curve.type <- match.arg(curve.type,c("probability","cumulative"))
  class(z) <- c("moprofile",class(z))
  invisible(z)
}

moprofile.kalordinal <- function(z,curve.type="probability") {
  z$curve.type <- match.arg(curve.type,c("probability","cumulative"))
  class(z) <- c("moprofile",class(z))
  invisible(z)
}

plot.moprofile <- function(z,nind=1,observed=T,main=NULL,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,lty=NULL,pch=NULL,add=F,axes=F,bty="n",at=NULL,touch=F,...) {
  if(max(nind)>length(nobs(z))||min(nind)<1)
    stop("No such individual.")
  nc <- length(unique(z$response$y))-1
  if(is.null(z$response$times))
    stop("Times are required to perform a marginal ordinal profile.")
  ii <- covind(z)
  ns <- length(nind)
  if(is.null(xlim))
    xlim <- range(z$response$times)
  if(axes)
    xaxt <- yaxt <- "s"
  else
    xaxt <- yaxt <- "n"
  if(touch) {
    xaxs <- yaxs <- "i"
    xaxt <- yaxt <- "n"
  }
  else
    xaxs <- yaxs <- "r"
  if(z$curve.type=="probability") {
    if(is.null(main))
      if(ns==1)
        main <- paste("Individual ",nind,".",sep="")
      else
        main <- paste("Individuals ",paste(nind,collapse=", "),".",sep="")
    else
      if(length(main)!=1) {
        warning("Dropping contents of option main as it must have lentgh 1.")
        if(ns==1)
          main <- paste("Individual ",nind,".",sep="")
        else
          main <- paste("Individuals ",paste(nind,collapse=", "),".",sep="")
      }
    if(is.null(xlab))
      xlab <- "Time"
    else
      if(length(xlab)!=1) {
        warning("Dropping contents of option xlab as it must have lentgh 1.")
        xlab <- "Time"
      }
    if(is.null(ylab))
      ylab <- "Highest probabilities"
    else
      if(length(ylab)!=1) {
        warning("Dropping contents of option ylab as it must have lentgh 1.")
        ylab <- "Highest probabilities"
      }
    if(is.null(ylim))
      ylim <- c(0,nc)
    if(!is.null(lty))
      if(length(lty)!=1&&length(lty)!=ns) {
        warning(paste("Dropping contents of option lty as it must have lentgh ",ns,".",sep=""))
        lty <- NULL
      }
      else
        if(length(lty)==1)
          lty <- rep(lty,ns)
    if(!is.null(pch))
      if(length(pch)!=1&&length(pch)!=ns) {
        warning(paste("Dropping contents of option pch as it must have lentgh ",ns,".",sep=""))
        pch <- NULL
      }
      else
        if(length(pch)==1)
          pch <- rep(pch,ns)
    first <- !add
    lt <- pc <- j <- 0
    for(i in nind) {
      j <- j+1
      lt <- if(is.null(lty)) lt%%6+1 else lty[j]
      if(first) {
        plot(z$response$times[ii==i],z$pred[ii==i],main=main,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,type="l",lty=lt,xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,
             bty=bty,...)
        if(!axes||touch) {
          if(!axes)
            if(is.null(at))
              axis(1,at=z$response$times[ii==i],...)
            else
              axis(1,at=at,...)
          else
            axis(1,at=pretty(xlim),...)
          if(touch)
            axis(2,at=pretty(ylim),...)
          else
            axis(2,...)
        }
        first <- F
      }
      else
        lines(z$response$times[ii==i],z$pred[ii==i],lty=lt,...)
      if(observed) {
        pc <- if(is.null(pch)) pc%%6+1 else pch[j]
        points(z$response$times[ii==i],z$response$y[ii==i],pch=pc,...)
      }
    }
  }
  else {
    if(ns>25)
      stop("Too many individuals were chosen to be plotted at once (maximum 25).")
    if(ns!=1&&!add)
      if(ns==2)
        oldpar <- par(mfrow=c(1,2))
      else {
        ng <- ceiling(sqrt(ns))
        ng <- c(ceiling(ns/ng),ng)
        oldpar <- par(mfrow=ng)
      }
    if(is.null(main))
      main <- paste("Individual ",nind,".",sep="")
    else
      if((length(main)!=1&&length(main)!=ns)) {
        warning(paste("Dropping contents of option main as it must have lentgh ",ns,".",sep=""))
        main <- paste("Individual ",nind,".",sep="")
      }
      else
        if(length(main)==1)
          main <- rep(main,ns)
    if(is.null(xlab))
      xlab <- rep("Time",ns)
    else
      if((length(xlab)!=1&&length(xlab)!=ns)){
        warning(paste("Dropping contents of option xlab as it must have lentgh ",ns,".",sep=""))
        xlab <- rep("Time",ns)
      }
      else
        if(length(xlab)==1)
          xlab <- rep(xlab,ns)
    if(is.null(ylab))
      ylab <- rep("Cumulative probabilities",ns)
    else
      if((length(ylab)!=1&&length(ylab)!=ns)) {
        warning(paste("Dropping contents of option ylab as it must have lentgh ",ns,".",sep=""))
        ylab <- rep("Cumulative probabilities",ns)
      }
      else
        if(length(ylab)==1)
          ylab <- rep(ylab,ns)
    if(is.null(ylim))
      ylim <- c(0,1)
    if(!is.null(lty))
      if(length(lty)!=1&&length(lty)!=nc) {
        warning(paste("Dropping contents of option lty as it must have lentgh ",nc,".",sep=""))
        lty <- NULL
      }
      else
        if(length(lty)==1)
          lty <- rep(lty,nc)
    pred <- matrix(z$cpred,nrow=nc,byrow=T)
    j <- 0
    for(i in nind) {
      j <- j+1
      lt <- if(is.null(lty)) 1 else lty[1]
      plot(z$response$times[ii==i],pred[1,ii==i],main=main[j],xlim=xlim,ylim=ylim,xlab=xlab[j],ylab=ylab[j],type="l",lty=lt,
           xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
      if(!axes||touch) {
        if(!axes)
          if(is.null(at))
            axis(1,at=z$response$times[ii==i],...)
          else
            axis(1,at=at,...)
        else
          axis(1,at=pretty(xlim),...)
        if(touch)
          axis(2,at=pretty(ylim),...)
        else
          axis(2,...)
      }
      if(nc>1)
        for(k in 2:nc) {
          lt <- if(is.null(lty)) lt%%6+1 else lty[k]
          lines(z$response$times[ii==i],pred[k,ii==i],lty=lt,...)
        }
      if(observed)
        mtext(z$response$y[ii==i],at=z$response$times[ii==i],...)
    }
    if(ns!=1&&!add)
      par(oldpar)
  }
}

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

ioprofile.default <- function(z,curve.type="mean") {
  if(inherits(z$response,"response"))
    stop("Object z must contain the response.")
  tmp <- c("mean","probability","both","cumulative")
  curve.type <- match.arg(curve.type,tmp)
  if(!inherits(z,"recursive"))
    stop("The object must have class, recursive")
  else {
    if(curve.type=="mean"||curve.type=="both") {
      if(is.null(z$rpred))
        stop("Recursive means not available (rpred).")
    }
    if(curve.type=="probability"||curve.type=="both") {
      if(is.null(z$ppred))
        stop("Recursive highest probabilities not available (ppred).")
    }
    if(curve.type=="cumulative") {
      if(is.null(z$cprob))
        stop("Recursive cumulative probabilities not available (cprob).")
    }
  }
  z$curve.type <- curve.type
  class(z) <- c("ioprofile",class(z))
  invisible(z)
}

ioprofile.kalordinal <- function(z,curve.type="mean") {
  z$curve.type <- match.arg(curve.type,c("mean","probability","both","cumulative"))
  class(z) <- c("ioprofile",class(z))
  invisible(z)
}

plot.ioprofile <- function(z,nind=1,observed=T,main=NULL,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,lty=NULL,pch=NULL,add=F,axes=F,bty="n",at=NULL,touch=F,...) {
  if(max(nind)>length(nobs(z))||min(nind)<1)
    stop("No such individual.")
  nc <- length(unique(z$response$y))-1
  if(is.null(z$response$times))
    stop("Times are required to perform an individual ordinal profile.")
  ii <- covind(z)
  ns <- length(nind)
  if(is.null(xlim))
    xlim <- range(z$response$times)
  if(axes)
    xaxt <- yaxt <- "s"
  else
    xaxt <- yaxt <- "n"
  if(touch) {
    xaxs <- yaxs <- "i"
    xaxt <- yaxt <- "n"
  }
  else
    xaxs <- yaxs <- "r"
  if(z$curve.type=="cumulative") {
    if(is.null(ylim))
      ylim <- c(0,1)
    if(!is.null(lty))
      if(length(lty)!=1&&length(lty)!=nc) {
        warning(paste("Dropping contents of option lty as it must have lentgh ",nc,".",sep=""))
        lty <- NULL
      }
      else
        if(length(lty)==1)
          lty <- rep(lty,nc)
  }
  else {
    if(is.null(ylim))
      ylim <- c(0,nc)
    if(z$curve.type=="both") {
      if(is.null(lty))
        lty <- 1:2
      else
        if((length(lty)!=1&&length(lty)!=2)) {
          warning("Dropping contents of option lty as it must have lentgh 2.")
          lty <- 1:2
        }
        else
          if(length(lty)==1)
            lty <- c(lty,lty+1)
      if(is.null(pch))
        pch <- 1
      else
        if(length(pch)!=1) {
          warning("Dropping contents of option pch as it must have lentgh 1.")
          pch <- 1
        }
    }
    else {
      if(is.null(main))
        if(ns==1)
          main <- paste("Individual ",nind,".",sep="")
        else
          main <- paste("Individuals ",paste(nind,collapse=", "),".",sep="")
      else
        if(length(main)!=1) {
          warning("Dropping contents of option main as it must have lentgh 1.")
          if(ns==1)
            main <- paste("Individual ",nind,".",sep="")
          else
            main <- paste("Individuals ",paste(nind,collapse=", "),".",sep="")
        }
      if(is.null(xlab))
        xlab <- "Time"
      else
        if(length(xlab)!=1) {
          warning("Dropping contents of option xlab as it must have lentgh 1.")
          xlab <- "Time"
        }
      if(!is.null(lty))
        if(length(lty)!=1&&length(lty)!=ns) {
          warning(paste("Dropping contents of option lty as it must have lentgh ",ns,".",sep=""))
          lty <- NULL
        }
        else
          if(length(lty)==1)
            lty <- rep(lty,ns)
      if(!is.null(pch))
        if(length(pch)!=1&&length(pch)!=ns) {
          warning(paste("Dropping contents of option pch as it must have lentgh ",ns,".",sep=""))
          pch <- NULL
        }
        else
          if(length(pch)==1)
            pch <- rep(pch,ns)
    }
  }
  if(z$curve.type=="mean") {
    if(is.null(ylab))
      ylab <- "Recursive means"
    else
      if(length(ylab)!=1) {
        warning("Dropping contents of option ylab as it must have lentgh 1.")
        ylab <- "Recursive means"
      }
    first <- !add
    lt <- pc <- j <- 0
    for(i in nind) {
      j <- j+1
      lt <- if(is.null(lty)) lt%%6+1 else lty[j]
      if(first) {
        plot(z$response$times[ii==i],z$rpred[ii==i],main=main,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,type="l",lty=lt,
             xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
        if(!axes||touch) {
          if(!axes)
            if(is.null(at))
              axis(1,at=z$response$times[ii==i],...)
            else
              axis(1,at=at,...)
          else
            axis(1,at=pretty(xlim),...)
          if(touch)
            axis(2,at=pretty(ylim),...)
          else
            axis(2,...)
        }
        first <- F
      }
      else
        lines(z$response$times[ii==i],z$rpred[ii==i],lty=lt,...)
      if(observed) {
        pc <- if(is.null(pch)) pc%%6+1 else pch[j]
        points(z$response$times[ii==i],z$response$y[ii==i],pch=pc,...)
      }
    }
  }
  else
    if(z$curve.type=="probability") {
      if(is.null(ylab))
        ylab <- "Recursive highest probabilities"
      else
        if(length(ylab)!=1) {
          warning("Dropping contents of option ylab as it must have lentgh 1.")
          ylab <- "Recursive highest probabilities"
        }
      first <- !add
      lt <- pc <- j <- 0
      for(i in nind) {
        j <- j+1
        lt <- if(is.null(lty)) lt%%6+1 else lty[j]
        if(first) {
          plot(z$response$times[ii==i],z$ppred[ii==i],main=main,xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,type="l",lty=lt,
               xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
          if(!axes||touch) {
            if(!axes)
              if(is.null(at))
                axis(1,at=z$response$times[ii==i],...)
              else
                axis(1,at=at,...)
            else
              axis(1,at=pretty(xlim),...)
            if(touch)
              axis(2,at=pretty(ylim),...)
            else
              axis(2,...)
          }
          first <- F
        }
        else
          lines(z$response$times[ii==i],z$ppred[ii==i],lty=lt,...)
        if(observed) {
          pc <- if(is.null(pch)) pc%%6+1 else pch[j]
          points(z$response$times[ii==i],z$response$y[ii==i],pch=pc,...)
        }
      }
    }
    else {
      if(ns>25)
        stop("Too many individuals were chosen to be plotted at once (maximum 25).")
      if(ns!=1&&!add)
        if(ns==2)
          oldpar <- par(mfrow=c(1,2))
        else {
          ng <- ceiling(sqrt(ns))
          ng <- c(ceiling(ns/ng),ng)
          oldpar <- par(mfrow=ng)
        }
      if(is.null(main))
        main <- paste("Individual ",nind,".",sep="")
      else
        if((length(main)!=1&&length(main)!=ns)) {
          warning(paste("Dropping contents of option main as it must have lentgh ",ns,".",sep=""))
          main <- paste("Individual ",nind,".",sep="")
        }
        else
          if(length(main)==1)
            main <- rep(main,ns)
      if(is.null(xlab))
        xlab <- rep("Time",ns)
      else
        if((length(xlab)!=1&&length(xlab)!=ns)){
          warning(paste("Dropping contents of option xlab as it must have lentgh ",ns,".",sep=""))
          xlab <- rep("Time",ns)
        }
        else
          if(length(xlab)==1)
            xlab <- rep(xlab,ns)
      if(z$curve.type=="both") {
        if(is.null(ylab))
          ylab <- rep("Recursive predictions",ns)
        else
          if((length(ylab)!=1&&length(ylab)!=ns)) {
            warning(paste("Dropping contents of option ylab as it must have lentgh ",ns,".",sep=""))
            ylab <- rep("Recursive predictions",ns)
          }
          else
            if(length(ylab)==1)
              ylab <- rep(ylab,ns)
        j <- 0
        for(i in nind) {
          j <- j+1
          plot(z$response$times[ii==i],z$rpred[ii==i],main=main[j],xlim=xlim,ylim=ylim,xlab=xlab[j],ylab=ylab[j],type="l",lty=lty[1],
               xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
          if(!axes||touch) {
            if(!axes)
              if(is.null(at))
                axis(1,at=z$response$times[ii==i],...)
              else
                axis(1,at=at,...)
            else
              axis(1,at=pretty(xlim),...)
            if(touch)
              axis(2,at=pretty(ylim),...)
            else
              axis(2,...)
          }
          lines(z$response$times[ii==i],z$ppred[ii==i],lty=lty[2],...)
          if(observed)
            points(z$response$times[ii==i],z$response$y[ii==i],pch=pch,...)
        }
      }
      else {
        if(is.null(ylab))
          ylab <- rep("Recursive cumulative probabilities",ns)
        else
          if((length(ylab)!=1&&length(ylab)!=ns)) {
            warning(paste("Dropping contents of option ylab as it must have lentgh ",ns,".",sep=""))
            ylab <- rep("Recursive cumulative probabilities",ns)
          }
          else
            if(length(ylab)==1)
              ylab <- rep(ylab,ns)
        pred <- matrix(z$cprob,nrow=nc,byrow=T)
        j <- 0
        for(i in nind) {
          j <- j+1
          lt <- if(is.null(lty)) 1 else lty[1]
          plot(z$response$times[ii==i],pred[1,ii==i],main=main[j],xlim=xlim,ylim=ylim,xlab=xlab[j],ylab=ylab[j],type="l",lty=lt,
               xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
          if(!axes||touch) {
            if(!axes)
              if(is.null(at))
                axis(1,at=z$response$times[ii==i],...)
              else
                axis(1,at=at,...)
            else
              axis(1,at=pretty(xlim),...)
            if(touch)
              axis(2,at=pretty(ylim),...)
            else
              axis(2,...)
          }
          if(nc>1)
            for(k in 2:nc) {
              lt <- if(is.null(lty)) lt%%6+1 else lty[k]
              lines(z$response$times[ii==i],pred[k,ii==i],lty=lt,...)
            }
          if(observed)
            mtext(z$response$y[ii==i],at=z$response$times[ii==i],...)
        }
      }
      if(ns!=1&&!add)
        par(oldpar)
    }
}

poprofile <- function(mu,...) UseMethod("poprofile")

poprofile.default <- function(mu,pintercept,preg,pinitial=NULL,depend="independence",times=NULL,distribution="multinomial",curve.type="probability") {
  call <- sys.call()
  tmp <- c("binary","multinomial","continuation-ratio","proportional-odds")
  mdl <- match(distribution <- match.arg(distribution,tmp),tmp)
  tmp <- c("independence","Markov","serial","frailty")
  dep <- match(depend <- match.arg(depend,tmp),tmp)-1
  tmp <- c("probability","cumulative")
  curve.type <- match.arg(curve.type,tmp)
  if(is.null(times)&&(depend=="serial"||depend=="Markov"))
    stop("No times. Serial and Markov dependence cannot be fitted.")
  init <- !is.null(pinitial)
  if(init)
    if(pinitial<=0)
      stop("Estimate of initial parameter must greater than 0")
    else
      pinitial <- log(pinitial)
  else
    if(depend!="independence")
      stop(paste("To introduce a ",depend," dependence, pinit must have an initial estimate.",sep=""))
  nc <- length(pintercept)+1
  b <- mu(preg)
  n <- length(b)
  p <- c(pintercept,preg,pinitial)
  x <- .C("pkord",
          p=as.double(p),
          nbs=as.integer(n),
          nc=as.integer(nc),
          init=as.integer(init),
          model=as.integer(mdl),
          dep=as.integer(dep),
          preg=as.integer(length(preg)),
          cpred=double(n*(nc-1)),
          pred=double(n),
          bbb=as.double(b),
          beta=double((nc-1)),
          bt=double((nc-1)),
          bt2=double((nc-1)),
          bet=double((nc-1)),
          cum=double((nc-1)),
          DUP=F)
  z <- list(call=call,
            curve.type=curve.type,
            distribution=distribution,
            mdl=mdl,
            mu=mu,
            npr=nc-1+length(preg),
            depend=depend,
            response=list(y=0:(nc-1),times=times,nobs=n),
            cpred=x$cpred,
            pred=x$pred,
            init=init,
            coefficients=p)
  class(z) <- "poprofile"
  invisible(z)
}

plot.poprofile <- function(z,main="Predicted profile",xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,lty=NULL,add=F,axes=F,bty="n",at=NULL,touch=F,...)
  plot.moprofile(z,nind=1,observed=F,main=main,xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,lty=lty,pch=NULL,add=add,axes=axes,bty=bty,at=at,touch=touch,...)

hist.repeated <- function(z,ccov=NULL,breaks=NULL,nclass=NULL,col=NULL,border=par("fg"),levels=NULL,side=F,main=NULL,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,add=F,
                          axes=F,bty="n",at=NULL,touch=F,...) {
  if(z$response$type!="ordinal")
    stop("To perform a histogram, the response in the `repeated' object must be of type `ordinal'.")
  hist.ordinal(z=z,ccov=ccov,breaks=breaks,nclass=nclass,col=col,border=border,levels=levels,side=side,main=main,xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,add=add,
               axes=axes,bty=bty,at=at,touch=touch,...)
}

hist.ordinal <- function(z,ccov=NULL,breaks=NULL,nclass=NULL,col=NULL,border=par("fg"),levels=F,side=F,main=NULL,xlab=NULL,ylab=NULL,xlim=NULL,ylim=NULL,add=F,
                         axes=F,bty="n",at=NULL,touch=F,...) {
  if(inherits(z,"response")) {
    if(!is.null(ccov))
      warning("No covariates available in response objects.")
    ccov <- NULL
    z <- rmna(z)
  }
  if(!inherits(z,"repeated"))
    stop("The object z must either be a response or a repeated object.")
  nc <- length(unique(z$response$y))
  if(is.null(times(z))) {
    all.times <- 1
    m <- 1
  }
  else {
    all.times <- sort(unique(times(z)))
    m <- length(all.times)
  }
  use.br <- !is.null(breaks)
  if(use.br) {
    if(!is.null(nclass))
      warning("`nclass' not used when `breaks' specified")
  }
  else
    if(!is.null(nclass)&&length(nclass)==1)
      breaks <- nclass
  rm(nclass)
  use.br <- use.br&&length(breaks)>1
  if(use.br) {
    if(any(is.na(breaks))|length(breaks)!=m+1)
      stop(paste("Invalid or wrong number of breaks, no NAs are allowed and ",m+1," breaks are required.",sep=""))
    breaks <- sort(breaks)
    if(length(unique(diff(breaks)))!=1)
      stop("The breaks are invalid as they must all delimit the same interval.")
    nB <- length(breaks)
  }
  else {
    if(is.null(breaks))
      breaks <- all.times
    else
      if(is.na(breaks)|breaks!=m)
        stop(paste("Invalid number of breaks, ",m," are required.",sep=""))
      else
        breaks <- all.times
    if(m>1)
      breaks <- c(breaks,max(breaks)+unique(diff(breaks)))
    else
      breaks <- c(all.times,all.times+1)
    breaks <- breaks-(unique(diff(breaks))/2)
    nB <- length(breaks)
  }
  if(is.logical(side)&&side)
    side <- "l"
  else
    if(!is.logical(side)&&side!="l"&&side!="r") {
      warning("Droping contents of option `side' as it must either be `l' (for left side) or `r' (for right side).")
      side <- "l"
    }
  if(is.null(levels))
    levels <- 1:nc
  else
    if(length(levels)!=nc) {
      warning(paste("Wrong number of levels, vector `levels' should be of length ",nc,".",sep=""))
      levels <- 1:nc
    }
  if(is.null(col))
    col <- gray(seq(1,0,length=nc+1))[-(nc+1)]
  else
    if(length(col)!=nc) {
      warning(paste("Wrong number of colors, vector `col' should be of length ",nc,".",sep=""))
      col <- gray(seq(1,0,length=nc+1))[-(nc+1)]
    }
  if(is.null(xlim))
    xlim <- c(min(breaks)-unique(diff(breaks))/2,max(breaks)+unique(diff(breaks))/2)
  if(is.null(ylim))
    ylim <- 0:1
  if(axes)
    xaxt <- yaxt <- "s"
  else
    xaxt <- yaxt <- "n"
  if(touch) {
    xaxs <- yaxs <- "i"
    xaxt <- yaxt <- "n"
  }
  else
    xaxs <- yaxs <- "r"
  if(!is.null(ccov)) {
    if(length(ccov)!=1) {
      warning("Only using first covariate.")
      ccov <- ccov[1]
    }
    cov.name <- ccov
    ccov <- covariates(z,names=ccov)
    if(is.data.frame(ccov))
      ccov <- codes(ccov)
    np <- unique(ccov)
    ns <- length(np)
    if(ns>25)
      stop("The covariate chosen has more than 25 levels (too many categories for plotting).")
    if(length(ccov)==length(nobs(z)))
      ccov <- rep(ccov,nobs(z))
    if(is.null(main))
      main <- paste(cov.name," ",np,".",sep="")
    else
      if(length(main)!=1&&length(main)!=ns) {
        warning(paste("Dropping contents of option main as it must have lentgh ",ns,".",sep=""))
        main <- paste(cov.name," ",np,".",sep="")
      }
      else
        if(length(main)==1)
          main <- rep(main,ns)
    if(is.null(xlab))
      xlab <- rep("Time",ns)
    else
      if(length(xlab)!=1&&length(xlab)!=ns) {
        warning(paste("Dropping contents of option xlab as it must have lentgh ",ns,".",sep=""))
        xlab <- rep("Time",ns)
      }
      else
        if(length(xlab)==1)
          xlab <- rep(xlab,ns)
    if(is.null(ylab))
      ylab <- rep("Cumulative probabilities",ns)
    else
      if(length(ylab)!=1&&length(ylab)!=ns) {
        warning(paste("Dropping contents of option ylab as it must have lentgh ",ns,".",sep=""))
        ylab <- rep("Cumulative probabilities",ns)
      }
      else
        if(length(ylab)==1)
          ylab <- rep(ylab,ns)
  }
  else {
    ccov <- np <- ns <- 1
    if(is.null(main))
      main <- "Cumulative histograms over time."
    if(is.null(xlab))
      xlab <- "Time"
    if(is.null(ylab))
      ylab <- "Cumulative probabilities"
  }
  if(ns!=1&&!add)
    if(ns==2)
      oldpar <- par(mfrow=c(1,2))
    else {
      ng <- ceiling(sqrt(ns))
      ng <- c(ceiling(ns/ng),ng)
      oldpar <- par(mfrow=ng)
    }
  i <- 0
  for(k in np) {
    i <- i+1
    times <- sort(unique(times(z)[ccov==k]))
    ntimes <- length(times)
    if(is.null(weights(z)))
      tab <- table(response(z)[ccov==k],times(z)[ccov==k])
    else
      tab <- tapply(weights(z)[ccov==k],list(as.vector(response(z)[ccov==k]),times(z)[ccov==k]),sum)
    plot.new()
    plot.window(xlim=xlim,ylim=ylim,xaxt=xaxt,yaxt=yaxt,xaxs=xaxs,yaxs=yaxs,bty=bty,...)
    title(main=main[i],xlab=xlab[i],ylab=ylab[i],...)
    if(!axes||touch) {
      if(!axes)
        if(is.null(at))
          axis(1,at=times,...)
        else
          axis(1,at=at,...)
      else
        axis(1,at=pretty(xlim),...)
      if(touch)
        axis(2,at=pretty(ylim),...)
      else
        axis(2,...)
    }
    bottom <-  rep(0,m)
    for(j in 1:nc) {
      x <- apply(matrix(tab[1:j,],ncol=ntimes),2,sum)/apply(tab,2,sum)
      if(side=="l")
        height <- (bottom[1]+x[1])/2
      else
        if(side=="r")
          height <- (bottom[m]+x[m])/2
      rect(breaks[-nB],bottom,breaks[-1],x,col=col[j],border=border)
      if(side=="l")
        text(x=xlim[1],y=height,label=levels[j])
      else
        if(side=="r")
          text(x=xlim[2],y=height,label=levels[j])
      bottom <- x
    }
  }
  if(ns!=1&&!add)
    par(oldpar)
}
#
#  Ordinal: A Library of Ordinal Models
#  Copyright (C) 1998, 1999, 2000, 2001 P.J. 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
#
#  print.lcr(z,digits=max(3,.Options$digits-3),print.gap=2)
#  print.kalordinal(z,digits=max(3,.Options$digits-3))
#  print.moprofile(z,nind=NULL,digits=max(3,.Options$digits-3),
#                  print.gap=2)
#  print.listmatrix(x,row.names=NULL,digits=max(3,.Options$digits-3),
#                   print.gap=2)
#
#  DESCRIPTION
#
#    Printing functions for linear categorical regression models,
#  marginal ordinal time profiles, and lists of matrices.
#

print.lcr <- function(z,digits=max(3,.Options$digits-3),print.gap=2) {
  links <- c("logit","probit","loglog","cloglog","uniform","log-normal","exponential","Pareto","Cauchy","Laplace","Levy","simplex",
             "gamma","Weibull","inverse Gauss","t","chi-square","generalized logistic","generalized extreme value","Box-Cox",
             "power exponential","Burr","Hjorth","beta","stable","generalized gamma","generalized Weibull","generalized inverse Gauss",
             "F","non-central t","non-central chi-square","Tukey","non-central beta","non-central F")
  cat("\nCall:",deparse(z$call),"",sep="\n")
  if(z$individual) {
    cat("Individual data.\n")
    cat("Total number of individuals: ",z$events,"\n",sep="")
    cat("Number of observations:      ",z$n,"\n\n",sep="")
  }
  else {
    cat("Frequency table.\n")
    cat("Number of non-empty cells: ",z$n,"\n",sep="")
    cat("Total number of events:    ",z$events,"\n\n",sep="")
  }
  switch(z$distribution,
         "binary" = {
           cat("Binary distribution.\n")
           cat("Transformation: ",z$transformation,".\n",sep="")
           cat("Link:           ",links[z$link],".\n",sep="")
         },
         "multinomial" = {
           cat("Multinomial distribution.\n")
           cat("Transformation: ",z$transformation,".\n",sep="")
         },
         "simplified multinomial" = {
           cat("Simplified multinomial distribution.\n")
           cat("Transformation: ",z$transformation,".\n",sep="")
         },
         "continuation-ratio" = switch(z$direction,
           "upwards" = {
             cat("Continuation-ratio distribution (upwards).\n")
             cat("Transformation: ",z$transformation,".\n",sep="")
           },
           "downwards" = {
             cat("Continuation-ratio distribution (downwards).\n")
             cat("Transformation: ",z$transformation,".\n",sep="")
           }
           ),
         "proportional odds" = {
           cat("Proportional odds distribution.\n")
           cat("Transformation: ",z$transformation,".\n",sep="")
           cat("Link:           ",links[z$link],".\n",sep="")
         },
         "adjacent categories" = {
           cat("Adjacent categories distribution.\n")
           cat("Transformation: ",z$transformation,".\n",sep="")
         }
         )
  if(z$distribution=="multinomial")
    np <- (z$nc-1)*(z$ncv+1)
  else
    np <- z$nc-1+z$ncv
  print.matrix(rbind(-round(z$likelihood,2),if(z$individual) as.character(z$n-np) else as.character(z$events-np),as.character(np),-round(z$likelihood,2)+np,as.character(z$iterations)),
               rowlab=c("-Log likelihood ","Degrees of freedom ","Number of parameters ","AIC ","Iterations "),collab="",quote=F,right=T)
  if(z$iterlim==z$iterations)
    cat("\nWarning: Number of iterations exceeded.\n")
  if(z$info!=0&&z$iterlim!=0)
    cat("\nWarning: Exact singularity in QR decomposition.\n")
  if(z$rank!=np&&z$iterlim!=0)
    cat("\nWarning: Second derivative matrix is singular.\n")
  res <- cbind(z$coefficients,z$se)
  if(z$distribution=="binary")
    rownames(res) <- c("(Intercept)",colnames(z$tvcov$tvcov))
  else
    if(z$distribution=="multinomial"&&z$ncv!=0)
      rownames(res) <- as.vector(t(cbind(matrix(paste("(Intercept",1:(z$nc-1),")",sep=""),ncol=1),
                                         matrix(rep(colnames(z$tvcov$tvcov),(z$nc-1)),nrow=(z$nc-1),byrow=T))))
    else
      if((z$distribution=="continuation-ratio"&&z$direction=="upwards")||z$distribution=="proportional odds")
        rownames(res) <- c(paste("(Intercept",0:(z$nc-2),")",sep=""),colnames(z$tvcov$tvcov))
      else
        rownames(res) <- c(paste("(Intercept",1:(z$nc-1),")",sep=""),colnames(z$tvcov$tvcov))
  colnames(res) <- c("estimate","s.e.")
  cat("\nRegression coefficients\n")
  print.default(res,digits=digits,quote=F,print.gap=print.gap)
  if(length(z$coefficients)>1&&z$iterations>0) {
    rownames(z$correlations) <- colnames(z$correlations) <- 1:length(z$coefficients)
    cat("\nCorrelation matrix\n")
    print.default(z$correlations,digits=digits,quote=F,print.gap=print.gap)
  }
  cat("\n")
}

print.kalordinal <- function(z,digits=max(3,.Options$digits-3)) {
  nc <- length(unique(z$response$y))
  if(!is.null(z$ccov))
    nccov <- dim(z$ccov$ccov)[2]
  else
    nccov <- 0
  deppar <- (z$depend=="serial"||z$depend=="Markov")
  cat("\nCall:",deparse(z$call),sep="\n")
  cat("\n")
  if(z$optimize&&z$code>2)
    cat("Warning: no convergence - error",z$code,"\n\n")
  if(length(nobs(z))!=length(z$response$y))
    cat("Number of subjects    ",length(nobs(z)),"\n")
  cat("Number of observations",length(z$response$y),"\n")
  if(!is.null(z$response$time))
    cat("Times centred at      ",mean(z$response$time),"\n")
  if(z$depend=="independence")
    cat(" with independence\n")
  else
    cat(" with",z$depend,"dependence\n")
  cat("\n-Log likelihood   ",z$maxlike,"\n")
  cat("Degrees of freedom",z$df,"\n")
  cat("AIC               ",z$aic,"\n")
  if(z$optimize)
    cat("Iterations        ",z$iterations,"\n")
  cat("\nLocation 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(nc==2)
    cname <- "(Intercept)"
  else
    cname <- paste("(Intercept",1:(nc-1),")",sep="")
  if(inherits(z$mu,"formulafn"))
    cname <- c(cname,if(is.character(attr(z$mu,"model"))) attr(z$mu,"model") else attr(z$mu,"parameters"))
  else {
    if(nccov)
      cname <- c(cname,colnames(z$ccov$ccov))
    if(z$torder) {
      cname <- c(cname,paste("t^",1:z$torder,sep=""))
      if(length(z$interaction)>0) {
        for(i in 1:nccov)
          if(z$interaction[i]>0) {
            cname <- c(cname,paste(colnames(z$ccov$ccov)[i],".t^",1:z$interaction[i],sep=""))
          }
      }
    }
    if(!is.null(z$tvcov))
      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$init||deppar) {
    cat("\nNonlinear parameters\n")
    if(z$init)
      coef <- exp(z$coef[(z$npt-deppar):z$npt])
    else
      coef <- exp(z$coef[z$npt])
    cname <- NULL
    if(z$init)
      cname <- c(cname,"initial")
    if(deppar) {
      if(z$init)
        coef[2] <- coef[2]/(1+coef[2])
      else
        coef <- coef/(1+coef)
      cname <- c(cname,"depend")
    }
    if(z$init)
      coef.table <- cbind(z$coef[(z$npt-deppar):z$npt],z$se[(z$npt-deppar):z$npt],coef)
    else
      coef.table <- cbind(z$coef[z$npt],z$se[z$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[z$npt-deppar]))
      cat("Correlation =",tmp/(tmp+trigamma(1)),"\n")
    }
  }
  if(z$optimize&&dim(z$corr)[1]>1) {
    cat("\nCorrelation matrix\n")
    print.default(z$corr,digits=digits)
  }
}

print.moprofile <- function(z,nind=NULL,digits=max(3,.Options$digits-3),print.gap=2) {
  if(!is.null(nind)&&!z$individual) {
    warning("Option `nind' is ignored as the response corresponse to a frequency table.")
    nind <- NULL
  }
  if(!is.null(nind)) {
    if(max(nind)>z$events||min(nind)<1)
      stop("No such individual.")
    ns <- length(nind)
  }
  if(z$curve.type=="probability")
    cat("\nCategory containing highest probability:\n\n")
  else
    cat("\nCumulative probabilities:\n\n")
  tim <- unique(z$response$times)
  nobs <- nobs(z)
  ii <- covind(z)
  if(is.null(nind)) {
    if(z$individual) {
      cat("Individual data.\n")
      cat("Total number of individuals: ",z$events,"\n",sep="")
      cat("Number of observations:      ",z$n,"\n",sep="")
      cat("Categories available:        ",paste(0:(z$nc-1),collapse=" "),"\n",sep="")
      if(length(tim)>1)
        cat("Number of time points:       ",length(tim),"\n",sep="")
    }
    else {
      cat("Frequency table.\n")
      cat("Number of non-empty cells: ",z$n,"\n",sep="")
      cat("Total number of events:    ",z$events,"\n",sep="")
      cat("Categories available:      ",paste(0:(z$nc-1),collapse=" "),"\n",sep="")
      if(length(tim)>1)
        cat("Number of time points:     ",length(tim),"\n",sep="")
    }
  }
  else {
    cat("Individual data.\n")
    if(ns==1)
      cat("Individual selected:    ",nind,"\n",sep="")
    else
      cat("Individuals selected:   ",paste(nind,collapse=",",sep=""),"\n",sep="")
    cat("Number of observations: ",z$n,"\n",sep="")
    cat("Categories available:   ",paste(0:(z$nc-1),collapse=" "),"\n",sep="")
    if(length(tim)>1)
      cat("Number of time points:  ",length(tim),"\n",sep="")
  }
  cat("\n")
  if(z$curve.type=="probability") {
    if(is.null(nind)) {
      if(z$individual) {
        if(length(tim)>1) {
          res <- matrix(NA,nrow=length(nobs),ncol=length(tim),dimnames=list(1:length(nobs),paste("Time ",tim,sep="")))
          for(i in tim)
            res[ii[z$response$times==i],i] <- z$pred[z$response$times==i]
          rm(i)
        }
        else
          res <- z$pred
        print.default(res,digits=digits,quote=F,print.gap=print.gap)
        cat("\n")
      }
      else {
        if(length(tim)>1) {
          res <- NULL
          for(i in tim) {
            tmp <- matrix(NA,nrow=length(nobs),ncol=2,dimnames=list(NULL,c("Level","Freq.")))
            tmp[ii[z$response$times==i],] <- cbind(z$pred[z$response$times==i],z$response$wt[z$response$times==i])
            res <- c(res,list(tmp))
          }
          rm(i,tmp)
          names(res) <- paste("Time ",tim,sep="")
          print.listmatrix(res,row.names=1:length(nobs),digits=digits,print.gap=print.gap)
          cat("\n")
        }
        else {
          res <- cbind(z$pred,z$response$wt)
          colnames(res) <- c("Level","Freq.")
          rownames(res) <- 1:length(nobs)
          print.default(res,digits=digits,quote=F,print.gap=print.gap)
          cat("\n")
        }
      }
    }
    else {
      if(length(tim)>1) {
        res <- NULL
        for(i in tim) {
          tmp <- rep(NA,length(nind))
          for(j in 1:length(nind))
            tmp[j] <- z$pred[z$response$times==i][ii[z$response$times==i]==nind[j]]
          res <- cbind(res,tmp)
          rm(tmp,j)
        }
        rm(i)
        colnames(res) <- paste("Time ",tim,sep="")
        rownames(res) <- nind
      }
      else
        res <- z$pred[ii=nind]
      print.default(res,digits=digits,quote=F,print.gap=print.gap)
      cat("\n")
    }
  }
  else {
    if(is.null(nind)) {
      if(z$individual)
        if(length(tim)>1) {
          res <- NULL
          j <- 0
          for(i in tim) {
            j <- j+1
            res <- c(res,list(cbind(matrix(z$cpred,nrow=z$n)[z$response$times==i,,drop=F],1)))
            colnames(res[[j]]) <- c("0",apply(cbind(0,1:(z$nc-1)),1,function(x) paste(x,collapse="-")))
            rownames(res[[j]]) <- 1:length(nobs)
          }
          names(res) <- paste("Time ",tim,":",sep="")
          print.listmatrix(res,row.names=1:length(nobs),digits=digits,print.gap=print.gap)
          cat("\n")
        }
        else {
          res <- cbind(matrix(z$cpred,nrow=z$n),1)
          colnames(res) <- c("0",apply(cbind(0,1:(z$nc-1)),1,function(x) paste(x,collapse="-")))
          rownames(res) <- rep(1:length(nobs),dim(res)[1]/length(nobs))
          print.default(res,digits=digits,quote=F,print.gap=print.gap)
          cat("\n")
        }
      else
        if(length(tim)>1) {
          res <- NULL
          j <- 0
          for(i in tim) {
            j <- j+1
            res <- c(res,list(cbind(matrix(z$cpred,nrow=z$n)[z$response$times==i,],1,z$response$wt[z$response$times==i])))
            colnames(res[[j]]) <- c("0",apply(cbind(0,1:(z$nc-1)),1,function(x) paste(x,collapse="-")),"Freq.")
          }
          names(res) <- paste("Time ",tim,":",sep="")
          print.listmatrix(res,row.names=1:length(nobs),digits=digits,print.gap=print.gap)
          cat("\n")
        }
        else {
          res <- cbind(matrix(z$cpred,nrow=z$n),1,z$response$wt)
          colnames(res) <- c("0",apply(cbind(0,1:(z$nc-1)),1,function(x) paste(x,collapse="-")),"Freq.")
          rownames(res) <- rep(1:length(nobs),dim(res)[1]/length(nobs))
          print.default(res,digits=digits,quote=F,print.gap=print.gap)
          cat("\n")
        }        
    }
    else {
      if(length(tim)>1) {
        res <- NULL
        j <- 0
        for(i in tim) {
          j <- j+1
          tmp <- NULL
          for(k in nind)
            tmp <- rbind(tmp,(matrix(z$cpred,nrow=z$n)[z$response$times==i,])[ii[z$response$times==i]==k,,drop=F])
          res <- c(res,list(cbind(tmp,1)))
          colnames(res[[j]]) <- c("0",apply(cbind(0,1:(z$nc-1)),1,function(x) paste(x,collapse="-")))
          rm(tmp,k)
        }
        rm(j)
        names(res) <- paste("Time ",tim,":",sep="")
        print.listmatrix(res,row.names=nind,digits=digits,print.gap=print.gap)
        cat("\n")
      }
      else {
        res <- cbind(matrix(z$cpred,nrow=z$n)[ii==nind,,drop=F],1)
        colnames(res) <- c("0",apply(cbind(0,1:(z$nc-1)),1,function(x) paste(x,collapse="-")))
        rownames(res) <- nind
        print.default(res,digits=digits,quote=F,print.gap=print.gap)
        cat("\n")
      }
    }
  }
}

print.listmatrix <- function(x,row.names=NULL,digits=max(3,.Options$digits-3),print.gap=2) {
  res <- NULL
  for(i in 1:length(x)) {
    if(is.null(colnames(x[[i]])))
      names <- rep("",ncol(x[[i]]))
    else
      names <- colnames(x[[i]])
    dimnames(x[[i]]) <- NULL
    x[[i]] <- round(x[[i]],digits)
    mat <- array("",dim=dim(x[[i]]))
    for(j in 1:ncol(x[[i]])) {
      n <- lapply(strsplit(x[[i]][,j],"\\."),nchar)
      r <- unlist(lapply(n,length))==2
      if(any(r)) {
        len.l <- len.r <- vector(mode="numeric",length=length(n))
        len.l[r] <- unlist(n[r])[(1:length(n[r]))*2-1]
        if(any(!r))
          len.l[!r] <- unlist(n[!r])
        len.r[r] <- unlist(n[r])[(1:length(n[r]))*2]
        tot <- max(nchar(names[j]),max(len.l)+1+max(len.r))
        names[j] <- paste(paste(rep(" ",tot-nchar(names[j])),collapse="",sep=""),names[j],collapse="",sep="")
        mat[,j] <- paste(substr(as.character(x[[i]][,j]),1,len.l),".",substr(as.character(x[[i]][,j]),len.l+2,len.l+1+len.r),sep="")
        for(k in (1:length(x[[i]][,j]))[nchar(mat[,j])!=tot])
          mat[k,j] <- paste(paste(rep(" ",max(len.l)-len.l[k]),collapse="",sep=""),mat[k,j],paste(rep(0,max(len.r)-len.r[k]),collapse="",sep=""),collapse="",sep="")
        rm(n,r,len.l,len.r,tot,k)
      }
      else {
        len <- unlist(n)
        tot <- max(c(nchar(names[j]),len))
        names[j] <- paste(paste(rep(" ",tot-nchar(names[j])),collapse="",sep=""),names[j],collapse="",sep="")
        mat[,j] <- as.character(x[[i]][,j])
        for(k in (1:length(x[[i]][,j]))[nchar(mat[,j])!=tot])
          mat[k,j] <- paste(paste(rep(" ",tot-len[k]),collapse="",sep=""),mat[k,j],collapse="",sep="")
        rm(n,r,len,tot,k)
      }
    }
    mat <- rbind(names,mat)
    mat <- as.matrix(apply(mat,1,function(x) paste(x,collapse=" ")))
    dimnames(mat) <- NULL
    cname <- names(x[i])
    len <- unique(nchar(mat))
    center <- len-nchar(cname)
    if(center>0)
      cname <- paste(paste(rep(" ",floor(center/2)),collapse="",sep=""),cname,paste(rep(" ",center-floor(center/2)),collapse="",sep=""))
    else
      if(center<0)
        mat <- as.matrix(paste(paste(rep(" ",floor(abs(center)/2)),collapse="",sep=""),mat,
                               paste(rep(" ",abs(center)-floor(abs(center)/2)),collapse="",sep=""),sep=""))
    if(i!=length(x))
      mat <- as.matrix(paste(mat,paste(rep(" ",3*print.gap),collapse="",sep=""),sep=""))
    colnames(mat) <- cname
    res <- cbind(res,mat)
    rm(names,mat,j,cname,len,center)
  }
  if(is.null(row.names)||(length(row.names)+1)!=dim(res)[1])
    rownames(res) <- rep("",dim(res)[1])
  else
    rownames(res) <- c("",row.names)
  print.default(res,quote=F,print.gap=print.gap)
  invisible(res)
}
