f.auc <- function(dat, time)
{
   ## Purpose: calculates the AUC (= area under curve)
   ## ------------------------------------------------
   ## Arguments: dat  : data frame or matrix
   ##            time : vector of times
   ## ------------------------------------------------
   ## Author: Christian Keller, Date:  4 Mar 97, 15:09
   ##         Peter Holzer, Date: 17. July 2002, 09:30

   n <- length(time)
   if (n != ncol(dat))
   stop("length of argument 'time' must equal the number of columns in argument 'data'")  
   auc <- numeric(nrow(dat))
   for (i in seq(nrow(dat))) {
     datrow <- unlist(dat[i, !is.na(dat[i,])])  # unlist: needed if dat is 
                                                # a data.frame
     y <- (datrow[-length(datrow)] + datrow[-1]) / 2
     tdiff <- diff(time[!is.na(dat[i,])])       # omit time points where 
                                                # value of subject i is NA
     auc[i] <- as.vector(y %*% tdiff)
     }
     auc[is.na(dat[,1]) | is.na(dat[,n])] <- NA # if border point of subject
                                                # i is NA -> auc is NA
     return (auc)
}


f.mult2uni <- function(data.mult, var)
{
  ## Purpose: transformes multivariate data set to univariate
  ## -------------------------------------------------------------------------
  ## Arguments: data.mult : multivariate data set (data frame)
  ##                  var : time variables
  ## -------------------------------------------------------------------------
  ## Author: Christian Keller, Date: 26 Feb 97, 14:04
  ## Update: R. Frisullo, Date: 22 Juli 2002 10:18

  p <- length(var)
  n <- dim(data.mult)[1]
  data.uni <- data.mult[rep(1:n, each=p), -var, drop=FALSE]
  ymat <- data.matrix(data.mult[,var])
  names.var <- names(data.mult)[var]
  data.uni <- cbind(data.uni,
                    period=ordered(rep(names.var,n),levels=names.var),
                    y=as.vector(t(ymat)))
  row.names(data.uni) <- 1:(p*n)
  data.uni
}

## Greenhouse-Geisser and Huynh-Feldt epsilons
## -------------------------------------------
f.epsi <- function(S,p,g,n)
   {
   ## Purpose: calculates the Greenhouse-Geisser and Huynh-Feldt epsilons
   ## -------------------------------------------------------------------
   ## Arguments: S pxp covariance matrix
   ##            p dimension of observation vectors
   ##            g number of groups
   ##            n number of subjects
   ## Lit:    E.F. Vonesh + V.M. Chinchilli (1997), p.84-86
   ##         M.J. Crowder and D.J. Hand (1990), p.54-55
   ## Author: H.-R. Roth
   ## Date:   23.07.2002
   ## -------------------------------------------------------------------

   # U is a matrix of (p-1) orthonormal contrasts
   U <- t(cbind(diag(p-1),0) - outer(1:(p-1),1:p,"<") / ((p-1):1))
   a <- 1/sqrt(colSums(U^2))
   U <- U%*%diag(a)
   V <- t(U) %*% S %*% U
   e <- (sum(diag(V)))^2/sum(diag(V%*%V))/(p-1)

   GGepsilon <- e
   HFepsilon <- min(1, (n*(p-1)*e - 2) / ((p-1)* (n-g-(p-1)*e) ))
   t.output <- c(GGepsilon,HFepsilon)
   names(t.output)  <- c("G-G-epsilon", "H-F-epsilon")
   t.output
   }