# ==========================================================================
f.table <- function(x1=NULL, x2=NULL, count="count", data=NULL,
  missing.cell=NA, na.exclude=FALSE, print.it="table", 
  row=TRUE, col=TRUE, test=TRUE, test.tables=test, 
  perc.digits=0, yates=TRUE)
{
  ## Purpose:   analysis of 2-dim. contingency tables
  ## -------------------------------------------------------------------------
  ## Arguments:
  ##   usually given: x1, x2 and possibly data.
  ##        x1, x2 : variables, or
  ##            x1 : matrix (already a table)
  ##         count : frequencies of [x1,x2] (if different from 1)
  ##          data : data frame. If it is given, x1 and x2 must be columns
  ##                 of data (they default to the first two columns)
  ##  missing.cell : if count is given and some combinations do not occur,
  ##                 they are filled with the value of missing.cell
  ##    na.exclude : should NAs be shown in the table?
  ##      print.it : what should be printed directly
  ##                 (regardless of an assignment of value)
  ##   Further arguments select the components of the value
  ##
  ## Value: list with components
  ##      col, row : column and row percentages
  ##          test : chisquare test
  ##   test.tables : expected counts and pearson residuals
  ##      print.it : print the results
  ##   perc.digits : number of digits for column and row percents
  ##         yates : Yates' continuity correction (only for 2 x 2 tables) 
  ## -------------------------------------------------------------------------
  ## Author: Werner Stahel, Date: 14 Oct 96
  #on.exit(browser())
  twodim <- TRUE
  if (!is.null(data)) {
    nm <- names(data)
    if (is.character(count)&length(count)==1) {
      jc <-match(count,nm)
      if(!is.na(jc)) count <- data[[jc]]
      nm <- nm[-jc] }
    if (is.null(x1)) x1 <- nm[1]
    if (is.null(x2)) x2 <- nm[2]
    nm <- names(data)
    if (!is.character(x1)) x1 <- nm[x1]
    if (!is.character(x2)) x2 <- nm[x2]
    nm <- c(x1,x2)
    x1 <- data[[x1]]
    x2 <- data[[x2]]
  } else 
    nm <- c(paste(substitute(x1),collapse=""),
      if (!is.null(x2)) paste(substitute(x2),collapse="") else "")
# ---
  if(is.matrix(x1)) {
    tb <- x1
    nm <- names(dimnames(x1))
    if (all(nm=="")) nm <- NULL }
    else {
      if (is.null(x2)) {
        tb <- rbind(table(x1))
        twodim <- col <- test <- test.tables <- FALSE }
      else {
      # --- cross two variables
      x1 <- factor(x1)
      x2 <- factor(x2)
    if (length(x1)!=length(x2)) stop("f.table: unsuitable arguments x1, x2")
    if (is.numeric(count)) {
	nm1 <- levels(x1)
	nm2 <- levels(x2)
	tb <- matrix(missing.cell, length(nm1),length(nm2)) 
        dimnames(tb) <- list(nm1,nm2)
	tb[cbind(x1,x2)] <- count
      # --- more than one row for a cell
	if (any(duplicated(outer(nm1,nm2,paste)))) {
 	  tb <- 0*tb
	  tb[cbind(x1,x2)] <- tb[cbind(x1,x2)]+count }
    } else {
      if (na.exclude) {
        tb <- table(x1,x2)
      } else {
        tb <- table(x1,x2, exclude=NULL)
      }
     if (is.numeric(x1))
       dimnames(tb)[[1]] <- c(format(sort(unique(x1))),
                              if((!na.exclude)&&any(is.na(x1))) "NA")
     if (is.numeric(x2))
       dimnames(tb)[[2]] <- c(format(sort(unique(x2))), 
                              if((!na.exclude)&&any(is.na(x2))) "NA")
    }
#    }
#    else  
  } }
  tb1 <- cbind(tb, total = sr <- apply(tb, 1, sum))
  sr <- c(sr, n <- sum(sr))
  if(twodim) tb1 <- rbind(tb1, total = sc <- apply(tb1, 2, sum))
  attr(tb1, "margin.names") <- nm
  names(dimnames(tb1)) <- names(dimnames(tb))
  tbc <- if(col) 100 * sweep(tb1, 2, sc, "/") else NULL
  tbr <- if(row) 100 * sweep(tb1, 1, sr, "/") else NULL
  ch <- if(test) chisq.test(tb, correct=yates) else NULL
  ex <- pr <- NULL
  if(test.tables){
    ex <- outer(last(sr,-1),last(sc,-1)/n)
    pr <- (tb-ex)/sqrt(ex)
  }
  rr <- list(table=tb1, margin.names=nm, colperc=tbc, rowperc=tbr,
             fit=ex, resid.pearson=pr, chisq=ch)
  class(rr) <- c("f.table","list")
  if (!(is.logical(print.it)&&!print.it)) {
    prt <- list(margin.names=nm,
      table=tb1[,], colperc=if(col) round(tbc[,],perc.digits),
      rowperc=if(row) round(tbr[,],perc.digits),
      fit=if(test.tables) round(ex[,],2),
      resid.pearson=if(test.tables) round(pr[,],2), chisq=if(test) ch)
    prt <- prt[!sapply(prt,is.null)]
    if (is.logical(print.it)) 
      if(print.it) print(prt) else return(rr)
    else if (is.character(print.it)) print(prt[print.it]) }
  #if (!Browse) on.exit()
  invisible(rr)
}
# ==========================================================================



last <- function(x, length.out = 1, na.rm = FALSE)
{
  ## Purpose: last element(s) of a vector
  ## Author: Werner Stahel, Date:  Tue Jan 21 17:29:42 1992
  ## ----------------------------------------------------------------
  ## Arguments:
  ##   x:          vector
  ##   length.out: if positive, return the  length.out last elements of x,
  ##               if negative, the last  length.out  elements are dropped
  ## ----------------------------------------------------------------
  if (na.rm) 
    x <- x[!is.na(x)]
  n <- length(x)
  x[sign(length.out)*(n-abs(length.out)+1):n]
}