###########################################################################
pletter.list <- function( plet, condense = T )
  {
    plet <- as.character( plet )
    nprod <- length( plet )
    x <- list( )
    i <- 1
    li <- letters[i]
    sprod <- seq( nprod )
    pplet <- paste( plet, collapse = "" )
    anychar <- ".*"
    if( is.null( version$language ))
      {
        if( pmatch( "win", casefold( version$platform )))
          anychar <- ".*"
      }
    else
      {
        if( version$language != "R" )
          anychar <- ".*"
      }
    
    while( length( grep( paste( anychar, "[", li, "-z]", anychar, sep = "" ), pplet )) > 0 )
      {
        tmp <- grep( paste( anychar, li, anychar, sep = "" ), plet )
        if( length( tmp ))
          x[[li]] <- tmp
        i <- i + 1
        li <- letters[i]
      }
    if( length( x ) == 0 )
      x <- list( sprod )
    # condense out singletons and subsets
    if( condense )
      {
        # singletons
        for( i in seq( length( x ), 1 ))
          if( length( x[[i]] ) <= 1 )
            x[[i]] <- NULL
        nx <- length( x )
        # subsets
        drop <- NULL
        if( nx > 1 ) for( i in seq( nx - 1 ))
          {
            tmp <- match( x[[i]], x[[i+1]] )
            if( any( is.na( tmp )))
              {
                tmp <- match( x[[i+1]], x[[i]] )
                if( !any( is.na( tmp )))
                  drop <- c( drop, i + 1 )
              }
            else
              drop <- c( drop, i )
          }
        if( !is.null( drop ))
          for( i in rev( drop ))
            x[[i]] <- NULL
      }
    x
  }
###########################################################################
pletter <- function( pdiff, level = .05, let = letters )
  {
    n <- nrow( pdiff )
    pdiff[ seq( from = 1, to = n*n, by = n+1 )] <- 1
    for( i in seq( n - 1 ))
      pdiff[ i, seq( i + 1, n ) ] <- 0

    pdiff <- pdiff > level
    pdiff[ is.na( pdiff ) ] <- T
    prev <- pdiff[,1]

    code <- character( n )
    code[ prev ] <- let[1]
    j <- 1
    for( i in seq( 2, n )) if( any( (!prev) & pdiff[,i] ))
      {
        j <- j + 1
        prev <- pdiff[,i]
        code[prev] <- paste( code[prev], let[j], sep = "" )
      }
    names(code) <- dimnames( pdiff )[[1]]
    code
  }
###########################################################################
pdiff <- function( fit, lsm = lsmean( fit ),
                   pred = outer( lsm$pred[order.pred], lsm$pred[order.pred], "-" ),
                   se = outer( lsm$se[order.pred], lsm$se[order.pred],
                     function( x, y ) sqrt( x^2+y^2 )),
                   rdf = df.resid( fit ))
  {
    order.pred <- order( -lsm$pred )
    tmp <- - abs( pred / se )
    2 * pt( tmp, rdf )
  }
###########################################################################
pvalue.ellipse <- function( plet, lsmeans, se, rdf,
                           level = .05, scale = 0.3, full.width = T )
  {
    pset <- pletter.list( plet )
    level <- 1 - level / 2
    ni <- length( plet )
    ns <- length( pset )
    if( ns == 0 )
      return( matrix( NA, 2, 0 ))
    elstuff <- matrix( NA, 2, ns, dimnames = list( c("x","sx"), 1:ns ))
    if( length( rdf ) == 1 )
      rdf <- array( rdf, ni )

    pnames <- names( lsmeans )
    if( is.null( pnames ))
      pnames <- as.character( seq( lsmeans ))

    tmp <- function( r, se, rdf, level, scale, full.width,
                    tmpfn = function( x, y ) sqrt( x^2 + y^2 ))
      {
        if( length( se ) > 1 )
          {
            se <- outer( se, se, tmpfn )
            se <- max( se[ row( se ) > col( se ) ] )
            radius <- scale * se + diff( r ) / 2
          }
        else
          {
            radius <- se <- 0
          }
        rdf <- min( rdf )

        if( full.width )
          c( mean( r ), max( radius, se * qt( level, rdf ) / 2 ))
        else
          c( mean( r ), max( radius ))
      }
    for( i in seq( ns ))
      {
        same <- pset[[i]]
        elstuff[,i] <- tmp( range( lsmeans[same] ), se[same], rdf[same],
                           level, scale, full.width )
      }
    as.matrix( elstuff )
  }
arrow <- function(x,y,head=3,tip=.025,width=tip/2,fine=.005,...)
{
lines(x,y,...)
n <- length(x)
u <- par()$usr
diag <- (u[1]-u[2])^2+(u[3]-u[4])^2

if (head%%2)
	dotip(x[1], x[2]-x[1], y[1], y[2]-y[1], diag, tip, width,...)
if (head>1)
	dotip(x[n], x[n-1]-x[n], y[n], y[n-1]-y[n], diag, tip, width,...)
}
dotip <- function(x,dx,y,dy,diag,tip,width,...)
{
	dd <- sqrt((dx^2+dy^2)/diag)
	a <- tip / dd
	b <- width / (2*dd)
	tx <- x + a*dx + c(1,-1)*b*dy
	ty <- y + a*dy + c(-1,1)*b*dx
	polygon(c(x,tx),c(y,ty),...)
}

arc.arrow <- function(diam,head=3,fine=.005,x=diam,y=.5+2*arc,
	arc=(.5-diam)/4,...)
{
	tmp <- circle(x,y,diam,show=F,fine=fine)[[1]]
	tmp <- tmp[tmp$y>=y-arc&tmp$y<=y+arc&tmp$x<x,]
	arrow(tmp$x,tmp$y,head,...)
	invisible(tmp)
}
path3 <- function(cell=rep("",4),xlab="",xlim=c(-.1,1),diam=.15,ylab="",...)
{
	tmp <- diam / sqrt(5)

	plot(xlim,0:1,type="n",axes=F,xlab="",ylab=ylab,...)
	mtext(xlab,1)

if(nchar(cell[1])){
	circle(1-diam,.5,diam,cen=0)
	text(1-diam,.5,cell[1])
}
if(nchar(cell[2])){
	circle(diam,1-diam,diam,cen=0)
	text(diam,1-diam,cell[2])
	arrow(c(diam+2*tmp,1-diam-2*tmp),c(1-diam-tmp,.5+tmp),2)
}
if(nchar(cell[3])){
	circle(diam,.5,diam,cen=0)
	text(diam,.5,cell[3])
	arrow(c(2*diam,1-2*diam),c(.5,.5),2)
}
if(nchar(cell[4])){
	circle(diam,diam,diam,cen=0)
	text(diam,diam,cell[4])
	arrow(c(diam+2*tmp,1-diam-2*tmp),c(diam+tmp,.5-tmp),2)
}
}
#bekk <- read.table("../src/bekk.dat",header=T)
#bekk$lab <- factor(bekk$lab)
#bekk$mat <- factor(bekk$mat)
# see help(trans.plot)
box.fig <- function(fit,xpos=1,width=.2,jitter=width)
{
	fit <- boxplot(fit,plot=F)
  if( is.null( fit$stats ))
    fit <- fit[[1]]
        tmp <- fit$stats
	polygon(xpos+width*c(-1,1,1,-1),rep(tmp[c(2,4)],rep(2,2)))
	lines(xpos+width*c(-1,1),rep(tmp[3],2))
	lines(xpos+width*c(-1,1),rep(tmp[1],2),lty=2)
	lines(xpos+width*c(-1,1),rep(tmp[5],2),lty=2)
	lines(rep(xpos,2),tmp[1:2],lty=2)
	lines(rep(xpos,2),tmp[4:5],lty=2)
	if (!is.null(fit$out) & length( fit$out ))
		points(xpos+runif(length(fit$out),-jitter,jitter),fit$out)
	fit
}
#choose <- function(n, k)
#{
#        if(k > n | k < 0)
#                return(0)
#        k1 <- min(n - k, k)
#        j <- n - k1
#        coef <- 1
#        for(i in 1:k1) {
#                j <- j + 1
#                coef <- (coef * j)/i
#        }
#        coef
#}
fact <- function(k)
{
        k[k == 0] <- 1
        f <- array(1, length(k))
        for(i in 1:length(f))
                for(j in 1:k[i])
                        f[i] <- f[i] * j
        f
}
perms <- function(x,k)
{
	n <- length(x)
	if (n == 1) n <- x
	if (n < k)
		stop("length of x must be greater than k")
	a <- matrix(c(rep(F,n-k),rep(T,k)),n,1)
	if (k < n & k > 0) for (i in (n-k-1):0) {
		b <- c(rep(F,i),T)
		d <- perms(n-i-1,k-1)
		a <- cbind(a,rbind(matrix(b,i+1,ncol(d)),d))
	}
	if (length(x) == 1)
		a
	else
		matrix(x,n,ncol(a))[a]
}
ci.plot <- function(object, ...) UseMethod("ci.plot")

ci.plot.default <- function(x.factor,response,group,
   xlab=xlabs,ylab=ylabs,bar.plot="ci",...)
{
   xlabs <- deparse(substitute(x.factor))
   ylabs <- deparse(substitute(response))

   if (!is.factor(x.factor))
      x.factor <- factor(x.factor)
   if( missing(group) ) {
      data <- data.frame(response,x.factor)
      names(data) <- c(ylabs,xlabs)
      formula.y <- formula( paste( ylabs, xlabs, sep="~" ))
      fit <- lm(formula.y,data)
   }
   else {
      glabs <- deparse(substitute(group))
      if (!is.factor(group))
         group <- factor(group)
      data <- data.frame(response,x.factor,group)
      names(data) <- c(ylabs,xlabs,glabs)
      formula.y <- formula( paste( ylabs, ".^2", sep="~" ))
      fit <- lm(formula.y,data)
   }
   invisible(int.plot.lm(fit,data,xlab=xlab,ylab=ylab,bar.plot=bar.plot,...))
}

ci.plot.lm <- function(fit,...,bar.plot="ci")
invisible(int.plot.lm(fit,...,bar.plot=bar.plot))
ci.width <- function(fit,data=eval(fit$call$data),
   factors=all.factors(fit,data),lsm=lsmean(fit,data,factors,...),
   rdf=df.resid(fit),level=.05,crit=qt(1-level/2,rdf),...)
{
   if (level > 0 & level < 1) crit * lsm$se else 0
}
dotna <- function(data)
{
   tmp <- function(x)
   {
      dot <- x == "."
      if(any(dot)) {
         tmpopt <- options(warn=-1)
         on.exit(options(tmpopt))
         x <- as.character(x)
         x[dot] <- NA
         if(!any(is.na(as.numeric(x[!dot]))))
            x <- as.numeric(x)
         else
            x <- factor(x)
      }
      x
   }
   data <- as.data.frame(data)
   for (i in seq(length(data)))
      data[[i]] <- tmp(data[[i]])
   data
}
effect.plot <- function(fit,effect=effect.fit(fit,adjust=adjust,...),
   adjust=T,boxplot.min=10,jitter=.1,width=jitter,xlim=xlims,ylim=ylims,xaxt="i",
   xlab="Terms",ylab=ylabs,...)
{
   nterms <- length(effect)
   # eliminate degenerate effects
   for (i in 1:nterms)
      if (all(is.na(effect[[i]])))
         effect[[i]] <- NULL
   nterms <- length(effect)

   if (adjust)
      ylabs <- "MS Adjusted Effects"
   else
      ylabs <- "Effects"

   xlims <- c(1-max(jitter,width),nterms+max(jitter,width))
   ylims <- range(unlist(effect))
   plot(xlim,ylim,type="n",xaxt="n",xlab=xlab,ylab=ylab,...)
   if( xaxt != "n" )
      axis(1,1:nterms,names(effect))
   abline(h=0,lty=2)
   par(lty=1)

   for (i in 1:nterms) {
      neff <- length(effect[[i]])
      if (neff < boxplot.min) {
         name1 <- names(effect[[i]])
         if(is.null(name1))
            points(i+runif(neff,-jitter,jitter),effect[[i]])
         else
            text(i+runif(neff,-jitter,jitter),effect[[i]],labels=name1)
            
      }
      else
         box.fig(effect[[i]],i,width,jitter)
   }
   invisible(effect)
}


# effects contrasts

effect.fit <- function(object, ...) UseMethod("effect.fit")

effect.fit.default <- function(fit,data=eval(fit$call$data),
   factors=all.factors(fit,data),adjust=T,influence=F,
   terms.fit=as.data.frame(predict(fit,type="terms")),...)
{
   if (adjust)
      tmp <- function(x,f=NULL) {
         if (is.null(f)) {
            nx <- tapply(x,x,length)
            ux <- as.numeric(names(nx))
            ff <- NULL
         }
         else {
            nx <- tapply(x,f,length)
            ux <- tapply(x,f,mean)
            ff <- tapply(f,f,unique)
            if (is.factor(f))
               ff <- levels(f)[ff]
         }
         lx <- length(ux)
         x <- ux * sqrt(nx*lx/(lx-1))
         list(x=x,f=ff)
      }
   else
      tmp <- function(x,f=NULL) {
         if(is.null(f)) {
            x <- unique(x)
            ff <- NULL
         }
         else {
            x <- tapply(x,f,mean)
            ff <- tapply(f,f,unique)
            if (is.factor(f))
               ff <- levels(f)[ff]
         }
         list(x=x,f=ff)
      }

   se <- .01*diff(range(terms.fit))
   eff <- list()
   for (i in names(terms.fit)) {
      tmpterm <- tmp(precision(terms.fit[[i]],se),data[[i]])
      eff[[i]] <- tmpterm$x
      if (!is.na(match(i,factors)))
         names(eff[[i]]) <- as.character(tmpterm$f)
      else
        {
          # for some reason this does not work in R
         names(eff[[i]]) <- NULL
         names( eff[[i]] ) <- rep( "o", length( eff[[i]] ))
       }
   }

   tmp <- resid(fit)
   if (influence)
      eff[["resid"]] <- tmp / sqrt(1-lm.influence(fit)$hat)
   else {
      if (adjust)
         eff[["resid"]] <- tmp * sqrt(length(tmp)/df.resid(fit))
      else
         eff[["resid"]] <- tmp
   }
   eff
}

effect.fit.listof <- function(fit,data=eval(attr(fit,"call")$data),
   factors=all.factors(fit,data),terms.fit=fit.proj,stratum,...)
{
   fit.proj <- as.data.frame(proj(fit)[[stratum]])
   fit.proj$Residuals <- NULL
   effect.fit.default(fit[[stratum]],data=data,factors=factors,
      terms.fit=terms.fit,...)
}
circle <- function(x=0,y=0,radius=1,sx=radius,sy=radius,...)
   invisible(ellipse(x,y,sx=sx,sy=sy,...))

ellipse <- function(x=0,y=0,sx=1,sy=1,rho=0,fine=.005,center=.1,show=T,add=T,
   xlab="",ylab="",bty="n",xaxt="n",yaxt="n",type="n",lty=1,pty="s",err=-1,
   ...) {
   if(min(sx) < 0 | min(sy) < 0 | max(abs(rho)) > 1)
      stop("sx,sy or rho out of bounds")
   if(fine <= 0 | fine > 1)
      stop("fine out of bounds")
   n <- max(length(x),length(y),length(sx),length(sy),length(rho))
   x <- array(x,n)
   y <- array(y,n)
   lty <- array(lty,n)
   sx <- array(sx,n)
   sy <- array(sy,n)
   rho <- array(rho,n)
   center <- array(center,n)
   prin1 <- sqrt((1+rho)/2)
   prin2 <- sqrt((1-rho)/2)
   xarc <- sin(pi * seq(0,.5,by=4*fine))
   yarc <- c(rev(xarc),-xarc)
   xarc <- c(xarc,yarc,-rev(xarc))
   yarc <- c(yarc,rev(yarc))
   arc <- list()
   for (i in 1:n)
      arc[[i]] <- data.frame(x=x[i]+sx[i]*(prin1[i]*xarc+prin2[i]*yarc),
         y=y[i]+sy[i]*(prin1[i]*xarc-prin2[i]*yarc))

   if(show) {
      if(!add)
        {
          tmp <- options( warn = -1 )
          plot(c(min(x-sx),max(x+sx)),c(min(y-sy),max(y+sy)),
               bty=bty,xaxt=xaxt,yaxt=yaxt,type=type,xlab=xlab,ylab=ylab,pty=pty,
               ...)
          options( tmp )
       }
      for (i in 1:n) {
         lines(arc[[i]]$x,arc[[i]]$y,lty=lty[i],err=err,...)
         if(center[i] > 0) {
            lines(x[i]+center[i]*c(-1,1)*sx[i]*prin1[i],
               y[i]+center[i]*c(-1,1)*sy[i]*prin1[i],err=err)
            lines(x[i]+center[i]*c(-1,1)*sx[i]*prin2[i],
               y[i]+center[i]*c(1,-1)*sy[i]*prin2[i],err=err)
         }
      }
   }
   invisible(arc)
}
norm.ellipse <- function(x, y, group, homosked=T, level=.05,
   radius=sqrt(qchisq(1-level,2)),...) {
   mx <- mean(x)
   my <- mean(y)
   sx <- std.dev(x)
   sy <- std.dev(y)
   rho <- if(length(x) > 1) cor(x,y,trim=.00001) else 1
   ellipse(mx,my,sx*radius,sy*radius,rho,...)
   invisible(list(x=c(mx,sx),y=c(my,sy),rho=rho))
}
df.resid <- function(object, ...) UseMethod("df.resid")

df.resid.default <- function(...)
{
   x <- list(...)
   if(length(x) == 1 && is.list(x[[1]]))
      x <- x[[1]]
   x <- tapply(x[[1]],x,length) - 1
   x[is.na(x)] <- 0
   sum(x)
}

df.resid.lm <- function(object,resid=residuals(object))
{
   rdf <- object$df.resid
   if(is.null(rdf)) {
      p <- object$rank
      if(is.null(p))
         p <- sum(!is.na(coefficients(object)))
      rdf <- length(resid) - p
   }
   rdf
}
se.coef <- function(object,summary=summary.lm(object))
{
   summary$coefficients[,2]
}

cor.coef <- function(object,summary=summary.lm(object))
{
   summary$correlation
}

cov.coef <- function(object,summary=summary.lm(object))
{
   summary$cov.unscaled
}
std.dev <- function(object)
{
   # adapted from summary.lm
   sd <- summary(object)$sigma
   if(is.null(sd))
   {
      resid <- residuals(object)
      if(is.null(resid))
         sd <- sqrt(var(object))
      else
      {
         rdf <- df.resid(object)
         if(rdf > 0)
            sd <- sqrt(sum(resid^2)/rdf)
         else
            sd <- NA
      }
   }
   if(!is.na(sd)) sd
      else 0
}
sample.size <- function(object, ...) UseMethod("sample.size")

sample.size.default <- function(...)
{
   x <- list(...)
   if(length(x) == 1 && is.list(x[[1]]))
      x <- x[[1]]
   x <- tapply(x[[1]],x,length)
   x[is.na(x)] <- 0
   x
}

sample.size.lm <- function(fit,data=eval(fit$call$data),
   factors=all.factors(fit,data),na.action=na.omit)
{
   n <- length(factors)
   if (n < 1)
      length(data[[model.response(fit)]])
   else {
      replications(formula(paste("~",paste(factors,collapse=":"),sep="")),
         data, na.action = na.action)[[1]]
   }
}
cor.coef <- function(object)
{
   summary.lm(object)$correlation
}
harmonic.mean <- function(x)
{
	length(x)/sum(1/x)
}
expt <- function(cell,effects,inter,grand=0,design=c("rcbd","crd","bibd"),
	reps=4,factors=c("A","B"),tukey=F,sd=c(1,.5),anova=T,bibd)
{
	if(missing(cell)){
		if(missing(effects))
			stop("must specify cell means or effects")
		if(length(effects)!=2)
			stop("only two factor effects allowed")
		marg <- expand.grid(effects)
		cell <- grand + apply(marg,1,sum)
		if(!missing(inter))
			cell <- cell + array(c(inter),length(cell))
		else {if(tukey!= 0)
			cell <- cell + tukey*apply(marg,1, prod)
		}
		cell <- array(cell,unlist(lapply(effects,length)))
		tmp <- function(x){
			n <- names(x)
			if(is.null(n))
				letters[seq(x)]
			else
				n
		}
		dimnames(cell) <- lapply(effects,tmp)
		factors <- names(effects)
	}
	nt <- length(cell)
	if (!missing(bibd)) {
		if (bibd >= nt | bibd < 2)
			stop(paste("bibd =",bibd,"must be between 2 and",nt))
		design <- "bibd"
	}

	sd <- array(c(sd,0),2)
	des <- pmatch(design[1],c("rcbd","crd","bibd"))
	if (des == 1 & reps > 1) {
		# rcbd
		y <- rep(rnorm(reps,0,sd[1]),rep(nt,reps)) + rep(cell,reps) +
			 rnorm(reps*nt,0,sd[2])
		block <- "reps +"
	}
	else if (des == 3) {
		# bibd
		numib <- choose(nt,bibd)
		y <- rep(perms(cell,bibd),reps)
		reps <- reps * numib
		nt <- bibd
		y <- y + rep(rnorm(reps,0,sd[1]),rep(nt,reps)) +
			rnorm(reps*nt,0,sd[2])
		block <- "reps +"
	}
	else {
		# crd
		y <- rep(cell,reps) + rnorm(reps*nt,0,sqrt(sd[1]^2+sd[2]^2))
		block <- ""
	}
browser()
	data <- data.frame(reps=factor(rep(1:reps,rep(nt,reps))))
	tmp <- list(ncol(cell),rep(nrow(cell),ncol(cell)))
browser()
	for (i in 1:2) {
		cellnames <- dimnames(cell)[[i]]
		if (is.null(cellnames))
			cellnames <- as.character(seq(dim(cell)[i]))
		data[[factors[i]]] <- factor(rep(rep(cellnames,tmp[[i]]),reps))
	}
browser()
	data$y <- as.numeric(y)
	if(anova) {
		red <- as.formula(paste("y ~",block,factors[1],"+",factors[2]))
		full <- as.formula(paste("y ~",block,factors[1],"*",factors[2]))
		fullfit <- aov(full,data)
		print(summary(fullfit))
		list(full=fullfit,reduced=aov(red,data),data=data)
	}
	else
		data
}
all.factors <- function (fit,data=eval(fit$call$data))
{
   predictors <- all.names(formula(fit),unique=T,functions=F)[-1]
   factors <- unlist(lapply(data,is.factor))[predictors]
   names(factors[factors])
}

all.covars <- all.covariates <- function (fit,data=eval(fit$call$data))
{
   predictors <- all.names(formula(fit),unique=T,functions=F)[-1]
   factors <- unlist(lapply(data,is.factor))[predictors]
   names(factors[!factors])
}

all.predictors <- function (fit,data=eval(fit$call$data))
{
   all.names(formula(fit),unique=T,functions=F)[-1]
}
#model.response <- function (fit,data=eval(fit$call$data))
#{
#   all.names(formula(fit),functions=F)[1]
#}
factor.labels <- function(object,data) {
   p <- pmatch(names(data),labels(object))
   labels(object)[p[!is.na(p)]]
}
unfactor <- function(x.factor,class="factor")
{# adapted from interaction.plot
   ok <- inherits(x.factor,substitute(class))
   if (ok)
   {
     xval <- as.character( x.factor )
     tmp <- options( warn = -1 )
     nval <- as.numeric( xval )
     options( tmp )
     ok <- !all( is.na( nval ))
     if( ok )
       xval <- nval
   }
   if (ok)
      xval
   else
      unclass(x.factor)
}
get.list <- function(data,factors,sortlist=F)
{
  r <- list()
  for( i in factors )
    r[[i]] <- data[[i]]
   if(sortlist)
   {
      # sort em!
   }
   r
}
ghostview <- function(file, background = T)
{
	if (background)
		background <- "&"
	else
		background <- ""
        invisible(unix(paste("ghostview", file, background), output = F))
}
lsd.plot <- function(object, ...) UseMethod("int.plot")
int.plot <- function(object, ...) UseMethod("int.plot")

int.plot.default <- function(x.factor,response,group,
   xlab=xlabs,ylab=ylabs,...)
{
   xlabs <- deparse(substitute(x.factor))
   ylabs <- deparse(substitute(response))

   if (!is.factor(x.factor))
      x.factor <- factor(x.factor)
   if( missing(group) ) {
      data <- data.frame(response,x.factor)
      names(data) <- c(ylabs,xlabs)
      formula.y <- formula( paste( ylabs, xlabs, sep="~" ))
      fit <- lm(formula.y,data)
   }
   else {
      glabs <- deparse(substitute(group))
      if (!is.factor(group))
         group <- factor(group)
      data <- data.frame(response,x.factor,group)
      names(data) <- c(ylabs,xlabs,glabs)
      formula.y <- formula( paste( ylabs, ".^2", sep="~" ))
      fit <- lm(formula.y,data)
   }
   invisible(int.plot.lm(fit,data,xlab=xlab,ylab=ylab,...))
}

int.plot.lm <- function(fit,data=eval(fit$call$data),
                        factors=all.factors(fit,data),
                        lsm=lsmean(fit,data,factors),
                        ci=ci.width(fit,data,factors,lsm,rdf,level),
                        offset=0.01,bar.plot="lsd",
                        xlim=xlims,ylim=ylims,type="b",
                        width=.02,sort.mean=F,edge=wi,
                        xlab=xlabs,ylab=as.character(formula(fit)[[2]]),
                        xaxt="s",level=.05,col=1,rdf = df.resid(fit),
                        fine=.01,...)
{
  bar.plot <- c("lsd","ci","test","ellipse","none")[
     pmatch( bar.plot, c("lsds","cis","tests","ellipses"), nomatch = 5 ) ]
  
  # make room for confidence intervals or test intervals
  if( pmatch( bar.plot, c("ci","test"), nomatch = 0 ))
    {
      wi <- .02
      if ( bar.plot == "test" )
        ci <- ci / sqrt(2)
      ylims <- c(min(lsm$pred-ci),max(lsm$pred+ci))
    }
  else {
    wi <- 0
    ylims <- range(lsm$pred)
  }

   # set up plotting groups if any
   xlabs <- factors[1]
   lf <- length(factors)
   if (lf == 1) { # only 1 factor
      group <- NULL
      x.factor <- lsm[[factors]]
   }
   else {
      if (lf > 2) { # x is 1st factor, all others have traces
         group <- interaction( get.list( lsm, factors[-1] ), drop = T )
         x.factor <- lsm[[factors[1]]]
      }
      else { # list of two vectors of factor names
         group <- interaction( get.list( lsm, factors[[2]] ), drop = T )
         x.factor <- interaction( get.list( lsm, factors[[1]] ), drop = T )
         xlabs <- paste(factors[[1]],collapse=".")
   }
}
   # set up horizontal axis -- factor levels or sorted by mean
   if (sort.mean) {
      if (is.null(group))
         x <- order(order(-lsm$pred))
      else {
         x <- order(order(-tapply(lsm$pred,x.factor,mean)))
         names(x) <- levels(x.factor)
         x <- x[x.factor]
      }
   }
   else
      x <- unfactor(x.factor)

  # set up ellipses; adjust ylims for ellipses
  if( pmatch( bar.plot, "ellipse", nomatch = 0 ))
    {
      width <- 3 * width
      plet <- pletter( pdiff( fit, lsm, rdf = rdf ), level = level )
      o <- order( order( - lsm$pred ))
      plet <- plet[o]
      
      xx <- unfactor(x,"ordered")
      ux <- unique( xx )
      ellipx <- list()
      for( i in seq( length( ux )))
        {
          ii <- xx == ux[i]
          ellipx[[i]] <- tmp <- pvalue.ellipse( plet[ii], lsm$pred[ii],
                                               lsm$se[ii], rdf,
                                               level = level )
          if( length( tmp ))
            ylims <- range( ylims, tmp[1,] - tmp[2,], tmp[1,] + tmp[2,] )
        }
    }
  
   # horizontal limits
   xlims <- range(x)
   if (width+edge > 0)
      xlims <- xlims * (1+width+edge) - rev(xlims) * (width+edge)

   # multiple group plot

   # KLUDGE: mplot is on its way out; also, ... causes problems with lsd.bar options
   tmpopt <- options( warn = -1 )
   mplot( x, lsm$pred, group, xlim = xlim, ylim = ylim, xaxt = "n", type = type,
      xlab=xlab,ylab=ylab,col=col,...)
   options( tmpopt )
   # add factor levels on horizontal axis
   if (xaxt != "n")
      mtext(as.character(x.factor),1,at=x)

   # intervals at each point or one LSD bar or none?
  switch( bar.plot,
         none =
         {
           invisible( )
         },
         lsd =
         {
           invisible( lsd.bar.lm( fit, data, factors, lsm = lsm,
                                 width = width, level = level,
                                 rdf = rdf, ... ))
         },
         test =, ci =
         {
           x <- unfactor(x,"ordered")
           group <- if (length(factors) == 2) lsm[[factors[2]]] else NULL
           if (!is.null(group) & offset > 0)
             x <- x + offset * diff(par("usr"))[1] *
               (unclass(group) - (1+length(levels(group)))/2)
           
           invisible(se.bar(x,lsm$pred,2*ci,width=width,cap="",...))
         },
         ellipse =
         {
           tmpfn <- function( y, origin, width )
             {
               if( !is.na( y[2] ) & y[2] > 0 )
                 ellipse( origin[1], origin[2] + y[1], width, y[2], center = 0,
                         fine = fine )
               return(0)
             }
           width <- diff( range( xx )) * width
           ux <- unique( xx )
           for( i in seq( length( ux )))
             {
               tmp <- ellipx[[i]]
               if( length( tmp ) & !is.na( tmp[1] ))
                 apply( tmp, 2, tmpfn, c(ux[i],0), width )
             }
           invisible(  )
         })
}
lpr <- function(file, printer = unix("printenv PRINTER"))
{
       invisible(unix(paste("lpr -P",printer,file),output=F))
}

lpq <- function(printer = unix("printenv PRINTER"))
{
        invisible(unix(paste("lpq -P", printer),output=F))
}
enscript <- function(file,options="-2r")
{
	invisible(unix(paste("enscript",options,file),output=F))
}
lsd.bar <- function(object, ...) UseMethod("lsd.bar")

lsd.bar.default <- function(x.factor,response,group,...)
{
   xlabs <- deparse(substitute(x.factor))
   ylabs <- deparse(substitute(response))

   if( missing(group) ) {
      data <- data.frame(response,x.factor)
      names(data) <- c(ylabs,xlabs)
      formula.y <- formula( paste( ylabs, xlabs, sep="~" ))
      fit <- lm(formula.y,data)
   }
   else {
      glabs <- deparse(substitute(group))
      if (!is.factor(group))
         group <- factor(group)
      data <- data.frame(response,x.factor,group)
      names(data) <- c(ylabs,xlabs,glabs)
      formula.y <- formula( paste( ylabs, ".^2", sep="~" ))
      fit <- lm(formula.y,data)
   }
   invisible(lsd.bar.lm(fit,data,...))
}

lsd.bar.lm <- function(fit,data=eval(fit$call$data),
   factors=all.factors(fit,data),lsm=lsmean(fit,data,factors),
   xpos=xposn,ypos=yposn,rdf=df.resid(fit),
   level=.05,crit=qt(1-level/2,rdf),
   cap=paste(100*level,"%\nLSD",sep=""),mod=mods,tol=1e-5,...)
{
   usr <- par("usr")
   xlen <- usr[2] - usr[1]
   xposn <- usr[2] - .2 * xlen
   yposn <- (2 * usr[3] + usr[4]) / 3

   lsd <- sqrt(mean(lsm$se[lsm$se>0]^2))
   mods <- if (any(abs(lsd-lsm$se)>tol*diff(range(lsm$pred)) & lsm$se>0))
      "*" else ""
   lsd <- crit * sqrt(2) * lsd

   se.bar(xpos,ypos,lsd,mod=mod,cap=cap,...)
   invisible(list(lsd=lsd,rdf=rdf,level=level))
}
lsmean <- function(object, ...) UseMethod("lsmean")

lsmean.default <- function(y,...,factors=all.factors(fit,data)[1:2],
   effects=F,se.fit=T,adjust.covar=T)
{
   data <- eval(substitute(data.frame(y,...)))
   y <- deparse(substitute(y))
   fit <- lm(y~.,data,singular.ok=T)
   lsmean.lm(fit,data,factors,effects=effects,se.fit=se.fit,adjust.covar=adjust.covar)
}

lsmean.lm <- function(fit,data=eval(fit$call$data),
   factors=names(predictors[predictors]),expr=formula(fit),
   contrast=fit$contrasts,effects=F,se.fit=T,adjust.covar=T)
{
   expr.labels <- attr(terms(expr),"term.labels")
   # get predictors from right side of equation
   predictors <- unlist(lapply(data, is.factor))[
      all.names(expr, unique = T, functions = F)[-1]]
   predictors <- predictors[!is.na(charmatch(names(predictors),expr.labels))]

   # set up factors of interest, auxiliary factors and covariates
   covars <- names(predictors[!predictors])
   aux.factors <- predictors[predictors]
   if (length(aux.factors)) {
      pm <- pmatch(unlist(factors),names(aux.factors))
      pm <- pm[!is.na(pm)]
   }
   else
      pm <- numeric(0)
   if(length(pm)==0) {
      factors <- NULL
      aux.factors <- names(aux.factors)
   }
   else {
      factors <- names(aux.factors[pm])
      aux.factors <- names(aux.factors[-pm])

      # check contrasts if using effects option
      if (effects) {
         for (i in factors) if (any(abs(apply(contrast[[i]],2,sum)>1e-12)))
            stop(paste("effects=T not implemented for factors",
               "without sum-to-zero contrasts (",i,")"))
      }
   }

   # separate auxiliary factors with respect to sum-to-zero contrasts
   cont.factors <- NULL
   for (i in aux.factors) if (any(apply(contrasts(data[[i]]),2,sum)))
      cont.factors[i] <- i
   if(length(cont.factors)) if(!is.null(aux.factors))
      aux.factors <- aux.factors[-pmatch(aux.factors,cont.factors,nomatch=0)]

   # ignore aux factors whose contrasts sum to zero
   fac <- attr(terms(expr),"factors")
   tmp <- c(factors,covars,cont.factors)
   term.factors <- fac[tmp,]
   if (length(tmp) > 1)
      term.factors <- apply(term.factors,2,any)
   not.ignore <- 1 - fac[aux.factors,]
   if (length(aux.factors) > 1)
      not.ignore <- apply(not.ignore,2,all)
   if (length(aux.factors) & length(not.ignore))
      term.factors <- term.factors & not.ignore
   term.factors <- expr.labels[term.factors]

   # get indexes for terms used in LS mean fit
   coef <- coefficients(fit)
   coef <- coef[ !is.na( coef ) ]
   # pos <- attr(coef,"assign")
   # kludge! (R does not have assign as attribute for coef
   mat <- model.matrix( fit )
   pos <- attr( mat, "assign" )
   if( !is.list( pos ))
     {
       posnames <- attr( terms( fit ), "term.labels" )
       posnum <- pos
       pos <- list()
       pos[["(Intercept)"]] <- 1

       spos <- seq( length( posnum ))
       for( i in seq( length( posnames )))
         pos[[ posnames[i] ]] <- spos[ i == posnum ]
     }
   index <- if (!effects) pos[["(Intercept)"]]
   for (i in term.factors)
      index <- c(index,pos[[i]])
   # make sure index values are not missing for mat and for coef
   is.coef <- match( dimnames( mat )[[2]], names( coef ), nomatch = 0 )
   mindex <- index[ is.coef[index] > 0 ]
   index <- is.coef[index]
   
   # simplify design: unique factor combinations, covariates at mean level

   # unique factor combinations
   if(is.null(factors)) {
      newdata <- data[1,]
      int.factors <- "mean"
      row.names(newdata) <- int.factors
   }
   else {
      int.factors <- interaction( data[,factors], drop = T )
      design <- !duplicated(int.factors)
      newdata <- data[design,]
   }
   # average covariates to mean value
   if (adjust.covar | is.null(factors))
     {
       for (i in covars)
         newdata[[i]] <- rep( mean( data[[i]], na.rm = T ), nrow( newdata ))
     }
   else
     {
       for (i in covars)
         newdata[[i]] <- c( tapply( data[[i]], int.factors, mean ) )
     }
   if(!is.null(factors)) {
      for (i in rev(factors))
         newdata <- newdata[order(newdata[[i]]),]
      row.names(newdata) <- int.factors <- interaction( newdata[,factors], drop = T )
   }

   # expand matrix for auxiliary factors without sum-to-zero contrasts
   if (!is.null(cont.factors)) {
     factornames <- list()
     for( i in names( data ))
       factornames[[i]] <- levels( data[[i]] )
      cont.levels <- expand.grid(get.list(factornames,cont.factors))
      cont.n <- nrow(cont.levels)
      newdata <- newdata[rep(int.factors,cont.n),]
      tmp <- rep(seq(len=cont.n),rep(length(int.factors),cont.n))
      newdata[,cont.factors] <- cont.levels[tmp,]
   }

   # THIS BREAKS FOR NESTED MODELS (in model.frame)
   # SEE lsmean.listof() below for a kludge fix
   mat <- t(as.matrix(model.matrix.default(expr,newdata,contrast)[,mindex]))
   # average over levels for auxiliary factors without sum-to-zero contrasts
   if (!is.null(cont.factors)) {
      tmp <- array(mat,c(nrow(mat),ncol(mat)/cont.n,cont.n))
      dimnames(tmp) <- list(dimnames(mat)[[1]],NULL,NULL)
      mat <- apply(tmp,c(1,2),mean)
   }
   mat <- as.matrix(mat)

   # predicted values (THIS IS BROKEN IF THERE ARE EMPTY CELLS!)
   pred <- coef[index] * mat
   if (length(index) == 1)
      pred <- unlist(pred)
   else
      pred <- apply( pred, 2, sum, na.rm = T )
   # standard error
   if(se.fit) {
      cov <- cov.coef(fit)[index,index]
      if(length(index) == 1) 
         se <- sqrt(cov) * abs(mat)
      else
         se <- sqrt(apply(mat*(cov%*%mat),2,sum))
      se <- std.dev(fit) * se
   }
   if(is.null(factors) & is.null(covars))
      newdata <- data.frame(pred=matrix(pred,1,1),row.names=int.factors)
   else {
      newdata <- as.data.frame(newdata[int.factors,c(factors,covars)])
      dimnames(newdata) <- list(int.factors,c(factors,covars))
      names(pred) <- as.character(int.factors)
      newdata$pred <- pred
   }
   if(se.fit) {
      if(is.null(factors) & is.null(covars))
         se <- matrix(se,1,1)
      newdata$se <- se
      names(newdata$se) <- row.names(newdata)
   }
   newdata
}

lsmean.listof <- function(fit,data=eval(attr(fit,"call")$data),
   factors=names(fit[[stratum]]$contrasts),stratum,expr=formula(fit),
   contrast=fit.cont,...)
{
   if(missing(stratum))
	stop("must specify stratum for fit of nested model")
   fit.cont <- NULL
   for (i in names(fit))
      if (i != "(Intercept)")
      fit.cont <- c(fit.cont,fit[[i]]$contrasts)
   lsmean.lm(fit[[stratum]],data=data,factors=factors,expr=expr,
      contrast=contrast,...)
}
margin.plot <- function(object, ...) UseMethod("margin.plot")

margin.plot.default <- function(x.factor,response,group,
   orient = "default", average = "full",
   xlab=xlabs,ylab=ylabs,...)
{
   xlabs <- deparse(substitute(x.factor))
   xlabs[2] <- deparse(substitute(group))
   ylabs <- deparse(substitute(response))

   data <- data.frame(response,x.factor,group)
   names(data) <- c(ylabs,xlabs)

   formula.y <- formula( paste(ylabs,".",sep="~"))
   reduced <- lm( formula.y, data)
   full <- update(reduced,~.^2,singular.ok=T)
 
  if(is.na(pmatch(average,"full")))
      lsm <- lsmean(reduced,data)
   else
      lsm <- lsmean(full,data)

   orient <- pmatch(orient,c("switch","both"),nomatch=0)
   if (orient == 1)
      ylabs <- c("",ylabs)
   if(orient != 1)
      ret <- margin.plot.lm(full,reduced,data,factors=xlabs,lsm=lsm,
	 xlab=xlab[1],ylab=ylab[1],...)
   if(orient > 0)
      ret <- margin.plot.lm(full,reduced,data,factors=rev(xlabs),lsm=lsm,
         xlab=xlab[2],ylab=ylab[2],...)
   invisible(list(ret=ret,reduced=reduced,full=full))
}

margin.plot.lm <- function(full,reduced=full,
   data=eval(full$call$data),factors=all.factors(full,data),
   lsm=lsmean(full,data,factors),margin=margin.labels,error="full",
   bar.plot="lsd",xlab=xlabs,ylab="",xlim=xlims,ylim=ylims,
   pty="s",xaxt="s",...)
{
   # allow for multi-factor horizontal axis
   x.factor <- factors[[1]]
   xlabs <- paste(x.factor,collapse=".")
   factors[[1]] <- xlabs
   x <- interaction( get.list( lsm, xlabs ), drop = T )
   ylims <- range(lsm$pred)
   xlims <- range(unfactor(lsm[[x.factor[1]]]))

   # use margin values if they correspond to levels of first factor
   # lazy evaluation: margin=margin.labels
   margin.labels <- precision(lsmean(reduced,data,x.factor))
   if (margin[1]) {
      names(margin) <- names(margin.labels)
      lsm[[xlabs]] <- ordered(margin[as.character(x)])
      tmpar <- par(pty=pty)
      xlims <- ylim
   }
   if(is.na(pmatch(error,"full"))) {
      error <- reduced
      lsm$se <- lsmean(reduced,data,factors)$se
   }
   else
      error <- full
   ret <- int.plot(error,data,factors,lsm=lsm,xlab=xlab,ylab=ylab,
      xlim=xlim,ylim=ylim,bar.plot=bar.plot,xaxt=xaxt,...)

   # add actual factor 1 levels if using marginal mean values 
   if (margin[1]) {
      if (xaxt != "n") {
         mtext(names(margin),1,1,at=margin)
         right <- par("usr")[1]
         mtext(xlabs,1,1,at=right)
         mtext("margin",1,0,at=right)
      }
      par(tmpar)
   }
   ret$par <- tmpar
   invisible(ret)
}
mplot <- function(a, b, group=NULL, type="p",
   xlim=range(a),ylim=range(b),
   xlab=deparse(substitute(a)),ylab=deparse(substitute(b)),
   main="",sub="",...)
{
   plot(a, b, type = "n",xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,main=main,
      sub=sub, ...)
   mlines(a, b, group, type=type, ...)
   invisible(par(lty=1))
}
mpoints <- function(a, b, group=NULL, type="p", ...)
{
   mlines(a, b, group, type, ...)
}
mlines <- function(a, b, group=NULL, type="l", lty=1:lu, col=1:lu,
   fulltext=F, ...)
{
   if(is.null(group)) {
      lu <- 1
      lines(sort(a), b[order(a)], type=type, lty=lty[1], col=col[1] )
   }
   else
   {
      if (is.factor(group))
         ugroup <- levels(group)
      else
         ugroup <- unique(sort(group))
      char <- as.character(ugroup)
      lu <- length(ugroup)
      lty <- array(lty,lu)
      col <- array(col,lu)
      if (fulltext) 
         text(a, b, label = group)
      else
         for(i in 1:lu) {
            u <- group == ugroup[i]
            if(any(u))
               lines(sort(a[u]), b[u][order(a[u])], lty=lty[i], col=col[i],
                  pch = char[i], type=type)
         }
   }
}
mqqnorm <- function(x,group = abbreviate( qq$labels, 8 ), plot.it = T,
                    xlab = "half-normal quantiles",
                    ylab = "effects", ...)
{
  if( inherits( x, "aov" ) & !exists( "qqnorm.aov" ))
    {
      qq <- list()
      qq$y <- sort( abs( coef(x)[-1] ))
      n <- length( qq$y )
      qq$y <- qq$y * sqrt( n )
      qq$x <- qnorm(( 1 + ppoints( n )) / 2 )
      tmpopt <- options(warn = -1)
      plot( qq$x, qq$y, xlab = xlab, ylab = ylab, ..., type = "n" )
      options(tmpopt)
#     qq <- qqnorm( y, ..., type = "n" )
    }
  else
    qq <- qqnorm(x,...,type="n")
  qq$labels <- abbreviate( names( qq$y ), 8 )
  mpoints(qq$x, qq$y, group, ...)
  ii <- !is.na(qq$x)
  quant <- round(c(1,3)*(1+sum(ii))/4)
  abline(lsfit(qq$x[ii][quant],qq$y[ii][quant]),lty=2)
  invisible(qq)
}
char <- function(...,pch=c(1:9,0,letters,LETTERS),nomatch=".")
{
   x <- c(...)
   pm <- pmatch(pch,x)
   x <- rep(NA,length(x))
   x[pm[!is.na(pm)]] <- pch[!is.na(pm)]
   x[is.na(x)] <- nomatch
   paste(x,collapse="")
}
nested <- function(object, ...) UseMethod("nested")

nested.default <- function( fit, data = eval( attr( fit, "call" )$data ),
   project = proj( fit ), ... )
   nested.aovprojlist( project, data, ... )

# this assumes that the right list of factors goes along with the nesting

nested.aovprojlist <- function( project, data, factors, response="y",
   nesting=c(1,ncol(ref)), ref = projref( project ), ... )
{
   newdata <- data[, factors]
   newdata[[response]] <- apply( ref[ , nesting ], 1, sum )
   projwt( newdata, factors, ... )
}

projwt <- function( data, factors, u = interaction( data[,factors], drop = T ),
   weight = 1, covars )
{
   weight <- array( weight, length( u ) )
   newdata <- as.data.frame( data[ !duplicated( u ), ] )
   if ( !missing( covars )) for (i in covars)
      newdata[[i]] <- tapply( data[[i]], u, mean )
   newdata$weight <- tapply( weight, u, sum )
   row.names( newdata ) <- unique( u )
   newdata
}

projref <- function( project )
{
   ref <- data.frame(Reference=c(project[[1]]))
   for (i in names(project)[-1])
      ref[[i]] <- apply( project[[i]], 1, sum )
   ref
}
nested.tree <- function( data, nesting=names(means),
   means=nested.means(data,nesting[-1]), offset = 0.05, small = 0.02,
   estimate="estimate",bar="stderr", jitter=offset/2, ...)
{
   nn <- length(nesting)

   plot( c(0, 1+nn), range( means[[nn]][[estimate]] ),
      type = "n", xaxt = "n", ... )
   mtext(nesting, 1, 0, at = seq(nn) )

   small <- rep(small,nn)
   offset <- rep(offset,nn)
   jitter <- rep(jitter,nn)

   mean2 <- means[[1]]  
   estim2 <- mean2[[estimate]]
   for ( i in seq(nn-1) ) {
      mean1 <- mean2
      mean2 <- means[[i+1]]
      estim1 <- estim2
      estim2 <- mean2[[estimate]]
      nester <- nesting[i]
      nest2 <- mean2[[nester]]
      nest1 <- mean1[[nester]]
      nr <- nrow(mean1)
      x <- rep(i,nr)
      text( x-offset[i], uncollide(estim1,small[i]), as.character(nest1), adj = 1 )
      for ( j in seq(nr) )
         for (k in estim2[ nest2==nest1[j] ])
            lines(seq(i,length=2),c(estim1[j],k))
      points( x, estim1 )
      if (jitter[i] > 0)
         x <- x + runif( nr, -jitter[i], jitter[i] )
      se.bar( x, estim1, 2 * mean1[[bar]], cap = "" )
   }
   nester <- nesting[nn]
   invisible(text( nn + offset[nn], uncollide( estim2, small[nn] ),
      as.character( mean2[[nester]] ), adj = 0 ))
}

nested.means <- function(data, nesting, parameter="parameter",
   carry=c("estimate","stderr","ddf"))
{
   x <- list( data[ 1, carry ] )
   root <- as.character(data[ 1, parameter ])
   x[[1]][[root]] <- names(x) <- root
   for (i in seq(nesting)) {
      x[[nesting[i]]] <- data[ data[[parameter]] == nesting[i],
         c(nesting[1:i],carry) ]
   }
      x[[2]][[root]] <- rep(root,nrow(x[[2]]))
   x
}
plotter <- function(..., term = unix("echo $TERM"), 
	host = unix("echo $SHOME"), height = 24, las = 1)
{
   switch(term,
      {
         printer(height = height, ...)
      },
      xterm = {
	 if(is.na(match(host,"/usr/stat/splus")))
            x11(...)
         else
            motif(...)
         par(las = las, ...)
      },
      tek4014 =,
      versaterm = {
         tek4014(...)
      })
   invisible(.Device)
}
power.curve <- function(n = 100, sigma = 1, alpha = 0.05, dstd = seq(-5, 5, length = 100),df)
{
if(missing(df)){
   cstd <- qnorm(1 - alpha/2)
   data.frame(dif = (dstd * sigma)/sqrt(n/2),
      pow = 1 - pnorm(cstd - dstd) + pnorm( - cstd - dstd))
}
else {
   cstd <- qt(1 - alpha/2,df)
   data.frame(dif = (dstd * sigma)/sqrt(n/2), pow = 1 - pt(cstd - dstd,df) +
      pt( - cstd - dstd,df))
	}
}
precision <- function(x, se, digits = 0)
{
        if(is.list(x)) {
                se <- x$se
                p <- x$pred
                if(is.data.frame(x))
                        names(p) <- names(se) <- dimnames(x)[[1]]
                x <- p
        }
        if(length(se) == 1 && se == 0) {
                x[abs(x) < 1e-05] <- 0
        }
        else {
                if(length(se) > 1 && length(se) != length(x))
                        stop("x and se must have the same length")
                scale <- 10^(digits + 1 - round(log10(se)))
                p <- round(scale * x)/scale
                p[is.na(p)] <- x[is.na(p)]
                p[abs(p) < 1e-05] <- 0
                names(p) <- names(x)
                p
        }
}

#precision <- function(value,se)
#{
#   shifts <- 10^(1-round(log10(2*se)))
#   round(value*shifts) / shifts
#}
   
se.bar <- function(xpos,ypos,bar,mod="",cap="SE",adj=0,width=.01,lty=1,horiz=F,...)
{
   if(!all(is.na(bar)) && max(bar) > 0 && min(bar) >= 0) {
      if(length(xpos) != length(ypos))
         stop("xpos and ypos must have the same length")
      if(length(xpos) != 1 && length(bar) != 1 && length(bar) > length(xpos))
         stop("xpos and bar must have the same length or one must be scalar")

      bar <- array(bar,length(xpos))
      xpos <- xpos[!is.na(bar)]
      ypos <- ypos[!is.na(bar)]
      bar <- bar[!is.na(bar)]
      lty <- rep( lty, length( bar ) )
      if (horiz) {
         top <- xpos+bar*.5
         bot <- xpos-bar*.5
         for (i in 1:length(ypos))
            lines(c(bot[i],top[i]),rep(ypos[i],2))
   
         left <- right <- ypos
         if (width > 0) {
            width <- width * diff(par("usr"))[3]
            left <- ypos-width
            right <- ypos+width
            for (i in 1:length(ypos)) {
               lines(rep(bot[i],2),c(left[i],right[i]))
               lines(rep(top[i],2),c(left[i],right[i]))
            }
         }
      }
      else {
         top <- ypos+bar*.5
         bot <- ypos-bar*.5
         for (i in 1:length(xpos))
            lines(rep(xpos[i],2),c(bot[i],top[i]),lty=lty[i])
   
         left <- right <- xpos
         if (width > 0) {
            width <- width * diff(par("usr"))[1]
            left <- xpos-width
            right <- xpos+width
            for (i in 1:length(xpos)) {
               lines(c(left[i],right[i]),rep(bot[i],2))
               lines(c(left[i],right[i]),rep(top[i],2))
            }
         }
         if (cap != "") {
            if(adj)
               text(left,ypos,paste(cap,mod,"=",sep=""),adj=1)
            else
               text(right,ypos,paste("=",cap,mod,sep=""),adj=0)
         }
      }
   }
   list(xpos=xpos,ypos=ypos,bar=bar,mod=mod,cap=cap,adj=adj,width=width)
}
burst <- function(root,branch,pos=c(2,1))
{
   for (i in branch)
      lines(c(root,i),pos)
}
branch <- function(root,limb,pos=1:2,flip=F)
{
   if (flip)
      for (i in limb)
         lines(pos,c(root,i))
   else
      for (i in limb)
         lines(c(root,i),pos)
}
split.plot <- function(response,nest.factors,covars,weights,data)
{
   nest.int <- interaction( get.list( data, nest.factors ), drop = T )
   # need to check that nest.factors are factors
   # how to add weighting to fit?
   formula.y <- formula( paste( response, "~",
                               paste( nest.factors, sep = "*" )))
   nest.fit <- lm( formula.y, data)
   # whole plot
   Within <- data
   if (missing(weights))
      weights <- "weights"
   if(is.na(pmatch(weights,names(data))))
      Within[[weights]] <- rep(1,nrow(data))
   Within[[response]] <- resid(nest.fit)
   # split plot
   Between <- data.frame(data[!duplicated(nest.int),],
      row.names=unique(nest.int))
   Between[[response]] <- predict(nest.fit)[!duplicated(nest.int)]

   Between[[weights]] <- tapply(Within[[weights]],nest.int,sum)
   browser()
   # covariates
   if(!missing(covars)) {
      covar.data <- data
      for (i in covars) {
         covar.data[[response]] <- data[[i]]
         Within[[i]] <- resid(nest.fit,covar.data)
         Between[[i]] <- predict(nest.fit,covar.data)
      }
   }
   list(Between=Between,Within=Within)
}
subsamp <- function(variance=c(1,4),diff=1,alpha=.05,
   sub=10^seq(0,log10(50),length=100),eu=tot/sub,tot=100)
{
ems <- variance[1]+variance[2]*sub
delta <- eu/(2*ems)
crit <- qt(1-alpha/2,eu-1)
pow <- pt(abs(diff)*sqrt(delta)+crit,eu-1)-pt(abs(diff)*sqrt(delta)-crit,eu-1)
crit <- qt(1-alpha/2,eu-1)/sqrt(delta)
delta <- delta*diff^2
data.frame(sub=sub,eu=eu,df=eu-1,ems=ems,delta=delta,pow=1-pow)
}

power.curve <- function(n=100,sigma=1,alpha=.05,dstd=seq(-5,5,length=100))
{
cstd <- qnorm(1-alpha/2)
data.frame(dif=dstd*sigma/sqrt(n/2),pow=1-pnorm(cstd-dstd)+pnorm(-cstd-dstd))
}
trans.plot <- function(fit,
	predictors=all.predictors(fit)[1:2],
	terms=predict(fit,type="terms"),
	xlab="product of margins",ylab="interaction",...)
{
	a <- terms[,predictors[1]]
	b <- terms[,predictors[2]]
	ab <- terms[,paste(predictors[1:2],collapse=":")]
	if (is.null(ab))
		stop(paste(paste(predictors[1:2],collapse=":"),
			"is not in list of terms"))
	m <- attributes(terms)$constant
	add <- a*b/m
	plot(add,ab,xlab=xlab,ylab=ylab,...)
	data <- data.matrix(list(x=add,y=ab))
	fit <- lm(y~x,data)
	lines(sort(add),predict(fit)[order(add)])
	c(power=1-coef(fit)["x"],se=se.coef(fit)["x"])
}
tukey.plot <- function(x.factor,response,group,...,orient="default",
   xlab=xlabs,ylab=ylabs)
{
   xlabs <- deparse(substitute(x.factor))
   xlabs[2] <- deparse(substitute(group))
   ylabs <- c(deparse(substitute(response)),"")

   data <- data.frame(y=response,a=x.factor,b=group)
   names(data) <- c(ylabs[1],xlabs)

   reduced <- lm(y~.,data)
   data$inter <- predict(reduced)^2
   formula.y <- formula( paste( paste( ylabs[1], xlabs[1], sep = "~" ),
                               xlabs[2], "inter", sep = "+" ))
   full <- lm( formula.y, data,singular.ok=T)
   lsm <- lsmean(full,data,xlabs,adjust.covar=F)

   orient <- pmatch(orient,c("switch","both"),nomatch=0)
   if(orient != 1)
      ret <- margin.plot.lm(full,reduced,data,factors=xlabs,lsm=lsm,
	 xlab=xlab[1],ylab=ylab[1],...)
   if(orient > 0)
      ret <- margin.plot.lm(full,reduced,data,factors=rev(xlabs),lsm=lsm,
         xlab=xlab[2],ylab=ylab[2],...)
   invisible(list(reduced=reduced,full=full,data=data,lsm=lsm,ret=ret))
}

mandel.plot <- function(x.factor,response,group,...,orient="default",
   xlab=xlabs,ylab=ylabs)
{
   xlabs <- deparse(substitute(x.factor))
   xlabs[2] <- deparse(substitute(group))
   ylabs <- c(deparse(substitute(response)),"")

   data <- data.frame(y=response,a=x.factor,b=group)
   names(data) <- c(ylabs[1],xlabs)

   reduced <- lm(y~.,data)
   data$inter <- unfactor(lsmean(reduced,data,se.fit=F)[[1]][codes(x.factor)])
   formula.y <- formula( paste( paste( ylabs[1], xlabs[1], sep = "~" ),
                               xlabs[2], paste( xlabs[2], "inter", sep = ":"),
                               sep = "+" ))
   full <- lm( formula.y, data, singular.ok = T )
   lsm <- lsmean(full,data,xlabs,adjust.covar=F)

   orient <- pmatch(orient,c("switch","both"),nomatch=0)
   if(orient != 1)
      ret <- margin.plot.lm(full,reduced,data,factors=xlabs,lsm=lsm,
	 xlab=xlab[1],ylab=ylab[1],...)
   if(orient > 0)
      ret <- margin.plot.lm(full,reduced,data,factors=rev(xlabs),lsm=lsm,
         xlab=xlab[2],ylab=ylab[2],...)
   invisible(list(reduced=reduced,full=full,data=data,lsm=lsm,ret=ret))
}
uncollide <- function(x, small = 0.02)
{
        if(length(x) > 1) {
                o <- order(x)
		n <- names(x)
                x <- sort(x)
                d <- (small * diff(range(x)) - diff(x))/2
                d <- d * (d > 0)
                x[o] <- (x - c(d, 0) + c(0, d))
		names(x) <- n
        }
	x
}