info      MACRO
) MacAnova macros that are useful in experimental design and analysis
)
) aberration2    Computes aberration in a fractioned 2 series
) aliases2       Gets aliases in fractioned 2 series
) aliases3       Gets aliases in fractioned 3 series
) allaliases2    Complete aliases structure in fractioned 2 series
)                factorial
) all3anova      Looks at all hierarchical models in a three factor
)                anova and sorts them by Cp
) all4anova      Looks at all hierarchical models in a four factor anova
)                and sorts them by Cp
) boxcoxvec      Gets the error SS for several boxcox transformations
) buildfactor    Creates factor for balanced multi-factor experiments
)                with data in canonical order
) choosedef2     Chooses defining contrasts for a blocked 2 series
)                factorial
) choosegen2     Chooses generators for a 2 series fractional factorial
) confound2      Confounds a 2 series factorial into blocks
) confound3      Confounds a 3 series factorial into blocks
) doconfound2    Improved front end to confound2()
) docontrast     Improved front end to contrast()
) doff2          Improved front end to choosegen2,alias2, etc
) dorandtest     Front end for randt2 and randsign
) ems            Computes the expected mean squares for the terms in the
)                ANOVA in a specified mode with specified random terms
) ffdesign2      Determines factor/level combinations to use for
)                fractioned 2 series
) findncp        Computes the noncentrality parameter in a one-way ANOVA
) findpower      Computes the power of the F-test in a one-way ANOVA
) findsampsize   Computes the sample size needed to attain a requested
)                minimum power in a one-way ANOVA
) interactplot   Makes interaction plots of marginal means in a
)                factorial
) interblock     Does interblock recovery for incomplete block designs
) mixed          Performs a mixed effects ANOVA, computing appropriate
)                error MS and approximate df as needed.
) pairwise       Does pairwise comparisons and generates an "underline"
)                diagram
) quadmax        Finds the location of the maximum of a quadratic
)                function, with optional linear equality and inequality
)                constraints on the solution
) randsign       Does a permutation paired t-test
) randt          Does a permutation t-test using values from a single
)                vector
) randt2         Does a permutation t-test combininng values from two
)                vectors
) reml           Does restricted maximum liklihood estimation (still
)                beging tested)
) rscanon        Does canonical analysis of 2nd order response surface
)                design
) sidebyside     Makes a side-by-side plot of model effects
) stdordlabels   Make letter labels for effects in standard order
) varcomp        Computes ANOVA estimates of variance components.
) yatesplot      Do plots of effects for 2 series designs
)
) Internal macros not called directly by users
)
) buildRFStr     Parses a mixed model and constructs matrices describing
)                random and fixed parts of the model.
) colproduct     Automatically loaded by ems and used to compute all
)                element-wise products of the columns of two matrices
) makemat        Automatically loaded by ems and used to compute various
)                basis matrices.
) quadmaxlin     Automatically read by quadmax and used to find the
)                maximum of x'Ax + b'x subject to Cx = y
)
) Written by Gary Oehlert, Applied Statistics, University of Minnesota,
) St. Paul, MN 55108 (gary@stat.umn.edu), with some
) modifications by Christopher Bingham (kb@stat.umn.edu)
) Version of 940907
) Version of 950313  Fixed bug in ffdesign2() and allaliases2()
) Version of 960319  Added all3anova(), all4anova();
)                    Added functionality to pairedcomp()
)                    Added blocking to rscanon()
) Version of 970224  Added error: keyword to pairedcomp()
)                    Added quadmax()
) Version of 970311  Fixed bugs in quadmax(), quadmaxlin()
) Version of 970520  Added boundedness checking to quadmax()
) Version of 970522  Add ems()
) Version of 971029  Corrections to ems() and addition of usage comments
) Version of 980129  Correction to makemat(), addition of varcomp() and mixed()
) Version of 980302  Added useneg:T to mixed()
)                    Added warnings to mixed() and varcomp()
) Version of 980306  Fixed bounds bug in quadmax()
) Version of 980318  Changed makemat() to make it faster
) Version of 980619  Changed to new format, not requiring nlines
) Version of 980803  Bug fix in mixed()
) Version of 980903  Fixes to aliases2() and aliases3()
) Version of 981210  Fix to all3anova()
) Version of 991118  Changed choosegen2() to do minimum aberration
)                    Added choosedef2()
)                    Added interactplot()
)                    Added interblock()
)                    Changed pairedcomp to pairwise()
)                    Removed the obsolete typeIIIss() macro
) Version of 000129  Stripped most '$$' and put DOLLARS on headers and
)                    made some conversion to using argvalue() and keyvalue()
) Version of 000208  Corrected errors, improved headers, modified
)                    interactplot()
) Version of 000430  Corrected more errors, made most macros out-of-line
) Version of 000620  Added sidebyside()
) Version of 000907  Modified ems() to accept variates
)                    Added reml() to do restricted maximum likelihood
)                    Modifed supporting macros such as makemat()
) Version of 000911  Cosmetic cleaning of code
) Version of 020322  A number of keywords are specified with a trailing
)                    '*' for added flexibility
)                    boxcoxvec() uses new functions popmodel(), pushmodel()
)                    Many expressions involving A' %*% B have been
)                    replaced by A %c% B and %c% has been used in other
)                    situations as well.
)                    Several cases of trace (A %*% B) have been replaced
)                    by sum(vector(A * B'))
) Version of 020829  Modified some prefatory comments to macros
)                    Changed some labeling in sidebyside()
)                    boxcoxvec() now can be used as
)                    boxcoxvec(model [,power:pow])
) Version of 021022  New macro buildfactor() to create factor for
)                    balanced multifactor design in canonical order
) Version of 021227  Fixed bug in allanova3() and made minor changes to
)                    error messages containing '$S'
)                    buildfactor() recognizes 'reverse:T'
) Version of 040827  Changed reml() to denote random terms by either
)                    random:termlist or Z:termlist; reml is also faster
)                    and returns predictions of random effects
)                    Updated makemat
) Version of 040903  Added findpower, findncp, findsampsize
) Version of 050730  Added docontrast
) Version of 050802  Added dorandtest
) Version of 050806  Added rcb:T and ncp1/ngrp to findsampsize
) Version of 050810  Added error bars and ls means to interactplot
) Version of 050814  Changed reml to deal with varcomps all zero
) Version of 050820  Added yatesplot and stdordlabels
) Version of 050820  Fixed bug in rscanon
) Version of 050821  Added aberration2
) Version of 050824  New, faster version of choosegen2, plus recuraber
%info%

====> buildfactor <====
buildfactor     MACRO DOLLARS
) Macro to create factors for a regular factorial treatment
) Usage:
)  a <- buildfactor(j,dims [,reverse:T])
)   j           positive integer identifying subscript for which to
)               compute a factor
)   dims        vector of positive integers with length(dims) >= j
)   reverse:T   first subscript changes slowest, last changes fastest
)   a           factor of length N = prod(dims) with dims[j]
)               levels
)
) Each of the N cases is assumed to be identified by k = length(dims)
) subscripts, with i-th subscript running from 1 to dims[i]
)
) Without reverse:T: subscript 1 changes fastest and subscript k slowest
) With reverse:T: subscript k changes fastest and subscript 1 slowest
)
)) Written by C. Bingham (kb@umn.edu) 021018
)) 021227 added keyword reverse
# a <- $S(j,dims), integer vector dims, integer j <= length(dims)
@j <- argvalue($1,"subscript number","positive count")
@dims <- argvalue($2,"dimensions","positive integer vector")
@reverse <- keyvalue($K,"revers*","TF",default:F)
@ndims <- length(@dims)
if (@j > @ndims){
	error("subscript number > number of subscripts")
}
if (delete(@reverse,return:T)) {
	@dims <- reverse(@dims)
	@j <- @ndims - @j + 1
}
@N1 <- if (@j == 1){
	1
} else {
	prod(@dims[run(@j-1)])
}
@N2 <- if (@j == @ndims){
	1
} else {
	prod(@dims[-run(@j)])
}
@result <- factor(rep(rep(run(@dims[@j]),rep(@N1,@dims[@j])), @N2))
delete(@j,@dims,@ndims,@N1,@N2)
delete(@result,return:T)
%buildfactor%

====> choosegen2 <====
choosegen2  MACRO   DOLLARS OUTLINE
) choosegen2(k,p, all:T or tries:m [,res:r]) finds a set of generators
) for a 2^(k-p) fractional factorial design.  Either all:T or tries:m
) must be an argument
)
) When all:T is an argument, all potential generators are searched.
)
) When tries:m is an argument, m random sets of generators are searched.
)
) When res:r is not an argument, a set of generators giving the greatest
) resolution among those searched is returned.  When res:r is an argument,
) the first set of generators found giving resolution r is returned; if
) no such set is found, a warning is printed and the generators giving
) the best resolution found are printed.
)
) The value returned is
)   structure(resolution:bestres,generators:bestgen,
)         aberration:abvec, basis:genmat)
)    bestres is the resolution of the best design found
)    bestgen is a CHARACTER vector with the defining contrasts of the
)      best design found
)    abvec is an integer vector of length k with abvec[i] containing the
)      numbers of interactions with i letters in the generating set
)    basis is a p by k matrix of 0's and 1's, with the 1's in row i
)      defining the factors in the bestgen[i]
))Written by Gary W. Oehlert (gary@stat.umn.edu)
))020307 Modified error messages; replace "res" by "res*" in keyvalue()
))Version 020307
))Version 050824
# usage: $S(k,p,all:T,tries:m,res:r)
if($v != 2){
	error("usage is $S(k,p[,all:T,tries:m,res:r])", macroname:F)
}
@k <- argvalue($1,"$1","positive count")
@p <- argvalue($2,"$2","positive count")
if(@k < 3 || @k > 32) {
    error("k must be an integer between 3 and 32")
}
if(@p >= @k ) {
	error(paste("p =",@p,"not a positive integer < k =",@k))
}
@res <- 0
@tries <- 0
if($k == 0) { # no keywords
    error("you must specify one of all:T or tries:m")
}
@res <- keyvalue($K,"res*","positive count")
if(alltrue(!isnull(@res),@res < 3)) {
    error("res must be integer >= 3")
}
@tries <- keyvalue($K,"tries","positive count")
@all <- keyvalue($K,"all","TF", default:F)

if(@all && !isnull(@tries)) {
    error("you can't use both all:T and tries:n")
}
if(isnull(@tries) && !@all) {
    error("you must specify one of all:T or tries:m")
}
@rhs <- 2^run(@k-@p,@k-1)
@bestab <- rep(1,@k+1)
@bestr <- 0
@besttry <- 0

@lhs <- run(0,2^(@k-@p)-1)
@gt <- @lhs > @lhs'
@bindigl <- (@lhs %& 2^run(0,@k-@p-1)') > 0
@bindig <- 1*@bindigl
for(@i,run(length(@lhs))) {
	@m <- max(run(@k-@p)[@bindigl[@i,]'])
	if(isnull(@m)) {
		@m <- 0
	}
	if(@m < @k-@p-1) {
		for(@j,run(@m+2,@k-@p)) {
			@out <- !@bindigl[,@j-1] && @bindigl[,@j]
			@gt[@out,@i] <- F
		}
	}
}
@thislist <- @lhs[@gt[,1]]
@gt <- @gt[-1,-1]
@lhs <- @lhs[-1,]
if(@all) {
	for(@i,run(length(@thislist))) {
		# now go through each value in the new list
		@start <- @thislist[@i]
		recuraber(@start,@p,@gt,@besttry,@bestab,@lhs,@rhs)
	}
} else {
    for(@i,run(@tries)) {
        @thistry <- @lhs[grade(runi(length(@lhs)))<=@p]+@rhs
        @thisalias <- vector(0,@thistry[1])
        if(length(@thistry) > 1) {
            for(@j,run(2,length(@thistry))) {
                @thisalias <- vector(@thisalias %^ vector(0,@thistry[@j])')
            }
        }
        @thisr <- min(nbits(@thisalias[-1]))
        @thisab <- tabs(@thisalias[-1],nbits(@thisalias[-1])+1,count:T)
        for(@j,run(length(@thisab))) {
            if(@thisab[@j] < @bestab[@j]) {
                @besttry <- @thistry
                @bestab <- @thisab
                @bestr <- @thisr
                break
            } else {
                if(@thisab[@j] > @bestab[@j]) {
                    break
                }
            }
        }
        if(alltrue(!isnull(@res), @thisr >= @res)) {
            break
        }
    }
}

@letters <- vecread(string:"ABCDEFGHJKLMNOPQRSTUVWXYZ1234567",\
	bychar:T)[run(@k)]

@out <- rep("A",@p)
@genmat <- NULL
for(@i,run(@p)) {
    @tmp <- @besttry[@i] %& 2^run(0,@k-1)
    @genmat <- vconcat(@genmat,1*(@tmp>0)')
    @out[@i] <- paste(@letters[@tmp>0],sep:"")
}
delete(@i,@j,@p,@tmp,@letters,@res,\
	@all,@besttry,@lhs,@rhs)
@bestr <- min(run(0,length(@bestab)-1)[@bestab > 0])

structure(resolution:delete(@bestr,return:T),generators:delete(@out,return:T),\
	aberration:padto(delete(@bestab,return:T)[-1],delete(@k,return:T)),\
	basis:delete(@genmat,return:T))
%choosegen2%

recuraber MACRO DOLLARS
# recuraber is a utility routine used in finding fractions
# of large design.  It recursively explores a minimal set of
# potential generators
# $1 is vector containing current vals
# $2 is desired final length
# $3 is gt matrix
# $4 is besttry
# $5 is bestaber
# $6 is lhs
# $7 is rhs
# if (length(currentvals) < desired val)
@curlen <- length($1)
if(@curlen < $2) {
	# if last curr val is already at the end, we can't go further
	# so we might as well give up on this one
	if($1[@curlen] == $6[length($6)]) {
		$1 <- $1[-@curlen]
		return(NULL)
	}
	# this list is combinations greater than current effect
	@thislist <- $6[$3[,$1[@curlen]]]
	for(@i,run(length(@thislist))) {
		# now go through each value in the new list
		$1 <- vector($1,@thislist[@i])
		recuraber($1,$2,$3,$4,$5,$6,$7)
	}
	$1 <- $1[-@curlen]
} else {
	# here we have a set of $2 (k) effects, so compute aberration
        #@thistry <- @lhs[these elements]+@rhs
        @thistry <- $1+$7
        @thisalias <- vector(0,@thistry[1])
        if(length(@thistry) > 1) {
            for(@j,run(2,length(@thistry))) {
                @thisalias <- vector(@thisalias %^ vector(0,@thistry[@j])')
            }
        }
        @thisab <- tabs(@thisalias[-1],nbits(@thisalias[-1])+1,count:T)
        for(@j,run(length(@thisab))) {
            if(@thisab[@j] < $5[@j]) {
                #@besttry <- @thistry
                $4 <- @thistry
                #@bestab <- @thisab
                $5 <- @thisab
		$1 <- $1[-@curlen]
                return(NULL)
            } else {
                if(@thisab[@j] > $5[@j]) {
			$1 <- $1[-@curlen]
                    return(NULL)
                }
            }
        }
	$1 <- $1[-@curlen]
	return(NULL)
}
%recuraber%

====> choosedef2 <====
choosedef2  MACRO  DOLLARS OUTLINE
) choosedef2(k,p,all:T or tries:m) finds a set of defining contrasts for
) a 2^k factorial design in 2^p blocks.  Either all:T or tries:m must be
) an argument
)
) When all:T is an argument, all potential defining contrasts are
) searched.
)
) When tries:m is an argument, m random sets of are searched.
)
) choosedef2 determines the set of defining contrasts for the minimum
) aberration design among those searched.
)
) The value returned is
)  structure(generators:bestgen, aberration:abvec,basis:genmat)
)    bestgen is a CHARACTER vector of length p containing the names
)      of the defining contrasts found
)    abvec is an integer vector of length k with abvec[i] containing the
)      numbers of interactions with i letters in the generating set
)    genmat is a p by k matrix of 0's and 1's, with the 1's in row i
)      defining the factors in the bestgen[i]
)) Written by Gary W. Oehlert (gary@stat.umn.edu)
)) Version of 000209
# usage: $S(k,p,all:T,tries:m)
if($v != 2){
	error("usage: $S(k,p[,all:T,tries:m])", macroname:F)
}
@k <- argvalue($1,"$1","positive count")
@p <- argvalue($2,"$2","positive count")
if(@k < 3 || @k > 32) {
    error("k must be an integer between 3 and 32")
}
if(@p >= @k ) {
	error(paste("p =",@p,"not a positive integer < k =",@k))
}
@tries <- 0
if($k == 0) { # no keywords
    error("you must specify one of all:T or tries:m")
}
@tries <- keyvalue($K,"tries","positive count")
@all <- keyvalue($K,"all","TF", default:F)
if(@all && !isnull(@tries)) {
    error("you cannot use both all:T and tries:n")
}
if(isnull(@tries) && !@all) {
    error("you must specify one of all:T or tries:m")
}
@pbase <- @p
while( (2^(@k-@pbase) - 1) < @p) {
    @pbase <- @pbase - 1
}
@lhs <- run(2^(@k-@pbase)-1,1)
@rhs <- 2^run(@k-@pbase,@k-1)
if(@pbase < @p) {
	@rhs <- vector(@rhs,rep(0,@p-@pbase))
}
delete(@pbase)
@bestab <- rep(1,@k+1)
@besttry <- 0
@k2 <- length(@lhs)

if(@all) {
    @basetry <- rep(1,@p)
    for(@i,run(0,(@k2^@p)-1)) {
        @tmp <- @i
        for(@j,run(@p)) {
            @basetry[@j] <- (@tmp %% @k2) + 1
            @tmp <- floor(@tmp / @k2)
        }
        @thistry <- @lhs[@basetry]+@rhs
        @thisalias <- vector(0,@thistry[1])
        if(length(@thistry) > 1) {
            for(@j,run(2,length(@thistry))) {
                @thisalias <- vector(@thisalias %^ vector(0,@thistry[@j])')
            }
        }
        @thisab <- tabs(@thisalias[-1],nbits(@thisalias[-1])+1,count:T)
		#print(@thistry,@thisab,@bestab,@thisalias,@besttry)
        for(@j,run(length(@thisab))) {
            if(@thisab[@j] < @bestab[@j]) {
                @besttry <- @thistry
                @bestab <- @thisab
                break
            } else {
                if(@thisab[@j] > @bestab[@j]) {
                    break
                }
            }
        }
    }
} else {
    for(@i,run(@tries)) {
        @thistry <- @lhs[floor(runi(@p)*@k2)+1]+@rhs
        @thisalias <- vector(0,@thistry[1])
        if(length(@thistry) > 1) {
            for(@j,run(2,length(@thistry))) {
                @thisalias <- vector(@thisalias %^ vector(0,@thistry[@j])')
            }
        }
        @thisab <- tabs(@thisalias[-1],nbits(@thisalias[-1])+1,count:T)
        for(@j,run(length(@thisab))) {
            if(@thisab[@j] < @bestab[@j]) {
                @besttry <- @thistry
                @bestab <- @thisab
                break
            } else {
                if(@thisab[@j] > @bestab[@j]) {
                    break
                }
            }
        }
    }
}

@letters <- vecread(string:"ABCDEFGHJKLMNOPQRSTUVWXYZ1234567",\
	bychar:T)[run(@k)]

@out <- rep("A",@p)
@genmat <- NULL
for(@i,run(@p)) {
    @tmp <- @besttry[@i] %& 2^run(0,@k-1)
    @genmat <- vconcat(@genmat,1*(@tmp>0)')
    @out[@i] <- paste(@letters[@tmp>0],sep:"")
}
delete(@i,@j,@tmp,@letters,@thisab,@besttry,@thistry,@thisalias,@tries,\
	@p, @k2, @lhs, @rhs)
structure(generators:delete(@out,return:T),\
	aberration:padto(delete(@bestab,return:T)[-1],delete(@k,return:T)),\
	basis:delete(@genmat,return:T))
%choosedef2%

====> aliases2 <====
aliases2        MACRO  DOLLARS OUTLINE
) aliases2(basis, effect:vec [,length:m]) finds aliases of vec in a
) 2^(k-p) fractional factorial design.  It returns a CHARACTER vector
) with typical elements "I", "AB" and "-ABF".
)
) basis, a p x k matrix of 0's, 1's and -1's, contains linearly
) independent generators (defining contrasts) for the aliasing, one row
) for each generator and one column for each factor in the design.  A 1
) or -1 indicates a factor is present in the alias, with -1 indicating a
) changed sign.  For example, 1 0 1 0 0 -1 means -ACF is a generator
) (alias of I).
)
) vec is a vector of 0's and 1's of length k with 1's indicating the
) presence of a factor in an effect.
)
) If length:m is an argument, only aliases of length m are returned.
)
) aliases2(basis [,length:m]) is equivalent to aliases2(basis,
) effect:rep(0,k) [,length:m]), returning aliases of I.
)
)) Written by Gary W. Oehlert, gary@stat.umn.edu
)) Version of 000209
# usage: $S(basis [,effect:vec] [,length:m])
if ($v != 1 || $k > 2){
	error("usage: $S(basis [,effect:vec] [,length:m])", macroname:F)
}
@basis <- matrix(argvalue($1,"basis","integer matrix"))
if (max(abs(vector(@basis))) > 1){
    error("basis not matrix with elements -1, 0 and 1")
}
if(isvector(@basis)){
	@basis <- @basis'
}#change column to row vector
@k <- ncols(@basis)
if(@k > 25) {
    error("$S() doesn't work with more than 25 factors", macroname:F)
}
@p <- nrows(@basis)
@vec <- rep(0,@k)

if($k > 0) { # process keywords
	@vec <- vector(keyvalue($K,"effect*","integer vector",default:@vec))
	if(nrows(@vec) != @k || max(abs(@vec)) > 1){
		error(paste("Value of 'effect' not vector of -1, 0, and 1 of length",\
			@k))
    }
	@length <- keyvalue($K,"length","nonnegative count",default:0)
}else{
	@vec <- rep(0,@k)
	@length <- 0
}
@letters <- vecread(string:"ABCDEFGHJKLMNOPQRSTUVWXYZ",bychar:T)[run(@k)]
@out <- "I"

if(sum(abs(@vec)) > 0) {
    @out <- paste(@letters[@vec!=0],sep:"")
    if((sum(@vec < 0) %% 2) != 0) {
        @out <- paste("-",@out,sep:"")
    }
}
@counter <- rep(0,@p)
for(@i,run(2^@p-1)) {
    @counter[1] <- @counter[1]+1
    for(@l,run(@p)) {
        if(@counter[@l] <= 1) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }


    @which <- (vector(sum(@counter*abs(@basis))) + @vec) %% 2
    if(@length <= 1 || sum(@which)==@length) {
        if(sum(@which)==0) {
            @tmp <- "I"
        } else {
            @tmp <- paste(@letters[@which!=0],sep:"")
        }
        if(((sum(vector(@counter*@basis) < 0) +\
			sum(@vec < 0)) %% 2) != 0) {
            @tmp <- paste("-", @tmp, sep:"")
        }
        @out <- vector(@out,@tmp)
    }
} #for(@i,run(2^@p-1))
delete(@i, @l, @p, @k, @vec, @length, @tmp, @letters, @counter,@basis, @which)
delete(@out, return:T)
%aliases2%

====> aliases3 <====
aliases3          MACRO  DOLLARS OUTLINE
) aliases3(basis, effect:vec [,length:m]) finds aliases of vec in a
) 3^(k-p) fractional factorial design.  It returns a CHARACTER vector
) with typical elements "I", "A^1 C^2" and "A^1 B^2 C^1".
)
) basis, a p x k matrix of 0's, 1's and 2's, contains linearly
) independent generators (defining contrasts) for the aliasing, one row
) for each generator and one column for each factor in the design.  A 1
) or 2 indicates a factor is present in the alias. For example, 1 0 2 0
) 0 1 means AC^2F is a generator (alias of I).
)
) vec is a vector of 0's and 1's of length k with 1's indicating the
) presence of a factor in an effect.
)
) If length:m is an argument, only aliases of length m are returned.
)
) aliases3(basis [,length:m]) is equivalent to aliases3(basis,
) effect:rep(0,k) [,length:m]), returning aliases of I.
)
)) Written by Gary W. Oehlert, gary@stat.umn.edu
)) Version of 000430
# usage: $S(basis,effect:vec,length:j)
@basis <- argvalue($1,"basis","nonnegative integer matrix")
if (max(vector(@basis)) > 2){
    error("basis not matrix with elements 0, 1 and 2")
}
if(isvector(@basis)){
	@basis <- @basis'
}#change column to row vector
@k <- ncols(@basis)
@p <- nrows(@basis)
@vec <- rep(0,@k)
@length <- 0
if(@k > 25) {
    error("$S() doesn't work with more than 25 factors",macroname:F)
}

if($k > 0) { # process keywords
	@vec <- keyvalue($K, "effect*", "nonnegative integer vector",default:@vec)
	if (nrows(@vec) != @k || max(@vec) > 2){
		error(paste("value for 'effect' not length",@k,"vector of 0, 1 and 2"))
	}
	@length <- keyvalue($K, "length", "positive count")
	if (isnull(@length)){
		@length <- 0
	}
}
@letters <- vecread(string:"ABCDEFGHJKLMNOPQRSTUVWXYZ",bychar:T)[run(@k)]
@out <- "I"

if(sum(@vec) > 0) {
    @out <- ""
    @first <- min(run(@k)[@vec > 0])
    if(@vec[@first] == 2) {
        @vec <- 2*@vec %% 3
    }
    for(@l,run(@first,@k)) {
        if(@vec[@l] > 0) {
            @out <- paste(@out,@letters[@l],"^",\
				@vec[@l]," ",sep:"")
        }
    }
	delete(@first)
}
@counter <- rep(0,@p)
for(@i,run(3^@p-1)) {
    @D <- @counter[1] <- @counter[1]+1
    for(@l,run(@p)) {
        if(@counter[@l] <= 2) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @which <- (vector(sum(@counter*@basis)) + @vec) %% 3

    if(@length <= 1 || sum(@which > 0) <= @length) {
        if(sum(@which)>0) {
            @this <- ""
            @first <- min(run(@k)[@which > 0])
            if(@which[@first] == 2) {
                @which <- 2*@which %% 3
            }
            for(@l,run(@first,@k)) {
                if(@which[@l] > 0) {
                    @this <- paste(@this,@letters[@l],"^",\
					@which[@l]," ",sep:"")
                }
            }
        } else {
            @this <- "I"
        }
		delete(@first)
        @out <- vector(@out,delete(@this, return:T))
    }
}
delete(@D,@length,@which,@i,@l,@p,@k,@counter,@basis,@vec)
delete(@out, return:T)
%aliases3%

====> allaliases2 <====
allaliases2    MACRO  DOLLARS OUTLINE
) allaliases2(basis) finds the full set of aliases in a 2^(k-p)
) fractional factorial design.  It returns a CHARACTER vector with
) typical elements "I", "AB" and "-ABF".
)
) basis, a p x k matrix of 0's, 1's and -1's, contains linearly
) independent generators (defining contrasts) for the aliasing, one row
) for each generator and one column for each factor in the design.  A 1
) or -1 indicates a factor is present in the alias, with -1 indicating a
) changed sign.  For example, 1 0 1 0 0 -1 means -ACF is a generator
) (alias of I).
)
)) Written by Gary W. Oehlert, gary@stat.umn.edu
)) Version of 000430
# usage: $S(basis)
@basis <- argvalue($1,"basis","integer matrix")
if (max(abs(vector(@basis))) > 1){
    error("basis not matrix with elements -1, 0 and 1")
}
@k <- ncols(@basis)
@p <- nrows(@basis)
if(@k > 25) {
    error("$S() doesn't work with more than 25 factors", macroname:F)
}

@letters <- vecread(string:"ABCDEFGHJKLMNOPQRSTUVWXYZ",bychar:T)[run(@k)]

 # set up all aliases of I
@Ialiases <- matrix(rep(0,@k*(2^@p - 1)),@k)
@counter <- rep(0,@p)
for(@i,run(2^@p-1)) {
    @D <- @counter[1] <- @counter[1]+1
    for(@l,run(@p)) {
        if(@counter[@l] <= 1) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @which <- vector(sum(@counter*abs(@basis))) %% 2
    if((sum(vector(@counter*@basis) < 0) %% 2) != 0) {
        @which[min(run(@k)[@which!= 0])] <- -1
    }
    @Ialiases[,@i] <- @which
}
 #@Ialiases

 #find seed variables for included factorial
@seedvars <- run(@k-@p)
for(@i,run( prod(run(@k))/prod(run(@p))/prod(run(@k-@p)) )) {
    @tmp <- rep(0,@k);@tmp[@seedvars] <- 1
    @tmp <- sum(abs( (1-@tmp) * @Ialiases))
    if(min(vector(@tmp)) > 0) {
        break
    }
    # next seeds
    if(@seedvars[@k-@p] == @k ) { # if end at max, find last not at max
        for(@j,run(@k-@p-1,1)) {
            if(@seedvars[@j] < @p+@j) {
                break;
            }
        }
        @D <- @seedvars[@j] <- @seedvars[@j]+1 # bump it up
        for(@m,run(@j+1,@k-@p)) {
            @D <- @seedvars[@m] <- @seedvars[@j]+@m-@j
        }
		@delete(@m)
    } else {
        @D <- @seedvars[@k-@p] <- @seedvars[@k-@p]+1
    }
}
 #@seedvars

 # do the included factorial in standard order
@counter <- rep(0,@k-@p)
@this <- "I"
for(@i,run(2^@p-1)) {
    @tmp <- paste(@letters[abs(@Ialiases[,@i])>0],sep:"")
    if((sum(@Ialiases[,@i] < 0) %% 2) != 0) {
        @tmp <- paste("-",@tmp,sep:"")
    }
    @this <- paste(@this,@tmp,sep:" = ")
}
@out <- @this

for(@i,run(2^(@k-@p)-1)) {
    @counter[1] <- @counter[1]+1
    for(@l,run(@k-@p)) {
        if(@counter[@l] <= 1) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @vec <- rep(0,@k)
    @vec[ (@seedvars[@counter>0])] <- 1
    @this <- paste(@letters[@vec>0],sep:"")
    for(@j,run(2^@p-1)) {
        @tmp <- paste(@letters[((@vec+abs(@Ialiases[,@j])) %% 2) >0],sep:"")
        if((sum(@Ialiases[,@j] < 0) %% 2) != 0) {
            @tmp <- paste("-",@tmp,sep:"")
        }
        @this <- paste(@this,@tmp,sep:" = ")
    }
    @out <- vector(@out,@this)
}
delete(@D,@this,@vec,@seedvars,@i,@j,@l, @k, @p, @counter,@Ialiases,@tmp,\
	@letters, @basis)
delete(@out, return:T)
%allaliases2%

====> ffdesign2 <====
ffdesign2    MACRO  DOLLARS OUTLINE
) ffdesign2(basis) finds the set of factor/level combinations to
) run in the 2^(k-p) fractional factorial corresponding to the given
) generators.   The p x k matrix basis contains the
) generators for the aliasing, one row for each generator and one
) column for each factor in the design.  The elements are 0, -1, and 1.
) A 1 indicates a factor is present in the alias, signs change sign; eg.
) 1 0 1 0 0 -1   means -ACF is a generator (alias of I).
) Results are returned in a character vector.
# usage: $S(basis)
@basis <- argvalue($1,"basis","integer matrix")
if (max(abs(vector(@basis))) > 1){
    error("basis not matrix with elements -1, 0 and 1")
}
@k <- ncols(@basis)
if(@k > 25) {
    error("$S() doesn't work with more than 25 factors", macroname:F)
}
@p <- nrows(@basis)

@letters <- vecread(string:"abcdefghjklmnopqrstuvwxyz",bychar:T)[run(@k)]

 # set up all aliases of I
@Ialiases <- matrix(rep(0,@k*(2^@p - 1)),@k)
@counter <- rep(0,@p)
for(@i,run(2^@p-1)) {
    @D <- @counter[1] <- @counter[1]+1
    for(@l,run(@p)) {
        if(@counter[@l] <= 1) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @which <- vector(sum(@counter*abs(@basis)))  %% 2
    if((sum(vector(@counter*@basis) < 0) %% 2) != 0) {
        @which[min(run(@k)[@which!= 0])] <- -1
    }
    @D <- @Ialiases[,@i] <- @which
}
 #@Ialiases

 #find seed variables for included factorial
@seedvars <- run(@k-@p)
for(@i,run( prod(run(@k))/prod(run(@p))/prod(run(@k-@p)) )) {
    @tmp <- rep(0,@k);@D <- @tmp[@seedvars] <- 1
    @tmp <- sum(abs( (1-@tmp) * @Ialiases))
    if(min(vector(@tmp)) > 0) {
        break
    }
    # next seeds
    if(@seedvars[@k-@p] == @k ) { # if end at max, find last not at max
        for(@j,run(@k-@p-1,1)) {
            if(@seedvars[@j] < @p+@j) {
                break;
            }
        }
        @D <- @seedvars[@j] <- @seedvars[@j]+1 # bump it up
        for(@m,run(@j+1,@k-@p)) {
            @D <- @seedvars[@m] <- @seedvars[@j]+@m-@j
        }
		delete(@m)
    } else {
        @D <- @seedvars[@k-@p] <- @seedvars[@k-@p]+1
    }
}
 #@seedvars

 # cut down @Ialiases so that we only leave other factors in
 # terms of the seedvars, ie, the generators

@tmp <- vector(sum(abs(@Ialiases[-@seedvars,]))==1)
@Ialiases <- @Ialiases[,@tmp]

 # find the generated factors for each generator
@genfac <- rep(0,@p)
for(@i,run(@p)) {
    @D <- @genfac[@i] <- (run(@k)[-@seedvars])[vector(@Ialiases[-@seedvars,@i]) != 0]
}


 # do the included factorial in standard order
@counter <- rep(0,@k-@p)
@this <- "(1)"
for(@i,run(@p)) {
    @tmp <- 2*(@counter-.5)*@Ialiases[@seedvars,@i]
    @tmp <- prod(@tmp[@tmp!=0])/ @Ialiases[@genfac[@i],@i]
    if(@tmp > 0) {
        if(@this == "(1)") {
            @this <- @letters[@genfac[@i]]
        } else {
            @this <- paste(@this,@letters[@genfac[@i]],sep:"")
        }
    }
}
@out <- @this

for(@i,run(2^(@k-@p)-1)) {
    @D <- @counter[1] <- @counter[1]+1
    for(@l,run(@k-@p)) {
        if(@counter[@l] <= 1) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @vec <- rep(0,@k)
    @D <- @vec[ (@seedvars[@counter>0])] <- 1
    @this <- paste(@letters[@vec>0],sep:"")
    for(@j,run(@p)) {
        @tmp <- 2*(@counter-.5)*@Ialiases[@seedvars,@j]
        @tmp <- prod(@tmp[@tmp!=0])/ @Ialiases[@genfac[@j],@j]
        if(@tmp>0) {
            @this <- paste(@this,@letters[@genfac[@j]],sep:"")
        }
    }
    @out <- vector(@out,@this)
}
delete(@D,@i,@j,@l,@counter,@vec,@this,@tmp,@Ialiases,@genfac,@letters,\
	@k, @p, @seedvars, @which, @basis)
delete(@out, return:T)
%ffdesign2%


====> confound2 <====
confound2        MACRO  DOLLARS OUTLINE
) confound2(basis) confounds the two series factorial into blocks based on
) the generators given in the matrix basis.  The p x k matrix basis contains
) the generators for the confounding, one row for each generator and one
) column for each factor in the design.  The elements are 0 and 1.  A 1
) indicates a factor is present in the generator.  The results are returned
) in a structure with component names block1, block2, etc.  Each component
) has a character vector of factor/level combinations for that block.
# usage: $S(basis)
@basis <- argvalue($1,"basis","nonnegative integer matrix")
if (max(vector(@basis)) > 1){
    error("basis not matrix with elements 0 and 1")
}
@k <- ncols(@basis)
if(@k > 25) {
    error("$S() doesn't work with more than 25 factors",macroname:F)
}
@p <- nrows(@basis)

@letters <- vecread(string:"abcdefghjklmnopqrstuvwxyz",bychar:T)[run(@k)]

 # do the factorial in standard order
 # first transpose basis to make things easier
@basis <- @basis'
@powers <- 2^run(0,@p-1)

@combos <- "(1)"
@blocks <- 1

@counter <- rep(0,@k)
for(@i,run(2^(@k)-1)) {
    @D <- @counter[1] <- @counter[1]+1
    for(@l,run(@k)) {
        if(@counter[@l] <= 1) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @blocks <- vector(@blocks,1+sum((sum(@counter*@basis)' %%2)*@powers))
    @combos <- vector(@combos,paste(@letters[@counter>0],sep:""))
}
@out <- structure(1)
for(@i,run(2^@p)) {
    @out <- changestr(@out,@i,<<paste("block",@i,sep:"")>>:@combos[@blocks==@i])
}
delete(@combos,@blocks,@i, @l,@k, @p,@D,@counter,@powers,@basis)
delete(@out,return:T)
%confound2%

====> confound3 <====
confound3        MACRO  DOLLARS OUTLINE
) confound3(basis) confounds the three series factorial into blocks based on
) the generators given in the matrix basis.  The p x k matrix basis contains
) the generators for the confounding, one row for each generator and one
) column for each factor in the design.  The elements are 0, 1, and 2,
) indicating the exponent of each factor in the generator. The results are
) returned in a structure with component names block1, block2, etc.  Each
) component has a character vector of factor/level combinations for that
) block.
# usage: $S(basis)
@basis <- argvalue($1,"basis","nonnegative integer matrix")
if (max(vector(@basis)) > 2){
    error("basis not matrix with elements 0, 1 and 2")
}
@k <- ncols(@basis)
if(@k > 25) {
    error("$S() doesn't work with more thatn 25 factors", macroname:F)
}
@p <- nrows(@basis)

 # do the factorial in standard order
 # first transpose basis to make things easier
@basis <- @basis'
@powers <- 3^run(0,@p-1)

@combos <- paste(rep(0,@k),sep:"")
@blocks <- 1

@counter <- rep(0,@k)
for(@i,run(3^(@k)-1)) {
    @D <- @counter[1] <- @counter[1]+1
    for(@l,run(@k)) {
        if(@counter[@l] <= 2) {
            break
        }
        @D <- @counter[@l] <- 0
        @D <- @counter[@l+1] <- @counter[@l+1]+1
    }

    @blocks <- vector(@blocks,1+sum((sum(@counter*@basis)' %% 3)*@powers))
    @combos <- vector(@combos,paste(@counter,sep:""))
}
@out <- structure(1)
for(@i,run(3^@p)) {
    @out <- changestr(@out,@i,<<paste("block",@i,sep:"")>>:@combos[@blocks==@i])
}
delete(@D,@i,@l,@combos,@blocks,@basis,@powers,@counter,@k,@p)
delete(@out, return:T)
%confound3%

===> randsign <===
randsign   MACRO  DOLLARS OUTLINE
) randsign(diffs) computes sum(s_i * diffs_i) for all 2^length(diffs)
) possible combinations of signs s_i.  The form randsign(diffs,trials:n)
) simulates a distribution of sum(s_i * diffs_i) based on n sets of random
) signs.  This is appropriate when diffs is long, as 2^length(diffs) grows
) quickly!!
) Results are returned in a vector.
) Version 000129
# usage: $S(diffs)
@diffs <- argvalue($1, "diffs", "real vector")
if (anymissing(@diffs)){
	print(paste("WARNING:",\
		sum(ismissing(@diffs)),"MISSING values in argument to $S() removed"))
	@diffs <- @diffs[!ismissing(@diffs)]
	if (isnull(@diffs)){
		error("all value are missing in argument")
	}
}
@ln <- length(@diffs)

@trials <- keyvalue($K, "trial*", "positive count")
if (!isnull(@trials)){
	@result <- rep(0,@trials)
	for(@i,run(@trials)) {
		@tmp <- 2*((runi(@ln) - 0.5) <= 0) - 1
		@D <- @result[@i] <- sum(@tmp*@diffs)
	}
	delete(@tmp)
}else{
	@plusminus <- vector(1,-1)'
    @result <- vector(@plusminus * @diffs[1])
    for(@i,run(2,@ln)) {
        @result <- vector(@result+@plusminus * @diffs[@i])
    }
	delete(@plusminus)
}
delete(@diffs, @i, @ln, @trials)
delete(@result,return:T)
%randsign%

===> randt <===
randt     MACRO  DOLLARS OUTLINE
) randt(x,m) computes xbar_1 - xbar_2 for all combinations with m data
) values from x in group 1 and the remainder n - m values in group 2,
) where n = length(x).
)
) The result is a vector of length n!/(m!(n-m)!)
)
) x must be a REAL vector with no MISSING values with length(x) > m.
)
) randt(x,m,trials:N) simulates the permutation distribution of
) xbar_1 - xbar_2 based on N random sets of size m and the complementary
) subset of size n - m.
)
) The result is a vector of length N
)) Version 000203
# usage: $S(x,m)
@x <- argvalue($1,"x","real vector nonmissing")
@n1 <- argvalue($2,"m","positive count")
@ln <- length(@x)
@n2 <- @ln - @n1
if(@n2 <= 0) {
    error("m not a positive integer < length(x)")
}
@trials <- keyvalue($K, "trial*", "positive count")
if (!isnull(@trials)){ # use random sampling?
	@result <- rep(0,@trials)
	for(@i,run(@trials)) {
		@J <- grade(runi(@ln)) <= @n1
		@result[@i] <- sum(@x[@J])/@n1 - sum(@x[!@J])/@n2;;
    }
}else{
    @N <- prod(run(@ln))/(prod(run(@n1))*prod(run(@n2)))
    @result <- rep(0,@N)
    @J <- run(@n1)
    @D <- @result[1] <- sum(@x[@J])/@n1-sum(@x[-@J])/@n2
    for(@i,run(2,@N)) {
        # advance J
        if(@J[@n1] == @ln) { # if end at max, find last not at max
            for(@j,run(@n1-1,1)) {
                if(@J[@j] < @n2+@j) {
                    break;
                }
            }
            @D <- @J[@j] <- @J[@j]+1 # bump up one
            for(@j2,run(@j+1,@n1)) {
                @D <- @J[@j2] <- @J[@j] + @j2 - @j
            }
        } else {
            @D <- @J[@n1] <- @J[@n1] + 1
        }

        @D <- @result[@i] <- sum(@x[@J])/@n1 - sum(@x[-@J])/@n2
    }
	delete(@N,@j,@j2,@D)
}
delete(@x, @n1, @n2, @ln, @trials, @i, @J)
delete(@result, return:T)
%randt%

===> randt2 <===
randt2     MACRO  DOLLARS OUTLINE
) randt2(y1, y2), where y1 and y2 are REAL vectors, computes all values
) xbar_1 - xbar_2 where xbar_1 is the mean of a subset of vector(y1, y2) of
) size n1 = length(y1) and xbar_2 is the mean of the complementary subset
) of size n2 = length(y2)
)
) Any MISSING values are removed from y1 and y2 before determining
) n1, n2 and doing the computation.
)
) The result is a REAL vector of length (n1+n2)!/(n1!n2!).
)
) randt2(y1, y2, trials:N) randomly sample N values from the complete
) permutation distribution of xbar_1 - xbar_2.
) The result is a vector of length N
)) Version 000203
# usage: $S(x1,x2 [,trials:N]), x1, x2 REAL vectors, N > 0 integer
if ($v != 2 || $k > 1){
	error("usage: $S(x1,x2 [,trials:N]), x1, x2 REAL vectors, N > 0 integer",\
		macroname:F)
}
@x1 <- argvalue($1,"x1","real vector")
@x2 <- argvalue($2,"x2","real vector")
if(anymissing(@x1)){
	@x1 <- @x1[!ismissing(@x1)]
	if (isnull(@x1)){
		error("x1 has no non-MISSING elements")
	}
}
if(anymissing(@x2)){
	@x2 <- @x2[!ismissing(@x2)]
	if (isnull(@x2)){
		error("x2 has no non-MISSING elements")
	}
}

if ($k > 0){
	if (isnull(keyvalue($K, "trial*", "positive count"))){
		error(paste("unrecognized keyword",compnames($K)))
	}
}
if (!ismacro(randt)){
	getmacros(randt, quiet:T,printname:F)
	if (isnull(randt)){
		error("cannot proceed without macro randt")
	}
}
@n1 <- length(@x1)
@x <- vector(delete(@x1,return:T),delete(@x2,return:T))
@result <- if($k == 0){
	randt(@x, @n1)
}else{
	randt(@x,@n1,$K)
}
delete(@x,@n1)
delete(@result, return:T)
%randt2%

===> rscanon <===
rscanon      MACRO DOLLARS
) rscanon(y,x1,x2,...,xk [,block:f1,block:f2,...block:f3]) performs the
) canonical analysis for the quadratic response surface model with
) response y and predictors x1, ... ,xk.
)
) The output is structure(b0, b, B, x0, y0, H, lambda), where
) the components are the intercept, linear coeffs, the quadratic/cross
) product coefficient matrix, the stationary point, the predicted response
) at the stationary point, the matrix of canonical directions and the
) eigenvalues respectively
)
) Keyword phrases of the form 'block:fj'
) Factors f1, f2, ... specified by any keywords phrases of the form
) block:fj are included in the model as blocking variables.
) Modified 3-18-96 by GWO to add blocks
)) Version 000131, strip $$, use argvalue, cosmetic changes
)) Version 020322, pushmodel() and popmodel() used to preserve
)) GLM information
# usage: $S(y,x1,x2,...,xk [,block:f1,block:f2,...block:f3])
if ($v < 2){
	error("$S() must have at least 2 REAL vector arguments", macroname:F)
}
@y <- argvalue($1,"$1","nonmissing real vector")
@n <- length(@y)
@xnms <- $A[-1]
for(@i,run(1,$v-1)) {
    @tmp <- argvalue(<<@xnms[@i]>>,@xnms[@i],"nonmissing real vector")
	if(length(@tmp) != @n){
        error(paste("x-variable", @xnms[@i], "not the same length as $1"))
    }
}
@allxvars <- hconcat($V)[,-1]
@nx <- $v-1

 #@yname <- paste(paste("@Y",$$,sep:""),$$,sep:"_")
@yname <- paste("@Y","_$$",sep:"")
<<@yname>> <- @y
@model <- paste(@yname,"=1")
if($k > 0) { # take care of blocks
	@keys <- structure($K)
	@keynames <- compnames(@keys)
	for(@i,1,ncomps(@keys)){
		if (match("block*",@keynames[@i],0,exact:F) != 0){
			@bname <- paste("@B",@i,"_$$",sep:"")

			<<@bname>> <- argvalue(@keys[@i], paste("blocking variable",@i),\
								   "positive integer vector")
			if (length(<<@bname>>) != @n || !isfactor(<<@bname>>)){
				error(paste("blocking variable",@i,"not a factor the same length as y"))
			}
			@model <- paste(@model,@bname,sep:"+")
		}else{
			error(paste("unrecognized keyword '",@keynames[@i],"'",sep:""))
		}
	}
}

 #linear terms
for(@i,run(@nx)) {
    @xname <- paste("@X",@i,"_$$",sep:"")
    <<@xname>> <- @allxvars[,@i]
    @model <- paste(@model,@xname,sep:"+")
}

 #quadratic terms
for(@i,run(@nx)) {
	@xname <- paste("@X",@i,"sq_$$",sep:"")
    <<@xname>> <- @allxvars[,@i]^2
    @model <- paste(@model,@xname,sep:"+")
}

 # cross products
for(@i,run(@nx-1)) {
    for(@j,run(@i+1,@nx)) {
		@xname <- paste("@X",@i,"X",@j,"_$$",sep:"")
        <<@xname>> <- @allxvars[,@i]*@allxvars[,@j]
        @model <- paste(@model,@xname,sep:"+")
    }
}
delete(@xname,@allxvars)
@canpush <- alltrue(isfunction(pushmodel),pushmodel(canpush:T))
if (@canpush){
	pushmodel()
}
anova(@model,silent:T)
@info <- modelinfo(aliased:T,coefs:T)
if (delete(@canpush,return:T)){
	popmodel()
}
@aliased <- @info$aliased
@coefs <- @info$coefs
@tmp <- length(@coefs)
@use <- vector(1,run(@tmp+1-(@nx+@nx*(@nx+1)/2),@tmp))
if(prod(@aliased[@use])==1) {
	error("singular columns in $S()", macroname:F)
}
@coefs <- @coefs[@use]
delete(@use,@tmp,@aliased,@info,@model)

@b0 <- @coefs[1]
@b <- @coefs[run(2,@nx+1)]
@B <- dmat(@coefs[run(@nx+2,2*@nx+1)])
@offset <- 1
for(@i,run(@nx-1)) {
    for(@j,run(@i+1,@nx)) {
        @D <- @B[@i,@j] <- @B[@j,@i] <- @coefs[2*@nx+@offset+1]/2
        @offset <-+ 1
    }
}
delete(@i,@j,@nx,@D,@coefs,@offset)
@x0 <- solve(@B,-@b)/2
@y0 <- @b0 + sum(@x0*@b)/2
@eig <- eigen(@B)
structure(b0:delete(@b0,return:T),b:delete(@b,return:T),\
	B:delete(@B,return:T),x0:delete(@x0,return:T),y0:delete(@y0,return:T),\
	H:@eig[2],lambda:@eig[1])
%rscanon%

===> boxcoxvec <===
boxcoxvec   MACRO  DOLLARS
) Macro to compute the error SS for an ANOVA of Box-Cox transformed
) response.
) Usage:
)  result <- boxcoxvec(Model [, powers:pows])
)  result <- boxcoxvec(Rhs_model, y [,powers:pows])
)   Model            CHARACTER scalar specifying a model that
)                    can be used as first argument to anova(), say
)                    "y=a+b"
)   Rhs_model        CHARACTER scalar specifying r.h.s. of model, say
)                    "a+b"
)   y                REAL vector whose Box-Cox transform is to be
)                    on the l.h.s. of the model
)   pows             REAL vector with no MISSING elements; default
)                    is run(-1,2,.25)
)   result           structure(power:pows, SS:ss), ss is a REAL vector
)                    the same length as pows with ss[i] the error SS when
)                    boxcox(y,pows[i]) is the response variable.
)
)) Written by G. Oehlert (gary@stat.umn.edu) and C. Bingham
))  (kb@stat.umn.edu)
)) 020320 Usage boxcoxvec("y=x" [powers:pows]) now legal
))        GLM model information preserved by pushmodel()
# usage: $S(rhs_model,y [,powers:pows]) or $S(model [,powers:pows])
@rhs <- argvalue($1, "rhs_model or model", "string")
if (match("*=*", @rhs, 0, exact:F) != 0){
	@y <- modelvars(y:T, @rhs)
	@modchars <- vecread(string:@rhs,bychar:T)
	@rhs <- paste(@modchars[-run(match("=",@modchars))],sep:"")
	delete(@modchars)
} else {
	@y <- argvalue($2, "y", "real vector")
}
@SS <- @pow <- keyvalue($K, "pow*", "real vector nonmissing",\
	default:run(-1,2,.25))

@canpush <- alltrue(isfunction(pushmodel),pushmodel(canpush:T))

if (@canpush){
	pushmodel()
}
for(@i,run(length(@pow))) {
    @BOXCOXY <- boxcox(@y,@pow[@i])
    anova(paste(nameof(@BOXCOXY),@rhs,sep:"="),silent:T)
    @SS[@i] <- SS[length(SS)]
}
if (delete(@canpush, return:T)){
	popmodel()
}
delete(@rhs, @y, @i, @BOXCOXY)
structure(power:delete(@pow, return:T),SS:delete(@SS, return:T))
%boxcoxvec%

====> all3anova <====
all3anova        MACRO DOLLARS
) macro to fit all hierarchical anova models with 3 predictor variables
) usage: all3anova(y, a, b, c [,mse:val])  where a, b, and c are factors
) and mse is an optional value to use for the error variance
) Adapted from the macanova macro all3() by Kit Bingham which was
) adapted from the S function all3() by Kinley Larntz
) You may add the argument keep:T; in that case, the model df,
) Cp, adjusted R^2, R^2, and model are returned in a structure.
) You may suppress the printing with print:F; keep:T implies
) print:F, but you may use both print:T and keep:T if you want
) both printing and returned values.
)) Written Jan 1996 by GWO
)) Version of 981210: Bug fix
)) Version of 000129 remove $$, use argvalue(), keyvalue()
)) Version of 020322, pushmodel(), popmodel() used to preserve
))  current GLM information
)) Version of 021227, incorrect placement of keyword
# usage: $S(y, a, b, c [,mse:val])
@usage <- "$S(y, a, b, c [,mse:val])"
if ($v != 4){
	error(paste("usage:",@usage), macroname:F)
}
@y <- argvalue($1,"y","real vector")
@a <- argvalue($2,"$2","positive integer vector")
@b <- argvalue($3,"$3","positive integer vector")
@c <- argvalue($4,"$4","positive integer vector")
if (!isfactor(@a) || !isfactor(@b) || !isfactor(@c)){
	error("factor arguments not all factors")
}
@N <- length(@y)
if (length(@a) != @N || length(@b) != @N || length(@c) != @N){
	error("not all factors are the same length as y")
}
@mse <- keyvalue($K,"mse","positive number")
if (isnull(@mse)){
	@mse <- ?
}
@keepit <- keyvalue($K, "keep", "TF", default:F)
@printit <- keyvalue($K, "print", "TF", default:!@keepit)

@canpush <- alltrue(isfunction(pushmodel),pushmodel(canpush:T))

if (@canpush){
	pushmodel()
}

if(ismissing(@mse)) {
	anova(paste(nameof(@y),"=",nameof(@a),".",nameof(@b),".",nameof(@c),\
			sep:""), silent:T)
    if(DF[3] == 0 || SS[3] == 0) {
		error("no MSE available for use in Cp")
	}
    @mse <- SS[3]/DF[3]
}
@mdlnames <- vector("a","b","c","a+b","a+c","b+c",\
    "a*b","a*c","b*c",\
    "a+b+c","a+b*c","b+a*c",\
    "c+a*b","a*b+a*c","a*b+b*c",\
    "a*c+b*c","a*b+a*c+b*c",\
    "a*b*c")

 # '$$' in the following are needed because they are quoted
@models <- vector("@a$$","@b$$","@c$$","@a$$+@b$$","@a$$+@c$$","@b$$+@c$$",\
    "@a$$.@b$$","@a$$.@c$$","@b$$.@c$$",\
    "@a$$+@b$$+@c$$","@a$$+@b$$*@c$$","@b$$+@a$$*@c$$",\
    "@c$$+@a$$*@b$$","@a$$*@b$$+@a$$*@c$$","@a$$*@b$$+@b$$*@c$$",\
    "@a$$*@c$$+@b$$*@c$$","@a$$*@b$$+@a$$*@c$$+@b$$*@c$$",\
    "@a$$*@b$$*@c$$")

@nmodels <- length(@models)
@errordf <- @Cp <- @SSE <- @R2 <- @adjR2 <- rep(0,@nmodels)
for (@i, run(@nmodels)){
	anova(paste(nameof(@y),@models[@i],sep:"="),silent:T)
    @errorterm <- length(SS)
    @tmp <- @errordf[@i] <- DF[@errorterm]
    @tmp <- @SSE[@i] <- SS[@errorterm]
}
@Cp <-  @SSE/@mse + @N - 2*@errordf
@R2 <- 1-@SSE/sum(SS[-1])
@adjR2 <- 1-(@N-1)/(@errordf)*(1-@R2)
@order <- grade(@Cp)
@mdlnames <- @mdlnames[@order]
@modeldf <- @N - @errordf[@order]
@Cp <- @Cp[@order]
@R2 <- @R2[@order]
@adjR2 <- @adjR2[@order]
if (delete(@canpush,return:T)){
	popmodel()
}
if(delete(@printit, return:T)) {
    print(paste(charwidth:4,"   p",\
		charwidth:10," C(p)","   Adj R^2","   R^2","   Model"))
    for(@i,run(@nmodels)){
        print(paste(intwidth:4,format:"10.5f",charwidth:14,\
            @modeldf[@i],@Cp[@i]+1e-12,@adjR2[@i]+1e-12,\
            @R2[@i]+1e-12,@mdlnames[@i]))
	}
}
if(!delete(@keepit,return:T)) {
	delete(@modeldf, @Cp, @adjR2, @R2, @mdlnames)
} else {
    structure(p:delete(@modeldf,return:T),cp:delete(@Cp,return:T),\
		adjrsq:delete(@adjR2,return:T),rsq:delete(@R2,return:T),\
		model:delete(@mdlnames,return:T))
}
%all3anova%

====> all4anova <====
all4anova        MACRO  DOLLARS
) macro to fit all hierarchical anova models with 4 predictor variables
) usage: all4anova(y, a, b, c, d [,mse:val])  where a, b, c and d are factors
) and mse is an optional value to use for the error variance
) Adapted from the macanova macro all4() by Kit Bingham which was
) adapted from the S function all4() by Kinley Larntz
) You may add the argument keep:T; in that case, the model df,
) Cp, adjusted R^2, R^2, and model are returned in a structure.
) You may suppress the printing with print:F; keep:T implies
) print:F, but you may use both print:T and keep:T if you want
) both printing and returned values.
) Written Jan 1996 by GWO
)) Version of 000129 remove $$, use argvalue(), keyvalue()
# usage: $S(y, a, b, c, d [,mse:val])
@usage <- "$S(y, a, b, c, d [,mse:val])"
if ($v != 5){
	error(paste("usage: ",@usage, macroname:F))
}
@y <- argvalue($1,"y","real vector")
@a <- argvalue($2,"$2","positive integer vector")
@b <- argvalue($3,"$3","positive integer vector")
@c <- argvalue($4,"$4","positive integer vector")
@d <- argvalue($5,"$5","positive integer vector")
if (!isfactor(@a) || !isfactor(@b) || !isfactor(@c) || !isfactor(@d)){
	error("factor arguments not all factors")
}

@N <- length(@y)
if (length(@a) != @N || length(@b) != @N || length(@c) != @N ||\
	length(@d) != @N){
	error("not all factors are the same length as y")
}

@mse <- keyvalue($K, "mse", "positive number")
if (isnull(@mse)){
	@mse <- ?
}
@keepit <- keyvalue($K,"keep","TF",default:F)
@printit <- keyvalue($K,"print","TF",default:!@keepit)

@canpush <- alltrue(isfunction(pushmodel),pushmodel(canpush:T))

if (@canpush){
	pushmodel()
}

if(ismissing(@mse)) {
	anova(paste(nameof(@y),"=",nameof(@a),".",nameof(@b),".",nameof(@c),\
			".", nameof(@d), sep:""), silent:T)
    if(DF[3] == 0 || SS[3] == 0) {
		error("no MSE available for use in Cp")
	}
    @mse <- SS[3]/DF[3]
}

 # '$$' are needed because they are quoted
@models <- vector("@a$$","@b$$","@c$$","@d$$",\
    "@a$$+@b$$","@a$$+@c$$","@b$$+@c$$","@a$$+@d$$","@b$$+@d$$","@c$$+@d$$",\
    "@a$$.@b$$","@a$$.@c$$","@b$$.@c$$","@a$$.@d$$","@b$$.@d$$","@c$$.@d$$",\
    "@a$$+@b$$+@c$$","@a$$+@b$$*@c$$","@b$$+@a$$*@c$$",\
    "@c$$+@a$$*@b$$","@a$$*@b$$+@a$$*@c$$","@a$$*@b$$+@b$$*@c$$",\
    "@a$$*@c$$+@b$$*@c$$","@a$$*@b$$+@a$$*@c$$+@b$$*@c$$",\
    "@a$$*@b$$*@c$$",\
    "@d$$+@a$$+@b$$","@d$$+@a$$*@b$$","@a$$+@d$$*@b$$",\
    "@b$$+@d$$*@a$$","@d$$*@a$$+@d$$*@b$$","@d$$*@a$$+@a$$*@b$$",\
    "@d$$*@b$$+@a$$*@b$$","@d$$*@a$$+@d$$*@b$$+@a$$*@b$$",\
    "@d$$*@a$$*@b$$",\
    "@d$$+@a$$+@c$$","@d$$+@a$$*@c$$","@a$$+@d$$*@c$$",\
    "@c$$+@d$$*@a$$","@d$$*@a$$+@d$$*@c$$","@d$$*@a$$+@a$$*@c$$",\
    "@d$$*@c$$+@a$$*@c$$","@d$$*@a$$+@d$$*@c$$+@a$$*@c$$",\
    "@d$$*@a$$*@c$$",\
    "@d$$+@b$$+@c$$","@d$$+@b$$*@c$$","@b$$+@d$$*@c$$",\
    "@c$$+@d$$*@b$$","@d$$*@b$$+@d$$*@c$$","@d$$*@b$$+@b$$*@c$$",\
    "@d$$*@c$$+@b$$*@c$$","@d$$*@b$$+@d$$*@c$$+@b$$*@c$$",\
    "@d$$*@b$$*@c$$",\
    "@a$$+@b$$+@c$$+@d$$","@a$$*@b$$+@c$$+@d$$",\
    "@a$$*@c$$+@b$$+@d$$","@a$$*@d$$+@b$$+@c$$","@b$$*@c$$+@a$$+@d$$",\
    "@b$$*@d$$+@a$$+@c$$","@c$$*@d$$+@a$$+@b$$",\
    "@a$$*@b$$+@a$$*@c$$+@d$$","@a$$*@b$$+@a$$*@d$$+@c$$",\
    "@a$$*@b$$+@b$$*@c$$+@d$$","@a$$*@b$$+@b$$*@d$$+@c$$",\
    "@a$$*@b$$+@c$$*@d$$","@a$$*@c$$+@a$$*@d$$+@b$$",\
    "@a$$*@c$$+@b$$*@c$$+@d$$","@a$$*@c$$+@b$$*@d$$",\
    "@a$$*@c$$+@c$$*@d$$+@b$$","@a$$*@d$$+@b$$*@c$$",\
    "@a$$*@d$$+@b$$*@d$$+@c$$","@a$$*@d$$+@c$$*@d$$+@b$$",\
    "@b$$*@c$$+@b$$*@d$$+@a$$","@b$$*@c$$+@c$$*@d$$+@a$$",\
    "@b$$*@d$$+@c$$*@d$$+@a$$","@a$$*@b$$+@a$$*@c$$+@a$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@b$$*@c$$+@d$$","@a$$*@b$$+@a$$*@c$$+@b$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@c$$*@d$$","@a$$*@b$$+@a$$*@d$$+@b$$*@c$$",\
    "@a$$*@b$$+@a$$*@d$$+@b$$*@d$$+@c$$","@a$$*@b$$+@a$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@b$$*@c$$+@b$$*@d$$","@a$$*@b$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$+@b$$*@d$$+@c$$*@d$$","@a$$*@c$$+@a$$*@d$$+@b$$*@c$$",\
    "@a$$*@c$$+@a$$*@d$$+@b$$*@d$$","@a$$*@c$$+@a$$*@d$$+@c$$*@d$$+@b$$",\
    "@a$$*@c$$+@b$$*@c$$+@b$$*@d$$","@a$$*@c$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@c$$+@b$$*@d$$+@c$$*@d$$","@a$$*@d$$+@b$$*@c$$+@b$$*@d$$",\
    "@a$$*@d$$+@b$$*@c$$+@c$$*@d$$","@a$$*@d$$+@b$$*@d$$+@c$$*@d$$",\
    "@b$$*@c$$+@b$$*@d$$+@c$$*@d$$+@a$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@b$$*@c$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@b$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@b$$*@c$$+@b$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@d$$+@b$$*@c$$+@b$$*@d$$",\
    "@a$$*@b$$+@a$$*@d$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@d$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@c$$+@a$$*@d$$+@b$$*@c$$+@b$$*@d$$",\
    "@a$$*@c$$+@a$$*@d$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@c$$+@a$$*@d$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@c$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@d$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@b$$*@c$$+@b$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@d$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@c$$+@a$$*@d$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$+@a$$*@c$$+@a$$*@d$$+@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@d$$","@a$$*@b$$*@c$$+@a$$*@d$$",\
    "@a$$*@b$$*@c$$+@b$$*@d$$","@a$$*@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@d$$+@b$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@d$$+@b$$*@d$$+@c$$*@d$$","@a$$*@b$$*@d$$+@c$$",\
    "@a$$*@b$$*@d$$+@a$$*@c$$","@a$$*@b$$*@d$$+@b$$*@c$$",\
    "@a$$*@b$$*@d$$+@c$$*@d$$","@a$$*@b$$*@d$$+@a$$*@c$$+@b$$*@c$$",\
    "@a$$*@b$$*@d$$+@a$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$*@d$$+@b$$*@c$$+@c$$*@d$$",\
    "@a$$*@b$$*@d$$+@a$$*@c$$+@b$$*@c$$+@c$$*@d$$","@a$$*@c$$*@d$$+@b$$",\
    "@a$$*@c$$*@d$$+@a$$*@b$$","@a$$*@c$$*@d$$+@b$$*@c$$",\
    "@a$$*@c$$*@d$$+@b$$*@d$$","@a$$*@c$$*@d$$+@a$$*@b$$+@b$$*@c$$",\
    "@a$$*@c$$*@d$$+@a$$*@b$$+@b$$*@d$$",\
    "@a$$*@c$$*@d$$+@b$$*@c$$+@b$$*@d$$",\
    "@a$$*@c$$*@d$$+@a$$*@b$$+@b$$*@c$$+@b$$*@d$$","@b$$*@c$$*@d$$+@a$$",\
    "@b$$*@c$$*@d$$+@a$$*@b$$","@b$$*@c$$*@d$$+@a$$*@c$$",\
    "@b$$*@c$$*@d$$+@a$$*@d$$","@b$$*@c$$*@d$$+@a$$*@b$$+@a$$*@c$$",\
    "@b$$*@c$$*@d$$+@a$$*@b$$+@a$$*@d$$",\
    "@b$$*@c$$*@d$$+@a$$*@c$$+@a$$*@d$$",\
    "@b$$*@c$$*@d$$+@a$$*@b$$+@a$$*@c$$+@a$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@b$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@b$$*@d$$+@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@c$$*@d$$+@b$$*@d$$",\
    "@a$$*@b$$*@c$$+@b$$*@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@b$$*@c$$*@d$$+@a$$*@d$$",\
    "@a$$*@b$$*@d$$+@a$$*@c$$*@d$$",\
    "@a$$*@b$$*@d$$+@a$$*@c$$*@d$$+@b$$*@c$$",\
    "@a$$*@b$$*@d$$+@b$$*@c$$*@d$$",\
    "@a$$*@b$$*@d$$+@b$$*@c$$*@d$$+@a$$*@c$$",\
    "@a$$*@c$$*@d$$+@b$$*@c$$*@d$$",\
    "@a$$*@c$$*@d$$+@b$$*@c$$*@d$$+@a$$*@b$$",\
    "@a$$*@b$$*@c$$+@a$$*@b$$*@d$$+@a$$*@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@b$$*@d$$+@b$$*@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@c$$*@d$$+@b$$*@c$$*@d$$",\
    "@a$$*@b$$*@d$$+@a$$*@c$$*@d$$+@b$$*@c$$*@d$$",\
    "@a$$*@b$$*@c$$+@a$$*@b$$*@d$$+@a$$*@c$$*@d$$+@b$$*@c$$*@d$$",\
    "@a$$*@b$$*@c$$*@d$$")

@mdlnames <- vector("a","b","c","d",\
    "a+b","a+c","b+c","a+d","b+d","c+d",\
    "a.b","a.c","b.c","a.d","b.d","c.d",\
    "a+b+c","a+b*c","b+a*c",\
    "c+a*b","a*b+a*c","a*b+b*c",\
    "a*c+b*c","a*b+a*c+b*c",\
    "a*b*c",\
    "d+a+b","d+a*b","a+d*b",\
    "b+d*a","d*a+d*b","d*a+a*b",\
    "d*b+a*b","d*a+d*b+a*b",\
    "d*a*b",\
    "d+a+c","d+a*c","a+d*c",\
    "c+d*a","d*a+d*c","d*a+a*c",\
    "d*c+a*c","d*a+d*c+a*c",\
    "d*a*c",\
    "d+b+c","d+b*c","b+d*c",\
    "c+d*b","d*b+d*c","d*b+b*c",\
    "d*c+b*c","d*b+d*c+b*c",\
    "d*b*c",\
    "a+b+c+d","a*b+c+d",\
    "a*c+b+d","a*d+b+c","b*c+a+d",\
    "b*d+a+c","c*d+a+b",\
    "a*b+a*c+d","a*b+a*d+c",\
    "a*b+b*c+d","a*b+b*d+c",\
    "a*b+c*d","a*c+a*d+b",\
    "a*c+b*c+d","a*c+b*d",\
    "a*c+c*d+b","a*d+b*c",\
    "a*d+b*d+c","a*d+c*d+b",\
    "b*c+b*d+a","b*c+c*d+a",\
    "b*d+c*d+a","a*b+a*c+a*d",\
    "a*b+a*c+b*c+d","a*b+a*c+b*d",\
    "a*b+a*c+c*d","a*b+a*d+b*c",\
    "a*b+a*d+b*d+c","a*b+a*d+c*d",\
    "a*b+b*c+b*d","a*b+b*c+c*d",\
    "a*b+b*d+c*d","a*c+a*d+b*c",\
    "a*c+a*d+b*d","a*c+a*d+c*d+b",\
    "a*c+b*c+b*d","a*c+b*c+c*d",\
    "a*c+b*d+c*d","a*d+b*c+b*d",\
    "a*d+b*c+c*d","a*d+b*d+c*d",\
    "b*c+b*d+c*d+a",\
    "a*b+a*c+a*d+b*c",\
    "a*b+a*c+a*d+b*d",\
    "a*b+a*c+a*d+c*d",\
    "a*b+a*c+b*c+b*d",\
    "a*b+a*c+b*c+c*d",\
    "a*b+a*c+b*d+c*d",\
    "a*b+a*d+b*c+b*d",\
    "a*b+a*d+b*c+c*d",\
    "a*b+a*d+b*d+c*d",\
    "a*b+b*c+b*d+c*d",\
    "a*c+a*d+b*c+b*d",\
    "a*c+a*d+b*c+c*d",\
    "a*c+a*d+b*d+c*d",\
    "a*c+b*c+b*d+c*d",\
    "a*d+b*c+b*d+c*d",\
    "a*b+a*c+a*d+b*c+b*d",\
    "a*b+a*c+a*d+b*c+c*d",\
    "a*b+a*c+a*d+b*d+c*d",\
    "a*b+a*c+b*c+b*d+c*d",\
    "a*b+a*d+b*c+b*d+c*d",\
    "a*c+a*d+b*c+b*d+c*d",\
    "a*b+a*c+a*d+b*c+b*d+c*d",\
    "a*b*c+d","a*b*c+a*d",\
    "a*b*c+b*d","a*b*c+c*d",\
    "a*b*c+a*d+b*d",\
    "a*b*c+a*d+c*d",\
    "a*b*c+b*d+c*d",\
    "a*b*c+a*d+b*d+c*d","a*b*d+c",\
    "a*b*d+a*c","a*b*d+b*c",\
    "a*b*d+c*d","a*b*d+a*c+b*c",\
    "a*b*d+a*c+c*d",\
    "a*b*d+b*c+c*d",\
    "a*b*d+a*c+b*c+c*d","a*c*d+b",\
    "a*c*d+a*b","a*c*d+b*c",\
    "a*c*d+b*d","a*c*d+a*b+b*c",\
    "a*c*d+a*b+b*d",\
    "a*c*d+b*c+b*d",\
    "a*c*d+a*b+b*c+b*d","b*c*d+a",\
    "b*c*d+a*b","b*c*d+a*c",\
    "b*c*d+a*d","b*c*d+a*b+a*c",\
    "b*c*d+a*b+a*d",\
    "b*c*d+a*c+a*d",\
    "b*c*d+a*b+a*c+a*d",\
    "a*b*c+a*b*d",\
    "a*b*c+a*b*d+c*d",\
    "a*b*c+a*c*d",\
    "a*b*c+a*c*d+b*d",\
    "a*b*c+b*c*d",\
    "a*b*c+b*c*d+a*d",\
    "a*b*d+a*c*d",\
    "a*b*d+a*c*d+b*c",\
    "a*b*d+b*c*d",\
    "a*b*d+b*c*d+a*c",\
    "a*c*d+b*c*d",\
    "a*c*d+b*c*d+a*b",\
    "a*b*c+a*b*d+a*c*d",\
    "a*b*c+a*b*d+b*c*d",\
    "a*b*c+a*c*d+b*c*d",\
    "a*b*d+a*c*d+b*c*d",\
    "a*b*c+a*b*d+a*c*d+b*c*d",\
    "a*b*c*d")

@nmodels <- length(@models)
@errordf <- @Cp <- @SSE <- @R2 <- @adjR2 <- rep(0,@nmodels)
for (@i, run(@nmodels)){
	anova(paste(nameof(@y),@models[@i],sep:"="), silent:T)
    @errorterm <- length(SS)
    @tmp <- @errordf[@i] <- DF[@errorterm]
    @tmp <- @SSE[@i] <- SS[@errorterm]
}
delete(@tmp,@errorterm,@i,@a,@b,@c,@d,@y)

@Cp <- @SSE/delete(@mse,return:T) + @N - 2*@errordf
@R2 <- 1 - @SSE/sum(SS[-1])
@adjR2 <- 1 - (@N - 1)*(1 - @R2)/@errordf
@order <- grade(@Cp)
@mdlnames <- @mdlnames[@order]
@modeldf <- @N - @errordf[@order]
@Cp <- @Cp[@order]
@R2 <- @R2[@order]
@adjR2 <- @adjR2[@order]
if (@canpush){
	popmodel()
}

delete(@order)
if(delete(@printit,return:T)) {
    print(paste(charwidth:4,"   p",\
		charwidth:10," C(p)","   Adj R^2","   R^2","   Model"))
    for(@i,run(@nmodels)){
        print(paste(intwidth:4,format:"10.5f",charwidth:14,\
            @modeldf[@i],@Cp[@i]+1e-12,@adjR2[@i]+1e-12,\
            @R2[@i]+1e-12,@mdlnames[@i]))
	}
	delete(@i)
}
delete(@nmodels)
if(delete(@keepit,return:T)) {
    structure(p:delete(@modeldf,return:T),cp:delete(@Cp,return:T),\
		adjrsq:delete(@adjR2,return:T),	rsq:delete(@R2,return:T),\
		model:delete(@mdlnames,return:T))
}else{
	delete(@modeldf,@Cp,@adjR2,@R2,@mdlnames)
}
%all4anova%

===> pairedcomp <===
pairedcomp   MACRO
) calls pairwise, the new name for pairedcomp
pairwise($0)
%pairedcomp%

===> pairwise <===
pairwise   MACRO  DOLLARS OUTLINE
) pairwise(termname,lev [,method:T] [,error:term] ) does pairwise comparisons
) of the means of the levels of the factor given in termname at level of
) signficance lev.  If only these two arguments are used, then the pairwise
) comparisons are done using the Bonferroni method.  Other techniques can be
) requested by using a method keyword.  Currently supported methods are: lsd
) (least significant difference), bsd (Bonferroni), snk (Student-Newman-
) Keuls), hsd (Tukey's honest significant difference, also called the
) studentized range procedure), regwb (regw with Bonferroni test), and regwr
) (regw with studentized range).
)
)) Thus, for example, pairwise("trt",.01,hsd:T) does pairwise comparisons
)) between the levels of trt at signficance .01 using the hsd method.
)) termname should be a single factor, not an interaction.
)) lev should be a number between 0 and 1
))
)) The printed output is one row for each level of the term, sorted from
)) smallest to largest effect, giving the "underlines", term number, and
)) effect.
)
) By default, the error in these tests is taken from the last error term of
) the current model (the last line of the ANOVA table).  You may specify a
) different error term with a keyword error:term.  term may be a number,
) indicating a line in the ANOVA table, or the name of a term in the ANOVA
) table, for example, error:4 or error:"ERROR1".
)
) An alternative form is pairwise(termname,critval:val)  In this form, all
) pairwise t-test statistics are compared to the critical value given in val.
) The same output is printed as before.
)
)) Version of 000430 made revisions on handling of temporary macros
# usage: $S(termname,lev [,method:T] [,error:term])
#        or $S(termname,critval:val)
if ($v > 2 || $v == 0 || $v == 2 && $k > 2 || $v == 1 && $k != 1){
	error("usage: $S(termname,lev [,method:T] [,error:term]) or $S(termname,critval:val)",\
		macroname:F)
}

if(!isvector(TERMNAMES) || !ischar(TERMNAMES)){
	error("apparently no previous GLM command")
}
@term <- argvalue($1, "term name", "string")
if (match(@term,TERMNAMES,0) == 0){
	error(paste("\"",@term,"\" is not a term in the current model",sep:""))
}
if (match("*.*", @term, 0, exact:F) > 0){
	error(paste("\"",@term,"\" is an interaction term",sep:""))
}
@cfs <- coefs(@term)
@ord <- grade(@cfs)
@nlvl <- length(@cfs)
if (anymissing(@cfs)) {
    error("coefficients for term \"$1\" contain missing values")
}

if($v > 1) {
    @alpha <- argvalue($2, "level", "positive number")
	if (@alpha >= 1){
        error("level not a scalar between 0 and 1 in $S()", macroname:F)
    }
	if (@alpha >= .5){
		print("WARNING: error rate >= .5", macroname:T)
	}
}else{
	@alpha <- .05
}

@error <- length(TERMNAMES)
if ($k > 0){
	@keynames <- compnames(structure($K))
	for(@i,1,length(@keynames)){
		if (match(@keynames[@i], vector("bsd","lsd","regwb","snk",\
			"regwr","hsd","error","critval"),0) == 0){
			error(paste("unrecognized keyword \"",@keynames[@i],"\"",sep:""))
		}
	}
}


if ($v == 2){
	if (!isnull(keyvalue($K,"critval*"))){
		error("keyword 'critval' cannot be used when level is specified")
	}

	@bsd <- T
	@lsd <- @regwb <- @snk <- @regwr <- @hsd <- F
	if($k > 0) {
		@lsd <- keyvalue($K,"lsd","TF",default:F)
		@regwb <- keyvalue($K,"regwb","TF",default:F)
		@snk <- keyvalue($K,"snk","TF",default:F)
		@regwr <- keyvalue($K,"regwr","TF",default:F)
		@hsd <- keyvalue($K,"hsd","TF",default:F)
		@bsd <- keyvalue($K,"bsd","TF",\
			default:!(@lsd || @regwb || @snk || @regwr || @hsd))

		@error <- keyvalue($K,"error")
	    if (isnull(@error)){
			if($k == 2){
				error("when 2 keywords are used, one must be 'error'")
			}
			@error <- length(TERMNAMES)
		}else{
			@error <- if (ischar(@error)){
				argvalue(@error,"error","string")
			}elseif(isreal(@error)){
				argvalue(@error,"error","positive count")
			}else{
				error("value for 'error' not CHARACTER scalar or positive integer")
			}
			if (ischar(@error)){
				@which <- match(@error,TERMNAMES,-1)
				if (@which < 0){
					error(paste("\",@error,\" not a term in the current model"))
				}
				@error <- delete(@which,return:T)
			}elseif(@error > length(TERMNAMES)){
				error("value for 'error' > 1 + number of terms in current model")
			}
		}
	}

	# arguments for critval(df,alpha,numgps,length_of_stretch)
	@critval <- if (@lsd){
		macro("invstu(1 - \\$2/2, \\$1)", inline:F)
	}elseif (@hsd){
		macro("invstudrng(1 - \\$2, \\$3, \\$1)/sqrt(2)",inline:F)
	}elseif (@regwb){
		macro("if(\\$3 - \\$4 <= 1){
			invstu(1-\\$2/\\$4/(\\$4-1), \\$1)
		} else {
			invstu(1-\\$2/\\$3/(\\$4-1),\\$1)
		}",inline:F)
	}elseif (@regwr){
		macro("if(\\$3 - \\$4 <= 1) {
			invstudrng(1-\\$2, \\$4, \\$1)/sqrt(2)
		} else {
			invstudrng(1-\\$2*\\$4/\\$3,\\$4,\\$1)/sqrt(2)
		}",inline:F)
	}elseif (@snk){
		macro("invstudrng(1-\\$2,\\$4,\\$1)/sqrt(2)",inline:F)
	}else{
		macro("invstu(1-\\$2/\\$3/(\\$3-1),\\$1)",inline:F)
	}
	delete(@bsd,@lsd,@regwb,@snk,@regwr,@hsd)
}elseif ($v == 1){
	@usercrit <- keyvalue($K,"critval*","positive number")
	if (isnull(@usercrit)){
		error("must have 'critval:val' when the level is not specified")
	}
	@critval <- macro(paste(delete(@usercrit,return:T),format:".17g"),inline:F)
}

 # test all pairs and record matches
 # a 1 indicates not significantly different; 0 indicates different
@pairs <- dmat(rep(1,@nlvl))
@lowzero <- 1*(run(@nlvl)-run(@nlvl)' <= 0)

for(@i,run(1,@nlvl-1)) {
    @ti <- @ord[@i];
    for(@j,run(@nlvl,@i+1)) {
        if(@pairs[@i,@j] == 1) { # if already not different
            break
        }
        @tj <- @ord[@j];
        @tmin <- min(@ti,@tj)
        @tmax <- max(@ti,@tj);
        @c <- rep(vector(0,1,0,-1,0),\
					vector(@tmin-1,1,@tmax-@tmin-1,1,@nlvl-@tmax))
        @ctr <- contrast(@term,@c,error:@error)
		@stretch <- @j - @i + 1
        if(abs(@ctr[1]/@ctr[3]) < @critval(DF[@error],@alpha,@nlvl,@stretch)) {
            @pairs[run(@i,@nlvl),run(@j)] <- 1
            @pairs <-* @lowzero
        }
    }
}
delete(@tj,@tmin,@tmax,@c,@ctr,@lowzero,@critval,@term, @alpha)

 # now coalesce
 # remove rows contained in earlier rows
@nrows <- @nlvl
for(@i,run(@nrows,2)) {
    for(@j,run(1,@i-1)) {
        if(sum(abs(vector(@pairs[@i,run(@i,@nlvl)]-\
            @pairs[@j,run(@i,@nlvl)]))) == 0) {
            @pairs <- @pairs[-@i,]
            break
        }
    }
}

 # remove singleton rows
@nrows <- dim(@pairs)[1]
if(@nrows > 1) {
    for(@i,run(@nrows,2)) {
        if(sum(vector(@pairs[@i,])) == 1) {
            @pairs <- @pairs[-@i,]
        }
    }
}

@nrows <- dim(@pairs)[1]
if(sum(vector(@pairs)) == 1) {
    @pairs <- 0
} else {
    if(sum(vector(@pairs[1,])) == 1) {
        @pairs <- @pairs[-1,]
    }
}
delete(@j,@nrows)

 #printout
for(@i,run(1,@nlvl)) {
    @out <- paste(format:"4.0f",@ord[@i],format:"8.3g",@cfs[@ord[@i]])
    if(length(@pairs) > 1) {
        @out <- paste(" ",vector(" ","|")[1+vector(@pairs[,@i])],@out)
    }
    print(@out)
}
delete(@i,@nlvl,@out,@pairs,@cfs,@ord)
%pairwise%

===> quadmax <===
quadmax     MACRO  DOLLARS OUTLINE
) quadmax(A,b[,eq:eqmat][,gte:gtemat][,ckbounds:F])
) quadmax finds the maximum of x'Ax + b'x; if there is no unique maximum
) then quadmax returns NULL A is a p by p real matrix and b is a p by 1
) vector
)
) You may specify linear equality constraints on the solution by using the
) eq:eqmat keyword phrase.  eqmat is an q by p+1 matrix partitioned [Q:y]
) where Q is q by p and y is q by 1.  When this keyword phrase is used,
) the solution is constrained to satisfy Qx = y.  If the constraints
) cannot be met, then quadmax returns NULL.  If the constrained problem is
) unbounded, quadmax returns an error.
)
) You may specify linear inequality constraints on the solution by using
) the gte:gtemat keyword phrase.  gtemat is an g by p+1 matrix partitioned
) [G:z] where G is g by p and z is g by 1.  When this keyword phrase is
) used, the solution is constrained to satisfy Gx >= z elementwise.  If
) the constraints cannot be met, then quadmax returns NULL.  If the
) constrained problem is unbounded, quadmax returns an error.
)
)) For example, in a three variable mixture problem you might have the
)) equality constraint that the sum of the x's is 1 and each element of x
)) is at least .05.  Then use eq:eqmat,gte:gtemat where eqmat is [1 1 1 1]
)) and gtemat is [ 1 0 0 .05 ]
))               [ 0 1 0 .05 ]
))               [ 0 0 1 .05 ]
)
) quadmax tries to determine when the problem is unbounded.  This can
) increase the computational time substantially, so you may tell quadmax
) not to check for boundedness by using the ckbounds:F keyword phrase.
)) This only seems reasonable when the problem is known to be bounded, for
)) example, in a problem where the range of the x's is totally bounded by
)) inequality constraints, or in a problem where the A matrix is known to
)) have a unique maximum.
)) 020322 Instances of a' %*% b changed to a %c% b
# usage: $S(A,b[,eq:eqmat][,gte:gtemat][,ckbounds:F])
 # Read in needed macro
if(!ismacro(quadmaxlin)) {
    getmacros(quadmaxlin,quiet:T,printname:F)
    if(!ismacro(quadmaxlin)) {
        error("cannot find macro quadmaxlin")
    }
}
@A <- argvalue($1, "argument A", "nonmissing real square")
@b <- argvalue($2, "argument b", "nonmissing real vector")

if(nrows(@b) != nrows(@A) ) {
    error("$S() arguments A and b don't have the same number of rows",\
			macroname:F)
}
@p <- nrows(@b)

 # need to process keywords

@tmp2 <- keyvalue($K,"eq","nonmissing real matrix")
if(isnull(@tmp2)) {
    # throw in a default null constraint to make programming easy
    @Q <- rep(0,@p)'
    @q <- 1
    @y <- 0
} else {
    if(ncols(@tmp2) != @p + 1) {
        error("value of 'eq' not REAL matrix with 1 more column than A")
    }
    @Q <- vconcat(@tmp2[,run(@p)])
    @q <- nrows(@Q)
    @y <- vector(@tmp2[,@p+1])
}

@tmp2 <- keyvalue($K,"gte","nonmissing real matrix")
if(isnull(@tmp2)) {
    @g <- 0
    @z <- @G <- NULL
} else {
    if(ncols(@tmp2) != @p + 1) {
        error("value of 'gte' not REAL matrix with 1 more column than A")
    }
    @G <- @tmp2[,run(@p)]
    @g <- nrows(@G)
    @z <- vector(@tmp2[,@p+1])
}
delete(@tmp2)
@ckbounds <- keyvalue($K,"ckbound*","TF", default:T)

@ckbmin <- 2^@g
if(delete(@ckbounds,return:T)) {
    # if we check bounds, we want to put maxima and minima on the
    # directions in which A has positive eigen values and add these
    # to our inequality constraints

    @tmpbig <- 1e50

    @tmpe <- eigen(@A)
    @tmpp <- sum(@tmpe$values > 0)
    if(@tmpp > 0) {
        # positive eigen values
        @Hp <- @tmpe$vectors[,run(@tmpp)]
        @tmp <- -@tmpbig/sqrt(@tmpe$values[run(@tmpp)])
        @G <- vconcat(@G,@Hp',-@Hp')
        @z <- vector(@z,@tmp,@tmp)
        @g <- nrows(@G)
		delete(@tmp,@Hp)
    }
	delete(@tmpp,@tmpe,@tmpbig)
}

 #print(@G,@z)
 # what we do is try to find the best minimally constrained, then add
 # increasing numbers of constraints

@out <- @bestcrit <- @bestset <- @boundsets <- NULL

for(@thispat,run(0,2^@g-1)) {
    @thisset <- (@thispat %& 2^run(0,max(1,@g)-1)) > 0
    # augment equality constraints
    if(anytrue(isnull(@boundsets),\
				sum( (@thispat %& @boundsets) == @boundsets) == 0 )) {
        # do this if we have no prior "interior" solutions or this constraint pattern
        # is not a strict superset of any previous interior solution pattern
        if(@thispat == 0) {
            @agmntdQ <- @Q
            @agmntdy <- @y
        } else {
            @rows <- run(@g)[@thisset]
            @agmntdQ <- vconcat(@Q,@G[@rows,])
            @agmntdy <- vector(@y,@z[delete(@rows,return:T)])
        }
        @thistry <- quadmaxlin(@A,@b,@agmntdQ,@agmntdy);
        if( alltrue(!isnull(@thistry),@g > 0, \
            min(@G %*% @thistry-@z+1e-7*(abs(@G) %*% (rep(1,@p) + abs(@thistry)))) < 0 )){
            # null it out if gt constraints present and not met
            @thistry <- NULL
        }
        # if not null now, we have a candidate solution; compare it to best so far
        if( !isnull(@thistry)) {
            @boundsets <- vector(@boundsets,@thispat)
            @thiscrit <- @thistry %c% @A %*% @thistry + @b %c% @thistry
            if( anytrue(isnull(@out),@thiscrit > @bestcrit )) {
                if(@thispat >= @ckbmin) {
                    error("$S() finds that there is probably an unbounded solution",\
						macroname:F)
                }
                # first candidate or better candidate
                @out <- @thistry
                @bestcrit <- @thiscrit;
				@bestset <- @thispat
            }
        }
    }
}
delete(@thistry,@thispat,@thisset,@bestset,@bestcrit, @ckbmin, @boundsets,\
	@agmntdQ,@agmntdy,@A,@b,@G,@z,@p,@Q,@y,@g)
delete(@out,return:T)
%quadmax%

===> quadmaxlin <===
quadmaxlin  MACRO  DOLLARS
) quadmaxlin(A,b,C,y) finds the maximum of x'Ax + b'x subject to Cx = y
)
) If Cx = y has no solution or the constrained problem has a saddle point
) or a mimumum, then quadmaxlin returns NULL otherwise, quadmaxlin returns
) the x
)) Version 000129 removed $$, use argvalue()
# usage: $S(A,b,C,y)
@A <- argvalue($1,"argument A","real square nonmissing")
@b <- argvalue($2,"argument b","real vector nonmissing")
@C <- argvalue($3,"argument C","real matrix nonmissing")
@y <- argvalue($4,"argument y","real vector nonmissing")

 # p is dimension of the problem
@p <- nrows(@b)

if(ncols(@A) != ncols(@C)) {
    error("arguments A and C must have the same number of columns")
}
if(nrows(@A) != @p ) {
    error("argument b must have the same number of rows as A")
}
if( nrows(@C) != nrows(@y)) {
    error("argument y must have the same number of rows as C")
}

 # check for a solution to Cx = y
@tmp <- hconcat(@C,@y)
@tmp <- @tmp %c% @tmp

@warn <- getoptions(warnings:T)
setoptions(warnings:F)
@tmp <- swp(@tmp,run(@p))
setoptions(warnings:delete(@warn,return:T))

@out <- ?
if( @tmp[@p+1,@p+1] > 1e-10 * sum(@y^2) ) {
    # no solution to constraint
    @out <- NULL
} else {
    # now decompose C
    if(nrows(@C) >= ncols(@C)) {
        @tmp <- svd(@C,all:T)
        @L <- @tmp$leftvectors
        @D <- @tmp$values
        @R <- @tmp$rightvectors
    } else {
        @tmp <- svd(@C',all:T)
        @L <- @tmp$rightvectors
        @D <- @tmp$values
        @R <- @tmp$leftvectors
    }
    # r is the dimension of the constraint
    @r <- sum(@D > 1e-10)
    if(@r == 0) {
        # no constraints
        if( sum(@y^2) != 0) {
            @out <- NULL;  # can't meet constraint
        } else {
            @F <- @A;
            @f <- @b'
        }
	} else {
		@J <- run(@r)
		if (@r == @p) {
			# completely constrained
			@R1 <- @R[,@J];
			@D1 <- @D[@J]
			@L1 <- @L[,@J]
			@out <- (@R1' / @D1) %c% (@L1 %c% @y)
		} else {
			# have some constraints
			@R1 <- @R[,@J]
			@D1 <- @D[@J]
			@L1 <- @L[,@J]
			@tmp <- dmat(rep(1,@p)) - @R1 %*% @R1'
			@R2 <- eigen(@tmp)$vectors[,run(@p-@r)]
			@F <- @R2 %c% @A %*% @R2
			@f <- 2*((@L1 / @D1') %c% @y) %c% @R1' %*% @A %*% @R2 \
				+ @b %c% @R2
		}
		delete(@J)
    }
    if(alltrue(!isnull(@out),anymissing(@out))) {
        @tmp <- eigen(@F)
        if(max(@tmp$values) >= 0) {
            # saddle or minimum
            @out <- NULL;
        } else {
            @out <- solve(@F,-@f'/2)
            if( @r > 0 && @r < @p) { # a constrained problem
				@out <- @R2 %*% @out + ((@L1' / @D1) %c% @R1') %c% @y
            }
        }
    }
}
delete(@tmp,@y,@p,@F,@f,@r,@R2,@R1,@D1,@L1, @b, silent:T)
delete(@out,return:T)
%quadmaxlin%

ems     MACRO  DOLLARS
) ems(Model,Randomvars) computes the expected mean squares for the terms
) in the ANOVA for the model given in CHARACTER scalar Model.  Randomvars
) is a CHARACTER vector specifying the names of factors in the model which
) are random.  Randomvars can also be REAL with integer elements
) specifying the index of a factor in the model.  If there are no random
) factors, Randomvars should be NULL.
)
) In this default use, ems() computes sequential (Type I) sums of squares
) for the the restricted (mixed effects add to zero across fixed factors)
) model, prints these expected means squares for each term, and returns no
) value.
)
) Any variates in the model must be fixed effects.
)
)) ems() assumes that if a factor first appears in an interaction, then
)) that factor is nested in the other terms of the interaction.  For
)) example, if the first appearance of factor c is in the term a.b.c, then
)) c is assumed nested in the a.b combinations.  This nesting is assumed in
)) the remainder of the model.  That is, continuing the example, if there
)) is a later term c.d, it will be interpreted as a.b.c.d even though
)) a.b.c.d is not specifically in the model.
))
)) ems() works for both balanced and unbalanced data.
))
) ems(Model,Randomvars,marg:T) computes expected mean squares based on
) adjusted (Type III) sums of squares.
)
) ems(Model,Randomvars,restrict:F) computes expected mean squares assuming
) no marginal restrictions on the random effects in the model.
)
) ems(Model,Randomvars,nonhier:T) computes expected mean squares for an
) analysis of variance that does not enforce the usual MacAnova hierarchy
) assumptions.
)) That is, for example, model "y=a+b+c+a.b.c" does not imply
)) that the two-way interaction degrees of freedom are part of the "a.b.c"
)) term.  You cannot use anova() to compute such an analysis although it
)) can be done (if you know how) using swp().
))
)) These keywords can be used together.  For example, ems(Model,Randomvars,
)) marg:T,restrict:F) provides answers equivalent to the EMS in SAS PROC GLM.
)
) ems(Model,Randomvars [,...],keep:T) suppresses printed output but
) returns a structure containing the results.  If you want the printed
) output too, use keep:T,print:T.  The structure returned has components
) 'df', 'ss', 'termnames', 'coefs', and 'rterms'.
)
)) Component   Description
))  df         REAL Vector of degrees of freedom for all terms in model
))  ss         REAL Vector of sums of squares for all terms in model
))  termnames  CHARACTER vector of labels for each term
))  coefs      REAL matrix with coefs[i,j] the coefficient for term j in
))             the EMS of term i
))  rterms     LOGICAL vector with T indicating that a term is random.
))
)) Components ss and df are just those computed from a MacAnova anova()
)) command (possibly with marg:T as needed), and may not be in conformance
)) with the model as used by ems for the following reasons:
)) 1. anova() computes only hierarchical models, while you may specify
))    nonhierarchical models in ems by using nonhier:T.
)) 2. ems enforces nesting.  If b first appears in a.b then b is nested
))    in a and any appearance of b in a later term implies the presence
))    of a.  anova() does no such enforcing.  For example, in
))    "y=a+a.b+c+b.c", b.c would be interpreted by ems as a.b.c while
))    anova() would not include a.b.c in the model.
))
)) ems() uses the "synthesis" method of Hartley, as explained in 10.5.2 of
)) Hocking's linear models book.
# usage: $S(Model,Randomvars [,marg:T, restrict:T or F, nonhier:T, keep:T])
if(!ismacro(colproduct)) {
	getmacros(colproduct,quiet:T,printname:F)
}
if(!ismacro(makemat)) {
	getmacros(makemat,quiet:T,printname:F)
}
if(!ismacro(buildRFStr)) {
	getmacros(buildRFStr,quiet:T,printname:F)
}

@model <- argvalue($1,"argument 1","string")
if (match("*=*",@model,0,exact:F) == 0){
	paste("\"",@model,"\" is not a model of the form \"lhs=rhs\"")
}

@keep <- keyvalue($K,"keep","TF",default:F)

@print <- keyvalue($K,"print","TF", default:!@keep)

@restrict <- keyvalue($K,"restrict*","TF",default:T)
@unrest <- !@restrict

@marg <- keyvalue($K,"marg*","TF", default:F)

@nonhier <- keyvalue($K,"nonhier*","TF", default:F)

if (@nonhier){
	if(@unrest) {
		error("don't use 'nonhier:T' with 'restrict:F'")
	}
	if(@marg) {
		error("don't use 'nonhier:T' with 'marg:T' in $S()", macroname:F)
	}
}

@canpush <- alltrue(isfunction(pushmodel),pushmodel(canpush:T))
if (@canpush){
	pushmodel()
}
 # collect information about the model in use
anova(@model,silent:T,unbal:T,marg:@marg)
@df <- DF
@bitmodel <- modelinfo(bitmodel:T)
@nterms <- length(@bitmodel)
@tnames <- modelinfo(termnames:T)
@colcount <- modelinfo(colcount:T)
@offsets <- autoreg(1,vector(0,@colcount))
@tmpopt <- getoptions(warning:T)
setoptions(warning:F)
@xmat <- xvariables(missing:?)
setoptions(warning:@tmpopt)
@userows <- !ismissing(@xmat[,1])
@vnames <- varnames()[-1]
@yname <- varnames()[1]
@nvars <- length(@vnames)

if(@marg) {
	@xtxinv <- modelinfo(xtxinv:T)
	@aliased <- modelinfo(aliased:T)
}

@random <- rep(F,@nvars)
@rvars <- argvalue($2, "argument 2")
if(!isnull(@rvars)) {
	if(isreal(@rvars)) {
		argvalue(@rvars, "argument 2", "positive integer vector")
		if (max(@rvars) > @nvars){
			error("element of argument 2 > number of variables in model")
		}
	} elseif( ischar(@rvars)) {
		@tmp <- match(@rvars,@vnames,0)
		if(min(@tmp) == 0) {
			@out <- paste("variable(s)",@rvars[@tmp == 0],\
				"named in $2\nare not in model")
			error(@out)
		}
		@rvars <- @tmp
	} else {
		error("$2 is not REAL, CHARACTER, or NULL.")
	}
	@random[@rvars] <- T
}

@nlevels <- rep(0,@nvars)
@variate <- rep(F,@nvars)
 # check for variates
for(@i,run(@nvars)) {
	@nlevels[@i] <- max(<<@vnames[@i]>>)
	if(!isfactor(<<@vnames[@i]>>)) {
		@variate[@i] <- T
		@nlevels[@i] <- 1
	#	error("model $1 contains 1 or more variates")
	}
	if(@random[@i] && @variate[@i]) {
		error("variates must be fixed effect")
	}
}

 # build up matrices to indicate levels of factors in the
 # original data.  These are used to construct random and
 # mixed effects.

 # base matrices for random effects are just indicators for those
 # levels.  Base matrices for fixed effects are columns from the
 # design matrix.  We will also need to deal with nesting and
 # grouping.

 # there is a random-base matrix and a fixed-base matrix for
 # every variable

 # we also need to enforce nesting and grouping

@tmp <- buildRFStr(@nvars,@vnames,@nlevels,@userows,@bitmodel,@random,\
  @offsets,@colcount,@xmat,@nterms,@rvars,@variate)
@Rstr <- @tmp$Rstr
@Fstr <- @tmp$Fstr
@nesting <- @tmp$nesting
@grouping <- @tmp$grouping
@bitmodel <- @tmp$bitmodel
@random <- @tmp$random
@rvars <- @tmp$rvars

# which terms are random

if(!isnull(@rvars)) {
	@rterms <- nbits(sum(2^(@rvars-1)) %& @bitmodel) > 0
} else {
	@rterms <- rep(F,@nterms)
}

@mmat <- NULL
@mccount <- NULL
@info <- structure(random:@random,bitmodel:@bitmodel,Rstr:@Rstr,\
	Fstr:@Fstr,nesting:@nesting,grouping:@grouping,variate:@variate,\
	unrest:@unrest,nonhier:@nonhier,nlevels:@nlevels)

for(@thisterm,run(@nterms)) {
	@thisout <- makemat(@info,@thisterm)
	@mmat <- hconcat(@mmat,@thisout$matrix)
	@mccount <- vector(@mccount,ncols(@thisout$matrix))
}

 # @coefs[i,j] has coefficient for effect j in model term i
@coefs <- matrix(rep(0,(@nterms+1)^2),@nterms+1)
@coefs[,@nterms+1] <- 1

@pivoted <- 0
@tol <- 1e-5
@getd <- 1*(run(@nterms) == rep(run(@nterms),@mccount)')

if(@marg) {
	# for marginal (type III) SS

	# first remove pivoted columns
	@xmat <- @xmat[@userows,!@aliased]
	@xtxinv <- @xtxinv[!@aliased,!@aliased]
	@tmp <- 1*(run(@nterms) == rep(run(@nterms),@colcount)')
	@colcount <-- @tmp %*% (1*@aliased)
	@mcols <- sum(@mccount)

	for(@thisterm,run(@nterms)) {
		if(@df[@thisterm] > 0) {
			@cols <- run(@pivoted+1,@pivoted+@colcount[@thisterm])
			@pivoted <-+ @colcount[@thisterm]
			@tmpxtx <- swp(@xtxinv,@cols)
			@tmpx <- @xmat[,@cols] -\
				@xmat[,-@cols] %*% (-@tmpxtx[-@cols,@cols])
			@tmpx <- hconcat(@mmat,@tmpx)
			@tmpxtx <- @tmpx %c% @tmpx
			@thisd <- @getd %*% diag(@tmpxtx)[run(@mcols)]
			@tmpxtx <- swp(@tmpxtx,\
				run(@mcols+1,@mcols+@df[@thisterm]))
			@newd <- @getd %*% diag(@tmpxtx)[run(@mcols)]
			@diffd <- @thisd - @newd
			@coefs[@thisterm,run(@nterms)] <- @diffd/@df[@thisterm]
		}
	}
} else {
	# for type I SS
	@xtx <- @mmat %c% @mmat

	@thisd <- @getd %*% diag(@xtx)

	for(@thisterm,run(@nterms)) {
		@thisdf <- 0
		for(@thiscol,run(@pivoted+1,@pivoted+@mccount[@thisterm])) {
			if(@xtx[@thiscol,@thiscol] > @tol) {
				@xtx <- swp(@xtx,@thiscol)
				@thisdf <-+ 1
			}
		}
		@df[@thisterm] <- @thisdf
		@pivoted <-+ @mccount[@thisterm]
		@newd <- @getd %*% diag(@xtx)
		@diffd <- @thisd - @newd
		@getd[@thisterm,] <- 0
		if(@df[@thisterm] > 0) {
			for(@j,run(@nterms,@thisterm)) {
					@cf <- if(@j == @thisterm){
						@thisd[@thisterm]
					} else {
						@diffd[@j]
					}
					@cf <- @cf/@df[@thisterm]
					@coefs[@thisterm,@j] <- @cf
			}
		}
		@thisd <- @newd
	}
}
if (@keep){
	@ss <- SS
}
if (@canpush){
	popmodel()
}


if(@print) {
	for(@i,run(@nterms+1)) {
		@out <- paste("EMS(",@tnames[@i],") = ",sep:"")
		if(@df[@i] == 0) {
			@out <- paste(@out,"cannot be estimated")
		} else {
			@out <- paste(@out,"V(",@tnames[@nterms+1],")",sep:"")
			for(@j,run(@nterms,1)) {
				@symbol <- if(@rterms[@j]) {
					"V("
				} else {
					"Q("
				}
				if(@coefs[@i,@j] > @tol) {
					@cf <- @coefs[@i,@j]
					@out <- paste(@out," + ",@cf,\
						@symbol,@tnames[@j],")",sep:"")
				}
			}
		}
		print(@out);;
	}
}
if(delete(@keep,return:T)) {
	structure(df:delete(@df,return:T),ss:delete(@ss,return:T),\
		termnames:delete(@tnames,return:T),coefs:delete(@coefs,return:T),\
		rterms:delete(@rterms,return:T))
} else {
	delete(@df,@tnames,@coefs,@rterms)
	NULL
}
%ems%

====> reml <====
reml     MACRO  DOLLARS
) reml(Model,random:Randomvars,Z:Zmatrices) performs a restricted maximum
) likelihood analysis for the model given in CHARACTER scalar Model.  At
) least one of random: or Z: must be included. Randomvars is a CHARACTER
) vector specifying the names of factors in the model which are random.
) Zmatrices is a CHARACTER vector specifying the names of REAL matrices with
) no MISSING values.  All matrices in Zmatrices should have the same number
) of rows as the response in Model.  Each column of every matrix in Zmatrices
) will be given an independent random effect.  Variances of the random
) effects can differ between matrices.
) 
) The return value of reml() is a structure with the following components:
)      theta:   estimates of the fixed effects
)        phi:   estimates of the variance components
)   thetavar:   variance matrix of the fixed effects
)     phivar:   variance matrix of the variance components
)      phidf:   equivalent degrees of freedom for the variance components
)      gamma:   estimates (predictions) of random effects
)   gammavar:   variances of predictions of random effects
)          L:   REML log likelihood
)
) Any variates in the model must be fixed effects.
)
)) reml() assumes that if a factor first appears in an interaction, then
)) that factor is nested in the other terms of the interaction.  For
)) example, if the first appearance of factor c is in the term a.b.c, then
)) c is assumed nested in the a.b combinations.  This nesting is assumed in
)) the remainder of the model.  That is, continuing the example, if there
)) is a later term c.d, it will be interpreted as a.b.c.d even though
)) a.b.c.d is not specifically in the model.
))
)) reml() works for both balanced and unbalanced data.
))
) reml(Model,random:Randomvars,restrict:F) performs the REML analysis
) assuming no marginal restrictions on the random effects in the model.
)
) reml(Model,random:Randomvars,nonhier:T) performs the REML analysis for an
) analysis of variance that does not enforce the usual MacAnova hierarchy
) assumptions.
)
)) That is, for example, model "y=a+b+c+a.b.c" does not imply
)) that the two-way interaction degrees of freedom are part of the "a.b.c"
)) term.  You cannot use anova() to compute such an analysis although it
)) can be done (if you know how) using swp().
))
) reml(Model,random:Randomvars,usemle:T) performs a maximum likelihood
) analysis instead of the REML analysis.
)
) reml(Model,random:Randomvars,retV:T) returns the estimated covariance
) matrix of the data as an additional component named V.
)
) reml(Model,random:Randomvars,tolerance:value) uses value as a tolerance for
) determining singularity and convergence (default is 1e-10).
)
) reml(Model,random:Randomvars,maxiter:value) uses value as the maximum
) number of iterations in the fitting process (default is 60).
)
))040827 Changed to specify random terms by either random:termlist
))       or Z:termlist; it is also faster and returns predictions of
))       random effects
))040902 Added a bunch of delete()s
# usage: $S(Model,random:Randomvars,Z:Zmatrices [,restrict:T or F, nonhier:T, ])
if(!ismacro(colproduct)) {
	getmacros(colproduct,quiet:T,printname:F)
}
if(!ismacro(makemat)) {
	getmacros(makemat,quiet:T,printname:F)
}
if(!ismacro(buildRFStr)) {
	getmacros(buildRFStr,quiet:T,printname:F)
}

@model <- argvalue($1,"argument 1","string")
if (match("*=*",@model,0,exact:F) == 0){
	paste("\"",@model,"\" is not a model of the form \"lhs=rhs\"")
}

@restrict <- keyvalue($K,"restrict*","TF",default:T)
@unrest <- !@restrict

@nonhier <- keyvalue($K,"nonhier*","TF", default:F)

@usemle <- keyvalue($K,"usemle","TF", default:F)

@tolerance <- keyvalue($K,"tol*","positive number",default:1e-10)

@kmax <- keyvalue($K,"maxit*","positive count",default:60)

@retV <- keyvalue($K,"retV","TF",default:F)

@rnames <- keyvalue($K,"random","character vector")

@Znames <- keyvalue($K,"Z","character vector")

if(isnull(@rnames) && isnull(@Znames)) {
	error("you must specify at least one of random: or Z:")
}

if (@nonhier){
	if(@unrest) {
		error("don't use 'nonhier:T' with 'restrict:F'")
	}
}

@canpush <- alltrue(isfunction(pushmodel), pushmodel(canpush:T))
if (@canpush) {
	pushmodel()
}
 # collect information about the model in use

anova(@model,silent:T,unbal:T)

if(TERMNAMES[1] != "CONSTANT") {
	error("model must include the constant")
}

@df <- DF
@bitmodel <- modelinfo(bitmodel:T)
@nterms <- length(@bitmodel)
@tnames <-	TERMNAMES
@colcount <- modelinfo(colcount:T)
@offsets <- autoreg(1,vector(0,@colcount))
@tmpopt <- getoptions(warning:T);setoptions(warning:F)
@xmat <- xvariables(missing:?)
setoptions(warning:@tmpopt)
@userows <- !ismissing(@xmat[,1])
@vnames <- varnames()[-1]
@yname <- varnames()[1]
@nvars <- length(@vnames)
@N <- sum(@userows)
@Y <- modelinfo(y:T)[@userows]

if (@canpush){
	popmodel()
}
@random <- rep(F,@nvars)
@rvars <- NULL
if(!isnull(@rnames)) {
	@tmp <- match(@rnames,@vnames,0)
	if(min(@tmp) == 0) {
		@out <- paste("variable(s)",@rnames[@tmp == 0],\
			"named in random:2\nare not in model")
		error(@out)
	}
	@rvars <- @tmp
	@random[@rvars] <- T
}

@tmpn <- length(@userows)
if(!isnull(@Znames)) {
	for(@i,run(length(@Znames))) {
		if( !alltrue( isdefined(<<@Znames[@i]>>), isreal(<<@Znames[@i]>>), \
			ismatrix(<<@Znames[@i]>>), !anymissing(<<@Znames[@i]>>),\
			nrows(<<@Znames[@i]>>) == @tmpn) ) {
			@out <- paste("variable",@Znames[@i],\
				"named in Z: is not a real matrix with no",\
				"missing values and number of rows equal to",\
				"length of",@yname)
			error(@out)
		}
	}
}

@nlevels <- rep(0,@nvars)
@variate <- rep(F,@nvars)
 # check for variates
for(@i,run(@nvars)) {
	@nlevels[@i] <- max(<<@vnames[@i]>>)
	if(!isfactor(<<@vnames[@i]>>)) {
		@variate[@i] <- T
		@nlevels[@i] <- 1
	}
	if(@random[@i] && @variate[@i]) {
		error("variates must be fixed effect")
	}
}

 # build up matrices to indicate levels of factors in the
 # original data.  These are used to construct random and
 # mixed effects.

 # base matrices for random effects are just indicators for those
 # levels.  Base matrices for fixed effects are columns from the
 # design matrix.  We will also need to deal with nesting and
 # grouping.

 # there is a random-base matrix and a fixed-base matrix for
 # every variable

 # we also need to enforce nesting and grouping

@tmp <- buildRFStr(@nvars,@vnames,@nlevels,@userows,@bitmodel,@random,\
  @offsets,@colcount,@xmat,@nterms,@rvars,@variate,scalefix:F)
@Rstr <- @tmp$Rstr
@Fstr <- @tmp$Fstr
@nesting <- @tmp$nesting
@grouping <- @tmp$grouping
@bitmodel <- @tmp$bitmodel
@random <- @tmp$random
@rvars <- @tmp$rvars

# which terms are random

if(!isnull(@rvars)) {
	@rterms <- nbits(sum(2^(@rvars-1)) %& @bitmodel) > 0
} else {
	@rterms <- rep(F,length(@bitmodel))
}
@Ulabs <- vector(@tnames[vector(@rterms,F)],@Znames,@tnames[length(@tnames)])
@nphisp1 <- length(@Ulabs)
@nphis <- @nphisp1 - 1

@phis <- rep(1,@nphis)
@phi0 <- 1
@U <- NULL
@ms <- rep(0,@nphis)

@X <- NULL
@Xlabs <- NULL

@tmpr <- 1
#@termfit <- rep(F,2^@nvars)
@info <- structure(random:@random,bitmodel:@bitmodel,Rstr:@Rstr,\
	Fstr:@Fstr,nesting:@nesting,grouping:@grouping,variate:@variate,\
	unrest:@unrest,nonhier:@nonhier,nlevels:@nlevels)

for(@thisterm,run(@nterms)) {
	@thisout <- makemat(@info,@thisterm)
	if(@rterms[@thisterm]) {
		@ms[@tmpr] <- ncols(@thisout$matrix)
		@U <- hconcat(@U,@thisout$matrix)
		@tmpr <-+ 1
	} else {
		@X <- hconcat(@X,@thisout$matrix)
		if(@thisterm == 1) {
			@tmp2 <- "CONSTANT"
		} else {
			@tmp <- ncols(@thisout$matrix)
			@tmp2 <- NULL
			for(@i,run(@tmp)) {
				@tmp2 <- vector(@tmp2,paste(@tnames[@thisterm],@i))
			}
		}
		@Xlabs <- vector(@Xlabs,@tmp2)
	}
}

if(!isnull(@Znames)) {
	for(@i,run(length(@Znames))) {
		@U <- hconcat(@U,<<@Znames[@i]>>[@userows,])
		@ms[@nphis-length(@Znames)+@i] <- ncols(<<@Znames[@i]>>)
	}
}

# remove singular columns of X
@R <- qr(@X,ronly:T)
@tmp <- diag(@R)^2/vector(sum(@R^2))
@X <- @X[,@tmp > @tolerance]
@Xlabs <- @Xlabs[@tmp > @tolerance]
@r <- ncols(@X)

@Usubs <- split(run(@nphis)')
@tmp <- autoreg(1,vector(0,@ms))
for(@i,run(2,length(@tmp))) {
	@Usubs[@i-1] <- run(@tmp[@i-1]+1,@tmp[@i])
}

@XtX <- @X %c% @X
@XtY <- @X %c% @Y
@UtX <- @U %c% @X
@UtU <- @U %c% @U
@YU <- @Y %c% @U

if(@usemle) {
	@XU <- @X %c% @U
	@YY <- sum(@Y^2)
} else {
	# this is the old way.  It involves computing the projection
	# matrix, and we'd like to avoid that if possible.
	#@tmpP <- @X %*% solve(@X %c% @X, @X')
	#@tmpe <- svd(@tmpP,right:T)
	#@tmpvals <- @tmpe$values
	#@tmpuse <- @tmpvals/@tmpvals[1] < @tolerance
	#@P11 <- @tmpe$rightvectors[,@tmpuse]'

	#@PU <- @P11 %*% @U
	#@UPPU <- @PU %c% @PU
	#@YPPU <- @Y %c% (@P11 %c% @PU)
	#@YPPY <- @Y %c% (@P11 %c% @P11) %*% @Y
	#@XPPU <- @X %c% (@P11 %c% @PU)

	# old way involved @P11, a basis for the null space
	# new way finds @P22, as basis for the nonnull space
	# and uses @P11' @P11 = I - @P22 @P22'

	@tmpe <- svd(@X,left:T,maxit:250)
	@tmpvals <- @tmpe$values
	@tmpuse <- @tmpvals/@tmpvals[1] > @tolerance
	@P22 <- @tmpe$leftvectors[,@tmpuse]
	# @P11' @P11 = I - @P22 @P22'
	@PPU <- @U -  @P22 %*% (@P22' %*% @U)
	@UPPU <- @U %c% @PPU
	@YPPU <- @Y %c% @PPU
	@XPPU <- @X %c% @PPU
	@YPPY <- @Y %c% (@Y - @P22 %*% (@P22 %c% @Y))
	delete(@PPU,@P22,@tmpuse,@tmpvals)
}


@Omega <- matrix(rep(0,@nphisp1^2),@nphisp1)
@rhos <- rep(0,@nphisp1)

# use Hemmerle and Hartley a la Hocking

for(@k,run(@kmax)) {

	@inQ <- vector(@Usubs[@phis > 0])

	# first get theta
	#we need to be careful if all varcomps are estimated 0, and
	# thus inQ is NULL
	if(isnull(@inQ)) {
		@XVXinv <- solve(@XtX)
		@theta <- @XVXinv %*% @XtY 
	} else {
		@Qinv <- solve(dmat(rep(@phi0/@phis[@phis>0],@ms[@phis>0])) +\
			@UtU[@inQ,@inQ])

		@XVXinv <- solve(@XtX - @UtX[@inQ,] %c% @Qinv %*% @UtX[@inQ,])
		@theta <- @XVXinv %*% (@XtY - @UtX[@inQ,] %c% @Qinv %*% @YU[,@inQ]')
	}

	@tres <- @Y - @X %*% @theta
	if(@usemle) {
		if(isnull(@inQ)) {
			@Bu <- (@UtU)/@phi0
		} else {
			@Bu <- (@UtU - @UtU[,@inQ] %*% @Qinv %*% @UtU[@inQ,])/@phi0
		}
	} else {
		if(isnull(@inQ)) {
			@Bu <- (@UPPU)/@phi0
		} else {
			@Qinv <- solve(dmat(rep(@phi0/@phis[@phis > 0],@ms[@phis > 0])) + \
				@UPPU[@inQ,@inQ])
			@Bu <- (@UPPU - @UPPU[@inQ,] %c% @Qinv %*% @UPPU[@inQ,])/@phi0
		}
	}

	for(@i,run(@nphis)) {
		for(@j,run(@i,@nphis)) {
	#		@Omega[@i,@j] <- trace(@Bu[@Usubs[@i],@Usubs[@j]] %*% @Bu[@Usubs[@j],@Usubs[@i]])

			@Omega[@j,@i] <- @Omega[@i,@j] <- \
				sum(vector(@Bu[@Usubs[@j],@Usubs[@i]])^2)
		}
	}


	if(@usemle) {
		@Omega[@nphisp1,@nphisp1] <- @N/@phi0^2
		if(!isnull(@inQ)) {
			@tmp <- @Qinv %*% @UtU[@inQ,@inQ]
		#	@Omega[@nphisp1,@nphisp1] <- (@N - 2*trace(@tmp) + trace(@tmp %*% @tmp))/@phi0^2
			@Omega[@nphisp1,@nphisp1] <- (@N - 2*trace(@tmp) + \
				sum(vector(@tmp	* @tmp')))/@phi0^2
		}
		for(@i,run(@nphis)) {
			@tmpO <- trace(@UtU[@Usubs[@i],@Usubs[@i]])
			if(!isnull(@inQ)) {
				@tmp <- @UtU[@Usubs[@i],@inQ]
				@tmpO <-- 2*sum(vector((@tmp %*% @Qinv) * @tmp))
		#		@tmpO <-- 2*trace(@tmp %*% @Qinv %*% @tmp')
				@tmpO <-+ sum(vector((@tmp %*% @Qinv %*% @UtU[@inQ,@inQ] %*% @Qinv) * @tmp))
		#		@tmpO <-+ trace(@tmp %*% @Qinv %*% @UtU[@inQ,@inQ] %*% @Qinv %*% @tmp')
			}
			@Omega[@nphisp1,@i] <- @tmpO/@phi0^2
			@Omega[@i,@nphisp1] <- @tmpO/@phi0^2
		}
		if(isnull(@inQ)) {
			@BYU <- (@YU)/@phi0
			@BXU <- (@XU)/@phi0
		} else {
			@tmp <- @Qinv %c% @UtU[@inQ,]
			@BYU <- (@YU - @YU[,@inQ] %*% @tmp)/@phi0
			@BXU <- (@XU - @XU[,@inQ] %*% @tmp)/@phi0
		}
		@tmp <- vector(@BYU - @theta %c% @BXU)
		for(@i,run(@nphis)) {
			@rhos[@i] <- sum(@tmp[@Usubs[@i]]^2)
		}
		@R <- @Y - @X %*% @theta
		if(isnull(@inQ)) {
			@tmp <- (@R)/@phi0
		} else {
			@tmp <- (@R - @U[,@inQ] %*% (@Qinv %c% (@U[,@inQ] %c% @R)))/@phi0
		}
		@rhos[@nphisp1] <- sum(@tmp^2)
	} else {
		@Omega[@nphisp1,@nphisp1] <- (@N - @r)/@phi0^2
		if(!isnull(@inQ)) {
			@tmp <- @Qinv %c% @UPPU[@inQ,@inQ]
		#		@Omega[@nphisp1,@nphisp1] <- (@N - @r - 2*trace(@tmp) + trace(@tmp %*% @tmp))/@phi0^2
			@Omega[@nphisp1,@nphisp1] <- (@N - @r - 2*trace(@tmp) + \
									  sum(vector(@tmp * @tmp')))/@phi0^2
		}
		for(@i,run(@nphis)) {
			@tmpO <- trace(@UPPU[@Usubs[@i],@Usubs[@i]])
			if(!isnull(@inQ)) {
				@tmp <- @UPPU[@Usubs[@i],@inQ]
				@tmpO <-- 2*sum(vector(@tmp' * (@Qinv %c% @tmp')))
		#		@tmpO <-+ trace(@tmp %*% @Qinv %*% @UPPU[@inQ,@inQ] %*% @Qinv %*% @tmp')
				@tmpO <-+ sum(vector((@tmp %*% (@Qinv %c% @UPPU[@inQ,@inQ] %*% @Qinv)) * @tmp))
			}
			@Omega[@nphisp1,@i] <- @tmpO/@phi0^2
			@Omega[@i,@nphisp1] <- @tmpO/@phi0^2
		}
		if(isnull(@inQ)) {
			@tmp <- vector(@YPPU)/@phi0
		} else {
			@tmp <- vector(@YPPU - @YPPU[,@inQ] %*% (@Qinv %c% @UPPU[@inQ,]))/@phi0
		}
		for(@i,run(@nphis)) {
			@rhos[@i] <- sum(@tmp[@Usubs[@i]]^2)
		}
		if(isnull(@inQ)) {
			@rhos[@nphisp1] <- (@YPPY)/@phi0^2
		} else {
			@rhos[@nphisp1] <- (@YPPY - 2*@YPPU[,@inQ] %*% (@Qinv %c% @YPPU[,@inQ]') + \
		   		@YPPU[,@inQ] %*% (@Qinv %c% @UPPU[@inQ,@inQ]) %*% (@Qinv %c% @YPPU[,@inQ]'))/@phi0^2
		}
	}

	@oldphis <- vector(@phis,@phi0)

	@phis <- solve(@Omega,@rhos)

	# only allow a gradual zeroing; give small ones a few
	# iterations before letting them zero completely
	if(@k <= 5) {
		@zeroed <- @phis < .01^@k
		while(min(@phis) < .01^@k) {
			@phis[!@zeroed] <- solve(@Omega[!@zeroed,!@zeroed],\
				(@rhos-@Omega[,@zeroed] %*% rep(.01^@k,sum(@zeroed)))[!@zeroed])
			@phis[@zeroed] <- .01^@k
			@zeroed$$ <- @zeroed$$ || @phis$$ < .01^@k
		}
	} else {
		@zeroed <- @phis < 0
		while(min(@phis) < 0) {
			@phis[!@zeroed] <- solve(@Omega[!@zeroed,!@zeroed],@rhos[!@zeroed])
			@phis[@zeroed] <- 0
			@zeroed$$ <- @zeroed$$ || @phis$$ < 0
		}
	}

	if(sum(abs(@oldphis - @phis))/sum(@oldphis+1e-30) < @tolerance) {
		break
	}

	@phi0 <- @phis[@nphisp1]
	@phis <- @phis[-@nphisp1]
}
delete(@UtU, @UtX, @XtX, @XtY, @YU)

if(@k == @kmax) {
	print("WARNING: reml may not have converged")
	@phis <- vector(@phis,@phi0)
}

@R <- @Y - @X %*% @theta

@Usig <- rep(@phis[-@nphisp1],@ms)
@W <- @U * sqrt(@Usig')
@UG <- @U * @Usig'
@sig <- @phis[@nphisp1]

# V = sig*I + WW' so eigenvalues of V are
# sig+diag(d^2) and sig, where d are svs of W
@eigenV <- rep(@sig,nrows(@W))
@use <- run(min(dim(@W)))
@eigenV[@use] <- @eigenV[@use] + svd(@W)[@use]^2

# V = sig*I + WW' so
# V^-1 = 1/sig*I - W(I+W'W/sig)^-1 W'/sig^2

@RtWsig <- @R %c% @W / @sig
@RVinvR <- sum(@R^2)/@sig - @RtWsig %*% \
	solve(dmat(rep(1,ncols(@W)))+@W %c% @W/@sig,@RtWsig')


@V <- @W %*% @W'
@V <-+ dmat(rep(@sig,@N))
 # these likelihoods usually match SAS PROC MIXED
# V = sig*I + WW' so
# V^-1 = 1/sig*I - W(I+W'W/sig)^-1 W'/sig^2

@XtWsig <- @X %c% @W / @sig
@XVinvX <- @X %c% @X/@sig - @XtWsig %*% \
	solve(dmat(rep(1,ncols(@W)))+@W %c% @W/@sig,@XtWsig')
delete(@XtWsig)
if(@usemle) {
	@L <- -.5*@N*log(2*PI) - .5*sum(log(@eigenV))
	@L <-- .5 * @RVinvR
} else {
#	@tmp <- solve(@V,@X)
#	@XVinvX <- @X %c% @tmp
#	@L <- -.5*sum(log(eigenvals(@V))) -.5*sum(log(eigenvals(@XVinvX))) \
#		-.5*@R' %*% solve(@V,@R) - (@N-@r)/2*log(2*PI)

	@L <- -.5*sum(log(@eigenV)) -.5*sum(log(eigenvals(@XVinvX))) \
		-.5*@RVinvR - (@N-@r)/2*log(2*PI)
	delete(@RVinvR)
}

@VinvUG <- @UG/@sig - @W %*% solve(dmat(rep(1,ncols(@W)))+@W %c% @W/@sig,\
	@W %c% @UG/@sig^2)
@tmp <- @X %*% solve(@XVinvX,@X %c% @VinvUG)
delete(@XVinvX)

@Vinvtmp <- @tmp/@sig - @W %*% solve(dmat(rep(1,ncols(@W)))+@W %c% @W/@sig,\
	@W %c% @tmp/@sig^2)
@lambda <- @VinvUG - @Vinvtmp
@gamma <- @lambda %c% @Y
@tmp <- @lambda %c% @W
@gammavar <- @lambda %c% @lambda*@sig + @tmp %*% @tmp'

@gammavar <- diag(@gammavar) + @Usig-2*vector(sum(@lambda*(@U*@Usig')))
delete(@Usig, @U, @tmp, @VinvUG, @Vinvtmp)

setlabels(@theta,structure(@Xlabs,""))
setlabels(@phis,structure(@Ulabs,""))
setlabels(@XVXinv,structure(@Xlabs,@Xlabs))
@Omegainv <- @Omega
@Omegainv[!@zeroed,!@zeroed] <- solve(@Omega[!@zeroed,!@zeroed])*2
@Omegainv[@zeroed,] <- 0
@Omegainv[,@zeroed] <- 0
setlabels(@Omegainv,structure(@Ulabs,@Ulabs))
@phidf <- 2*@phis^2/diag(@Omegainv)
@XVXinv <- @XVXinv * @phis[length(@phis)]
delete(@Xlabs)
if(delete(@retV,return:T)) {
        structure(theta:delete(@theta,return:T),phi:delete(@phis,return:T),\
		thetavar:delete(@XVXinv,return:T),phivar:delete(@Omegainv,return:T),\
		phidf:delete(@phidf,return:T), gamma:delete(@gamma,return:T),\
		gammavar:delete(@gammavar,return:T),loglike:delete(@L,return:T),\
		residuals:delete(@R,return:T),V:delete(@V,return:T))
} else {
        structure(theta:delete(@theta,return:T),phi:delete(@phis,return:T),\
		thetavar:delete(@XVXinv,return:T),phivar:delete(@Omegainv,return:T),\
		phidf:delete(@phidf,return:T), gamma:delete(@gamma,return:T),\
		gammavar:delete(@gammavar,return:T),\
		loglike:delete(@L,return:T),residuals:delete(@R,return:T))
}
%reml%

====> colproduct <====
colproduct      MACRO  DOLLARS
) this macro is used in ems and makemat
# usage: $S(n1,n2)
@n1 <- ncols($1)
@n2 <- ncols($2)
@tmpopt <- getoptions(warning:T);setoptions(warning:F)
@out <- $1[,rep(run(@n1),@n2)]*$2[,rep(run(@n2),rep(@n1,@n2))]
setoptions(warning:delete(@tmpopt,return:T))
delete(@n1,@n2)
delete(@out,return:T)
%colproduct%


buildRFStr     MACRO  DOLLARS

# build up matrices to indicate levels of factors in the
# original data.  These are used to construct random and
# mixed effects.

# base matrices for random effects are just indicators for those
# levels.  Base matrices for fixed effects are columns from the
# design matrix.  We will also need to deal with nesting and
# grouping.

# there is a random-base matrix and a fixed-base matrix for
# every variable

# we also need to enforce nesting and grouping

@nvars <- $1
@vnames <- $2
@nlevels <- $3
@userows <- $4
@bitmodel <- $5
@random <- $6
@offsets <- $7
@colcount <- $8
@xmat <- $9
@nterms <- $10
@rvars <- $11
@variate <- $12
@scalefix <- keyvalue($K,"scalefix","TF",default:T)


@tmpopt <- getoptions(warning:T);setoptions(warning:F)
@Rstr <- structure(NULL)
@Fstr <- structure(NULL)
setoptions(warning:@tmpopt)

 # @nesting[i,j] == T if variable i is nested in variable j
 # @grouping[i,j] == T if variable i is in a group led by variable j
@nesting <- matrix(rep(F,@nvars*@nvars),@nvars)
@grouping <- matrix(rep(F,@nvars*@nvars),@nvars)


for(@i,run(@nvars)) {
	if(@variate[@i]) {
		@Rstr <- changestr(@Rstr,@i,(<<@vnames[@i]>>)[@userows,])
	} else {
		@Rstr <- changestr(@Rstr,@i,1*(<<@vnames[@i]>> == run(@nlevels[@i])')[@userows,])
	}
}

@i <- 1
while(@i <= @nvars) {
	# we are going to have to enforce nesting later when it occurs
	# earlier, so we have to look for nesting (a variable ocurring
	# first in an interaction).  We also have to fix up bitmodels
	# so that we don't come up with too many columns later.

	@whichtrms <- run(@nterms)[(@bitmodel %& 2^(@i-1)) > 0]

	# find any nesting and/or grouping of terms

	# variables in the first term where var @i appears
	@whichvars <- (@bitmodel[@whichtrms[1]] %& 2^run(0,@nvars-1)) > 0

	# last variable in the first term where var @i appears
	@ilast <- max(run(@i,@nvars)[@whichvars[run(@i,@nvars)]])

	# fix up Rstr to deal with any grouping
	@thismat <- rep(1,sum(@userows))
	for(@j,run(@i,@ilast)) {
		if(@whichvars[@j]) {
			@thismat <- colproduct(@thismat,@Rstr[@j])
		}
	}
	@Rstr <- changestr(@Rstr,@i,@thismat)

	# treat variables from @i to @ilast as one group fused together
	# if any are random, treat whole bunch as random
	if(sum(@random[run(@i,@ilast)]) > 0) {
		# here is random
		@random[run(@i,@ilast)] <- T # enforce that @i will be random
		@rvars <- run(@nvars)[@random]
	}

	# now do Fmat, just a column of 1's unless all variables fixed
	@thismat <- rep(1,sum(@userows))
	if(sum(@random[run(@i,@ilast)]) == 0) {
		# here is fixed, just get cols from original anova
		@thismat <- @xmat[,@offsets[@whichtrms[1]]+run(@colcount[@whichtrms[1]])]
		# rescale so that multiplier for Q term matches usual definitions
		@thismat <- @thismat[@userows,]
		if(@scalefix) {
			@thismat <- @thismat/sqrt(2^nbits(@bitmodel[@whichtrms[1]]))
			if(@ilast>@i) {
				# more rescaling for fused variables
				@tmp1 <- run(ncols(@thismat))
				@tmpdiv <- prod(@nlevels[run(@i+1,@ilast)])
				for(@j,run(@i,@ilast)) {
					@thisind <- floor(@tmp1/@tmpdiv) %% @nlevels[@j]
					@thisind <- @thisind == 0
					# if variable j is missing from this column, scale by 2 (for lost var)
					# over nlevels of j
					@thismult <- sqrt(2/@nlevels[@j])
					@thismat[,@thisind] <- @thismat[,@thisind]*@thismult
					if(@j<@nvars) {
						@tmpdiv <- @tmpdiv/@nlevels[@j+1]
					}
				}
			}
		}
	}
	@Fstr <- changestr(@Fstr,@i,@thismat)

	# for later vars in group, just give a column of ones
	if(@i < @ilast) {
		for(@j,run(@i+1,@ilast)) {
			@Fstr <- changestr(@Fstr,@j,rep(1,sum(@userows))'')
			@Rstr <- changestr(@Rstr,@j,rep(1,sum(@userows))'')
		}
	}

	# now fix nesting/grouping indicators and bit model to enforce
	# nesting and grouping

	# vars we nest into (may be empty)
	@nestbits <- @bitmodel[@whichtrms[1]] %& (2^(@i-1) - 1)
	if(@i > 1) {
		@nesting[@i,run(@i-1)] <- (@nestbits %& 2^run(0,@i-2)) > 0
	}

	# vars in this group (may be just this var)
	@groupbits <- 2^@ilast - 2^(@i-1)
	if(@ilast > @i) {
		for(@j,run(@i+1,@ilast)) {
			@grouping[@j,@i] <- T
		}
	}

	# now fixup bitmodel so that nesting is enforced
	# this means to zero out the nested-in bits, and make sure bit @i
	# is on if any in the group is on, and take out higher bits in
	# grouping bits
	for(@j,run(@nterms)) {
		if((@bitmodel[@j] %& @groupbits) > 0) {
			@bitmodel[@j] <- @bitmodel[@j] %|  @nestbits
			@bitmodel[@j] <- @bitmodel[@j] %& (%! @groupbits)
			@bitmodel[@j] <- @bitmodel[@j] %| (2^(@i-1))
		}
	}

	@i <- @ilast + 1
}
#cleanup
delete(@tmpdiv,@thisind,@thismult,@nestbits,@groupbits,silent:T)
delete(@nvars,@vnames,@nlevels,@userows,@offsets,@colcount,@xmat,@nterms)
delete(@tmpopt,@whichtrms,@whichvars,@ilast,@thismat,@tmp1,@i,@j,@ilast,silent:T)
structure(Rstr:delete(@Rstr,return:T),Fstr:delete(@Fstr,return:T),\
   nesting:delete(@nesting,return:T),grouping:delete(@grouping,return:T),\
   bitmodel:delete(@bitmodel,return:T),random:delete(@random,return:T),\
   rvars:delete(@rvars,return:T))
%buildRFStr%

====> makemat <====
makemat  MACRO  DOLLARS
) makemat is called from ems to produce various basis matrices
# usage: $S(info,term)
@random <- $1$random
@bitmodel <- $1$bitmodel
@Rstr <- $1$Rstr
@Fstr <- $1$Fstr
@nesting <- $1$nesting
@grouping <- $1$grouping
@unrest <- $1$unrest
@nonhier <- $1$nonhier
@nlevels <- $1$nlevels
@variate <- $1$variate
@thisterm <- $2
@bits <- @bitmodel[@thisterm]
@nvars <- length(@random)
@nterms <- length(@bitmodel)
@scalefix <- keyvalue($K,"scalefix","TF",default:T)

@tol <- 1e-5

@matout <- NULL
@which <- (2^(run(0,@nvars-1)) %& @bits) > 0
@varsntrm <- sum(@which)
@varsindx <- 0
@od <- setodometer(lower:0,upper:0)
if(@varsntrm > 0) {
	@varsindx <- (2^(run(0,@nvars-1)) %& @bits)[@which]
	@od <- setodometer(ndigits:@varsntrm,lower:0,upper:1,place:0)
}

if(sum(@which) == 0) {
	@fixed <- T
	@randomT <- F
} else {
	@fixed <- sum(@random[@which])==0
	@randomT <- sum(@random[@which]) == nbits(@bits)
}
@mixed <- !(@fixed || @randomT)
if(@mixed || @fixed) {
	@maxlvls <- 1
	for(@i,run(@nvars)) {
		if(@which[@i]) {
			@maxlvls <-* @nlevels[@i]
		}
	}
}

for(@kk,run(0,2^@varsntrm - 1)) {
	@k <- sum(@od$digits*@varsindx)
	@od <- setodometer(@od,step:1)
	@out <- rep(1,nrows(@Rstr[1]))
	@which <- (2^(run(0,@nvars-1)) %& @k) > 0
	@included <- (@k %& @bits) == @k
	@strictin <- nbits(@k %& @bits) < nbits(@bits)
	if(@thisterm == 1) {
		@termfit <- F
	} else {
		if(@nonhier) {
			@termfit <- sum(@k == @bitmodel[run(@thisterm-1)]) > 0
		} else {
			@termfit <- sum((@k %& @bitmodel[run(@thisterm-1)]) == @k) > 0
		}
	}

	# we don't need to include an earlier term literally if
	# 1, the term is not included in thisterm, or
	# 2, the term has already been fit, or
	# 3, thisterm is pure random, or
	# 4, thisterm has any random vars and the model is unrestricted, or
	# 5, the  model is nonhierarchical
	@skipterm <- !@included || @termfit ||\
			@randomT  && @strictin||\
			(@unrest && !@fixed && @strictin ) ||\
			(@nonhier && @k != @bits)

	if(@skipterm) {
		next
	}

	# how we actually include term @k depends a bit
	# if all vars in thisterm are fix or all are random, no
	# problem.  If thisterm is mixed, but @k is pure random,
	# then just treat it as random.
	if(@fixed) {
		for(@i,run(@nvars)) {
			if(@which[@i]) {
				@out <- colproduct(@out,@Fstr[@i])
			}
		}
		if(@scalefix && (@tmpdiv <- sqrt(@maxlvls/\
			if(sum(@which) > 0){
			prod(@nlevels[@which])
			} else {
				1
			})) > 1) {
			@out <- @out/@tmpdiv
		}
	} elseif (@randomT || (@unrest && !@fixed) ||\
		sum(!@random[@which]) == 0) {
			# treat term as random
			for(@i,run(@nvars)) {
				if(@which[@i]) {
					@out <- colproduct(@out,@Rstr[@i])
				}
			}
			if(@mixed) {
				# here we have a pure random subterm of a mixed term
				# rescale for matching with other subterms to be joined
				@out <- @out/sqrt(@maxlvls/ncols(@out))'
			}
	} else {
		# here mixed
		# the strategy is to build an N by abcd-? matix that
		# implies the right covariance structure (ie, coefs add
		# to zero across any fixed coordinate) for coefs
		# this matrix is the colproduct of all the Rstr matrices
		# for the vars in the term times an abcd by h orthog matrix
		# spanning the nonzero eigen vectors of the cov matrix
		# of the coefficients

		# get colproduct of all Rstr matrices for this term
		# this is a "too big" matrix that we will shrink if we
		# have a restricted mixed term
		@J <- rep(1,nrows(@Rstr[1]))
		for(@i,run(@nvars)) {
			if(@which[@i]) {
				@J <- colproduct(@J,@Rstr[@i])
			}
		}

 		# build up matrices to indicate levels of factors in an
 		# a x b x c ... matrix of interaction coefficients
 		# these are used in the construction of mixed effects bases
		@tmpopt <- getoptions(warning:T);setoptions(warning:F)
		@coefind <- structure(NULL)
		setoptions(warning:@tmpopt)
		@nlvl2 <- @nlevels
		for(@i,run(@nvars)) {
			if(!@which[@i]) {
				@nlvl2[@i] <- 1
			}
		}
		for(@i,run(@nvars)) {
			@rep1 <- prod(@nlvl2[run(@i)])/@nlvl2[@i]
			@rep2 <- prod(@nlvl2[run(@i,@nvars)])/@nlvl2[@i]
			@thislvl <- @nlvl2[@i]
			@tmp <- rep( rep(run(@thislvl),rep(@rep1,@thislvl)),@rep2 )
			@coefind <- changestr(@coefind,@i,1*(@tmp==run(@thislvl)'))
		}


		# get colproduct of all @coefind matrices for random vars in this term
		# and all vars grouped with a random var in this term
		@JR <- rep(1,nrows(@coefind[1]))
		for(@i,run(@nvars)) {
			if(@which[@i] && @random[@i]) {
				@JR <- colproduct(@JR,@coefind[@i])
				if(sum(@grouping[,@i]) > 0) {
					@j <- @i+1
					while(@grouping[@j,@i]) {
						@JR <- colproduct(@JR,@coefind[@j])
						@j <- @j+1
					}
				}
			}
		}


		# find variables in this term that are fixed and not variates
		@tmp <- @which %& !@random %& !@variate
		@fixedbits <- sum(2^(run(@nvars)-1) * @tmp)

		# now we loop through all the fixed vars in this
		# term, leaving out one at a time
		@thisout <- NULL
		for(@ii,run(@nvars)) {
			@partout <- @JR
			if( (2^(@ii-1) %& @fixedbits) > 0) {
				# var @ii is fixed in this term
				for(@j,run(@nvars)) {
					if((2^(@j-1) %& @fixedbits) > 0 && @j != @ii) {
						@partout <- colproduct(@partout,@coefind[@j])
					}
				}
				@thisout <- hconcat(@thisout,@partout)
			}
		}
		@out <- @thisout

		# OK, now we have a set of vectors across which our mixed effect
		# coefficients must add to zero
		# now find orthog basis for orthogoncal complement of this space
		if(!isnull(@out)) {
			@tmpe <- eigen(@out %*% @out')
			@H <- @tmpe$vectors[,@tmpe$values<@tol]
			@out <- @J %*% @H
		} else {
			@out <- @J
		}
		# rescale so all subterms of this mixed term match for variance
		@out <- @out/sqrt(@maxlvls/ncols(@J))'

	}
	@matout <- hconcat(@matout,@out)
}
 # send out just enough columns to span what we need
if(@mixed && !@unrest) {
	if(nrows(@matout) >= ncols(@matout)) {
		@tmpsvd <- svd(@matout,left:T,maxit:250)
		@matout <- @tmpsvd$leftvectors[,@tmpsvd$values>@tol]
	} else {
		@tmpsvd <- svd(@matout',right:T,maxit:250)
		@matout <- @tmpsvd$rightvectors[,@tmpsvd$values>@tol]
	}
	#@tmpe <- eigen(@matout %C% @matout) # @matout %*% @matout'
	#@matout <- @tmpe$vectors[,@tmpe$values>@tol]
	@matout <- @matout*sqrt(nrows(@matout)/@maxlvls)
	for(@i,run(@nvars)) {
		if(@which[@i] && @variate[@i]) {
			@matout <- @matout * sqrt(sum(@Rstr[@i]^2)/nrows(@matout))
		}
	}
}
structure(matrix:delete(@matout,return:T))
%makemat%

====> varcomp <====
varcomp    MACRO  DOLLARS
) This macro estimates variance components via the ANOVA method
) (ie, equate expected and observed mean squares).  varcomp can take
) two kinds of arguments.  First, it can take a single (structure)
) argument EMS, where EMS is the output from an
) ems(model,randomvars,keep:T[,...]) command.  Alternatively, you
) may give varcomp all the arguments you would give to ems, and
) varcomp will call ems internally and keep any needed output.
) The value of varcomp is a matrix giving the estimated variance
) components, their standard errors, and (approximate) degrees
) of freedom.
#usage $S(EMSoutput) or $S(model,randomvars[,ems options])
@EMS <- argvalue($1,"argment 1")
if(anytrue(!isstruc(@EMS), length(@tn <- compnames(@EMS)) != 5,\
	prod(@tn==vector("df", "ss", "termnames", "coefs", "rterms")) < .5)) {
	if(!ismacro(ems)) {
    		getmacros(ems,quiet:T,printname:F)
    		if(!ismacro(ems)) {
        		error("cannot find macro ems")
    		}
	}
	@EMS <- ems($0,keep:T)
}

@df <- @EMS$df
@ss <- @EMS$ss
@termnames <- @EMS$termnames
@coefs <- @EMS$coefs
@rterms <- vector(@EMS$rterms,T)
@finr <- sum( vector( @coefs[@rterms,!@rterms] ) > 0 ) > 0
if(@finr) {
	print("WARNING: fixed effects contribute to some random terms")
}
@vcfs <- solve(@coefs[@rterms,@rterms])
@msr <- (@ss/@df)[@rterms]
   #@vc <- solve(@coefs[@rterms,@rterms],(@ss/@df)[@rterms])
@vc <- @vcfs %*% @msr
@msrv <- 2*@msr^2/@df[@rterms]
@vcv <- (@vcfs^2) %*% @msrv
@vcdf <- @vc^2/((@vcfs)^2 %*% (@msr^2/@df[@rterms]))
matrix(hconcat(delete(@vc, return:T),sqrt(delete(@vcv, return:T)),\
	delete(@vcdf, return:T)),\
labels:structure(delete(@termnames, return:T)[delete(@rterms, return:T)],\
	vector("Estimate","SE","DF")))
%varcomp%

====> mixed <====
mixed    MACRO  DOLLARS
) This macro does mixed effects analysis of variance estimating
) approximate error MS where needed via the ANOVA method with a
) Satterthwaite approximate error df.
) There are two ways to use mixed.
)    EMSoutput <- ems(model, randomvars, keep:T [, ems options])
)    mixed(emsOutput [,mixed options])
) and
)    mixed(model, randomvars [,ems options] [,mixed options])
)
) The mixed options are useneg:T and/or keepmixed:T.
)
) Option useneg:T
)  By default, only linear combinations of mean squares with positive
)  coefficients are used in approximate tests.  When useneg:T is an
)  argument, the original numerator MS is used, but negative
)  multipliers may be used in the denominator.  Use of all nonnegative
)  coefficients tends to give better approximate chisquare
)  distributions.
)
) Option keepmixed:T
)  By default, mixed prints out a table giving each line of
)  the ANOVA with its df and MS, the error MS and approximate
)  error df, the F statistic and its p value.  In this case,
)  mixed returns a NULL value.  If the kewword phrase keepmixed:T
)  is given to mixed, the table will not be printed, but will
)  instead be returned as the value of mixed.  (The table is
)  returned as a labeled matrix.)
)) 980803 bug fix
)) 000207 stripped $$, made more use of argvalue() and keyvalue()
#usage $S(EMSoutput[,keepmixed:T][,useneg:T]) or
#      $S(model,randomvars[,ems options][,keepmixed:T][,useneg:T])
@EMS <- argvalue($1,"argument 1")
if(anytrue(!isstruc(@EMS), length(@tn <- compnames(@EMS)) != 5,\
	prod(@tn==vector("df", "ss", "termnames", "coefs", "rterms")) < .5)) {
	if(!ismacro(ems)) {
		getmacros(ems,quiet:T,printname:F)
		if(!ismacro(ems)) {
			error("cannot find macro ems")
		}
	}
	@EMS <- ems($0,keep:T)
}

@mkeep <- keyvalue($K,"keepmixed","TF", default:F)
@useneg <- keyvalue($K,"useneg*","TF",default:F)

@df <- @EMS$df
@ss <- @EMS$ss
@termnames <- @EMS$termnames
@coefs <- @EMS$coefs
@rterms <- vector(@EMS$rterms,T)

@nterms <- length(@df)
@finr <- sum( vector( @coefs[@rterms,!@rterms] ) > 0 ) > 0
if(@finr) {
	print("WARNING: fixed effects contribute to some random terms")
}

@cf2 <- @coefs
@cf2[run(@nterms)+@nterms*run(0,@nterms-1)] <- 0 # take out self
@cf2 <- @cf2'
@dcoef <- matrix(rep(0,@nterms^2),@nterms)
@dcoef[,@rterms] <- solve(@coefs[@rterms,@rterms]',@cf2[@rterms,])'
@ncoef <- dmat(rep(1,@nterms));
if(!@useneg && (sum(@use <- vector(@dcoef) < 0) > 0) ) {
	@ncoef[@use] <- abs(@dcoef[@use])
	@dcoef[@use] <- 0
}
@numer <- @ncoef %*% (@ss/@df)
@denom <- @dcoef %*% (@ss/@df)
@addf <- @denom^2/sum((@dcoef'*(@ss/@df))^2/@df)'
@andf <- @numer^2/sum((@ncoef'*(@ss/@df))^2/@df)'

@oldwarn <- getoptions(warning:T);setoptions(warning:F)
@result <- hconcat(@andf,@numer,@addf,@denom,@numer/@denom,\
1-cumF(@numer/@denom,@andf,@addf))
setoptions(warning:@oldwarn)

@result <- matrix(@result,labels:structure(@termnames,\
	vector("DF","MS","Error DF","Error MS","F","P value")))
if(@mkeep) {
	@result
} else {
	print(@result,format:"10.4g",header:F)
	NULL
}
%mixed%


====> interactplot <====
interactplot MACRO  DOLLARS
) interactplot() draws an interaction plot of marginal means of a
) response for all combinations of values of one or more factors
) interactplot(y,a [,b,c ...] [errors:T or errors:x,pool:T|F,errorvar:v]
)   [graphics keywords]), y a REAL vector, a, b, c, ... vectors of positive
)   integers the same length as y, x and v positive scalars.
) interactplot(means [errors:T or errors:x,errormat:s] [graphics
)   keywords]), means a REAL matrix or array,x a positive scalar, s a
)   positive matrix the same shape as means.
) interactplot(a [,b,c ...],frommodel:T,[errors:T or errors:x] [graphics
)   keywords]), a, b, c, ... factors in the current model, x a positive
)   scalar.
)) 021204 modified title when there are more than 1 non-keyword args
)) 050810 modified to do error bars and take from models
#$S(y,a,b,...[,graphics keywords])
@frommodel <- keyvalue($K,"frommodel","logic scalar",default:F)
@errors <- keyvalue($K,"errors")
if(isnull(@errors)) {
	@errorx <- 0
	@doerrors <- F
} elseif(isreal(@errors)) {
	if(!isscalar(@errors)) {
		error("when real, @errors must be a nonnegative scalar")
	}
	if(@errors < 0) {
		error("when real, @errors must be a nonnegative scalar")
	}
	@errorx <- @errors
	@doerrors <- T
} elseif(islogic(@errors)) {
	if(!isscalar(@errors)) {
		error("@errors must be a scalar")
	}
	@doerrors <- @errors
	@errorx <- 2
} else {
	error("errors must be T|F or a positive real scalar")
}
@pool <- keyvalue($K,"pool","logic scalar",default:F)
@errorvar <- keyvalue($K,"errorvar","positive real scalar")
@errormat <- keyvalue($K,"errormat","positive real array")
@frommodel <- keyvalue($K,"frommodel","logic scalar",default:F)
if(@doerrors) {
	@offwidth <- .15
} else {
	@offwidth <- 0
}

@mns <- argvalue($1,"argument 1","real")

if (!isvector(@mns)){
	if ($v > 1){
		error("more than 1 nonkeyword argument with non-vector argument 1")
	}
	@xlab <- "Dimension 1 of $1"
	@ylab <- "$1"
	@title <- "Elements of $1 against dimension 1"
	@xtickmax <- dim(@mns)[1]
	@mndim <- dim(@mns)
	@mns <- matrix(@mns, @xtickmax)
	if(@doerrors && isnull(@errormat)) {
		error("you must specify errormat:s when errors:T is used for a matrix of means")
	}
	if(isnull(@errormat) || !@doerrors) {
		@errormat <- 0*@mns
	} else {
		@errormat <- @errorx * @errormat
	}
	if(anytrue(\
		length(dim(@mns)) != length(dim(@errormat)),\
	   sum(dim(@mns) != dim(@errormat)) > 0)) {
		error("dimensions of mean matrix and error matrix do not match")
	}
} elseif(@frommodel) {
	if($v < 1) {
		error("$S() needs at least one model variable argument with frommodel:T",\
			macroname:F)
	}

	#check factor arguments
	@input <- structure($V)
	@names <- $A[run($v)]
	@term <- paste(@names,sep:".")
	@usenames <- $v > 1
	@mns <- glmtable(@term,estimate:T,seest:F)
	@mndim <- dim(@mns)
	@xtickmax <- @mndim[1]
	@mns <- matrix(@mns,@mndim[1])
	if(@doerrors) {
		@errormat <- glmtable(@term,seest:T,estimate:F)
		@errormat <- matrix(@errormat,@mndim[1]) * @errorx
	} else {
		@errormat <- 0*@mns
	}
	@xlab <- "Levels of $1"
	@ylab <- "Least squares means"
	@title <- paste("LS means of",@term,"against $1")
	if(@usenames) {
		@title <- paste(@title,"by",paste(@names[-1],sep:"."))
	}
} else {
	if($v < 2) {
		error("$S() needs at least one factor argument with vector argument 1",\
			macroname:F)
	}

	#check factor arguments
	@input <- structure($V)
	@names <- $A
	@usenames <- $v > 2
	for(@i,run(2,$v)) {
		argvalue(@input[@i],paste("factor",@i-1), "positive integer vector")
		if (@i > 2) {
			@usenames <- @usenames && isname(@names[@i])
		}
	}
	@names <- @names[-1]
	@mns <- tabs($V,mean:T)
	@mndim <- dim(@mns)
	@mns <- matrix(@mns,@mndim[1])
	if(@doerrors) {
		@cts <- matrix(tabs($V,count:T),max(@input[2]))
		if(isnull(@errorvar)) {
			@vars <- matrix(tabs($V,var:T),max(@input[2]))
			if(@pool) {
				@tmpv <- vector(@vars)
				@tmpn <- vector(@cts)
				@errorvar <- sum(@tmpv*(@tmpn-1))/sum(@tmpn-1)
				@errormat <- sqrt(@errorvar/@cts)
			} else {
				@errormat <- sqrt(@vars/@cts)
			}
		} else {
			@errormat <- sqrt(@errorvar/@cts)
		}
		@errormat <- @errorx * @errormat
	} else {
		@errormat <- 0*@mns
	}
	@title <- "Interaction plot of $1 vs $2"
	if (delete(@usenames,return:T)) {
		@title <- paste(@title,"by",paste(@names[-1],sep:"."))
	}
	@xlab <- "Levels of $2"
	@ylab <- "Means of $1"
	@xtickmax <- max(@input[2])
	delete(@input,@i,@names)
}

@offsets <- run(ncols(@mns))/ncols(@mns)
@offsets <- @offsets - sum(@offsets)/length(@offsets)
@offsets <- @offsets*@offwidth

@d <- @xtickmax/40

@basex <- run(nrows(@mns))+@offsets'
@hval <- vector(vconcat(@basex,rep(?,ncols(@mns))'))
@vval <- vector(vconcat(@mns,rep(?,ncols(@mns))'))

@tmpx <- vector(@basex)'
@tmpx <- vconcat(@tmpx,@tmpx,rep(?,length(@tmpx))')
@tmpy <- vconcat(vector(@mns-@errormat)',vector(@mns+@errormat)',\
	rep(?,length(@mns))')

@hval <- vector(@hval,@tmpx)
@vval <- vector(@vval,@tmpy)

chplot(@hval, @vval, lines:T, $K, symbols:"\7",\
	   xlab:@xlab, ylab:@ylab, title:@title,\
	   xticks:run(@xtickmax),xmin:1-@d,xmax:@xtickmax+@d,\
	   xticklen:-.5, yticklen:-.5, xaxis:F,yaxis:F,ticks:"LB")

@alllabs <- NULL
if(length(@mndim) > 1) {
	@nr <- nrows(@mns)
	@od <- setodometer(lower:1,upper:@mndim[-1])
	for(@i,run(ncols(@mns))) {
		@lab <- paste(@od$digits,sep:".")
		@alllabs <- vector(@alllabs,rep(@lab,@nr))
		@od <- setodometer(@od,step:1)
		#addstrings(@basex[,@i],@mns[,@i],rep(@lab,@nr))
	}
	addstrings(vector(@basex),vector(@mns),@alllabs)
}

delete(@xlab,@ylab,@title,@xtickmax,@mns,@d)
%interactplot%


====> interblock <====
interblock MACRO  DOLLARS
) Macro to perform interblock recovery of information in incomplete
) block designs.  The usage is
) interblock(y,block,trt,[contrast:coefs])
) where y is the vector of responses, block is the factor of block
) levels, trt is the factor of treatment levels and coefs is an
) optional vector of contrast coefficients (it must be the same
) length as the number of levels of trt).  If a contrast is not
) specified, the output is the intra, inter, and combined treatment
) effects and their standard errors.  If a contrast is specified,
) the output is the intra, inter, and combined estimates of the
) contrasts with their standard errors.
) 11-09-1999 GWO
))Version 000207 uses argvalue, keyvalue, strip $$
# $S(y,block,trt,[contrast:coefs])
if($v < 3) {
	error("Usage: $S(y,block,trt,[contrast:coefs])",macroname:F)
}
@y <- argvalue($1,"argument 1","real vector")
@block <- argvalue($2,"argument 2","positive integer vector")
@trt <- argvalue($3,"argument 3","positive integer vector")

@ccfs <- keyvalue($K,"con*","nonmissing real vector")
@docon <- !isnull(@ccfs)
if (@docon){
    if(abs(sum(@ccfs)) > 1e-7*sum(abs(@ccfs)) ) {
		error("Contrast argument not a REAL vector that sums to 0")
    }
}

if(!isfactor(@block) || !isfactor(@trt)) {
	error("Arguments 2 and 3 not both factors")
}
if(nrows(@y) != nrows(@block) || nrows(@block) != nrows(@trt)) {
	error("Arguments 1, 2, and 3 must have the same length")
}

@aovmodel <- paste(nameof(@y),"=",nameof(@block),"+",nameof(@trt))
anova(@aovmodel,silent:T)

@intracfs <- coefs(3) #coefs(nameof(@trt))
@intraV1 <- secoefs(3)$se^2 #secoefs(nameof(@trt))$se^2
@g <- length(@intracfs)

if(@docon) {
    if(@g != nrows(@ccfs)) {
		error("Length of contrast argument != number of levels of argument 3")
    }
    @intrac <- contrast(nameof(@trt),@ccfs)
    @intraest <- @intrac$estimate
    @intraV2 <- @intrac$se^2
}

@vout <- varcomp(@aovmodel, nameof(@block), marg:T)

@ct <- tabs(@y,@block,count:T)
@V <- dmat(@vout[1,1]*@ct^2 + @vout[2,1]*@ct)
@btot <- tabs(@y,@block,mean:T)*@ct
@inc <- tabs(@y,@block,@trt,count:T)
@incg <- @inc[,@g]
@inc <- @inc[,-@g]
@use <- @incg == 1
@inc[@use,] <- @inc[@use,]-1
@inc <- hconcat(rep(1,nrows(@inc)),@inc)
@cfmat <- hconcat(rep(0,@g-1),dmat(rep(1,@g-1)))
@cfmat <- vconcat(@cfmat,vector(0,rep(-1,@g-1))')
	#print(@vout[1,1],@vout[2,1])
@xtxinvt <- solve(@inc %c% @inc,singok:T) 
if(isnull(@xtxinvt)) {
	error("cannot compute all interblock estimates")
}
@xtxinvt <- @xtxinvt %*% @inc'
@intercfs <- @cfmat %*% @xtxinvt %*% @btot
@interV <- @cfmat %*% @xtxinvt %*% @V %*% (@xtxinvt %c% @cfmat')
@interV1 <- diag(@interV)

if(@docon) {
    @interest <- sum(@intercfs*@ccfs)
    @interV2 <- @ccfs %c% @interV %*% @ccfs
    @intraw <- (1/@intraV2)/(1/@intraV2 + 1/@interV2)
    @combest <- @intraw*@intraest + (1-@intraw)*@interest
    @combse <- sqrt(1/(1/@intraV2 + 1/@interV2))
    @out <- matrix(vector(@intraest,@interest,@combest,@intrac$se,\
		@interV2^.5,@combse),3,\
	    labels:structure(vector("intra","inter","combined"),\
			vector("estimate","se")))
} else {
    @intraw <- (1/@intraV1)/(1/@intraV1 + 1/@interV1)
    @combest <- @intraw*@intracfs + (1-@intraw)*@intercfs
    @combse <- sqrt(1/(1/@intraV1 + 1/@interV1))
    @out <- matrix(vector(@intracfs, sqrt(@intraV1), @intercfs, \
		sqrt(@interV1), @combest, @combse), @g, \
		labels:structure(rep("", @g), \
		vector("intra est", "intra se", "inter est", "inter se", \
			"combined est", "combined se")))
}
@out
%interblock%

sidebyside MACRO DOLLARS
) Produce a sidebyside plot.  There must be an active ANOVA model.
) sidebyside([termlaby:y,labels:c,rescale:tf,showconst:tf,boxcut:int])
) Optional arguments are
) termlaby:real      Specify height for term labels. This can be a single
)                    value or a vector of length equal to number of terms.
) labels:charvector  Specify your own term labels.
) rescale:logic      Should effects be divided by their standard errors.
)                    Default is F.  This uses the standard errors as
)                    reported by secoefs() and may not be correct for
)                    all mixed models (secoefs() depends on terms labeled
)                    ERRORX).  Residuals are divided by root MSE.
) showconst:logic    Should the coefficient of the CONSTANT be shown?
)                    Default is F.
) boxcut:integer     Cutoff for using a boxplot for a term instead of
)                    plotting individual effects.  Default is 20.
)) 020829 modified by kb; default labels are now drawn as xticklabs
))        ymin and ymax are determined slightly differently
# usage: $S([termlaby:y,labels:c,rescale:tf,showconst:tf,boxcut:int])
@labyval <- keyvalue($K,"termlaby","real vector")
@userlabs <- keyvalue($K,"labels","char")
@stndrdz <- keyvalue($K,"rescale","logic",default:F)
@showconst <- keyvalue($K,"showconst","logic",default:F)
@boxcut <- keyvalue($K,"boxcut*","real",default:20)

if(@stndrdz) {
	@cfs <- coefs()/secoefs(coefs:F)
	@mse <- SS[length(SS)]/DF[length(DF)]
} else {
	@cfs <- coefs()
}
@nt <- ncomps(@cfs)
@names <- compnames(@cfs)
@ilvl <- 0
@outx <- NULL
@outy <- NULL
@outc <- NULL
@labs <- NULL
for(@i,run(@nt)) {
	if(@showconst || (@names[@i] != "CONSTANT")) {
		@ilvl <-+ 1
		@this <- vector(@cfs[@i])
		@outx <- vector(@outx,rep(@ilvl,length(@this)))
		@outy <- vector(@outy,@this)
		@outc <- vector(@outc,run(length(@this)))
		@labs <- vector(@labs,@names[@i])
	}
}
@ilvl <-+ 1
@outx <- vector(@outx,rep(@ilvl,length(RESIDUALS)))
@outc <- vector(@outc,run(length(RESIDUALS)))
@labs <- vector(@labs,"RESIDUALS")
if(!isnull(@userlabs)) {
	if(length(@userlabs) != length(@labs)) {
		error(paste("you must supply",length(@labs),"labels"))
	}
	@labs <- @userlabs
}
if(@stndrdz) {
	@outy <- vector(@outy,RESIDUALS/sqrt((1-HII)*@mse))
} else {
	@outy <- vector(@outy,RESIDUALS)
}

@tmp <- tabs(@outy,@outx,count:T)
if(max(@tmp) < @boxcut) {
	chplot(@outx,@outy,@outc,xlab:" ",xmin:0.5,xmax:@ilvl+.5,xticks:NULL,xaxis:F,\
		ylab:{
			if(@stndrdz){
				"Std. Effect or Resid"
			}else{
				"Effect or Residual"
			}
		},show:F)
} else {
	@boxn <- run(length(@tmp))[@tmp>=@boxcut]
	@boxuse <-  vector(sum(@boxn == @outx') > 0)
	@oldwarn <- getoptions(warnings:T)
	setoptions(warnings:F)
	boxplot(split(@outy[@boxuse],@outx[@boxuse]),show:F,vertical:T,xticks:NULL)
	setoptions(warnings:delete(@oldwarn,return:T))
	addchars(@outx[!@boxuse],@outy[!@boxuse],@outc[!@boxuse],show:F)
	showplot(xmin:0.5,xmax:@ilvl+.5,title:" ",xlab:" ",\
		ylab:{
			if(@stndrdz){
				"Std. Effect or Resid"
			}else{
				"Effect or Residual"
			}
		},show:F)
	delete(@boxuse,@boxn)
}
@ymin <- min(@outy) - .05*(max(@outy) - min(@outy))
@ymax <- max(@outy) + .05*(max(@outy) - min(@outy))
if(isnull(@labyval)) {
	showplot(xticks:run(@ilvl),xticklab:@labs,show:F)
} else {
	if(length(@labyval)==1) {
		@labyval <- rep(@labyval,length(@labs))
	}
	if(length(@labyval) != length(@labs)) {
		error(paste("length of label height must be either 1 or",length(@labs)))
	}
	addstrings(run(length(@labs)),@labyval,@labs,show:F)
}
@model <- vecread(string:STRMODEL,bychars:T)
@J <- match("=",@model) + run(2)
if (paste(@model[@J],sep:"") == "1+") {
	@model <- @model[-@J]
}
delete(@outx,@outy,@outc,@labyval,@labs,@ilvl)
showplot(notakey:NULL,$K,ymin:delete(@ymin,return:T),\
		 ymax:delete(@ymax,return:T),\
		 title:paste("Side by side plot for \"",\
					 delete(@model,return:T),"\"",sep:""))
%sidebyside%



====> findpower <====
findpower     MACRO  DOLLARS
) findpower(means,nis,sigma2,alpha) computes the power for the F-test
) in the one-way anova testing the null hypothesis of no treatment
) differences when the means, sample sizes, error variance, and type I
) error rate are as given.
if($v != 4){
	error("usage is $S(means,nis,sigam2,alpha)", macroname:F)
}
@rcb <- keyvalue($K,"rcb","logic scalar",default:F)
if(length($1) > 1) {
	@means <- argvalue($1,"means","real vector nonmissing")
	@nis <- argvalue($2,"sample sizes","integer vector positive")
	@sigma2 <- argvalue($3,"error variance","real scalar positive")

	if(length(@means) != length(@nis)) {
		error("lengths of mean and sample size arguments differ")
	}
	if(@rcb && length(unique(@nis)) > 1) {
		error("sample sizes must be equal for rcb:T")
	}
	@wtmean <- sum(@nis*@means)/sum(@nis)
	@noncen <- sum(@nis*(@means-@wtmean)^2)/@sigma2
	@df1 <- length(@means)-1
	@df2 <- sum(@nis)-length(@means)
	if(@rcb) {
		@df2 <- @df2 - @nis[1] + 1
	}
} else {
	@noncen <- argvalue($1,"noncentrality parameter","real scalar nonnegative")
	@df1 <- argvalue($2,"numerator df","positive scalar")
	@df2 <- argvalue($3,"denominator df","positive scalar")
	if(@rcb) {
		error("cannot use rcb with noncentrality form")
	}
}
@alpha <- argvalue($4,"argument 4","real scalar positive")
if(@alpha >= 1) {
	error("type 1 error rate must be between 0 and 1")
}

power2(@noncen,@df1,@df2,@alpha)
%findpower%

====> findsampsize <====
findsampsize     MACRO  DOLLARS
) findsampsize(means,sigma2,alpha,reqpower) computes the sample size
) required to get power for the F-test in the one-way anova at least
) as large as reqpower.  The returned value is a structure with the
) vector of sample sizes and the attained power.  In this usage, the
) sample sizes are assumed to be equal.
)
) findsampsize(ncp1,ngrps,alpha,reqpower) computes the sample size
) required to get power for the F-test in the one-way anova at least
) as large as reqpower.  Here, ncp1 is the n=1 noncentrality 
) parameter and ngrps is the number of groups.
)
) findsampsize([means,sigma2] or [ncp1,ngrps],alpha,reqpower,rcb:T) 
) computes the sample size required to get power for the F-test in
) a randomized complete block at least as large as reqpower.  
)
) findsampsize(means,sigma2,alpha,reqpower,prop:propvec) does similarly
) except that the group sample sizes are assumed to be (approximately)
) proportional to propvec.  Sample sizes are actually round(X*propvec),
) so X will be at least 1/min(propvec)/2 to get at least one
) observation in each group.  

if($v != 4){
	error("usage is $S(means,sigma2,alpha,reqpower) or $S(ncp1,ngrps,alpha,reqpower)", macroname:F)
}
if(length($1) > 1) {
	@means <- argvalue($1,"means","real vector nonmissing")
	@sigma2 <- argvalue($2,"variance","real scalar positive")
	@g <- length(@means)
} else {
	@ncp1 <- argvalue($1,"noncentrality parameter","real scalar nonmissing")
	@g <- argvalue($2,"number of groups","count positive")
}
@alpha <- argvalue($3,"alpha","real scalar positive")
@reqpower <- argvalue($4,"power","real scalar positive")

if(@alpha >= 1) {
	error("type 1 error rate must be between 0 and 1")
}
if(@reqpower >= 1) {
	error("requested power must be between 0 and 1")
}
@propvec <- keyvalue($K,"prop","vector real positive nonmissing",\
	default:rep(1,@g))
@rcb <- keyvalue($K,"rcb","logic scalar",default:F)

@df1 <- @g-1

@iter <- 1
@rightX <- 1/min(@propvec)/3.999
@rightnis <- round(@rightX*@propvec)
@rightp <- 0
while(@rightp < @reqpower && @iter < 20) {
	@leftp <- @rightp
	@leftX <- @rightX
	@leftnis <- @rightnis
	@rightX <- 2*@rightX
	@rightnis <- round(@rightX*@propvec)
	if(isdefined(@ncp1)) {
		@noncen <- @ncp1*@rightnis[1]
	} else {
		@wtmean <- sum(@rightnis*@means)/sum(@rightnis)
		@noncen <- sum(@rightnis*(@means-@wtmean)^2)/@sigma2
	}
	@df2 <- sum(@rightnis)-@g
	if(@rcb) {
		@df2 <- @df2 - @rightnis[1] + 1
	}
	if(@df2 > 0) {
		@rightp <- power2(@noncen,@df1,@df2,@alpha)
	}
}
if(@rightp < @reqpower) {
	errorout <- paste("sample sizes truncated at ", @rightnis,\
		" with attained power ", @rightp)
	error(errorout)
}
@lasttryp <- 0
while(sum((@rightnis-@leftnis)^2) > 0.5) {
	@tryX <- (@rightX + @leftX)/2
	@trynis <- round(@tryX*@propvec)
	if(isdefined(@ncp1)) {
		@noncen <- @ncp1*@trynis[1]
	} else {
		@wtmean <- sum(@trynis*@means)/sum(@trynis)
		@noncen <- sum(@trynis*(@means-@wtmean)^2)/@sigma2
	}
	@df2 <- sum(@trynis)-@g
	if(@rcb) {
		@df2 <- @df2 - @trynis[1] + 1
	}
	if(@df2 > 0) {
		@tryp <- power2(@noncen,@df1,@df2,@alpha)
	}
	if(@tryp > @reqpower) { # replace right side
		@rightX <- @tryX
		@rightnis <- @trynis
		@rightp <- @tryp
	} else {
		@leftX <- @tryX
		@leftnis <- @trynis
		@leftp <- @tryp
	}
	if(abs(@lasttryp-@tryp) < 1e-6) {
		break
	}
	@lasttryp <- @tryp
}
structure(nis:@rightnis,power:@rightp)
%findsampsize%

====> findncp <====
findncp     MACRO  DOLLARS
) findncp(means,nis,sigma2) computes the noncentrality for the F-test
) in the one-way anova testing the null hypothesis of no treatment
) differences when the means, sample sizes, and error variance
) are as given.

if($v != 3){
	error("usage is $S(means,nis,sigam2)", macroname:F)
}
@means <- argvalue($1,"argument 1","real vector nonmissing")
@nis <- argvalue($2,"argument 2","integer vector positive")
@sigma2 <- argvalue($3,"argument 3","real scalar positive")

if(length(@means) != length(@nis)) {
	error("lengths of mean and sample size arguments differ")
}

@wtmean <- sum(@nis*@means)/sum(@nis)
sum(@nis*(@means-@wtmean)^2)/@sigma2
%findncp%

====> docontrast <====
docontrast  MACRO DOLLARS
if($v != 2) {
	error("$S requires two non-keyword arguments")
}
@term <- argvalue($1,"term name","char scalar")
@coefs <- argvalue($2,"contrast coefficients","real")
@by <- keyvalue($K,"by*","char scalar")
@error <- keyvalue($K,"error*","char scalar")
if(isnull(@by) && isnull(@error)) {
	@cout <- contrast(@term,@coefs)
} elseif(isnull(@error)) {
	@cout <- contrast(@term,@coefs,<<@by>>)
} elseif(isnull(@by)) {
	@cout <- contrast(@term,@coefs,error:@error)
} else {
	@cout <- contrast(@term,@coefs,<<@by>>,error:@error)
}
if(isnull(@error)) {
	@df <- DF[length(DF)]
} else {
	@df <- DF[match(@error,TERMNAMES)]
}
@nc <- length(@cout$estimate)
if(@df < .05) {
	@out <- hconcat(@cout$estimate',@cout$ss,rep(?,@nc)',rep(0,@nc)')
	setlabels(@out,structure("(",vector("estimate","ss","se","df")))
	return(@out)
}
@cover <- keyvalue($K,"cover*","real scalar",default:.95)
if(@cover <= 0 || @cover >= 1) {
	error("Coverage must be between 0 and 1")
}
@out <- hconcat(@cout$estimate'',@cout$se'',rep(@df,@nc)'')
if(!isnull(keyvalue($K,"twosided*","logic scalar"))) {
	if(isnull(@cover)) {
		error("You must supply cover:xx for a confidence interval")
	}
	@pm <- invstu((1-@cover)/2,@df)*vector(1,-1)'*@cout$se''
	@out <- hconcat(@out,@out[,1]+@pm)
	setlabels(@out,structure("(",vector("estimate","se","df",\
		"lowerb","upperb")))
	return(@out)
}
if(!isnull(keyvalue($K,"lowerb*","logic scalar"))) {
	if(isnull(@cover)) {
		error("You must supply cover:xx for a confidence interval")
	}
	@pm <- invstu((1-@cover),@df)*@cout$se''
	@out <- hconcat(@out,@out[,1]+@pm)
	setlabels(@out,structure("(",vector("estimate","se","df",\
		"lowerb")))
	return(@out)
}
if(!isnull(keyvalue($K,"upperb*","logic scalar"))) {
	if(isnull(@cover)) {
		error("You must supply cover:xx for a confidence interval")
	}
	@pm <- invstu(@cover,@df)*@cout$se''
	@out <- hconcat(@out,@out[,1]+@pm)
	setlabels(@out,structure("(",vector("estimate","se","df",\
		"upperb")))
	return(@out)
}
@out <- hconcat(@cout$estimate'',@cout$se'',@cout$ss'',rep(@df,@nc)'')
if(!isnull(keyvalue($K,"twotail*","logic scalar"))) {
	@t <- @out[,1]/@out[,2]
	@pval <- twotailt(@t,@df)
	@out <- hconcat(@out,@pval'')
	setlabels(@out,structure("(",vector("estimate","se","ss","df",\
		"pvalue")))
	return(@out)
}
if(!isnull(keyvalue($K,"lowertail*","logic scalar"))) {
	@t <- @out[,1]/@out[,2]
	@pval <- cumstu(@t,@df)
	@out <- hconcat(@out,@pval'')
	setlabels(@out,structure("(",vector("estimate","se","ss","df",\
		"pvalue")))
	return(@out)
}
if(!isnull(keyvalue($K,"uppertail*","logic scalar"))) {
	@t <- @out[,1]/@out[,2]
	@pval <- 1-cumstu(@t,@df)
	@out <- hconcat(@out,@pval'')
	setlabels(@out,structure("(",vector("estimate","se","ss","df",\
		"pvalue")))
	return(@out)
}
%docontrast%


====>dorandtest<====
dorandtest MACRO DOLLARS
@type <- keyvalue($K,"type","char scalar")
@trials <- keyvalue($K,"trials","positive integer scalar")
if(@type == "twosample") {
	if($v != 2) {
		error("You must specify two real vector arguments when type is "twosample"")
	}
	@x <- argvalue($1,"x","real vector")
	@y <- argvalue($2,"y","real vector")
	if(isnull(@trials)) {
		@dist <- randt2(@x,@y)
	} else {
		@dist <- randt2(@x,@y,trials:@trials)
	}
	@val <- describe(@x,mean:T)-describe(@y,mean:T)
	@xlab <- "Difference of means"
} elseif(@type == "sign") {
	if($v != 1) {
		error("You must specify one real vector argument")
	}
	@x <- argvalue($1,"x","real vector")
	if(isnull(@trials)) {
		@dist <- randsign(@x,@y)
	} else {
		@dist <- randsign(@x,@y,trials:@trials)
	}
	@xlab <- "Sum of differences"
	@val <- sum(@x)
} else {
	error("Type must be one of twosample or sign")
}
@nd <- length(@dist)
@lowerpv <- sum(@dist - @val < 1e-7)/@nd
@upperpv <- sum(@dist - @val > -1e-7)/@nd
@twotpv <- sum(abs(@dist)-abs(@val) > -1e-7)/@nd
@upper <- keyvalue($K,"upper*","logic scalar",default:F)
@lower <- keyvalue($K,"lower*","logic scalar",default:F)
@twot <- keyvalue($K,"twot*","logic scalar",default:F)
if((@upper+@lower+@twot) != 1) {
	error("You must specify exactly one of upper:T, lower:T, or twot:T")
}
if(@upper) {
	@out <- vector(@val,@upperpv)
} elseif (@lower) {
	@out <- vector(@val,@lowerpv)
} else {
	@out <- vector(@val,@twotpv)
}
setlabels(@out,vector("statistic","pvalue"))
@hist <- keyvalue($K,"hist*","logic scalar",default:F)
if(@hist) {
	@w <- 3.5*describe(@dist,var:T)^.5/length(@dist)^.333
	@k <- ceiling(max(abs(@dist))/@w)
	@tmp <- hist(@dist,run(-@k*@w,@k*@w,@w),save:T)
	@ym <- max(@tmp$y)*1.05
	@xp <- vector(@tmp$x,?,rep(@val,2))
	@yp <- vector(@tmp$y,?,vector(0,@ym))
	lineplot(@xp,@yp,xlab:" ",ylab:"Density",yaxis:F,\
		title:"Randomization distribution")
}
@out
%dorandtest%


====>stdordlabels<====
stdordlabels MACRO DOLLARS
# stdordlabels(k:val) produces A, B, AB, C, ... for val factors
# stdordlabels(letters:"def") produces d, e, de, f, ...
@k <- keyvalue($K,"k","positive integer scalar")
@letters <- keyvalue($K,"letters","character scalar")
if(isreal(@k)) {
	if(@k > 15) {
		error("A factor count must be a positive integer scalar less than 16")
	}
	@letters <- vector("A","B","C","D","E","F","G","H","K","J","L","M","N","O","P")[run(@k)]
} elseif(ischar(@letters)) {
	@letters <- vecread(string:@letters,byword:T)
	@letters <- paste(@letters,sep:"")
	@letters <- vecread(string:@letters,bychar:T)
	@k <- length(@letters)
	if(@k < 1 || @k > 15) {
		error("Between 1 and 15 factor letters can be used")
	}
} else {
	error("No factor count or letter list given")
}
if(length(unique(@letters)) < length(@letters)) {
	error("Factor letters must be unique")
}
@labs <- NULL
for(@i,run((2^@k)-1)) {
	@which <- (@i %& 2^run(0,@k-1)) > 0
	@labs <- vector(@labs,paste(@letters[@which],sep:""))
}
@labs
%stdordlabels%

yatesplot MACRO DOLLARS
# yatesplot(y) does a halfnormal plot of factor effects for 
# data y from a 2 series design in standard order
# yatesplot(y,fullnorm:T) uses a fullnormal plot rather than
# a halfnormal
# yatesplot(y,letters:"defg..") uses the letters to signify
# factors rather than the default A, B, C, ...
@y <- argvalue($1,"data","real vector nonmissing")
@k <- log2(length(@y))
if(abs(@k-round(@k)) > 1e-7) {
	error("Length of data must be a power of 2")
}
@yts <- yates(@y)
@dofull <- keyvalue($K,"full*","logic scalar",default:F)
if(@dofull) {
	@scores <- rankits(@yts)
	@xlab <- "Normal scores"
	@ylab <- "Effects"
	@title <- "Normal scores plot of effects"
} else {
	@scores <- halfnorm(@yts)
	@yts <- abs(@yts)
	@xlab <- "Half-normal scores"
	@ylab <- "Absolute effects"
	@title <- "Half-Normal plot of absolute effects"
}
if(@k > 10.5) {
	chplot(@scores,@yts,"\27",$K,xlab:@xlab,ylab:@ylab,title:@title)
	return(NULL)
}
@lts <- keyvalue($K,"let*","character scalar")
if(!isnull(@lts)) {
	@labs <- stdordlabels(letters:@lts)
} else {
	@labs <- stdordlabels(k:@k)
}
if(length(@labs) != length(@scores)) {
	error("wrong number of factor letters")
}
@xmax <- (1+(1+@dofull)*@k/75)*max(@scores)
if(@dofull) {
	@xmin <- (1+(1+@dofull)*@k/75)*min(@scores)
} else {
	@xmin <- 0
}
chplot(@scores,@yts," ",$K,xlab:@xlab,ylab:@ylab,title:@title,\
	xmax:@xmax,xmin:@xmin)
addstrings(@scores,@yts,@labs,wind:0)
return(NULL)
%yatesplot%


aberration2  MACRO   DOLLARS 
) aberration2(basis) where basis is a k by p matrix
) of 0s and 1s (or -1s) indicating generators for the
) design.  aberration2 compute the vector of alias and
) then the aberration, the number of aliases of each
) length
@basis <- argvalue($1,"basis","integer matrix")
@basis <- abs(@basis)
if (max(vector(@basis)) > 1){
    error("basis not matrix with elements 0, -1, +1")
}
@k <- ncols(@basis)
@p <- nrows(@basis)
@basis <- @basis %*% 2^(run(@k-1,0))

@thisalias <- vector(0,@basis[1])
if(@p > 1) {
	for(@i,run(2,@p)) {
		@thisalias <- vector(@thisalias %^ vector(0,@basis[@i])')
	}
}
@thisalias <- @thisalias[-1]
@ab <- tabs(,nbits(@thisalias))
padto(@ab,@k)
%aberration2%



doconfound2 MACRO DOLLARS
@p <- keyvalue($K,"p","positive integer scalar")
@k <- keyvalue($K,"k","positive integer scalar")
@cfdout <- keyvalue($K,"cfdout","structure")
@basis <- keyvalue($K,"basis","nonnegative integer matrix")
@doce <- keyvalue($K,"confeff","TF",default:F)
@doaber <- keyvalue($K,"aber*","TF",default:F)
@doassign <- keyvalue($K,"assign*","TF",default:F)
@_cfdout <- NULL
if (!ismacro(aliases2)){
	getmacros(aliases2, quiet:T,printname:F)
	if (isnull(aliases2)){
		error("cannot proceed without macro aliases2")
	}
}
if (!ismacro(confound2)){
	getmacros(confound2, quiet:T,printname:F)
	if (isnull(confound2)){
		error("cannot proceed without macro confound2")
	}
}
if (!ismacro(aberration2)){
	getmacros(aberration2, quiet:T,printname:F)
	if (isnull(aberration2)){
		error("cannot proceed without macro aberration2")
	}
}
if(!isnull(@cfdout)) {
	@ce <- aliases2(@cfdout$basis)[-1]
	@aber <- @cfdout$aberration
	@assign <- confound2(@cfdout$basis)
} elseif(!isnull(@basis)) {
	@ce <- aliases2(@basis)[-1]
	@aber <- aberration2(@basis)
	@assign <- confound2(@basis)
} elseif(!isnull(@p) && !isnull(@k)) {
	if(@k - @p < 1) {
		error("You must have fewer blocks than treatments")
	}
	if(@k > 8) {
		print("Warning: non-cataloged design; this could take a while.")
		@tries <- keyvalue($K,"tries","count")
		if(isnull(@tries)) {
			@_cfdout <- choosedef2(@k,@p,all:T)
		} else {
			@_cfdout <- choosedef2(@k,@p,tries:@tries)
		}
	} else {
		@_cfdout <- read("Design.dat.txt",paste("cfd",@k,@p,sep:"_"),silent:T)
	}
	@ce <- aliases2(@_cfdout$basis)[-1]
	@aber <- @_cfdout$aberration
	@assign <- confound2(@_cfdout$basis)
} else {
	error("must supply basis, or confounding output, or p and k")
}
if(@doce) {
	print(name:"Confounded effects",@ce)
}
if(@doaber) {
	print(name:"Aberration",@aber)
}
if(@doassign) {
	print(name:"Block assignments",@assign)
}
@_cfdout
%doconfound2%



doff2 MACRO DOLLARS
@p <- keyvalue($K,"p","positive integer scalar")
@k <- keyvalue($K,"k","positive integer scalar")
@alout <- keyvalue($K,"alout","structure")
@basis <- keyvalue($K,"basis","nonnegative integer matrix")
@doIal <- keyvalue($K,"Ialiases","TF",default:F)
@doallal <- keyvalue($K,"allal*","TF",default:F)
@doaber <- keyvalue($K,"aber*","TF",default:F)
@dofrac <- keyvalue($K,"showfr*","TF",default:F)
@dorand <- keyvalue($K,"random*","TF",default:F)
@_alout <- NULL
if(!isnull(@alout)) {
	@Ialiases <- aliases2(@alout$basis)
	@allalias <- allaliases2(@alout$basis)
	@aber <- @alout$aberration
	@frac <- ffdesign2(@alout$basis)
} elseif(!isnull(@basis)) {
	@Ialiases <- aliases2(@basis)
	@allalias <- allaliases2(@basis)
	@aber <- aberration2(@basis)
	@frac <- ffdesign2(@basis)
} elseif(!isnull(@p) && !isnull(@k)) {
	if(@k - @p < 2) {
		error("You must have at least four runs")
	}
	if(@k > 12) {
		print("Warning: non-cataloged design; this could take a while.")
		@tries <- keyvalue($K,"tries","count")
		if(isnull(@tries)) {
			@_alout <- choosegen2(@k,@p,all:T)
		} else {
			@_alout <- choosegen2(@k,@p,tries:@tries)
		}
	} else {
		@_alout <- read("Design.dat.txt",paste("al",@k,@p,sep:"_"),\
			silent:T)
	}
	if(@dorand) {
		@basis <- @_alout$basis
		for(@i,run(@k-@p+1,@k)) {
			@basis[,@i] <- @basis[,@i]*(2*(runi(1) < .5) - 1)
		}
		@_alout <- changestr(@_alout,"basis",@basis)
	}
	@Ialiases <- aliases2(@_alout$basis)
	@allalias <- allaliases2(@_alout$basis)
	@aber <- @_alout$aberration
	@frac <- ffdesign2(@_alout$basis)
} else {
	error("must supply basis, or confounding output, or p and k")
}
if(@doIal) {
	print(name:"Aliases of I",@Ialiases)
}
if(@doallal) {
	print(name:"Aliases",@allalias)
}
if(@doaber) {
	print(name:"Aberration",@aber)
}
if(@dofrac) {
	print(name:"Treatments",@frac)
}
@_alout
%doff2%


