info    MACRO
) Regression related macros
)   anovapred)   Compute anova cell means and standard errors.
)   betalimits   Compute confidence limits for a regression coefficient
)   testbeta     Compute t-statistic and optionally DF and
)                P-value for test of H0: betaj = hypValue
)   estimlimits  compute confidence limits for E(y | x)
)   predlimits   compute confidence limits for y for specified x
)   testestim    compute t-statistic and optional DF and P-value
)                for test of H0: E(y | x) = hypValue
)   stepsetup    Initialize stepwise regression
)   entervar     Enter variable in stepwise regression
)   regcoefs     Compute labelled matrix of regression coefficients,
)                with standard errors, t-statistics, P-values
)   regs         Carry out regression of y on columns of a matrix
)   nlreg        Non linear regression
)   removevar    Remove variable in stepwise regression
)   resid        Compute various case statistics related to residuals
)   _resvsxxxx   Back end for next three macros
)   resvsindex   Plot studentized residuals against case number
)   resvsrankits Plot studentized residuals against normal scores
)   resvsyhat    Plot studentized residuals against predicted values
)   stepstatus   Print report on current state of stepwise regression
)   steplook     Retrieve information on current state of stepwise
)                regression
)   yhat         Compute various case statistics related to fitted
)                values
) The following are macros defining test functions to be fit by nlreg
)   linear       Linear function f <- x %*% b
)   asymptot     f = b[1] + b[2]*(b[3])^x & derivatives
)   DandSFunc    f = b[1] + (.49 - b[1])*exp(-b[2]*(x - 8))   & derivatives
)   testfun      f = b[1]+b[2]*x+b[3]*x^2+b[4]*sin(b[5]*x)+b[6]*exp(-b[7]*x)
)
) The following are test data for nlreg() (read using read())
)   SandCT19.8.1 Data from Snedecor and Cochran for testing asymptot()
)   DandS_T10.2  Data from Draper & Smith for testing DandSFunc()
)
) Help for most macros follows them in the file
)) Copyright (C) 2003 by Christopher Bingham and Gary W. Oehlert
)) Version of 001004
)) Version of 001210 added subtopics to help
)) Version of 001216 edited topics and subtopics
)) Version of 010516 fixed bug in stepwise macros so that they work when
))       model has temporary variables
)) Version of 010519 fixed bug introduced by bug fix
)) Version of 010710 changed default symbol in resvsrankits,resvsindex,
))       resvsyhat back to "\6"
)) Version of 010721, resvsrankits and resvsindex work after arima
)) Version of 010724, keyword usehii added to resvs{rankits,index,yhat}
)) Version of 010725  resvsrankits(), resvsindex(), resvsyhat() now
))       differ only in their names and are now merely front ends to a
))       new invisible macro _resvsxxxx()
)) Version of 011212 fixed bug in _resvsxxxx()
)) 020806 Amended some uses of cumstu() and cumF() to use 'upper:T'
)) 030320 Added nlreg() and related macros and data sets (moved from
))   arima.mac)
)) 030329 Improved help for nlreg()
)) 030430 Corrected and simplified help for resvsxxx() macros.
)) 030814 added subtopic titles to help
%info%

regcoefs      MACRO  DOLLARS
) Macro to compute labeled matrices of regresson coefficients,
) their standard errors, t-statistics and optionally P-values,
) in a form similar to regress() output
) Usage:
)  regcoefs(model [,pvals:T] [,byvar:F])
)  regcoefs([pvals:T] [,byvar:F])
)   model     CHARACTER scalar, a GLM model
) If the response is multivariate, the result is a structure of
) of matrices, one for each variable.
) Without pvals:T, no P-values are computed.
) byvar:F is meaningful only if the response is multivariate when
) the result is one big matrix instead of a structure of matrices.
)) 990309 Converted to using argvalue()
)) 990901 minor change in arguments to keyvalue()
)) 000218 Stripped $$, replaced changestr() by str[i] <- comp
)) 020806 use keyword 'upper' on cumstu()
# $S([pvals:T][,byvar:T]) or $S(model[,pvals:T][,byvar:T])
if ($v > 1 || $k > 2){
	error("usage: $S([pvals:T][,byvar:T]) or $S(model[,pvals:T][,byvar:T])",\
		macroname:F)
}
@pvals <- keyvalue($K,"pval*","TF",default:getoptions(pvals:T))
@byvar <- keyvalue($K,"byvar*","TF")
if ($v == 1){
	@model <- argvalue($1,"model","string")
	manova(delete(@model, return:T),silent:T)
}
@p <- nrows(SS[1,,])
@mterms <- dim(SS)[1]
if (isnull(@byvar)){
	@byvar <- (@p != 1)
}elseif(@p == 1 && @byvar){
	print("WARNING: byvar:T ignored by $S with univariate response variable")
	@byvar <- F
}
@stuff <- secoefs(byterm:F)
if (max(vector(length(@stuff[1]))) > @p){
	error("$S does not allow factors in the model", macroname:F)
}
if (isstruc(@stuff[1])){
	@rowlabels <- compnames(@stuff[1])
}else{
	@rowlabels <- TERMNAMES[1]
}
@collabels <- if (@p == 1 || @byvar){
	vector("Coef","StdErr","t")
}else{
	@makelabs <- macro("getlabels(vector(run(\\$1),labels:\\$2))",inline:F)
	vector(@makelabs(@p,"Coef "),\
	@makelabs(@p,"StdErr "),@makelabs(@p,"t "))
}
if (DF[@mterms] == 0){
	@pvals <- F
}
if (@pvals){
	@collabels <- if (@p == 1 || @byvar){
		vector(@collabels,"P-Value")
	}else{
		vector(@collabels,@makelabs(@p,"P-Value "))
	}
}
delete(@makelabs,silent:T)
@coefs <- vector(@stuff[1])
@stderrs <- vector(@stuff[2])
if (@p > 1){
	@coefs <- matrix(@coefs,@p)'
	@stderrs <- matrix(@stderrs,@p)'
}
@tstats <- if (DF[@mterms] > 0){
	@coefs/@stderrs
}else{
	@pvals <- F
	array(rep(?,(@mterms-1)*@p),dim(@coefs))
}
if (!@byvar || @p == 1){
	@o <- hconcat(@coefs,@stderrs,@tstats)
	if (@pvals){
		@o <- hconcat(@o,2*cumstu(abs(@tstats),DF[@mterms],upper:T))
	}
	@o <- matrix(@o,labels:structure(@rowlabels,@collabels))
}else{
	@o <- split(rep(0,@p)', compnames:"Variable_") #create structure
	for (@i,run(@p)){
		@oi <- hconcat(@coefs[,@i],@stderrs[,@i],@tstats[,@i])
		if (@pvals){
			@oi <- hconcat(@oi,\
				2*cumstu(abs(@tstats[,@i]),DF[@mterms],upper:T))
		}
		@o[@i] <- matrix(@oi,labels:structure(@rowlabels,@collabels))
	}
	delete(@oi,@i)
}
delete(@stuff,@p,@mterms,@pvals,@byvar,@coefs,@stderrs,@tstats)
delete(@o,return:T)
%regcoefs%

betalimits   MACRO DOLLARS
) Macro to compute confidence limits for a regression coefficient
) Its use must be preceded by regress("y=x") or
) regress("y=x1 + x2 + ... xk"))
) It also works after anova()
) Usage:
)  betalimits(term, level)
)   term      positive integer, CHARACTER scalar, variable name
)             like the first argument to secoefs()
)   level     positive REAL scalar < 1, confidence level (1 -
)             confidence level when level > .5)
)
)  After regress("y=x1 + x2 + x3"), any of the following will work
)   betalimits(x2, .95), betalimits("x2",) or betalimits(3)
)   (counting the constant term, x2 is the third term)
) If level < .5, a warning message is printed and 1 - level is used
)
) After anova(), when the term is a main effect, the result is a matrix
) with two columns and a row for each level
) When the term is an interaction of k factors, the result is a k+1
) dimensional array, the first k dimensions matching the dimensions of
) the interaction and the last being of length 2
)) Version 010721
)) 020806 use 'upper:T' on invstu()
# $S(term,level)
if (isnull(modelinfo(strmodel:T,nomodelok:T)) || !isdefined(DF)) {
	error("Apparently no active model")
}
@stuff <- secoefs($1)
@level <- argvalue($2,"confidence level", "positive scalar")
if (@level >= 1){
	error("Confidence level >= 1")
}
if (@level < .5){
	@level <- 1 - @level
	print(paste("WARNING: Confidence level < .5; will use 1 - level =",\
		@level))
}
@se <- @stuff[2]
@errmargin <- invstu((1 - delete(@level,return:T))/2,\
			reverse(DF)[1],upper:T)*@se
@plusminus <- if (isscalar(@se)) {
	vector(-1,1)
} else {
	vector(-1,1)'
}
@retval <- vector(@stuff[1]) + @plusminus * vector(@errmargin)
if (!isscalar(@se)){
	@retval <- array(@retval,dim(@se),2)
}
delete(@stuff,@errmargin, @plusminus, @se)
delete(@retval,return:T)
%betalimits%

testbeta     MACRO DOLLARS
) Macro to compute a t-statistic for a hypothesis of the form
) H0: betaj = hypValue, where betaj is a regression coefficient
) Its use must be preceded by regress("y=x") or
) regress("y=x1 + x2 + ... xk")
) Usage:
)   testbeta(term, hypValue [,df:T] [,pval:T])
)   term      positive integer, CHARACTER scalar, variable name
)             like the first argument to secoefs()
)   hypValue  REAL scalar, hypothesized value of coeficient
)
) If neither df:T or pval:T is an argument, the t-statistic is
) returned.  Otherwise a structure with 2 or elements is returned
) with components tstat and one or both of df (error DF from regression)
) or pvalue (two-tail P-value for t-statistic).
)) Version 010721
)) 020806 use 'upper:T' on cumstu()
# $S(term, hypValue [,df:T] [,pval:T])
if (!isdefined(COEF) || !isdefined(SS) || !isdefined(DF) || !isdefined(XTXINV)){
	error("Apparently no active model")
}
@hypValue <- argvalue($2, "hypothesized value","number")
@stuff <- secoefs($1)
@df <- keyvalue($K,"df", "TF",default:F)
@pval <- keyvalue($K,"pval*", "TF",default:F)
@tstat <- (@stuff[1] - @hypValue)/@stuff[2]
if (@df || @pval){
	@tstat <- structure(tstat:@tstat)
	@fe <- reverse(DF)[1]
	if (@df){
		@tstat <- strconcat(@tstat,df:@fe)
	}
	if (@pval){
		@tstat <- strconcat(@tstat,pvalue:2*cumstu(abs(@tstat[1]),@fe,upper:T))
	}
	delete(@fe)
}
delete(@stuff,@hypValue,@df)
delete(@tstat,return:T)
%testbeta%

estimlimits  MACRO DOLLARS
) Macro to compute confidence limits for the expected
) value E(y|x) after a regression
) Its use must be preceded by regress("y=x") or
) regress("y=x1 + x2 + ... xk")
) or more generally by anova(Model) where Model may contain both
) one or more predictors (covariates) and/or one or more factors
)
) Usage after regress() or anova() with no factors
)  estimlimits(x, level), where x is like the
)  argument to regpred() and level is the confidence level
)
) Usage after anova() with covariates and factors
)  estimlimits(x,factorValues, level), where x and factorValues
)  are like arguments to glmpred()
)
) Usage after anova() with no covariates
)  estimalimits(NULL,factorValues, level)
)
)  After regress("y=x1 + x2 + x3"), for example
)   estimlimits(vector(2,3,4), .95) returns vector(lower,upper)
)   where lower and upper are the limits when x1=2, x2=3 and x3=4
)
)   estimlimits(hconcat(x01, x02, x03),.95), returns
)   hconcat(lower,upper), where lower[I] and upper[I]
)   are limits whenx1 = x01[I], x2 = x02[I] amd x3 = x03[I].
)   x01, x02, and x03 must all be vectors of the same length.
)
)  After anova("y = a + b")
)   estimlimits(NULL,vector(1,2),.95) returns vector(lower,upper), where
)   lower and upper are limits when a = 1 and b = 2
)
)  After anova("y = x + a + b")
)   estimlimits(3,vector(1,2),.95) returns vector(lower,upper), where
)   limits are those for x = 3, a = 1 and b = 2
) If level < .5, a warning message is printed and 1 - level is used
)) Version 000904
)) 020806 use 'upper:T' on invstu()
# $S(xvalues [,factorLevels], level)
if ($v < 2 || $v > 3){
	error("Usage: $S(xvalues [,factorLevels] [,confidenceLevel]")
}
if (isnull(modelinfo(strmodel:T,nomodelok:T)) || !isvector(DF)){
	error("Apparently no active model")
}
@nvariates <- modelvars(nvariates:T)
@nfactors <- modelvars(nfactors:T)
if (@nfactors > 0 && $v == 2){
	error("you must supply factor level(s) when there are factors in model")
}
@x <- argvalue($01,"x value(s)")
if (@nvariates == 0) {
	if (!isnull(@x)){
		error("argument 1 must be NULL when there are no variates in model")
	}
} else {
	if (isnull(@x)){
		error("you must supply variate value(s) when there are variates in model")
	}
	argvalue(@x,"x value(s)","nonmissing real")
	if (!ismatrix(@x)){
		error("x values must be vector or matrix")
	}
	if (@nvariates == 1){
		if (!isvector(@x)){
			error("x values must be vector when there is one variate")
		}
	} elseif (isvector(@x) && length(@x) != @nvariates || \
		!isvector(@x) && ncols(@x) != @nvariates) {
		error("number of variate values doesn't match model")
	}
}

@level <- if ($v == 2){
	@factors <- NULL
	$2
} else {
	@factors <- argvalue($2,"factor level(s)","positive integer")
	if (@nfactors == 0){
		error("Don't supply factor levels when model has no factors")
	}
	if (!ismatrix(@factors)){
		error("Factor levels must be vector or matrix")
	}
	if (isvector(@factors) && length(@factors) != @nfactors ||\
		!isvector(@factors) && ncols(@factors) != @nfactors) {
		error(paste("number of factor levels !=",@nfactors))
	}
	$3
}
delete(@nfactors,@nvariates)
argvalue(@level,"confidence level", "positive scalar")
if (@level >= 1){
	error("Confidence level >= 1")
}
if (@level < .5){
	@level <- 1 - @level
	print(paste("WARNING: Confidence level < .5; will use 1 - level = ",@level))
}

@stuff <- glmpred(delete(@x,return:T),delete(@factors,return:T))
@critval <- invstu((1-@level)/2,reverse(DF)[1],upper:T)
@lower <- @stuff[1] - @critval*@stuff[2]
@upper <- @stuff[1] + @critval*@stuff[2]
delete(@critval,@stuff)
@result <- if (isscalar(@lower)){
	vector(@lower,@upper)
}else{
	hconcat(@lower,@upper)
}
delete(@lower,@upper)
delete(@result,return:T)
%estimlimits%

predlimits   MACRO DOLLARS
) Macro to compute prediction limits for y = E(y|x) + epsilon
) Its use must be preceded by regress("y=x") or
) regress("y=x1 + x2 + ... xk")
) Usage:  predlimits(x, level), where x is like the
)  argument to regpred() and level is the confidence level
)
)  After regress("y=x1 + x2 + x3"), for example
)   predlimits(vector(2,3,4), .95) returns vector(lower,upper)
)   where lower and upper are the limits when x1=2, x2=3 and x3=4
)
)   predlimits(hconcat(x01, x02, x03)), returns
)   hconcat(lower,upper), where lower[I] and upper[I]
)   are limits whenx1 = x01[I], x2 = x02[I] amd x3 = x03[I].
)   x01, x02, and x03 must all be vectors of the same length.
)
) If level < .5, a warning message is printed and 1 - level is used
)) Version 000904
)) 020806 use 'upper:T' on invstu()
# $S(x, level)
if (isnull(modelinfo(strmodel:T,nomodelok:T)) || !isvector(DF, real:T)) {
	error("Apparently no active regression model")
}
if (modelvars(nfactors:T) > 0){
	error("Previous model contains one or more factors")
}
@x <- argvalue($1,"x value(s)","nonmissing real")
if (!ismatrix(@x)){
	error("x values must be vector or matrix")
}
@level <- argvalue($2,"confidence level", "positive scalar")
if (@level >= 1){
	error("Confidence level >= 1")
}
if (@level < .5){
	@level <- 1 - @level
	print(paste("WARNING: Confidence level < .5; will use 1 - level = ",@level))
}

@stuff <- regpred(delete(@x,return:T))
@critval <- invstu((1-@level)/2,reverse(DF)[1],upper:T)
@lower <- @stuff[1] - @critval*@stuff[3]
@upper <- @stuff[1] + @critval*@stuff[3]
delete(@level,@critval,@stuff)
@result <- if (isscalar(@lower)){
	vector(@lower,@upper)
}else{
	hconcat(@lower,@upper)
}
delete(@lower,@upper)
delete(@result,return:T)
%predlimits%

testestim    MACRO DOLLARS
) Macro to test a hypothesis H0: E(y|x) = hypValue
) Its use must be preceded by regress("y=x") or
) regress("y=x1 + x2 + ... xk")
) Usage:  testestim(x, hypValue), where hypValue is the
) hypothesized value, x is a scalar (simple linear
) regression) or a vector of length k, and level is
) the confidence level
)
) Return value:
)  After regress("y=x1 + x2 + x3"), for example
)  testestim(vector(2,3,4), 17) returns the t-statistic for
)  testing H0: E(y|x) = 17
)  testestim(vector(2,3,4), 17, pval:T, df:T) returns
)  structure(tstat:t_statistic, df:ErrorDF, pval:p_value)
)) Version 000904
)) 020806 use 'upper:T' on cumstu()
# $S(x, hypValue [,df:T] [,pval:T])
if (isnull(modelinfo(strmodel:T,nomodelok:T)) || !isvector(DF,real:T)){
	error("Apparently no active model")
}
if (modelvars(nfactors:T) > 0){
	error("Model has at least one factor")
}
@x <- argvalue($1,"x value","nonmissing real vector")
if (!ismatrix(@x)){
	error("x values must be vector or matrix")
}
@hypValue <- argvalue($2, "hypothesized value","number")
@df <- keyvalue($K,"df", "TF",default:F)
@pval <- keyvalue($K,"pval*", "TF",default:F)
@stuff <- regpred(delete(@x,return:T))
@tstat <- (@stuff[1] - delete(@hypValue,return:T))/@stuff[2]
if (@df || @pval){
	@tstat <- structure(tstat:@tstat)
	@fe <- reverse(DF)[1]
	if (@df){
		@tstat <- strconcat(@tstat,df:@fe)
	}
	if (@pval){
		@tstat <- strconcat(@tstat,\
			pvalue:2*cumstu(abs(@tstat[1]),@fe,upper:T))
	}
	delete(@fe)
}
delete(@stuff,@df, @pval)
delete(@tstat,return:T)
%testestim%

_resvsxxxx         MACRO DOLLARS
) Back end to macros resvsindex(), resvsrankits() and resvsyhat(),
) all of which are now identical and use this common macro.
) Usage:
)   _resvsxxxx(which ,usehii , standres ,var ,C [,graphics keyword phrases])
)   name          unquoted name of calling macro (resvsindex,resvsrankits,
)                 or resvsyhat)
)   var           Positive integer <= ncols(RESIDUALS)
)   usehii        LOGICAL scalar or NULL
)   standres      LOGICAL scalar
)   C             CHARACTER or non-negative integer scalar or vector
) In case of error, a CHARACTER scalar error message is returned;
) otherwise NULL is returned
)
) See resvsindex(), resvsrankits(), resvsyhat() for more info
) Graphics keywords are passed by $K in calling macro; this may include
) keywords 'usehii' and 'standres' but they are ignored here
) Care has been taken so that error messages don't mention _resvsxxxx.
) In particular all use of argvalue() and keyvalue() is in front end macros
)) 010725 new macro
)) 011212 fixed bug.  Several uses of length() changed to nrows() to
))        take into account the multivariate case.
#$S(which, usehii, standres ,var ,C [,graphics keywords])
@name <- "$1"
@usehii <- $2
@standres <- $3

if (match("resvsindex*",@name,0,exact:F) != 0){
	@which <- 1
	@xname <- "Case Numbers"
} elseif (match("resvsrankit*",@name,0,exact:F) != 0) {
	@which <- 2
	@xname <- "Normal Scores"
} elseif (match("resvsyhat*",@name,0,exact:F) != 0) {
	@which <- 3
	@xname <- "Fitted Values (Yhat)"
} else {
	error(paste("\"",@name,"\" improper first argument",sep:""))
}
@name <- paste(@name,"()",sep:"")

@var <- $4

if (!isscalar(@var,positive:T,integer:T)){
	return("Argument 1 not positive integer")
}

@chars <- $5 # already checked

if (!isdefined(RESIDUALS)){
	return("Apparently no linear, GLM or ARIMA model has been fitted")
}

@n <- nrows(RESIDUALS)
if (alltrue(isdefined(ALLRESIDUALS),nrows(ALLRESIDUALS) >= @n)){
	@J <- run(nrows(ALLRESIDUALS) - @n + 1, nrows(ALLRESIDUALS))
	@arma <- sum(RESIDUALS != ALLRESIDUALS[@J]) == 0
	delete(@J)
} else {
	@arma <- F
}

if (@var > ncols(RESIDUALS)) {
	return("argument 1 > ncols(RESIDUALS)")
}

if (isnull(@usehii)) {
	@usehii <- alltrue(!@arma, @standres, isvector(HII,real:T,nonneg:T))
} elseif (@usehii && !@standres) {
	print(paste("WARNING: usehii:T ignored with stand:F in",@name))
	@usehii <- F
}

if (@usehii && anytrue(!isvector(HII,real:T,nonneg:T),length(HII) != @n)){
	return("no HII or HII is not a REAL vector => 0 of right length")
}

if (@arma){
	if (@which == 3){
		return(paste(@name, "can't be used after arima()"))
	}
	if (@standres) {
		if (!isdefined(NPAR)){
			return("Need NPAR to be defined after arima() with standres:T")
		}
		@sd <- sqrt(sum(RESIDUALS^2)/(@n - NPAR))
	}
} elseif (@standres) {
	if(!isdefined(DF) || !isdefined(SS)) {
		return("SS and/or DF not defined.  No linear model fitted?")
	}
	@m <- length(DF)
	if (DF[@m] == 0){
		return("no degrees of freedom for error; can't standardize")
	}
	@sd <- sqrt(SS[@m,@var,@var]/DF[@m])
	delete(@m)
}

if (isnull(@chars)) {
	@chars <- "\6"
} elseif (alltrue(isscalar(@chars, real:T), @chars == 0)){
	@chars <- "###"
}

@r <- if(!@arma && isdefined(WTDRESIDUALS)){
	WTDRESIDUALS[,@var]
}else{
	RESIDUALS[,@var]
}

@w <- getoptions(warnings:T)
setoptions(warnings:F)

@kind <- if (@standres) {
	if (@usehii) {
		@sd <-* sqrt(1 - HII)
	}
	@r <-/ @sd
	delete(@sd,@usehii)
	"Standardized"
} else {
	""
}

@x <- if (@which == 1) { #resvsindex
	1
} elseif (@which == 2) { # resvsrankits
	rankits(@r)
} else { # resvsyhat
	modelvars(y:T)[,@var] - RESIDUALS[,@var]
}
setoptions(warnings:@w)

chplot(delete(@x,return:T),@r,symbols:@chars,$K,yaxis:F,\
	title:paste(@kind,"Residuals vs",@xname),\
	xlab:@xname,ylab:paste(@kind,"Resids"))
delete(@w,@var,@r,@chars,@kind,@standres)
%_resvsxxxx%

resvsindex         MACRO DOLLARS
) Macro to plot standardized residuals computed from variable
) RESIDUALS or WTDRESIDUALS against case number
) Usage:
)  resvsindex(var [,usehii:T or F] [, stand:F] [,graphics keyword phrases])
)   var       positive integer
)   stand:F   don't standardize
)  resvsindex([keyword phrases])
)   is equivalent to resvsindex(1 [,keyword phrases])
)
)  The default for usehii is T after GLM commands and F after arima.
)  usehii:T is ignored with stand:F
)
) It may be used after all OLS GLM commands (regress(), anova(), and
) manova()), most non-linear GLM commands such as logistic() and
) poisson(), and an ARIMA fit using macro arima().  The only time you
) need to provide var is after manova() when you want a residual plot of
) column var > 1 of RESIDUALS, the matrix of residuals.
)
) With stand:F, raw unstandardized residuals are plotted.
)
) Without stand:F, when usehii is T (default after GLM commands), residuals
) are standardized by dividing RESIDUALS by sqrt(mse*(1 - HII)), where HII
) are the leverages computed by the GLM command or arima().
)
) Without stand:F, when usehii is F (default after arima()), residuals are
) standardized simply by dividing RESIDUALS by sqrt(mse).
)
) In both cases, mse is the estimated residual variance after a linear GLM
) command, the mean error deviance after a nonlinear GLM command or the
) estimated innovation variance after arima().
)
) You can use any of the usual graphics keyword phrases.  In particular,
) you can change the plotting symbol using 'symbols' and change default
) axis labels and title using 'xlab', 'ylab' and 'title'.
)
) symbols:0 is special.  It specifies case numbers be used as plotting
) symbol.
)
)) Version of 951227, works properly after weighted analysis.
)) 000220 stripped $$, use argvalue()
)) 010710 changed default symbol to "\6"; make symbol:chars work
))        also it now works after arima()
)) 010724 added keyword usehii
)) 010725 drastically shortened; now it's a frontend for _resvxxxxx
))        and should be identical to resvsrankits and resvsyhat
#$S([usehii:T or F] [,standres:T][,graphics keywords]) or $S(var [,keywords])

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

@standres <- keyvalue($K,"stand*","TF",default:T)
@usehii <- keyvalue($K,"usehii","TF")

@chars <- keyvalue($K,"symbol*","vector")

@var <- if($v == 0){
	1
}else{
	argvalue($1, "argument 1", "positive integer scalar")
}

if ($v > 2) {
	error("at most 2 non-keyword arguments allowed")
} elseif ($v == 2){ #resvsindex(i,c) allowed for backward compatibility
	if (!isnull(@chars)){
		error("at most 1 non-keyword argument with 'symbols'")
	}
	@chars <- argvalue($2, "argument 2", "vector")
	if (!ischar(@chars) && !isreal(@chars)){
		error("argument 2 not REAL or CHARACTER")
	}
}elseif (!isnull(@chars) && !isreal(@chars)  && !ischar(@chars)){
	error("value of 'symbols' not REAL or CHARACTER")
}

@msg <- _resvsxxxx($S, @usehii, @standres, @var, @chars, $K)
if (!isnull(@msg)){
	error(@msg, macroname:match("resvs*",@msg,0,exact:F) == 0)
}
delete(@usehii,@standres,@var,@chars,@msg)
%resvsindex%

resvsrankits    MACRO DOLLARS
) Macro to plot standardized residuals computed from variable
) RESIDUALS or WTDRESIDUALS computed by a GLM comand or arima()
) against normal scores (rankits)
) Usage:
)  resvsrankits(var [,usehii:T or F] [, stand:F] [,graphics keyword phrases])
)   var       positive integer
)   stand:F   don't standardize
)  resvsrankits([keyword phrases])
)   is equivalent to resvsrankits(1 [,keyword phrases])
)
)  The default for usehii is T after GLM commands and F after arima.
)  usehii:T is ignored with stand:F
)
) It may be used after all OLS GLM commands (regress(), anova(), and
) manova()), most non-linear GLM commands such as logistic() and
) poisson(), and an ARIMA fit using macro arima().  The only time you
) need to provide var is after manova() when you want a residual plot of
) column var > 1 of RESIDUALS, the matrix of residuals.
)
) With stand:F, raw unstandardized residuals are plotted.
)
) Without stand:F, when usehii is T (default after GLM commands), residuals
) are standardized by dividing RESIDUALS by sqrt(mse*(1 - HII)), where HII
) are the leverages computed by the GLM command or arima().
)
) Without stand:F, when usehii is F (default after arima()), residuals are
) standardized simply by dividing RESIDUALS by sqrt(mse).
)
) In both cases, mse is the estimated residual variance after a linear GLM
) command, the mean error deviance after a nonlinear GLM command or the
) estimated innovation variance after arima().
)
) You can use any of the usual graphics keyword phrases.  In particular,
) you can change the plotting symbol using 'symbols' and change default
) axis labels and title using 'xlab', 'ylab' and 'title'.
)
) symbols:0 is special.  It specifies case numbers be used as plotting
) symbol.
)
)) Version of 951227, works properly after weighted analysis.
)) 000220 stripped $$, uses argvalue()
)) 010710 changed default symbol to "\6"; make symbol:chars work
)) 010720 it now works after arima()
)) 010724 added keywords stand and usehii
)) 010725 drastically shortened; now it's a frontend for _resvxxxxx
))        and should be identical to resvsindex and resvsyhat
#$S([usehii:T or F] [,standres:T][,graphics keywords]) or $S(var [,keywords])

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

@standres <- keyvalue($K,"stand*","TF",default:T)
@usehii <- keyvalue($K,"usehii","TF")

@chars <- keyvalue($K,"symbol*","vector")

@var <- if($v == 0){
	1
}else{
	argvalue($1, "argument 1", "positive integer scalar")
}

if ($v > 2) {
	error("at most 2 non-keyword arguments allowed")
} elseif ($v == 2){ #resvsindex(i,c) allowed for backward compatibility
	if (!isnull(@chars)){
		error("at most 1 non-keyword argument with 'symbols'")
	}
	@chars <- argvalue($2, "argument 2", "vector")
	if (!ischar(@chars) && !isreal(@chars)){
		error("argument 2 not REAL or CHARACTER")
	}
}elseif (!isnull(@chars) && !isreal(@chars)  && !ischar(@chars)){
	error("value of 'symbols' not REAL or CHARACTER")
}

@msg <- _resvsxxxx($S, @usehii, @standres, @var, @chars, $K)
if (!isnull(@msg)){
	error(@msg, macroname:match("resvs*",@msg,0,exact:F) == 0)
}
delete(@usehii,@standres,@var,@chars,@msg)
%resvsrankits%

resvsyhat          MACRO DOLLARS
) Macro to plot standardized residuals computed from variable
) RESIDUALS or WTDRESIDUALS against fitted values
) Usage:
)  resvsyhat(var [,usehii:F] [, stand:F] [,graphics keyword phrases])
)   var       positive integer
)   stand:F   don't standardize
)  resvsyhat([keyword phrases])
)   is equivalent to resvsyhat(1 [,keyword phrases])
)
)  usehii:T is ignored with stand:F
)
) It may be used after all OLS GLM commands (regress(), anova(), and
) manova()) and most non-linear GLM commands such as logistic() and
) poisson().  It can't be used following an ARIMA fit using macro
) arima().  The only time you need to provide var is after manova() when
) you want a residual plot of column var > 1 of RESIDUALS, the matrix of
) residuals.
)
) With stand:F, raw unstandardized residuals are plotted.
)
) Without stand:F, when usehii is T (default), residuals are standardized
) by dividing RESIDUALS by sqrt(mse*(1 - HII)), where HII are the leverages
) computed by the GLM command or arima().
)
) Without stand:F, when usehii is F, residuals are standardized simply by
) dividing RESIDUALS by sqrt(mse).
)
) In both cases, mse is the estimated residual variance after a linear GLM
) command, or the mean error deviance after a nonlinear GLM command.
)
) You can use any of the usual graphics keyword phrases.  In particular,
) you can change the plotting symbol using 'symbols' and change default
) axis labels and title using 'xlab', 'ylab' and 'title'.
)
) symbols:0 is special.  It specifies case numbers be used as plotting
) symbol.
)
)) Version of 951227, works properly after weighted analysis.
)) 000220 stripped $$, uses argvalue()
)) 010710 changed default symbol to "\6"; make symbol:chars work
)) 010724 added keyword usehii
)) 010725 drastically shortened; now it's a frontend for _resvxxxxx
))        and should be identical to resvsrankits and resvsindex
#$S([usehii:T or F] [,standres:T][,graphics keywords]) or $S(var [,keywords])

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

@standres <- keyvalue($K,"stand*","TF",default:T)
@usehii <- keyvalue($K,"usehii","TF")

@chars <- keyvalue($K,"symbol*","vector")

@var <- if($v == 0){
	1
}else{
	argvalue($1, "argument 1", "positive integer scalar")
}

if ($v > 2) {
	error("at most 2 non-keyword arguments allowed")
} elseif ($v == 2){ #resvsindex(i,c) allowed for backward compatibility
	if (!isnull(@chars)){
		error("at most 1 non-keyword argument with 'symbols'")
	}
	@chars <- argvalue($2, "argument 2", "vector")
	if (!ischar(@chars) && !isreal(@chars)){
		error("argument 2 not REAL or CHARACTER")
	}
}elseif (!isnull(@chars) && !isreal(@chars)  && !ischar(@chars)){
	error("value of 'symbols' not REAL or CHARACTER")
}

@msg <- _resvsxxxx($S, @usehii, @standres, @var, @chars, $K)
if (!isnull(@msg)){
	error(@msg, macroname:match("resvs*",@msg,0,exact:F) == 0)
}
delete(@usehii,@standres,@var,@chars,@msg)
%resvsyhat%

resid      MACRO DOLLARS
) Macro to mimic residual command in Multreg when used after regress(),
) anova() or manova().
) Usage:
)  resid()  or  resid(model)
)   model           CHARACTER scalar, a GLM model
) If model includes factors to be treated as variates, use resid(model,T)
) Successive columns of result are y (p cols), Studentized residual (p cols),
) HII, Cook's distance (p cols), t-values (p cols) where p is the number of
) dependent variables.
)) 951227, works correctly after weighted analyses & nonlinear GLMs
)) 971107 rows and columns of output are labeled.
)) 000220 stripped $$, use argvalue()
)) 010721 mainly cosmetic changes
# $S() or $S(model) or $S(model,T), where model is a linear model
if($v > 2){
	error("usage: $S() or $S(model) or $S(model,T)", macroname:F)
}

@reg <- if($v <= 1){
	F
}else{
	argvalue($02,"argument 2","TF")
}#T forces regress()
if($v >= 1){
	@model <- argvalue($01,"model","string")
	if(@reg){
		regress(@model,silent:T)
	}else{
		manova(@model,silent:T)
	}
	delete(@model)
}
if(!isdefined(RESIDUALS)||!isdefined(SS)||!isdefined(HII)||!isdefined(DF)){
	error("no active model; try $S(model)", macroname:F)
}
@m <- dim(SS)[1]
if (DF[@m] == 0){
	error("no degrees of freedom for error; can't studentize")
}
@r <- if(isdefined(WTDRESIDUALS)){
	WTDRESIDUALS
}else{
	RESIDUALS
}
@w <- getoptions(warnings:T)
setoptions(warnings:F)
@r <-/ sqrt((1-HII)*(diag(SS[@m,,])'/DF[@m]))
@p <- nrows(SS[1,,])
if (@p == 1){
	@labels <- vector("Depvar","StdResids","HII","Cook's D",\
		"t-stats")
}else{
	@J <- run(@p)
	@dlabels <- getlabels(vector(@J,labels:"Depvar "))
	@rslabels <- getlabels(vector(@J,labels:"StdResids "))
	@rclabels <- getlabels(vector(@J,labels:"Cook's D "))
	@rtlabels <- getlabels(vector(@J,labels:"t-stats "))
	@labels <- vector(@dlabels,@rslabels,"HII",@rclabels,@rtlabels)
	delete(@dlabels,@rslabels,@rtlabels,@rclabels,@J)
}
@o <- matrix(vector(\
	modelvars(0),\
	@r,\
	HII,\
	@r^2*HII/((1-HII)*sum(DF[-@m])),\
	@r*sqrt((DF[@m]-1)/(DF[@m]-@r^2))),nrows(@r),\
	labels:structure("(",@labels))
setoptions(warnings:@w)
delete(@w,@r,@m,@p,@labels)
delete(@o,return:T)
%resid%

yhat               MACRO DOLLARS
) Macro to mimic yhat command in Multreg when used after regress(),
) anova() or manova().
) Usage:
)  yhat()  or  yhat(model)
)   model      CHARACTER scalar, a GLM model
) If model includes factors to be treated as variates, use yhat(model,T)
) Successive columns of output are y (p cols), yhat (p cols), predictive
) residuals (p cols), SE[yhat] (p cols), and SE[pred] (p cols), where p is
) the number of dependent variables.
)) 951222 fixed to work properly after weighted analyses
)) 971108 informative labels added to output
)) 000220 stripped $$, use argvalue()
)) 010721 cosmetic changes
# $S() or $S(model) or $S(model,T), where model is a linear model
if($v > 2){
	error("usage: $S() or $S(model) or $S(model,T)", macroname:F)
}
@reg <- if($v <= 1){
	F
}else{
	argvalue($02,"argument 2","TF")
}
if($v >= 1){
	@model <- argvalue($01, "model", "string")
	if(@reg){
		regress(@model,silent:T)
	}else{
		manova(@model,silent:T)
	}
	delete(@model)
}
if(!isdefined(RESIDUALS) || !isdefined(SS) ||\
	!isdefined(HII) || !isdefined(DF)){
	error("apparently no active linear model; try $S(model)", macroname:F)
}
@m <- dim(SS)[1]
if (DF[@m] == 0){
	error("no degrees of freedom for error; can't studentize")
}

@p <- nrows(SS[1,,])
if (@p == 1){
	@labels <- vector("Depvar","Pred","Pred Resid","SE Est",\
		"SE Pred")
}else{
	@ylabels <- @plabels <- @prlabels <- @seelabels <- @seplabels <- NULL
	@J <- run(@p)
	@ylabels <- getlabels(vector(@J,labels:"Depvar "))
	@plabels <- getlabels(vector(@J,labels:"Pred "))
	@prlabels <- getlabels(vector(@J,labels:"Pred Resid "))
	@seelabels <- getlabels(vector(@J,labels:"SE Est "))
	@seplabels <- getlabels(vector(@J,labels:"SE Pred "))
	@labels <- \
		vector(@ylabels,@plabels,@prlabels,@seelabels,@seplabels)
	delete(@ylabels,@plabels,@prlabels,@seelabels,@seplabels,@J)
}
@wts <- if(isdefined(WTDRESIDUALS)){
	modelinfo(weights:T)
}else{
	rep(1,nrows(RESIDUALS))
}
@mse <- diag(SS[@m,,])'/DF[@m]
@depvar <- modelvars(0)
delete(@reg, @m)

@w <- getoptions(warnings:T)
setoptions(warnings:F)
@o <- matrix(vector(\
	@depvar ,\
	@depvar-RESIDUALS,\
	RESIDUALS/(1-HII),\
	sqrt(HII*@mse/@wts),\
	sqrt((1+HII/@wts)*@mse)),nrows(@depvar),\
	labels:structure("(",@labels))
setoptions(warnings:@w)
delete(@w,@wts,@depvar,@mse)
delete(@o, return:T)
%yhat%

anovapred     MACRO DOLLARS
) Macro to compute anova cell means and standard errors.
) usage:
)  anovapred(a,b, ... ) , where a, b, ..., are all the factors in model
)
) The result is structure(estimate:predvals,SEest:seest, SEpred:sepred)
) where the components are vectors, matrices or arrays with a dimension for
) each factor.  If the response is multivariate, each component has an
) extra dimension, with the last dimension indexing variables.
)
) predvals contains estimated cell means, seest contains standard errors
) of the estimated cell means and sepred contains standard errors of
) predicted values for each cell.
)
) NOTE: Values for any empty cells are set to MISSING, even if they are
) mathematically well defined. Use glmpred() in that case.
)
) Results are based on the side effect variables from the most recent
) anova() or manova() command.
)) 000220 stripped $$
)) 000520 recognize case when error DF = 0 and set SE's to MISSING
)) 010514 uses modelinfo(y:T) instead of <<DEPV>>
# $S(a,b,...) , where a, b, ..., are all the factors in the model
if($v == 0){
	error("usage is $S(a,b,...) where a,b,... are all the factors in model")
}
if (isnull(modelinfo(strmodel:T,nomodelok:T))){
	error("apparently no active linear model")
}
if(!isdefined(RESIDUALS) || !isdefined(HII) || !isdefined(DF) || !isdefined(SS)){
	error("$S needs RESIDUALS, HII, DF and SS; not all defined",macroname:F)
}
@error <- dim(SS)[1]
if (DF[@error] == 0){
	@seest <- @sepred <- tabs(,$V,count:T)
	@seest[] <- @sepred[] <- ?
}else{
	@mse <- diag(SS[@error,,])'/DF[@error]
	@w <- getoptions(warnings:T)
	setoptions(warnings:F)
	@seest <- sqrt(tabs(@mse*HII,$V,mean:T))
	@sepred <- sqrt(tabs(@mse*(1+HII),$V,mean:T))
	setoptions(warnings:delete(@w,return:T))
	delete(@mse)
}
delete(@error)
@estimate <- tabs(modelinfo(y:T) - RESIDUALS,$V,mean:T)
structure(estimate:delete(@estimate,return:T),\
	SEest:delete(@seest,return:T), SEpred:delete(@sepred,return:T))
%anovapred%

regs               MACRO  DOLLARS
) regs(x,y [,GLM keywords]), REAL matrix x, REAL matrix or vector y.
) regs(x,y, T [,GLM keywords])
) regression of y on columns of x (manova if y has more than 1 column)
) When argument 3 is T, no constant is fit
)) Version of 990901 uses argvalue(), can use GLM keywords
)) 990901 fixed some bugs; added ability to use GLM keywords
)) 991012 modified so that '$$' could be removed and DOLLARS pu
))        on the header; it now uses argvalue() to check arguments
)) 000217 added capability of fitting model without intercept
)) 000218 stripped $$
#  $S(x,y [,T] [,GLM keyword phrases]),matrix or vector y, matrix x
@xname <- "@X"
@yname <- "@Y"
@Xvars <- matrix(argvalue($1,"argument 1","real matrix nonmissing"))
<<@yname>> <- matrix(argvalue($2,"argument 2","real matrix nonmissing"))
@nocon <- if($v > 2){
	argvalue($3,"argument 3","TF")
}else{
	F
}
if (nrows(<<@yname>>) != nrows(@Xvars)){
	error("number of rows of x and y are different")
}
@p <- ncols(@Xvars)

<<paste(@xname,1,sep:"")>> <- @Xvars[,1]
@model <- paste(@yname,"=",@xname,1,sep:"")
if(@p > 1){
	for(@i,run(2,@p)){
		<<paste(@xname,@i,sep:"")>> <- @Xvars[,@i]
		@model <- paste(@model,"+",@xname,@i,sep:"")
	}
	delete(@i)
}
if (delete(@nocon,return:T)){
	@model <- paste(@model,"-1",sep:"")
}
if(isvector(<<@yname>>)){
	if ($k > 0){
		regress(@model,$K)
	}else{
		regress(@model)
	}
}else{
	if ($k > 0){
		manova(@model,$K)
	}else{
		manova(@model)
	}
	if (getoptions(warnings:T)){
		print("NOTE: use secoefs() to obtain coefficients and standard errors.")
	}
}
delete(@xname,@yname,@p,@model,@Xvars)
%regs%

===> nlreg <===
nlreg         MACRO DOLLARS
) Nonlinear least squares regression
) Usage:
)  nlreg(b,x,y [,func],param [,resid:resmac,crit:vec,active:active,\
)        maxit:itmax,minit:itmin,print:T, keep:T, quiet:T])
) b      REAL vector of starting values for the iterative fitting
) x      REAL matrix
) y      REAL vector
) func   macro called as fit <- func(x,b,param) (not allowed with resid)
) param  a vector or structure of additional parameters for func or NULL
) vec    vector(numsig, nsigsq, delta), 3 criteria for convergence
)        numsig = number of digits of accuracy in coefficients
)        nsiqsq = number of digits of accuracy in residual SS
)        delta  = norm of gradient threshhold
) active LOGICAL vector the same length as b
) resmac a macro called as resmac(b,x,y [,param]) to compute residuals
)        from a function defined or directly referenced in resmac.
)        func should be omitted when resid:resmac is an argument.
)        Conversely, when resid:resmac is not an argument, argument
)        func is required
) itmin  the minimum number >= 0 of iterations performed
) itmax  the maximum number >= itmin of iterations allowed
) print  If T, partial results printed at each iteration
) keep   If T, nlreg returns the structure returned by levmar()
)        with components, coefs, hessian, jacobian, gradient,
)        rss, residuals, nobs, iter, iconv plus component edf
) quiet  If F (default, unless keep:T), no summary results are
)        printed.  quiet:T is illegal without keep:T
)) Other macros used (should be loaded automatically if available)
))   levmar()
))    _cgrad()
))    _lmout() (only with print:T)
))
) Written by C. Bingham, December 1998
) Version 000909
)) 011129 fixed bug with keep:T (@edf wasn't being created)
)) 030401 modified usage comments
# $S(b,x,y,f,param [,deriv:deriv,crit:vec,active:active,maxit:itmax,\
#       minit:itmin,print:T, keep:T, quiet:T])
# or
# $S(b,x,y, param ,resid:res [,deriv:deriv, crit:vec,\
#       active:active,maxit:itmax,minit:itmin,print:T, keep:T, quiet:T])
# b      REAL vector of starting values for the iterative fitting
# x      REAL matrix with nrows(x) = nrows(y) when f is an argument 
# y      REAL vector
# f      macro called as fit <- f(b,x,param); omitted if resid:res
#        is an argument
# res    macro called as res(b,x,y [,param]) to compute residuals
#        from a function defined or directly referenced in res();
#        when resid:res is an argument, f must not be an argument
# param  a vector or structure of additional parameters for f or NULL
# vec    vector(numsig, nsigsq, delta), 3 criteria for convergence
#        numsig = number of digits of accuracy in coefficients
#        nsiqsq = number of digits of accuracy in residual SS
#        delta  = norm of gradient
# active LOGICAL vector the same length as b
# itmin  the minimum number >= 0 of iterations performed
# itmax  the maximum number >= itmin of iterations allowed
# print  If T, partial results printed at each iteration
# keep   If T, $S returns the structure returned by levmar()
#        with components, coefs, hessian, jacobian, gradient,
#        rss, residuals, nobs, iter, iconv
# quiet  If F (default, unless keep:T), summary results are
#        printed.  quiet:T is illegal without keep:T
@keep <- keyvalue($K, "keep", "TF", default:F)

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

@maxit <- keyvalue($K,"maxit*","count",default:30)

@crit <- keyvalue($K,"crit*", "nonmissing real vector",\
	default:vector(5, 8, -1))
@crit <- padto(@crit,3)
@active <- keyvalue($K,"active", "nonmissing logical vector")

if (@quiet && !@keep){
	error("quiet:T illegal without keep:T on $S", macroname:F)
}

if(!ismacro(levmar)){
	getmacros(levmar,silent:T)
}
@result <- levmar($0)
@coefs <- @result$coefs
if (isnull(@active)){
	@active <- (@coefs == @coefs)
}
COEF <- @coefs[@active]
XTXINV <- solve(@result$hessian)
RESIDUALS <- @result$residuals

@nobs <- @result$nobs

@extra <- nrows(RESIDUALS) - @nobs
@jacobian <- if (@extra > 0){
	RESIDUALS <- RESIDUALS[-run(@extra)]
	@result$jacobian[-run(@extra),]'
}else{
	@result$jacobian'
}

HII <- vector(sum(@jacobian * (XTXINV %*% @jacobian)))
delete(@jacobian, @extra)

@edf <- @nobs - length(COEF)
if (!@quiet){
	@iter <- @result$iter
	@iconv <- @result$iconv
	@npar <- length(@coefs)
	@mse <- sum(RESIDUALS^2)/@edf

	@fmt <- getoptions(format:T)
	@width <- floor(vecread(string:@fmt,silent:T))

	@labw <- if (@npar < 10){
		3
	}else{
		4
	}
	print(paste(charwidth:@labw," ",charwidth:@width,\
		"Coef","StdErr","t","P Value", justify:"r"))

	@se <- sqrt(@mse*diag(XTXINV))
	@tstat <- COEF/@se
	@k <- 0
	for (@i,1,@npar){
		@line <- paste(sep:"",charwidth:@labw-2,"B",format:"2.0f",@i,\
			sep:" ",format:@fmt,@coefs[@i])
		if (@active[@i]){
			@k <-+ 1
			print(paste(sep:" ",format:@fmt,@line,@se[@k], @tstat[@k],\
			2*(1-cumstu(abs(@tstat[@k]),@edf))))
		}elseif(@coefs[@i] != 0){
			print(@line)
		}
	}
	print(paste(rep("-",@labw + 4 + 4*@width),sep:""))
	print(paste("N: ",@nobs, ", MSE: ", @mse, ", DF: ", @edf,sep:""))

	if(@iconv <= 0){
		print(paste("Did not converge in",@maxit,"iterations"))
	}elseif(@iconv <= 3){
		@msg <- if (@iconv == 1){
			paste("relative change in all coefs <", 10^-@crit[1])
		}elseif (@iconv == 2){
			paste("relative change in RSS <", 10^-@crit[2])
		}else{
			paste("norm of gradient =",sqrt(sum(@result$gradient^2)),"<",\
				@crit[3])
		}
		print(paste("Converged with",@msg,"in", @iter,"iterations"))
	}else{
		print(paste("Halving step did not reduce RSS on",@iter,"iteration"))
	}
	delete(@fmt,@se,@tstat,@i,@iter,@npar,@iconv, @width)
}
delete(@nobs, @quiet, @maxit)
if (delete(@keep,return:T)){
	strconcat(delete(@result,return:T),edf:delete(@edf,return:T))
}else{
	delete(@result)
}
%nlreg%

===> linear <===
linear        MACRO
) Sample function linear(b,x) for use with nlreg and levmar()
) It computes x %*% b which is linear in b
# fit <- linear(b,x)
($2) %*% ($1)
%linear%

===> asymptot <===
asymptot      MACRO DOLLARS
) Sample function asymptot(b,x [,j]) for use with nlreg and levmar()
) It computes b[1] + b[2]*b[3]^x
) asymptot(b,x,j) computes derivative with respect to b[j]
) An example is in Sec. 19.8, p. 409ff of Snedecor and Cochran
) Statistical Methods, 7th Edition
) Test data from Snedecor and Cochran are in data set SandCT19.8.1
) version of 990111
# fit <- asymptot(b,x,[y,param,j])
@b <- $1
@tmp <- if ($v < 5){
	@b[1] + @b[2]*@b[3]^($2)
}else{
	@k <- $5
	if (@k == 1){
		rep(1,nrows($2)) # derivative w.r.t. b[1]
	}elseif (@k == 2){
		@b[3]^($2) # derivative w.r.t. b[2]
	}else{
		@x <- $2
		@tmp1 <- @b[2]*@x*@b[3]^(@x-1) # derivative w.r.t. b[3]
		delete(@x)
		delete(@tmp1,return:T)
	}
}
delete(@b)
delete(@tmp,return:T)
%asymptot%

===> DandSFunc <===
DandSFunc  MACRO DOLLARS
) Function to compute nonlinear function and derivatives
) with respect to its parameters defined in equation
) 10.3.1 on p. 475 in Draper and Smith,
) Applied Regression Analysis, 2nd Edition, Wiley 1981
) Function is f(b,x) = b[1] + (.49 - b[1])*exp(-b[2]*(x - 8))
) Test are in data set DandS_T10.2, with suggested starting
) values vector(.30, .02)
@b <- $1
@tmp <- if ($v < 5){
	@b[1] + (.49 - @b[1])*exp(-@b[2]*(($2) - 8))
}else{
	@k <- $5
	if (@k == 1){#df(b,x)/d b[1]
		1 - exp(-@b[2]*(($2) - 8))
	}else{#df(b,x)/d b[2
		-(.49 - @b[1])*(($2) - 8)*exp(-@b[2]*(($2) - 8))
	}
}
delete(@b)
delete(@tmp,return:T)
%DandSFunc%

===> testfun <===
testfun    MACRO DOLLARS
) Sample function testfun(b,x) for use with nlreg() and levmar()
) It computes b[1]+b[2]*x+b[3]*x^2+b[4]*sin(b[5]*x,radian:T)+b[6]*exp(-b[7]*x)
# fit <- testfun(b,x)
@tmp <- getoptions(angles:T)
setoptions(angles:"radians")
@x <- $2
@b <- $1
@f <- @b[1] + @b[2]*@x + @b[3]*@x^2 +\
	@b[4]*sin(@b[5]*@x) + @b[6]*exp(-@b[7]*@x)
setoptions(angles:delete(@tmp,return:T))
delete(@b,@x)
delete(@f,return:T)
%testfun%

===> SandCT19.8.1 <===
SandCT19.8.1     6     2 ENDED
) Data from Table 19.8.1, p. 411 in Snedecor and Cochran
) Statistical Methods, 7th Edition, Iowa State University
) Press (1980)
) For testing nonlinear fitting of function in Snedecor &
) Cochran eq. 19.8.1, defined in macro asymptot()
) Suggested starting values are vector(30,25,.55)
) Col. 1: x
) Col. 2: y
)"%lf %lf"
  0.0 57.5
  1.0 45.7
  2.0 38.7
  3.0 35.3
  4.0 33.1
  5.0 32.2
%SandCT19.8.1%

===> DandS_T10.2 <===
DandS_T10.2    44     2 ENDED
) Data from Table 10.2, p. 476 in Draper & Smith
) Applied Regression Analysis, 2nd Edition, Wiley 1981
) For use with function
) Sample data for use with nlreg() with macro DandSFunc()
) Suggested starting values are b0 <- vector(.30, .02)
) Col. 1: x
) Col. 2: y
)"%lf %lf"
  8  0.49
  8  0.49
 10  0.48
 10  0.47
 10  0.48
 10  0.47
 12  0.46
 12  0.46
 12  0.45
 12  0.43
 14  0.45
 14  0.43
 14  0.43
 16  0.44
 16  0.43
 16  0.43
 18  0.46
 18  0.45
 20  0.42
 20  0.42
 20  0.43
 22  0.41
 22  0.41
 22  0.40
 24  0.42
 24  0.40
 24  0.40
 26  0.41
 26  0.40
 26  0.41
 28  0.41
 28  0.40
 30  0.40
 30  0.40
 30  0.38
 32  0.41
 32  0.40
 34  0.40
 36  0.41
 36  0.38
 38  0.40
 38  0.40
 40  0.39
 42  0.39
%DandS_T10.2%

stepsetup        MACRO  DOLLARS
) Setup stepwise regression
) Usage:
)   stepsetup(Model [,allin:T or in:logvec] [,silent:T]), where
)     Model is a regression model similar to
)        "y=x1+x2+...+xk" or "y=x1+x2+...+xk-1"
)   stepsetup([allin:T] [,silent:T]), with no model, is equivalent
)   to stepsetup(STRMODEL [,allin:T] [,silent:T])
)
) stepsetup initializes an invisible structure _STEPSTATUS and,
) except with silent:T, prints F-to-enter with P-values
)
) With allin:T, stepsetup starts with all variables in.
) With in:logvec, logvec must be a LOGICAL vector of length k with
) logvec[j] True if and only if xj is in the initial model.  allin:T
) and in:rep(T,k) have the same effect.
)
) Adapted from a macro by Gary Oehlert
)) 010515 uses varnames() instead of TERMNAMES and DEPVNAME to get variable
))        names
)) Version 010515
#$S(model [,silent:T] [,in:logVector  or  allin:T])
@model <- if ($v > 0){
	$01
}else{
	NULL
}
@model <- if (!isnull(@model)){
	argvalue(@model,"model","string")
}elseif(isscalar(STRMODEL,char:T)){
	STRMODEL
}else{
	error("No model provided and STRMODEL does not exist")
}
if(match("*=*",@model, 0, exact:F) == 0){
	error(paste("Illegal model \"",@model,"\"",sep:""))
}
@allin <- keyvalue($K,"allin","TF",default:F)
@instart <- keyvalue($K,"in","logical vector nonmissing")
if (@allin && !isnull(@instart)){
	error("keyword 'in' illegal with 'allin:T'")
}
@nocon <- match("*-1",@model,0, exact:F) != 0 || \
	match("*- 1",@model,0, exact:F) != 0
if (!ismacro(stepstatus)){
	getmacros(stepstatus,silent:T)
}
regress(@model,silent:T)
@varnames <- varnames()
@depv <- length(@varnames)
@varnames <- @varnames[rotate(run(@depv),-1)]
@nvars <- @depv - 1
if (!isnull(@instart)){
	if(length(@instart) != @nvars){
		error("Length value of 'instart' different from number of variables")
	}
	@allin <- (sum(@instart) == @nvars)
}else{
	@instart <- rep(F,@nvars)
}
@Fs <- matrix(rep(0, 2*@nvars), 2, labels:structure(vector("F", "P"),\
	@varnames[-@depv]))
if (!@allin){
	@X <- modelvars(x:T)
	@y <- modelvars(y:T)
	@sscp <- if(!@nocon){
		bcprd(@X,@y)[-1,-1]
	}else{
		hconcat(@X, @y) %c% hconcat(@X, @y)
	}
	@dfe <- nrows(RESIDUALS) - 1*(!@nocon)
	if (sum(@instart) > 0){
		@sscp <- swp(@sscp, run(@nvars)[@instart])
		@dfe <-- sum(@instart)
	}
	delete(@X,@y)
}else{
	@sscp <- dmat(@depv, 0)
	@sscp[-@depv,-@depv] <- if(@nocon){
		XTXINV
	}else{
		XTXINV[-1,-1]
	}
	@sscp[@depv, -@depv] <- if(@nocon){
		COEF'
	}else{
		COEF[-1]'
	}
	@sscp[-@depv, @depv] <- if(@nocon){
		-COEF
	}else{
		-COEF[-1]
	}
	@sscp[@depv, @depv] <- reverse(SS)[1]
	@dfe <- reverse(DF)[1]
	@instart <- rep(T,@nvars)
}

@history <- if(sum(@instart) > 0){
	run(@nvars)[@instart]
}else{
	NULL
}

 # build model for _STEPSTATUS$model
@nin <- sum(@instart)
if (@nin == 0){
	@model <- "y=1"
}else{
	@model <-\
		paste(@varnames[@depv],"=",@varnames[-@depv][@instart][1],sep:"")
	if (@nin > 1){
		for(@i,2,@nin){
			@model <- paste(@model,@varnames[-@depv][@instart][@i],sep:"+")
		}
	}
}
if (@nocon){
	@model <- paste(@model,"-1",sep:"")
}

setlabels(@sscp, structure(@varnames,@varnames))
_STEPSTATUS <- structure(model:delete(@model,return:T),\
	sscp:delete(@sscp,return:T),\
	in:vector(delete(@instart,return:T),labels:@varnames[-@depv]),\
	F:delete(@Fs,return:T), dfe:delete(@dfe,return:T),\
	fullmse:reverse(SS)[1]/reverse(DF)[1],history:delete(@history,return:T) )
delete(@nocon,@nvars,@varnames,@allin,@depv,@nin)
stepstatus($K)
%stepsetup%

steplook  MACRO DOLLARS
) Macro to acess elements of _STEPSTATUS
) Usage:
)  steplook(name1, name2, ...)   names  model, sscp, in, F, dfe or
)                                fullmse, history
) If there is one argument, steplook returns a scalar, vector or matrix
) Otherwise, step look returns a structure.
)) Version 991010
# $(name1, name2, ...) names model, sscp, in, F, dfe, history
if ($k > 0){
	error("no keywords allowed as arguments")
}
@compnames <- compnames(_STEPSTATUS)
@names <- $A
for(@i,run($v)){
	if (match("\"*\"",@names[@i], 0,exact:F) != 0){
		@names[@i] <- <<@names[@i]>>
	}
}

@result <- split(rep(0,$v)',compnames:@names)
for(@i,run($v)){
	@k <- match(@names[@i],@compnames,0)
	if(@k == 0){
		error(paste(@names[@i],"not a legal name"))
	}
	@result[@i] <- _STEPSTATUS[@k]
}
delete(@compnames,@names,@i,@k)
if (ncomps(@result) == 1){
	@result <- @result[1]
}
delete(@result,return:T)
%steplook%

stepstatus       MACRO DOLLARS
) Macro to report status of current stage of stepwise regression
) Usage:
)  stepstatus([silent:T])
) Adapted from a macro by Gary Oehlert
)) Version 991011
)) 020806 use 'upper:T' on cumF()
# $S([silent:T])
@silent <- keyvalue($K,"silent","TF",default:F)
@sscp <- _STEPSTATUS$sscp
@depv <- nrows(@sscp)
@nvars <- @depv - 1
@Fs <- _STEPSTATUS$F
@thisSSE <- @sscp[@depv, @depv]
@in <- _STEPSTATUS$in
@ss <- @sscp[@depv,-@depv]'*@sscp[-@depv,@depv]/(diag(@sscp)[-@depv])
@dfe <- _STEPSTATUS$dfe
for(@i, run(@nvars)) {
	if(@in[@i]) {#swp just put variable out
		@Fs[1, @i] <- -@ss[@i]/(@thisSSE/@dfe)
		@Fs[2, @i] <- cumF(@Fs[1, @i], 1, @dfe, upper:T)
	} else {
		@Fs[1, @i] <- @ss[@i]/((@thisSSE - @ss[@i])/(@dfe-1))
		@Fs[2, @i] <- cumF(@Fs[1, @i], 1, @dfe - 1, upper:T)
	}
}
@ii <- match("F",compnames(_STEPSTATUS))
_STEPSTATUS[@ii] <- delete(@Fs,return:T)

@model <- _STEPSTATUS$model
@varnames <- getlabels(@sscp,1)
@nocon <- match("*-1",@model,0, exact:F) != 0
@nin <- sum(@in)
if (@nin == 0){
	@model <- paste(@varnames[@depv],"=1",sep:"")
}else{
	@model <- \
		paste(@varnames[@depv],"=",@varnames[-@depv][@in][1],sep:"")
	if (@nin > 1){
		for(@i,2,@nin){
			@model <- paste(@model,@varnames[-@depv][@in][@i],sep:"+")
		}
	}
}
if (@nocon){
	@model <- paste(@model,"-1",sep:"")
}
@ii <- match("model",compnames(_STEPSTATUS))
_STEPSTATUS[@ii] <- delete(@model,return:T)

delete(@i,@ii,@thisSSE,@ss)

if (!delete(@silent,return:T)){
	print(paste("Current model: \"",_STEPSTATUS$model,"\"",sep:""))
	if (@nin > 0){
		@p <- @nin + 1*(!@nocon)
		@n <- @dfe + @p
		@in <- vector(@in,F)
		@ssreg <- -@sscp[@depv,@in] %*% solve(@sscp[@in,@in], @sscp[@in,@depv])
		@sse <- @sscp[@depv,@depv]
		@rsq <- @ssreg/(@ssreg + @sse)
		@adjrsq <- 1 - (@n - 1*(!@nocon))*(1 - @rsq)/@dfe
		@fstat <- (@ssreg/@nin)/(@sse/@dfe)
		@cp <- @sse/_STEPSTATUS$fullmse + 2*@p - @n
		print(paste("Model F(",@nin,",",@dfe,") = ",@fstat,", P = ",\
			cumF(@fstat,@nin,@dfe, upper:T),", MSE = ", @sse/@dfe,sep:""))
		print(paste("p = ",@p,", C(p) = ",@cp,", Adj R^2 = ",@adjrsq,\
			", R^2 = ",@rsq,sep:""))
		delete(@p,@n,@ssreg,@sse,@rsq,@adjrsq,@fstat,@cp)
	}
	print("")
	if(sum(_STEPSTATUS$in) > 0) {
		print(paste("F(1, ",@dfe,") to delete",sep:""))
		print(_STEPSTATUS$F[, _STEPSTATUS$in], header:F)
	}else{
		print("No variables are \"in\"")
	}
	print("")
	if(sum(_STEPSTATUS$in) < @nvars) {
		print(paste("F(1, ",@dfe - 1,") to enter",sep:""))
		print(_STEPSTATUS$F[, !_STEPSTATUS$in], header:F)
	}else{
		print("All variables are \"in\"")
	}
}
delete(@dfe,@in,@nin)
_STEPSTATUS # invisible so it won't print but can be assigned
%stepstatus%

entervar     MACRO DOLLARS
) Macro to take one or more steps forward, entering a variable or
) several variables into the current stepwise regression model
) Usage:
)   entervar(var1 [,var2 ...] [,silent:T]), var1, var2 ... names
)       or numbers of variables to be entered
) Adapted from a macro by Gary Oehlert
)) 010514 Bug fix so it can work with models with temporary variables
)) 010519 fixed bug in bug fix
)) Version 010519
# $S(var1 [,var2 ...] [,silent:T]), varj a variable name or number
if ($v == 0 || $k > 1){
	error("usage is $S(indepvar [,silent:T])")
}
if (!isstruc(_STEPSTATUS)){
	error("You must run stepsetup(Model) before using $S")
}
@varnames <- getlabels(_STEPSTATUS$in)
@vars <- $A
@nv <- length(@vars)
@silent <- keyvalue($K,"silent","TF",default:T)
for (@j, 1, @nv){
	@var <- @vars[@j]
	if(match("*:*",@var,0,exact:F) == 0){
		if (match("\"*\"",@var,0,exact:F) != 0){
			@var <- <<@var>> #strip off quotes if present
		}

		if (alltrue(isdefined(<<@var>>),\
					isscalar(@newvar <- <<@var>>, real:T))){
			if (@newvar != floor(@newvar) || @newvar <= 0 ||\
				@newvar > length(@varnames)){
				error(paste(@var,"not positive integer <= nvars"))
			}
			@var <- @varnames[@newvar]
		}else{
			@newvar <- match(@var, @varnames, 0)
			if(@newvar == 0) {
				error(paste("Unrecognized variable \"",@var,"\"",sep:""))
			}
		}
		if (_STEPSTATUS$in[@newvar]){
			error(paste("Variable \"",@var,"\" already in model",sep:""))
		}
		delete(@var)

		@ii <- match("sscp",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- swp(_STEPSTATUS$sscp, @newvar)
		@in <- _STEPSTATUS$in
		@in[@newvar] <- T
		@ii <- match("in",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- @in
		@ii <- match("dfe",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- _STEPSTATUS$dfe - 1
		@ii <- match("history",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- vector(_STEPSTATUS$history,@newvar)

		delete(@ii,@newvar)
		stepstatus($K)
		if (!@silent && @j < @nv){
			print("")
		}
	}
}
delete(@vars, @j, @nv, @varnames, @silent)
_STEPSTATUS
%entervar%

removevar    MACRO  DOLLARS
) Macro to take one or more steps backward, removing a variable or
) several variables from the current stepwise regression model
) Usage:
)   removevar(var1 [,var2 ...] [,silent:T]), var1, var2 ... names
)       or numbers of variables to be removed
) Adapted from a macro by Gary Oehlert
)) 010514 fixed bug so it works with models with temporary predictors
)) 010519 fixed bug in bug fix
)) Version 010519
# $S(var1 [,var2 ...] [,silent:T]), varj a variable name or number
if ($v == 0 || $k > 1){
	error("usage is $S(var1 [,var2 ...] [,silent:T])")
}
if (!isstruc(_STEPSTATUS)){
	error("You must run stepsetup(Model) before using $S")
}
@silent <- keyvalue($K, "silent", "TF",default:F)
@varnames <- getlabels(_STEPSTATUS$in)
@vars <- $A
@nv <- length(@vars)
for (@j, 1, @nv){
	@var <- @vars[@j]
	if(match("*:*",@var,0,exact:F) == 0){
		@outvar <- match("$1", @varnames, 0)
		if (match("\"*\"",@var,0,exact:F) != 0){
			@var <- <<@var>> #strip off quotes if present
		}
		if (alltrue(isdefined(<<@var>>),\
					isscalar(@outvar <- <<@var>>, real:T))){
			# it's a number
			if (@outvar != floor(@outvar) || @outvar <= 0 ||\
				@outvar > length(@varnames)){
				error(paste(@var,"not positive integer <= nvars"))
			}
			@var <- @varnames[@outvar]
		}else{
			@outvar <- match(@var, @varnames, 0)
			if(@outvar == 0) {
				error(paste("Unrecognized variable \"",@var,"\"",sep:""))
			}
		}
		if (!_STEPSTATUS$in[@outvar]){
			error(paste("Variable \"",@var,"\" not currently in model",sep:""))
		}
		delete(@var)

		@ii <- match("sscp",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- swp(_STEPSTATUS$sscp, @outvar)
		@in <- _STEPSTATUS$in
		@in[@outvar] <- F
		@ii <- match("in",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- @in
		@ii <- match("dfe",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- _STEPSTATUS$dfe + 1
		@ii <- match("history",compnames(_STEPSTATUS))
		_STEPSTATUS[@ii] <- vector(_STEPSTATUS$history,-@outvar)
		stepstatus($K)
		if (!@silent && @j < @nv){
			print("")
		}
	}
}
delete(@ii,@outvar, @varnames, @vars, @j, @silent)
_STEPSTATUS
%removevar%

regresshelp   MACRO
) Macro to get help on macros in regress.mac
# usage $S(topic1 [, topic2 ...] [help keywords])
if(!ismacro(_gethelp)){
	getmacros(_gethelp,quiet:T,printname:F)
}
___HELPFILE_ <- "regress.mac"
___INDEXTOP_ <- "regress_index"
___MACRO_ <- "$S"
_gethelp($0)
%regresshelp%

_E_N_D_O_F_M_A_C_R_O_S_

Help file for regress.mac for MacAnova
(C) 2003 by Gary W. Oehlert and Christopher Bingham
Updated 030814 CB

Note: Help entries for resvsyhat(), resvsindex(), resvsrankits(),
resid(), yhat() and regs() no longer duplicate entries in MacAnova.hlp

!!!! Starting marker for message of the day
!!!! Ending marker for message of the day

???? Starting marker for list of up to 32 comma/newline separated keys
ANOVA
Confidence limits
General
GLM
Hypothesis test
Nonlinear fitting
Plotting
Prediction limits
Regression
Residuals
Standard error
Stepwise regression
???? Ending marker for keys

====anovapred()#glm,anova,standard error,confidence limits,prediction limits
%%%%
anovapred(a,b,...), a, b, ... all the factors in STRMODEL
%%%%
@@@@usage#Usage
anovapred(a,b,...), where a, b, ... are all the factors in the most
recent GLM model, computes the fitted (predicted) value, the standard
error of estimation, and the standard error of prediction for each
cell. The result is a structure with components 'estimate', 'SEest' and
'SEpred', each of which is a vector, matrix, or array with dimensions
derived from the number of levels of a, b, ....  It uses side effect
variables DEPVNAME, RESIDUALS, and HII.

When the most recent GLM model was manova() with a p-dimensional
dependent variable, each component will have an extra dimension of size
p.

When the most recent GLM model included variates (non-factors), or when
you do not include all factors in the argument list, the results will
probably be wrong, although no warning message will be printed.

anovapred() is implemented as a pre-defined macro.

@@@@see_also#Cross references
See also predtable(), regpred(), 'glm'.
@@@@______

====betalimits()#regression,confidence limits
%%%%
betalimits(Term, level), Term a CHARACTER scalar, a positive integer or
  a variable in the most recent regression model, REAL positive scalar
  level < 1.
%%%%
@@@@usage#Usage
You use betalimits() to compute confidence limits for a regression
coefficient.  Its use must be preceded by a GLM command, usually
regress() or anova().

betalimits(Term, Level), where Term specifies a variable or a term in
the ANOVA table, and Level is a REAL scalar between 0.5 and 1, returns
upper and lower confidence limits for a regression coefficeint or ANOVA
main effects or interactions.

Term can be a quoted or unquoted predictor variable or term name, or
CONSTANT or a positive integer specifying the number of the term.  If
the term is an iteraction it must be quoted.

Level is the desired confidence coefficient.  If level < .5, a warning
message is given and the confidence coefficient is assumed to be 1 -
level.

After regress(), the result is vector(lower, upper), where lower and
upper are the confidence limits.

After anova(), when the term is a main effect with k levels, the result
is k by 2 matrix hconcat(lower, upper).  When the term is an
interaction of factors with k1, k2, ... levels, the result is
array(vector(upper,lower),k1,k2,...,2).

@@@@example
Example:
After regress("y=x1 + x2 + x3"), the following are all equivalent
  Cmd> betalimits(x2, .95)
  Cmd> betalimits("x2",.95)
  Cmd> betalimits(3,.95)  #counting the constant, x2 is the third term

They all return a vector of length 2.

If a and b are factors with 3 and 4 levels, respectively, after
anova("y = a*b") (equivalent to anova("y = a + b + a.b"))
  Cmd> betalimits("a.b",.95) # quotes are required
and
  Cmd> betalimits(4, .95) # a.b is term 4
are equivalent and return a 3 by 4 by 2 array with a[i,j,] the upper
and lower limits for the i,j interaction effect.

In all these replacing .95 by .05 gives the same result except a
warning message is printed.

@@@@see_also#Cross references
See also coefs(), secoefs(), regress()
@@@@______

====entervar()#stepwise regression,regression
%%%%
entervar(var1 [,var2 ...] [,silent:T]), var1, var2, ... the names or
  numbers of independent variables in the complete stepwise model but
  not in the current stepwise model
%%%%
@@@@usage#Usage
entervar(NewVar) enters independent variable Var to the current stepwise
regression model.  NewVar can either be a quoted name ("z3"), an
unquoted name (z3) or the number of a variable in the complete stepwise
model but not in the current stepwise model. Thus if the full model is
"y=x1+x2+x3+x4+x5", entervar(x2), entervar("x2") and entervar(2) are
equivalent.

It is an error if NewVar is not an independent variable in the complete
stepwise model or if it has already been entered in.

Invisible variable _STEPSTATUS is updated to reflect the changed model.
See topic '_STEPSTATUS'.

@@@@printed_output#Printed output
The F-to-remove statistics with P-values are printed for all the
variables in the model, including NewVar, and F-to-enter statistics
with P-values are printed for any variables in the full model that have
not yet been entered.

In addition, entervar() prints an overall F statistic and its P-value,
Mallow's Cp statistic, adjusted R^2 and R^2.  The F-statistic tests
the null hypothesis that the coefficients of the "in" variables are 0.
Because the "in" variables have been selected because they appear to
contributed importantly to the regression, the P-value should not be
interpreted literally.

@@@@value_returned#Value returned
The value, which can be assigned (stuff <- entervar(x3)) but is not
printed, is a copy of the updated invisible variable _STEPSTATUS.

@@@@entering_several_variables#Entering several variables
entervar(Var1, Var2 ...) does the same except that more than one
variable is entered.  The model and other statistics are printed after
each variable is entered.  The value returned is _STEPSTATUS after all
have been entered.  An example when the full model is "y=x1+x2+x3+x4+x5"
might be entervar(x1,"x2",4).  This would enter x1, x2 and x4 in that
order.

@@@@silent_keyword#Keyword silent
entervar(Var1 [, Var2 ...], silent:T) does the same, except that the
model and summary statistics are not printed.  You can then use
stepstatus() to print the summary statistics for the new current state
of the stepwise regression process.

@@@@see_also#Cross references
See also removevar(), stepsetup(), stepstatus()
@@@@______

====estimlimits()#regression,confidence limits
%%%%
estimlimits(x, confLevel), REAL vector or matrix x with no MISSING
  elements, 0 < confLevel < 1 a REAL scalar
estimlimits(NULL,factorValues, confLevel), REAL vector or matrix
  factorValues with no MISSING elements
estimalimits(x,factorValues, confLevel)
%%%%
Macro estimlimits() computes confidence limits for the expected value
E(y|x) after a regression.  You can also use it to compute limits for
the expectation of y after anova() when the model includes one or more
factors.

Before using estimlimits(), you must have run regress("y=x") or
regress("y=x1 + x2 + ...  xk") or more generally anova(Model) where
Model may contain both one or more predictors (covariates) and/or one
or more factors

When the arguments to estimlimits() specify a single condition under
which limits are wanted for E(y), the result is vector(lower,upper),
where lower and upper are the limits.

When the arguments specify several conditions (several values of x
and/or several factor levels), the result is hconcat(lower,upper),
where lower and upper are vectors with length(lower) = length(upper) =
number of conditions.

@@@@usage_after_regress
                      Usage after regress()
estimlimits(x, confLevel), where REAL variable x is a scalar (simple
linear regression) or a vector (multiple regression) returns
vector(lower, upper), the confidence limits for E(y|x) for the
specified value of of the predictor variable(s).  x can contain no
MISSING elements and confLevel must be a real scalar between 0 and 1.

When there is more than one predictor variable, length(x) = number of
predictors.

When there is only one predictor (simple linear regression), x can be a
vector.  When there is more than one predictor, x can be a matrix with
ncols(x) = number of predictors.  In both cases, the result is
hconcat(lower,upper).

When confLevel < .5, a warning message is printed and the confidence
level is assumed to be 1 - confLevel.

@@@@usage_after_anova
            Usage after anova() with factors in the model
estimlimits(NULL,FactorValues,confLevel) is appropriate after anova()
with a model with no variates (covariates).

When there is just one factor ("y = a"), FactorValues should be a
scalar or vector of permissible factor levels.  When there are nFactors
> 1 factors, FactorValues should be a vector with length(FactorValues) =
nFactors, or a matrix with ncols(FactorValues) = nFactors.

estimlimits(x,FactorValues,confLevel) is appropriate after anova() with
a model containing both factors and (co)variates.  x and FactorValues
are as in the no factor case and no variate case, respectively.

@@@@examples_after_regress#Examples after regress()
After regress("y=x1 + x2 + x3"):
 estimlimits(vector(2,3,4),.95) returns where lower and upper are the
 limits when x1 = 2, x2 = 3 and x3 = 4.

 estimlimits(hconcat(x01, x02, x03),.95), returns hconcat(lower,upper),
 where lower[I] and upper[I] are limits whenx1 = x01[I], x2 = x02[I]
 and x3 = x03[I]. x01, x02, and x03 must all be vectors of the same
 length.

@@@@examples_after_anova#Examples after anova()
After anova("y = a + b")
  estimlimits(NULL,vector(1,2),.95) returns vector(lower,upper), where
  lower and upper are limits when a = 1 and b = 2

 After anova("y = x + a + b")
  estimlimits(vector(3,4),vconcat(vector(1,2)',vector(1,3)'),.95)
  returns hconcat(lower,upper), where lower[1] and upper[1] are limits
  when x = 3, a = 1 and b = 2 and lower[2] and upper[2] are limits
  when x = 4, a= 1 and b = 3.

@@@@see_also#Cross references
See also predlimits(), regpred(), glmpred().
@@@@______

====news*#general
%%%%
help(news) or help(news:yymmdd1) or help(news:vector(yymmdd1,yymmdd2)),
  where yymmdd1 and yymmdd2 are integers like 990416 and 19990416
  (April 16, 1999) or 010416 and 20010416 (April 16, 2001)
%%%%
@@@@2003_march#2003 march
030320 Moved nlreg() and related macros and data sets from arima.mac to
this file.

@@@@2001_july#2001 july
010725 resvsrankits(), resvsindex(), resvsyhat() now differ only in
their names.  They are now merely front ends to a new invisible macro
_resvsxxxx().

010723 resvsindex(), resvsrankits() and resvsyhat() have been modified.

The default plotting symbol is now a drawn asterisk ("\6") instead of a
printed asterisk ("*").  You specify a different plotting symbol using
keyword 'symbols' just as on other plotting commands and macros.  In
addition, resvsindex() and resvsrankits() (but not resvsyhat()) work
after an ARIMA time series fit using macro arima().

New keyword phrase usehii:T or usehii:F affects how residuals are
standardized.

Default labels and titles have been modified.
@@@@______

====nlreg()#nonlinear fitting,regression
%%%%
nlreg(b,x,y,f,param [,deriv:deriv,crit:vec,active:active,maxit:itmax,\
  minit:itmin, print:T, keep:T, quiet:T])
nlreg(b,x,y param ,resid:res [,deriv:deriv,crit:vec,active:active,\
  maxit:itmax,minit:itmin,print:T,keep:T, quiet:T])
 b      REAL vector of starting values for the iterative fitting
 x      REAL matrix; with argument f, nrows(x) = nrows(y); with
        'resid:res' nrows(x) = nrows(y) can be omitted; if x is not
        used by res(), it should be 0
 y      REAL vector
 f      macro: f(b,x,param) computes a vector of fitted values for y
        (not allowed with 'resid:res'; required without
        'resid:res')
 param  a vector or structure of additional parameters for f() or NULL
 res    macro: res(b,x,y [,param]) computes a vector of residuals from
        a function defined or directly referenced in res() (not allowed
        when f is an argument; required when f is not an argument).
 deriv  optional macro: deriv(b,x,y,param,j) computes derivative of
        f(x,b,param) or -res(b,x,y) with respect to b[j], returning a
        vector the same length as y
 active LOGICAL vector the same length as b specifies which parameters
        are included in iteration
 vec    vector(numsig, nsigsq, delta): 3 criteria for convergence
        numsig = number of digits of accuracy in coefficients
        nsiqsq = number of digits of accuracy in residual SS
        delta  = threshhold for norm of gradient
 itmin  the minimum number >= 0 of iterations performed
 itmax  the maximum number >= itmin of iterations allowed
 print  If T, partial results printed at each iteration
 keep   If T, nlreg() returns the structure returned by levmar()
        with components, coefs, hessian, jacobian, gradient,
        rss, residuals, nobs, iter, iconv plus component edf
 quiet  If F (default, unless keep:T), no summary results are printed.
        quiet:T is illegal without keep:T
%%%%
@@@@introduction#Introduction
nlreg() is a macro for carrying out nonlinear least squares regression.
It uses macro levmar() to do the actual minimization.  The function
fitted and/or residuals from the function are specified by macros in the
argument list.  Derivatives may be computed by a macro.  If no macro to
compute derivatives is provided, derivatives are computed numerically by
differencing.

@@@@usage_and_arguments#Usage and arguments
 nlreg(b,x,y,f,param [,deriv:deriv,crit:vec,active:active,maxit:itmax,\
       minit:itmin,print:T, keep:T, quiet:T])
or
 nlreg(b,x,y, param ,resid:res [,deriv:deriv, crit:vec,\
       active:active,maxit:itmax,minit:itmin,print:T, keep:T, quiet:T])

  b      REAL vector of starting values for the iterative fitting
  x      REAL matrix; when f is an argument, nrows(x) = nrows(y) is
         required; when 'resid:res' is an argument, nrows(x) !=
         nrows(y) is allowed; if res() does not use x, x should be 0
  y      REAL vector
  f      macro; f(b,x,param) computes a vector with length nrows(y) of
         fitted values for y (not allowed with 'resid:res')
  param  a vector or structure of additional parameters for func or NULL
  res    macro; res(b,x,y [,param]) computes a vector of residuals from
         a function defined or directly referenced in res().  f must not
         be an argument when resid:res is an argument.  Conversely,
         when resid:res is not an argument, f is a required  argument.
         res() must return a vector of length nrows(y) and need not use
         x.
  deriv  optional macro: deriv(b,x,y,param,j) computes derivative of
         f(x,b,param) or -res(b,x,y,param) with respect to b[j],
         returning a vector the same length as y
  active LOGICAL vector the same length as b.  active[i] = F means b[i]
         is kept constant and does not participate in iteration
  vec    vector(numsig, nsigsq, delta), 3 criteria for convergence
         numsig = number of digits of accuracy in coefficients [8]
         nsiqsq = number of digits of accuracy in residual SS [5]
         delta  = norm of gradient threshhold [-1]
  itmin  the minimum number >= 0 of iterations performed
  itmax  the maximum number >= itmin of iterations allowed
  print  If T, partial results printed at each iteration
  keep   If T, nlreg() returns the structure returned by levmar() with
         components, coefs, hessian, jacobian, gradient, rss, residuals,
         nobs, iter, iconv plus component edf (error degrees of freedom)
  quiet  If F (default, unless keep:T), no summary results are
         printed.  quiet:T is illegal without keep:T

@@@@value_returned#Value returned
With 'keep:T', nlreg() returns structure(coefs:b_hat,hessian:hes,
jacobian:jac,gradient:g,rss:Rssmin,residuals:resids,nobs:n,iter:niter,
iconv:convflag,edf:errordf).  The components are the same as returned by
levmar() with the addition of 'edf'.

The component values are as follows:

  b_hat     REAL vector which minimizes Rss
  hes       jac' %*% jac, an approximation to the Hessian matrix H,
            where H[i,j] = 2nd order partial derivative of Rss/2
            with respect to b_hat[i] and. b_hat[j]
  jac       the nrows(y) by nrows(b) Jacobian matrix; -jac[,j] = vector
            of partial derivatives of the elements of the residual
            vector with respect to b_hat[j].
  g         the nrows(b) REAL gradient vector with g[j] = partial
            derivative of Rss/2 with respect to b_hat[j].  g should be
            close to a vector of zeros at convergence
  Rssmin    the minimized value of Rss
  resids    vector of residuals of length nrows(y)
  n         positive integer = nrows(y) = nrows(resids) = nrows(jac)
  niter     positive integer = the number of iterations
  conflag   Convergence status flag; 0 = not converged, 1 = met relative
            change in b_hat criterion, 2 = met relative change in Rss
            criterion, 3 = met norm of gradient vector criterion, 4 =
            failed to reduce Rss on a step of the iteration.
  edf       Error degrees of freedom = nrows(y) - length(b_hat).

When any parameters are inactive as specified by keyword 'active' (see
below), jac is nrows(y) by p, hes is p by p, g has length p and edf =
nrows(y) - p, where p = number of active parameters.

@@@@keyword_arguments#Keyword phrase arguments
There are several of keywords that can be used to control the iteration
and hold parameters at fixed values.

  Keyword Phrase   Value and Explanation
  crit:crvec       vector(numsig, nsigsq, delta), 3 criteria for
                   convergence (default = vector(8,5,-1).  See below.
  active:act       LOGICAL vector the same length as b.  b[j]
                   "participates" in the iteration only if act[j] is
                   True (default = rep(T,length(b))).  When act[j] is
                   False, b_hat[j] remains at the starting value
  maxit:itmax      Non-negative integer specifying the maximum number of
                   iterations (default = 30).  When itmax = 0, no iterations
                   are done and the quantities returned are computed at
                   the starting value b
  minit:itmin      Non-negative integer < itmax specifies the mininum
                   number of iterations
  print:T          When T, partial results are printed on each iteration

Iteration stops when (i) itmax is exceeded, (ii) any of the three
convergence criteria are satisfied, or (iii) when there has been no
reduction in Rss on an iteration after 10 halvings of the initial step
size.

@@@@use_of_argument_param#Use of argument param
One use for param might be to provide a choice between several functions
to be fit.  For example, the code for macro f might be the following:

  @b <- $1; @x <- $2; @p <- $3
  @val <- @b[1]
  for(@i,1,@p){
      @P <- @b[3*@i+1]
      @val <- @val + @b[3*@i-1]*cos(@x/@P,cycles:T) +\
        @b[3*@i]*sin(@x/@P,cycles:T)
  }
  @val # return

Here is the equivalent code for a macro res():

  @b <- $1; @x <- $2; @y <- $3; @p <- $4
  @res <- @y - @b[1]
  for(@i,1,@p){
      @P <- @b[3*@i+1]
      @res <- @res - @b[3*@i-1]*cos(@x/@P,cycles:T) -\
        @b[3*@i]*sin(@x/@P,cycles:T)
  }
  @res # return

Either might be used in fitting a curve with p cosine components with
unknown periods, amplitudes or phases, with the number of terms
specified by params.

@@@@convergence_criteria#Convergence criteria
There are 3 possible convergence criteria, at least one of which must be
enabled.  nlreg() terminates iteration when at least itmin iterations
have been completed and any convergence criterion is satisfied or when
itmax iterations have been completed.

The criteria are specified by optional argument crit:vec, where vec is
vector(numsig [, nsigsq [, delta]]) with length <= 3.  The default values
for numsig, nsiqsq and delta are 8, 5 and -1, respectively.

A negative value for a criterion means it is not active.

  numsig   Desired number of significant digits in every active
           coefficient.  Specifically, the criterion is satisfied if,
           for every active coefficient b[j], the change d_j satisfies
           abs(d_j) < 10^-numsig*max(.5,abs(b_j)) where b_j is the
           updated value of b[j].  Component 'iconv' of the return value
           is 1 when satisfied.

  nsigsq   Desired number of significant digits in Rss = the residual
           sum of squares.  Specifically, the criterion is satisfied
           when abs(Rss_new - Rss_old) < 10^-nsigsq*max(.5,Rss_new).
           Component 'iconv' of the return value is 2 when satisfied.

  delta    Desired maximum norm ||g|| of the gradient vector g.
           Specifically, the criterion is satisfied when
           sqrt(sum(g^2)) <= delta.  Component 'iconv' of the return
           value is 3 when satisfied.

For numsig and nsigsq, the returned coefficients are the values updated
on that iteration.  For delta, the returned coefficients are the values
found on the previous iteration.

Criteria are checked in the order delta, numsig and nsigsq.

@@@@test_macros_and_data#Test macros and data
Included in this file are the following four macros used in debugging
and testing nlreg():

 linear()    Linear function f <- x %*% b
 asymptot()  f = b[1] + b[2]*(b[3])^x & derivatives
 DandSFunc() f = b[1] + (.49 - b[1])*exp(-b[2]*(x - 8)) & derivatives
 testfun()   f = b[1]+b[2]*x+b[3]*x^2+b[4]*sin(b[5]*x)+b[6]*exp(-b[7]*x)

In addition, the following data sets are included in this file:
 DandS_T10.2  Data from Draper and Smith to be used with DandSFunc()
 SandCT19.8.1 Data from Snedecor and Cochran to be used with
              asymptot()
@@@@example
Example:
Fit the function b1 + b2*b3^x to data from Snedecor and Cochran with
starting values b1 = b2 = 40 and b3 = 1.

The data are also available in data set SandCT19.8.1 and the macro is
available as macro asymptot(), both in this file.

  Cmd> x <- vector(0, 1, 2, 3, 4, 5)

  Cmd> y <- vector(57.5, 45.7, 38.7, 35.3, 33.1, 32.2)

  Cmd> f <- macro("@b <- $1; @b[1]+@b[2]*@b[3]^($2)", dollars:T)

  Cmd> startVal <- vector(40,40,1) # starting values for iteration

  Cmd> nlreg(startVal,x,y,f)
             Coef      StdErr           t     P Value
  B 1      30.724     0.23099      133.01  9.3704e-07
  B 2      26.821      0.2577      104.08  1.9555e-06
  B 3     0.55184    0.008448      65.322  7.9056e-06
  ---------------------------------------------------
  N: 6, MSE: 0.032416, DF: 3
  Converged with relative change in all coefs < 1e-05 in 7 iterations

  Cmd> resfunc <- macro("($3) - f($1,$2)") # computes residuals

  Cmd> nlreg(startVal,x,y,resid:resfunc) # identical
             Coef      StdErr           t     P Value
  B 1      30.724     0.23099      133.01  9.3704e-07
  B 2      26.821      0.2577      104.08  1.9555e-06
  B 3     0.55184    0.008448      65.322  7.9056e-06
  ---------------------------------------------------
  N: 6, MSE: 0.032416, DF: 3
  Converged with relative change in all coefs < 1e-05 in 7 iterations

  Cmd> stuff <- nlreg(startVal,x,y,f,keep:T)

  Cmd> compnames(stuff)
   (1) "coefs"
   (2) "hessian"
   (3) "jacobian"
   (4) "gradient"
   (5) "rss"
   (6) "residuals"
   (7) "nobs"
   (8) "iter"
   (9) "iconv"
  (10) "edf"

  Cmd> stuff[vector(1,2,5,8,9,10)]
  component: coefs
  (1)       30.724       26.821      0.55184
  component: hessian
  (1,1)            6       2.1683       111.39
  (2,1)       2.1683       1.4367       30.242
  (3,1)       111.39       30.242       2675.8
  component: rss
  (1)     0.097248
  component: iter
  (1)            7
  component: iconv
  (1)            1
  component: edf
  (1)            3

  Cmd> # compute approximate standard errors

  Cmd> sqrt((stuff$rss/stuff$edf)*diag(solve(stuff$hessian)))
  (1)      0.23099       0.2577     0.008448

@@@@see_also#Cross reference
See also levmar().
@@@@______

====predlimits()#regression,prediction limits
%%%%
predlimits(x, confLevel), x REAL scalar, vector or matrix with no
  MISSING elements, 0 < confLevel < 1 scalar
%%%%
@@@@usage#Usage
You can use macro predlimits() to compute prediction limits for y =
E(y|x) + epsilon or y = E(y | x1, x2 ...) + epsilon after running
regress("y=x") or regress("y=x1+x2+..+xk").  These are limits on the a
future value of y for specified values of the predictor variable or
variables.

predlimits(x, confLevel), where x is a REAL vector with length(x) =
number of predictors (1 for simple linear regression), returns
vector(lower,upper), where lower and upper are prediction limit with
confidence level confLevel.  Argument confLevel must be a REAL scalar
between 0.5 and 1.

When confLevel < .5, a warning message is printed and 1 - confLevel is
used.

@@@@example_of_single_prediction#Example of single prediction
Example:
  After regress("y=x1 + x2 + x3"), predlimits(vector(2,3,4), .95)
  returns vector(lower,upper) where lower and upper are the limits when
  x1=2, x2=3 and x3=4

@@@@limits_for_several_predictions#Limits for several predictions
You can use predlimits() to get limits for several values at once,
returning hconcat(lower,upper), where lower and upper are vectors of
limits for the various values.

For simple linear regression, x should be a vector containing the
values for which you want prediction limits.

When there are k predictors, x should be a matrix with k columns, and
lower and upper will have length(nrows(x)).

@@@@example_of_several_predictions#Example of several predictions
Example:
  After regress("y = x1 + x2"), predlimits(vconcat(vector(1,2)',
  vector(2,1.5)',vector(3,3.2)'), .95) returns hconcat(lower,upper)
  where, for example, lower[2] and upper [2] are the limits when x1 = 2
  and x2 = 1.5.

@@@@see_also#Cross references
See also estimlimits(), regpred(), glmpred().
@@@@______

====regcoefs()#glm,anova,regression,confidence limits,standard error
%%%%
regcoefs(Model [,pvals:T] [,byvar:F]) or regcoefs([pvals:T] [,byvar:F]),
  where Model is a CHARACTER scalar
%%%%
@@@@usage#Usage
regcoefs(Model) returns a matrix with appropriately labeled rows and
columns of the regression coefficients, their standard errors and
t-statistics from a least squares fit to the regression model specified
by Model.  There can be no factors in Model.  If Model is omitted, the
most recent GLM model is used.

regcoefs(Model,pvals:T) or regcoefs(pvals:T) also computes two-tail P
values corresponding to the t-statistics on the basis of Student's
t-distribution with degrees of freedom from the last element of side
effect variable DF.

Because of the row and column labels, after any GLM command with a model
withut factors, typing regcoefs([pvals:T]) produces a table similar to
that produced by regress().  After non-linear fits such as logistic() or
poisson(), the P-values will not necessarily be appropriate.

@@@@multivariate_response#Multivariate response
If the response variable is multivariate, the result is a structure,
each of whose components is a labeled matrix of coefficients, standard
errors and t-statistics.  regcoefs(Model,byvar:F) or regcoefs(byvar:F)
returns a single labeled matrix, with separate columns for the
coefficients, standard errors, ... for each variable.

@@@@see_also#Cross references
See also topics 'glm', regress(), secoefs().
@@@@______

====regresshelp()#general
%%%%
regresshelp(topic1 [, topic2 ...] [,usage:T])
regresshelp(index:T)
%%%%
regresshelp(topicname) prints help on a topic related to file
regress.mac.  Usually topicname is the name of a macro in the file.

When quoted, topicname may contain "wildcard" characters "*" and "?".
You can also use help keyword 'key'.  See help() for details.

regresshelp(topicname1, topicname2, ...) prints help on more than one
topic.

regresshelp(topicname1 [, topicname2 ...], usage:T) prints just a brief
summary of usage for the each topic.
@@@@______

====regress_index*#
%%%%
Topics in this file:
  anovapred(), betalimits(), entervar(), estimlimits(), nlreg(),
  predlimits(), regcoefs(), regresshelp(), regress_index, regs(),
  removevar(), resid(), resvsindex(), resvsrankits(), resvsyhat(),
  steplook(), stepsetup(), stepstatus(), testbeta(), testestim(), yhat(),
  _STEPSTATUS
%%%%
@@@@list_of_help_entries#List of help entries
Help entries in this file
anovapred()    Macro to compute anova cell means and standard errors.
betalimits()   Macro to compute confidence limits for a regression
               coefficient betaj
entervar()     Macro to update the current stepwise regression model
               to include a specified variable
estimlimits()  Macro to compute confidence limits for E(y | x)
nlreg()        Macro to carry out non linear regression
predlimits()   Macro to compute confidence limits for y for specified x
regcoefs()     Macro to compute labelled matrix of regression
               coefficients, with standard errors, t-statistics,
               P-values
regresshelp()  Macro to get help from this file
regs()         Macro to carry out regression of y on columns of a matrix
removevar()    Macro to delete a specified variable from the current
               stepwise regression model
resid()        Macro to compute various case statistics related to
               residuals
resvsindex()   Macro to plot studentized residuals against case number
resvsrankits() Macro to plot studentized residuals against normal scores
resvsyhat()    Macro to plot studentized residuals against predicted
               values
steplook()     Macro to retrieve various aspects, such as error DF,
               the "in" variables and the regression model, of the
               current stepwise regression status
stepsetup()    Macro to initialize the stepwise regression process
stepstatus()   Macro to print report on current state of stepwise
               regression process
testbeta()     Macro to Compute t-statistic and optionally DF and
               P-value for test of H0: betaj = hypValue
testestim()    Macro to compute t-statistic and optional DF and P-value
               for test of H0: E(y | x) = hypValue
yhat()         Macro to compute various case statistics related to
               fitted values
_STEPSTATUS    "Invisible" Variable containing current state of stepwise
               regression process
@@@@______

====regs()#glm,regression
%%%%
regs(x,y [,T] [GLM keywords]), x and y REAL matrices with the same
  number of rows; T means no intercept
%%%%
@@@@usage#Usage
regs(x,y) computes the regression of y on the columns of x.  For
example, when data is a n by 7 matrix, say, you can compute the
regression of the last column on the first 6 by regs(data[,-7],
data[,7]).

regs(x,y,T) does the same, except the model fit has no constant term
(intercept).

You can use GLM keyword phrases such as 'pvals:T', 'silent:T',
'wts:weights', and 'marginal:T' as additional arguments.

Macro regs() creates temporary variables @X1, @X2, ... and @Y and then
invokes regress() or, if y has more than 1 column, manova().

When y is univariate, immediately follow regs() by anova() to see the
ANOVA table.  When y is multivariate, follow regs() by regcoefs() to see
coefficients and standard errors.

Because the model for regress() or manova() uses temporary variables,
STRMODEL cannot be used as a model for a subsequent regress(), anova(),
or manova() command.  However, these variables can be retrieved, by
modelvars().   For example, when x has 3 columns,

  Cmd> makecols(modelvars(x:T),X1,X2,X2)

creates variables X1, X2 and X3.  Of course, if x still exists,
makecols(x,X1,X2,X3) does the same.

@@@@see_also#Cross references
See also regress(), anova(), modelvars(), regcoefs(), and 'glm_keys'.
@@@@______

====removevar()#stepwise regression,regression
%%%%
removevar(var1 [,var2 ...] [,silent:T]), var1, var2, ... the names or
  numbers of independent variables in the current stepwise model
%%%%
@@@@usage#Usage
removevar(Var) removes independent variable from the current stepwise
regression model.  Var can be either a quoted name ("z3"), an unquoted
name (z3) or the number of a variable in the current stepwise model (an
'in' variable).  Thus if the full model is "y=x1+x2+x3+x4+x5",
removevar(x2), removevar("x2") and removevar(2) are equivalent.

It is an error if the variable is not an independent variable in the
full stepwise model or if it is not in the current model.

Invisible variable _STEPSTATUS is updated to reflect the changed model.
See topic '_STEPSTATUS'.

@@@@printed_output#Printed output
The F-to-remove statistics with P-values are printed for all the
variables in the model, and F-to-enter statistics with P-values are
printed for any variables not in the model, including Var.

In addition, if there are any variables left in the model, removevar()
prints an overall F statistic and its P-value, Mallow's Cp statistic,
adjusted R^2 and R^2.  The F-statistic tests the null hypothesis that
the coefficients of the "in" variables are 0.

@@@@value_returned#Value returned
The value returned is the updated invisible variable _STEPSTATUS.  It
can be assigned (stuff <- removevar(x3)), but is not printed.

@@@@removing_several_variables#Removing several variables
removevar(Var1, Var2 ...) does the same except that more than one
variable is removed.  All variable must be in the current stepwise
model.  The model and other statistics are printed after each variable
is removed.  The value returned is _STEPSTATUS after all have been
removed.  An example when the full model is "y=x1+x2+x3+x4+x5"
might be removevar(x1,"x2",4).  This would remove x1, x2 and x4 in that
order.

@@@@silent_keyword#Keyword silent
removevar(Var1 [,Var2 ...], silent:T) does the same, except that the
model and summary statistics are not printed.

@@@@see_also#Cross references
See also entervar(), stepsetup(), stepstatus()
@@@@______

====resid()#glm,residuals,anova,regression
%%%%
resid() or resid(Model)
%%%%
@@@@usage#Usage
resid(), with no argument, computes a REAL matrix of various quantities
useful in the analysis of residuals.  It uses side effect variables
RESIDUALS, HII, etc. produced by the most recent GLM (generalized linear
or linear model) command such as regress(), anova(), or poisson().

It is an error if any of the needed side effect variables do not exist.

resid(Model) first executes manova(Model, silent:T) to compute the
required side effect variables before computing the residual-related
quantities.  Model should be a CHARACTER variable or string specifying a
linear ANOVA or MANOVA model.  Any factors in the model will be treated
as factors.  If you want them treated as variates, use resid(Model,T).

@@@@description_of_output#Description of output
Each row of the result corresponds to a case.  When the dependent
variable Y is univariate, there are 5 columns, as follows:
   Col. 1   Y = observed response
   Col. 2   Studentized residuals = RESIDUALS/SE(RESIDUALS)
   Col. 3   HII = leverage
   Col. 4   Cook's distance
   Col. 5   t-statistics = externally studentized residuals =
            RESIDUALS/SE*(RESIDUALS), where SE* for each case is
            a standard error based on the model fit excluding that case.

When Y is multivariate of dimension p, there are 4*p + 1 columns -- the
p values for Y, the p standardized residuals, HII, the p Cook's
distances, and the p externally studentized residuals.

If a case has missing values, most entries for that case will be MISSING
and there are no useful numbers.

@@@@after_nonlinear_GLM_commands#After nonlinear GLM commands
After non-linear GLM commands such as poisson() and logistic(), the
results are based on the last stage of the iteratively reweighted least
squares algorithm used to fit the model.  Residuals are standardized by
the error mean square in the linear scale.  They should still be valid
for diagnosing departures from the model.
@@@@______

The output of resid() is modeled on what is printed by the resid()
command in program Multreg.

resid() is implemented as a pre-defined macro.

@@@@see_also#Cross references
See also topics 'glm', yhat(), resvsindex(), resvsrankits(),
resvsyhat().
@@@@______

====resvsindex()#plotting,glm,residuals,anova,regression
%%%%
resvsindex([varNo,] [usehii:T or F] [,standres:F]\
  [,graphics keyword phrases]), 1 <= varNo <= ncols(RESIDUALS)
%%%%
@@@@usage#Usage
You use resvsindex() to plot standardized (default) or non-standardized
residuals vs case numbers.

resvsindex([graphics keyword phrases]) plots standardized residuals
against case number.

resvsindex(usehii:T [,graphics keyword phrases]) does the same using
leverages HII in standardizing.  This is the default after a GLM
command.

resvsindex(usehii:F [,graphics keyword phrases]) does the same without
using leverages HII.  This is the default after arima().

resvsindex(standres:F [,graphics keyword phrases]) does the same without
any standardization.

The residuals are from variable RESIDUALS or WTDRESIDUALS produced by
the most recent GLM (generalized linear or linear model) command such as
regress(), anova(), or poisson(), or from an ARIMA fit computed by macro
arima().

If the most recent command was manova(), only column 1 of the residual
matrix is plotted, but see the next usage for plotting other columns.

@@@@after_manova#After manova()
resvsindex(varNo [, usehii:T or F] [, standres:F] [,graphics keyword
phrases]), where varNo is an integer between 1 and ncols(RESIDUALS),
plots residuals associated with variable varNo against case numbers.
varNo > 1 is legal only when RESIDUALS was computed by manova().

@@@@plotting_symbols#Plotting symbols
The default plotting symbol is the same as for plot(), a drawn asterisk
or star ("\6").  You can change it by including 'symbols:c' as an
argument, where c is a CHARACTER or integer scalar or vector.  c = 0 is
special: it is equivalent to c = "###" and results in points being
labeled with case number.  See chplot(), subtopic 'symbols_used'.

@@@@graphics_keywords#Graphics keywords
You can use all the usual graphics keywords to modify the default plot
characteristics.  These include 'title', 'xlab', 'ylab', 'symbols'
'impulse' and 'lines'.  See topics 'graphs', 'graph_keys',
'graph_border' and 'graph_ticks'.

When you have set option 'dumbplot' to False (see 'options'), the plot
will be a low resolution plot unless 'dumb:F' is an argument.

@@@@what_is_plotted#What is plotted
Without standres:T, the quantities plotted are r[i]/sd[i] where r[i] is
RESIDUALS[i] or WTDRESIDUALS[i] and sd[i] is the estimated standard
deviation.  WTDRESIDUALS[i] is used after regress(), anova(), or
manova() with 'weights:wts' or after nonlinear GLM commands such as
logistic() and poisson().

When usehii is True (the default after GLM commands), sd[ii] =
sqrt(mse*(1-HII[i])), where mse is the residual mean square after
regress(), anova() or manova(), the mean error deviance after non-linear
GLM commands or the estimated innovation variance after arima().

When usehii is False (the default after arima()), sd[i] = sqrt(mse).

With standres:F, the quantities plotted are r[i].

The values on the X-axis are 1, 2, ..., nrows(RESIDUALS).

@@@@______

resvsindex() is implemented as macro.

@@@@see_also#Cross references
See also topics resvsrankits(), resvsyhat(), resid().
@@@@______

====resvsrankits()#plotting,glm,residuals,anova,regression
%%%%
resvsrankits([varNo,] [usehii:T or F] [,standres:F]\
  [,graphics keyword phrases]), 1 <= varNo <= ncols(RESIDUALS)
%%%%
@@@@usage#Usage
resvsrankits([graphics keyword phrases]) plots standardized residuals
against normal scores as computed by function rankits().

resvsrankits(usehii:T [,graphics keyword phrases]) does the same using
leverages HII in standardizing.  This is the default after a GLM
command.

resvsrankits(usehii:F [,graphics keyword phrases]) does the same without
using leverages HII.  This is the default after arima().

resvsrankits(standres:F [,graphics keyword phrases]) does the same
without any standardization.

The residuals are from variable RESIDUALS or WTDRESIDUALS produced by
the most recent GLM (generalized linear or linear model) command such as
regress(), anova(), or poisson(), or from an ARIMA fit computed by macro
arima().

If the most recent command was manova(), only column 1 of the residual
matrix is plotted, but see the next usage for plotting other columns.

@@@@after_manova#After manova()
resvsrankits(varNo [, usehii:T or F] [, standres:F] [,graphics keyword
phrases]), where varNo is an integer between 1 and ncols(RESIDUALS),
plots residuals associated with variable varNo against case numbers.
varNo > 1 is legal only when RESIDUALS was computed by manova().

@@@@plotting_symbols#Plotting symbols
The default plotting symbol is the same as for plot(), a drawn asterisk
or star ("\6").  You can change it by including 'symbols:c' as an
argument, where c is a CHARACTER or integer scalar or vector.  c = 0 is
special: it is equivalent to c = "###" and results in points being
labeled with case number.  See chplot(), subtopic 'symbols_used'.

@@@@graphics_keywords#Graphics keywords
You can use all the usual graphics keywords to modify the default plot
characteristics.  These include 'title', 'xlab', 'ylab', 'symbols'
'impulse' and 'lines'.  See topics 'graphs', 'graph_keys',
'graph_border' and 'graph_ticks'.

When you have set option 'dumbplot' to False (see 'options'), the plot
will be a low resolution plot unless 'dumb:F' is an argument.

@@@@what_is_plotted#What is plotted
Without standres:T, the quantities plotted are r[i]/sd[i] where r[i] is
RESIDUALS[i] or WTDRESIDUALS[i] and sd[i] is the estimated standard
deviation.  WTDRESIDUALS[i] is used after regress(), anova(), or
manova() with 'weights:wts' or after nonlinear GLM commands such as
logistic() and poisson().

When usehii is True (the default after GLM commands), sd[ii] =
sqrt(mse*(1-HII[i])), where mse is the residual mean square after
regress(), anova() or manova(), the mean error deviance after non-linear
GLM commands or the estimated innovation variance after arima().

When usehii is False (the default after arima()), sd[i] = sqrt(mse).

With standres:F, the quantities plotted are r[i].

The values on the X-axis are normal scores computed by rankits(r).  See
rankits() for more information.

@@@@______

resvsrankits() is implemented as macro.

@@@@see_also#Cross references
See also topics resvsindex(), resvsyhat(), resid().
@@@@______

====resvsyhat()#plotting,glm,residuals,anova,regression
%%%%
resvsyhat([varNo,] [usehii:T or F] [,standres:F]\
  [,graphics keyword phrases]), 1 <= varNo <= ncols(RESIDUALS)
%%%%
@@@@usage#Usage
resvsyhat([graphics keyword phrases]) plots standardized residuals
against fitted or predicted values.

resvsyhat(usehii:F [,graphics keyword phrases]) does the same without
using leverages HII in standardizing.  The default is to use HII.

resvsyhat(standres:F [,graphics keyword phrases]) does the same without
any standardization, that is, the residuals y - y_hat are plotted.

The residuals are from variable RESIDUALS or WTDRESIDUALS produced by
the most recent GLM (generalized linear or linear model) command such as
regress(), anova(), or poisson().

If the most recent command was manova(), only column 1 of the residual
matrix is plotted, but see below for plotting other columns.

Unlike resvsrankits() and resvsindex(), resvsyhat() cannot be used to
make a residual plot after arima() was used to estimate an ARIMA time
series model.

@@@@after_manova#After manova()
resvsyhat(varNo [, usehii:T or F] [, standres:F] [,graphics keyword
phrases]), where varNo is an integer between 1 and ncols(RESIDUALS),
plots residuals associated with variable varNo against case numbers.
varNo > 1 is legal only when RESIDUALS was computed by manova().

@@@@plotting_symbols#Plotting symbols
The default plotting symbol is the same as for plot(), a drawn asterisk
or star ("\6").  You can change it by including 'symbols:c' as an
argument, where c is a CHARACTER or integer scalar or vector.  c = 0 is
special: it is equivalent to c = "###" and results in points being
labeled with case number.  See chplot(), subtopic 'symbols_used'.

@@@@graphics_keywords#Graphics keywords
You can use all the usual graphics keywords to modify the default plot
characteristics.  These include 'title', 'xlab', 'ylab', 'symbols'
'impulse' and 'lines'.  See topics 'graphs', 'graph_keys',
'graph_border' and 'graph_ticks'.

When you have set option 'dumbplot' to False (see 'options'), the plot
will be a low resolution plot unless 'dumb:F' is an argument.

@@@@what_is_plotted#What is plotted
Without standres:T, the quantities plotted are r[i]/sd[i] where r[i] is
RESIDUALS[i] or WTDRESIDUALS[i] and sd[i] is the estimated standard
deviation.  WTDRESIDUALS[i] is used after regress(), anova(), or
manova() with 'weights:wts' or after nonlinear GLM commands such as
logistic() and poisson().

When usehii is True (the default after GLM commands), sd[ii] =
sqrt(mse*(1-HII[i])), where mse is the residual mean square after
regress(), anova() or manova() or the mean error deviance after
non-linear GLM commands.

When usehii is False, sd[i] = sqrt(mse).

With standres:F, the quantities plotted are r[i].

The values on the X-axis are the estimated means of the response
variable.  After a nonlinear GLM command, they are in the original
scale, not the transformed scale.

@@@@______

resvsyhat() is implemented as macro.

@@@@see_also#Cross references
See also topics resvsindex(), resvsrankits(), resid(), yhat().
@@@@______

====steplook()#stepwise regression,regression
%%%%
steplook(item1, item2, ...), item1, item2, ... unquoted words selected
  from 'model', 'sscp', 'in', 'F', 'dfe', 'fullmse' and 'history'
%%%%
@@@@usage#Usage
steplook(item1, item2, ... ) returns the values of components of
hidden variable _STEPSTATUS (see topic '_STEPSTATUS').  The arguments
must match component names of _STEPSTATUS, that is they must be one
or more of 'model', 'sscp', 'in', 'F', 'dfe', 'fullmse' and 'history'.

If only one item is requested, the value is a scalar, vector or
matrix, depending on the item.  Otherwise, the value is a structure.

@@@@examples
Examples:
  Cmd> steplook(model, history, in)
  returns structure(model:_STEPSTATUS$model,history:_STEPSTATUS$history,
    in:_STEPSTATUS$in)

  Cmd> steplook(in)
  returns _STEPSTATUS$in

  Cmd> steplook(F)[1,steplook(in)]
  returns F-to-remove for variables in current stepwise model

  Cmd> steplook(F)[2,!steplook(in)]
  returns P-values for F-to-enter for variables not in current stepwise
  model

  Cmd> regress(steplook(model),pvals:T)

  estimates the coefficients for the current stepwise model.

@@@@see_also#Cross references
See also topics stepsetup(), entervar(), removevar(), stepstatus()
@@@@______

====stepsetup()#stepwise regression,regression
%%%%
stepsetup([Model] [,silent:T] [,allin:T or in:logVec]), Model of form
  "y = x1+x2+...+xk" or "y = x1+x2+...+xk - 1", logVec a LOGICAL
  vector of length k
%%%%
@@@@usage#Usage
stepsetup(Model), where Model is a CHARACTER scalar specifying a
regression model, initializes a stepwise regression process.

It creates invisible variable _STEPSTATUS and prints the F-to-enter
statistics and P-values for all the variables in the model.

It returns _STEPSTATUS as value.  This can be assigned (stuff <-
stepsetup("y=x1+x2+x3+x4")) but is not printed.  See topic
'_STEPSTATUS'.

stepsetup(), with no model specified, does the same, except it uses
STRMODEL, the model for the most recent GLM command, as model.

@@@@silent_in_and_allin_keywords#Keywords silent, in and allin
stepsetup([Model,] silent:T) does the same except nothing is printed.

stepsetup([Model,] allin:T [,silent:T]) does the same, except all the
variables are entered in the model immediately so that backward
stepwise regression can be done.

stepsetup([Model,] in:In [,silent:T]) does the same, except the
variables specified by In are entered in the model immediately.  In
must be a LOGICAL vector the same length as the number of independent
variables in Model.  Variable j is considered in the model if In[j] is
True.

@@@@printed_output#Printed output
When variables are initially entered in the model, stepsetup() prints an
overall F statistic and its P-value, Mallow's Cp statistic, adjusted
R^2 and R^2.  The F-statistic tests the null hypothesis that the
coefficients of the "in" variables are 0.

@@@@examples
Examples:
  Cmd> stepsetup("y=x1+x2+x3+x4+x5")
starts stepwise regression process with no variables in the model.

  Cmd> stepsetup("y=x1+x2+x3+x4+x5+x6-1",in:vector(F,F,F,T,T,T))
starts stepwise regression process with variables x4, x5 and x6 in the
model.  Because of "-1" in the model, no intercept is included.

@@@@see_also#Cross references
See also stepstatus(), entervar(), removevar(), steplook().
@@@@______

====stepstatus()#stepwise regression,regression
%%%%
stepstatus([silent:T])
%%%%
@@@@usage#Usage
stepstatus() prints the current stepwise regression model, F-to-remove
statistics with P-values for all the variables currently in the model,
and F-to-enter statistics with P-values for any variables not in the
model.

In addition, if there are any variables in the model, stepstatus()
prints an overall F statistic and its P-value, Mallow's Cp statistic,
adjusted R^2 and R^2.  The F-statistic tests the null hypothesis that
the coefficients of the "in" variables are 0.

It returns _STEPSTATUS as value.  This can be assigned (stuff <-
stepstatus()) but is not printed.  See topic '_STEPSTATUS'.

@@@@silent_keyword#Keyword silent
stepstatus(silent:T) prints nothing but returns _STEPSTATUS as value.
@@@@______

====testbeta()#regression,hypothesis test
%%%%
testbeta(Term, hypValue [,df:T] [,pval:T]), Term a quoted or unquoted
  variable name or an integer term number, hypValue a non-MISSING real
  scalar
%%%%
@@@@usage#Usage
You can use testbeta() to compute a t-statistic to test a null
hypothesis of the form H0: betaj = hypValue, where betaj is a regression
coefficient.  There must be an active regression model, that is you need
to have executed a command of the form regress("y = x" or regress("y=x1
+ x2 + ... xk").

testbeta(Term, betaj0), where Term is the quoted or unquoted name of a
predictor in the regression (for example, "x3" or x3) and betaj0 is a
REAL scalar whose value is the hypothesized value of the regression
coefficient of the predictor.  Term can also be the number of the term
associated with the variable in an ANOVA.

The value is the test statistic Tstat = (betajHat-betaj0)/SE[betajHat].

@@@@pvalue_and_df_keywords#Keywords pvalue and df
testbeta(Term, betaj0, pvalue:T) does the same except the value returned
is structure(tstat:Tstat, pvalue:Pvalue), where Pvalue is the two-tail
P-Value associated with Tstat.

testbeta(Term, betaj0, df:Df [, pvalue:T]), returns
structure(tstat:Tstat, df:Df [,pvalue:Pvalue]), where Df = error
degrees of freedom.

@@@@example
Example:
  Cmd> regress("y = x", silent:T) # output suppressed

  Cmd> testbeta(x, 1, pvalue:T, df:T) # tests H0: beta = 1

@@@@see_also#Cross references
See also secoefs(), twotailt().
@@@@______

====testestim()#regression,hypothesis test
%%%%
testestim(x, hypValue), x REAL scalar or vector with no MISSING
  elements, hypValue a non MISSING REAL scalar
%%%%
@@@@usage#Usage
You can use macro testestim() to test a null hypothesis of the form
H0: E(y|x) = hypValue or E(y|x1,x2,...,xk) = hypValue, where hypValue
is a hypothesized value for the expectation of y for given x values.
You must previously have run regress("y=x") or regress("y=x1+...+xk").

testestim(x, hypValue), where hypValue is the hypothesized value, x is
a scalar (simple linear regression) or a vector of length k (multiple
regression) returns a t-statistic of for testing H0.  It is an error if
there is not an active regression model or if length(x) != k.

@@@@df_and_pvalue_keywords#Keywords df and pvalue
testestim(x, hypValue, df:T) does the same except the result is
structure(tstat:tvalue, df:errorDF), where errorDF is the error degrees
of freedom needed to compute a P-value or find a critical value.

@@@@pvalue_keyword#Keyword pvalue
testestim(x, hypValue [,df:T] , pvalue:T) does the same except the
result is structure(tstat:tvalue [,df:errorDF], pvalue:Pvalue), where
Pvalue is the two-tail P-value associated with the test statistic.

@@@@example
Example:
After regress("y=x1 + x2 + x3"),
  testestim(vector(2,3,4), 17) returns the t-statistic for testing H0:
     E(y|x) = 17
  testestim(vector(2,3,4), 17, pval:T, df:T) returns
     structure(tstat:t_statistic, df:ErrorDF, pval:p_value)

@@@@see_also#Cross references
See also estimlimits(), predlimits(), regpred().
@@@@______

====yhat()#glm,regression,anova
%%%%
yhat() or yhat(Model [,T])
%%%%
@@@@usage#Usage
yhat(), with no argument, computes a REAL matrix of various quantities
useful in making predictions from a regression or analysis of variance
model.

yhat() uses side effect variables RESIDUALS, HII, etc. produced by the
most recent GLM (generalized linear or linear model) command such as
regress() or anova().  If weights were supplied, it also uses the result
of modelinfo(weights:T).

NOTE: If the most recent GLM was not linear, that is, not regress(),
anova(), manova() or their weighted variants, only the first two columns
computed by yhat() are meaningful.

It is an error if any of the needed side effect variables do not exist.

@@@@specified_model#Specified model
yhat(Model) first executes manova(Model, silent:T) to compute the side
effect variables and then computes its usual output.  Model should be a
CHARACTER variable or string specifying a linear model.  Any factors in
the model will be treated as factors.

yhat(Model,T) does the same, except any factors in the model are treated
as variates.

@@@@value_returned#Value returned
Each row of the result corresponds to a case.  When the dependent
variable Y is univariate (has one column), the result has the following
5 columns:
   Col. 1   Y = observed response
   Col. 2   Yhat = predicted or fitted value computed using all data
   Col. 3   Predictive residuals = Y - (Yhat computed excluding the
            case)
   Col. 4   SE[Yhat] = estimated standard error of Yhat as estimate of
            E[Y | x]
   Col. 5   SE[pred] = estimated s.e. of prediction error

When the all the independent variables in the model "y=x1+x2+x3+...+xk"
are variates and not factors, columns 2, 4, and 5 of the output
correspond to components 'estimate', 'SEest', and 'SEpred" in the output
of regpred(hconcat(x1,x2,...)) following regress("y=x1+x2+...+xk").  See
regpred().

@@@@multivariate_model#Multivariate model
When Y is multivariate, with p columns, there are 5*p columns in groups
of p -- the p columns of Y, the p columns of YHat, and so on.

If a case has missing values, most entries will be MISSING and there are
no useful numbers.

@@@@example
Example:
  Cmd> yhat("y=x1+x2+x3+x4")
@@@@______

The output of yhat() is modelled on output printed by the yhat() command
in program Multreg.

yhat() is implemented as a pre-defined macro.

@@@@see_also#Cross references
See also topics regpred(), 'glm', resid().
@@@@______

====_STEPSTATUS*#stepwise regression,regression
%%%%
Type help(_STEPSTATUS) for information about hidden variable
  _STEPSTATUS
%%%%
@@@@description#Description
The macros used for doing stepwise regression communicate among
themselves through "invisible" structure _STEPSTATUS (it's invisible
because its name starts with '_').

In the following description of the components of _STEPSTATUS, k is
the total number of variables in the "full" model with independent
variables x1, x2, ..., xk, and dependent variable y.

_STEPSTATUS is a structure with the following components
  model      CHARACTER scalar which is the current stepwise regression
             model.  It is a legal MacAnova regression model except when
             there are no variables "in" to model and the model does not
             include an intercept
  sscp       k+1 by k+1 matrix of sums of squares and products of x1,
             x2, ..., xk and y after swp() has swept the rows and
             columns of the variables currently in the model
  in         LOGICAL vector of length k with in[j] = T if and only if
             xj is in the current model
  F          2 by k REAL matrix.  If xj is in the current model, column
             j contains the F-to-remove statistic and P-value for xj.
             If xj is not in the current model, column j contains the
             F-to-enter statistic and P-value for xj
  dfe        integer scalar containing the error degrees of freedom for
             the current model
  fullmse    mean square error for full model
  history    vector of integers jj, with 1 <= abs(jj) <= k, representing
             the history of variable selection.  If j = history[i] is
             positive, xj was entered at step i; if j = -history[i] is
             positive, xj was removed at step i.  If keywords 'allin' or
             'in' are used on setupstep to start up with one or more
             variables already in the model, the initial elements of
             history are the numbers of these variables.  Otherwise,
             history is initialized to NULL.

@@@@retrieving_components#Retrieving components
You can retrieve components of _STEPSTATUS using macro steplook().  For
instance, you can get regression estimates from the current model by
  Cmd> regress(steplook(model))

Macros stepsetup(), stepstatus(), entervar() and removevar() all return
_STEPSTATUS as value.  With keyword argument 'silent:T', they print
nothing.

@@@@see_also#Cross references
See also stepsetup(), stepstatus(), entervar(), removevar(), steplook().
@@@@______

_E_O_F_#This should be the last line, an internal End Of File marker
