meansanova     macro     dollars
) Macro to do multi-way anova using means, SD's and n's
) Usage:
)  meansanova(xbar, s, n [,Model] [,anova keywords])
)   xbar        REAL vector, matrix or array of means with no
)               MISSING values
)   s           vector, matrix or array of positive std deviations,
)               the same size and shape as xbar
)   n           positive integer or vector, matrix or array of positive
)               integers the same size and shape as xbar 
)               a scalar is equivalent to
)               array(rep(n,length(xbar)),dims(xbar))
)   Model       CHARACTER scalar in the form of an ANOVA model, all of
)               whose variable names start with "@".  The default is
)               "@Y=@GROUPS" (xbar a vector) or "@Y=@A*@B..." otherwise
)               There can be no more than 8 factors
)
)  meansanova() 
)     Creates a non-random dependent variable with group sizes,
)       means and standard deviations matching n, xbar and s
)     Creates ndims(xbar) factors
)     Runs anova(Model ...)
)     Deletes RESIDUALS and WTDRESIDUALS
)  The response and factor names are those in Model
)  meansanova can be followed by coefs(), secoefs(), contrast()
)  and modelinfo(), but not resvsrankits(), resvsyhat() and resvsindex()
)  The usual side effect variables are created except RESIDUALS and
)  WTDRESIDUALS
)) 030503 Written by C. Bingham, kb@umn.edu
@ybars <- argvalue($1,"means","nonmissing real")
if (isvector(@ybars)) {
	@ybars <- vector(@ybars)
}

@sds <- argvalue($2,"std devs","positive real")
if (isvector(@sds)) {
	@sds <- vector(@sds)
}

@n <- argvalue($3,"sample sizes","positive integer")
if (isvector(@n)) {
	@n <- vector(@n)
}
@usemodel <- $v >= 4
@model <- if (!@usemodel) {
	NULL
} else {
	argvalue($4,"model","string")
}
@dims <- dim(@ybars)
@ndims <- ndims(@ybars)
if (anytrue(ndims(@sds) != @ndims,sum(dim(@sds) != @dims) != 0)) {
	error("Sizes or shapes of means and sd's differ")
}
if (isscalar(@n)) { # equal sample sizes
	@n <- array(rep(@n,length(@ybars)),@dims)
} elseif (anytrue(ndims(@n) != @ndims,sum(dim(@n) != @dims) != 0)) {
	error("Sizes or shapes of n's and sd's differ")
}
if (@ndims > 8) {
	error("More than 8 dimensions are illegal")
}
if (@usemodel) {
	@names <- varnames(@model)
	if (length(@names) != @ndims+1) {
		error("Wrong number of variables in the model")
	}
	for (@i,1,length(@names)) {
		if (match("@*",@names[@i],0,exact:F) == 0) {
			error("All names in the model must start with '@'")
		}
	}
	@respname <- @names[1]
	@facnames <- @names[-1]
	delete(@names)
} else {
	@respname <- "@Y"
	@facnames <- if (@ndims == 1) {
		"@GROUPS"
	} else {
		vector("@A","@B","@C","@D","@E","@F","@G","@H")[run(@ndims)]
	}
}

if (!ismacro(buildfactor)) {
	getmacros(buildfactor,silent:T)
}

<<@respname>> <- rep(0,sum(vector(@n)))
@op <- "*"
if (!@usemodel) {
	@model <- paste("@Y=",@facnames[1],sep:"")
	for(@i,1,@ndims) {
		@model <- if (@i == 1) {
			paste(@respname,@facnames[@i],sep:"=")
		} else {
			paste(@model,@op,@facnames[@i],sep:"")
		}
	}
}
@N <- vector(@n)
for(@i,1,@ndims) {
	<<@facnames[@i]>> <- factor(rep(buildfactor(@i,@dims),@N))
}
@ncells <- prod(@dims)
@start <- 0
for(@i,1,@ncells) {
	@J <- @start + run(@N[@i])
	<<@respname>>[@J] <- @ybars[@i] +\
		  @sds[@i]*(run(@N[@i]) - (@N[@i]+1)/2)/sqrt(@N[@i]*(@N[@i]+1)/12)
	@start <-+ @N[@i]
}
if ($k == 0) {
	anova(@model)
} else {
	anova(@model,$K)
}
delete(RESIDUALS,WTDRESIDUALS,silent:T)
for (@i,1,@ndims) {
	delete(<<@facnames[@i]>>)
}

delete(@ybars,@sds,@n,@ncells,<<"@Y">>,@facnames,@start,@i,@J)
%meansanova%