R: fun.R

From MathWiki

"

This collection of functions is now available as the spida package at http://r-forge.r-project.org/. 'spida' and a companion package for 3-d plotting called 'p3d' can be installed with:

> install.packages(c('spida','p3d'),repos="http://r-forge.r-project.org")

and used with:

> library(spida)
> library(p3d)

Both are in development.

If you are interested in more information on developing packages in R-Forge, consult this document (http://cran.r-project.org/doc/contrib/Graves+DoraiRaj-RPackageDevelopment.pdf).

fun.R: A collection of utility functions and functions for multilevel modeling kept by Georges Monette. These file can be downloaded or sourced in R from http://www.math.yorku.ca/~georges/R/fun.R

NOTE 1: Please send any suggestions or problems to georges@yorku.ca.

NOTE 2: There is a copy of this file at //wiki.math.yorku.ca/index.php/R:fun.R which can be edited. Changes will be incorporated regularly into the the downloadable file.

NOTE 3: THIS FILE USES QUOTES AND TAGS SO IT CAN BE RENDERED AS A WIKI FILE AS WELL AS R SOURCE

BE SURE TO LEAVE WIKI TEXT IN DOUBLE QUOTES AND AVOID DOUBLE QUOTES WITHIN WIKI TEXT

Last uploaded to http://www.math.yorku.ca/~georges/R/fun.R : May 28, 2006

Decribe modifications here
Modified: --Georges Monette 14:30, 28 May 2006 (EST)
added functions from funRDC
Modified: --Georges Monette 05:56, 22 Feb 2006 (EST)
added Lmat
Modified: --Georges Monette 09:45, 2 Nov 2005 (EST)
added class 'cat' to coursefun
defined print.cat to print with cat
Modified:
plot3d, identify3d by John Fox so they work with matrices, data.frames or variable arguments
ell - corrected error when using radius


Table of contents

General description

"
##
##
##  Some R functions for PSYC 6140 and MATH 6630
##  2005-2006
##
## Last update: October 27, 2005
## Summary:
##
##
## Tables:
##
##    atotal: border an array with sums
##    abind : glue two arrays together on a selected dimension
##
## General Linear Hypothesis
##
#     glh   : glh( fit, L )
##    Lmat  : generates L matrix for anova in library(lmer) or lht in library(car)
##
## Graphics:
##
##    td    : easy front end to trellis.device and related functions
##    xqplot: extended quantile plots
##
## 3D graphics by John Fox:
##
##    scatter3d
##    identify3d
##    ellipsoid
##    plot3d     - wrapper for scatter3d
##
## Inference
##    cell   - a modified version of car::confidence.ellipse.lm that
##             creates a confidence ellipse to plot with lines(...)
##    dell   - data ellipse
##
##            
# Note that some functions that were transferred and improved in coursefun.R
# have been 'disabled' by appending '.rdc' to function name
#
#  Splus functions written in RDC
#  2005:
#  April 19       Modified for R
#  May 3          copied atotal, abind
#                 wrote acond
#  May 10         new function: cap1 to capitalize initial letters and turn underscores to blanks
#  May 12         new function adapted from gm: td 
#  May 13         new function: write.sas to write data frame to SAS without truncating variable names
#  May 16         added Poly and centering to splines
#  June 13        added constant to check if values are constant within levels of id
#                 added varLevel( data.frame, ~lev1/lev2) to report level of each variable
#                 added Lmat to generate L matrix with 1's for effects containing a string
#                 modified anova.lme to use min DFs in 'denominator' DFs
#  August 3       modified Lag to accept 'at' argument
#  August 5       changed 'arep' to 'apct' in order to parallel atotal and acond
#                 changed acond to aprop
#  August 15      getFix, glh and and print.glh, Q, Vcov, Vcor
#  August 25      Contrasts
#  2006
#  May 19         fill and capply




fun <- "
fun.R                         Last update: May 28, 2006
An R script containing functions and some datasets for the courses
PSYC6140 and MATH6630 in 2005-06 and MATH6643 in Summer 2006

Help on some of following functions is available by typing the name
of the function. This text is available by typing
   > fun
   
A current version of this file can be sourced or downloaded from
   http://www.math.yorku.ca/~georges/R/fun.R
   
A copy is kept at:
   http://wiki.math.yorku.ca/R:_fun.R

Please make corrections, changes and additions to the version on the wiki.
They will be periodically transferred to the downloadable version.

Contact: Georges Monette <georges@yorku.ca> 

Partial list of functions:

Tables:
    atotal: border an array with sums
    abind : glue two arrays together on a selected dimension

Graphics:

    td    : easy front end to trellis.device and related functions
    xqplot: extended quantile plots

3D graphics by John Fox:

    scatter3d
    identify3d
    ellipsoid

Inference
    glh    - General Linear Hypothesis of fixed effects using a Wald test
             - should work on lm, lme, nlme and lmer objects
             - glh( fit, L ) where 'L' is a hypothesis matrix performs a Wald test
             - glh( fit, 'pat' ) where 'pat' is a regular expression (see ?regex)
               uses an L matrix that tests all fixed effects whose names match
               'pat'.  For example 
                  > glh( fit, ':.*:' )
               will test all 2- and higher-way interactions together.
             - glh( fit, c(2,3,5) )  will test the 2nd, 3rd and 5th coefficients.
    Lmat   - Lmat ( fit, 'pat' ) constructs a hypothsis matrix.
    Lmu
    Ldiff
               
    cell   - a modified version of car::confidence.ellipse.lm that
             can add to a plot

Graphics for linear algebra

    vplot  - plots column vectors adding to current plot
    vell   - ellipse as a 2 x n matrix
    vbox   - box around unit circle
    orthog - 2 x 2 rotation
    orthog.prog 2 x 2 matrix of orthog projection

File utilities

    Read.spss  : uses spss.get in library (Hmisc)
    Read.dta   : uses read.dta in library (foreign) with trimming and conversion 
                 to factors
    trim       : removes trailing blanks in character variables and factor levels
                 - useful to post-process files from SPSS read with read.dta


Multilevel utilities:

    varLevel( dd, fmla )   shows the level of each variable in data.frame dd
           - e.g. varLevel( dd, ~Board/School/Class/Subject ) or just
           - varLevel( dd, ~ID ) to identify inner and outer variables
           
    capply ( x, by, FUN, ...) creates 'contextual' variables 
           - like tapply but replicating the original shape

    fill ( x, by, FUN , ... )   uses capply to create contextual variable by selecting
             a value. Works intelligently with factors and 'Date's.
           - the default FUN selects a unique value and can be used when
             a long data frame has been created with values for macro variables
             only in the first position.
"
 

class(fun) <- "cat"
print.cat <- function(x,...) {
     cat(x,...)
     invisible(x)
}

print(fun)

"

General Linear Hypothesis

"
##
##
##   Functions to perform a GLH on lm, lme or lmer models
##   August 13, 2005
##
##
##
##   Lmat: generate a hypothesis matrix based on a pattern
##
##   glh
##   Lmat
##   Ldiff
##   getFix
##
##   print.glh
##

Lmat <- function( fit, pat, ... ) {
     if ( is.character(fit)) {
          x <- fit
          fit <- pat
          pat <- x
     }
     # generates L matrix for anova(fit, L = ...) with all terms
     # containing a pattern
     gf <- getFix(fit)
     fe <- gf$fixed
     L.indices <- grep(pat, names(fe))
     ret <- diag(length(fe))[ L.indices,]
     rownames(ret) <- names(fe) [L.indices]
     ret
}


Ldiff <- function( fit, pat, levnames = c(substring(rownames(L),cut+1),reflevel),
         reflevel ="<ref>", cut=nchar(pat)) {
      L <- Lmat(fit, pat)
      nam <- rownames(L)
      n <- nrow(L)
      zm <- matrix(1:n,nrow=n,ncol=n)
      plus <- zm[ col(zm) < row(zm)]
      minus <- rep(1:(n-1), (n-1):1)
      Lp <- L[plus,]
      Lm <- L[minus,]
      Lret <- rbind( L, Lp - Lm)
      rn <- paste( levnames[ c(1:n,plus)+1], minus[ c(rep(0,n),minus) + 1], sep = " - ")
      rownames(Lret) <- rn
      Lret
}

getFix <- function(fit,...) UseMethod("getFix")

getFix.lm <- function(fit,...) {
       ss <- summary(fit)
       ret <- list()
       ret$fixed <- coef(fit)
       ret$vcov <- ss$sigma^2 * ss$cov.unscaled
       ret$df <- rep(ss$df[2], length(ret$fixed))
       ret
}

getFix.lme <- function(fit,...) {
       require(nlme)
       ret <- list()
       ret$fixed <- nlme::fixef(fit)
       ret$vcov <- fit$varFix
       ret$df <- fit$fixDF$X
       ret
}

getFix.lmer <- function(fit,...) {
       ret <- list()
       ret$fixed <- fixef(fit)
       ret$vcov <- vcov(fit)
       ret$df <- Matrix:::getFixDF(fit)
       ret
}

Vcov <- function(fit) {
     getFix(fit)$vcov
}

Vcor <- function(fit) {
     vc <- cov2cor(getFix(fit)$vcov)
     svds <- svd(vc)$d
     attribute(vc,'conditionNumber') <- svds[1]/svds[length(svds)]
     vc
}


# fit <- lmer( Yield ~ Location *   Family  + (1|Block), data = Genetics)
# getFix(fit)
# fit <- lme( Yield ~ Location *   Family  , data = Genetics, random = ~1|Block)

# L <- rbind( c(1,1,1,1), c(0,1,0,0), c(1,0,1,1))


tfun <- function( x) {
      ret <- NULL
      if ( is.character(x)) ret <- "it's character"
      else { 
         if ( is.numeric(x)) {
            if ( is.null(dim(x))) {
                if ( length(x) != 4 ) ret <- diag(4)[x,]
                else ret <- rbind( x )
            }
         }
      }
      ret 
}

glh <- function(fit, Llist,clevel=0.95,debug=F) {
help <- "
glh:  General Linear Hypothesis with Wald test
      for lm, lme, nlme and lmer objects.

      Can be extended to other objects (e.g.) 'glm' by writing 'getFix.glm'
      
      Usage:
         glh ( fit, L ) where L is hypothesis matrix
         glh ( fit, 'pat' ) where 'pat' a regular expression (see ?regex)
           used to match names of coefficients of fixed effects.
           e.g. glh( fit, ':.*:') tests all 2nd and higher order interactions.
         glh ( fit, c(2,5,6)) to test 2nd, 5th and 6th coefficients. 
"

           unique.rownames <- function(x) {
                ret <- c(tapply(1:length(x), x , function(xx) {
                    if ( length(xx) == 1) ""
                    else 1:length(xx)
                })) [tapply(1:length(x),x)]
                ret <- paste(x,ret,sep="")
                ret
           }
           if(!is.list(Llist)) Llist <- list(Llist)
           ret <- list()
           fix <- getFix(fit)
           if(debug) print(fix)
           beta <- fix$fixed
           vc <- fix$vcov


           dfs <- fix$df
           for (ii in 1:length(Llist)) {
               ret[[ii]] <- list()
               Larg <- Llist[[ii]]
               L <- NULL
               if ( is.character(Larg)) {
                  L <- Lmat(fit,Larg)
               } else { 
                  if ( is.numeric(Larg)) {
                     if ( is.null(dim(Larg))) {
                        if(debug) cat("is.null(dim(Larg))\n")
                        if ( (length(Larg) < length(beta)) && all(Larg>0) ) {
                              L <- diag(length(beta))[Larg,]
                              dimnames(L) <- list( names(beta)[Larg], names(beta))
                        } else ret <- rbind( Larg )
                     } 
                     else L <- Larg
                  }
               }
               if (debug) {
                  cat("Larg:", Larg, "\nL:\n")
                  print(L)
                  cat("\n")
               }
               ## Anova
               qqr <- qr(t(L))
               L.rank <- qqr$rank
               L.full <- t(qr.Q(qqr))[ 1:L.rank,,drop=F]
               if (debug) {
                  cat("\nqr(L)::::::::::\n")
                  print(qqr)
                  cat("\nL.full:::::::::\n")
                  print(dim(L.full))
                  print(L.full)
                  cat("\nvc:::::::::::\n")
                  print(dim(vc))
                  print(vc)

               }

               if (debug) print(L.full)
               vv <-  L.full %*% vc %*% t(L.full)
               eta.hat <- L.full %*% beta
               Fstat <- (t(eta.hat) %*% qr.solve(vv,eta.hat,tol=1e-10)) / L.rank
               included.effects <- apply(L,2,function(x) sum(abs(x))) != 0
               denDF <- min( dfs[included.effects])
               numDF <- L.rank
               ret[[ii]]$anova <- list(numDF = numDF, denDF = denDF,
                               "F-value" = Fstat,
                               "p-value" = pf(Fstat, numDF, denDF, lower.tail = F))
               ## Estimate
               etahat <- L %*% beta
               etavar <- L %*% vc %*% t(L)
               etasd <- sqrt( diag( etavar ))

               denDF <- apply( L , 1 , function(x,dfs) min( dfs[x!=0]), dfs = dfs)


               aod <- cbind( Estimate=c(etahat),
                   Std.Error = etasd,
                   DF = denDF,
                   "t-value" = c(etahat/etasd),
                   "p-value" = 2*pt(abs(etahat/etasd), denDF, lower.tail =F))
                colnames(aod)[ncol(aod)] <- 'p-value'
             if (debug ) {
                cat("\n\naod::::::::::::::::::\n")
                print(aod)
             }
             if ( !is.null(clevel) ) {
                #print(aod)
                #print(aod[,'DF'])
                #print(aod[,'etasd'])
                 hw <- qt(1 - (1-clevel)/2, aod[,'DF']) * aod[,'Std.Error']
                 #print(hw)
                 aod <- cbind( aod, LL = aod[,"Estimate"] - hw, UL = aod[,"Estimate"] + hw)
                 #print(aod)
                 if (debug ) {
                    cat("\ncolnames::::::::::\n")
                    print(colnames(aod))
                 }
                 labs <- paste(c("Lower","Upper"), format(clevel))
                 colnames(aod) [ ncol(aod) + c(-1,0)] <- labs
             }
             aod <- as.data.frame(aod)

             rownames(aod) <- rownames(as.data.frame(L))
             ret[[ii]]$estimate <- aod
        }
        names(ret) <- names(Llist)
        attr(ret,"class") <- "glh"
        ret
}

print.glh <- function(x,round = 6, pround = 5,...) {
    pformat <- function(x, digits = pround) {
        x <- format(xx <- round(x,digits))
        x[ as.double(xx) == 0 ] <- paste(c("<.",rep('0',digits-1),'1'),collapse="")
        x
    }
    rnd <- function(x,digits) {
        if (is.numeric(x)) x <- round(x,digits=digits)
        format(x)
    }
             for( ii in 1:length(x)) {
                  nn <- names(x)[ii]
                  tt <- x[[ii]]
                  ta <- tt$anova
                  #cat("\n\n",nn,"\n")
                  #cat("Overall test:\n")
                  ta[["p-value"]] <- pformat(ta[["p-value"]])
                  print(as.data.frame(ta,row.names=nn))
                  te <- tt$estimate
                  #cat("Estimates:\n")
                  te[,'p-value'] <- pformat( te[,'p-value'])
                  if ( !is.null(round)) {
                     for ( ii in 1:length(te)) {
                         te[[ii]] <- rnd(te[[ii]],digits=round)
                     }
                  }
                  print(as.data.frame(te),digits=round,...)
             }
        invisible(x)
}


"

Trellis graphics

"




###
###  td
###



td <- function(
        new = F,
	col=c(3,5,4,6,7,8,2), lty=1:7, lwd=1,
	pch = 1:7, cex = 0.8, font = 1,
	len = 7,
	long = F,
        record = T,
        basecol = NULL,
	colsets = c('plot.symbol','plot.line','dot.symbol',
			'dot.line','cloud.3d','box.dot'),...) {

help <- "
td                   coursefun.R     for PSYC 6140/MATH 6630 05/06

Easier way to call 'trellis.device'

Description:

     'td' calls 'trellis.device' and sets graphical parameters for
     'superpose.line' and 'superpose.symbol'. 'td' also initializes
     a new trellis device with a white background by default.

Usage:

     td( new = F,
         col=c(3,5,4,6,7,8,2),
         lty=1:7,
         lwd=1,
	 pch = 1:7, cex = 0.8, font = 1,
	 len = 7,
	 long = F,
         record = T,
         basecol = NULL,
	 colsets = c('plot.symbol','plot.line','dot.symbol',
			'dot.line','cloud.3d','box.dot'),...)

Arguments:

     new : open a new device with: 'td(T)'.

     col, lty, lwd, pch, cex: graphical parameters for superpose.line and superpose.symbol

     len : extend the length of graphical parameters by recycling

     long: if TRUE generate a default combination of col, lty and pch with length 42

     record : initiate the history mechanism for the graphical device

     colsets: additional graphical parameter lists to be modified

Details:

     By using col and lty/pch with lengths that are relatively prime and
     by using the len argument, one can generate unique combinations,
     e.g. for len = 42 with col of length 6 and pch of length 7

Value:

References:

Contributed by:  G. Monette  2005-10-10

Modifications:

"

        # Modified for R: Oct. 10, 2004
	#
	# reset superpose.symbol and superpose.line so they have consistent length
	# equal to the max of input parameters:
	#			sps: cex, col, font, pch
	#			spl: lty, col, lwd
	# or to len
	#  This allows distinctive line styles, for example for 12 objects by using
	#  good lty's:     1,4,8,6,3,5,7
	#  and good col's: 1,8,6,3,4,5,2
	#

  require(lattice)
  if ( long ) {
    col <- c(3,5,4,6,8,2)   # drop yellow
    len <- 42 # generates 42 unique combinations of pch/lty and col
  }
  if(new) trellis.device(theme = col.whitebg, record = record, new = new, ...)
                                        # NOTE: fixed panel.superpose so lty argument
                                        # is passed to points for type = 'b'
  len <- max(len,length(col),length(lty),length(lwd),length(pch),length(cex),length(font))
  spl <- trellis.par.get("superpose.line")
  spl$lty <- rep(lty,length=len)
  spl$col <- rep(col,length=len)
  spl$lwd <- rep(lwd,length=len)
  trellis.par.set("superpose.line",spl)
  sps <- trellis.par.get("superpose.symbol")
  sps$pch <- rep(pch, length=len)
  sps$col <- rep(col, length=len)
  sps$cex <- rep(cex, length=len)
  sps$font <- rep(font, length=len)

  trellis.par.set("superpose.symbol",sps)
  list(superpose.symbol = sps, superpose.line = spl)
  if ( !is.null(basecol)) {
    for ( ii in colsets ) {
      tt <- trellis.par.get(ii)
      tt$col <- basecol
      trellis.par.set(ii,tt)
    }
  }
  invisible(attr(.Device,"trellis.settings"))
}


"

Quantile plots

"



###
###  xqplot
###


xqplot <- function(x, labels = dimnames(x)[[2]], ..., ask = F,
                   ptype = "quantile", mfrow = findmfrow ( ncol(x)),
                   mcex = 0.8, maxlab = 12 , debug = F,
                   mar = c(5,2,3,1),
                   text.cex.factor = 1 ,
                   left.labs = F,
                   maxvarnamelength = 20)
{
help <- "
xqplot                   coursefun.R     for PSYC 6140/MATH 6630 05/06

Extended Quantile Plots and Barcharts

Description:

     'xqplot' produces uniform quantile plots of numeric variables and
     barcharts of factor variables.  The display is tuned to provide a
     quick view of a data frame at a glance.

Usage:
     xqplot( x, mcex = 0.8, mar = c(5,2,3,1), text.cex.factor = 1,
             mfrow,
             maxlabs = 12,
             left.labs = F)

Arguments:

     x    : a data frame or list of variables to plot

     mcex : character expansion factor for marginal text

     mar  : size of margins

     text.cex.factor : character expansion factor for barchart labels

     left.labs : determines placement of barchart labels

     maxlab : maximum number of categories to label in barcharts

     mfrow :  number of rows and columns per page. If missing, an attempt
            is made to choose a reasonable number.

     maxvarnamelength : maximum length of variable name without splitting
            on two lines.
Details:

Value:

Bugs:

     'mfrow' should take the total number of variables into account if they will
     fill more than one page so the last page is close to being full.

     The current version of the function could be made much simpler and
     more transparent. Some code is redundant.

References:

Contributed by:  G. Monette  2005-10-10

Modifications:

"
  ## Adapted from myplot.data.frame for R by G. Monette, Oct. 25, 2004
  ##    maxlab is maximum number of labels
  if (! is.list(x)) x <- as.data.frame(x) 
  left.labs <- rep( left.labs, length = length(x))
  findmfrow <- function( x ) {
	if ( x > 9) c(3,4)
	else cbind( '1'=c(1,1),'2'=c(1,2),'3'=c(2,2),
                   '4'=c(2,2),'5'=c(2,3),'6'=c(2,3),
                   '7'=c(3,3), '8'=c(3,3), '9'=c(3,3)) [, x]
      }

  opt <- par( mfrow = mfrow, mar = mar + 0.1 )
  on.exit(par(opt))

  iscat <- function( x ) is.factor(x) || is.character(x)

  compute.cex <- function( x ) {
    ll <- length(x)
    cex <- 2 * ifelse( ll < 5, 2,
                      ifelse( ll < 10, 1,
                             ifelse( ll < 20, .7, .5)))/mfrow[1]
  }
  for ( ii in 1: dim(x)[2]) {
    vv <- x[[ii]]
    nam <- labels[[ii]]
    Nmiss <- sum(is.na(vv))
    N <- length(vv)
    if ( iscat(vv) ){
      tt <- table(vv)

      xlab <- paste("N =", N )
      if ( Nmiss > 0 ) {
        tt <- c( "<NA>" = sum(is.na(vv)), tt)
        xlab <- paste(xlab, "  Nmiss =", Nmiss)
      }
      ll <- names(tt)
      nn <- length(ll)
      if ( left.labs[ii] ) barplot( tt, horiz = T,
                                   xlab = xlab,
                                   cex.names = text.cex.factor * compute.cex(nn) )
      else {
        zm <- barplot( tt, names = rep("",nn), horiz = T, xlab = xlab)
        ## If nn > maxlab drop labels for smaller frequencies
        sel <- rep( T, length(tt))
        tt.sorted <- rev(sort(tt))
        if ( nn > maxlab ) sel <- tt > tt.sorted[maxlab]
        if (debug) {
          print(sel)
          print(nam)
          print(tt)
          print(tt.sorted)
          print(maxlab)
          print(tt.sorted[maxlab])
          print(sel)
          print(zm[sel])
          print(rep(max(tt),nn)[sel])
          print( ll[sel])
        }
        if ( any(sel) ) text( rep( max( tt ), nn)[sel]  ,
                             zm[sel], ll[sel], adj = 1, cex = text.cex.factor * compute.cex( nn ))
      }
    }
    else {
      vv <- vv[!is.na(vv)]
      xxvar <- ppoints(length(vv))
      xlab <- paste("Fraction of",N-Nmiss,"(NA:",Nmiss,")")
      if ( pmatch( ptype, 'normal', nomatch = 0) == 1 ) {
        xxvar <- qnorm( ppoints(length(xxvar)) )
        xlab <- paste("Normal quantile for",N-Nmiss,"obs. (NA:",Nmiss,")")
      }

      ## Plot continuous
      if ( N == Nmiss ) {
        xxvar <- 1
        vv <- 1
        plot( xxvar, vv, xlab = xlab, ylab="", type = 'n')
        text( 1, 1, "NA")
      }
      else {
        plot(xxvar, sort(vv), xlab = xlab, ylab = "Data", ...)
        xm <- mean(vv)
        xs <- sqrt(var(vv))
        abline( h= xm,lty=1)
        abline( h= c(xm-xs,xm+xs),lty=2)
      }
    }
    ## Titles for all plots
    vlab <- labels[ii]
    line.offset <- 1.0
    if ( nchar( vlab ) > maxvarnamelength) {
        vlab <- paste( substring(vlab,1,maxvarnamelength), "\n",substring(vlab, maxvarnamelength + 1))
        line.offset <- 0.2
    }
    mtext(vlab, 3, line.offset , cex = mcex)
  }
  invisible(0)
}



"

Contingency tables

"




### atotal

atotal <- function(arr, FUN = sum, ...) {
help <- "
atotal                coursefun.R     for PSYC 6140/MATH 6630 05/06

Adds border of sums to an array

Description:

     'atotal' adds by default a border of sums to an array.
     The funtion FUN is used instead of 'sum' is supplied. Additional
     arguments to FUN can also be given.

Usage:

     atotal( arr , FUN = sum, ...)

Arguments:

     arr: array, matrix or vector

     FUN: function to be applied to cross sections of arr

     ...: additional arguments to FUN

Details:

Value:

     An array with dimension dim(arr) + 1

References:

Contributed by:  G. Monette  2005-10-10

Modifications:

"
	d <- dim(arr)
	if (length(d) == 1) {
		arr <- c(arr)
		d <- dim(arr)
	}
	if(is.character(FUN))
		FUN <- get(FUN, mode = "function")
	else if(mode(FUN) != "function") {
		farg <- substitute(FUN)
		if(mode(farg) == "name")
			FUN <- get(farg, mode = "function")
		else stop(paste("\"", farg, "\" is not a function", sep = ""))
	}

	if (is.null(d)) return (c( arr, Total=FUN(arr)))
	n <- length(d)
	ret <- arr
	ind <- 1:n
	for ( i in n:1) {
		new <- apply(ret,ind[-i],FUN)
		ret <- abind( ret, new, i, "Total")
	}
	ret
}






###
###  abind
###


abind <- function(arr1,arr2,d,facename="") {

help <- "
abind                coursefun.R     for PSYC 6140/MATH 6630 05/06

Binds two conformable arrays along a dimension

Description:

     'abind' binds two conformable arrays, arr1 and arr2 along
     the dimension d.

Usage:

     abind( arr1, arr2 , d ,  facename )

Arguments:

     arr1, arr2: arrays, matrices or vectors

       d: dimension along which arr1 and arr2 are joined

     ...: Name of extended index if arr2 has one fewer dimensions than arr1

Details:

      dim( arr1 ) and dim( arr2 ) must be equal except in the dth dimension.
      If the length of dim( arr2 ) is 1 less than that of dim( arr1 ), then
      'arr2' is treated as if it had dropped the dth dimension with size 1.

Value:

     The returned value, ret, is an array with dimension dim(arr1)
     except for the dth dimension where dim(ret)[d] == dim(arr1)[d] + dim(arr2)[d].
     If length(dim(arr2)) == length(dim(arr1)) - 1, then arr2 is treated as if
     it dropped a dimension of size 1 in the dth dimension. 'facename' is used
     as the name of the dimnames list for this dimension.

Notes:

     abind is used by atotal

References:

Contributed by:  G. Monette  2005-10-10

Modifications:

"

	d1 <- dim(arr1)
	n1 <- length(d1)
	dn1 <- dimnames(arr1)
	d2 <- dim(arr2)
	n2 <- length(d2)
	dn2 <- dimnames(arr2)

	arenull <- is.null(dn1) & is.null(dn2)
	if (is.null(dn1)){
		dn1 <- lapply( as.list(d1), function(x) seq(1,x))
		dimnames(arr1) <- dn1
	}

	if ( n1 != n2 ) {
		d2 <- d1
		d2[d] <- 1
		dn2 <- dn1
		dn2[[d]] <- facename
		dimnames(arr2) <- NULL
		dim(arr2) <- d2
		dimnames(arr2) <- dn2
		n2 <- n1
	}
	if (is.null(dn2)){
		dn2 <- lapply( as.list(d2), function(x) seq(1,x))
		dimnames(arr2) <- dn2
	}

	# check input for consistency
	# ... later
	#

	perm <- 1:n1
	perm[c(d,n1)] <- c(n1,d)	# perm is an involution

	#
	# permute arr1
	#

	arr.p1 <- aperm(arr1,perm)

	#
	# permute arr2 if length of dimension same as arr1
	#

	arr.p2 <- aperm(arr2,perm)
	dret <- d1[perm]
	dret[n1] <- dret[n1] + d2[d]
	dnret <- dn1
	dnret[[d]] <- c(dnret[[d]],dn2[[d]])

	ret <- c(arr.p1, arr.p2)
	dim(ret) <-  dret

	#
	# permute response back
	#

	ret <- aperm(ret, perm)

	dimnames(ret) <- dnret
	ret
}





###
###   tab
###


tab <- function(...) {
help <- "
abind                fun.R     for PSYC 6140/MATH 6630 05/06

Cross Tabulation and Table Creation Including Missing Values

Description:

     'tab' does the same thing as 'table' except that it includes
     missing values for factors.  The argument 'exclude = NULL' to 'table'
     results in the inclusion of missing values for numeric variable but
     excludes missing values for factors. 'tab' is intended to remedy this
     deficiency of 'table'.

Usage:

     tab(...)

Arguments:
     ...: objects which can be interpreted as factors (including
          character strings), or a list (or data frame) whose
          components can be so interpreted.

Details:

Value:

Notes:

References:

Bugs:
      Does not use argument name as a dimension name, in contrast with 'table'.

Contributed by:  G. Monette  2005-10-10

Modifications:

"
  args <- list(...)
  if( is.list(args[[1]])) args <- args[[1]]
  # for ( ii in 1:length(a)) if ( is.factor( a[[ii]])) a[[ii]] <- factor(a[[ii]],exclude = NULL)
  for ( ii in 1:length(args))  args[[ii]] <- factor(args[[ii]],exclude = NULL)
  do.call("table", args)
}


"

3d graphics from John Fox

"


###
###    3d functions by J. Fox
###    with a few modifications by GM
###



# last modified 17 Sept 05 by J. Fox

ellipsoid <- function(center=c(0, 0, 0), radius=1, shape=diag(3), n=30){
# adapted from the shapes3d demo in the rgl package
  degvec <- seq(0, 2*pi, length=n)
  ecoord2 <- function(p) c(cos(p[1])*sin(p[2]), sin(p[1])*sin(p[2]), cos(p[2]))
  v <- t(apply(expand.grid(degvec,degvec), 1, ecoord2))
  v <- center + radius * t(v %*% chol(shape))
  v <- rbind(v, rep(1,ncol(v)))
  e <- expand.grid(1:(n-1), 1:n)
  i1 <- apply(e, 1, function(z) z[1] + n*(z[2] - 1))
  i2 <- i1 + 1
  i3 <- (i1 + n - 1) %% n^2 + 1
  i4 <- (i2 + n - 1) %% n^2 + 1
  i <- rbind(i1, i2, i4, i3)
  qmesh3d(v, i)
  }

scatter3d <- function(x, y, z, xlab=deparse(substitute(x)), ylab=deparse(substitute(y)),
    zlab=deparse(substitute(z)), revolutions=0, bg.col=c("white", "black"),
    axis.col=if (bg.col == "white") "black" else "white",
    surface.col=c("blue", "green", "orange", "magenta", "cyan", "red", "yellow", "gray"),
    neg.res.col="red", pos.res.col="green", point.col="yellow",
    text.col=axis.col, grid.col=if (bg.col == "white") "black" else "gray",
    fogtype=c("exp2", "linear", "exp", "none"),
    residuals=(length(fit) == 1), surface=TRUE, fill=TRUE, grid=TRUE, grid.lines=26,
    df.smooth=NULL, df.additive=NULL,
    sphere.size=1, threshold=0.01, speed=1, fov=60,
    fit="linear", groups=NULL, parallel=TRUE, ellipsoid=FALSE, level=0.5,
    model.summary=FALSE){
    require(rgl)
    require(mgcv)
    summaries <- list()
    if ((!is.null(groups)) && (nlevels(groups) > length(surface.col)))
        stop(sprintf(
            gettextRcmdr("Number of groups (%d) exceeds number of colors (%d)."),
            nlevels(groups), length(surface.col)))
    if ((!is.null(groups)) && (!is.factor(groups)))
        stop(gettextRcmdr("groups variable must be a factor."))
    bg.col <- match.arg(bg.col)
    fogtype <- match.arg(fogtype)
    if ((length(fit) > 1) && residuals && surface)
        stop(gettextRcmdr("cannot plot both multiple surfaces and residuals"))
    xlab  # cause these arguments to be evaluated
    ylab
    zlab
    rgl.clear()
    rgl.viewpoint(fov=fov)
    rgl.bg(col=bg.col, fogtype=fogtype)
    valid <- if (is.null(groups)) complete.cases(x, y, z)
        else complete.cases(x, y, z, groups)
    x <- x[valid]
    y <- y[valid]
    z <- z[valid]
    if (!is.null(groups)) groups <- groups[valid]
    x <- (x - min(x))/(max(x) - min(x))
    y <- (y - min(y))/(max(y) - min(y))
    z <- (z - min(z))/(max(z) - min(z))
    size <- sphere.size*((100/length(x))^(1/3))*0.015
    if (is.null(groups)){
        if (size > threshold) rgl.spheres(x, y, z, color=point.col, radius=size)
            else rgl.points(x, y, z, color=point.col)
            }
    else {
        if (size > threshold) rgl.spheres(x, y, z,
            color=surface.col[as.numeric(groups)], radius=size)
        else rgl.points(x, y, z, color=surface.col[as.numeric(groups)])
        }
    rgl.lines(c(0,1), c(0,0), c(0,0), color=axis.col)
    rgl.lines(c(0,0), c(0,1), c(0,0), color=axis.col)
    rgl.lines(c(0,0), c(0,0), c(0,1), color=axis.col)
    rgl.texts(1, 0, 0, xlab, adj=1, color=text.col)
    rgl.texts(0, 1, 0, ylab, adj=1, color=text.col)
    rgl.texts(0, 0, 1, zlab, adj=1, color=text.col)
    if (ellipsoid) {
        dfn <- 3
        if (is.null(groups)){
            dfd <- length(x) - 1
            radius <- sqrt(dfn * qf(level, dfn, dfd))
            ellips <- ellipsoid(center=c(mean(x), mean(y), mean(z)),
                shape=cov(cbind(x,y,z)), radius=radius)
            if (fill) shade3d(ellips, col=surface.col[1], alpha=0.1, lit=FALSE)
            if (grid) wire3d(ellips, col=surface.col[1], lit=FALSE)
            }
        else{
            levs <- levels(groups)
            for (j in 1:length(levs)){
                group <- levs[j]
                select.obs <- groups == group
                xx <- x[select.obs]
                yy <- y[select.obs]
                zz <- z[select.obs]
                dfd <- length(xx) - 1
                radius <- sqrt(dfn * qf(level, dfn, dfd))
                ellips <- ellipsoid(center=c(mean(xx), mean(yy), mean(zz)),
                    shape=cov(cbind(xx,yy,zz)), radius=radius)
                if (fill) shade3d(ellips, col=surface.col[j], alpha=0.1, lit=FALSE)
                if (grid) wire3d(ellips, col=surface.col[j], lit=FALSE)
                coords <- ellips$vb[, which.max(ellips$vb[1,])]
                if (!surface) rgl.texts(coords[1] + 0.05, coords[2], coords[3], group,
                    col=surface.col[j])
                }
            }
        }
    if (surface){
        vals <- seq(0, 1, length=grid.lines)
        dat <- expand.grid(x=vals, z=vals)
        for (i in 1:length(fit)){
            f <- match.arg(fit[i], c("linear", "quadratic", "smooth", "additive"))
            if (is.null(groups)){
                mod <- switch(f,
                    linear = lm(y ~ x + z),
                    quadratic = lm(y ~ (x + z)^2 + I(x^2) + I(z^2)),
                    smooth = if (is.null(df.smooth)) gam(y ~ s(x, z))
                        else gam(y ~ s(x, z, fx=TRUE, k=df.smooth)),
                    additive = if (is.null(df.additive)) gam(y ~ s(x) + s(z))
                        else gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
                            s(z, fx=TRUE, k=(rev(df.additive+1)[1]+1)))
                    )
                if (model.summary) summaries[[f]] <- summary(mod)
                yhat <- matrix(predict(mod, newdata=dat), grid.lines, grid.lines)
                if (fill) rgl.surface(vals, vals, yhat, color=surface.col[i],
                    alpha=0.5, lit=FALSE)
                if(grid) rgl.surface(vals, vals, yhat, color=if (fill) grid.col
                    else surface.col[i], alpha=0.5, lit=FALSE, front="lines", back="lines")
                if (residuals){
                    n <- length(y)
                    fitted <- fitted(mod)
                    colors <- ifelse(residuals(mod) > 0, pos.res.col, neg.res.col)
                    rgl.lines(as.vector(rbind(x,x)), as.vector(rbind(y,fitted)),
                        as.vector(rbind(z,z)), color=as.vector(rbind(colors,colors)))
                    }
                }
            else{
                if (parallel){
                    mod <- switch(f,
                        linear = lm(y ~ x + z + groups),
                        quadratic = lm(y ~ (x + z)^2 + I(x^2) + I(z^2) + groups),
                        smooth = if (is.null(df.smooth)) gam(y ~ s(x, z) + groups)
                            else gam(y ~ s(x, z, fx=TRUE, k=df.smooth) + groups),
                        additive = if (is.null(df.additive)) gam(y ~ s(x) + s(z) + groups)
                            else gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
                                s(z, fx=TRUE, k=(rev(df.additive+1)[1]+1)) + groups)
                        )
                    if (model.summary) summaries[[f]] <- summary(mod)
                    levs <- levels(groups)
                    for (j in 1:length(levs)){
                        group <- levs[j]
                        select.obs <- groups == group
                        yhat <- matrix(predict(mod, newdata=cbind(dat, groups=group)),
                            grid.lines, grid.lines)
                        if (fill) rgl.surface(vals, vals, yhat, color=surface.col[j],
                            alpha=0.5, lit=FALSE)
                        if (grid) rgl.surface(vals, vals, yhat, color=if (fill) grid.col
                            else surface.col[j], alpha=0.5, lit=FALSE,
                                front="lines", back="lines")
                        rgl.texts(0, predict(mod, newdata=data.frame(x=0, z=0,
                            groups=group)), 0,
                            paste(group, " "), adj=1, color=surface.col[j])
                        if (residuals){
                            yy <- y[select.obs]
                            xx <- x[select.obs]
                            zz <- z[select.obs]
                            fitted <- fitted(mod)[select.obs]
                            rgl.lines(as.vector(rbind(xx,xx)),
                                as.vector(rbind(yy,fitted)), as.vector(rbind(zz,zz)),
                                col=surface.col[j])
                            }
                        }
                    }
                else {
                    levs <- levels(groups)
                    for (j in 1:length(levs)){
                        group <- levs[j]
                        select.obs <- groups == group
                        mod <- switch(f,
                            linear = lm(y ~ x + z, subset=select.obs),
                            quadratic = lm(y ~ (x + z)^2 + I(x^2) + I(z^2),
                                    subset=select.obs),
                            smooth = if (is.null(df.smooth)) gam(y ~ s(x, z),
                                    subset=select.obs)
                                else gam(y ~ s(x, z, fx=TRUE, k=df.smooth),
                                    subset=select.obs),
                            additive = if (is.null(df.additive)) gam(y ~ s(x) + s(z),
                                    subset=select.obs)
                                else gam(y ~ s(x, fx=TRUE, k=df.additive[1]+1) +
                                    s(z, fx=TRUE, k=(rev(df.additive+1)[1]+1)),
                                    subset=select.obs)
                            )
                        if (model.summary)
                            summaries[[paste(f, ".", group, sep="")]] <- summary(mod)
                        yhat <- matrix(predict(mod, newdata=dat), grid.lines, grid.lines)
                        if (fill) rgl.surface(vals, vals, yhat, color=surface.col[j],
                            alpha=0.5, lit=FALSE)
                        if (grid) rgl.surface(vals, vals, yhat, color=if (fill) grid.col
                            else surface.col[j], alpha=0.5, lit=FALSE,
                                front="lines", back="lines")
                        rgl.texts(0, predict(mod,
                            newdata=data.frame(x=0, z=0, groups=group)), 0,
                            paste(group, " "), adj=1, color=surface.col[j])
                        if (residuals){
                            yy <- y[select.obs]
                            xx <- x[select.obs]
                            zz <- z[select.obs]
                            fitted <- fitted(mod)
                            rgl.lines(as.vector(rbind(xx,xx)), as.vector(rbind(yy,fitted)),
                                as.vector(rbind(zz,zz)),
                                col=surface.col[j])
                            }
                        }
                    }
                }
            }
        }
    if (revolutions > 0) {
        for (i in 1:revolutions){
            for (angle in seq(1, 360, length=360/speed)) rgl.viewpoint(-angle, fov=fov)
            }
        }
    if (model.summary) return(summaries) else return(invisible(NULL))
    }

identify3d <- function (x, y, z, groups = NULL, labels = 1:length(x), col = c("blue",
    "green", "orange", "magenta", "cyan", "red", "yellow", "gray"),
    offset = ((100/length(x))^(1/3)) * 0.02) {
    ## added by GM Nov 9
       if ( ( is.matrix(x) && (ncol(x) > 1) ) || is.data.frame(x) ) {
            dat <- x
            x <- dat[,1]
            y <- dat[,2]
            z <- dat[,3]
            labels <- rownames(dat)
       }
    select3d <-function (...) {
    # adapted from select3d and rgl.selecte3d in the rgl package, but passes
    #  through arguments to rgl.select
        rgl:::.check3d()
        rect <- rgl:::rgl.select(...)
        llx <- rect[1]
        lly <- rect[2]
        urx <- rect[3]
        ury <- rect[4]
        if (llx > urx) {
            temp <- llx
            llx <- urx
            urx <- temp
            }
        if (lly > ury) {
            temp <- lly
            lly <- ury
            ury <- temp
            }
        proj <- rgl:::rgl.projection()
        function(x, y, z) {
            pixel <- rgl.user2window(x, y, z, proj = proj)
            apply(pixel, 1, function(p) (llx <= p[1]) && (p[1] <=
                urx) && (lly <= p[2]) && (p[2] <= ury) && (0 <= p[3]) &&
                (p[3] <= 1))
            }
        }
    valid <- if (is.null(groups)) complete.cases(x, y, z)
    else complete.cases(x, y, z, groups)
    x <- x[valid]
    y <- y[valid]
    z <- z[valid]
    labels <- labels[valid]
    x <- (x - min(x))/(max(x) - min(x))
    y <- (y - min(y))/(max(y) - min(y))
    z <- (z - min(z))/(max(z) - min(z))
    rgl.bringtotop()
    identified <- character(0)
    groups <- if (!is.null(groups)) as.numeric(groups[valid])
    else rep(1, length(x))
    repeat {
        f <- select3d(button="middle")
        which <- f(x, y, z)
        if (!any(which)) break
        rgl.texts(x[which], y[which] + offset, z[which], labels[which],
            color = col[groups][which])
        identified <- c(identified, labels[which])
        }
    unique(identified)
    }

plot3d <- function(x, y, z,  ...) {
   help = "
   plot3d()   by GM Nov 9, 2005
   A wrapper for scatter3d by John Fox: added to coursefun.R
   "
   if ( is.matrix(x) || is.data.frame(x)) {
      if ( is.null( colnames( x ))) {
          colnames(x) <- paste( deparse(substitute(x)), 1:ncol(x))
      }
      scatter3d(x[,1],x[,2],x[,3],xlab = colnames(x)[1],
        ylab=colnames(x)[2], zlab = colnames(x)[3],...)
   } else {
      scatter3d(x,y,z, xlab = deparse(substitute(x)), ylab = deparse(substitute(y)), zlab = deparse(substitute(z)), ... )
   }
}


"

Ellipses: data and confidence

"


ell <- function(center = rep(0,2) , shape = diag(2) , radius = 1, n = 100) {
       angles <- (0:n)*2*pi/n
       circle <- radius * cbind( cos(angles), sin(angles))
       t( c(center) + t( circle %*% chol(shape)))
}


cell <-
function (model, which.coef, levels = 0.95, Scheffe = FALSE, dfn = 2,
    center.pch = 19, center.cex = 1.5, segments = 51, xlab, ylab,
    las = par("las"), col = palette()[2], lwd = 2, lty = 1,
    add = FALSE, ...)
{
help <- "
See help for car::confidence.ellipse.lm
except that 'cell' returns the points to form the ellipse
which must be plotted with plot(...,type='l') or lines(...)
-- Use dfn to determine Sheffe dimension, i.e. dfn = 1 to generate ordinary CIs, dfn = 2 for 2-dim CE, etc.
"
    require(car)
    which.coef <- if (length(coefficients(model)) == 2)
        c(1, 2)
    else {
        if (missing(which.coef)) {
            if (has.intercept(model))
                c(2, 3)
            else c(1, 2)
        }
        else which.coef
    }
    coef <- coefficients(model)[which.coef]
    xlab <- if (missing(xlab))
        paste(names(coef)[1], "coefficient")
    ylab <- if (missing(ylab))
        paste(names(coef)[2], "coefficient")
    if(missing(dfn)) {
        if (Scheffe) dfn <- sum(df.terms(model))
        else 2
    }
    dfd <- df.residual(model)
    shape <- vcov(model)[which.coef, which.coef]
    ret <- numeric(0)
    for (level in rev(sort(levels))) {
        radius <- sqrt(dfn * qf(level, dfn, dfd))
        ret <- rbind(ret, c(NA,NA), ell( coef, shape, radius) )
    }
    colnames(ret) <- c(xlab, ylab)
    ret
}

dell <- function( x, y, radius = 1, ...) {
    if ( (is.matrix(x) && (ncol(x) > 1))|| is.data.frame(x)) mat <- as.matrix(x[,1:2])
    else if (is.list(x)) mat <- cbind(x$x, x$y)
    else mat <- cbind( x,y)
    ell( apply(mat,2,mean), var(mat), radius = radius, ...)

}



"

Diagnostics

"


diags <- function(x, ...) UseMethod("diags")


diags.lm <- function(x, y, ..., ask, labels = names(residuals(x)), showlabs = text)
{
# diags.lm
# graphical diagnostics for lm, locally first-order for glm
# enlarged version of plot.lm with emphasis on diagnostics
# G. Monette, Dec. 94
# modified Nov. 97, May 98
# Slight modification to pairs adding labels, Jan 03
	if(!missing(ask)) {
		op <- par(ask = ask)
		on.exit(par(op))
	}
	form <- formula(x)
	f <- predict(x)
	r <- residuals(x)
	nams <- names(r)
	if(!missing(labels)) {
		nams <- names(residuals(x))	#
# if labels not same length as residuals assume it's a vector
# of len == original data and select elements included in residuals
		if(length(nams) != length(labels))
			labels <- labels[nams]
	}
	ret <- NULL
	if(missing(y)) {
		y <- f + r
		yname <- deparse(form[[2]])
	}
	else yname <- deparse(substitute(y))
	fname <- paste("Fitted:", deparse(form[[3]]), collapse = " ")
	plot(f, y, xlab = fname, ylab = yname, main = "Dependent var. vs. Predicted",
		...)
	abline(0, 1, lty = 1)
	lines(supsmu(f,y))
	showlabs(f, y, labels,...)
#
# get influence diags and model matrix while looking at first plot
#
	lmi <- lm.influence(x)
	hat <- lmi$hat
	sigma <- lmi$sigma	# drop 1 sigma
	mm <- scale(model.matrix(x), scale = F)	# centres each column
	mp <- predict(x, type = "terms")
	comp.res <- mp + r	# effect + residual
#
# Absolute residual vs. predicted
#
	plot(f, abs(r), xlab = fname, ylab = deparse(substitute(abs(resid(x)))),
		main = "Absolute Residual vs. Predicted", ...)
	showlabs(f, abs(r), labels, ...)	#
#
# Normal quantile plot
#
	zq <- qqnorm(r, main = "Normal Quantile Plot", ylab = "Residual", sub
		 = fname)
	qqline(r)
	showlabs(zq, labels,...)	#
#
# Symmetry plot of residuals (Lawrence C. Hamilton, Regression with
#       Graphics, Duxbury, 1992)
	n <- length(r)
	r.o <- sort(r)
	half <- (n + 1)/2
	if(n %% 2 == 1) {
# n is odd
		med <- r.o[half]
		below <- med - r.o[half:1]
		above <- r.o[half:n] - med
	}
	else {
# n is even
		med <- sum(r.o[c(half, half + 1)])/2
		below <- med - r.o[(n/2):1]
		above <- r.o[(n/2 + 1):n] - med
	}
	opt <- par(pty = "s")
	ran <- range(c(below, above))
	plot(below, above, main = "Symmetry plot of residuals", xlab =
		"Distance below median", ylab = "Distance above median", xlim
		 = ran, ylim = ran)
	abline(0, 1, lty = 2)
	par(opt)	#
#
# Studentized residual vs. leverage
#
	std.r <- r/(sigma * sqrt(1 - hat))
	plot(hat, std.r, xlab = "Leverage (hat)", ylab = yname, sub = fname,
		main = "Studentized residual vs. Leverage", ...)
	showlabs(hat, std.r, labels,...)	#	plot(lmi$sig, std.r)	#

#
# effect of dropping one observation DFBETA
#
	nams <- dimnames(lmi$coefficients)[[1]]
	pairs(lmi$coefficients)
	pairs(lmi$coefficients, panel = function(x,y,nams){
		points(x,y)
		text(x,y,nams)
	}, nams = nams)

	# main = "Effect of dropping one case", sub = fname)
	invisible(0)
}

"

vplot -- a plot function for matrix algebra

"

##
##  vplot  plots columns of a 2 x n matrix
##  Transferred to coursfun: Nov. 15, 2005

vplot <- function( mat , type = 'p', new = F,  pch = 16, pop = 0, ...) {
"
vplot    - plots the columns of a 2 x n matrix or a vector of length 2
         - vplot adds to the current plot resizing it to include all plotted
           objects in a 'euclidean' frame
         - to start a new plot, use 'new = T'
         - to remove the last element added use 'vplot(pop=1)'
         Associated functions:
         - vell( mean, var) generates an ellipse, default = unit circle
         - vbox() generates a box
         - vobj() generates a circle in a box
         - orthog(theta) generates an orthog matrix rotating through angle theta
         - orthog.proj generates the matrix of an orthog. projection into span (x)
         - vmat( .... ) generates a 2 by n matrix
         Examples:
           vplot( new = T )
           vplot( vell(), 'l' )
           vplot( cbind(c(3,1),c(1,4)) %*% vell())
           vplot( pop = 1)
           vplot( cbind(c(3,1),c(1,4)) %*% vell(), type = 'l', col = 'red')
"
     if (  new || !exists(".vplot")) assign(".vplot", list(list(x=0,y=0,type='n')),pos=1)
     a <- .vplot
     if ( ! missing(mat) ) {
        mat <- cbind(mat)
        if ( type == 'v' ) {
           zz <- rbind( 0*mat, mat, mat/0)
           mat <- matrix( c(zz), nrow = 2)
           type = 'b'
        }
        d <- dim(mat)
        if ( d[1] != 2 && d[2] == 2){
           mat <- t(mat)
           warning("mat is n x 2 and has been transposed")
        }
        a <- c(a,list( list(x=mat[1,],y = mat[2,],type=type, pch = pch, ...)))
     }
     dat <- NULL
     for ( i in seq( along = a )) {
         dat <- c( dat, a[[i]]$x, a[[i]]$y)
     }
     # print(a)
     par ( pty = 's')
     plot( range(na.omit(dat)), range(na.omit(dat)), type = 'n', xlab = '', ylab ='')
     if ( pop > 0 ) {
        keep <- 1:max(1,(length(a)-(pop+1)))
        a <- a[keep]
     }
     abline( h = 0, v = 0)
     for ( i in seq( along = a)) do.call('points', a[[i]])
     assign(".vplot", a, pos = 1)
     invisible(a)
}

vell <- function(...) t( ell(...))
vbox <- function(...) cbind( c(-1,-1), c(-1,1), c(1,1), c(1,-1), c(-1,-1))
vobj <- function(...) {
     cbind( vell(), NA, vbox(), NA, c(0,-1),c(0,1), NA, c(-1,0), c(1,0))
}
vsquare <- function(...) vmat( 0,0,0,1,1,1,1,0,0,0)

vmat <- function(...) {
help <- "
vmat creates a matrix entering data column by column
"
     aa <- list(...)
     aa <- do.call('c', aa)
     matrix(aa, nrow = 2)
}

orthog <- function( theta ) cbind( c( cos(theta), sin(theta)), c( - sin(theta), cos(theta)))
orthog.proj <- function ( x ) {
       x <- cbind(x)
       x %*% solve(t(x) %*% x , x)
}


"

Read.spss and Read.dta

"

###
###  trim
###

  trim <- function(x) {
       help <- "
trim in fun.R
  removes trailing blanks from character variables or from
  factor levels
  Use mainly to trim .dta files produced with SPSS
" 
       UseMethod("trim")
}

  trim.data.frame <- function(x) {
      for ( nn  in names(x)) x[[nn]] <- trim(x[[nn]])
      x
  }
  trim.factor <- function( x ) {
      levels(x) <- trim(levels(x))
      x
  }
  trim.character <- function( x ) {
      x[] <- sub(" +$", "", x )
      x
  }
  trim.default <- function(x) x
  
  Read.spss <- function( ... ) {
            require("Hmisc")
            dd <- spss.get ( ... )
            trim( dd )
  }
  
  Read.dta <- function ( ... ) {
help <- "
  Read.dta reads Stata files using 'read.dta' in 'library(foreign)'
  This appears to be an ideal way of importing spss files in order
  to keep full variable names. Direct use of 'read.spss' on a SPSS
  '.sav' file abbreviates variable names to 8 characters.
"
           require("foreign")
           dd <- read.dta(... , missing.type = T)
           cls <- lapply(dd,class)
           ch.nams <- names(dd) [ cls == "character" ]
           for ( nn in ch.nams ) dd[[nn]] <- trim(factor(dd[[nn]]) )
           dd
  }
  
  
"

RDC functions -- needs to be organized

"

   capply <- function( x , by, FUN, ... ) tapply( x, by, FUN, ...) [ tapply( x, by ) ]

   ch <- as.character

   select.first <- function( x , ... ) {
       ret <- unique( x [ !is.na(x) ])
       if ( length(ret) > 1 ) warning(paste("Multiple values:", paste(ret, collapse = ", ")))
       ret [1]
   }

   fill <- function( x, ... ) UseMethod("fill")

   fill.factor <- function( x, by, FUN = select.first, ...) {
      levs <- levels(x)
      ret <- capply( ch( x), by, FUN , ... )
      factor( ret, levels = levs)
   }

   fill.Date <- function( x, by, FUN = select.first, ...) {
      as.Date( capply(ch(x), by, FUN, ...))
   }

   fill.default <- function( x , by , FUN = select.first, ...) {
      capply( x, by, FUN, ...)
   }

   fill.data.frame <- function ( x, by ,... ) {
          ret <- list()
          for ( nn in names(x)) {
              ret[[nn]] <- fill( x[[nn]], by)
          }
          as.data.frame(ret)
   }



na20 <- function(x) {
     x[is.na(x)] <- 0
     x
}

describe.vector <-
function (x, descript, exclude.missing = TRUE, digits = 4, weights = NULL,
    normwt = FALSE, ...)
{
    # GM: modified by rounding n and missing
    oldopt <- options(digits = digits)
    on.exit(options(oldopt))
    if (length(weights) == 0)
        weights <- rep(1, length(x))
    special.codes <- attr(x, "special.miss")$codes
    labx <- attr(x, "label")
    if (missing(descript))
        descript <- as.character(sys.call())[2]
    if (length(labx) && labx != descript)
        descript <- paste(descript, ":", labx)
    un <- attr(x, "units")
    if (length(un) && un == "")
        un <- NULL
    fmt <- attr(x, "format")
    if (length(fmt) && (is.function(fmt) || fmt == ""))
        fmt <- NULL
    if (length(fmt) > 1)
        fmt <- paste(as.character(fmt[[1]]), as.character(fmt[[2]]))
    present <- if (all(is.na(x)))
        rep(FALSE, length(x))
    else if (is.character(x))
        (if (.R.)
            x != "" & x != " " & !is.na(x)
        else x != "" & x != " ")
    else !is.na(x)
    present <- present & !is.na(weights)
    if (length(weights) != length(x))
        stop("length of weights must equal length of x")
    if (normwt) {
        weights <- sum(present) * weights/sum(weights[present])
        n <- sum(present)
    }
    else n <- round(sum(weights[present]),2)
    if (exclude.missing && n == 0)
        return(structure(NULL, class = "describe"))
    missing <- round(sum(weights[!present], na.rm = TRUE),2)
    atx <- attributes(x)
    atx$names <- atx$dimnames <- atx$dim <- atx$special.miss <- NULL
    atx$class <- atx$class[atx$class != "special.miss"]
    isdot <- testDateTime(x, "either")
    isdat <- testDateTime(x, "both")
    x <- x[present, drop = FALSE]
    x.unique <- sort(unique(x))
    weights <- weights[present]
    n.unique <- length(x.unique)
    attributes(x) <- attributes(x.unique) <- atx
    isnum <- (is.numeric(x) || isdat) && !is.category(x)
    timeUsed <- isdat && testDateTime(x.unique, "timeVaries")
    z <- list(descript = descript, units = un, format = fmt)
    counts <- c(n, missing)
    lab <- c("n", "missing")
    if (length(special.codes)) {
        tabsc <- table(special.codes)
        counts <- c(counts, tabsc)
        lab <- c(lab, names(tabsc))
    }
    if (length(atx$imputed)) {
        counts <- c(counts, length(atx$imputed))
        lab <- c(lab, "imputed")
    }
    if (length(pd <- atx$partial.date)) {
        if ((nn <- length(pd$month)) > 0) {
            counts <- c(counts, nn)
            lab <- c(lab, "missing month")
        }
        if ((nn <- length(pd$day)) > 0) {
            counts <- c(counts, nn)
            lab <- c(lab, "missing day")
        }
        if ((nn <- length(pd$both)) > 0) {
            counts <- c(counts, nn)
            lab <- c(lab, "missing month,day")
        }
    }
    if (length(atx$substi.source)) {
        tabss <- table(atx$substi.source)
        counts <- c(counts, tabss)
        lab <- c(lab, names(tabss))
    }
    counts <- c(counts, n.unique)
    lab <- c(lab, "unique")
    x.binary <- n.unique == 2 && isnum && x.unique[1] == 0 &&
        x.unique[2] == 1
    if (x.binary) {
        counts <- c(counts, sum(weights[x == 1]))
        lab <- c(lab, "Sum")
    }
    if (isnum) {
        xnum <- if (.SV4.)
            as.numeric(x)
        else oldUnclass(x)
        if (isdot) {
            dd <- sum(weights * xnum)/sum(weights)
            fval <- formatDateTime(dd, atx, !timeUsed)
            counts <- c(counts, fval)
        }
        else counts <- c(counts, format(sum(weights * x)/sum(weights),
            ...))
        lab <- c(lab, "Mean")
    }
    if (n.unique >= 10 & isnum) {
        q <- if (any(weights != 1))
            wtd.quantile(xnum, weights, normwt = FALSE, na.rm = FALSE,
                probs = c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95))
        else quantile(xnum, c(0.05, 0.1, 0.25, 0.5, 0.75, 0.9,
            0.95), na.rm = FALSE)
        fval <- if (isdot)
            formatDateTime(q, atx, !timeUsed)
        else format(q, ...)
        counts <- c(counts, fval)
        lab <- c(lab, ".05", ".10", ".25", ".50", ".75", ".90",
            ".95")
    }
    names(counts) <- lab
    z$counts <- counts
    counts <- NULL
    if (n.unique >= 20) {
        if (isnum) {
            r <- range(xnum)
            xg <- pmin(1 + floor((100 * (xnum - r[1]))/(r[2] -
                r[1])), 100)
            z$intervalFreq <- list(range = as.single(r), count = as.integer(tabulate(xg)))
        }
        lo <- x.unique[1:5]
        hi <- x.unique[(n.unique - 4):n.unique]
        fval <- if (isdot)
            formatDateTime(c(oldUnclass(lo), oldUnclass(hi)),
                atx, !timeUsed)
        else format(c(format(lo), format(hi)), ...)
        counts <- fval
        names(counts) <- c("L1", "L2", "L3", "L4", "L5", "H5",
            "H4", "H3", "H2", "H1")
    }
    if (n.unique > 1 && n.unique < 20 && !x.binary) {
        tab <- wtd.table(if (isnum)
            format(x)
        else x, weights, normwt = FALSE, na.rm = FALSE, type = "table")
        pct <- round(100 * tab/sum(tab))
        counts <- t(as.matrix(tab))
        counts <- rbind(counts, pct)
        dimnames(counts)[[1]] <- c("Frequency", "%")
    }
    z$values <- counts
    structure(z, class = "describe")
}


na.include  <- function (obj)
{
     # from library(Hmisc)
    if (inherits(obj, "data.frame"))
        for (i in seq(along = obj)) obj[[i]] <- na.include(obj[[i]])
    else {
        if (length(levels(obj)) && any(is.na(obj)))
            obj <- factor(obj, exclude = NULL)
    }
    obj
}


summ <- function(x,...) UseMethod("summ")

summ.lmer <- function(x, ...) {
               ret <- c(AIC = AIC(x@logLik), BIC= BIC(x@logLik), logLik=x@logLik)
               ret
}

pr <- function(x,...) {
   # print to cut and paste as input
   UseMethod("pr")
}
pr.default <- function(x,pre="\t",...) {
          # cat('\nc(')
          for ( xx in x) cat(pre,'"',xx,'",\n',sep='')
          invisible(x)
}


reorder.factor <- function (x, v, FUN = mean, ...) {
        warning("gm version produces ordinary factor -- not ordered factor")
               factor(x, levels = levels(x)[order(tapply(v, x, FUN, ...))])
}



Q <- function(x, verbose = 0) {
    # find the Q matrix of a qr decomposition with possible NAs
    miss <- apply(x, 1, function(xx) any(is.na(xx)))
    xc <- x[!miss,]
    qf <- qr.Q(qr(xc))
    
    if(verbose > 0) {
               cat("xc:", dim(xc),'\n')
               cat("qf:", dim(qf), '\n')

               print(qf)
    }
    if( ncol(xc) > ncol(qf)) xc <- xc[,1:ncol(qf)]
    ip <- sign(apply(qf*xc,2,sum))   # sign of reln between Q and X
    qf <- qf * rep(ip, rep( nrow(qf), length(ip)))   # change sign of each row
    ret <- array(NA, dim = c(nrow(x),ncol(qf)))
    rownames(ret) <- rownames(x)
    colnames(ret) <- colnames(xc)
    ret[!miss,] <- qf
    ret
}

Contrasts <- function(x) UseMethod("Contrasts")

Contrasts.default <- function(x) contrasts(x)

Contrasts.lmer <- function(x) {
       dd <- x@frame[1,]
       ret <- list()
       for ( nn in names(dd)) {
           vv <- dd[[nn]]
           if (is.factor(vv)) ret[[nn]] <- contrasts(vv)
       }
       vv
}


###  getFix: function designed to be used internally to get coef, var(coef) and df.resid

getFix.rdc <- function(fit) UseMethod("getFix")

getFix.default <- function(fit) stop("not yet implemented")

getFix.rdc.lmer <- function(fit) {
            ret <- list()
            ret$fixed <- fixef(fit)
            ret$vcov <- vcov(fit)
            ret$df <- Matrix:::getFixDF(fit)
            ret
}

getFix.rdc.lm <- function(fit) {
          ret <- list()
          ret$fixed <- coef(fit)
          ret$vcov <- vcov(fit)
          ret$df <- fit$df.residuals
          ret
}

Lmat <- function(fit, pattern) {
     # pattern can be a character used as a regular expression in grep
     # of a list with each component generating  a row of the matrix
     umatch <- function( pat, x ) {
            ret <- rep(0,length(pat))
            for ( ii in 1:length(pat)) {
                imatch <- grep(pat[ii], x)
                if ( length(imatch) != 1) {
                   cat("\nBad match of:\n")
                   print(pat)
                   cat("in:\n")
                   print(x)
                   stop("Bad match")
                }
                ret[ii] <- imatch
            }
            ret
     }
     if ( is.character(fit)) {
        x <- pattern
        pattern <- fit
        fit <- x
     }
     fe <- getFix(fit)$fixed
     ne <- names(fe)
     if (is.character(pattern)) {
        L.indices <- grep(pattern,names(fe))
        ret <- diag( length(fe)) [L.indices,,drop = F]
        rownames(ret) <- names(fe) [L.indices]
     } else if (is.list(pattern)){
        ret <- matrix(0, nrow = length(pattern), ncol = length(fe))
        colnames(ret) <- ne
        for ( ii in 1:length(pattern)) {
            Lcoefs <- pattern[[ii]]
            pos <- umatch(names(Lcoefs), ne)
            if ( any( is.na(pos))) stop("Names of L coefs not matched in fixed effects")
            ret[ii, pos] <- Lcoefs
        }
        rownames(ret) <- names(pattern)
      }
      ret
}

Ldiff.old <- function(fit, pat, levnames = c(reflevel,substring(rownames(L),cut+1)),
                       reflevel = "<ref>", cut = nchar(pat)) {
         L <- Lmat(fit, pat)
         nam <- rownames(L)
         n <- nrow(L)
         if(n < 2) return(L)
         plus <- unlist( apply( rbind( 2:n), 2, seq, n))
         minus <- rep(1:(n-1), (n-1):1)
         Lp <- L[ plus, ]
         Lm <- L[ minus, ]
         Lret <- rbind( L, Lp - Lm)
         rn <- paste( levnames[ c(1:n,plus) + 1], levnames[ c(rep(0,n), minus)+1], sep = " - ")
         rownames(Lret) <- rn
         Lret
}

Ldiff.rdc <- function( fit, nam , ref = "no longer used") {
      # based on Lmu
      # Tests all pairwise difference in factor with model with Intercept term
      Lm <- Lmu(fit, nam)
      levs <- rownames(Lm)
      n <- nrow(Lm)
      if (n < 2) return (Lm)
      plus <- unlist( apply ( rbind(2:n), 2, seq, n))
      minus <- rep(1:(n-1), (n-1):1)
      Lret <- Lm[plus,] - Lm[minus,]
      rn <- paste( levs [plus], levs[minus] , sep = " - ")
      rownames(Lret) <- rn
      Lret
}




Lmu <- function(fit, nam, verbose = 0) {
       ## "Works only if 'nam' is a factor and a main effect and model has Intercept")
       if ( class(fit) != 'lmer' ) stop( "only implemented for lmer")
       v <- fit@frame[[nam]]
       if( !is.factor(v)) stop ("nam needs to specify the name of a factor")
       levs <- levels(v)
       if( verbose > 0) print(levs)
       cmat <- contrasts(v)
       if( verbose > 0) print(cmat)
       #  print(cmat)
       fe <- getFix(fit)$fixed
       if( verbose > 0) print(fe)
       if ( substring(nam,1,1) != '^') nam <- paste("^",nam,sep="")
       L.indices <- grep(nam,names(fe))
       if( verbose > 0) print(L.indices)
       L <- matrix(0,nrow=length(levs), ncol = length(fe))

       colnames(L) <- names(fe)
       if( verbose > 0) print(L)
          rownames(L) <- levs
       L[,L.indices] <- cmat
       if('(Intercept)' %in% colnames(L)) L[,'(Intercept)'] <- 1
       L
}



Lc <- function(fit, nam, ref = 1, verbose = 0) {
       ## Comparisons with one level
       ## Use Lmu
       ## "Works only if 'nam' is a factor and a main effect and model has Intercept?")
       if ( class(fit) != 'lmer' ) stop( "only implemented for lmer")
       L <- Lmu( fit, nam)
       Lref <- L[ ref,,drop = F]
       index <- 1:nrow(L)
       names(index) <- rownames(L)
       refind <- index[ref]
       if (length(refind) != 1) stop( paste( ref, "does not refer to a single level"))
       Lret <- L[-refind,]
       Lret <- Lret - cbind( rep(1,nrow(Lret))) %*% Lref
       attr(Lret,"heading") <- paste("Comparisons with reference level:", rownames(L)[refind])
       Lret
}

Lrm <- function(fit, nam, vals = 1:nrow(L.mu)) {
    ## Repeated measures polynomial contrasts
    ## Uses Lmu
    ## "Works only if 'nam' is a factor and a main effect and model has Intercept?")
    ##
    L.mu <- Lmu(fit, nam)
    # print(L.mu)
    pp <- cbind( 1, Poly(vals, nrow(L.mu) -1))
    # print(pp)
    ortho <- Q(pp)[,-1] # (linear, quad, etc.)
    # print(ortho)
    ortho <- ortho[,-1]
    maxp <- max( 5, nrow(L.mu))
    colnames(ortho) <- c('linear','quadratic','cubic', paste("^",4:maxp,sep=''))[1:ncol(ortho)]
    L <- t(ortho) %*% L.mu
    L
}
# Lrm(fit, "SRH94")




Lequal <- function(fit, pat) {
       # Test for equality within levels of pat using all differences
         L <- Lmat(fit, pat)
         nam <- rownames(L)
         n <- nrow(L)
         if(n < 2) return(L)
         plus <- unlist( apply( rbind( 2:n), 2, seq, n))
         minus <- rep(1:(n-1), (n-1):1)
         Lp <- L[ plus, ]
         Lm <- L[ minus, ]
         Lret <- rbind( Lp - Lm)
         rn <- paste( nam[plus], nam[minus], sep = " - ")
         rownames(Lret) <- rn
         Lret
}




#Lc <- function(fit, vec ){
#   fe <- getFix(fit)$fixed
#   ret <- 0 * fe
#   if ( is.null(names(vec))) ret[] <-
#}
# Lmu(fit,"SRH")

Lall <- function( fit , nam ) {
        if ( class(fit) != 'lmer' ) stop( "only implemented for lmer")
        v <- fit@frame[[nam]]
        if( !is.factor(v)) stop ("nam needs to specify the name of a factor")
        lev0 <- levels(v)[1]
        ret <-list()
        namf <- nam
        if ( substring(namf,1,1) != "^") namf <- paste("^", namf, sep ="")
        ret[[ nam ]] <- Lmat( fit, namf)
        ret[[ paste(nam,"mu",sep = '.') ]] <- Lmu(fit,nam)
        ret[[ paste(nam,"diff",sep = '.') ]] <- Ldiff( fit , nam)
        ret

}


comp.old <- function(fit, form, varname, varpattern = vname, data = fit@frame, ...) {
     ## Computing regression components
     ## currently works form lmer only
     ## idea is to return a data frame with value a variable and corresponding fitted values
     ## for the 'component' of that variable or variables
     ##
     ## This can be used to fit a model to a prediction data.frame with only
     ## the necessary variables defined!!
     ## However BEWARE not to use 'unsafe' transformations (e.g. Q, poly)
     ##
     model.mat <- model.matrix(form, data,...)
     #print(dim(model.mat))
     ret <- data[rownames(model.mat),varname,drop = F]
     fe <- fixef(fit)
     pos.fe <- grep( varpattern, names(fe))
     pos.mat <- grep( varpattern, colnames(model.mat))
     ret$comp <- c( model.mat[,pos.mat] %*% fe[pos.fe] )
     attr(ret,"predictors") <-  names(fe)[pos.fe]
     ret
}



comp <- function(fit, form, varname, varpattern = vname, data = fit@frame, ...) {
     ## Computing regression components
     ## currently works form lmer only
     ## idea is to return a data frame with value a variable and corresponding fitted values
     ## for the 'component' of that variable or variables
     ##
     ## This can be used to fit a model to a prediction data.frame with only
     ## the necessary variables defined!!
     ## However BEWARE not to use 'unsafe' transformations (e.g. Q, poly)
     ##
     ## Arguments:
     ##
     ## fit : model for which component to be estimated
     ##
     ## form : portion of model formula to generate model variables needed for fitting
     ##
     ## varname: variable names to be included in output data frame
     ##
     ## varpattern:  regular expression so select effects in component
     ##
     ## data:  the 'parent' data.frame -- often a prediction data frame.
     ##        When using the original data frame, it is often necessary to
     ##        'na.action = na.omit' for '...'
     ##
     model.mat <- model.matrix(form, data,...)
     #print(dim(model.mat))
     ret <- data[rownames(model.mat),varname,drop = F]
     fe <- fixef(fit)
     effect.names <- grep( varpattern, names(fe), value = T)
     ret$comp <- c( model.mat[,effect.names] %*% fe[effect.names] )
     attr(ret,"predictors") <-  effect.names
     ret
}


#  glh( fit, Lall(fit, 'Wave'))

cond <- function(x) {
     # reporting on conditioning of matrix
     Svd <- svd(x)
     ret <- list(svd = Svd, dim = dim(x), rank = qr(x)$rank, condition = Svd$d[1]/Svd$d[length(Svd$d)])
     if(nrow(x) == ncol(x)) ret <- c(ret, cor= list(cor(x)))
     ret
}



Vcov.rdc <- function( fit, L  = Lmat(fit,"") ) {
     # variance of etahat = L beta.hat
     vc <- getFix(fit)$vcov
     L %*% vc %*% t(L)
}

Vcor.rdc <- function(fit, L = Lmat(fit,"")) {
     ret <- cov2cor(vc <- Vcov(fit, L))
     sv <- svd(vc)$d

     attr(ret,'condition') <- sv[1]/sv[length(sv)]
     attr(ret,'class') <- "correl"
     ret
}


print.correl <- function(x) {
      correl <- format(round(x, 3), nsmall = 3,
                  digits = 4)
      correl[!lower.tri(correl)] <- ""
      if (!is.null(cond <- attr(correl,"condition"))) {
         attr(correl,"condition") <- NULL
      }
      print(correl, quote = FALSE)
      if(!is.null(cond)) cat("Condition:",cond,"\n")
      invisible(x)
}


glh.rdc <- function(fit, Llist, help = F, clevel = 0.95, debug = F) {
    if(help) {
         cat("help!!!")
         return(0)
    }
    if( !is.list(Llist) ) Llist <- list(Llist)
    ret <- list()
    fix <- getFix(fit)
    beta <- fix$fixed
    vc <- fix$vcov
    dfs <- fix$df
    for ( ii in 1:length(Llist)) {
        ret[[ii]] <- list()

        L <- rbind(zz <-  Llist[[ii]])
        heading <- attr(zz, "heading")
        nam <- names(Llist)[ii]
        
        ## Anova

        qqr <- qr(t(L))
       #if(debug) print(qqr)
        L.rank <- qqr$rank
        L.full <- t(qr.Q(qqr)) [ 1:L.rank,,drop=F]
        if(F) {
                  cat("\n+++++++++++++++++++\n")
                  print(nam)
                  print(L)
                  print(svd(L)$d)
                  print(L.full)
                  print(svd(L.full)$d)
         }
         if (debug) {
            print(L)
            print(dim(L.full))
            print(dim(vc))
         }
        vv <- L.full %*% vc %*% t(L.full)
        eta.hat <- L.full %*% beta
        Fstat <- (t(eta.hat) %*% qr.solve(vv, eta.hat))/L.rank
        Fstat2 <- (t(eta.hat) %*% solve(vv) %*% eta.hat)/L.rank
        included.effects <- apply(L, 2, function(x) sum(abs(x))) != 0
        denDF <- min(dfs[included.effects])
        numDF <- L.rank

        ret.anova <- rbind(c(numDF, denDF, Fstat,Fstat2, pf(Fstat, numDF, denDF, lower.tail = F)))
        colnames(ret.anova) <- c("numDF","denDF","F value","F2","Pr(>F)")
        rownames(ret.anova) <-  nam
        ret[[ii]]$anova <- ret.anova
        
        ## Estimate
        
        etahat <- L %*% beta
        etavar <- L %*% vc %*% t(L)
        etasd <- sqrt(diag(etavar))
        denDF <- apply( L , 1, function(x,dfs) min(dfs[x!=0]), dfs = dfs)
        aod <- cbind(
                 c(etahat),
                 etasd,
                 denDF,
                 c(etahat/etasd),
                 2*pt(-abs(etahat/etasd), denDF))
        colnames(aod) <- c("Estimate","Std.Error",'DF','t value','Pr(>|t|)')
        if( !is.null(clevel) ) {
            hw <- qt( 1-(1-clevel)/2, denDF) * etasd
            aod <- cbind(aod, LL = etahat - hw, UL = etahat + hw)
            labs <- paste( c("Lower","Upper"), format(clevel))
            colnames(aod)[ ncol(aod) + c(-1,0)] <- labs

        }
        #aod <- as.data.frame(aod)
        rownames(aod) <- rownames(L)
        ret[[ii]]$estimate <- aod
        
        ## Vcov
        
        ret[[ii]]$vcov <- Vcov( fit, L)
        ret[[ii]]$vcor <- Vcor(fit,L)
        ret[[ii]]$L <- L
        ret[[ii]]$L.full <- L.full
        if ( !is.null(heading)) attr(ret[[ii]],"heading") <- heading
    }
    names(ret) <- names(Llist)
    attr(ret,"class") <- "glh"
    ret
}





formatCoefmat <- function(x ,digits = 6, pdigits = digits-1 ) {
     pformat <- function(x, digits) {
             x <- format(xx <- round(x,digits))
             x[ as.double(xx) == 0 ] <- paste(c("<.",
                              rep('0',digits-1),"1"), collapse = "")
             x
     }
     xx <- array("",dim=dim(x), dimnames = dimnames(x))
     for ( i in 1:ncol(xx)) {
         xx[,i] <- format(round(x[,i],digits),digits = digits)
     }
     if ( length( isp <- grep("^Pr\\(",colnames(x))) > 0) {
         xx[,isp[1]] <- pformat( x[,isp[1]], digits = pdigits)
     }
     xx
}

print.glh.rdc <- function(x, round = 6, pround = round - 1, L  = T, cov = T, ...) {
     # should round by SD, i.e. keep 3 sig digits for sd and others rounded accordingly


     rnd <- function(x, digits) {
             if ( is.numeric( x)) x <- round(x, digits = digits)
             format(x)
     }
     for ( ii in 1:length(x)) {
         nn <- names(x)[[ii]]
         tt <- x[[ii]]
         ta <- tt$anova
         tap <- array("", dim = dim(ta), dimnames = dimnames(ta))

         # ta[["p-value"]] <- pformat( ta[["p-value"]], digits = pround)
         ## print(as.data.frame(ta,row.names = nn))
         cat("\n",nn,"\n",sep='')

         print(formatCoefmat(ta,digits=round,pdigits=pround),quote=F,right=T)
         cat("\n")
         te <- tt$estimate
         #tret <- te
         #mode(tret) <- 'character'
         #tret[,'p-value'] <- pformat( te[,'p-value'], digits = pround)
         #if( !is.null(round)) {
         #    for (i in 1:ncol(te) ) {
         #        tret[,i] <- rnd(te[,i], digits = round)
         #    }
         #print(tret,quote=F)
         if(!is.null(zhead <- attr(tt,'heading'))) cat(zhead,"\n")
         print(formatCoefmat(te,digits=round,pdigits=pround),quote=F,right=T)
         if (L == T ) {
            cat("\nL:\n")
            print(tt$L)
            if( dim(tt$L.full)[1] < dim(tt$L) [1]) {
             cat("\nL (full rank):\n")
             print(tt$L.full)
            }
         }
         if ( cov == T ) {
            cat("\nVar-Cov of estimates:\n")
            print(tt$vcov)
            cat("\nCorrelations:\n")
            print(tt$vcor)
         }

     }
     invisible(x)
}


##z <- glh(fit, list('language'=Lmat(fit, "language")))
##z

##z <- glh(fit, list('Language | Year=1998'=Ldiff(fit, "language",ref="Qc French Multi")))
##z


####################################################
####################################################  ABOVE is OKAY
xanova <- function(object,...) UseMethod("xanova")

#xanova.lmer <- function(fit, Llist ) {
#     if ( is.list(Llist) )
#}


xanova.lmer <- function( fit, Llist , df = NULL, clevel = .95) {
       # performs a Wald test on an object that has a fixef and a vcov methods
       warning("xanova.lmer uses Chi-Square tests")
       ret <- list()
       for ( ii in 1:length(Llist) ) {
           L <- rbind(Llist[[ii]])

           # anova step

           QR <- qr(L)
           R <- qr.R(QR)
           dfH <- QR$rank
           eta <- R %*% fixef(fit)
           vv <- R %*% vcov(fit) %*% t(R)
           chisq <- t(eta) %*% qr.solve(vv, eta)
           test <- list(ChiSquare = chisq, DF = dfH, "p-value" = 1-pchisq(chisq,dfH))
           ret[[ii]]$anova <- test
           
           # estimation
           
           eta <- L %*% fixef(fit)
           vv <- diag( L %*% vcov(fit) %*% t(L))
           etasd <- sqrt(vv)
           zval <- c(eta/etasd)
           aod <- cbind(Estimate=c(eta), Std.Error = etasd,
               "z-value" = zval, "p-value" = 2*pnorm(-abs(zval)))
           if( !is.null(clevel) ) {
               hw <- qnorm(1-(1-clevel)/2) * etasd
               aod <- cbind( aod, LL = eta - hw, UL = eta + hw)
               labs <- paste(c("Lower","Upper",format(clevel)))
               colnames(aod) [ ncol(aod) + c(-1,0)] <- labs
           }
           aod <- as.data.frame(aod)
           class(aod) <- c('estimate.lme','data.frame')
           
           ret[[ii]]$estimate <- aod
        }
}

####################################   NEED TO TEST ABOVE

            
           
           





enc <- function(x) {
		## this function will use the coding table in Stats Can documentation
		## to create a factor with the right names
		#
		tonum <- function(x) as.numeric(gsub(",","",as.character(x)))
		ff <- function(x) format(x, big.mark=',')
		tran.table <- scan(what=list('character','integer','character','character'),
			flush = T)
		# Brief report
		sam <- tonum(tran.table[[3]])
		pop <- tonum(tran.table[[4]])
		samp.pct <- 100 * (sam / pop) / ( sum(sam,na.rm=T)/sum(pop,na.rm=T)) 
		print( data.frame(
			Code = tran.table[[2]], Content = tran.table[[1]], Sample = ff(sam), Popn = ff(pop),
				Sampling.Pct = round(samp.pct,1))) 
		tran( tran.table[[2]], tran.table[[1]], x , tofactor = T)
	}

apct <- function(x,MARGIN=1) {
		if( length(dim(x)) == 1) {
			ret <- cbind( N = x, pct = 100*x/sum(x,na.rm=T))
			ret <- rbind( ret, Total = apply(ret,2,sum,na.rm=T))
			print( round(ret,1))
			return(invisible(ret))
		}	
		# report a table
		ret <- list( N=atotal(x), pct = 100 * acond(x,MARGIN) )		
		cat("\nN:\n")
		print( ret[[1]])
		cat("\nPercentage:\n")
		print( round( ret[[2]],1))
		invisible(ret)
	}

  arep <- apct   # old names

## From library gm

write.sas <- function( df , datafile="G:/SAS/R2SAS.txt", codefile = "G:/SAS/R2SAS.sas"){
	debug <- F
	pr <- function(x) {
		cat(deparse(substitute(x)),"\n")
		print(x)
		cat("==========\n")		
		cat(x)
		cat("\n=============\n")
	}
	
	lrecl <- 256
	if(!debug) { 
		write.table(df, file = datafile, row = FALSE, col = FALSE, 
        		sep = ";", na = ".")
		# compute lrecl
		lines <- scan(file = datafile, what = list("character"), sep = "\n")
		lrecl <- max(nchar(lines))
	}
    	nms <- names(df)
	nms.sas <- gsub("\\.","_",nms)
	## Check for duplicate SAS names

	if ( length( unique(toupper(nms.sas))) != length(nms.sas)) {
		ind <- duplicated( toupper(nms.sas))
		ind.rev <- duplicated( rev(toupper(nms.sas)))
		cat("Warning:\n")
		cat("The following R names may yield duplicate SAS names",
			"\n", paste(nms[ind | rev(ind.rev)],collapse=" "),"\n")
		warning("Possible duplicate SAS names")
	}
	
	factors <- sapply(df, is.factor) | sapply(df, is.character)
	## check for odd types
	classes <- sapply(df, class)
	odd.classes <- setdiff( sapply(df,class), c('numeric','factor','character') )
	if ( length(odd.classes) > 0 ) {
		cat("Warning:\n")
		cat("The following variables have classes that might not be handled properly by SAS\n")
		print( classes[ grep( odd.classes, classes)])
		cat("\n")
	}

	factor.names <- nms[factors]
	factor.names.sas <- nms.sas[factors]	
	dollarsign <- ifelse( factors, "$","")
	factor.lengths <- sapply( df[factor.names], function(x) {
		if(is.factor(x)) max(nchar(levels(x) ) ) else 
			max(nchar(x))
		})
	length.stmt <- paste(paste( "   ", factor.names.sas, "$", factor.lengths,"\n"),collapse = "")
	length.stmt <- paste( "LENGTH\n", length.stmt, ";\n")
	if (debug) pr(length.stmt)
	input.stmt <- paste(paste("    ", nms.sas, dollarsign,"\n"), collapse = "")
	input.stmt <- paste( "INPUT\n", input.stmt, ";\n")
	if (debug) pr(input.stmt)
	
	code <- paste("filename r2sas \'",datafile,"\';\n",   # might have to convert to backslashes
		"libname to \'G:/SAS\';\n",
		"data to.r2sas;\n",
		"infile r2sas delimiter=\';\' dsd LRECL =", lrecl+100, ";\n", sep = "")
	code <- paste(code,length.stmt, input.stmt, "\nrun;\n")
	if (debug) pr(code)
	if(!debug) cat(code, file = codefile)
	invisible(0)
}


## date()
## write.sas(dd[dd$wave > 1,setdiff(sort(names(dd)),c('qday','bday')) ]) # approx 1 min.
## date()
	 
if(F) { # current version in /R/coursefun.R
td <- function( basecol = NULL, col = c(3,5,4,6,7,8,2), lty = 1:7,
	lwd = 1, pch = 1:7, cex = 0.8, font = 1, len = 7, long = F,
	new = F, record = T,  theme = col.whitebg,
	col.symbol = col, col.line = col,
	colsets = c( "plot.symbol","plot.line", "dot.symbol","dot.line","cloud.3d","box.dot"),...) {
	require(lattice)
	if (long) {
		col <- c(3,5,4,6,8,2)
		len <- 42
	}
	if (new) trellis.device( theme = theme, record = record, new = new, ...)
	len <- max( len, length(col.symbol), length(col.line), length( lty), length(lwd), length(pch),
		length(cex), length(font))
	spl <- trellis.par.get("superpose.line")
	spl$lty <- if ( is.null(lty)) spl$lty else lty
	spl$col <- if ( is.null(col.line)) spl$col else col.line
	spl$lwd <- if ( is.null(lwd)) spl$lwd else lwd
	trellis.par.set("superpose.line", Rows(spl, 1:len))
	sps <- trellis.par.get("superpose.symbol")
	sps$pch <- if ( is.null(pch)) sps$pch else pch
	sps$col <- if ( is.null(col.symbol)) sps$col else col.symbol
	sps$cex <- if ( is.null(cex)) sps$cex else cex
	sps$font <- if ( is.null(font)) sps$font else font
	trellis.par.set("superpose.symbol", Rows(sps, 1:len))
	if( !is.null(basecol) ) {
		for(ii in colsets) {
			tt <- trellis.par.get(ii)
			tt$col <- basecol
			trellis.par.set(ii,tt)
		}
	}
	invisible( attr( .Device, "trellis.settings"))
}
} # end of F


cap1 <- function(x, tofactor = is.factor(x),
		stop=c(" The"," Of"," By"," To"," And"),
		blanks=c(" ","(","\"","/","+")) {
	# capitalizes first letters 
	under2blank <- T
	if ( is.factor(x)) {
		ret <- cap1(levels(x))
		if ( length(unique(ret)) != length(ret)) warning("factor levels have been shortened")
		levels(x) <- ret
		return(x)
	}
	ret <- as.character(x)
	for ( ii in 1:length(ret)) {
		z <- ret[[ii]]
		if(under2blank) z <- gsub("_"," ",z)
		n <- nchar(z)
		z <- substring( z, 1:n, 1:n)
		zu <- toupper(z)
		zl <- tolower(z)
		zb <- c(" ",zu[-n])
		z <- 	paste( ifelse(zb %in% blanks, zu, zl), collapse ="")	
		for ( ss in stop) z <- gsub(ss,tolower(ss), z)
		ret[[ii]] <- z
	}
	ret
}

tran <- function( from, to, x, tofactor = is.factor(x)) { 
	# Note that this uses the same argument order as grep and sub
	# which conflicts with my versions of tr, ergo this is called tran
	if(is.factor(from)) from <- as.character(from)
	if(is.factor(to)) to <- as.character(to)
	to <- rep(to, length=length(from))
	ret <- x
	if( is.factor(x) ) {
		ret <- as.character(ret)
		levs <- levels(x)
	}
	ret <- c(to,unique(ret)) [ match( ret, c(from, unique(ret)))]
	# DEBUG: print(ret)
	if (tofactor) {
		if(is.factor(x)) tolevs <- tran(from,to,levs)
		else tolevs <- to
		tolevs <- c(tolevs,unique(ret))
		tolevs <- intersect( tolevs, unique( ret ))
		ret <- factor(ret,levels = tolevs)
	}
	ret
} 





abind.rdc <- function( arr1, arr2, d, facename = "") {
	# glue arr1 to arr2 along dimension d (i.e. face 'not d')
	# copied from library gm 05 05 03
	d1 <- dim(arr1)
	n1 <- length(d1)	
	d2 <- dim(arr2)
	n2 <- length(d2)
	dn1 <- dimnames( arr1 )
	dn2 <- dimnames( arr2 )
	
	arenull <- is.null(dn1) & is.null(dn2)
	if ( is.null(dn1)) {
		dn1 <- lapply( as.list(d1), function(x) seq(1,x))
		dimnames(arr1) <- dn1
	}
	if ( n1 != n2 ) {
		d2 <- d1
		d2[d] <- 1
		dn2 <- dn1
		dn2[[d]] <- facename
		dimnames(arr2) <- NULL
		dim(arr2) <- d2
		dimnames(arr2) <- dn2
		n2 <- n1
	}
	if ( is.null(dn2)) {
		dn2 <- lapply( as.list(d2), function(x) seq(1,x))
		dimnames(arr2) <- dn2
	}
	perm <- 1:n1
	perm[c(d,n1)] <- c(n1,d)  # perm is an involution
	
	arr.p1 <- aperm(  arr1, perm )

	arr.p2 <- aperm(  arr2, perm )
	dret <- d1[perm]
	dret[n1] <- dret[n1] + d2[d]
	dnret <- dn1
	dnret[[d]] <- c(dnret[[d]],dn2[[d]])
	ret <- c(arr.p1, arr.p2)
	dim(ret) <- dret

	ret <- aperm( ret, perm )
	dimnames( ret ) <- dnret
	ret
}



# acond( with(dds, table( lang, prov)),2)

atotal.rdc <- function( arr, FUN = sum, name = "Total",...) {
	# copied from library gm  05 05 03
	# 05 05 03: added option for name
	d <- dim(arr)
	if ( length(d) == 1) {
		arr <- c(arr)
		d <- dim(arr)
	}
	if ( is.character( FUN ) ) FUN <- get(FUN,mode='function')
	else if( mode(FUN) != "function") {
		farg <- substitute(FUN)
		if (mode(farg) == 'name' ) FUN <- get(farg,mode='function')
		else stop(paste("\"", farg, "\" is not a function", sep=""))
	}
	if ( is.null(d) ) {
		ret <- c(arr, FUN(arr))
		names(ret)[length(ret)] = name
		return(ret)
	}
	n <- length (d)
	name <- rep(name, length= n)
	ret <- arr
	ind <- 1:n
	for ( i in n:1 ) {
		new <- apply(ret, ind[ -i ], FUN,...)
		ret <- abind( ret, new, i, name[i])
	}
	ret
}

aprop <- function( x, MARGIN = NULL ) {
	debug <- F
	## written 05 05 03: conditional proportions given a margin
	## like atotal but summing to 1 on selected margin
	x <- atotal(x)
	if ( is.null (MARGIN ) ) return( x / x[length(x)])
	# divide through by appropriate marginal total
	d <- dim(x)
	if (debug ) { cat("dim(x)\n"); print(d)}
	dn <- dimnames(x)
	if (debug ) { cat("dimnames(x)\n"); print(dn)}

	n <- length(d)
	m <- length(MARGIN)
	perm <- c( MARGIN, setdiff( 1:n, MARGIN))
	perm.inverse <- order(perm)
	if (debug ) print(perm)
	if (debug ) print(perm.inverse)
	x <- aperm( x, perm)
	d <- dim(x)
	zl <- list(x)
	for ( ii in 1:length(d)) zl[[ii+1]] <- seq( if(ii<=m) 1 else d[ii], d[ii])
	tots <- do.call("[",zl)
	ret <- array( c(x) / c(tots), dim = d)
	ret <- aperm( ret, perm.inverse)
	if(debug) print( dim(ret))
	if(debug) print( dn)
	dimnames(ret) <- dn
	ret
}

acond <- aprop    # old names



# Lagging
#    - interesting if 'index' is continuous and/or if there
#      are missing waves.
#    - need some kind of intrapolation to compute a derivative
#      at the time of observation

# The first function here only works well with complete data
# and an integer index ranging from 1 to n. (n can vary from 
# subject to subject)


Lag.0 <- function(x,id,idx,lag = 1,at= NULL, check=T) {
	# computes Lagged values but without intrapolation
	# adds 'at'

	if (check) {
		comb <- paste(id,idx)
		if (any(duplicated(comb))) stop("id not unique in each level of idx")
	}
	if(any(is.na(idx))) stop("NAs in idx")
	if(any( round(idx) != idx)) stop("Non integers in idx")
	ret <- x
	id <- as.character(id)
	names(x) <- id
	for ( i in max(idx):min(idx)) {
		to.pos <- idx == i
		if( is.null(at) ) from <- x[idx == (i-lag)]
		else from <- x[idx == at ]
		ids <- names( x[to.pos])
		ret[to.pos] <- from[ids]
	}
	ret
}

Lag <- function(x, id, idx, lag=1, at = NULL, check = T) {
   ## Takes approx. 30% of time of version Lag.0 on a small problem
  if (check) {
		comb <- paste(id,idx)
		if (any(duplicated(comb))) stop("id not unique in each level of idx")
	}
  ret <- x
  names(ret) <- paste(id,idx,sep='::')
  retnn <- if(is.null(at)) paste(id,idx - lag,sep='::')  else paste(id,at,sep="::")
  ret [ retnn ]
}
  
 if(F) {
## small test of Lag
zx <- c(2,3,2,4,2,4, 5,3,4,5,7,8,9)
zid<- c(1,1,2,2,3,3 ,3,4,4,4,4,4,4)
cbind( zid, zx, Lag(zx,zid,zx))
}


LagI <- function(x,id,time,lag=1,delta=.01,check=T) {
	# lags by intrapolating 
	# with complete data at each value of index, this does the same thing
	# as Lag
	# If values of Lag are skipped then we get linear intrapolations.
	# Note that 'delta' should be small enough so that values of x are
	# at least delta apart. However, too small a value for delta introduces
	# numerical error
	#  	
	if (check) {
		comb <- paste(id,time)
		if (any(duplicated(comb))) stop("id not unique in each level of idx")
	}
	ret <- x
	id <- as.character(id)
	names(x) <- id
	names(time) <- id
	for (nn in unique(id)){
		pos <- id == nn
		xx <- x[pos]
		tt <- time[pos]
		topred <- tt-delta
		drop <- is.na(xx)|is.na(tt)
		xxc <- xx[!drop]
		ttc <- tt[!drop]
		nl <- length(xxc)
		if ( nl > 0) {
			if ( nl > 1 ) xx.p <- approx(ttc,xxc,topred)$y
			else xx.p <- NA 
			xx.lag <- xx - lag*(xx - xx.p)/delta
			ret[pos] <- xx.lag
		}	
	}
	ret
}

DiffI <- function(xx,...) {
	xx - LagI(xx,...)
}

# tt <- c(1,2,3,1,2,3,1,2,3)
# xx <- c(1,2,3,1,2,3,1,2,3)
# id <- rep(letters[1:3],rep(3,3))
# xx-LagI(xx,id,tt)
# DiffI(xx,id,tt,delta=.1)

###
###  Splines and polynomial: quadratic and cubic
###


qs <- function(x, knots=quantile(x,pc), exclude = 0, pc = c(.25,.75)) {
	# simple quadratic spline
	ret <- cbind(x,x^2)
	nam <- deparse(substitute(x))

	for ( kk in knots ) {
		z <- x - kk
		z [z<0] <- 0
		ret <- cbind(ret, z^2)
	}

	if ( is.null(knots) ) dimnames(ret)[[2]] <- paste(nam,c('','^2'), sep='')
	else dimnames(ret)[[2]] <- paste(nam,c('','^2',paste('.',knots,sep='')),sep='')
	if (exclude > 0) ret <- ret[, -( 1:exclude)]
	ret

	ret
}



lsp <- function(x, knots=quantile(x,pc), exclude = 0, pc = c(.25,.75)) {
	# linear spline
	ret <- cbind(x)
	nam <- deparse(substitute(x))

	for ( kk in knots ) {
		z <- x - kk
		z [z<0] <- 0
		ret <- cbind(ret, z)
	}

	if ( is.null(knots) ) dimnames(ret)[[2]] <- paste(nam,c(''), sep='')
	else dimnames(ret)[[2]] <- paste(nam,c('',paste('.',knots,sep='')),sep='')
	if (exclude > 0) ret <- ret[, -( 1:exclude)]
	ret

	ret
}

cs <- function(x, knots=quantile(x,pc), exclude = 0, pc = c(.25,.75))  {
	# simple cubic spline
	ret <- cbind(x,(x)^2,(x)^3)
	nam <- deparse(substitute(x))
	for ( kk in knots ) {
		z <- (x) - kk
		z [z<0] <- 0
		ret <- cbind(ret, z^3)
	}
	if ( is.null(knots) ) dimnames(ret)[[2]] <- paste(nam,c('','^2','^3'), sep='')
	else dimnames(ret)[[2]] <- paste(nam,c('','^2','^3',paste('.',knots,sep='')),sep='')
	if (exclude > 0) ret <- ret[, -( 1:exclude)]
	ret
}

Poly <- function(x, order=2, exclude = 0)  {
	# polynomial
	# exclude: NULL to include intercept, otherwise exclude order up to 
	# including exclude
	ret <- cbind(rep(1,length(x)))
	for ( i in 1:order) ret <- cbind(ret, x^i)
	nam <- deparse(substitute(x))
	powers <- paste(nam,"^",0:order,sep="")
	powers[1] <- "Intercept"
	powers[2] <- nam
	dimnames(ret)[[2]] <- powers
	if (!is.null(exclude)) ret <- ret[, -( 1:(exclude+1))]
	ret
}






if (F) {
	## shows that columns of cs span same space as columns of bs	

	zz <- 0:20
	cs(zz, c(5,15))

	#search()
	#  detach(9)
	### comparison with ns

	zd <- data.frame( x = 0:20, y = (-10:10)^3 + rnorm(21))

	bss <- bs(zd$x,knots = c(5,15))
	css <- cs(zd$x,knots = c(5,15))

	fit <- lm(bss ~ css - 1)
	round(coef(fit),7)
	summary(fit)
	
	fit <- lm(css ~ bss - 1)
	summary(fit)
	round(coef(fit),7)

	fit <- lm(y ~ 1+cs(x,c(5,15)), zd)
	summary(fit)
	anova(fit)
	
	fit <- lm(y ~ x + I(x^2) + cs(x,c(5,15),2), zd)
	summary(fit)
	anova(fit)

	fit <- lm(y ~ x + I(x^2) + I(x^3) + cs(x,c(5,15),3), zd)
	summary(fit)
	anova(fit)
	
	fit <- lm(y ~ x + I(x^2) + cs(x,pc=c(.05,.95)), zd, singular.ok=T)
	summary(fit)
	anova(fit)
	
	qqr <- qr.Q(qr(cs(0:20,c(5,15))))
	fit <- lm( css ~ qqr-1)
	summary(fit)
	round(coef(fit),5) == round(qr.R(qr(cs(0:20, c(5,15)))),5)

}


## much better approach to stack:

mergec <- function( ... ,vname = '.type') {
	
	## stacks data frames adding a new variable .type to identify each source data frame
	## good way to combine data and prediction data frames to show data and fitted
	## values in a common plot
	
	z <- list(...)
	for ( ii in 1:length(z)) {
		z[[ii]][,vname] <- rep(ii, nrow( z[[ii]]))
	}
	ret <- z[[1]]
	for ( ii in 2:length(z)) ret <- merge(ret,z[[ii]], all = T)
	ret
}

if(F) {
	## test mergec
	zd1 <- data.frame(x=1:4, y = letters[1:4], z = 1:4,d=LETTERS[1:4])
	zd2 <- data.frame(x=5:7, y = letters[5:7])
	zd3 <- data.frame(x=8:9,  w=1:2, d=1:2)
	sapply(zz <- mergec(zd1,zd2,zd3),class)
	zz
	levels(zz$d)
}


##
##  constant: to check whether something is constant
##

constant <- function(x,...) UseMethod("constant")

constant.default <- function(x) length(unique(x)) == 1

constant.data.frame <- function( x, id  , all = FALSE) {
   ## Description:    (G. Monette, June 13, 2005)
   ## check if values of variables are constant within levels of id
   ## if no id, check each variable
   ## if all == TRUE, report overall instead of within levels of id
   ## Possible improvements:
   ##    allow nested formulas: ~id1/id2 for id and report level of
   ##    each variable
   ##    [see varLevel()]
   
   ## note that the following code allows id to be given as a name or as
   ## a formula
   
   if (missing(id)) ret <- sapply(x, constant)
   else {
        id <- eval(substitute(id), x, parent.frame())
        if( inherits(id,"formula") ) id <- c( model.frame(id,x) )
        ret <- sapply(x, function(xx) tapply(xx, id, constant))
        if ( all ) ret <- apply( ret, 2, all)
   }
   ret
}

varLevel <- function(x, form) {
         ## Description:    (G. Monette, June 13, 2005)
         ## shows levels of each variable with respect to grouping formula
         ## of form ~id or nested ids ~id1/id2
         ## Level 0 is a constant for the whole data frame
         ## Level 1 varies between levels of id1 but not within
         ## Level 2 varies between levels of id2 but not within, etc.
         ##
        sel <- model.frame( form, x )
        z <- list()
        idx <- rep('',nrow(x))
        z[[1]] <- constant(x)
        for ( ii in 1:ncol(sel)) {
            idx <- paste(idx, as.character(sel[[ii]]), sep = ";")
            z[[ii+1]] <- constant( x, idx, all = T)
        }
        # print(z)
        ret <- do.call("rbind", z)
        # print(ret)
        ret <- length(z) - apply( ret*1, 2 , sum)
        ret
}



###
###  Trellis additions
###


#Key <- function(text,type='l') {
#	## perhaps should use simpleKey instead
#	list(text=list(text),lines=Rows(trellis.par.get("superpose.line"),1:(length(text)))) 
#}

##

Key <- function(...) simpleKey(...)



###
###   Miscellaneous utility functions
###


val2na <- function( x, val) {
	## val2na(1:10, c(3,5))
	## val2na(factor(c('a','b','b','c','a')), c('b','a'))
	x[match(x,val,0)>0] <- NA
	x
}

tab <- function( ... ) {
	aa <- list(...)
	# print(aa)
	for ( ii in 1:length(aa) ) aa[[ii]] <- factor(aa[[ii]], exclude = NULL)
	ret <- do.call("table",aa)
	ret
}

# zf <- factor(c('a','a','b',NA,'c','a',NA))


grepv <- function(...) grep( ..., value = T)

p <- function(...) paste(..., sep ="")

ch <- function(x) as.character(x)

#########
#########   SCRAPS
#########

###
### Stack
###

## NOTE THAT R HAS A FUNCTION CALLED stack
if (F) {
stack <- function(...) {
	# simple version of stack, uses names of first df.
      # 05-04-19: THis function might be obsolete d/t mergec below
	ll <- list(...)
	llclass <- sapply(ll,function(x) inherits(x,'data.frame'))
	if(!all(llclass)){
		dfr <- ll[[1]]
		keep <- ll[[2]]
		add <- ll[[3]]
		ll <- list(0)
		for ( i in 1:length(add) ) ll[[i]] <- dfr[,c(keep,add[i])]
	} 
	nam <- names(ll[[1]])
	for ( ii in 2:length(ll)) names(ll[[ii]]) <- nam
	ret <- do.call('rbind',ll)
	nrows <- sapply(ll, function(x) dim(x)[1])
	ret$Type. <- rep(1:length(ll),nrows)
	ret
}
}



#####
#####  anova.lme
#####

anova.lme <- function (object, ..., test = TRUE, type = c("sequential", "marginal"),
    adjustSigma = TRUE, Terms, L, verbose = FALSE)
{
    warning("This is a modified version of anova.lme that uses min dfs for the denominator")
    Lmiss <- missing(L)
    dots <- list(...)
    if ((rt <- (length(dots) + 1)) == 1) {
        if (!inherits(object, "lme")) {
            stop("Object must inherit from class \"lme\" ")
        }
        vFix <- attr(object$fixDF, "varFixFact")
        if (object$method == "ML" && adjustSigma == TRUE) {
            vFix <- sqrt(object$dims$N/(object$dims$N - ncol(vFix))) *
                vFix
        }
        c0 <- solve(t(vFix), fixef(object))
        assign <- attr(object$fixDF, "assign")
        nTerms <- length(assign)
        if (missing(Terms) && Lmiss) {
            type <- match.arg(type)
            Fval <- Pval <- double(nTerms)
            nDF <- integer(nTerms)
            dDF <- object$fixDF$terms
            for (i in 1:nTerms) {
                nDF[i] <- length(assign[[i]])
                if (type == "sequential") {
                  c0i <- c0[assign[[i]]]
                }
                else {
                  c0i <- c(qr.qty(qr(vFix[, assign[[i]], drop = FALSE]),
                    c0))[1:nDF[i]]
                }
                Fval[i] <- sum(c0i^2)/nDF[i]
                Pval[i] <- 1 - pf(Fval[i], nDF[i], dDF[i])
            }
            aod <- data.frame(nDF, dDF, Fval, Pval)
            dimnames(aod) <- list(names(assign), c("numDF", "denDF",
                "F-value", "p-value"))
            attr(aod, "rt") <- rt
        }
        else {
            nX <- length(unlist(assign))
            if (Lmiss) {
                if (is.numeric(Terms) && all(Terms == as.integer(Terms))) {
                  if (min(Terms) < 1 || max(Terms) > nTerms) {
                    stop(paste("Terms must be between 1 and",
                      nTerms))
                  }
                }
                else {
                  if (is.character(Terms)) {
                    if (any(noMatch <- is.na(match(Terms, names(assign))))) {
                      stop(paste("Term(s)", paste(Terms[noMatch],
                        collapse = ", "), "not matched"))
                    }
                  }
                  else {
                    stop("Terms can only be integers or characters")
                  }
                }
                dDF <- unique(object$fixDF$terms[Terms])
                if (length(dDF) > 1) {
                  # stop("Terms must all have the same denominator DF")
                  warning("Terms do not all have the same denominator DF -- using the minimum")
                  dDF <- min(dDF)
                }
                lab <- paste("F-test for:", paste(names(assign[Terms]),
                  collapse = ", "), "\n")
                L <- diag(nX)[unlist(assign[Terms]), , drop = FALSE]
            }
            else {
                L <- as.matrix(L)
                if (ncol(L) == 1)
                  L <- t(L)
                nrowL <- nrow(L)
                ncolL <- ncol(L)
                if (ncol(L) > nX) {
                  stop(paste("L must have at most", nX, "columns"))
                }
                dmsL1 <- rownames(L)
                L0 <- array(0, c(nrowL, nX), list(NULL, names(object$fixDF$X)))
                if (is.null(dmsL2 <- colnames(L))) {
                  L0[, 1:ncolL] <- L
                }
                else {
                  if (any(noMatch <- is.na(match(dmsL2, colnames(L0))))) {
                    stop(paste("Effects", paste(dmsL2[noMatch],
                      collapse = ", "), "not matched"))
                  }
                  L0[, dmsL2] <- L
                }
                L <- L0[noZeroRowL <- as.logical((L0 != 0) %*%
                  rep(1, nX)), , drop = FALSE]
                nrowL <- nrow(L)
                if (is.null(dmsL1)) {
                  dmsL1 <- 1:nrowL
                }
                else {
                  dmsL1 <- dmsL1[noZeroRowL]
                }
                rownames(L) <- dmsL1
                dDF <- unique(object$fixDF$X[noZeroColL <- as.logical(c(rep(1,
                  nrowL) %*% (L != 0)))])
                if (length(dDF) > 1) {
                  ## stop("L may only involve fixed effects with the same denominator DF")
                  warn <- paste( "L involves fixed effects with the different denominator DF:",
                          paste(dDF, collapse=" "), collapse = " ")
                  warning(warn)
                  dDF <- min(dDF)
                }
                lab <- "F-test for linear combination(s)\n"
            }
            nDF <- sum(svd(L)$d > 0)
            c0 <- c(qr.qty(qr(vFix %*% t(L)), c0))[1:nDF]
            Fval <- sum(c0^2)/nDF
            Pval <- 1 - pf(Fval, nDF, dDF)
            aod <- data.frame(nDF, dDF, Fval, Pval)
            names(aod) <- c("numDF", "denDF", "F-value", "p-value")
            attr(aod, "rt") <- rt
            attr(aod, "label") <- lab
            if (!Lmiss) {
                if (nrow(L) > 1)
                  attr(aod, "L") <- L[, noZeroColL, drop = FALSE]
                else attr(aod, "L") <- L[, noZeroColL]
            }
        }
    }
    else {
        ancall <- sys.call()
        ancall$verbose <- ancall$test <- NULL
        object <- list(object, ...)
        termsClass <- unlist(lapply(object, data.class))
        if (!all(match(termsClass, c("gls", "gnls", "lm", "lmList",
            "lme", "nlme", "nlsList", "nls"), 0))) {
            stop(paste("Objects must inherit from classes \"gls\", \"gnls\"",
                "\"lm\",\"lmList\", \"lme\",\"nlme\",\"nlsList\", or \"nls\""))
        }
        resp <- unlist(lapply(object, function(el) deparse(getResponseFormula(el)[[2]])))
        subs <- as.logical(match(resp, resp[1], FALSE))
        if (!all(subs))
            warning(paste("Some fitted objects deleted because",
                "response differs from the first model"))
        if (sum(subs) == 1)
            stop("First model has a different response from the rest")
        object <- object[subs]
        rt <- length(object)
        termsModel <- lapply(object, function(el) formula(el)[-2])
        estMeth <- unlist(lapply(object, function(el) {
            val <- el[["method"]]
            if (is.null(val))
                val <- NA
            val
        }))
        if (length(uEst <- unique(estMeth[!is.na(estMeth)])) >
            1) {
            stop("All fitted objects must have the same estimation method.")
        }
        estMeth[is.na(estMeth)] <- uEst
        REML <- uEst == "REML"
        if (REML) {
            aux <- unlist(lapply(termsModel, function(el) {
                aux <- terms(el)
                val <- paste(sort(attr(aux, "term.labels")),
                  collapse = "&")
                if (attr(aux, "intercept") == 1) {
                  val <- paste(val, "(Intercept)", sep = "&")
                }
                val
            }))
            if (length(unique(aux)) > 1) {
                warning(paste("Fitted objects with different fixed effects.",
                  "REML comparisons are not meaningful."))
            }
        }
        termsCall <- lapply(object, function(el) {
            if (is.null(val <- el$call)) {
                if (is.null(val <- attr(el, "call"))) {
                  stop("Objects must have a \"call\" component or attribute.")
                }
            }
            val
        })
        termsCall <- unlist(lapply(termsCall, function(el) paste(deparse(el),
            collapse = "")))
        aux <- lapply(object, logLik, REML)
        if (length(unique(unlist(lapply(aux, function(el) attr(el,
            "nall"))))) > 1) {
            stop("All fitted objects must use the same number of observations")
        }
        dfModel <- unlist(lapply(aux, function(el) attr(el, "df")))
        logLik <- unlist(lapply(aux, function(el) c(el)))
        AIC <- unlist(lapply(aux, AIC))
        BIC <- unlist(lapply(aux, BIC))
        aod <- data.frame(call = termsCall, Model = (1:rt), df = dfModel,
            AIC = AIC, BIC = BIC, logLik = logLik, check.names = FALSE)
        if (test) {
            ddf <- diff(dfModel)
            if (sum(abs(ddf)) > 0) {
                effects <- rep("", rt)
                for (i in 2:rt) {
                  if (ddf[i - 1] != 0) {
                    effects[i] <- paste(i - 1, i, sep = " vs ")
                  }
                }
                pval <- rep(NA, rt - 1)
                ldf <- as.logical(ddf)
                lratio <- 2 * abs(diff(logLik))
                lratio[!ldf] <- NA
                pval[ldf] <- 1 - pchisq(lratio[ldf], abs(ddf[ldf]))
                aod <- data.frame(aod, Test = effects, L.Ratio = c(NA,
                  lratio), "p-value" = c(NA, pval), check.names = FALSE)
            }
        }
        row.names(aod) <- unlist(lapply(as.list(ancall[-1]),
            deparse))
        attr(aod, "rt") <- rt
        attr(aod, "verbose") <- verbose
    }
    class(aod) <- c("anova.lme", "data.frame")
    aod
}










if( F) {


#### NOT RUN: ####


glmmPQL <-
function (fixed, random, family, data, correlation, weights,
    control, niter = 10, verbose = TRUE, ...)
{
    if (!require("nlme"))
        stop("package 'nlme' is essential")
    if (is.character(family))
        family <- get(family)
    if (is.function(family))
        family <- family()
    if (is.null(family$family)) {
        print(family)
        stop("'family' not recognized")
    }
    m <- mcall <- Call <- match.call()
    nm <- names(m)[-1]
    keep <- is.element(nm, c("weights", "data", "subset", "na.action"))
    for (i in nm[!keep]) m[[i]] <- NULL
    allvars <- if (is.list(random))
        allvars <- c(all.vars(fixed), names(random), unlist(lapply(random,
            function(x) all.vars(formula(x)))))
    else c(all.vars(fixed), all.vars(random))
    Terms <- if (missing(data))
        terms(fixed)
    else terms(fixed, data = data)
    off <- attr(Terms, "offset")
    if (length(off <- attr(Terms, "offset")))
        allvars <- c(allvars, as.character(attr(Terms, "variables"))[off +
            1])
    m$formula <- as.formula(paste("~", paste(allvars, collapse = "+")))
    environment(m$formula) <- environment(fixed)
    m$drop.unused.levels <- TRUE
    m[[1]] <- as.name("model.frame")
    mf <- eval.parent(m)
    off <- model.offset(mf)
    if (is.null(off))
        off <- 0
    w <- model.weights(mf)
    if (is.null(w))
        w <- rep(1, nrow(mf))
    mf$wts <- w
    fit0 <- glm(formula = fixed, family = family, data = mf,
        weights = wts, ...)
    w <- fit0$prior.weights
    eta <- fit0$linear.predictor
    zz <- eta + fit0$residuals - off
    wz <- fit0$weights
    fam <- family
    nm <- names(mcall)[-1]
    keep <- is.element(nm, c("fixed", "random", "data", "subset",
        "na.action", "control"))
    for (i in nm[!keep]) mcall[[i]] <- NULL
    fixed[[2]] <- quote(zz)
    mcall[["fixed"]] <- fixed
    mcall[[1]] <- as.name("lme")
    mcall$random <- random
    mcall$method <- "ML"
    if (!missing(correlation))
        mcall$correlation <- correlation
    mcall$weights <- quote(varFixed(~invwt))
    mf$zz <- zz
    mf$invwt <- 1/wz
    mcall$data <- mf
    for (i in 1:niter) {
        if (verbose) {
            cat("iteration", i, "\n")
            print(names(mcall))
        }
        fit <- eval(mcall)
        etaold <- eta
        eta <- fitted(fit) + off
        if (sum((eta - etaold)^2) < 1e-06 * sum(eta^2))
            break
        mu <- fam$linkinv(eta)
        mu.eta.val <- fam$mu.eta(eta)
        mf$zz <- eta + (fit0$y - mu)/mu.eta.val - off
        wz <- w * mu.eta.val^2/fam$variance(mu)
        mf$invwt <- 1/wz
        mcall$data <- mf
    }
    attributes(fit$logLik) <- NULL
    fit$call <- Call
    fit$family <- family
    oldClass(fit) <- c("glmmPQL", oldClass(fit))
    fit
}

#log <-function (x, base = exp(1)) {
#    x <- pmax(x,.00000001)
#    if (missing(base)) .Internal(log(x)) else .Internal(log(x, base))
#}




logLik.reStruct  <-
function (object, conLin, ...)
{
    if (any(!is.finite(conLin$Xy)))
        return(-1000000)
    .C("mixed_loglik", as.double(conLin$Xy), as.integer(unlist(conLin$dims)),
        as.double(pdFactor(object)), as.integer(attr(object,
            "settings")), loglik = double(1), double(1), PACKAGE = "nlme")$loglik
}

}

"

"