Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-05 - decus/20-0149/mulreg.doc
There are 2 other files named mulreg.doc in the archive. Click here to see a list.


 
First page


 










          A USER PROGRAM FOR MULTIPLE LINEAR REGRESSION ANALYSIS
          ======================================================







                            Marten van Gelderen


                              Version 5H(246)


                                 16-Feb-80


 























                     Copyright (C) 1975, 1979, 1980 by
                 Foundation Mathematical Centre, Amsterdam
             Institute for Nuclear Physics Research, Amsterdam

    General permission to make fair use in teaching or research of  all
    or  part  of  this material is granted to individual readers and to
    nonprofit organizations, provided  that  the  copyright  notice  of
    either  the  Foundation  Mathematical  Centre  or the Institute for
    Nuclear Physics Research is given and that  reference  is  made  to
    this  publication  and  to the fact that reprinting privileges were
    granted by permission of one of the above mentioned organizations.


CONTENTS                                                               Page


Abstract                                                                 1

Introduction                                                             2

Chapter 1 - Multiple linear regression analysis                          4
            1.   The regression model                                    4
            2.   Least squares                                           5
            2.1  Weighted least squares                                  8
            2.2  Residual analysis                                       9
            3.   Tests of hypotheses                                    10
            3.1  A particular regression coefficient                    10
            3.2  Analysis of variance                                   10

Chapter 2 - The input to the program                                    14
            1.   The model specification                                14
            1.1  Transformations                                        15
            2.   The input specification                                17
            3.   The option specification                               19
            4.   The data specification                                 21
            5.   The user program                                       22
            6.   Examples                                               23

Chapter 3 - The output from the program                                 25
            1.   Standard and optional printed output                   25
            2.   Standard and optional data output                      27
            3.   Error messages                                         28
            4.   Examples                                               30

Appendices  1.   Technical remarks                                      31
            2.   Definition of the syntax of a user program             34
            3.   Technical description of the program                   37

References                                                              41
                                                                     Page 1


ABSTRACT

     Performing  multiple  linear  regression  analysis  on  a   electronic
computer   can   be  very  laborious  when  preparation  of  the  input  is
complicated.   The  program  introduced  in  this  document  is  especially
designed   to   have   flexible   and   comprehensible   model   and  input
specifications.  It accepts a "Model" formula which resembles the  notation
of  regression  models  in common statistical literature quite closely.  An
accompanying "Input" formula provides the program  with  information  about
the  arrangement of the observations in the input "Data", which consists of
a series of numbers in freefield format.  A  "Run"  command  activates  the
program,  while  an  "Exit"  command causes the program to stop.  Extensive
runtime "Help" information is available.  The following  piece  of  program
may serve as an example of some of these ideas:

         "Model"  y  =  alfa0  +  alfa1 * x  +  alfa2 * x ^ 2;
         "Input"  5 * ([x], n, n * [y]);
         "Options" Transformed data matrix, Process submodels (1);
         "Data"   1   4    1.1    0.7    1.8    0.4
                  3   5    3.0    1.4    4.9    4.4    4.5
                  5   3    7.3    8.2    6.2
                 10   4   12.0   13.1   12.6   13.2
                 15   4   18.7   19.7   17.4   17.1
         "Run"
         "Exit"
                                                                     Page 2


INTRODUCTION

     Many computer programs or subroutines  exist,  which  in  one  way  or
another  calculate  estimates  for  the  parameters  of  a  multiple linear
regression model, but most of them  can  hardly  be  used  by  the  layman.
Subroutines  have  to  be  embedded  in  a higher level language program to
perform the necessary input and output;  as a layman however,  one  is  not
always  acquainted  with these languages.  The programs usually require the
input data to be presented in a specific standard form, which  in  practice
means that input data that are already available in a machine readable form
must be transformed into that standard form or that the input  data  should
be  punched  exactly in that standard form.  More difficult still and often
more confusing, is the way in which the program is  told  which  regression
model  the  user  wants  to  consider.   For  instance, for the program for
multiple polynomial regression analysis of the Mathematical Centre  a  code
matrix  of  zeros and ones had to be given in order to indicate the form of
the regression polynomial.  Transformations of one or  more  variables  are
hardly ever automatically possible and require a separate program to be run
in advance to prepare the transformed data matrix.

     In this document a program is described that will make it possible for
almost  everyone  to  obtain  his  desired  results  with  no  more  than a
superficial knowledge of the underlying programming system.   It  has  been
recognized  recently  that standard programs should not bother the user too
much with awkward input specifications.  A statistician is more  interested
in  the  results  of  a  program  and  in  how  he  can obtain them without
interference of software specialists of the computer centre, than in a  few
seconds gain in actual computing time.  Therefore one of the objectives was
to develop a program which would enable the user to specify the form of the
regression model and the arrangement of the numbers in his input data, in a
straightforward manner, even if he has hardly any knowledge of  programming
languages at all.  The program thus accepts a model formula which resembles
the notation of regression models in common  statistical  literature  quite
closely;   the accompanying input formula indicates which numbers or series
of numbers in the input data belong to which variable in the model formula.
                                                                     Page 3


     This input scheme gives the user the opportunity to work with existing
data  and to process possible transformations of the variables in the model
formula without the need for external data adjustment.  However, the  model
formula  must  be  given explicitly and completely and the structure of the
input  data  must  be  known  exactly  before  an  input  formula  can   be
constructed.  In the output the results are identified with the names given
by the user in the model formula.  Any technical (machine dependent) action
required to run the program is described in appendix 1.


ACKNOWLEDGEMENTS

     The basic ideas for the input arrangement presented in  this  document
originated  from  A.P.B.M. VEHMEYER.  The article 'ALGOL 60 translation for
everybody' by F.E.J. KRUSEMAN ARETZ [6] provided  the  foundation  for  the
translator  and  the execution section of the program, while the techniques
used in the storage section were designed by L.J. OOSTRIJK.  The matrix and
least  squares  routines  used  in  the  program were copied from 'ALGOL 60
procedures in numerical algebra, part 1', by T.J. DEKKER [2].  The  classic
book  'Applied Regression Analysis' by N.R. DRAPER & H. SMITH [1] served as
a frame of reference  throughout  the  design  and  implementation  of  the
program.   Numerous  recently  incorporated  statistical  improvements were
suggested by R.D. GILL.
                                                                     Page 4


                                 CHAPTER 1

                    MULTIPLE LINEAR REGRESSION ANALYSIS


1.1  THE REGRESSION MODEL

     In a regression problem the researcher postulates a certain  relation-
ship  between a random variable y (the realizations of which are subject to
some form of disturbance) on  the  one  side  and  a  number  of  variables
x1,...,xp  (which  are  without or at least almost without disturbances) on
the other side.  This relationship is expressed by a mathematical  formula,
which is called the (linear) regression model, for instance:

                  y = a0 + a1 * x1 + ... + ap * xp + e                  (1)

in which a0,...,ap represent unknown regression  coefficients  (parameters)
which  are to be estimated and e represents the disturbance.  If a constant
term is present in the model formula (in (1) the a0), the model is said  to
be  an  'intercept model',  if  no  constant  term is present, the model is
called a 'no-intercept model'.
The variables x1,...,xp and the  variable  y  can  also  represent  (other)
transformed  variables.  The researcher might have reasons to believe (from
background information concerning the experiment) that transformations  are
necessary, for instance:
1) to obtain normally distributed disturbances,
2) to obtain a greater homogeneity of the variances of the disturbances,
3) to linearize non-linear regression models (if possible).
The transformed regression model can be written as:

     G(y) = a0 + a1 * F1(x1,...,xm) + ... + ap * Fp(x1,...,xm) + e      (2)

in which G, F1,...,Fp represent the transformations,
            a0,...,ap represent the parameters to be estimated,
                    y represents the dependent variable,
            x1,...,xm represent the independent variables,
                    e represents the disturbance.
                                                                     Page 5


     The choice of a transformation by means of 'trial and error' is rather
time  consuming and costly.  The importance of the location parameter makes
for the difficulty.  It is not unusual that Log (x) yields no  improvement,
but  that  Log (c+x)  gives  better  results  for a particular choice of c.
Because this holds for almost any transformation  of  some  importance,  we
must  actually  solve  in  each case a nonlinear adjustment problem.  Often
though, a simple form of the transformation is suggested by the  researcher
who is better acquainted with the peculiarities of the experiment.


1.2  LEAST SQUARES

     Regression analysis consists in fact of the adjustment of a hyperplane
of the required dimension to the data.  The fitting is done with the method
of least  squares,  which  means  that  the  sum  of  the  squares  of  the
differences  between the observed values for y and the estimated values for
the expectation of y, are minimized.  This sum of squares  is  also  called
the  residual  sum of squares.  In matrix notation the regression model can
be written as (cf. DRAPER & SMITH [3] pp. 58-62):

                              Y = Xa + e                                (3)

in which Y is a (n*1) random vector of observations,
         X is a (n*p) matrix of known (fixed) values,
         a is a (p*1) vector of (unknown) parameters,
     and e is a (n*1) random vector of disturbances.

It is supposed that E(e) = 0 and var(e) = Isigma^2, in which I is the  unit
matrix, thus:

                               E(Y) = Xa                                (4)

     The sum of squares of the differences between the observed values of Y
and the estimated values for the expectation of Y thus equals:

                 (Y-Xa)'(Y-Xa) = Y'Y - 2a'X'Y + a'X'Xa                  (5)

(for a'X'Y is a scalar and therefore equal to Y'Xa).
                                                                     Page 6


     Choosing as least squares estimator b that value of a which  minimizes
(5),  involves  differentiating  with  respect  to  the  elements  of a and
equating the result to zero:

                   -2X'Y + 2X'Xb = 0,  thus:  X'Y = X'Xb                (6)

This system is called the normal equations.  If the rank of X equals p, X'X
is nonsingular and the inverse of X'X exists.  In that case the solution of
the normal equations can be written as:

                             b = inv(X'X)X'Y                            (7)

Observe that p <= n must hold, in order that the rank of X can be p at all.
Therefore  at  least  as  many  observations  must  be  made,  as there are
parameters in the model.  Also observe that E(b) = inv(X'X)X'E(Y) = a, thus
b is an unbiased estimator of a.

The least squares estimator has the following properties:
1. It is an estimator which minimizes the sum  of  squares  of  deviations,
   irrespective  of  any  distribution properties of the disturbances.  The
   assumption that the disturbances are normally distributed is, of course,
   necessary  for  tests  which depend on this assumption, such as t- or F-
   tests,  or  for  obtaining  confidence  intervals  based  on  t-  or  F-
   distributions.
2. According to the Gauss-Markov theorem, the elements of  b  are  unbiased
   estimators,  which  have minimum variance (of any linear function of the
   Y's which provides  unbiased  estimators),  again  irrespective  of  the
   distribution properties of the disturbances.
3. If the disturbances are mutually independent  and  normally  distributed
   (with  E(e) = 0  and  var(e) = Isigma^2),  then  b  is  also the maximum
   likelihood estimator.

The variance-covariance matrix of b is:

                         var(b) = inv(X'X)sigma^2                       (8)

The variances  are  the  diagonal  and  the  covariances  the  off-diagonal
elements.
                                                                     Page 7


An unbiased estimator for sigma^2 is given by:

                       s^2 = (Y'Y - b'X'Y) / (n-p)                      (9)

The square root of this estimator is frequently called 'standard  error  of
estimate'.   In  the  printed  output  of  the program it is indicated more
properly as 'standard deviation of the error term'.

     Let vij be the element in the i-th row and j-th  column  of  inv(X'X),
then  sdi = s * Sqrt(vii) estimates the standard deviation of bi, and cij =
vij / Sqrt(vii * vjj) gives the correlation coefficient between bi  and  bj
for i = 1,...,p and j = 1,...,p.  Thus:

                           vii = (sdi / s)^2                           (10)
and
          vij = cij * Sqrt(vii * vjj) = cij * (sdi * sdj) / s          (11)


     A frequently used statistical measure for evaluating regression models
is the multiple correlation coefficient R which is defined in the intercept
model as the square root of the proportion of the corrected  total  sum  of
squares accounted for by the model.  If the correction for means is denoted
by nu^2, with u = Sum(i,1,n,yi)/n, then R can be defined by:

       R^2 = (b'X'Y-nu^2)/(Y'Y-nu^2) = 1 - (Y'Y-b'X'Y)/(Y'Y-nu^2)      (12)

However, we must divide Y'Y-b'X'Y by n-p, not by n, to obtain  an  unbiased
estimator  of  sigma^2, moreover it is customary to divide Y'Y-nu^2 by n-1,
not by n.  If we adopt both modifications we obtain the  adjusted  multiple
correlation coefficient, which can thus be defined by:

           adj(R)^2 = 1 - (n-1)/(n-p) * (Y'Y-b'X'Y)/(Y'Y-nu^2)         (13)

     In the no-intercept model the correction for means is ignored,  giving
as definition of R^2:  b'X'Y/Y'Y = 1 - (Y'Y-b'X'Y)/Y'Y,  while the adj(R)^2
is defined correspondingly as:  1 - n/(n-p) * (Y'Y-b'X'Y)/Y'Y.  R^2  itself
is  often  called  the  'proportion of variation explained' (cf. THEIL [11]
pp. 178-179).
                                                                     Page 8


1.2.1  WEIGHTED LEAST SQUARES

     It sometimes happens that some of the observations for  the  dependent
variable  are  'less  reliable'  than  others.  This usually means that the
variances of the observations are not all equal;  in other words the matrix
V = var(e)  is  not  of  the  form  Isigma^2,  but is diagonal with unequal
diagonal elements.  The basic idea to solve this problem is, to transform Y
to  other  variables,  which do appear to satisfy the usual tentative model
assumptions,  and  then  apply  the  usual  (unweighted)  analysis  to  the
variables  so obtained.  The estimates can then be re-expressed in terms of
the original variables Y (cf. DRAPER & SMITH [3] pp. 77-81).

     Let the original regression model be:  Y = Xa + e, with  E(e) = 0  and
var(e) = Vsigma^2,  with V diagonal with unequal diagonal elements, and let
P = inv(V).  Premultiplying the original regression model with  Q = Sqrt(P)
gives as transformed regression model:

                             QY = QXa + Qe                             (14)

with E(Qe) = 0 and var(Qe) = Isigma^2.  The normal equations then become:

                           (QX)'QY = (QX)'QXa                          (15)

giving as solution if the indicated inverse matrix exists:

                b = inv((QX)'QX)(QX)'QY = inv(X'PX)X'PY                (16)

with variance-covariance matrix:

                       var(b) = inv(X'PX)sigma^2                       (17)

     In practical situations it  is  often  difficult  to  obtain  specific
information  on  the  form  of V at first.  For this reason it is sometimes
necessary to make the (known to be erroneous)  assumption  V = I  and  then
attempt  to  discover  something  about  the  form  of  V  by examining the
residuals from the regression analysis.
                                                                     Page 9


1.2.2  RESIDUAL ANALYSIS

     The vector of residuals D is defined as  the  difference  between  the
vector  of  observations  Y  and the vector of fitted values Z, obtained by
using the regression equation  Z = Xb.  So D = Y - Z  or  di = yi - zi  for
i = 1,...,n.   If  the model is correct, the residual mean square MSE = s^2
estimates sigma^2, and the estimated standard deviation of the fitted value
zi at xi = (xi1,...,xip)' is:

                     sd(zi) = s * Sqrt(xi'inv(X'X)xi)                  (18)

which can be used to construct a confidence interval for the expected value
of  yi: E(yi) at xi = (xi1,...,xip)', or to construct a prediction interval
for the mean of h new observations at this  point  (cf. DRAPER & SMITH  [3]
pp. 121-122).  In the first case the confidence interval is:

           zi +- t(n-p-1,1-alpha/2) * s * Sqrt(xi'inv(X'X)xi)          (19)

and in the second case the prediction interval is:

        zi +- t(n-p-1,1-alpha/2) * s * Sqrt(1/h + xi'inv(X'X)xi)       (20)

     Researchers often divide the residuals  di  by  s,  resulting  in  the
standardized residuals, which can be examined to see if they make it appear
that the assumption ei/sigma ~ N(0,1) is violated  (cf. DRAPER & SMITH  [3]
pp. 86-97).  It might be expected that roughly 95% of the di/s were between
the limits (-2,2).

     However, the variances  of  the  residuals  are  not  constant  but  a
function of the X matrix (see (18)), which suggests as standardization:

                   ti = di / s / Sqrt(1 - xi'inv(X'X)xi)               (21)

giving the studentized residual.  The maximum studentized residual  can  be
used  in  a  test for detecting outliers, as follows:  let t^2 = max(ti^2),
then  min(1, n * (1-Fisher(1, n-p-1, t^2*(n-p-1)/(n-p-t^2)))) is an  'upper
bound  for  the  right tail probability of the largest absolute studentized
residual' (cf. LUND [9] pp. 473-474).
                                                                    Page 10


1.3  TESTS OF HYPOTHESES

1.3.1  A PARTICULAR REGRESSION COEFFICIENT

     If the disturbances are mutually independent and normally  distributed
(with E(e) = 0 and var(e) = Isigma^2) and with a (preset) level of signifi-
cance alpha, a significance test for a  particular  regression  coefficient
can be performed, or more specifically:  the null hypothesis is:

          H0:  ai = 0 (given that all other aj are in the model),

which is tested against the alternative hypothesis:

          H1:  ai is not equal to zero,

by treating FRi = bi^2/var(bi) as a realization of a Fisher(1,n-p) variate.

     However, this test  must  be  used  with  caution,  because  with  the
(preset)  level  of  significance alpha, only one coefficient can be tested
properly, while the computer output lists statistics for all  coefficients.
It seems very tempting to test the coefficients serially one at a time, but
one must keep in mind that in doing so the level  of  significance  of  the
whole test rises above the nominal value (cf. DRAPER & SMITH [3] p. 65).


1.3.2  ANALYSIS OF VARIANCE

     In the analysis of variance table the different contributions  to  the
total uncorrected sum of squares Y'Y (which is the first part in the table)
are given (cf. DRAPER & SMITH [3] pp. 57 & 62).

     The second part of the table assumes  the  presence  of  an  (unknown)
constant  term  in  the  model;   if  this  term is absent, the 'mean'-line
disappears and in the 'regression'-line p-1 changes into p  and  b'X'Y-nu^2
changes into b'X'Y.

     The third part of the table is only present when repeated observations
for the dependent variable are available, in which case:
                                                                    Page 11


k  is the number of groups of replications,
mi is the number of replications in group i, and
W = (w1,...,wk)', with wi = Sum(j,1,mi,yij) / Sqrt(mi), for i = 1,...,k.

     The fourth part of the  table  is  only  present  if  a  reduction  is
requested  and  possible.   SSQ then stands for the residual sum of squares
from a regression analysis with the first p-q  out  of  the  p  independent
variables (1 <= q <= p-1), while SSE stands for the residual sum of squares
from a regression analysis with the original p independent variables.

                           Analysis of variance
source of                                                       right  tail
variation        df   sum of squares    mean square   F-ratio   probability
---------------------------------------------------------------------------
total            n          Y'Y
---------------------------------------------------------------------------
mean             1         nu^2          MSM = nu^2     FRM      P(FM>=FRM)
regression      p-1    b'X'Y - nu^2         MSR         FRR      P(FR>=FRR)
residual        n-p     Y'Y - b'X'Y      MSE = s^2
---------------------------------------------------------------------------
 lack of fit    k-p     W'W - b'X'Y         MSL         FRL      P(FL>=FRL)
 pure error     n-k      Y'Y - W'W          MSP
---------------------------------------------------------------------------
 reduction       q       SSQ - SSE          MSQ         FRQ      P(FQ>=FRQ)
---------------------------------------------------------------------------

     The column 'mean square' is  obtained  by  division  of  the  sums  of
squares by their corresponding degrees of freedom.  The column 'F-ratio' is
obtained by division of the mean  squares  by  the  residual  mean  square,
except  for  the  lack of fit F-ratio, which is obtained by division of the
lack of fit mean square by the pure error mean square, thus:

   MSM = nu^2/1, MSR = (b'X'Y-nu^2)/(p-1), MSE = s^2 = (Y'Y-b'X'Y)/(n-p),

     MSL = (W'W-b'X'Y)/(k-p), MSP = (Y'Y-W'W)/(n-k), MSQ = (SSQ-SSE)/q,

       FRM = MSM/MSE, FRR = MSR/MSE, FRL = MSL/MSP and FRQ = MSQ/MSE.
                                                                    Page 12


     If the disturbances are mutually independent and normally  distributed
(with E(e) = 0 and var(e) = Isigma^2) and with a (preset) level of signifi-
cance alpha, a significance tests can be performed for:

1.  The mean of the observations for the dependent variable, or more speci-
    fically:  the null hypothesis is:

          H0:  E(u) = 0,

    which is tested against the alternative hypothesis:

          H1:  E(u) is not equal to zero,

    by treating FRM as a realization of a Fisher(1,n-p) variate.

2.  The regression equation, or more specifically:  the null hypothesis is:

          H0:  a1 = ... = ap = 0, except for the ai that denotes
                                  the constant term (if present),

    which is tested against the alternative hypothesis:

          H1:  at least one of a1,...,ap is not equal to zero,

    by treating FRR as a realization of a Fisher(p-1,n-p) variate.

3.  The adequacy (linearity) of the model, or more specifically:  the  null
    hypothesis is:

          H0:  the linear model is adequate (that is, no model signifi-
               cantly improves the prediction of Y over the linear model),

    which is tested against the alternative hypothesis:

          H1:  the linear model is not adequate,

    by treating FRL as a realization of a Fisher(k-p,n-k) variate.
                                                                    Page 13


4.  A subset of regression coefficients,  or  more  specifically:   suppose
    without  loss  of  generality  that  the  subset consists of the last q
    coefficients, then the null hypothesis is:

          H0:  ar = ... = ap = 0, with r = p-q+1,

    which is tested against the alternative hypothesis:

          H1:  at least one of ar,...,ap is not equal to zero,

    by treating FRQ as a realization of a Fisher(q,n-p) variate.

5.  A linear combination of the regression coefficients, or  more  specifi-
    cally:  the null hypothesis is:

          H0:  c'a = m, in which c is a vector of constants with order q+1,

    which is tested against the alternative hypotheses:

          H1:  c'a is not equal to m,

    by substituting c'a = m in the original model, shifting the known terms
    to  the  left hand part, combining the corresponding terms in the right
    hand part, and testing the thus derived so called  'reduced  model'  by
    treating FRQ as a realization of a Fisher(q,n-p) variate.

In each case the right tail probability P(F >= FR) can be found in the last
column of the analysis of variance table (cf. DRAPER & SMITH [3] pp. 63-64,
68 & 72-75).
                                                                    Page 14


                                 CHAPTER 2

                         THE INPUT TO THE PROGRAM


2.0  The purpose of the input system is to  give  the  user  a  simple  and
adequate  formalism to tell the program what he wants.  In order to specify
the regression model and the correspondence between the  variables  in  the
model  and the numbers in the input data, the user must provide a so called
user program, which consists of one or more jobs which in  their  turn  are
made  up out of one or more specifications.  Each specification starts with
a keyword and terminates at the next keyword.  The keywords are:   "Model",
"Input", "Options", "Data", "Run", "Exit" and "Help" (quotes included).


2.1  THE MODEL SPECIFICATION

     To let the program  know  between  which  variables  the  statistician
expects  a  certain  kind of relationship, he must provide a model specifi-
cation, which consists of the keyword "Model" followed by  a  formula  (the
model  statement),  which  resembles  the  notation of regression models in
common statistical literature quite closely.  For instance:

             "Model"  y = alpha0 + alpha1 * x1 + alpha2 * x2;

     A model formula consists of an  identifier  to  denote  the  dependent
variable  (the left hand part), followed by an '=' (equal), followed by the
sum of a number of terms (the right hand part), while it is terminated with
a  ';' (semicolon).   Each  term  must  be  the product of an identifier to
denote the parameter (which is to be estimated) and an identifier to denote
the  independent  variable.  An exception is made for the optional constant
term, which is given as a single identifier denoting  that  constant  term,
and which may be placed anywhere in the model.

     Each identifier must start with a letter and is allowed to contain any
number  of  letters,  digits and blanks.  As most peripheral equipment of a
computer is unable to process sub- or superscriptions or Greek letters,  we
write alpha0, alpha1 and alpha2.  Identifiers have no inherent meaning, but
                                                                    Page 15


serve for the identification of variables, parameters and functions.   They
may  be chosen freely (except for the twentyone standard function names and
the ten option names, cf. section 2.1.1 and section 2.3).   It  is  advised
not  to  use  the  same  identifier  to  denote  two  (or  more)  different
quantities;  for regression parameters, however, it will not lead to  fatal
errors, whereas for the dependent and independent variables distinguishable
identifiers must be used indeed.  Correct model formulae are for instance:

          "Model" y variable = constant + parameter * x variable;
and
          "Model" depvar = const + beta1 * xvar1 + beta2 * xvar2;


2.1.1  TRANSFORMATIONS

     Almost all transformations a user would like to perform on  his  input
data  fit  quite  naturally  in  the model formula:  each transformation is
expressed as a formula itself.  If, for instance, the user wants to include
in  the  model  formula as an independent variable the natural logarithm of
the sum of two other variables, called xvar1 and xvar2, he writes:

                            Ln (xvar1 + xvar2)

In model formulae the operators  '+' (plus),  '-' (minus),  '*' (asterisk),
and  '/' (slash)  are  allowed,  all  with  their  conventional  meaning of
addition, subtraction, multiplication and division respectively.  Of course
the  normal  operator  precedence rules are obeyed.  Special operators are:
':' (colon), integer division and  '^' (uparrow), exponentiation.
The operation term : factor is defined  only  for  operands  both  of  type
integer  and  will  yield  a  result of type integer, with the same sign as
would be obtained by normal division,  while  the  magnitude  is  found  by
dividing  the  two quantities and taking the whole part;  mathematically it
can be defined as:  a : b = Sign (a / b) * Entier (Abs (a / b)),
for instance:  5 : 2 = 2 and -7 : 2 = -3.
The operation factor ^ primary denotes exponentiation, where the factor  is
the base and the primary is the exponent,
for instance:  5 ^ 2 = 25 and 2 ^ 3 ^ 2 = 64 but 2 ^ (3 ^ 2) = 512.  
                                                                    Page 16


Also the following twentyone standard functions are allowed:

   Abs (E),  Sign (E),  Sqrt (E),  Sin (E),  Cos (E),  Tan (E),  Ln (E),
   Log (E),  Exp (E), Entier (E), Round (E), Mod (E1, E2), Min (E1, E2),
   Max (E1, E2), Arcsin (E), Arccos (E), Arctan (E), Sinh (E), Cosh (E),
   Tanh (E) and Indicator (E1, E2, E3)

in which E, E1, E2 and E3 are expressions in terms of variables,  operators
and standard functions.  Round (E) is defined as:  Entier (E + 0.5) and
Indicator (E1, E2, E3) is defined as:  IF E1 <= E2 <= E3 THEN 1 ELSE 0.

     The dependent variable may be transformed in a similar way  and  as  a
consequence the model formula in its most general form looks like:  

   "Model" G (y) = b0 + b1 * F1 (x1,...,xm) + ... + bp * Fp (x1,...,xm);

Some examples of transformed model formulae are:

          "model" y = a0 + a1 * Sqrt (x1 + x2) + a2 * Sqrt (x3);
and
          "MODEL" Arcsin (Sqrt (Y)) = A0 + A1 * X + A2 * X ^ 2;

     A user can specify model formulae in which terms with known regression
coefficients  appear, by subtracting those terms from the left hand part of
the model formula, for instance:

           "Model" y - 5.4321 * x3 = a0 + a1 * x + a2 * x ^ 2;

This applies especially to the constant term; if this term is known it must
be shifted to the left hand part.

     If weights are present in the input data (or can be  computed  out  of
the input data), to indicate that the variances of the observations are not
all equal (cf. section 1.2.1), the left hand part of the model formula  can
be  expanded  with  a  so  called weight part (which can be an expression),
preceeded by a '&' (ampersand), for instance:

   "Model" Depvar & Max (Abs (Weight), 10) = Const + Param * Indepvar;
                                                                    Page 17


2.2  THE INPUT SPECIFICATION

     To indicate which numbers or series  of  numbers  in  the  input  data
belong  to  which  variable  in  the model formula and which numbers can be
skipped, the program expects an input specification.  It  consists  of  the
keyword  "Input",  followed  by  a  formula  (the  input  statement)  which
describes the arrangement of the observations in the input data,  while  it
is  terminated with a ';' (semicolon).  The basic idea is that numbers from
the input data are identified with the names from the input formula in such
a  way  that (in order of entry) numbers belonging to the same name are put
in a queue appended to that name.  For instance:

         "Input" 100 * (codenr, 10 * [yvar], [xvar1, xvar2], -1);

means that one hundred series of numbers (each, as a check,  terminated  in
this example by -1) are present in the input data.  Each series consists of
fourteen numbers:  first one value which is read and assigned to  the  name
codenr,  then  ten  values  for  the name yvar, then one value for the name
xvar1, followed by one value for the name xvar2 and finally the value -1.

     The basic constituent of an input formula is a  variable  enclosed  in
square brackets, in the example: [yvar].  The corresponding number from the
input data will be appended to the queue for that name.  Several  variables
can  be  put  together  in a variable list by separating them by commas and
enclosing them in square brackets, in  the  example: [xvar1, xvar2].   This
only serves to save the writing of several opening and closing brackets.

     Separate numbers, series or  blocks  of  numbers  can  be  treated  by
putting a repetition factor (control) followed by an asterisk in front of a
variable list (or in front of an input formula which must then be  enclosed
in parentheses), in the example:  100 * and 10 * .

     If a repetition factor is 1, it  may  be  omitted  together  with  the
asterisk  and  a  parentheses  pair,  but square bracket pairs must remain.
When a name is used as a repetition factor, a value must already have  been
assigned  to it, which is done by giving that name, without square brackets
and followed by a comma (or closing  parenthesis),  earlier  in  the  input
formula  than  the  use  of  that name as a repetition factor.  The corres-
                                                                    Page 18


ponding number from the input data is then assigned  as  a  value  to  that
name.   If  such  names  are  used  repeatedly  in  the  input formula, the
corresponding numbers from the input data are compared with the  first  one
and,  in  the  case  of inequality, an error message is supplied.  This may
serve as a check against shifted data reading.   A  similar  check  can  be
obtained  by  giving  an  explicit  number  followed by a comma (or closing
parenthesis), in the example: the -1.  The corresponding  number  from  the
input  data  is  then  compared  with that given number and, in the case of
inequality, an error message is produced.

     Also an expression is allowed as a  repetition  factor,  or  for  that
matter,  as  a check value, provided that it is enclosed in angle brackets,
for instance:  <k+n>.  As in the case of single names used as a  repetition
factor  each  (non-standard function  and  non-option)  name used in such a
(special) expression must have been given, followed by a comma, earlier  in
the input formula than the use of that name in the expression.

     The linkage between  the  model  formula  and  the  input  formula  is
established  by  using  the  same names in the model terms and in the input
variable lists.  Numbers from the input data  that  belong  to  such  input
names  will  be  treated  as  observations  for  the model variables, while
numbers that belong to input names between square  brackets  which  do  not
appear in the model formula, are skipped.

     Often, repeated observations for the dependent variable are available.
In  order  to  be  able  to process these observations automatically, it is
necessary that a variable list consisting entirely of  dependent  variables
is  preceeded  by  a repetition factor (followed by an asterisk) indicating
the number of repetitions.  If a variable list contains independent as well
as  dependent  variables, the number of replications is assumed to be 1.  A
series  of  (say  100)  observations  for  a  dependent  variable  with  no
replications is denoted as:

                             100 * ([dep var])

The repetition factor in front of the opening  square  bracket  is  omitted
(because  it  is  1),  although  the  parentheses  are  not.   Without  the
parentheses it would mean 100 replications of [dep var].
                                                                    Page 19


EXAMPLE    "Input" k, n, <k+n> * (c, m, m * [y], [x1,x2,x3,x4], c), -99;

means that: the first number is read and its value assigned to k,
            the next number is read and its value assigned to n,
            then k+n times the following happens:
                 a number is read and its value assigned to c,
                 the next number is read and its value assigned to m,
                 then the m replications for y are read,
                 next the observations for x1, x2, x3 and x4 are read,
                 then a number is read and its value compared with c,
            finally a number is read and its value compared with -99.

If the comparisons fail, an error message is supplied and execution of  the
job  is terminated, otherwise (k+n) observations for x1, x2, x3, x4 and for
each quadruple m replications for y, have been identified.


2.3  THE OPTION SPECIFICATION

     It is possible to have the program perform some  tasks  optionally  by
providing  an  option  specification  in a job.  It consists of the keyword
"Options" followed by a list of option identifiers or corresponding  option
numbers  (the  option statement), separated by commas and terminated with a
';' (semicolon).  The following ten options are available:

               option number              option name

                     1                  Transformed data matrix
                     2                  Correlation matrix
                     3                  Residual analysis
                     4                  No regression analysis
                     5                  Process submodels
                     6                  Print input data
                     7                  No input data rewind
                     8                  Save original model
                     9                  Test reduced model
                    10                  Missing values
                                                                    Page 20


     Options 1, 2, 3 and 6 cause the corresponding piece of information  to
be  printed.   However,  option  1  lists only those (possibly transformed)
variables that are present in the model formula in  a  neat  tabular  form,
while  option  6 lists all the original input data serially (eleven numbers
per line) without any special layout, because the input data  consists  (by
definition) of an unstructured series of numbers (cf. section 2.4).

     Option 4 suppresses the regression analysis;  it is meant to  be  used
in combination with option 1 and/or 2.

     Option 5 causes the program to process submodels, which are formed  by
a  form  of  backward  elimination:  each time the last term from the right
hand part from the model formula is omitted, by deleting  the  last  column
from  the  design  matrix,  and a regression analysis is performed with the
reduced design matrix.   Messages  are  generated  about  which  terms  are
omitted,  while  further  processing  of  the job ceases when the resulting
model formula is of the form:  y = c.  Moreover a test is made  (under  the
usual  assumptions)  whether the omitted terms did contribute significantly
to the regression sum of squares (cf. section 1.3.2.4).

     To option 5 a specifier list may be  appended,  to  prevent  the  pro-
duction of waste output for unwanted submodels.  In this list the number of
terms to be omitted from the model formula (counting backwards, starting at
the  end)  must  be given enclosed in parentheses.  For example the option:
process submodels (6,  10)  instructs  the  program  to  process  only  two
submodels,  one  with  the last six terms omitted and one with the last ten
terms omitted (from the original model formula).  If the user asks for more
terms to be omitted than are present in the model formula, an error message
is supplied and the execution of that job is terminated.  Moreover,  if  no
explicit  specifier list is appended to option 5, the options 2 and 3 yield
no effect (even if specified), which is also to prevent the  production  of
waste output for the submodels.

     Option 7 gives the user the opportunity to process consecutive  pieces
of  input  data  in consecutive jobs.  Normally the processing of the input
data for each job starts with the first number in  the  data  specification
(or  with  the  first  number  in  the datastream), and the program gives a
(warning) message if the input  formula  does  not  match  the  input  data
                                                                    Page 21


precisely.   This  option  disengages the message and causes the program to
continue processing input data where the previous job had finished.  

     Option 8 causes the residual degrees of freedom and  residual  sum  of
squares  from  the current job to be saved, in order to be able in the next
job, by means of specifying option 9,  to  test  whether  the  model  under
consideration  in  that  next job, shows a significant increase in residual
sum of squares in comparison with the model in the previous job.  In effect
this  gives  the  possibility  of  testing a hypothesis concerning a linear
combination of the parameters  from  a  model  (cf. section  1.3.2.5),  for
instance (cf. SEARLE [11] pp. 121-122):

             "Model 1"  y = b1 * x1 + b2 * x2 + b3 * x3;
             "Options"  Save original model;
             "Run"
             "Model 2"  y - 4 * x1 = b2 * (x1 + x2) + b3 * x3;
             "Options"  Test reduced model;
             "Run"

causes the null hypothesis: b1 = b2 + 4 to be tested (in the second job).

     Option 10 may be used to identify  some  observations  or  repetitions
as 'missing'.   In  a  specifier list, appended to this option, the missing
values must be given enclosed in parentheses.  When a repetition equal to a
missing  value  is  encountered in the input data, the corresponding set of
repetitions for the dependent variable(s) is not  included  in  the  design
matrix.  When an observation equal to a missing value is encountered in the
input data, or when none of the repetitions  are  included  in  the  design
matrix, the corresponding set of observations for the independent variables
together with the (possibly empty) set of  repetitions  for  the  dependent
variable(s) (i.e. the 'case') is not included in the design matrix.


2.4  THE DATA SPECIFICATION

     The data specification consists of an unstructured series  of  numbers
(the data statement), preceeded by the keyword "Data" and terminated at the
next keyword.  The structure is imposed onto it by the input formula.  
                                                                    Page 22


     A sequence of symbols is considered a number  when  it  satisfies  the
definition  of  number  in  appendix 2, together with the machine dependent
restrictions imposed by the underlying  SIMULA  programming  system.   Note
that  the  definition  of number does not allow FORTRAN-like numbers as for
instance  10. or 1# (at least one digit must follow).   It  is  recommended
always  to  use  blanks as delimiters between numbers, but the use of other
non-numerical symbols as delimiters will not lead to fatal errors, only  to
a slight increase in processing time.

     If a datafile is specified in response to the datastream  request,  or
is  appended to the "Run" keyword (cf. appendix 1), the program will try to
read a record of input data, which then do not have to be preceeded by  the
keyword "Data", from that file.  The record ends at the next "Eor" keyword.
However, a non empty data  specification  in  the  user  program  will  get
priority  over  reading  input data from the specified file, while an empty
data specification causes the program to start reading the next  record  of
input  data  from the specified file.  If the data specification as well as
the next nonempty record in the data file does not  contain  any  numerical
information, an error message is supplied.

EXAMPLES            real number             value

                       1.234                  1.234
                        .98                   0.98
                      -0.5673#2             -56.73
                       +.02#-1                0.002
                        #+3                1000.0


2.5  THE USER PROGRAM

     In a user program several jobs can be submitted to the program.   Each
job  is  separated  from  its preceding one by the keyword "Run", while the
entire user program is terminated with the keyword "Exit".   In  the  first
job,  the  model, input and data specification must be given in some order.
The option specification is, of course, optional.  In each following job  a
specification which is not changed may be omitted, the program then retains
the last given specification.  If options have been specified in a previous
                                                                    Page 23


job  and  one  wants to delete them, this is done by providing a new option
specification which may be empty if no options are to be executed (that  is
by only providing:  "Options";).

     In front of each job or in front of the keyword "Exit" a text  may  be
given for further identification of the output of a job or of the output of
the entire user program.  The use of quotes in that text should be  avoided
in  view  of  confusion  with  the  keywords.  The program starts reading a
(possibly empty) text at the beginning of the next line after  the  keyword
"Run"  of  the previous job (with the first job the program starts with the
first line in the inputstream).


2.6  EXAMPLES

     The following user program can be submitted without  any  modification
to  the  Multiple  Linear Regression Analysis program.  It consists of four
jobs, three of which are preceeded by  an  identifying  header,  while  the
whole user program ends with an identifying trailer:


**********************************
*  Example  1  originates from:  *
*  DE JONGE [4], pp. 472 & 479.  *
**********************************

"Model"   y  =  c * Log (x)  +  a  +  b * x;
"Input"   5 * ([x], 10 * [y]);
"Options" Transformed data matrix, Correlation matrix,
          Residual analysis, Process submodels (1, 2);

"Data"    25  0.67  0.70  0.75  0.76  0.78  0.80  0.83  0.84  0.88  0.89
          50  0.88  0.92  0.93  0.96  0.98  1.00  1.01  1.03  1.06  1.07
          80  0.96  0.98  0.99  1.03  1.05  1.06  1.08  1.11  1.15  1.17
         130  1.07  1.09  1.11  1.13  1.14  1.14  1.19  1.22  1.25  1.29
         180  1.10  1.13  1.17  1.19  1.20  1.21  1.23  1.25  1.28  1.33
"Run"
                                                                    Page 24


**********************************
*  Example  2  originates from:  *
*  SEARLE [11], pp. 121-123      *
**********************************

"Input"  5 * [y, x1, x2, x3];
"Data"   8    2    1   4
        10   -1    2   1
         9    1   -3   4
         6    2    1   2
        12    1    4   6

"Model 1"  y  =  a3 * x3  +  a2 * x2  +  a1 * x1;
"Options"  Save original model, Process submodels (1);
"Run"
"Model 2"  y - 4 * x1  =  b2 * (x1 + x2) + b3 * x3;   (eqn. 118, p. 121)
"Options"  Test reduced model, Transformed data matrix;
"Run"

****************************************
*  Example  3  originates from:        *
*  AFIFI & AZEN [1], pp. 88 & 93-100.  *
****************************************

"Model"  y  =  alfa0  +  alfa1 * x;
"Input"  5 * ([x], n, n * [y]);
"Option" Transformed data matrix, Print input data;

"Data"   1   4    1.1    0.7    1.8    0.4
         3   5    3.0    1.4    4.9    4.4    4.5
         5   3    7.3    8.2    6.2
        10   4   12.0   13.1   12.6   13.2
        15   4   18.7   19.7   17.4   17.1
"Run"

***  Marten van Gelderen; Mathematisch Centrum  ***

"Exit"
                                                                    Page 25


                                 CHAPTER 3

                        THE OUTPUT FROM THE PROGRAM


3.1  STANDARD AND OPTIONAL PRINTED OUTPUT

     After having read the keyword "Run", the  processing  of  the  job  is
initiated.   First  the  model,  input and option texts are printed in this
order.  Next an attempt is made at translating the specifications.   Errors
against  syntax  or semantics cause error messages to be printed below each
specification, while further processing of that job ceases.  Note that  the
processing  of the next job, if present, will be of little or no use unless
the specification which  developed  the  error(s)  is  changed.   Next  the
(transformed)  data matrix is formed and passed to the regression routines,
which supply the following printed output in the order indicated:

    1) a listing of the original input data (option 6),
    2) the (transformed) data matrix (option 1),
    3) per (transformed) variable the:
       mean, standard deviation, minimum and maximum,
    4) the correlation matrix of the (transformed) variables (option 2),
    5) the multiple correlation coefficient (with adjustment),
    6) the proportion of variation explained (with adjustment),
    7) the standard deviation of the error term,
    8) the estimates for the regression parameters with
       estimated standard deviation, F-ratio and right tail probability,
    9) the correlation matrix of the estimates (option 2),
   10) the analysis of variance table,
   11) the residual analysis (option 3).

 Ad 1) cf. section 2.4.
 Ad 2) The transformed data matrix gives  the  input  data  after  possible
       transformations  according  to  the  model  specifications have been
       applied.  If the model  formula  contains  no  transformations,  the
       original input data are given.  The dependent variable is given as a
       separate column.  In the case  of  replications  for  the  dependent
       variable,  the  mean  value  of  them  is  given,  and the number of
                                                                    Page 26


       replications is given as an  extra  (last)  column.   If  a  weight-
       variable  (or  -expression)  is  specified in the model formula, the
       (transformed) data comprising the weights  are  given  as  an  extra
       (last) column.  Each (transformed) independent variable is indicated
       by its corresponding parameter.  This originated from the fact  that
       it is not  obvious  how  to  denote  a variable which is transformed
       like:  Arcsin (Sqrt (y+25)), with 'Arcsin', with 'Sqrt'  or  perhaps
       with 'y' itself.  The dependent variable is indicated by 'dep.var.'.
 Ad 4) and 9) The matrix of the estimated correlation coefficients  of  the
       variables  and  of  the  estimates  are  both  supplied depending on
       whether option 2 is specified or not.
 Ad 5), 6) and 7) cf. section 1.2.
 Ad 8) The F-ratio and right tail probability give the user the opportunity
       to  test  the  significance  of  a particular regression coefficient
       (cf. section 1.3.1).
Ad 10) The layout of the table closely  resembles  that  of  the  table  in
       section  1.3.2.   The F-ratios and right tail probabilities give the
       user the opportunity to test the significance of all the  regression
       coefficients  or  of  a subset or combination thereof or to test the
       adequacy of the (linear) model (cf. section 1.3.2).
Ad 11) A table of observations, fitted values, standard deviations  of  the
       fitted  values,  residuals,  standardized  residuals and studentized
       residuals is provided (cf. section 1.2.2).  As a check  on  computa-
       tions,  the  sum  of  the  residuals  is  also given.  If an unknown
       constant term is present in the model formula, this  sum  should  be
       zero.   Furthermore the upperbound for the right tail probability of
       the largest absolute studentized residual is given.

     Without  options  specified,  the  printed  output  from  the  program
consists  of  3),  5),  6),  7), 8) and 10).  If option 5 is specified, the
output for the model itself is given as specified by the other options, but
for  the  submodels  it  depends  on  the use of a submodel specifier list.
Without that list the output from the options 1,  2  and  3  is  suppressed
(even if those options are specified).  With that list only the superfluous
parts  of  the  output  (that  is  the  transformed  data  matrix  and  the
correlation matrix of the variables) are suppressed.
                                                                    Page 27


3.2  STANDARD AND OPTIONAL DATA OUTPUT

     If an outputfile is specified in response to the outputstream request,
or  is  appended  to the "Run" keyword (cf. appendix 1), the program writes
the following pieces of information in one record to that file:
1) if option 1 is specified:  the transformed data matrix, preceded by  the
   number of rows and columns respectively,
2) if option 2 is specified:  the corrrelation  matrix  of  the  variables,
   preceded by its order,
3) the number of submodels specified in the list appended to option  5,  if
   that  option  is  specified  at  all;   otherwise  the  number 1.  It is
   followed by (for each (sub) model):  the number of estimated parameters,
   the estimates for the parameters of the (sub) model, and:
   a) if option 2 is specified:   the  variance-covariance  matrix  of  the
      estimates,  preceded  by  its  order,  (BE  CAREFUL:  this is NOT the
      correlation matrix of the estimates, which is printed;  however,  the
      correspondence  between  the  two  matrices  is  established  by  the
      relations (10) and (11) in "Help"/Theory),
   b) if option 3 is specified:  the number of respondents, followed by for
      each  respondent the:  observation, fitted value, standard deviation,
      residual, standardized residual and studentized residual,
and finishes by writing an "Eor" keyword.

     As in the case of printed output, the output described in 3)  is  only
effected  for submodels, if an explicit submodel specifier list is appended
to option 5.

     An input specification to describe one record of data written  to  the
outputstream  when  options  1, 2, 3  and 5 (with a submodel specifier list
appended to it) are specified, could read:

"Input"  n, m, n * (m * [transformed data element]),
         p, <p * (p+1) : 2> * [correlation element],
         s, s * (t, t * [estimate],
                 q, <q * (q+1) : 2> * [covariance element],
                 r, r * (6 * [residual element]) );
                                                                    Page 28


     For the original model the following relations  hold: q = t,  t = m-1,
r = n  and p = m (or p = m-1 if replications and/or weights are specified);
s is the number of processed (sub)models; for each submodel  t  and  q  are
decreased  with  the  number  of  terms  that are omitted from the original
model.

     Real numbers in the printed output are given  in  fixed  point  format
with  a  six decimal fractional part, the only exceptions are the estimates
for the regression parameters with their standard deviations, which have  a
ten  decimal  fractional  part and the numbers in the listings of the input
data and the transformed data matrix, which have a three decimal fractional
part.   Real  numbers in the data output are given in floating point format
with a sixteen decimal mantissa and a two decimal exponent part.


3.3  ERROR MESSAGES

Error messages against syntax or semantics have the following layout:  

                Error : <error text> or <error number>

The error text corresponding to the error numbers is:

 1   No input data given.
 2   All input data has been skipped.
 3   Attempt to process more input data than provided.
 4   Number in the input data is incorrect or too large.
 5   In a number '.' is not followed by a digit.
 6   In a number '#' is not followed by '+', '-' or a digit.

10   No model formula given.
11   Left hand part is not followed by '='.
12   Expression is not followed by ')'.
13   Option name used in a primary in an expression.
14   Incorrect primary in a factor in an expression.
15   Incorrect (control) identifier in an expression.
16   Parameter list of a standard function is not followed by ')'.
17   Standard function call with incorrect number of parameters.
                                                                    Page 29


20   No input formula given.
21   Expression in a control is not followed by '>'.
22   Option name used in a control in an input statement.
23   Input statement in a description is not followed by ')'.
24   Variable list in a description is not followed by ']'.
25   Incorrect description in an input statement.
26   Incorrect identifier in a variable list.
27   Item in a variable list is not an identifier.

30   Incorrect option number in an option statement.
31   Incorrect option name in an option statement.
32   Specifier list is not followed by ')'.
33   Number in a specifier list is incorrect or too large.
34   Specifier list is appended to incorrect option.
35   Specifier is not a number.
36   Specification is not properly continued.
37   Specification is not terminated with ';'.

40   No defined (independent) identifier to the right of '='.
41   Incorrect use of a parameter in a regression term.
42   Undefined (weight) identifier to the left of '='.
43   Undefined (dependent) identifier to the left of '='.
44   Number in a regression term is incorrect or too large.
45   Term does not have the form: param * factor or factor * param.
46   Undefined (independent) identifier in a regression term.
47   No regression parameter in a regression term.

50   Division by zero.
51   Integer division by zero.
52   Observation for dependent variable is in absolute value too large.
53   Observation for independent variable is in absolute value too large.
54   Exponentiation with zero base and non positive exponent.
55   Exponentiation with negative base and real exponent.
56   Weight factor is not positive.

60   Argument of  'Sqrt'  is negative.
61   Argument of   'Ln'   is not positive.
62   Argument of   'Log'  is not positive.
                                                                    Page 30


63   Argument of   'Exp'  is too large.
64   Argument of 'Arcsin' is in absolute value larger than one.
65   Argument of 'Arccos' is in absolute value larger than one.
66   Argument of  'Sinh'  is in absolute value too large.
67   Argument of  'Cosh'  is in absolute value too large.

70   Number of observations for the first dependent variable is zero.
71   Numbers of observations for the dependent variables are not equal.
72   Number of observations for the first independent variable is zero.
73   Numbers of observations for the independent variables are not equal.
74   Control reads an incorrect number in the input data.
75   Numbers of replications for the dependent variables are not equal.
76   Given, read or computed replication factor is not integral.

If the error number lies between:
 5 and 37, it is followed by the most recently processed identifier, number
           and  symbol.   Only  the first eight characters of each name are
           displayed.
41 and 47, it is followed by the number of the right hand  part  regression
           term  which causes the error, or a zero if the left hand part is
           at fault.
50 and 67, it is followed by the wrong value and the number of the line  in
           the  transformed data matrix which causes the error.  Instead of
           the wrong value, the number of the right  hand  part  regression
           term  which  causes the error is displayed when the error number
           lies between 50 and 53.
70 and 76, it is followed by the check value and the wrong value.   Instead
           of the wrong value, the value of the controlling variable of the
           next enclosing repetition  loop  is  displayed  when  the  error
           number is 76.


3.4  EXAMPLES

     An impression of the printed output of the  Multiple Linear Regression
Analysis program may be obtained by actually submitting the user program in
section 2.6 (which resides on:  Hlp: Mulexa.hlp) to the MULREG program.
                                                                    Page 31


                                APPENDIX 1

                             TECHNICAL REMARKS

     The following technical remarks reflect the  SIMULA implementation  of
the Multiple Linear Regression Analysis program as of version 5H(246).

     Any comments or queries concerning the  functioning  of  the  software
described in this document should be addressed to:  
      Marten van Gelderen, IKO Computer Systems Group, Postbox 4395,
    1009 AJ  Amsterdam, The Netherlands.  (telephone: 31-(0)20-930951).

     The program resides on a device called USR:.  It is started by typing:
.R MULREG,  and  responds to the standard output device (usually TTY:) with
an identifying header and requests to specify four files, as follows:

Multiple Linear Regression Analysis

Enter file specifications
Inputstream  :
Printstream  :
Datastream   :
Outputstream :

     The inputstream serves to read the user program from;  the printstream
receives  the  printed  output  from the program;  the datastream serves to
read the separate input data records from;  the outputstream  receives  the
data output records from the program.

     If the default carriage return is responded to the data-  and  output-
stream  requests,  the  program  assumes  that  no  separate input data are
present and that no data output is required.

     If the default carriage return is responded to the input-  and  print-
stream requests, the program connects the standard input and output devices
(usually TTY:) to the input- and printstream respectively.  To notify  that
the inputstream is connected to the standard input device, the program dis-
plays the prompting character '*' (asterisk).
                                                                    Page 32


     If both the input- and printstream are connected to the standard input
and  output  devices,  the  program  echoes  every text it cannot interpret
properly, preceeded by the error character '?' (question  mark).   However,
in response to a single carriage return the text 'For help type: "Help"' is
displayed.

     To the "Run" keyword a file specification list may be  appended,  pre-
ceeded by the character '/' (slash), with the following general format:

               /print-spec;output-spec=input-spec;data-spec

Each of the specifications will be connected to the  corresponding  program
stream  respectively.   Also,  each  of  the specifications may be omitted,
which means that nothing will be changed to  the  corresponding  stream  at
all.   However,  if  the  character  ' ' (blank) is substituted for one (or
more) of the specifications, the defaults - as described previously -  will
be connected to the corresponding streams respectively.

     If the specified files do not exist (for  input-  or  datastreams)  or
cannot  be created (for print- or outputstreams), the corresponding streams
are displayed again, followed by the  character  '?'  (question  mark),  to
indicate  the  erroneous situation and to enable the specification of other
files (or defaults).

     If the program encounters a premature  end-of-file  condition  in  the
inputstream,  it  will connect the inputstream to the standard input device
and thus respond with the prompting character.  New specifications for  the
"Model",  "Input",  "Options"  or "Data" (or keywords like "Run" or "Exit")
may then be entered.  

     In the inputstream, the program does not  discriminate  between  upper
and  lower  case  letters.   Identifiers  may contain any number of blanks,
however, a carriage return is not permitted, restricting the maximum length
of  identifiers  to the maximum number of characters in one input line.  In
front of and following the opening  quote  of  keywords  only  non-printing
ASCII  characters  like  tabs  and/or  blanks  are permitted, otherwise the
keyword (and the whole line following it) is not recognized.
                                                                    Page 33


     Two implementation dependent restrictions are  imposed  on  user  pro-
grams:   the  maximum number of differently spelled identifiers and numbers
is 789 and the maximum number of nested parentheses is 62.  In addition two
machine   dependent  restrictions  are  imposed:   the  maximum  number  of
characters in one input line is 132 and the number of significant digits in
computations is 18.

     In response to the keyword "Help", the information  in  this  appendix
(which  resides  on  HLP:)  is  copied  to  the  printstream.   Each of the
following switches may be appended to the "Help" keyword, preceeded by  the
character  '/' (slash), in order to obtain more detailed information (which
also resides on HLP:).

     /Theory   Regression model & least squares (section 1.1 & 1.2),
     /Tests    Possible tests of hypotheses (section 1.3),
     /Model    Specification of the model formula (section 2.1),
     /Input    Specification of the input formula (section 2.2),
     /Options  Possible options and their effects (section 2.3),
     /Data     Acceptable numbers and their delimiters (section 2.4),
     /User     Setup of a complete user program (section 2.5),
     /Example  Example of a complete user program (section 2.6),
     /Print    Standard and optional printed output (section 3.1),
     /Output   Standard and optional data output (section 3.2),
     /Errors   Meaning of the error numbers (section 3.3),
     /Syntax   Definition of the syntax of a user program (appendix 2).

If no help information is available an appropriate message is displayed.
                                                                    Page 34


                                APPENDIX 2

                DEFINITION OF THE SYNTAX OF A USER PROGRAM


     The syntax of a user program is defined in an extended  version  of  a
notation  known  as  the  Backus Naur Form, in short:  BNF (cf. NAUR [10]).
The extensions comprise an explicit repetition  and  optionality  construct
together with the possibility of factorization.

     The BNF may be regarded as a metalanguage for  the  description  of  a
user  program.   In  addition  to the symbols that are admissible in a user
program, the metalanguage  requires  a  number  of  extra  symbols,  called
metasymbols.   The ten metasymbols used in extended BNF are:  ::=, |, <, >,
{, }, [, ], ( and ).  The , and . are part of the metalanguage  English  in
which we are describing BNF.  We write:

          <expression>::= ['+' | '-'] <term> { ('+' | '-') <term> }

     The metasymbols < and > are used as delimiters to enclose the name  of
a class.  The metasymbol ::= may be read as 'is defined as' or as 'consists
of'.  The metasymbol | is read as 'or'.  Repetition  is  denoted  by  curly
brackets,  i.e. { a } stands for e | a | aa | ...  Optionality is expressed
by square brackets, i.e.  [ a ] stands for e | a.  Parentheses merely serve
for  grouping (factorization) i.e.  (a | b) c stands for ab | ac.  Terminal
symbols appear enclosed in single apostrophes.

     The above phrase defines an expression as a term, optionally preceeded
by  a  '+'  or a '-' and followed by an arbitrary repetition of terms, each
preceeded by a '+' or a '-'.

The syntax of a user program can thus be defined as follows:
                                                                    Page 35


<letter>::= 'A'|'B'|'C'|'D'|'E'|'F'|'G'|'H'|'I'|'J'|'K'|'L'|'M'|
            'N'|'O'|'P'|'Q'|'R'|'S'|'T'|'U'|'V'|'W'|'X'|'Y'|'Z'|
            'a'|'b'|'c'|'d'|'e'|'f'|'g'|'h'|'i'|'j'|'k'|'l'|'m'|
            'n'|'o'|'p'|'q'|'r'|'s'|'t'|'u'|'v'|'w'|'x'|'y'|'z'

<digit>::= '0' | '1' | '2' | '3' | '4' | '5' | '6' | '7' | '8' | '9'

<model keyword>::= '"Model"' | '"MO"'
<input keyword>::= '"Input"' | '"IN"'
<option keyword>::= '"Options"' | '"OP"'
<data keyword>::= '"Data"' | '"DA"'
<run keyword>::= '"Run"' | '"RU"'
<exit keyword>::= '"Exit"' | '"EX"'

<function name>::= 'Abs' | 'Sign' | 'Sqrt' | 'Sin' | 'Cos' | 'Tan' |
                   'Ln' | 'Log' | 'Exp' | 'Entier' | 'Round' | 'Mod' |
                   'Min' | 'Max' | 'Arcsin' | 'Arccos' | 'Arctan' |
                   'Sinh' | 'Cosh' | 'Tanh' | 'Indicator'

<option name>::= 'Transformed data matrix' | 'Correlation matrix' |
                 'Residual analysis' | 'No regression analysis' |
                 'Process submodels' | 'Print input data' |
                 'No input data rewind' | 'Save original model' |
                 'Test reduced model' | 'Missing values'

<option number>::= '1' | '2' | '3' | '4' | '5' |
                   '6' | '7' | '8' | '9' | '10'

<number>::= ['+' | '-'] <unsigned number>
<unsigned number>::= <decimal number> | <exponent part> |
                     <decimal number> <exponent part>
<decimal number>::= <unsigned integer> | <fractional part> |
                    <unsigned integer> <fractional part>
<exponent part>::= '#' <integer>
<fractional part>::= '.' <unsigned integer>
<integer>::= ['+' | '-'] <unsigned integer>
<unsigned integer>::= <digit> { <digit> }
                                                                    Page 36


<identifier>::= <letter> { <letter> | <digit> }

<data specification>::= <data keyword> [ <input data> ]
<input data>::= <number> { <number> }

<option specification>::= <option keyword> [ <option statement> ] ';'
<option statement>::= <option> { ',' <option> }
<option>::= <simple option> [ '(' <specifier list> ')' ]
<simple option>::= <option name> | <option number>
<specifier list>::= <specifier> { ',' <specifier> }
<specifier>::= <number>

<input specification>::= <input keyword> <input statement> ';'
<input statement>::= <input part> { ',' <input part> }
<input part>::= <control> | <description> | <control> '*' <description>
<control>::= <number> | <identifier> | '<' <expression> '>'
<description>::= '(' <input statement> ')' | '[' <variable list> ']'
<variable list>::= <variable> { ',' <variable> }
<variable>::= <identifier>

<model specification>::= <model keyword> <model statement> ';'
<model statement>::= <left hand part> '=' <right hand part>
<left hand part>::= <expression> [ '&' <weight part> ]
<weight part>::= <expression>
<right hand part>::= ['+'] <term> { '+' <term> }
<expression>::= ['+' | '-'] <term> { ('+' | '-') <term> }
<term>::= <factor> { ('*' | '/' | ':') <factor> }
<factor>::= <primary> { '^' <primary> }
<primary>::= <unsigned number> | <identifier> |
             <function designator> | '(' <expression> ')'
<function designator>::= <function name> [ '(' <parameter list> ')' ]
<parameter list>::= <parameter> { ',' <parameter> }
<parameter>::= <expression>

<user program>::= { <job> } <exit keyword>
<job>::= { <specification> } <run keyword>
<specification>::= <model specification> | <input specification> |
                   <option specification> | <data specification>
                                                                    Page 37


                                APPENDIX 3

                   TECHNICAL DESCRIPTION OF THE PROGRAM


     In this appendix a more or less technical description of  the  program
is given in terms of variables, procedures and control flow.

Basically the program logic is as follows:

                 Init program
                 Enter & open files
            Job: Read job
            Run: Init compiler tables
                 Compile model & input
                 Init data buffers
                 Execute (to produce design matrix)
                 Regression analysis
                 Print results
                 GOTO Job
           Exit: Close files

In case of an errorsituation in one of the sections, no further  action  is
undertaken and control is transferred to Job.  

1.  The program is initialized  by  setting  various  declared  system  and
compiler  constants  to  their appropriate values.  Some of them are highly
machine dependent constants, others are internal codings or  table  limits.
The  four  basic  informationstreams  in  the program:  the input-, print-,
data- and outputstream are connected to the external  file  specifications,
entered by the user, in the procedure enterfiles.

2.  A job is  read  specification  by  specification  (after  printing  and
skipping  leading text by a call on echotext) by means of calls on readtext
or readdata, depending on which specification from the job is  read.   Both
procedures  first  set  up  a buffer administration as will be described in
section 3, then process an input line or input numbers and finally call  on
readline  to  obtain  further information.  Their task is finished when the
                                                                    Page 38


variable endtext becomes true, that is when the next input line starts with
a '"' (quote), which is the beginning of what can be a keyword, or when the
end-of-file condition is met.  If an input data record is to be read from a
file, the same administration is setup and the same routine for reading the
actual numbers is used, only information is obtained  from  the  datastream
rather  than  the inputstream.  The echoing of input text and input data is
done via calls on the procedures printtext and printdata respectively.

3.  The text and data storage section provides two kinds of SIMULA classes,
one  for  text storage and another for data storage.  The class textstorage
provides a procedure  to  store  lines  of  text  into  a  linked  list  of
textbuffers  (starting  at  base).   The procedure nextline retrieves those
lines (after a reset).  The class datastorage provides a procedure to store
real  numbers  into a databuffer which is linked into a linked list of such
buffers (starting at base).  In case of buffer overflow  a  new  buffer  is
created  and linked into the list.  A super class inputstorage provides the
reset procedure and a procedure nextnumber to retrieve numbers out  of  the
buffers  in  a  FIFO  (first  in,  first  out) manner.  Another super class
rightstorage provides a procedure lastnumber, which  is  quite  similar  to
nextnumber,  except  that it retrieves numbers out of the buffers in a LIFO
(last in, first out) manner.

4.  The compiler uses two tables: one called  program  and  another  called
hashtable.   The  first table accepts the (macro) instructions generated by
the compiler via calls of load.  The second table is used by the  procedure
nextatom  to  identify  the  various items in the model and input formulae.
The recognition method used is twin-prime-hashing described  in  KNUTH  [5]
p. 522.   Alphanumerical  items can belong to one of the following classes:
identifier, number, functionname or optionname, which are all  superclasses
of  atomcell.   The  procedure  nextatom  calls  on nextchar and delivers a
pointer to an atomcell of the appropriate type  containing,  among  others,
the actual text of the item and its index in the hashtable.

     The compiler assumes that the running system has  at  its  disposal  a
(programmed  pseudo)  register  F which is capable of handling both integer
and real numbers.  Furthermore a (programmed  pseudo)  memory  organization
known  as  a stack must be available to the running system.  A stackpointer
refers to the first free position in the stack.  All binary operations will
                                                                    Page 39


take  place with the top of the stack as first operand and F as the second.
The result is delivered in F and as  a  side  effect  the  stackpointer  is
decreased  by  1.   When  the  contents  of  F  is  saved in the stack, the
stackpointer is increased by 1.

     The fundamental idea behind most procedures for translating the  model
and  input  formulae  is, that the first atom of the syntactical unit to be
processed by that procedure  has  been  read  already  (its  'value'  being
assigned to lastatom).  The procedure considers itself to have finished its
task after reading the first atom that no longer can belong  to  that  unit
syntactically.   Meanwhile  the translation of that unit has been produced.
A  more  elaborate  description  of  the  procedure   system   (arithmetic)
expression  can  be  found  in  KRUSEMAN  ARETZ [6] and [7].  Reference [8]
provides the description of a complete ALGOL 60 compiler.  We only  mention
here  that  every  expression  is  transformed  into  a  macro program that
corresponds to the reversed polish form, thus:

               (a+b) * (c-d) ^ e  becomes:  ab+ cd- e ^ *.

The procedure  for  translating  the  input  formula  must,  among  others,
generate  instructions  to perform the linkage between identifiers from the
model formula and numbers from the input data.

     While translating a model formula, identifiers to  the  right  of  the
equal sign are assigned type 1, those to the left of the equal sign type 2.
If these identifiers appear in a variable list the types are changed into 3
and  5 respectively.  While translating an input specification, identifiers
in a variable list not appearing in the model formula are assigned type  4,
those  in  the  input formula not appearing in a variable list are assigned
type 6.  Meanwhile instructions are generated to put the next  number  from
the  input data in the appropriate column of the (yet untransformed) design
matrix, or to skip that number.  For variable lists that consists  entirely
of  identifiers not appearing in the model formula, special instructions to
skip the corresponding numbers all in one, are generated.

5.  In the procedure check model a check is  made  if  the  model  formula,
after  the linkage to the input data (by means of the input formula), still
satisfies some elementary statistical conditions, like:
                                                                    Page 40


a) each term must be the product of a parameter and a factor,
b) in that factor no identifier  may  appear  that  is  not  present  in  a
   variable  list  in the input formula.  (An attempt to perform regression
   analysis with variables for which no  input  data  is  present  may  not
   succeed.)

     In the procedure check input a check is made if for each  variable  in
the  model  formula an equal amount of numbers is present, moreover a check
is made if all numbers in the  input  data  have  actually  been  processed
(option 7 disengages this check).

6.  The execution section of the program is activated  by  a  call  of  the
procedure execute which, among other things, simulates the basic cycle of a
computer:

         next: get the instruction indicated by the programcounter
               increase the programcounter with 1
               isolate the instruction and address part
               execute the instruction
               GOTO next

This cycle ends when the programcounter  tries  to  leave  the  (translated
macro)  program.   The (macro) instructions itself are coded via the switch
lists macro and macro2.

7.   After  the  execution  of  the  input-  and   modelinstructions,   the
(transformed) design matrix is delivered to the regression routine(s).  The
actual computation of the regression coefficients is done  via  a  call  of
lsqdec  followed  by  a  call of lsqsol and lsqinv.  The first two of these
procedures are described extensively in DEKKER [2] pp. 65-69.   The  vector
and  matrix  multiplication is done via calls of vecvec, matvec, tamvec and
tammat, described in DEKKER [2] pp. 8-9.  The algorithms for phi and Fisher
are copied from CACM: algorithms 209 and 322 respectively.

All other computations are straightforward.
                                                                    Page 41


REFERENCES


 [1]  AFIFI, A.A. & S.P. AZEN, Statistical Analysis;  A  Computer  Oriented
            Approach;  Academic Press, (1972).

 [2]  DEKKER, T.J., ALGOL 60 procedures in numerical algebra, part  1;   MC
            Tract 22, Mathematisch Centrum, Amsterdam.

 [3]  DRAPER, N.R. & H. SMITH, Applied Regression Analysis;  John  Wiley  &
            Sons, (1966).

 [4]  JONGE, H. DE,  Inleiding  tot  de  Medische  Statistiek,   deel   II;
            Nederlands Instituut voor Praeventieve Geneeskunde, (1960).

 [5]  KNUTH, D.E., The Art of Computer  Programming,  Vol. 3,  Sorting  and
            Searching;  Addison Wesley, (1973).

 [6]  KRUSEMAN ARETZ, F.E.J.,   ALGOL   60   translation   for   everybody;
            Elektronische Datenverarbeitung, Vol. 6 (1964), 6, p. 233-244.

 [7]  KRUSEMAN ARETZ, F.E.J.,  Programmeren  voor  rekenautomaten,  (De  MC
            ALGOL 60 vertaler voor de EL X8);  MC Syllabus 13, Mathematisch
            Centrum, (1972).

 [8]  KRUSEMAN ARETZ, F.E.J., P.J.W. TEN HAGEN &  H.L. OUDSHOORN, An  ALGOL
            60  compiler  in  ALGOL 60;  MC Tract 48, Mathematisch Centrum,
            (1973).

 [9]  LUND, R.E., Tables for an approximate test  for  outliers  in  linear
            models;  Technometrics, Vol. 17 (1975), 4, p. 473-476.

[10]  NAUR, P. (ed.), Revised report on the algorithmic language ALGOL  60;
            Regnecentralen, Copenhagen, (1964).

[11]  SEARLE, S.R., Linear Models;  John Wiley & Sons, (1971).

[12]  THEIL, H., Principles of Econometrics;  John Wiley & Sons, (1971).
                                                                    Page 42