Intro Overview
OVERVIEW

  This statistics helpfile has been developed for use as an aid to those wishing to learn about statistical analysis.  The experienced Windows user will find familiar procedures which will allow them to move through the helpfile with relative ease.

  New windows users or users with no experience of helpfiles can find more help and information about this helpfile by clicking on the beginners below.

BEGINNERS	CONTENTS


Intro Beginners Help
BEGINNERS HELP

  The helpfile is constructed as a series of topics or pages containing text or graphics explaining the statistical procedures covered.  Using the buttons displayed the user is taken through each of the steps that are needed to perform the statistical tests.

  'Continue' buttons take the user from step to step, further buttons such as 'Details', 'Minitab', 'Protocol' and 'Example' will expand the particular topic on view.  The 'Finished' button  will return the user to the first page concerned with that topic.

  A summary or Protocol of each of the procedures is provided and this can be printed and kept for future reference.


.$ABOUT WINDOWS HELPFILES

  When text is highlighted in green (with broken underlining) these 'hot-spots' can be clicked on to obtain either more information about a particular point or a definition.  A second click will remove the information.  Try clicking now.

  Detailed information about how to use Help can be found by clicking on Help on the menu bar.  (The print option, in File on the menu bar will print the page on view).

Here are some useful hints:-

1.   The maximise button [image] in the top right hand corner will maximise the screen this help file is best viewed at maximum size.

2.   The minimise button [image] will reduce the screen to an icon, double clicking on the icon will restore the screen.

3.   The scroll bar to the right of the screen allows the user to scroll up or down the page;  the arrow button at the top of the scroll box moves the screen down one line and the arrow at the bottom moves the screen up one line.

4.   Clicking anywhere on the bar below the position marker button (now at the top) will move the screen up by one page.  Clicking above this marker button will move the screen down by one page.

5.   The Contents button will return to the contents page.

6.   The Search button will search for a particular topic in the helpfile.

7.   The Back button will return to the last screen viewed.

8.   The History button displays the last 40 topics viewed. (Double click to re-visit).

9.   The Browse buttons (<< >>) move through an ordered set of topics.


Now click again to remove


Intro Contents
CONTENTS

Oneway analysis of variance (ANOVA) Overview Route through a one-way ANOVA Computing the one-way table Checking assumptions of the ANOVA Transformations Testing for outliers Comparisons of means	Factorial analysis of variance (ANOVA) Introduction Route through a factorial ANOVA Computing the factorial table Checking assumptions of the ANOVA Transformations Testing for ouliers Comparisons of means
Bivariate Regression Introduction Route through a bivariate regression Checking assumptions of bivariate regression Transformations Computing the regression Residual analysis Influence statistics Robust linear regression Regression through the origin Regession with more than one y for each x Curved regression lines	Multiple Regression Introduction Route through a multiple regression Checking assuptions of multiple regression Transfomations Computing the regression Residual analysis Influence statistics Multicollinearity Selection of variables Calculation of confidence or prediction intervals


Oneway ANOVA
ONE-WAY ANALYSIS OF VARIANCE(ANOVA)

  This is an extremely useful statistical approach in all branches of biology.  It can be used to test whether there are any significant differences between a set of means from an experiment or survey and to identify where those differences lie.

A typical data structure:-

Definitions	a-number of groups in the factor being tested
n-number of replicates
x-represents an individual factor (one group indicated by a box)


Oneway EXAMPLE DATA SET
EXAMPLE DATA SET

  Four replicate samples were taken from a large river in Northern England at each of six sites.  The numbers of Escherichia coli (per 100ml) were estimated in each sample by a membrane filtration method.  The following results were obtained.

A major tributary joined the river between sites 4 and 5.

  An analysis of this data allows any differences in E. coli numbers at different sites above and below the tributary to be investigated.

  This example will be referred to throughout one-way analysis of variance.


Oneway ROUTE
ROUTE THROUGH A ONE-WAY ANOVA

1.	Enter the data

2.	Compute the one-way ANOVA table

3.	Check assumptions of ANOVA

4.	Transform data if necessary

5.	Examine data for outliers

6.	Compare group means

7.	What to do if the assumptions do not hold

For more details see One-way Protocol


Oneway ENTERING DATA
1. ENTERING THE DATA

  Most packages require all data to be in a single column or variable.  These data are then labelled with another set of whole, integer numbers in another column or variable.  This column is called the subsript, labelling, level or category variable.


Oneway EXAMPLE DATA ENTRY
EXAMPLE DATA ENTRY

	RAW	SITES
The file is organised as 25 rows, the top row contains variable labels for the data (RAW; number of E. coli per ml) and a labellng variable (SITES;1-6) for the six sites sampled.  The top row may have to be removed before the file can be read into some packages.	7 4 0 5 61 13 0 42 50 155 106 100 11 6 2 13 0 0 0 1 1 4 1 0	1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6


Oneway PROTOCOL1
ONE-WAY ANOVA	PROTOCOL 1

Route through a one-way ANOVA*

1.	Compute one-way ANOVA table.

2.	Check assumptions of ANOVA

(a)	Tests for homogeneity of variances

 Bartlett's

 Hartley's Fmax

 Cochran's C

(b)	Tests for normality of residuals

 normal probability plots

 Shapiro-Wilk test

 boxplots

3.	Transform data if assumptions do not hold (see Transfomation Protocol ) and repeat steps 1 and 2.  (If assumptions do not hold after transformation, firstly check step 4, to see if outliers are responsible for some large variances, and secondly follow Trouble Protocol ).

4.	Examine residuals of transformed data and test for outliers, consider removing outliers from groups if there is justification.  Tests for outliers are

 Dixon's test (n < 25)

 Grubb's test (n > 25)

5.	Comparisons of group means (see Trouble Protocol  if assumptions do not hold).

(a)	A Priori tests (could be planned before results obtained)
 calculate LSD for pairs of means
 orthogonal decomposition for more than two means; sometimes called contrasts.

(b)	A Posteriori tests (could not be planned before results obtained)
 calculate MSD by Tukey-Kramer method for pairs of means
 use Gabriel's SS-STP procedure for more than two means

* Abbreviations are: ANOVA, analysis of variance;  LSD, least significant difference;  MSD, minimum significant difference;  SS-STP, sum of squares simultaneous test procedure.


Oneway ANOVA TABLE
2. COMPUTING THE ONE-WAY ANOVA TABLE

  The next step is to compute the basic ANOVA table by using the appropriate command from the package being used.  It is also useful to save the residuals at this stage for use later.  The residuals are our best estimate of the error in the analysis and are calculated by

where x are the data values and [image] are the groups means.

Click to view Minitab commands and output


Oneway MTB ONEWAY
MINITAB COMMANDS AND OUTPUT	page 1 of 3

The Minitab command for the one-way table is:-

ONEWAY  C1  C2  C3  C4

C1  is the variable containing the data

C2  is the variable containing the labels

C3  and  C4  (empty) are where Minitab stores the residuals and group means.

The output will appear

ANALYSIS OF VARIANCE ON....

SOURCE	DF	SS	MS	F	P
LEVEL			ms
ERROR			ms
TOTAL
Minitab also produces a table containing the means and standard deviations of the groups.

Click for a more detailed explanation and an example for you to try


Oneway TABLE DETAILS
ONE-WAY ANOVA TABLE DETAILS 	page 2 of 3

  The ANOVA table from the E. coli data shows important features which should be noted.  The mean squares for between groups (MSgroup) and within groups (MSerror) are calculated by dividing the appropriate sum of squares by the corresponding degrees of freedom.  These mean squares are measures of the variation between sites (MSgroup) and within the sites (MSerror).  The latter is an average variance over all the groups in the ANOVA and is hence an estimate of the error variance, so it is designated MSerror; it is sometimes also called the residual mean square.  The F-ratio is calculated as

and is compared with the table FP(a-1,dferror) to find the probability level at which Fcalc is significant.

  This is normally given by the package but it is useful to know that Fcalc should be large for the result to be significant.


ONEWAY FCALC PROBLEM
ONE-WAY ANOVA PROBLEM 	page 3 of 3

The one-way ANOVA table for the E. coli data would appear

SOURCE	DF	SS	MS	F	P
LEVEL	5	31810.000	6362.000	?	.0000
ERROR	18	7930.500	440.5833
TOTAL	23	39740.500

Since the equation is

the answer is?

142.700	0.069	14.440	0.000
  Thus for the E. coli data Fcalc should be large for the result to be significant at P<<0.001 as
 Fcalc>FP(a-1,dferror) because

F0.1(5,18)     =  2.20
F0.05(5,18)   =  2.77
F0.01(5,18)   =  4.25
F0.001(5,18) =  6.81

The result is therefore?

 SIGNIFICANT	NOT SIGNIFICANT

  This procedure is logical, a moment's thought will show that a high error variance results in a low
F-ratio and a non significant F-test.

  A significant F-test means that there are significant differences somewhere between the group means, providing the variances are homogeneous and the errors are normally distributed (step 2 Protocol 1).  It does not indicate which means are different from each other, this is examined by the tests for comparisons of means (step 5 Protocol 1).


Degrees of freedom


Sum of squares


Mean square(variance)=SS\DF  used to make decisions about the data


MSgroup this is an estimate of variance between the groups


MSerror this is an estimate of variance within the groups


F-ratio   calculated as MSgroup\MSerror


P-value   calculated by Minitab(less than 0.05 for significance)


Source of variation


WRONG!  Check your calculation


WRONG!  Try the other way round


CORRECT!  Now check this answer against F tables


WRONG!  Try again


CORRECT!


WRONG!  14.440 is larger than any of the F table values.  This means the result is signficant.


Oneway ASSUMPTIONS
3. CHECKING THE ASSUMPTIONS OF THE ANOVA

The main assumptions which underpin one-way ANOVA are

1.	Randomness of replicates within a group

2.	Independence of errors

3.	Homogeneity (equality) of group variances

4.	Normality of residuals(errors)

  The first two assumptions must be carefully considered when designing the experiment.  The remaining two assumptions are very important but can only be tested once the data is obtained.

  Homogeneity of variances can be tested directly and normality of errors by using the residuals as they are our best estimates of the real errors.


Oneway HOM OF VAR
HOMOGENEITY OF VARIANCES	page 1 of 3

  This is the most important assumption for ANOVA and must always be tested.  The three most useful tests for variance equality are Bartlett's [image], Hartley's Fmax and Cochran's C.  These are also the most commonly applied and available in many computer packages.  All three tests are calculated in SPSS/PC+ and Statgraphics, Systat only calculates Bartlett's [image], but provides an alternative F-statistic for small sample sizes.  Minitab does not have commands for any of these statistics.  However, a set of Minitab macros that will calculate all three statistics and store group means and variances, can be used.

MACROS

  Bartlett's test is best calculated by computer.  The two other tests can be calculated more easily from computer output if the group variances [image] or standard deviations (s) are provided (Minitab, Systat).


Oneway HARTLEYS+COCHRANS
HOMOGENEITY OF VARIANCE	page 2 of 3

Hartley's Fmax and Cochran's C

These tests are calculated by

where [image] is the largest group variance, [image] is the smallest group variance and [image] are all the group variances.


Oneway HARTLEYS+COCHRANS2
HOMOGENEITYOF VARIANCE	page 3 of 3

Hartley's Fmax and Cochran's C

  For the E. coli. data the group standard deviations are

 2.94, 27.6, 42.9, 4.97, 0.500 and 1.73 so

These values are compared with tabulated values at P = 0.05.

So with Fmax at a and n-1 degrees of freedom the tabulated value for Fmax(a,n-1) =Fmax(6,3)=62.0 and this is below the calculated Fmax so the variances are heterogenous.

  Similarly with the tabulated value for C(a,n-1)=C(6,3)=0.532 is smaller than the calculated value so the variances are also declared heterogenous.

  In all three of these tests calculated values smaller than the table values are required for the variances to be homogeneous and this assumption of ANOVA to be upheld.

  The Fmax and C tables are designed for equal group sizes.  If the data has unequal n  then the tables can still be used to give an approximate result.

  Care must be taken when interpreting the results of computer output with these tests.  The probability values often given (SPSS/PC+, Systat, Statgraphics) should be greater than P=0.05 for the variances to be significantly homogeneous.  Notice that for these tests the desirable outcome is for the P value to be large and not small as is the case for many other tests.  We can then be confident about proceeding with the one-way ANOVA assuming equal variances.


 For example,some data has n = 5-15 with Fmax = 148.8 and C = 0.651.  The tabulated values for highest and lowest n are Fmax(4,15) = 4.01, Fmax(4,5) = 13.7, C(4,16) = 0.437 (the nearest value in the table, to C(4,15)) and C(4,5) = 0.590.  Thus the calculated values are larger than all possibilities from the tables and so it is clear that the variances are once again declared heterogenous.  If the calculated value for either of these test was between the two tabulated values the result would be uncertain, but if less than the smallest tabulated value (that for the largest n) then the variances would be homogeneous.


Oneway NORM OF RESIDS
NORMALITY OF RESIDUALS	page 1 of 5

  One-way ANOVA and the tests for comparisons of means carried out subsequently
(step 6; Protocol 1) are not as sensitive to lack of normality of errors as they are to heterogeneity of variances.  However, all the tests for homogeneity of variances are susceptible to lack of normality and so can give misleading results if errors are badly non-normal.  For this reason it is essential to carry out some sort of test for this condition.

  Some of the best tests for normality are cumulative normal probability plots against the residuals.  Most statistical packages allow such plots to be prepared easily (Systat, Statgraphics) with simple plotting commands.


Oneway MTB NSCORE
MINITAB COMMANDS AND OUTPUT	page 2 of 5

  Minitab allows the calculation of normal scores with the NSCORE command.  Plotting the residuals against the normal scores gives the usual cumulative normal probability plot.  These plots will be straight lines if the residuals are normally distributed and curved if not.

  The Minitab commands are

NSCORE  C3  C5

C3 contains the stored residuals from the ONE-WAY command

C5 (empty) is where Minitab stores the normal scores

PLOT  C3  C5


These are a set of expected points from a normal distribution with a mean of zero and one value for each of the residuals.


Oneway MTB PLOT
MINITAB COMMANDS AND OUTPUT	page 3 of 5

The plotted data will appear:-

nscores-

          -                                                       *

          -

     1.5  +                                             *

          -                                    *

          -                               **

          -                              *2

          -                             2*

     0.0  +                             3

          -                             2

          -                            3

          -                          **

          -                     *

    -1.5  +               *

          -

          -   *

          -

            --------+---------+---------+---------+---------+--------raw residuals

                  -40       -20         0        20        40

The residuals from the E. coli data are clearly not on a straight line.Thus the data needs to be transformed.


Oneway MTB CORR
MINITAB COMMANDS AND OUTPUT	page 4 of 5

  Although the linearity of the probability plot can be assessed by eye a specific test for linearity is useful.  Correlation of the normal scores and residuals gives a coefficient which can be compared with table values to see if the hypothesis of normality is acceptable.  The Minitab command is

CORR  C3  C5

C3  are the residuals

C5  are the nscores

  The correlation coefficient from the E. coli data is 0.871.  The table value for N (total number of data points in the one-way ANOVA = 24) is 0.956.  As a high calculated value indicates normality the raw data is very clearly non-normal. (Failure by a narrow margin would lead to very few errors in the final stages of the analysis of variance and so would probably be an acceptable deviation from normality).  This test is essentially the same as the powerful Shapiro-Wilk test for normality.


Oneway MTB BOXPLOT
MINITAB COMMANDS AND OUTPUT	page 5 of 5

  Another straightforward graphical way to examine residuals is by using the Box-and-Whiskers plot or boxplot.  This is really a graphical representation of a letter value display.  The Minitab commands are

BOXPLOT  C1

or

GBOXPLOT  C1      (for high resolution graphics plot)

C1  contains the data to be plotted (raw or transformed)

The output will appear
  The boxplot shows the median (middle data point), and hinges (middle points between the median and maximum and minimum values), adjacent values (the end of the whiskers) and the more extreme values outside these values (outside and far outside values).  Thus residuals that are evenly distributed will show symmetrical boxplots and uneven residuals will give non-symmetrical boxplots.  As the normal distribution is even, normally distributed residuals will give an even boxplot with a central median and slightly longer whiskers than the median-hinge distance.  Boxplots are part of a useful, wider group of statistical procedures called exploratory data analysis.(see Transformations for more details)


Oneway TRANSFORMATIONS
4. TRANSFORMATIONS

  If analysis of the raw data shows that the variances are not homogeneous and the residuals are not normally distributed then none of the results of the one-way ANOVA can be relied upon.

  The usual action is to search for a suitable transformation which will modify the scale upon which the data is expressed so that the two critical assumptions are valid.  Biologists often think that transformation somehow corrupts the original data, but this is not so.  It is a perfectly valid procedure and will often make understanding the data easier as well as enabling ANOVA to be carried out.

It is also possible that the critical assumptions might not be valid because of severe outliers within the data.  Hence this possibility should be considered whilst the data is being examined (step 5; Protocol 1).

  After transformation the critical assumptions must be examined again to check that the transformation has worked (step 3; Protocol 1).

HOW TO CHOOSE TRANSFORMATIONS


Oneway CHOOSING TRANSF
HOW TO CHOOSE AN APPROPRIATE TRANSFORMATION

The approaches outlined in steps 2-5 below are given in order of ease of computation with most computer packages.

1.	Select an appropriate method according to the data structure and ease of computation
 any data; steps 2 and 5
 several groups with small numbers of replicates;  steps 2 and 3
 at least one large group;  step 4

2.	Choose a transformation based on theory.

3.	Taylor's power law.

4.	Letter value display method.

5.	Minimization of homogeneity of variance statistics

For more details see Transformation Protocol

Note

  There is no need to use all the methods suggested here to choose a transformation for unfamiliar data sets.  It is best to select one of these methods and to base the choice on the arrangement of the data and the ease with which the computations can be done with the computer package being used.


Oneway TRASF PROTOCOL
TRANSFORMATION PROTOCOL

1.	Select an appropriate method according to the data structure and ease of computation
 any data; steps 2 and 5
 several groups with small numbers of replicates;  steps 2 and 3
 at least one large group;  step 4

2.	Choose a transformation based on theory, some examples follow
 population data;  log10 x or log10 (x+1)
 low counts in quadrats; [image] or
 percentages or proportions; the arcsine transformation,
    or

3.	Taylor's power law

(a) Calculate group means [image] and variances [image] for raw data.

(b) Plot log10 [image] against log10 [image] for raw data.

(c) Calculate slope (b) of the regression line from the above plot where

log10{bm ssq.bmp}  =  a + b log10

(d) Calculate appropriate power (p) for transformation where

(e) Use power directly to transform by [image] or choose nearest power transformation from the ladder of powers .

4.	Letter value display method

(a) Transform data according to ladder of powers .

(b) Compute midsummaries for all transformations.

(c) Calculate se from spreads/normal spreads for all transformations.

(d) The best transformation has the most even range of midsummaries and se values.

5.	Minimization of homogeneity of variance statistics.

(a) Transform all data according to ladder of powers

(b) Calculate Hartley's Fmax, Cochran's C and Bartlett's [image] statistics for all transformations.

(c) Use transformation giving the minimum value for most statistics.


Oneway BASED ON THEORY
CHOICES BASED ON THEORY OR EXPERIENCE

  Biologists who regularly analyse their data with correctly applied statistical procedures will soon find the best transformations for the different types of data they generate.  Thus the more time consuming methods of searching for a good transformation will rarely need to be routinely applied.

  The most commonly used are

Logarithmic Transformations

Square Root Transformations

The Arcsine Transformation


Oneway LOG TRANSF
LOGARITHMIC TRANSFORMATIONS

  A log10 transformation is often suitable for population data.  If the data contains zero values then a small constant number should be added to all data points.  This constant should represent the smallest number that could sensibly be found, that is the limit of detection, for example log10(x+1).Many population related variables can be sensibly transformed by logarithms.  Enzyme activities of bacteria, primary production in seawater and grazing rates of zooplankton are all examples of this.  Data collected from survival experiments often naturally follow a logarithmic decline and hence can be log10 transformed.

  Using log10 the base of logarithms is easy to understand as the number of orders of magnitude; 1 million insects is [image] or 6.0 when transformed.  However, logarithms to any base can be used and some may prefer loge.  These two forms can be easily interconverted in statistical packages without one or the other transformations as log10x is equal to logex multiplied by loge10.


Oneway SQ ROOT TRANSF
SQUARE ROOT TRANSFORMATIONS

  Data collected in quadrats are distributed according to the Poisson distribution.  Haemocytometer counts of blood cells, counts of bacteria in petri dishes and plants in quadrats are examples of this type of data collection.  When the counts are low (<20) the distributions are very skewed and are best transformed by the square root of x or x+1.  When counts are high (>20) the Poisson and Normal distributions are so similar that transformation might not be necessary.


Oneway ARCSINE TRANSF
THE ARCSINE TRANSFORMATION

  This is a special transformation which is for proportions or percentages and can be applied almost automatically to data of this type.  The transformation for proportions.

and has been designed to make the variances homogeneous.  With percentages the values should be converted to proportions by division by 100.  It can be looked up in tables but is most conveniently calculated in degrees by computer.  The Minitab command for this conversion assumes the raw data are percentages in column 2 (C2) and the transformed data are to be stored in column 3 (C3).

LET C3 = (ASIN(SQRT(C2/100))) * 57.296

Other packages make use of similar commands.


Oneway TAYLOR'S POWER LAW
TAYLORS POWER LAW	page 1 of 2

  This provides a very convenient method for finding the best transformation for a set of data.  It works best when there are several groups (> 6) and is effective even when the group size is small.

  The method (Transformation Protocol , step 3) works by calculating the slope (b) of the regression line of the log10 of the group variances upon log10 group means.  The power (p) for the transformation [image] is then calculated by p = 1 - b/2.

  For the resulting power transformation to be effective the plot of log10 [image] versus log10 [image] should be linear.  That is all the points should lie close to a straight line.

  The slope from this regression of the Ecoli data was 1.56 hence the power suggested is
1 - 1.56/2 = 0.22 and would seem reliable as the points on the corresponding graph were close to a straight line.  Thus the transformation suggested is


Oneway TAYLOR'S POWER LAW2
TAYLORS POWER LAW	page 2 of 2

  The power transformations indicated by the Taylor's power law method will be slightly different for each new set of data collected and so do not provide continuity of statistical analysis within a research programme.  It would be better to use the calculated power as an indication of the best transformation to use.  The ladder of powers provides a logical way of doing this.

	Transformation		Addition of small  number needed for zero values
Power	Name	Formula
3 2 1 0.5 0 -0.5 -1 -2	cube square raw data square root logaritm reciprocal root reciprocal reciprocal square	x log10x 1/ 1/x 1/	no no no yes yes yes yes

Based on a table in Velleman and Hoaglin

a Not needed but [image] or [image] is often used.

  The ladder puts the raw data and a range of equivalent transformations into a logical sequence from which a suitable transformation can be selected from the power suggested from Taylor's power law.

  The power from the E. coli data was 0.22 which is midway between the log10 and the square root transformation.  This uncertainty is perhaps to be expected because the counts are essentially population data collected from samples with many low counts.


Reference

Velleman, P.F. and Hoaglin, D.C. (1981). Applications, Basics, and Computing of Exploratory Data Analysis. Duxberry Press, Boston


Oneway LETTER VALUE DISPLAYS1
LETTER VALUE DISPLAYS	page 1 of 2

This is another technique of exploratory data analysis and Minitab will allow calculation of the associated statistics.  This is a diagrammatic representation of how the letter values are assigned and the midsummaries and spreads calculated.

  Briefly the calculations are done as follows.

  When the data is ordered the median is the middle data point and is assigned the letter value M.  The two middle data points between the median and the lowest and highest values are called hinges and assigned H.  This principle continues through eighths (E) and letter values D,C,B,A,Z etc. which do not have names.  When there is not a single central value the central two data points are averaged to calculate the number required.  The midsummaries are then calculated from averages of the pairs of letter values; hence, for example, mid H = (upper H + lower H)/2.

  For a perfectly symmetrical distribution these midsummaries will all be equal but for a distribution skewed to the left they will increase and will decrease for one skewed right.


Oneway LETTER VALUE DISPLAYS2
LETTER VALUE DISPLAYS	page 2 of 2

  Spreads are calculated as the differences between the letter values and so are the magnitudes of ranges between the letter values.  Hence, for a normal distribution they will be entirely predictable.  The values of the spreads from a standard normal distribution (mean=0, standard deviation=1) are called
normal spreads.  Thus the relationship between a normal distribution and the data being analysed can be examined from the ratios of the calculated spreads and the normal spreads.  This ratio is an estimate of the standard deviation of the data (se) so

and will be constant if the distribution of the data is close to normal.

  As with the midsummaries the normal spreads will increase or decrease for skewed distributions.

  These ideas have been used to make a test to select the best transformation for a set of data.

  Click for details


Oneway LADDER TEST
TEST TO SELECT THE BEST TRANSFORMATION

  This test should be done on data from a single group in the one-way ANOVA and only works well for data sets with at least one large group (>ca.12).  The test is done by transforming the data from the largest group by all the transformations in the ladder of powers, taking care to add an appropriate constant (perhaps 1 or 0.5) to all the values if zeros are present.  The midsummaries and se values are then calculated and the best transformation is the one in which these values remain constant and do not increase or decrease in a regular manner.

  A Minitab macro (LADDER.MTB) can be used which calculates the midsummaries and se values.  Although this approach to selecting a transformation might not appear very precise it works remarkably well.


Spreads for letter values from a standard normal distribution

Letter value	Normal Spread
H E D C B A Z	1.349 2.301 3.068 3.726 4.308 4.836 5.320


Oneway MINIMIZATION
MINIMIZATION OF HOMOGENEITY OF VARIANCE STATISTICS

  The three statistics used to check data for homogeneity of variance (Bartlett's [image], Fmax, C) can also in principle be used to select a good transformation.  This can be done by calculating these statistics for data transformed according to the ladder of powers.  By comparing the critical values from tables at
P = 0.05 with the calculated values for the different transformations a suitable transformation can be chosen.


Oneway MINIMIZ TABLE
RESULTS FROM TESTS FOR HOMOGENEITY OF VARIANCE

Transformation	Bartlett's	Hartley's Fmax	Cochran's C
Square Raw Data Square Root Log10 Reciprocal Root Reciprocal TABLE VALUE	127.9 45.9 19.09 9.80 13.99 23.5 11.07	> 7374 218.4 29.0 263.6 7075 62.0	0.964 0.697 0.587 0.620 0.482 0.309 0.532

  Hence for the E. coli data Bartlett's [image] and Fmax both reach acceptable minima at the logarithmic transformation, but C is below the critical value for both types of reciprocal transformation.

Additionally Bartlett's [image] indicates that log10x and reciprocal are acceptable and C suggests reciprocal is satisfactory.  Clearly with different data sets these tests are likely to have different discriminatory powers.  Once again the overall conclusion is that log10x works best for the E. coli data.


Oneway TESTING FOR OUTLIERS
5. TESTING FOR OUTLIERS

  Sometimes data sets might contain one or more very extreme values that may not really be part of the same data.  Such values are outliers and can occur for a variety of reasons, such as transcription or experimental errors.

  Outliers can sometimes be seen when residuals are being tested for normality.  For example the probability plot of the log10 E. coli residuals shows one value with the lowest residual that might be an outlier.

  Inspection of a table of residuals shows that the suspect point from the E. coli data is the zero value at site 2.  This is the seventh point and can also be detected easily by plotting the residuals against the order of the data.  Minitab has a time series plot for this purpose.

The most useful test for the significance of outliers is DIXON'S TEST


Oneway MTB OUTLIERS
MINITAB COMMANDS AND OUTPUT

The Mintab command is  TSPLOT  C10   < contains the residuals from the transformed data

The output will appear

         -

     0.60+     5

         -        8

 resids  - 1                    2

         -    4     0  3  6   0

         -  2

     0.00+      6    12 4      1 3

         -                 789

         -         9

         -               5        4

         -

    -0.60+   3

         -

         -

         -

         -

    -1.20+       7

          +---------+---------+---------+

          0        10        20        30

This confirms that the seventh point is suspect.


OnewayDIXONS TEST1
TESTING FOR OUTLIERS	page 1 of 2

Dixon's Test

  The most useful test for the significance of outliers is Dixon's test which is designed for small group sizes (n<25).  However, Grubb's test is available for data sets with large group sizes (n>25).  To test the possible outlier in the E. coli data from site 2 the log10(x+1) transformed values must be used because Dixon's test only works with normally distributed data.The correct formula for a group size of n = 4 is

(The correct Dixon's test formula for a particular group size must be used).  The data is then ordered and values corresponding to x1, x2 and xn determined.  The outlier of interest is assigned x1 whether it is the smallest or largest value and the rest of the order maintained accordingly.


Oneway DIXONS TEST2
TESTING FOR OUTLIERS	page 2of 2

Dixon's Test

Hence for the E. coli data

	x	x + 1	log10(x + 1)
	0	1	0.00	x1
	13	14	1.15	x2
	42	43	1.63	x3
	61	62	1.79	x4
  Consulting a table of critical values at n = 4 shows that r10 is greater than the table values at all three probabilities quoted, so the outlier is not significantly (P > 0.1) lower than the other values in the group.  Thus the value should not be considered for exclusion from the analysis.  If an outlier is significant it should not automatically be excluded from the data set.  Firstly, field and/or laboratory notebooks should be consulted to see whether the data point is correct or whether there is a good scientific reason for exclusion.


OnewayCOMPARISONS OF MEANS
6. COMPARISONS OF MEANS

  The result of the basic one-way ANOVA only says whether or not there are any significant differences between the means of the groups in the data set.  It does not tell which means are different from each other and a highly significant ANOVA might contain only a few large differences.  Thus we need tests to compare means after the basic analysis is complete.  Note that all these tests must be carried out on the data transformed in the most appropriate way to ensure that the assumptions are upheld.  The tests are divided into two types

 A priori tests or planned comparisons

 A posteriori tests or unplanned comparisons.


Oneway A PRIORI TESTS1
A PRIORI TESTS	page 1 of 2

  The a priori tests sometimes called contrasts also include tests for similarity between means as opposed to differences between means.  For pairwise comparisons a version of the tests that work for both equal and unequal sample sizes can be used.

  These tests are very patchily represented in computer packages.  Minitab does not have any commands for these tests in release 7.1.  SPSS/PC+, Systat and Statgraphics have several tests of this type.  All three of these packages currently allow computation of the following pairwise comparisons and SPSS/PC+ and Systat allow a priori contrasts to be performed as well.


Oneway A PRIORI TESTS2
A PRIORI TESTS 	page 2 of 2

  These tests can only be done for groups which could sensibly have been tested before the experiment was done.  Thus for the E. coli data set the researcher could sensibly have planned to compare the sites immediately above and below the tributary site 3 and 4; (see site details) and so these are candidates for an a priori pairwise comparison.  Similarly the rocky site above tributary station (site 1) could be compared a priori with the other two less rocky sites above the tributary (sites 2 and 3).  These tests are more powerful at detecting differences than the a posteriori tests and so stringent rules must be followed.  These rules are

tests can be applied regardless of significance of basic ANOVA

the maximum number of comparisons allowed is a-1 (where a = total number of groups in the ANOVA).

Because these tests are so powerful they should be used in preference to a posteriori tests wherever validly possible.

Least significant difference (LSD)

Orthogonal decomposition


Oneway LEAST SIGNIFICANT DIFFERENCE
LEAST SIGNIFICANT DIFFERENCE (LSD)

  This pairwise comparison is valid only as an a priori test, despite being linked with a posteriori tests in the SPSS/PC+ command structure.  The LSD is the smallest difference between two means that is significant, one value is calculated for all the comparisons to be made.  Calculation of LSD is straightforward.

where the two-tailed t-statistic is quoted at a probability P and for within groups degrees of freedom (dferror) from the one-way ANOVA table; ni and nj are the number of replicates in the largest and smallest groups.  So for the Ecoli data at P = 0.05

where 0.176 is the MSerror from the log10(x+1) transformed data.  As the difference between the group means at the sites above and below the tributary is 1.096 there is a significant difference between them.


Oneway ORTHOGONAL DECOMPOSITION1
ORTHOGONAL DECOMPOSITON	page 1 of 4

  These tests are to allow comparisons among groups of group means.  The method essentially performs a new one-way ANOVA on a restructured set of data.  A test sum of squares (SStest) is calculated according to the following type of formula

where g1 to g4 are the sums of the values in the four groups used in this test and n1 to n4 are all the numbers of replicates in the corresponding groups.  This formula can be extended (indicated by dots) or contracted to test more or less group means in each of the terms, or by adding more terms.


Oneway ORTHOGONAL DECOMPOSITION2
ORTHOGONAL DECOMPOSITION	page 2 of 4

  The next step is to calculate a test mean square (MStest) and F (Ftest)

where dftest is one less than the number of terms used to calculate SStest.  Ftest is then compared with a table value (FP(dftest, dferror)) and if Ftest > Ftable the result is significant.

  Thus to compare the E. coli data at site 1 with those at sites 2 and 3 (groupsums,g, ) the calculations are

From F-tables,  FP(1,18) is	4.41    for  P = 0.05

8.28    for  P = 0.01

15.4    for  P = 0.001

Hence, there is a significant difference in E. coli numbers due to rockiness above the tributary at the
 P=0.01 but not the P=0.001 level.


Oneway ORTHOGONAL DECOMPOSITION3
ORTHOGONAL DECOMPOSITION 	page 3 of 4

  Particularly with larger data sets these calculations can be tedious and can easily be done by computer by re-labelling the groups appropriately.

	log10	sites	rocky	below
Using the log10(x+1) transformed E. coli data, the labels in the column headed 'rocky' compare means at sites 2and3.   Those in the column headed 'below' compare the means at sites 4, 5 and 6.   Misssing values are labelled with an asterisk	0.903 0.699 0.000 0.778 1.792 1.146 0.000 1.633 1.707 2.193 2.029 2.003 1.079 0.845 0.477 1.146 0.000 0.000 0.000 0.301 0.301 0.699 0.301 0.000	1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6	1 1 1 1 2 2 2 2 2 2 2 2 * * * * * * * * * * * *	* * * * * * * * * * * * 1 1 1 1 2 2 2 2 3 3 3 3
  Click for details of Minitab commands


Oneway ORTHO MTB COMMANDS
ORTHOGONAL DECOMPOSITION	page 4 of 4

Minitab Commands

The command is the same

ONEWAY  C1  C2  C3  C4

C1  contains the transformed data

C2  contains the labels(in this case the 'rocky' or 'below'column labels)

C3  and  C4 (empty)  are where Minitab stores the residuals and group means.

  Performing a one-way ANOVA using either of these new sets of subscripts will show the SStest in the normal SSgroups position.  The Ftest however, must be calculated with the MSerror from the entire original one-way ANOVA table.

  By this method the SStest below the tributaries is 1.382 and so MStest=1.382/2=0.691 and hence
Ftest=0.691/.176=3.926.  In this test dftest is 2 because there were 3 terms compared and can be found as dfgroups in th modified ANOVA table used to calculate SStest .  As F0.05(2,18) is 3.55 there is just a significant difference between thes means

  Four means comparisons have now been done by a priori tests (above and below tributary=1 ; 'rocky'=1; 'below'=2; total=4) and this is very close to the upper limit of 5 (a-1), so it would be unwise to search artificially for a fifth a priori comparison.


 Site details, and group means and sums for the log10 (x+1) transformed E. coli data (n=4)

Site	Group Mean	Group  Sum	Site details
1 2 3 4 5 6	0.595 1.143 1.983 0.887 0.075 0.325	2.380 4.517 7.993 3.547 0.310 1.301	<<<most rocky site tributary enters between 3 and 4


Oneway A POSTERIORI
A POSTERIORI TESTS

  These tests are less stringent than the a priori tests but allow us to look for means and groups of means which appear to be different in the data after it has collected.  This feature makes them very powerful and they can be used extensively.  The basic rules for most a posteriori tests are as follows

the initial ANOVA on all the data must be significant before these tests are attempted

large numbers of comparisons are allowed, always more than (a (a - 1))/2.

  The limit set above is a lower limit and essentially allows all possible pairs to be compared.  As we are rarely interested in absolutely all pairwise comparisons, there are essentially comparisons available for some groups of means to be compared also.

The common methods used are

Tukey-Kramer

Sum of squares, simultaneous test proceedure


Oneway TUKEY-KRAMER
TUKEY-KRAMER METHOD

  The best a posteriori pairwise comparison is the Kramer modification of the Tukey method.  With this test a minimum significant difference (MSD) is calculated which is analogous to the LSD described in the a priori tests.

where Q is the studentized range statistic and so for the E. coli data Q0.05(6,18) is 4.495.

Hence

and so none of the group means below the tributary are different from each other.  We can also use this statistic to investigate the E. coli numbers above the tributary in more detail.  It is clear that site 1 is only individually significantly different from site 3 and not from site 2.  This does not invalidate our previous
a priori conclusion because it is perfectly possible for a priori and a posteriori tests to give conflicting results, because they have markedly different powers.


Oneway SS-STP
SUM OF SQUARES, SIMULTANEOUS TEST PROCEDURE(SS-STP)

  This is a very powerful a posteriori test for comparing groups of means.  The SStest is calculated and then compared with a critical value (SScrit) calculated as follows

SScrit = (a - 1).MSerror.FP(a-1,dferror)

If SStest > SScrit then the comparison being made is significant.  Thus for the E. coli data at P = 0.05.

SScrit = 5 x 0.176 x 2.77 = 2.438

this statistic can be used to further investigate groups of means in the E.col data.  A SStest is calculated longhand or by computer in an identical manner to that previously described for orthogonal decomposition .


Oneway MTB TABLE
SUM OF SQUARES SIMULTANEOUS TEST PROCEDURE

Minitab Commands

The table shows the method of labelling groups for calculating sums of squares by computer for the a posteriori SS-STP test used to compare multiple groups of means.  The log10 E. coli data are used.

	log10	sites	lowest
Missing values are labelled with an asterisk.   The SStest calculated from the column headed 'lowest' would compare the group with the lowest group mean (site 5) below the tributary with the other two below tributary sites (4 and 6).	0.903 0.699 0.000 0.778 1.792 1.146 0.000 1.633 1.707 2.193 2.029 2.003 1.079 0.845 0.477 1.146 0.000 0.000 0.000 0.301 0.301 0.699 0.301 0.000	1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6	* * * * * * * * * * * * 1 1 1 1 2 2 2 2 1 1 1 1
Again the  ONEWAY  command this time using the 'lowest' labels is used.

  In this case SStest was 0.751 and so the difference was not significant.  However, using the same principles a test between sites 1-4 and sites 5 and 6 was significant (SStest = 4.831).


?????

6.3	Computer output for pairwise comparisons

Several statistical computer packages give tables as part of the output for one-way ANOVA which indicate results for pairwise means comparisons.  This output can take several forms.  Systat gives two matrices to indicate differences, one contains absolute mean differences and the other probabilities for the differences being significant, both are organised in an a x a  triangular matrix.  SPSS/PC+ gives either a similar type of matrix but with significant differences, marked with an asterisk and/or a table of homogeneous subsets, which is a list of groups which do not have mean values that are significantly different. Statgraphics indicates homogeneous subjects by labelling the groups in a list with a system of asterisks.


Oneway WHAT TO DO IF
7. WHAT TO DO IF THE ASSUMPTIONS DO NOT HOLD

  In almost all cases the above procedures will find a transformation, if one is required, that ensures that the assumptions of ANOVA are upheld.  However, occasionally data will be collected that cannot be made to fit the assumptions by transformation.

  A summary of what to do under these circumstances is provided in the following Trouble Protocol  which has already been referred to earlier (Protocol 1).


Oneway TROUBLE PROTOCOL
TROUBLE PROTOCOL   What to do if ANOVA assumptions do not hold

1.	If ANOVA shows no significant differences with the best transformation accept result.

2.	If ANOVA shows large differences (P<0.001) and assumptions are only slightly violated, continue with ANOVA but use a more stringent probability level to detect differences
(e.g.P<0.01).  Check results with non-parametric tests (see step3 and 4).

3.	Use a non-parametric test, that is equivalent to one-way ANOVA, to detect a significant overall ANOVA
 Kruskal-Wallis
 Mood median

4.	Use non-parametric tests to compare pairs of means
 notched boxplots of grouped data
 multiple Mann-Whitney tests

5.	Consider using parametric tests to compare pairs of means when variances are unequal
 Games and Howell procedure
 T3 method

Robustness of one-way ANOVA

Non-parametric test equivalent to one-way ANOVA

Comparisons of means

Non-parametric tests

Parametric tests


Oneway ROBUSTNESS OF ANOVA
ROBUSTNESS OF ONE-WAY ANOVA

  One-way ANOVA is a fairly robust procedure and so small deviations from homogeneity of variances, which is the most important assumption, rarely give large errors.  This is particularly true if the basic ANOVA is highly significant and the group sizes are equal.  However, it is wise to confirm the results obtained with the best transformation which can be obtained with the results obtained by non-parametric procedures.


Oneway EQUIVALENT NON-PARAMETRIC TESTS
NON-PARAMETRIC TESTS EQUIVALENT TO ONE-WAY ANOVA

  There are two tests which are widely available, and both determine whether the group medians are significantly different.  Both tests work best with fairly large numbers of observations in each group (n > 6), but the major statistical packages normally give results for smaller n and print warnings.  Both tests perform calculations on the ranks of the numbers in the data and so depend on the order of the data not the value.  For this reason transformation is normally not needed.  Data is usually submitted as a data column and a labelling column as described for one-way ANOVA.

  The Kruskal-Wallis test is the most commonly used procedure of this type.  The result for the E. coli data is P = 0.009 and so is similar to the parametric one-way ANOVA.  The Mood median test is at present only available in Minitab and gives a more conservative result for the E. coli data (P = 0.016).  However, this test also allows calculation of confidence intervals which can be used to make an initial assessment of the differences between means.  These intervals are value sensitive and so should be done on the best data transformation.

  Click to view some example output from these two tests


Oneway EQUIV NON PARA eg1
NON-PARAMETRIC TESTS EQUIVALENT TO ONE-WAY ANOVA

Example Output	page 1 of 2

Kruskal-Wallis analysis of raw by sites (produced by statgraphics)

Test statistic = 15.7839  Significance level = 7.48902E-3


Onewy EQUIV NON-PARA eg2
NON-PARAMETRIC TESTS EQUIVALENT TO ONE-WAY ANOVA

Example Output	page 2 of 2

Mood median test of log10 x +1

The Minitab command is

MOOD  C4  C2

C4       << contains the data (in this case log transformed)

C2       << contains the labels

The output appears


Oneway COMPARISONS OF MEANS_NON-PARA TESTS
COMPARISONS OF MEANS

  Non-parametric tests

  One of the most useful non-parametric test of this type is the notched boxplot, which is yet another exploratory data analysis technique.

The Minitab command and subcommands are

BOXPLOT  C3;   	<< data column

BY  C4;                 	<< labels column

NOTCH [K% confidence]sign ci.         << for confidence intervals

  The 'notches' are confidence intervals provided as extra information to a standard boxplot of all the groups placed side by side.  If the notches of two groups do not overlap then the group medians are significantly different.  This test is value sensitive and so the best transformation should be used.  Notched boxplots are produced by Minitab, Statgraphics and Systat.  The former two packages print the notches in the boxplots as brackets whilst Systat uses spikes to indicate the notches when they lie outside the hinges.  Comparison of the results from the notched boxplot with the matrix of a posteriori comparisons(Systat) for the log10 E. coli data, show that the results are very similar.


Oneway COMPARISONS OF MEANS_PARA TESTS
COMPARISONS OF MEANS

  Parametric tests

  Some authorities have recommended that specially designed parametric tests should be used when variances are not homogeneous.  Two of the best tests of this type are the Games and Howell procedure and the T3-test.  However, computation is in both cases difficult as a separate estimate of the variance and the degrees of freedom must be calculated for each comparison made.  For this reason computers are almost essential for these tests.  Unfortunately they are not provided as standard features in many of the major statistical packages.  There is a suite of programs (BIOM) written in Fortran to accompany the book by Sokal and Rohlf ' which is available free from W.H. Freeman.  This contains a program which performs the Games and Howell procedure (MCHETV).  Although using this program is probably the best approach for those wishing to undertake the best possible test on data with heterogeneous variances the notched boxplot will probably be more cost effective in time.

( ' Sokal,R.R. and Rolf, F.J. (1981). Biometry. W.H. Freeman, San Francisco.)


Bivariate Intro
BIVARIATE REGRESSION

  All biologists will use regression analysis between two variables at one time or another and all statistical packages have commands for performing bivariate or simple linear regression.

  Any data set which has two variables in which the two sets of numbers are paired is a candidate for regression.  This is the type of data which normally plotted on a simple two dimensional graph to examine the relationship between the two variables.

  One variable y would normally be plotted on the vertical axis and is called the dependent variable.  The other variable x is plotted horizontally and is the independent or predictor variable.  Each point on the graph is represented by a yx data pair and there are n points.


Bivariate Reression Example
BIVARIATE REGRESSION EXAMPLE

  Regression calculates a line of best fit through the data and aims to minimize in some way the distance between the data points and the line.

  These distances are called residuals (e) and are calculated from the difference between the y value of the data point, yi, and the value of y at the corresponding value of xi that lies on the line of best fit, called  or the predicted y value.

plot {fig 1}?

  The basic equation for the linear regression of y upon x is

y  =  bo + bx + e

where bo is the [image] value when x = 0, or the intercept of the line of best fit on the y axis, and b is the slope of the line.  bo and b are constants which are normally estimated as the regression coefficients bo and b respectively, our best estimate of the error term [image] comes from the residuals (e1 .... en).


Bivariate Regression Uses
USES OF REGRESSION

  There are two main uses of bivariate regression

 to provide a prediction equation for y based on x

 to understand the relationship between y and x

  Prediction equations are useful for many diverse purposes, for example for calibration curves or for predicting simple relationships like the rate of an enzymic reaction given a certain concentration of substrate.  Such equations use all the data and so provide summary statistics.

  Regression analysis is straightforward with most statistical packages and will provide much better estimates than those obtained from lines fitted by eye.  Regression will also help biologists understand whether variables are linearly related or not, or whether the relationship is linear only after transformation of x or y by, say, log10, square root or reciprocal.  There are also techniques for fitting data which show nonlinear or curved relationships.  Click 'Example' button for some example data sets.


Bivariate Regression Example Data Sets
EXAMPLE DATA SETS	page 1 of 2

           Optical Density	    Child

Dry Wt (ug/ml)	OD (660nm)		Age (months)	Height (cm)
1775                          3.92            1280                          3.26            1110                          2.20              900                          2.00              800                          2.40              388                          0.84            1088                          2.84              480                          1.15            2050                          4.13            1110                          2.46            1210                          2.90            1775                          3.80			109 113 115 116 119 120 121 124 126 129 130 133 134 135 137 139 141 142	137.6 147.8 136.8 140.7 132.7 145.4 135.0 133.0 148.5 148.3 147.5 148.8 133.2 148.7 152.0 150.6 165.3 149.9
  The optical density (OD) data is an example of a calibration curve in which the dry weight of bacteria in a culture could be estimated from the true optical density of the suspension.  This was measured after dilution to less than OD = 0.6 and then back-calculated to allow for the dilution used;  this should theoretically ensure a linear relationship between the two variables.  The child data consists of children's height over a narrow range of ages( [image]).


EXAMPLE DAT SETS 	page 2 of 2

Enzyme

s	vtreat	v-untr
20 20 60 60 110 110 220 220 560 560 1100 1100	76 47 97 107 123 139 159 152 191 201 207 200	67 51 84 86 98 115 131 124 144 158 160	s - substrate concentration (mg/ml)    vtreat - reaction velocity when treated with puromycin (counts [image])    vuntr - untreated reaction velocity               (counts [image])

  The enzyme data are from an experiment to investigate the effect of puromycin treatment on the relationship between substrate concentration and reaction velocity for a simple enzymic reaction ( [image]).  The OD and Child data have one value of y for each x value and the enzyme data have two values of y for all but one x value.  This latter type of data structure is common because biologists often perform replicate estimations.


ROUTE THROUGH A BIVARIATE REGRESSION

  1.	Decide which variable is to be y and which is to be x (N.B. y is predicted from x).

  2.	Plot the data.

  3.	Check assumptions of the regression.

  4.	Transform x and/or y if not even.

  5.	Compute the regression.

  6.	Analyze residuals.

  7.	Compare influence satistics with critical values.

  8.	If outliers or leverage points seem a problem and all the data is accurate consider using a robust method.

  9.	Regression through the origin.

10.	Regression with more than one y for each x (Replicate Protocol).

11.	Curved regression lines (Curve Protocol).

For more details click


Bivar Protocol
BIVARIATE REGRESSION PROTOCOL

1.	Decide which variable is to be y and which is to be x (N.B. y is predicted from x).

2.	Plot data; y as vertical axis, x as horizontal axis.  If line looks curved rather than approximately linear see Curve Protocol .

3.	Check evenness of x and y variables with
 boxplots
 letter value displays
 histograms

4.	Transform x and/or y if not even;  use ladder of powers to find suitable transformation
(e.g. , [image], log10y, 1/x) or the transformation plot for symmetry.

5.	Compute regression, save residuals, fitted y values ( [image]) and influence statistics (leverage, Cook's D and Dfits), calculate Durbin-Watson statistic if data is in a logical order.

6.	Plot studentized deleted or standardized residuals against
 fitted values,
 x variable
 the sequence of observations, if in a logical order.

Examine residual plots for outliers, data points a long way from the rest.  Consider rejection of outliers with studentized deleted or standardized residuals >2 and return to step 5.

7.	Compare influence statistics with critical values
 leverage	>  2p'/n
 Cook's D	>  F0.5(p',n-p') or 4/n
 Dfits	>  absolute value of

Where p' = the number of parameters in the model (bo = 1 parameter; b = 1 parameter);
n = number of data points in the regression.  If two or more influence statistics are greater than the critical values consider rejecting points and return to step 5.

8.	If outliers or leverage points seem a problem and all the data is accurate consider using a robust regression method.
 three-group resistant line
 linear regression using least absolute residuals or minimized [image].


Bivar Plotting the Data
2.PLOTTING THE DATA

  Once it has been decided which is the x variable and which the y variable (step1), it is essential to plot the data with the y variable on the vertical axis and the x variable on the horizontal axis.  In this way it can easily been seen if the variables are linearly related and the degree of variability about the straight line.  This step should never be omitted.

  These plots can easily be performed in Minitab


Bivar Plotting data Mtb pg1
PLOTTING THE DATA	page 1 of 2

Minitab commands and output

For the OD data

 PLOT  C1  C2        [C1  contains the y variable,  C2  contains the x variable]

 Dry Wt  -

         -                                                      *

         -

     1800+                                                 * *

         -

         -

         -

         -                                          *

     1200+                                    *

         -                          *   *     *

         -                        *

         -                             *

         -

      600+

         -           *

         -       *

         -

           ----+---------+---------+---------+---------+---------+--OD

            0.70      1.40      2.10      2.80      3.50      4.20

  The data is clearly linearly related with little scatter.


Bivar Plotting data Mtb pg2
PLOTTING THE DATA	page 2 of 2

Minitab commands and output

For the enzyme data a multiple plot is used

MPLOT C2 C1 C3 C1           [C1 contains the x variable, C2 & C3 contain the y variables]

         -

         -                             A                          2

         -                             A

      180+

         -

         -            2                B                          B

         -       A                     B

         -            B

      120+       2    B

         -    A

         -    A  B

         -    2

         -  2

       60+

         -  2

         -

         -

           +---------+---------+---------+---------+---------+------

           0       200       400       600       800      1000

  The data is nonlinear and shows two clear curves as velocity (v) reaches a maximum value with increasing substrate concentration (s);  this appears to be a typical Michaelis-Menten relationship.

  Both the plots were made with simple text based output and illustrate that this style of plot is quite sufficient for most purposes during routine regression analysis.  It is only necessary to use higher-resolution graphics output for special purposes (e.g. for confidence intervals and lines of best fit) or for final presentation of the data.


Bivar Assumptions
3.ASSUMPTIONS OF BIVARIATE REGRESSION

  The basic assumptions upon which the calculations done by the standard least squares regression analysis available in most statistical packages are

 x and y are linearly related by equation (1)

 x and y are measured for each of n observations

 x is measured without error

 y is a set of random observations measured with error

 the variance of the y observations is constant

 the errors ( [image]) of y are normally distributed

 the errors of y are mutually independant

  The assumption of normality of errors is only required when significance tests are done on the regression statistics, confidence or prediction intervals are calculated or when regression lines are compared, however a check is needed for non-constant variance.  The errors do not need to be normally distributed to calculate a line of best fit nor to calculate predicted values from this line.

Missing x or y values

Checking the evenness of x and y


Missing x or y values
  Most computer packages will produce output from data with missing x or y values.  It can sometimes be useful to provide extra x values without corresponding y values.  If this is done the data with both x and y values will be used for calculating the basic regression statistics, such as the line of best fit and the residuals.  However, [image], some influence statistics (e.g. leverage ) and the confidence intervals will be calculated for all x values.


Errors of x and y

  In theory only y should be measured with error, if x also has error then the regression calculations will be wrong.  Although it is possible to perform the regression calculations when x has error.  However, provided the errors of y are much greater than the errors of x few problems are likely to be encountered.


Bivar Evenness of x & y pg1
EVENNESS OF X AND Y	page 1 of 2

  Regression analysis is far more likely to calculate an accurate line of best fit if the values of x and y are evenly spread through the range.  This is easily tested with simple graphical procedures (Bivariate Regression Protocol) like boxplots and histograms.  Letter value displays can also be used in which the midsummaries should be equal in a perfectly symmetrical distribution (more detailed information about letter value displays can be found in One-way ANOVA Assumptions)


bivar Evenness of x & y pg2
EVENNESS OF X AND Y	page 2 of 2

Boxplots of some variables from the example data sets. (Using the BOXPLOT Minitab command)

                                  ---------------------

                ------------------I       +           I---------

                                  ---------------------

          ----+---------+---------+---------+---------+---------+--OD

           0.70      1.40      2.10      2.80      3.50      4.20

                            ---------------------

               -------------I       +           I---------------

                            ---------------------

          ----+---------+---------+---------+---------+---------+--Dry Wt

            350       700      1050      1400      1750      2100

                              --------------------------------

            ------------------I              +               I----

                              --------------------------------

          ------+---------+---------+---------+---------+---------+vtreat

               60        90       120       150       180       210

  All are acceptably even.


Bivar Transformations
4.TRANSFORMATIONS

  If a variable was very uneven with several outside and/or far outside values on one side of the boxplot only then transformation might be necessary.  Power transformations from the ladder of powers can be used and the best one found by trial and error, moving along the ladder according to the shape of the distribution obtained.  Data skewed to the right (with a long right tail) will use lower powers and left skewed data (with a long left tail) will need higher powers.

Ladder of Powers

  Emerson and Soto describe a useful method for choosing a power transformation that will make data symmetrical.  This is an exploratory data analysis technique and uses letter value displays.

  It is called the Transformation plot for symmetry


Emerson and Soto ref?


Bivar Ladder of Powers
TRANSFORMATIONS

Ladder of powers

  The ladder puts the raw data and a range of equivalent transformations into a logical sequence from which a suitable transformation can be selected.

	Transformation		Addition of small  number needed for zero values
Power	Name	Formula
3 2 1 0.5 0 -0.5 -1 -2	cube square raw data square root logarythm reciprocal root reciprocal reciprocal square	x log10x 1/ 1/x 1/	no no no   no* yes yes yes yes

Based on a table in Velleman and Hoaglin

* Not needed but [image] or [image] is often used.


Reference

Velleman P.F. and Hoaglin, D.C.  (1981).  Applications, Basics, and Computing of Exploratory Data Analysis.  Duxberry Press, Boston.


Bivar Transformation Plot for Symmetry pg1
TRANSFORMATIONS -- The Transfomation Plot for Symmetry       page 1 of 2

  From the upper ( LU ) and lower ( LL ) letter values and the median (M) of the variable, two sets of statistics are calculated for each letter value.

Letter  Value	LL	LU
H E D	102 76 47	195.5 201 207	3.25 -7.0 -18.5	7.55 13.59 23.17

M=145.5

These are then plotted against each other with

on the vertical axis and    [image]  on the horizontal axis.


Bivar Transformation Plot for Symmetry pg2
TRANSFORMATIONS -- The Transfomation Plot for Symmetry           page 2 of 2

For the puromycin treated velocities for the enzyme data (which is slightly skewed left)

  The slope (b) of the transformation plot is -1.37 and the appropriate power (p) transformation is calculated by  p  =  1 - b  which in this case is 2.37 (p = 1 + 1.37).  This is close to the square transformation ( [image] [image] [image]) which gives an almost perfectly symmetrical boxplot.  This plot works better and is more useful with larger data sets.


Bivar Computing The Regression
5.COMPUTING THE REGRESSION

  Almost all standard regression commands in computer packages perform least squares regression as the main primary procedure.  In this type of regression the residuals are optimized by minimising the sum of the squares of the residuals, and so the line of best fit calculated by this method is very sensitive to extreme data points.

  Although other methods are more robust, least squares regression is a very powerful approach, from which a lot of useful information can be calculated, which cannot be estimated using the more robust approaches.  So it is probably the best method to try first.


Bivar Computing the Regression Mtb
COMPUTING THE REGRESSION

Minitab Commands and Output

  With Minitab the variables are defined by name or column number (assume 'Age' in column 2 and 'Height' in column 1) and the number of x variables must also be specified;  hence the command is one of the following:

     REGRESS  'Height'  1  'Age'
REGRESS    C2   1   C1

  As can be seen in the Bivariate Regression Protocol it is useful to save the residuals, y values and influence statistics (leverage, Cook's D and Dfits).  This is normally done in the same command set as the basic regression definition, with subcommands and/or command line extensions.

The standard output from Minitab is

The regression equation is

OD = 79.7 + 0.511 Dry Wt

Predictor Constant Dry Wt	Coef 79.701 0.5113	Stdev 21.25 0.1670	t-ratio 3.75 3.06	p 0.002 0.007
s= 7.029	R-sq = 36.9%	R-sq(adj) = 33.0%

Minitab also outputs an analysis of variance table and any unusual observations.

Calculating Confidence Intervals

Systat and SPSS/PC+ Commands


This is the value of the intercept (bo) of the regression equation which calculates new values of Height ( [image]) from Age.


This is the value of the slope (b)


This column gives the values of the intercept and slope


This column gives the standard errors of the two coefficients (called Stdev here; sb for b and sbo for bo).  From these the 95% confidence intervals of the coefficients can be calculated.


The P values indicate whether the coefficients are significantly different from zero;  P<0.05 indicating significance.  Thus for these data both the intercept and slope are significant.


s  is the estimated standard deviation of the errors in y and represents the deviation about the regression line.


R-sq is more commonly called [image] which is the coefficient of determination and indicates the percentage of the variation in the data explained by the straight line relationship.  The square root of [image] is the linear correlation coefficient (r), which in this case when expressed as a proportion is 0.67 [image]. This figure can be compared with table values as the absolute value of the coefficient.  Here rtable = 0.468 (n=18) which is greater than the calculated value and so there is significant linear correlation between Age and Height.


The adjusted [image] can be ignored.


The analysis of variance table for the regression model used (bivariate linear regression) appears below the regression output.  This is similar to the one-way ANOVA table.  The first line gives the F and corresponding P value for the overall fit of the regression model, again significant.  The error mean square (49.4) is the variance about the regression line and will be designated MSrerror here.  Note [image] [image] is the same as s defined earlier in the output.


The unusual observations are given as warnings of points which are possibly extreme outliers or have a strong influence.


Bivar Calculating Confidence Intervals
COMPUTING THE REGRESSION

Calculating Confidence Intervals

  The 95% confidence intervals can be calculated from the standard errors (Stdev column) of the two coefficients.  This is done by multiplying by t0.05(n-2) from tables.  Note that the degrees of freedom for these data is 16.  More generally n-2 is derived from n-p', where p' is the number of regression coefficients.  In this case p' = 2, the coefficients being the intercept and the slope.  So in this case the coefficients and their 95% confidence intervals are

b ± sbt0.05(10)  =  0.0019715 ± (0.167 x 2.12) = 0.5113 ± 0.354

bo ± sbot0.05(16)  =  79.70 ± (21.25 x 2.12) = 79.70 ± 45.05

the large range in these intervals reflects the wide variation observed in the original plot.


Bivar Systat & SPSS/PC Commands
COMPUTING THE REGRESSION

Systat and SPSS/PC Commands

  With Systat the regression model is specified in a similar form to the regression equation, hence the command is

MODEL OD = CONSTANT + Age

which is part of the MGLH commands.  Alternatively equivalent boxes can be filled in the menu mode of the program.  Statgraphics relies entirely on filling in boxes for dependent variable (y), independent variable (x) and toggling the model box to linear.  With SPSS/PC+ the corresponding command line would be

REGRESSION VARIABLES = AGE, HEIGHT \DEPENDENT = HEIGHT.

  Again it is useful to save the residuals, y values and influence statistics (leverage, Cook's D and Dfits).  This is normally done in the same command set as the basic regression definition, with subcommands and/or command line extensions.


Bivar Residual Analysis
6.ANALYZING RESIDUALS

Residual plots and outliers

  Examination of residuals is essential for checking that the linear model assumed by the regression calculations and other assumptions are valid for the data set being analysed.  Three types of residuals are commonly provided by computer packages.

 ordinary residuals

 standardized residuals

 studentized residuals

All will give the same pattern of points in the residual plots.  The ordinary residuals are expressed in the original y units whilst the other two are in units of standard deviation.  It is normally best to examine the studentized residuals, although the standardized residuals are a good substitute but are more conservative for outlier detection.  Studentized residuals are sometimes called studentized deleted or just deleted residuals.  The ordinary residuals should be examined to see the magnitude of the errors in the original units; this is particularly useful for calibration curves (e.g. the OD data).


Bivar Residual Plots and Outliers pg1
RESIDUAL PLOTS AND OUTLIERS	page 1 of 2

  The studentized or standardized residuals should be plotted against the

 fitted values,

 x variable

 the sequence of the observations if they are in a logical order

  These plots are then examined along with the absolute magnitude of the ordinary residuals.  A well fitted model will give a random distribution of residuals around zero.

fig6??

  Outliers with studentized or standardized residuals >+2 or <-2 (OD data) should be considered for rejection (see Multiple Regression for a more accurate test).


Bivar Residual Plots and Outliers pg2
RESIDUAL PLOTS AND OUTLIERS	page 2 of 2

  When considering points for rejection it is important to

 examine the residual plots to ensure they are not just the extreme values of an unusually variable data set

 examine laboratory records and field note books to see if there are any basic problems with the recording of the data, the experimentation or whether a transcription error has occurred.

  Exclusion of data points from the analysis that are not real errors is very dangerous.

  Extreme outliers may tell you something important about the data.  For example if the outlier is the highest value of x or y it might indicate a region of nonlinearity just outside the range of observations (the OD data has an example of this).

  Fitting a linear model to a curve is evident from a sweeping arc of residuals (Enzyme data).  This is evident even if the magnitude of the ordinary residuals is very small in relation to the y values ( [image]).

Other patterns:    Sloping line of residuals

  Funnel shaped residuals

Example residual plot

Durbin-Watson test for autocorrelation of residuals


Bivar Sloping line of Residuals
RESIDUAL PLOTS AND OUTLIERS

Sloping Line of Residuals

  If one or a few data points have excessive influence on the slope of the regression line, this is easily spotted as a sloping line of residuals.

fig6b?


Bivar Funnel Shaped Residuals
RESIDUAL PLOTS AND OULIERS

Funnel Shaped Residuals

  Violation of the assumption of constant variance is seen by a funnel shaped pattern.  This is best corrected by transformation of y by square root, log10 or reciprocal.

fig6d?


Bivar Example Residual Plot
RESIDUAL PLOTS AND OUTLIERS

Example Residual Plot

  The two residual plots of the data sets gave essentially the same pattern.

Standardized residuals

  This shows that two points are possible outliers (ringed), they were also flagged as unusual observations.  There is a possibility of a poorly fitted line, indicated by the slope on the remaining residuals.  There is also the suspicion of funnelling this however, is probably due to the outliers and so can be discounted.


Bivar Durbin-Watson Test
THE DURBIN-WATSON TEST

  This test is commonly available in computer packages and tests for the significance of autocorrelation between residuals.  Hence it will pick out as significant the problems identified in the residual plots (outliers or curves).  It should only be used if the data is in a logical order.Thus it can be done on the Child data which is ordered by age, but the OD data would need sorting by OD before the result would be meaningful.

  The Minitab subcommand is

REGRESS  C2  1  C1  C3  C4;     (C3 and C4 empty for residuals and fits)

DW.

  The Durbin-Watson statistic (DW) for the Child data is 2.21.  This is compared with table values.  For the Child data the table value is (n = 18, k = 1) are dL = 1.158 and dU = 1.391, hence there is no doubt as 2.21> dU and so there is no significant autocorrelation in the residuals.


Bivar Influence Statistics
7. INFLUENCE STATISTICS

  There are statistics which test whether data points have a large influence on the slope of the regression line.  The most common statistics of this type calculated by statistical packages are leverage (also called HI in Minitab, and LEVER in SPSS/PC+), Cook's D and Dfits.

  Leverage (hii) is the diagonal of the hat matrix and is very useful for the calculation of
Confidence and Prediction Intervals.


Bivar Influence Stats Mtb pg1
INFLUENCE STATISTICS 	page 1 of 2

Minitab Commands and Output

  All three statistics should be calculated and compared with critical values.

The Minitab commands are

REGRESS  C2 1 C1 C3 C4;

COOKD  C5;	}

HI  C6;	} >>  empty columns

DFITS  C7.	}

  It is common for several different data points to be identified as having high influence by these statistics.  Thus at least two of the statistics should be above the critical value before they are considered as truly influential.  Once identified such influence points should be considered carefully for possible exclusion.  The same cautionary notes as given earlier for elimination of points with high residuals apply equally here.


Bivar Influence Stats Mtb pg2
INFLUENCE STATISTICS	page 2 of 2

Minitab Commands and Output

  The values for the three influence statistics from the child data (influencial points above the critical values are underlined).

Child number	leverage	Cook's D	D fits	Critical values and Formulae
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18	0.235 0.164 0.135 0.122 0.090 0.082 0.075 0.060 0.056 0.058 0.061 0.077 0.085 0.093 0.114 0.139 0.169 0.185	0.019 0.253 0.005 0.005 0.068 0.019 0.038 0.070 0.012 0.005 0.001 0.001 0.230 0.000 0.007 0.000 0.452 0.016	0.191  0.751 -0.099  0.093 -0.373  0.189 -0.275 -0.390  0.153  0.094  0.049  0.046 -0.791 -0.001  0.119 -0.010  1.083 -0.176	leverage   = 0.222  ( >2p'/n)      Cook's D = (i) 0.222  (>4/n)                    (ii) 0.724  (>F0.5(p',n-p') Dfits         = 0.667  (2 for Dfits compare the absolute value (ignore the sign) p' = the number of parameters in the model(bo = 1; b = 1) n = number of data points in the    regression

  No data point is identified by all the statistics, but the points for child 13 and 17 are identified for the lower critical value of Cook's D and for Dfits as well.  These points are clearly suspicious because they were identified as outliers from the residuals earlier.  Additionally, the point for child 2 is indicated by two of the statistics and so is a possible influence point.  The point for child 1 is only high for leverage and was not an outlier so it is not considered further.

  Note also that for these data the two critical values which have been suggested for Cook's D give different results.  The value for 4/n will always be lower than that for F0.5(p',n-p').  So the latter statistic will only indicate points that are very strongly influential.

  This is quite a large number of influence points and outliers for such a small data set.  So we should be suspicious and in the real world probably collect more data.  However, without that option the line of best fit is probably best described by a more robust type of regression (Bivariate Regression Protocol).


Bivar Confidence & Prediction Intervals pg1
CONFIDENCE AND PREDICTION INTERVALS	page 1 of 4

  Most packages allow calculation of confidence and prediction intervals and some allow these values to be stored for subsequent use in graphs (e.g. Statgraphics).  If they are not stored they are easily calculated from the leverage values (hii) which are more commonly available as a stored variable.  The formulae for these intervals

Statistic	Symbol for standard error	Formula for standard error	Formula for interval
Confidence intervals Regression coefficient, b Intercept, bo           One estimated [image] from a given x, CIt All [image] values, CIF Predicion intervals One new value of [image] from a new x value, PIt All new values of [image] from any number of new x values, PIF		_b _b


 These formulae can be applied to all x values in the regression and the result stored as a column for plotting.

MSrerror = mean square from the regression ANOVA; hii = leverage values (one hii valufor each x value)


 These standard errors are usually on the standard regession output


Bivar Confidence & Prediction Intervals pg2
CONFIDENCE AND PREDICTION INTERVALS	page 2 of 4

Minitab Commands and Output

  The Minitab PREDICT subcommand is used for all x values for the whole OD data with a constant term.  The output appears

Fit  Stdev.Fit	95% C.I.	95% P.I.
1758.7           65.5 1447.5           47.4  947.7           44.7  853.4           48.6 1042.0           41.9  306.4           84.6 1249.5           41.2  452.6           73.7 1857.8           72.5 1070.3           41.7 1277.8           41.7 1702.2           61.8	( 1612.7,  1904.8 ) ( 1341.9,  1553.2 ) (  848.2,  1047.2 ) (  745.0,   961.8 ) (  948.6,  1135.4 ) (  117.9,   494.9 ) ( 1157.6,  1341.4 ) (  288.3,   616.9 ) ( 1696.2,  2019.3 ) (  978.2,  1162.5 ) ( 1184.8,  1370.8 ) ( 1564.5,  1839.8 )	( 1413.3,  2104.2) ( 1117.1,  1778.0) (  619.2,  1276.2) (  522.1,  1184.7) (  715.3,  1368.7) (  -59.0,   671.9) (  923.2,  1575.8) (   99.0,   806.2) ( 1505.5,  2210.1) (  744.0,  1396.7) (  951.2,  1604.4) ( 1360.2,  2044.2)
  Note that in this example the column called Stdev.Fit contains the [image] values used to calculate the two types of confidence interval for [image] given in the previous table.  Also the 95% confidence and prediction intervals (called CI and PI respectively) are those based on t (CIt, PIt) and so only properly apply to one [image].  It is better to use the F or Scheffé's intervals (CIF, PIF) which apply to all [image] values ( [image]) and so are wider.


Bivar Confidence & Prediction Intervals pg3
CONFIDENCE AND PREDICTION INTERVALS	page 3 of 4

Minitab Commands and Output

  These four intervals have been calculated with a Minitab macro and are plotted with the line of best fit

  It can be seen that they are bi-concave about the mean values of x and y.  This indicates greater certainty for [image] near the middle of the data and less certainty close to the extremities.  When presenting data graphically with confidence and prediction intervals it is better to use the wider F intervals, as a continuous line implies all possible points have been considered.


Bivar Confidence & Prediction Intervals pg4
CONFIDENCE AND PREDICTION INTERVALS 	page 4 of 4

Minitab Commands and Output

  The shape of these intervals when the line is forced through the origin is different

  In this case the confidence intervals converge towards the origin as the error at this point must be zero if it is known with certainty.


Bivar Robust Linear Regression
8. ROBUST LINEAR REGRESSION

  When there are severe problems with large residuals or influence statistics it could be considered better to use a robust regression method which is less susceptible to extreme data points than the least squares method.

Two methods of robust regression are

The three-group resistant line

Least absolute residuals and minimum


Bivar Three Group Resistant Line pg1
THREE GROUP RESISTANT LINE	page 1 of 2

  This technique is perhaps the most resistant method to outliers and influence points.  Basically the method divides the data points into three groups and takes the median of each group and then calculates a straight line to fit through these medians.

  Minitab calculates the slope and intercept (called level) by this procedure with the RLINE command which has a similar form to the REGRESS command.  As with least squares regression, residuals and fitted values can be calculated, stored and examined to check the efficacy of the model as described previously.


Bivar Three Group Resistant Line pg2
THREE GROUP RESISTANT LINE 	page 2 of 2

  The table illustrates the robustness of this form of regression, which compares the slope and intercepts for both least squares and resistant line methods for the child data with outliers(child 13 and 7) and influence points(child 2 and 17) removed.

	Data set used
Method and coefficients	Full	Without influence points	Without outliers	Without outliers and influence points
Least squares slope intercept Resistant line slope intercept	0.5113 79.70 36.9 0.4286 91.01	0.4757 82.72 40.2 0.4296 91.12	0.4563 86.69 45.1 0.4333 90.40	0.5525 73.90 60.1 0.4333 90.37

  It can be seen that the coefficients hardly vary with the resistant line but they vary greatly for the least squares calculations.  Also these latter values never quite reach those for the resistant line, despite large increases in the [image] value, indicating improved fit.


Bivar Least Absolute Residuals & Minimum pg1
LEAST ABSOLUTE RESIDUALS AND MINIMUM 	page 1 of 2

  These two methods use approaches other than least squares to minimize the distance between the data points and the line of best fit.  The values minimized are the absolute value of y - [image] for the least absolute residuals (LAR) method and [image], which is

for the minimum [image] method.  Many statistical packages cannot do this type of calculation and of the four used here only Systat was capable of doing such regressions.  Systat's nonlinear regression procedure

(NONLIN) must be used because with this it is possible to minimize any function with the LOSS statement.  The regression model must also be given which in this case is the normal regression equation y = bo + bx.  It is not possible to find a nonlinear solution for the child data set with the algorithms used in Systat, so the OD data are used in the following examples.


Bivar Least Absolute Residuals & Minimum pg2
LEAST ABSOLUTE RESIDUALS AND MINIMUM 	page 2 of 2

  The Systat commands needed for LAR regression were

MODEL DRYWT = CONSTANT + B1 * OD

LOSS = ABS(DRYWT - ESTIMATE)

ESTIMATE / SIMPLEX

and for the minimum [image] method were

MODEL DRYWT = CONSTANT + B1 * OD

LOSS = ((DRYWT - ESTIMATE)^2) / DRYWT

ESTIMATE / QUASI

  In these command sets bo is called CONSTANT and b1 B1, but [image] must be called ESTIMATE. The variables are called DRYWT and OD as before. Both these methods produce results different from the least squares regression approach, as do the three-group resistant line procedure from Minitab.  However, for these data and regression methods, the differences between the lines of best fit and the residuals are small.  Table of Results

  Hoaglin et al. ( [image]) give detailed arguments about the robustness of several different regression methods. They conclude that the three-group resistant line method is best, followed by the LAR and least squares methods, in that order.  They do not consider the [image] method.  Thus only with data having very extreme outliers or influence points would it be worthwhile considering the LAR or [image] approaches when the three-group resistant line method is available.


Regression coeficients for the full OD data calculated by four different methods.

	Regression method used
Coefficient	least squares	least absolute residuals	minimum	three-group resistant line
Intercept Slope	-87.7 472 93.0	-11.5 456	-23.9 441	-36.1 464


Bivar Regression Through the Origin
9.REGRESSION THROUGH THE ORIGIN

  With some data sets it is known with certainty that the regression line must pass through the origin
(y = 0;  x = 0).  The OD data is an example where a zero OD must mean there were no bacterial cells in the culture, or at least a vanishingly small number on a linear scale.

Forcing the regression line through the origin is straightforward in most computer packages.

  Minitab uses the NOCONSTANT command placed before the REGRESS statement.

  SPSS/PC+ uses the NOORIGIN subcommand.

  Systat uses a modified regression model statement MODEL DRYWT = OD


Bivar Example Regression Through the Origin
REGRESSION THROUGH THE ORIGIN

Example

  For the OD data the regression equation with a constant is

DRYWT = -89.7 + (472 x OD)

   and without a constant (i.e. line forced through origin) is

DRYWT = 442 x OD

without the constant, [image] is not calculated but all the other regression statistics are given.  It is clear from the equations above that the slope in the model without a constant is modified to allow the line of best fit to go through the origin.  The model without a constant also gives point 9
(OD = 4.13;  DRYWT = 2050) as influential with all three statistics.  This indicates possible lack of linearity at the extremes of the data, so more data in this region is required to test this possibility further.

  Care should be taken when forcing regression lines through the origin.  For example, if the true model is curved or changes slope close to the origin, the normal regression equation will probably be more appropriate.  The chance of this problem occurring undetected can be minimized if data points are collected very close to zero.  The residuals should also be examined for any sign of curvature, so departures from linearity can be spotted.  The Durbin-Watson test for autocorrelation must be modified for regression through the origin and different tables for DW are required ( [image]).


Bivar Regression with more than 1 y
10.REGRESSION WITH MORE THAN ONE Y FOR EACH X

  This can be thought of as regression with replicates of the y values.  The replicates provide extra information about the errors in the data.  Thus there are now two sources of information valuable for regression analysis which are the

 standard regression ANOVA table

 oneway ANOVA table calculated with the replicates of y for each x assigned to a different group

  The oneway ANOVA provides the MSerror which is an estimate of the average variance of the replicate y values.  This has been called the pure or experimental error.  The differences between the degrees of freedom and the sums of squares between the error lines in the two ANOVA tables provide information for a lack of fit test for the linear regression model.  This test has also been called a test for deviations from regression and misspecification of the regression model.


Bivar More than 1 y  pg1
REGRESSION WITH MORE THAN ONE Y FOR EACH X	page 1 of 2

  Many suggest comparison of the regression mean square with the pure error mean square (MSerror from the one-way ANOVA).  This is certainly the best test for the significance of the regression when the lack of fit test is significant.

  However, when lack of fit is not significant (i.e. when the linear model is acceptable), it is often better to use the ordinary regression error mean square (MSrerror from the standard regression ANOVA), which usually has more degrees of freedom and so is a more accurate estimate of the variance ( [image]).

  These tests allow a corrected regression ANOVA table to be constructed.  Unfortunately few statistical packages allow this table to be calculated.  Minitab, however, does give the F and P values for the lack of fit test with the PURE subcommand.

  For these reasons the following explains how to calculate the complete, corrected table.


Bivar More than 1 y  pg2
REGRESSION WITH MORE THAN ONE Y FOR EACH X	page 2 of 2

Calculation of the full, corrected ANOVA table.

1.	Calculate a One-way ANOVA table with the replicate y values in each group.  This will normally involve labelling the groups with integers (see One-way ANOVA).

2.	Calculate a standard regression ANOVA table .

3.	Construct a corrected regression ANOVA table using values from the first two tables.

4.	Test for lack of fit.

5.	Test for linear regression.

Click to view full Protocol, (test and table details)

Comparison of Regression Lines


Bivar Replcate Protocol
REPLICATE PROTOCOL

1.	Calculate a one-way ANOVA table with the replicate y values in each group.  This will normally involve labelling the groups with integers (see One-way ANOVA).  Note that the number of replicate y values need not be the same for each x value. Assign letters* to the different parts of the ANOVA table as below.

Source Between groups Within groups Total	df a e h	SS b f i	MS c g	F d
2.	Calculate standard regression ANOVA table
Assign letters* to the different parts of the ANOVA table as below.

Source Regression Error Total	df A E H	SS B F I	MS C G	F D
3.	Construct a corrected regression ANOVA table as below.

Source Between groups of y	df a	SS b	MS c	F d
Linear Regression Lack of fit Pure error Total	A E-e = Z e H	B F-f = Y f I	C Y/Z = X g	C/g or C/G X /g
4.	Test for lack of fit as follows
   if Ftable < Fcalc (X/g)  then lack of fit is significant
   if Ftable > Fcalc (X/g)  then lack of fit is not significant
(i.e. the model is acceptably linear).

5.	Test for linear regression as follows
   use (C/g) when lack of fit is significant
   use (C/g) when lack of fit is not significant and E >> e.

* symbols used in the ANOVA tables in this protocol do not apply elsewhere.


Bivar Replicate Regression Example
REPLICATE Y REGRESSION EXAMPLE

  The puromycin treated enzyme data has replicate y values, although it is not linear.  The Hane's plot of such data obeying Michaelis-Menten kinetics is the best linearization (see Curved Regression)( [image];).  This plots s/v (where s is the substrate concentration and v is the reaction velocity) as the y variable against s as the x variable.  The enzyme data is acceptably linear when plotted in this way ( [image] = 99.8; residuals random).

The One-way ANOVA Table

Standard Regression ANOVA Table

Corrected Regression ANOVA Table

  The between groups of y line shows there are significant differences between the replicate y values.  The linear regression line shows the overall linear model is highly significant.  The lack of fit line shows that none of the points deviate significantly from a straight line.

  The confidence and prediction intervals for regressions with replicate y values are calculated as usual.  Most authors are agreed on this; however, Sokal and Rohlf ( [image]) give different formulae for these intervals with replicate y's, but as they are in a minority the formulae in the confidence and prediction interval table are satisfactory.


Oneway ANOVA table

Source Between groups Within groups Total	df  5  6 11	SS 37.041520 0.05012 37.46532	MS 7.48304 0.00835	F 895.87	P <<0.001


Standard Regression ANOVA table

Source Linear regression Error Total	df  1 10 11	SS 37.388  0.077 37.465	MS 37.388  0.008	F 4863.08	P <<0.001


Corrected Regression ANOVA table

Source Between groups of y Linear regression Lack of fit Pure error Total	df 5 1 4 6 11	SS 37.415 37.388 0.0269 0.05012 37.46532	MS 7.483 37.388 0.00673 0.00835	F 895.87 4477.6 0.805	P <<0.001 <<0.001 >0.5


Bivar Comparison of Regression Lines pg1
COMPARISON OF REGRESSION LINES	page 1 of 3

  This is a common requirement amongst biologists and analysis of covariance is best used for this, which works for any number of lines and is not affected by replicated y values.  Most packages use model specifications for this type of analysis.  Using Minitab, with the Hane's linearization of the enzyme data(s/v against s) described earlier, the data are entered

x 20 20 60 60 110 110 220 220 560 560 1100 1100 20 20 60 60 110 110 220 220 560 560 1100	y 0.26326 0.42553 0.61856 0.45075 0.89431 0.79137 1.38365 1.44737 2.93194 2.78607 5.31401 5.50000 0.29851 0.39216 0.71429 0.69767 1.12245 0.95652 1.67939 1.77419 3.88889 3.54430 6.87500	lines   1                              x = s, substrate concentration   1   1   1   1                              y = s/v, substrate concentration divided   1                                               by velocity   1   1   1   1                              lines = labels for the two regression lines   1                                            (treated = 1, untreated = 2)   2   2   2   2   2   2   2   2   2   2   2


Bivar Comparison of Reg Lines pg2
COMPARISON OF REGRESSION LINES	page 2 of 3

The Minitab command used to carry out the analysis is

GLM Y = lines ¦ X;	(Note Y, lines and X are column names)

COVARIATE  X.

and the output appears

Factor   Levels     Values

lines	        2         1          2

Analysis of Variance for y

Source lines x lines * x Error Total	DF 1 1 1 19 22	Seq SS  0.042 76.543  1.281  0.202 78.067	Adj SS  0.022 77.202  1.281 0 0.202	Adj MS  0.002 77.202  1.281  0.011	F 0.18 7260.43 120.42	P 0.672 0.000 0.000
Term Constant x x * lines         1	Coeff 0.32644 0.005309 -0.000684	Stdev 0.02871 0.000062 0.000062	t-value 11.37 85.21 -10.97	P 0.000 0.000 0.000


Bivar Comparison of Reg Lines pg3
COMPARISON OF REGRESSION LINES	Page 3 of 3

  The three top lines in the ANOVA table are interpreted from their P values as follows

(a)  lines - significance of the differences between the intercepts

(b)  x - whether the slopes, if parallel are significantly different from zero

(c)  lines * x  -  significance of difference between the slopes

  In this case the intercepts are not significantly different but the slopes are significantly different from each other and zero.  If the intercepts and/or the slopes are not significantly different, then their common values are given as the coefficients for constant and x in the table below the ANOVA table.  So for these data the common intercept is 0.326, but the two seperate slopes must be obtained from two separate regression analyses, one for each line.

  If many regression lines need to be compared then an overall test should be done first, followed by multiple tests comparing seperate pairs of lines until it is clear which lines are different.  For such multiple tests the acceptable level of significance needs to be reduced, because with multiple comparisons the chance of an error increases.  A suitable significance level can be calculated by using Bonferroni's approximation (a/m), where a is the normally acceptable probability and m is the number of comparisons made.  Thus for making 5 comparisons from 10 lines the significance level used should beP < 0.05/5 = 0.01.


Bivar Curved Regression Lines
11. CURVED REGRESSION LINES

  If the initial plot of y against x looks curved or ther is a curved pattern to the residual plots and a significant Durbin -Watson statistic; then linear regression is inappropriate and another method is needed.  There are several methods suitable for curvilinear regression and these can be chosen using the following steps

1.	Select the most appropriate approach from the following guidelines:

(a) Polynomial regression (step 2)
(b) Transformations (steps 3 and 4)

(c) Direct nonlinear solutions (go to step 5)

2.	Perform polynomial regression using best set of predictor variables

3.	Transformations for 'one-bend' curves

4.	Transformations for 'two-bend' curves

5.	Direct nonlinear solutions

6.	Calculate regression statistics and check regression diagnostics

7.	Calculate and plot confidence and prediction intervals as required

For full Details

Choosing the best approach


Bivar Curved Protocol
CURVILINEAR REGRESSION PROTOCOL

1.	Select the most appropriate approach from the following guidelines:
(a) Polynomial regression (step 2)
  use as a purely empirical method of fitting a curve
  gives no idea of functional relationships.
(b) Transformations (steps 3 and 4)
  use for known relationships (e.g. log10y against x) or when a knowledge of the relationship is required
(c) Direct nonlinear solutions (go to step 5)
  use when an equation relating y and x is known from theory.

2.	Perform polynomial regression as follows:
(a) Transform x to [image], [image], [image].... as required to make a set of predictor variables
(b) Pick best set of predictor variables to fit a polynomial equation.  (N.B. See Multiple Regression for methods; the Minitab BREG command works well).

3.	Transformations for 'one-bend' curves
(a) Use the power transformations(ladder of powers; one-wayANOVA) for y and /or x
(b) Choose powers based on
  theory
  shape of curve
  comparisons of plots selected from a correlation matrix of [image] against [image]; where p = 3 to -2.

4.	Transformations for 'two-bend' curves
Choose appropriate transformation from the following:
  logit y = loge [y/(1-y)]
  generalised 'two bend' y = loge [y/(K-y)]
  probit y

5.	Direct nonlinear solutions
Use correct equation from biological theory as regression model.

6.	Calculate regression statistics and check regression diagnostics carefully as in Bivariate Regression Protocol (steps 5-7).  (N.B. Durbin-Watson test and influence statistics are not applicable to nonlinear regression).

7.	Calculate and plot confidence and prediction intervals as required (not applicable to nonlinear regression).


bivar Polynomial Regression pg1
POLYNOMIAL REGRESSION	page 1 of 3

  This type of regression analysis is simply a multiple regression of y against a suitable set of powers of x.  So users with no knowledge of multiple regression should view Multiple Regression before attempting polynomial regression.  The basic polynomial equation is

which is simply a logical extension of the normal bivariate regression equation.  There are still only two variables as [image], [image] and [image] are all transformations of x.  Firstly a set of transformed x variables are calculated, these can be called predictors.  Then the best set of predictors is chosen by a suitable variable selection procedure.  It is rarely neccessary to exceed [image] and beyond this collinearity is often a problem.

  The BREG command in Minitab is one of the best selection methods for polynomial regression.  The best set of predictors is simply indicated by the lowest value for Mallow's Cp, which is normally very close to the number of predictors.


Bivar Poly.Reg Mtb pg2
POLYNOMIAL REGRESSION	page 2 of 3

Minitab Commands and Output

For the Enzyme data

BREG  'vtreat'  's'  's^2'  's^3'  's^4'  's^5'

* ERROR * Predictor colums are highly correlated.

* ERROR * Use REGR command to find correlated vars.

* ERROR * Completion of computation impossible

BREG  'vtreat'  's'  's^2'  's^3'

Best  Subsets Regression of vtreat

		Adj			s s    ^ ^
Vars	R-sq	R-sq	C-p	s	s  2  3
1 1 2 2 3	69.1 48.3 90.8 88.1 95.9	66.0 43.1 88.7 85.4 94.4	53.0 93.9 12.2 17.5  4.0	30.898 39.944 17.808 20.237 12.514	x    x x  x x     x x  x  x

With the puromycin treated enzyme data a three predictor model is clearly optimal, although in this case collinearity prevented predictors past [image] being used.


Bivar Poly.Reg Mtb pg3
POLYNOMIAL REGRESSION	page 3 of 3

Minitab Commands and Output

  After selecting the predictors the regression constants are calculated, they give the following prediction equation

vtreat = 56.28 + 0.717s - 0.00116 [image] = 0.00000057

  This provides a resonable fit, particularly at lower concentrations.

Reaction velocity c.p.m.

With this solution all predictors were significant, [image] = 96.0 and the residuals were random, with no outliers nor important influence points.


 (DW = 2.2; dU - dL = 0.824 - 1.320, n = 9, k = 1).  Despite this two of the three influence statistics highlighted both the extreme values as possible influence points and the standardized residual of the 100% mortality point was >2.0.  Hence, it is possible that the calculated values at these extemities were poor approximations and so were unreasonably affecting the regression.  This problem has prompted some authors to recommend careful comparison of regressions, with and without the incl is best treated in two parts.


Bivar Transformations for 1 bend curves pg1
TRANSFORMATIONS FOR 'ONE-BEND' CURVES	page 1of 2

  'One-bend' curves are simple, with essentially only one curve (e.g. Enzyme data).  Whilst 'two-bend' curves are more complex, showing two distinct parts to the curve, and are often sigmoidal in shape.

  The use of transformations for the linearization of curved relationships has been discussed frequently and many different approaches have been suggested ( [image]).


Bivar Transformations for 1 bend curves pg2
TRANSFORMATIONS FOR 'ONE-BEND' CURVES	page 2of 2

  The series of power transformations from the ladder of powers provide an excellent series for 'one bend' curves.  Either y or x or both can be transformed as necessary.  It is not a trivial problem to select the best transformations for linearization and all the effective approaches take considerable time and effort.

  Whichever method is used, plots of the transformations must be examined to check their efficacy.

1.	Theoretical considerations

2.	Curve shape

3.	Correlation matrix method


Bivar Theoretical Considerations
THEORETICAL CONSIDERATIONS

  The analyst might know that a certain transformation can be expected to give a straight line from biological theory.  For example, populations grow logarithmically and so a plot of change in numbers of an animal (y) with time (x) might be best analysed as log10y against x whilst the population is growing exponentially.

  The enzyme data would be expected to follow Michaelis-Menten kinetics; hence the best simple pair of power transformations for linearization should be the reciprocals of y and x as this is the well known Lineweaver-Burk plot (1/v against 1/s ; [image]).


Bivar Curve Shape pg1
CURVE SHAPE	page 1 of 2

  This second approach aims to guess the power required from the shape of the curve.  However, a great deal of experience is required before powers can be guessed effectively.  Hence, for the beginner this normally results in large numbers of plots having to be examined before a good straight line is obtained.

  In order to select possible powers a large arrow is drawn on the plot pointing in the direction of the bulge in the curve.  Two vectors are then drawn against the diagonal arrow.  One vector is horizontal, parallel to the x axis, and one is vertical, parallel to the y axis.  When the vectors point to increasing values a higher power transformation is suggested and vice versa.

  Appropriate guesses of powers can then be made and the resulting plot examined.  If it is still markedly curved the steps are repeated until an approximately linear graph results.  It should be noted that this is a crude method and provides only an indication of what powers might be suitable.


Bivar Curve Shape pg2
CURVE SHAPE 	page 2 of 2

  Thus for a curve like the puromycin treated enzyme data the arrow and vectors indicate increasing the power for y and decreasing it for x.(A)

B,C and D give the rules for other shaped curves.

fig14???


bivar Correlation Matrix Method pg1
CORRELATION MATRIX METHOD	page 1 of 3

  This involves calculating a correlation matrix of the possible power transformations of both y and x.  The largest absolute correlation coefficient should give the best pair of transformations.  This will be illustrated with the puromycin treated enzyme data.  The correlation matrices from the full data with replicate y values and from the data with only the mean y values are given in the Table.

  The full data gives lower absolute values for the correlation coefficients because of the extra variability introduced by the replicate y values.  The highest correlation coefficient for the full data was for y against logx (0.982) and for the mean y data was [image] against 1/x (1.000).  Note that this latter transformation pair is very similar to the best theoretical simple power transformation (1/ y against 1/ x, the Lineweaver-Burk plot), as are the others with the very highest absolute correlation coefficients (>0.995 or <0.995) in the mean y correlation matrix.


Bivar Corr Matrix Table
Correllation Matrices for the Power Transformations in the Ladder of Powers for y and x (using the Enzyme data)

Power and transformation for  y			Power and transformation for
			3	2	1	0.5	0	-0.5	-1	-2
Power		Trans-forma-tion			x		log10x	1/	1/x	1/
Full data(n=12)
	3		0.734	0.807	0.922	0.973	0.962	-0.858	-0.730	-0.585
	2		0.690	0.764	0.891	0.961	0.981	-0.910	-0.801	-0.668
	1	y	0.623	0.695	0.831	0.922	0.982	-0.956	-0.878	-0.766
	0.5		0.578	0.647	0.786	0.887	0.968	-0.968	-0.910	-0.813
	0	log10y	0.525	0.591	0.729	0.839	0.941	-0.968	-0.930	-0.852
	-0.5	1/	-0.468	-0.529	-0.664	-0.779	-0.898	0.952	0.936	0.878
	-1	1/y	-0.409	-0.465	-0.595	-0.712	-0.845	0.921	0.925	0.886
	-2	1/	-0.308	-0.353	-0.467	-0.582	-0.727	0.833	0.868	0.861
Mean y data (n = 6)
	3		0.739	0.812	0.928	0.980	0.969	-0.866	-0.737	-0.592
	2		0.694	0.768	0.896	0.967	0.989	-0.918	-0.809	-0.675
	1	y	0.630	0.703	0.841	0.933	0.994	-0.967	-0.888	-0.775
	0.5		0.590	0.661	0.802	0.904	0.987	-0.985	-0.924	-0.825
	0	log10y	0.546	0.615	0.757	0.869	0.971	-0.996	-0.955	-0.872
	-0.5		-0.500	-0.566	-0.708	-0.827	-0.949	0.999	0.978	0.912
	-1	1/y	-0.455	-0.517	-0.657	-0.782	-0.920	0.994	0.992	0.944
	-2	1/	-0.376	-0.431	-0.565	-0.696	-0.858	0.969	1.000	0.982


Bivar Correlation Matrix Method pg2
CORRELATION MATRIX METHOD	page 2 of 3

  Inspection of the plots (using the full data) for the two pairs of 'best' transformations from the correlation matrices shows a great difference in the spread of variability

  The replicates in (A) are much more evenly spaced than in (B).  This is because use of the whole data set suggested transformations which were a compromise between even variability of replicates and linearity, whilst with the mean y data the suggested transformation ignored variability in favour of linearity.


Bivar Correlation Matrix Method pg 3
CORRELATION MATRIX METHOD	page 3 of 3

  It is well known that the Lineweaver-Burk transformation (C) emphasises variability at low substrate concentration and minimises variability at high substrate concentration ( [image]).  So it is not suprising that the similar transformation suggested from the mean y correlation matrix did the same.

  These general points should be carefully considered when selecting transformations with replicated data.  The ideal compromise would be a transformation that gave a straight line and had even variability.  The Hane's plot used earlier in the comparison of lines example provides this compromise and is clearly satisfactory (D).  However, it is not a simple power transformation and so could not have been found by any of the methods described in this section.


Bivar Transfomations for 2 bend curves pg1
TRANSFORMATIONS FOR 'TWO-BEND' CURVES	page 1 of 2

  When there are two bends in a curve the power series will not linearize the relationship and special transformations are needed.  Mortality, survival and dose-reponse curves are often sigmoidal.  In these mortality, survival or response is on the y axis and concentration, dose or time is on the x axis.  These are typical 'two-bend' curves and the following Rotenone data ( [image]) will be used as an example in this section.

rotenone concentration (mg/l)	mortality %	This data set is a typical dose-response
12.7 10.2 9.0 7.7 6.4 5.1 3.8 2.6 1.2 0.0	100 88 92 86 66 52 33 12 0 0	experiment often used in pharmacology, toxicology and other branches of biolgy. The data are derived from an experiment to test the effect of the pesticide, rotenone, the chrysanthemum aphis, in batches of 50 insects for each of the pesticide concentrations used.


Bivar Transfomations for 2 bend curves pg2
TRANSFORMATIONS FOR 'TWO-BEND' CURVES	page 2 of 2

  Apart from the commonly used 'two-bend' transformations below (and in Curved Protocol ) the arcsine transformation (One-way Transformations) can also be used to linearize 'two-bend' curves ( [image]).

Logit method

A generalized 'two-bend' transformation

Probit analysis

Calculating a predicted x ( [image])from y and its confidence intervals

Comparison of the logit and probit methods


Bivar Logit Method pg1
LOGIT METHOD	page1 of 5

  The logit transformation is derived from the logistic equation ( [image]) and essentially fits a logistic model to the data.  To do this with the rotenone data the log10 of concentration is used as the x variable and the logit of the proportion killed (y) is used as the y axis where

  With this type of data the mortality values at the extremes are often poorly estimated by experiment; because the tails of the distribution are very flat and so it is hard to get accurate values when relatively small numbers of animals are tested.  For this reason, the mortality values at 0% and 100%, which are closest to the centre of the data should be recalculated ( [image])


Bivar Logit Method pg2
LOGIT METHOD	page 2 of 5

For the rotenone data

Rotenone concentration (mg/l)	Percenatage mortality	Proportion killed'	Log10  rotenone concentration	Logit proportion killed
12.7 10.2   9.0   7.7   6.4   5.1   3.8   2.6   1.2   0.0	100   88   92   86   66   52   33   12     0     0	0.99' 0.88 0.92 0.86 0.66 0.52 0.33 0.12 0.01' 0.00	1.104 1.009 0.954 0.886 0.806 0.708 0.580 0.415 0.079 *	4.595 1.992 2.442 1.815 0.663 0.080 -0.708 -1.992 -4.595 *

' To include extreme data points (0% and 100% values for percentage mortality) the proportion killed has been recalculated (1) as 1/2n1   for the 0% point and 1-1/2n1   for the 100% point.  Where n1 is the number of individuals tested at the maximum concentration for which the proportion killed is being recalculated.  A maximum of one 0% point and one 100% point should be included.

The control data without rotenone cannot be included as log100.0 is incalculable


Bivar Logit Method pg3
LOGIT METHOD	page 3 of 5

Minitb Commands and Output

Regression of the logit transformed rotenone data

REGRESS 'logit-m' 1 'lg-conc'.

The regression equation is

logit-m = -5.46 + 8.17 lg-conc

Predictor Constant lg-conc	Coef -5.4596 8.1685	Stdev 0.4629 0.5870	t-ratio -11.79 13.92	p 0.000 0.000
s = 0.5391   R-sq = 96.5%   R-sq(adj) = 96.0%

Analysis of Variance

SOURCE Regression Error Total	DF 1 7 8	SS 56.296 2.035 58.330	MS 56.296 0.291	F 193.67	p 0.000
  The linear model fits these data fairly well as [image] was high and both t ratios and the F value were highly significant; also there was no significant autocorrelation in the residuals (DW = 2.2; dU - dL = 0.824 - 1.320, n = 9, k = 1).  Despite this two of the three influence statistics highlighted both the extreme values as possible influence points and the standardized residual of the 100% mortality point was >2.0.  Hence, it is possible that the calculated values at these extemities were poor approximations and so were unreasonably affecting the regression.


Bivar Logit Method pg4
LOGIT METHOD	page 4 of 5

  This problem has prompted some authors to recommend careful comparison of regressions, with and without the inclusion of the exteme points, before any decision to use them is made.  The regression without the extreme points gave similar results, but in this case the 88% mortality point was identified as a possible influence point.  This was not unexpected as it looks a bit lower than the other points on the plotted data.

Logit proportion killed


Bivar Logit Method pg5
LOGIT METHOD	page 5 of 5

  Some authors recommend weighted least squares regression with the logit method as the variances are not equal at all x values ( [image]).  If required this is simply done by using a column of weights (w), one for each data point.  These are calculated as

w = n1 y(1 - y)

where n1 is the number of individuals tested for that data point and y is the proportion killed.

.  In Minitab this type of regression is done by referring to the column containing the weights in the WEIGHTS subcommand of the REGRESSION statement.


Bivar A Generalized 2 bend Transformation
A GENERALIZED 'TWO-BEND' TRANSFORMATION

  There is a widely suggested, more universally applicable transformation for 'two-bend' curves.  This is similar to the logit method and is defined as

where K is unknown and must be estimated before the transformation can be applied.  This is best done by nonlinear regression.  However, plotting the transformed y against x with a variety of values for K will also prove satisfactory.  With practice K can be estimated crudely from the shape of the untransformed relationship.


Bivar Probit Analysis pg1
PROBIT ANALYSIS	page 1 of 2

  As the culmulative curve of the normal distribution is sigmoidal in shape, normal density functions can be used as the basis for another 'two-bend' linearizing transformation.  The popular probit transformation is based on this idea because probits are normal equivalent deviates with a mean of 5.0 and a standard deviation of 1.0.  Finney ( [image]) extensively discusses the application of the probit method to dose-response curves and these texts should be consulted for more details.


Bivar Probit Analysis pg 2
PROBIT ANALYSIS	page 2 of 2

Minitab Commands and Output

  The Minitab command to obtain probits is

INVCDF  C3  C4;
NORMAL  5  1.

where C3 contains the values for the proportion killed and C4 is where the probits for these proportions are stored.  The same rules for extreme values as suggested for the logit method can also be used here and these values are also recalculated.(see Table)

  The regression is then straightforward with the probits as the y variable and log10 rotenone concentrations as the x variable.  Some authors have suggested using a minimized [image] regression technique in place of least squares for probit analysis.  However, the three-group resistant line method would almost certainly be better as a robust technique to minimize the effect of outliers (see Bivariate Regression Protocol).


Bivar Probit Analysis table
PROBIT ANALYSIS Table

Rotenone concentration (mg/l)	Percenatage mortality	Proportion killed'	Log10  rotenone concentration	Probit proportion killed
12.7 10.2   9.0   7.7   6.4   5.1   3.8   2.6   1.2   0.0	100   88   92   86   66   52   33   12     0     0	0.99' 0.88 0.92 0.86 0.66 0.52 0.33 0.12 0.01' 0.00	1.104 1.009 0.954 0.886 0.806 0.708 0.580 0.415 0.079 *	7.326 6.175 6.405 6.080 5.412 5.050 4.560 3.825 2.674 *

' To include extreme data points (0% and 100% values for percentage mortality) the proportion killed has been recalculated (1) as 1/2n1   for the 0% point and 1-1/2n1   for the 100% point.  Where n1 is the number of individuals tested at the maximum concentration for which the proportion killed is being recalculated.  A maximumof one 0% point and one 100% point should be included.

The control data without rotenone cannot be included as log100.0 is incalculable


Bivar Calc of a Pred x from y & C.I.pg1
CALCULATING A PREDICTED X FROM Y AND ITS CONFIDENCE    page 1 of 4
INTERVALS

  Biologists often want to calculate the dose or concentration giving a 50% reponse together with the confidence interval for this value.  To do this it is necessary to predict a new value of x from a specific y value; the logit and probit values for 50% are 0.0 and 5.0 respectively.

  As this is the reverse of the normal way in which predictions are made another approach is needed.


Bvar Calc of a Pred x from y & C.I. pg2
Method for Calculating Confidence Intervals for a New  Predicted Value of x ( [image]) from a New Value of y

1.

2.	Calculate statistics used later
 mean of all the x and y values

3.	Calculate the first intermediate statistic, D

4.	Calculate the second intermediate statistic, H

or when y is known

5.	Calculate the third intermediate statistic, J

6.	Calculate confidence intervals   CIU = J + H  and  CIL = J - H

7.	Apply confidence intervals to [image]  CIL < [image] < CIU  note that the confidence intervals are not symmetrical about [image].

(based on a method given in Sokal & Rohlf; [image] ).


Bivar Calc of a Pred x from y & C.I. pg3
CALCULATING A PREDICTED X FROM Y AND ITS CONFIDENCE    page 3 of 4
INTERVALS

  As an example of this method the calculations can be applied to the regression of the rotenone data using the logit method with the inclusion of the recalculated extreme mortality values.  Some of the regression statistics needed for this are in the ouput obtained when regression of the logit transformed data is carried out.  In this case [image] is the lethal concentration for  50% of the aphids (LC50) and the confidence intervals calculated will be the 95% intervals for the LC50.  Hence, we need to find [image] when logit y = 0.0 (mortality = 50%).

  First the predicted x value, the mean x and y values and the [image], which for these data are

= (0.0 + 5.460)/8.169 = 0.668

 = 0.477

 = 0.727

 = 0.8437


Bivar Calc of a pred x from y & C.I. pg 4
CALCULATING A PREDICTED X FROM Y AND ITS CONFIDENCE   page 4 of 4
INTERVALS

Next the three intermediate statistics, D, H and J, must be calculated sequentially

D  =

the LC50 is known by definition so in this case we use the second formula for H in step 4

So the confidence interval for [image] (0.668) is  0.613  <  0.668  <  0.721

as CIL = 0.667 - 0.0539 = 0.613 and CIU = 0.667 + 0.0539 = 0.721; which converted into the original units is

4.10  <  4.66  <  5.26


Bivar Comparison of Logit & Probit pg1
COMPARISON OF THE LOGIT AND PROBIT METHODS	page 1 of 2

  The logit and probit methods normally produce very similar results.  They both work well for dose-response type of data and other data with a symmetrical 'two-bend' curve.  Finney ( [image]) discusses these methods in relation to the analysis of the age of the onset of menstruation with data from 3918 Warsaw girls and concludes that the models fit the data well and give almost identical results.  The table summarizes regression results from both methods with and without the recalculation of extreme points.

Regression method'	Extreme values''	bo	b	(%)	Most extreme standardized residual'''	LC50 (mg/l)
Logit	+ -	-5.46 -5.07	8.17 7.43	96.6 96.5	2.27 -1.63	0.668 0.683
Probit	+ -	2.12  2.01	4.35 4.38	97.4 96.9	-1.80 -1.70	0.663 0.684
Probit, robust	+ -	1.84  2.03	4.78 4.33	* *	-1.80   1.46	0.661 0.686

'  All regressions were of transformed proportions killed against the log10 rotenone concentration.  Robust indicates the three-group resistant line method, all other results were for least squares regression.

'' + = extreme values included (n=9);  - = extrame points excluded (n=7).

''' values for the three-group resistant line method were calculated by dividing the residuals by thier standard deviation.

* values not calculated by Minitab


Bivar Comparison of Logit & Probit pg2
COMPARISON OF THE LOGIT AND PROBIT METHODS	page 2 of 2

  It is clear that the two approaches give very similar results.  Only the logit method with the extreme points gave an outlier >2.0, but this is probably not a significant difference as all the other statistics were very similar.  The three-group resistant line method had only a very slight effect on the regression statistics, further indicating that neither outliers nor the extreme points have an important effect on the regression solution for these data.

  The intercept and slope values appear different, but this is due to the logits and probits having means of 0.0 and 5.0 respectively.  Although the results are not given here, all the regressions gave highly significant t ratios for slope and intercept (P < 0.001) and F values from the regression ANOVA
(F > 140; P << 0.001).  The LC50 values were consistently slightly higher for the regressions without exteme values, but this difference is probably of minor significance as the 95% confidence intervals were much larger than this difference.

  So we can conclude that for these data it makes no practical difference whether the logit or probit method is used; this conclusion is probably generally applicable to most examples of this type of data.


Bivar Direct Nonlinear solns pg1
DIRECT NONLINEAR SOLUTIONS	page 1 of 4

  Nonlinear regression provides a method of estimating a line of best fit directly from accepted, biologically useful equations ( [image]).  This means that the empirical regression equations do not have to be used and biologically important constants can be directly estimated.  The disadvantages of nonlinear regression are that the calculation of confidence intervals and regression diagnostics is not as well developed as for the linear, least squares methods.

  Systat, SPSS/PC+ and Statgraphics can all do nonlinear regression but Minitab cannot.

  Using the Enzyme data the variable names will be as used previously and the constant names will be VMAX = Vmax = maximum reaction velocity and KM = Km = Michaelis or half saturation constant.  The Michaelis-Menten equation

where v = reaction velocity (vtreat) must be used as the regression model.


Bivar Direct Nonlinear solns pg2
DIRECT NONLINEAR SOLUTIONS	page 2 of 4

  The commands needed to execute this nonlinear regression in Systat are

NONLIN

MODEL VTREAT = (VMAX*S)/(KM+S)

ESTIMATE / QUASI

and for SPSS/PC+ are

MODEL PROGRAM VMAX=1 KM=1.

COMPUTE PRED=(VMAX*S)/(KM+S).

NLR VTREAT WITH S.

  Note that for Systat there is no LOSS statement, so the loss function will default to least squares, whilst with SPSS/PC+ only least squares is available.  Also starting values of 1.0 are used by default in Systat but must always be defined in SPSS/PC+.

  Sometimes more realistic starting values are required and should be roughly estimated directly from the data.  In this case they could be estimated directly from the simple plots of the data, sensible guesses being Vmax = 200 and Km = 100; but the use of starting values other than 1.0 did not effect the answers here.


Bivar Direct Nonlinear solns pg3
DIRECT NONLINEAR SOLUTIONS	page 3 of 4

  The calculation procedures are iterative and with both packages the key results at each step are printed, followed by some statistics for the final solution.  These statistics are an ANOVA table, [image] value and the estimates of VMAX and KM with their standard errors and 95% confidence intervals.

The SPSS/PC+ output using the larger starting values

Iteration	Residual SS	VMAX	KM
1 1.1 2 2.1 3 3.1 4 4.1 5 5.1 6 6.1	7964.185201 1593.159441 1593.159441 1201.034526 1201.034526 1195.508871 1195.508871 1195.449383 1195.449383 1195.448820 1195.448820 1195.448814	200.000000 212.023790 212.023790 211.772798 211.772798 212.563259 212.563259 212.671588 212.671588 212.682567 212.682567 212.683630	100.000000 54.2873641 54.2873641 62.3244628 62.3244628 63.9264826 63.9264826 64.1022790 64.1022790 64.1194489 64.1194489 64.1211054
Run stopped after 12 model evaluations and 6 derivative evaluations. Iterations have been stopped because the relative reduction between successive residual sums of squares is at most SSCON = 1.00E-08

Nonlinear Regression Summary Statistics   Dependent Variable VTREAT

Source	DF	Sum of Squares	Mean Square
Regression Residual Uncorrected Total (Correctd Total)	2 10 12 11	270213.55119   1195.44881 271409.00000  308598.91667	135106.77559    119.54488
R squared =1-Residual SS/Corrected SS =  .96126

			Asymtotic 95%
Parameter	Estimate	Asymptotic Std. Error	Confidence Lower	Interval Upper
VMAX RM	212.68362994 64.121105356	6.94709058 8.280753406	197.20454734 45.670436968	228.16271255 82.571773744


Bivar Direct Nonlinear solns pg4
DIRECT NONLINEAR SOLUTIONS	page 4 of 4

  Output with starting values of 1.0 was similar but had more iterations.  As expected this solution fitted the original data well.

Reaction velocity, c.p.m.

  The alogarithms for these two packages gave slightly different results.  The constants given by Systat were VMAX = 204.9 and KM = 59.04, both being well within the SPSS/PC+ confidence limits.  Such differences between packages are not unexpected for nonlinear regression.


Bivar Best Method pg1
THE BEST METHOD	page 1 of 3

  It is clear that there are normally several approaches that can be used for regression analysis where the relationship between y and x is curved.  Choosing the best approach is not always straightforward.

A summary of the regression results for the puromycin treated enzyme data.

Regression method	(%)	Most extreme standardized residual	Correlation coefficient of normal scores and standardized residuals	Vmax (counts/ml)	Km (mg/ml)
Polynomial	95.9	-2.27	0.975*	-a	-
Lineweaver-Burk	85.6	-3.10	0.859	195.8	48.4
Hane's plot	99.8	-1.77	0.977*	216.2	67.9
v against log10s	96.5	-1.94	0.986*	-	-
against 1/s	75.3	3.12	0.802	-	-
Nonlinear by Systat	99.5	1.97	0.977*	204.9	59.0
Nonlinear by SPSS/PC+	99.3	2.36	0.959*	212.7	64.1
 {v = reaction velocity, vtreat; s = substrate concentration}

 a not easily calculable from these regressions


uncorrected value given; calculated from sums of squares, where

Note that because y varies in the different regression methods the [image] values are not strictly comparable.


values for nonlinear regression were calculated by dividing the ordinary residuals by their standard deviation; the ordinary residuals themselves were calculated from the predicted values of v, obtained from the nonlinear regression estimates of Vmax and Km and the Michaelis-Menten equation.


critical value = 0.927; * = residuals normally distributed, P < 0.05;

(see One-way ANOVA for further details).


obtained by multiple regression of v against s, [image] and [image].


1/v against 1/s; Vmax and Km calculated from slope and intercept;

( [image])


s/v against s; Vmax and Km calculated from slope and intercept,
 ( [image])


Bivar Best Method pg2
THE BEST METHOD	page 2 of 3

  It is clear that the Lineweaver-Burk method and [image] against 1/s were the least satisfactory; they both had the lowest [image] values, strongest outliers and did not have normally distributed residuals.  All the other regressions gave acceptable and fairly similar results.  So if a purely empirical prediction equation was required they would all be satisfactory.

  However, if values of the biologically important constants (Vmax and Km in this example) were required, only the nonlinear methods or Hane's plot could be used.  Neither the polynomial solution nor v against log10s could provide estimates of these constants.  Hane's plot depends on a knowledge of the rearranged Michaelis-Menten equation, so it is likely that nonlinear regression provides the most generally useful approach when the underlying biological equation is known.


Bivar Best Method pg3
THE BEST METHOD	page 3 of 3

  When a suitable equation is not known, it is likely that polynomial regression would be both easier and less prone to problems, such as uneven variances, than the methods which select suitable power transformations by examining curve shape or by studying correlation coefficients.  However, polynomial regression does have its own problems.

  The polynomial solution gives an unlikely dip in the fitted line between the data points at substrate concentrations of 560 and 1100 mg/ ml.  Also, the y values below the lowest substrate concentration soon became negative and those over the highest concentration rapidly became unnaturally high.

  Such problems can be minimized by using more evenly spaced x values and never extrapolating outside the limits of the data.

  Users should not be tempted to try to find an arbitary transformation of the data where the dependent variable, y, is calculated from the predictor variable, x.  When this is done it can sometimes happen that an apparent, but spurious, relationship can be demonstrated purely because y was calculated from x.  The Hane's plot used here is satisfactory because it is derived from a good theoretical basis.


Polynomial solution showing an 'unlikely dip'

Reaction velocity c.p.m.


Factorial Intro
FACTORIAL ANALYSIS OF VARIANCE (ANOVA)

  Experiments where the response is measured at all combinations of levels of several factors are called factorial experiments.  The value of such experiments is in the opportunity of assessing the inter-relationships or interactions between factors;  it is possible that the response to one factor may differ according to the levels of a second factor.

  Also, a single well-designed factorial experiment is usually more efficient than a series of single factor experiments.

  An important pre-requisite of any factorial experiment is a consideration of the experimental design so that the experiment satisfies the objectives of the investigation, and also that the experiment is efficient in its use of information.  To make the best possible use of experimental resources the design should be considered before data one collected and analysis attempted.

Experimental Design


Factorial Experimental Design
EXPERIMENTAL DESIGN

The following list gives some guidance on the considerations of experimental design.

1.	State the objectives of the experiment and make sure that the experiment will answer the objectives.

2.	Consider all available relevant information, using a pilot experiment where necessary.

3.	Choose the response variable, the observation, so that it genuinely provides useful information.

4.	Choose the factors to investigate so that the experiment is representative and informative, paying particular attention to the control of the levels of the factors.

5.	Consider the nature of the experimental units or plots carefully, ensuring these are representing.  What out for split plots.

6.	Ensure that the plots are homogeneous and, if necessary, divide them into blocks and include blocks as a factor in the analysis.

7.	Make certain that sufficient replicates are used to assess all effects and interactions, using prior information on variability.

8.	Allocate combinations of factors to plots at random, within the constraints of the design, so as to eliminate bias.

9.	Monitor the conduct of the experiment and keep a detailed accurate record of all relevant facts.

10.	Use statistical analysis to show objectively those differences between factors that are genuinely significant and those that are merely random fluctuations.


Factorial Data Structure pg 1
FACTORIAL ANOVA DATA STRUCTURE

  The two forms of Factorial ANOVA data structure are Crossed and Hierarchical (nested).

The Crossed layout:

  Each observation is represented by an x and here there are two replicated data for each combination of levels of the two factors one combination of levels, called a cell, is indicated with a box.


Each measurement of a response at the indicated levels of the two factors needs an amount of experimental material, often called a plot or experimental unit.  Two replicate data in a cell require two independent plots to be made available for each cell.  A genuine replicate cannot be obtained by simply duplicating the measurement made on a single plot.  In some factorial experiments it may be expedient to have more than two replicates in each cell.


Factorial Data Structure pg 2
FACTORIAL ANOVA   DATA STRUCTURE

The hierarchical layout

The distinction between a crossed layout and an hierarchical layout is most easily clarified by an example.

Example Data Sets

Balanced Experiments


Factorial Differences pg 1
DIFFERENCE BETWEEN CROSSED AND HIERARCHICAL LAYOUT	   page 1 of 2

Yields (tonnes per hectare) of four varieties of wheat at three farms.

	Variety 1	Variety 2	Variety 3	Variety 4
Farm 1 Farm 2 Farm 3	0.327      0.280 0.532      0.526 0.269      0.277	0 500      0.510 0.599      0.637 0.308      0.268	0.442      0.360 0.516      0.499 0.241      0.328	0.471      0.460 0.638      0.455 0.305      0.301

  So as to obtain enough replicate data to enable comparisons to be made three farms were used in the experiment.  At each farm eight similar plots of land were selected, each variety was sown on two replicate plots.

This experiment is described as crossed because each variety is grown at each farm.

  All the data in the first row have in common that they are measurements from farm 1.  Similarly the data in the first two columns have in common that they are measurements of variety 1.  The three farm averages and the four variety averages can be compared to determine if there are differences between the three farms or the four varieties.


Factorial differences
DIFFERENCE BETWEEN CROSSED AND HIERARCHICAL LAYOUTS	page 2 of 2

  Suppose the experiment had been conducted differently, although still comparing four varieties by measuring the yields from 24 plots of land.

  If the farms listed for variety 1 differ from those for variety 2 and those for variety 3 and 4 then the experiment is described as hierarchical or nested, with differences between farms nested within varieties.

  The comparison of the four variety means is still made, averaging each variety across farms.  However, there is no meaning in the comparison of the three row averages since the farms differ from variety to variety.  Instead, the farm averages (of two data) are compared within varieties.

  In this context the crossed experiment is preferred.  However, in other contexts hierarchical experiments give useful information.

Conventions


Although it is a convention to arrange the data as seen, the layout can readily be re-arranged to give the appearance of either type of experiment, hierarchical or crossed.  Thus care needs to be taken when deciding which form the experiment takes.  It is the inter-relationships of the factors themselves that determine whether they are crossed or hierarchical.


Factorial balanced experiments
BALANCED EXPERIMENTS

  In factorial experiments with two or more factors it is an advantage to have a balanced experiment.  In such experiments the number of levels of one factor for which measurements are made is the same at every level of all other factors, and the number of replicate data is the same for all cells.  There are three advantages of balanced experiments in preference to unbalanced ones.

 they are easier to analyse

 they are usually much easier to interpret

 they are more robust to departures from statistical assumptions used to justify the analysis.


Factorial example data sets pg 1
EXAMPLE DATA SETS	page 1 of 2

Three sets of data are used for illustration of ANOVA techniques.

  The first data set is the one used to describe the differences between the two types of data layout (wheat data).  The second data set is the experiment describing the numbers of seedlings produced from loblolly pine seeds.

	Number of seedlings on
Sowing date  in march	Burnt ground		Conventionally  cleared ground
2 6 10 14 18 22	900' 760 880' 1060 1530 1390' 1970' 1820 1960 1310 830 570	810 1040 1170 910' 1160 1540 1890 2140' 1670 1480' 420 760	880 960' 1050' 1110 1140 1320 1360' 1490' 1270 1500 150 420	1100 1040 1240 1120 1270 1080 1510 1270 1380 1450 380' 270	observations marked with a ( ' )  are assumed to be missing  in some examples
  A fixed number of seeds were sown randomly on six dates in spring into four replicate plots on land cleared by burning and, on the same dates, in another four plots cleared by conventional means.


Factorial example data sets pg 2
EXAMPLE DATA SET	page 2 of 2

  The data for the third experiment were obtained by students in Cardiff.  Twelve students were involved, six postgraduates (PG) and six undergraduates (UG).  Each student made three replicate determinations of the number of viable heterotrophic bacteria present in water samples (by a plate count determination).

Water type	Student	Postgraduates	Undergraduates
Sewage	1 2	2700   2800   1700 2600   3000   3200	4900   3600   1300 3500   4100   2500
Polluted River	1 2	52      49       61 68      75       83	60     140       75 40       60       80
Clean Stream	1 2	5.9     7.6    16.0 5.6     5.9      6.3	6.4      4.2      4.1 8.3      3.1      5.2
  Four students, two undergraduates and two postgraduates, worked on water from sewage, four worked on water from a polluted river and four on water from a clean stream.  Notice that the students have been labelled 1 or 2 in pairs for each water source/PG-UG combination.  Thus convention for labelling is the one used in texts on ANOVA, and also for data entry in statistical computer packages.

Catagorical data

Sample population data


The analysis of variance is not suitable for the analysis of categorical data where a population is categorised, possibly by several classifications simultaneously, and the number in a sample from the population are counted in each category.  A [image] contingency table may be suitable for the analysis of such data.


The analysis of variance is not suitable for the analysis of experiment where measurements of a number of variables are made on a sample from a population.  Such data may be analysed by a method of multivariate analysis.


General Principles of Factorial ANOVA
GENERAL PRINCIPLES OF ANOVA

  The idea of ANOVA is that the overall variability of the data, measured by a sum of squares (of deviations of the data from the overall mean), can be partitioned into several components.  Associated with each sum of square is a number called the degrees of freedom, df.  A sum of squares is divided by its degrees of freedom to give a mean square, and it is these that are tested against each other to determine the significance of effects.  The mean squares also enable variances to be estimated for random effects.

  For balanced experiments both Minitab and SAS allow the expected mean squares column to be printed.  For unbalanced experiments, however, the calculations are computationally more difficult and can be time consuming.  Only SAS evaluates the expected mean squares column for unbalanced data.


Route trough fact.Anova
ROUTE THROUGH A FACTORIAL ANOVA

1.	Enter the data for computer solutions

2.	Compute the factorial ANOVA table

3.	Check the assumptions of the factorial ANOVA

4.	Transform daa if assumptions do not hold

5.	Test for outliers

6.	Comparison of meansof main effects

For more details click


Fact ANOVA Protocol
FACTORIAL ANOVA PROTOCOL

1.	Enter data for computer solutions of factorial ANOVA

2.	Compute the factorial ANOVA table , this can be done in Minitab by a one-way ANOVA with cells as groups

3.	Check assumptions of factorial ANOVA

(a)  Tests for homogeneity of variances

 Bartlett's

 Hartley's F max

 Cochran's C

(b)  Tests for normality of residuals

 normal probability plots

 Shapiro-Wilk test

 boxplots

4.	Transform data if assumptions do not hold and repeat the checks

5.	Examine residuals of transformed data and test for outliers and for patterns in the residuals.

 Dixon's test (n < 25)

 Grubb's test (n > 25)

6.	Comparison of means of main effects

(a)  A posteriori tests

(b)  A priori tests


Factorial data entry
1. ENTERING THE DATA

  The form in which data have to be entered for statistical packages is not usually that of a table as seen on the example data sets.

  In Minitab, SPSS/PC+ and SAS all the data are entered in one column.  A further column is needed for each factor in the model giving the levels of the factors.  These columns are called indicator (labelling or subscript) variables.


Factorial Expl data entry pg 1
EXAMPLE DATA ENTRY	page 1 of 2

Row	Yield	Farm	Variety
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24	0.327 0.280 0.00 0.510 0.442 0.360 0.471 0.460 0.532 0.526 0.599 0.637 0.516 0.499 0.638 0.455 0.269 0.277 0.308 0.286 0.241 0.328 0.305 0.314	1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3	1 1 2 2 3 3 4 4 1 1 2 2 3 3 4 4 1 1 2 2 3 3 4 4	Here the four types of variety are 'sub-labelled within' the three farms.  Data entry for the loblolly seeds is essentially the same with the six dates 'sub-labelled within' the type of ground (burnt/conventionally cleared).


Factorial Exple data entry pg 2
EXAMPLE DATA ENTRY	page 2 of 2

Row	Count	Waters	Students	PG-UG
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36		1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3	1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 2 2 2 2 2	1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2 1 1 1 2 2 2	Here the heterotroph data shows the convention for the numbering of the nested students factor.  Although in the experiment there were twelve students, the numbering for entry into the computer is the students within water types, hence 1 and 2 are repeated six times.


Factorial Computing Fact table
2. COMPUTING THE FACTORIAL ANOVA TABLE

  Before computing the ANOVA table it must be decided which factors are fixed and which are random.  Any random factors can then be specified using a Minitab subcommand.

Fixed effects models

Random effects models

Minitab commands and output for:

Balanced experiments with two factors

1.  Two crossed factors, no replicates

2.  Two crossed factors with replicates

3.  Two factors with one hierarchically nested with the other

Balanced experiments with more than two factors

Unbalanced experiments


Factorial fixed effects models
FIXED EFFECTS MODELS

  Fixed effects models are appropriate for the analysis of experiments where the groups represent specific levels of a factor, determined in advance by the experiment.

  The detection of significant differences in the means for the levels (groups) and, if identified, the isolation of those levels whose means differ significantly from the others is the appropriate approach for analysis.


Factorial random effect models
RANDOM EFFECTS MODELS

  Random effects models are appropriate for experiments in which the levels of the factor represent different individuals in some populations, either finite or infinite, which are sampled by a random selection.  The means for the particular levels observed would not be of any interest, if the experiment were repeated different levels would be observed.  (In the fixed effects case the same levels would be observed in any repetition of the experiment).

  Random effects models are often wrongly used.  They should be reserved for those experiments in which the identification of the levels of a factor involves random sampling from some population.  Deciding the levels of measurement of a factor, or those treatments to use as levels of a factor in an experiment, does not represent random sampling and such factors should be regarded as fixed effects.


Factorial 2 crossed, no replicates fig 1
COMPUTING THE FACTORIAL ANOVA TABLE	page 1 of 3

Two Crossed Factors, no replicates

  This is the simplest multi-factorial experiment demonstrated by a modification of the wheat data.  In the absence of replications it is not possible to assess whether the effect of factor B is consistent from one level of factor A to another.  Such an inconsistency is called an interaction between factor A and B, and the assessment of interactions is done only when there are replicates.

  The absence of replicate data for any cells also precludes the calculation of cell variances, so that diagnostic checking for homogeneity of variance can only be done crudely by an examination of the residuals.  If they increase markedly with the levels of one or other of the factors, then the homogeneity of variance assumption may be called into question.

  Two factor experiments with no replicates are limited therefore and in general it is advisable to replicate measurements wherever this is possible.


Yields (tonnes per hectare) of four varieties of wheat at three farms.

	Variety 1	Variety 2	Variety 3	Variety 4
Farm 1 Farm 2 Farm 3	0.327      0.280 0.532      0.526 0.269      0.277	0.500      0.510 0.599      0.637 0.308      0.268	0.442      0.360 0.516      0.499 0.241      0.328	0.471      0.460 0.638      0.455 0.305      0.301

Assume that only the first of each pair of data in each cell is available


Factorial 2 crossed, no replicates pg 2
COMPUTING THE FACTORIAL ANOVA TABLE	page 2 of 3

Two Crossed Factors, no replicates

  A particular case of the two-Factor crossed experiment with no replicates is the randomized block experiment.  There is also a non-parametric equivalent; called Friedman's test, based on the ranks of the data within the blocks.

  For the wheat data, the factor representing varieties is regarded as fixed but the interpretation of the farms factor is open to question.  In this case it will be assumed that the farms are random.  (The farms factor could have been included as blocks, and the experiment is then a randomised block).


Randomized block experiment

  Here a treatments (or levels of a factor A) are compared with b replicates per level are required.  This requires ab observations to be made.  The randomized block experiment is used where the experimental material or plots on which measurements are made can be naturally divided into b sets called blocks, each block containing a plots.  the purpose in designing the experiment in this way is that differences between blocks can be allowed for in assessing the effect of factor A.  In each block one plot is used for each level of factory A.  the efficiency of the design depends on the plots within the blocks being roughly similar.


 Friedman's Test

  This test is useful if the errors in the data are so markedly different from normally distributed that no simple transformation enables a standard ANOVA to be done.

  The test is available on many statistical packages including Minitab and SPSS/PC+.  For more complicated factorial experiments non-parametric methods are still being developed.


Factorial 2 crossed, no replicates pg 3
TWO FACTORS CROSSED, NO REPLICATES	page 3 of 3

Minitab Commands and Output

For the Wheat data:

ANOVA C1 = C2 C3;	[C1 - 'yield', C2 - 'farm', C3 - 'variety']

RANDOM C2;	['farm' is a random factor]

EMS.	[gives the expected mean square]

The output appears:

Factor	Type	Levels	values
Farm	random	3	1	2	3
Variety	fixed	4	1	2	3	4
Analysis of Varience for Yield

Source	DF	SS	MS	F	P	Click on the green
Farm	2	0.168996	0.084498	56.86	0.000	'hotspots' for
Variety	3	0.021185	0.007062	4.75	0.050	explanation of
Error	6	0.008917	0.001486			output
Total	11	0.199098
Source	Variance component	Error term	Expected Mean Square (using unrestricted model)
1 Farm	0.02075	3	(3) + 4(1)
2 Variety		3	(3) + Q[2]
3 Error	0.00149		(3)

Friedman's Test, Commands and Output

SAS Commands


 The F-ratio for varieties is 4.75 (with 3 and 6 df) and this is just significant, as indicated by the p value of 0.050.  Thus a marginally significant difference in the varieties has been indicated bt these data.

  The F-ratio for farms is56.86 (with 2 and 6df) and this is highly significant, with a P value of <0.001.


 The Expected Mean Square column uses as the symbol for the error variance [image], as the symbol for the farms variance [image] and Q[2] for the quadratic expression in the adjustments for the fixed varieties effects.

  The use of the unrestricted model refers to assumptions made about a mixed model with replicate data and can be ignored here.


 The two components of variance are 0.02075 for farms and 0.00149 for error. (Notice that no variance component is attributed to varieties.  The varieties are a fixed effect and they are assessed by comparing their means).  Since a linear model is assumed, the variance of any observation is the sum of the variance components.  Thus, in this case the variability in farms repesents an estimated proportion of 0.02075/(0.02075 + 0.00149) = 0.933 of the total variance of any observation (for any particular variety).

  This division of variability into proportions attributable to the random components of the model is more easily understood than a simple reporting of the variances.


Factorial Friedman's Test
FRIEDMAN'S TEST

Minitab Commands and Output

  Using the same data for the purposes of comparison, the varieties are regarded as the factor under investigation and the farms are treated as blocks.  Each variety is measured on three plots but these plots are on different farms.  All four varieties are measured on each farm. Differences between varieties is to be investigated, making allowance for differences between the farms.  Differences between the farms are allowed for but not themselves assessed by Friedman's test.

Using the farm data:

   FRIEDMAN C1 C3 C2   	[C1 - 'yield', C2 - 'farm', C3 - 'variety']

Friedman test of Yield by Variety blocked by Farm

S = 7.40  d.f. = 3  p = 0.601

Variety	N	Est. Median	Sum of RANKS
1 2 3 4	3 3 3 3	0.40825 0.48175 0.41075 0.48925	5.0 11.0 4.0 10.0
  The P value of 0.061 associated with the [image] test of the hypothesis that the median yield is the same for each variety confirms the ANOVA result, that there is only marginally significant difference in the yields for the four varieties.

  The similarity of the median yields is evident from the column of medians.  For these data the conclusions drawn from Friedman's test are similar to those of the parametric analysis of variance.


SAS uses the VARCOMP procedure to give similar results, with the METHOD = type1 option and specifying only the class variable for farms as random on the MODEL line.


Factorial 2 crossed with replicated pg 1
COMPUTING THE FACTORIAL ANOVA TABLE	page 1 of 4

Two Crossed Factors, with replicates

  Here the interaction between A and B can be assessed.  Either or both of the factors can be regarded as random and if one is random then both the main effect for that factor and the interaction are regarded as independent random variables with associated variances.  For the mixed model where one effect is fixed and the other random two alternatives are possible;  the restricted model or the unrestricted model.  Both models give sensible statistical estimates of the parameters of the model, and it is an unresolved question of which is to be preferred.

  SAS assumes the unrestricted mode, but Minitab enables a choice to be made.  The main difference is that in the unrestricted model both effects are tested against the interaction.  In a restricted model the random effect is tested against the error.


In the restricted model, as well as assuming that the fixed effects sum to zero, it is assumed that the interaction effects sum to zero over the levels of the fixed effect.

  This assumption is not made in the unrestricted model.

The unrestricted model is the one usually assumed for unbalanced data.


Factorial 2 crossed with replicates pg 2
TWO CROSSED FACTORS WITH REPLICATES	page 2 of 4

Minitab Commands and Output

For the Wheat data:

ANOVA C1=C2 C3 C2*C3;	[C1 - 'yield', C2 - 'farm', C3 - variety']

RANDOM C2;	['farm' - random factor]

RESTRICTED;	[specifies the restricted model]

EMS.	[gives the expected mean square]

Factor	Type	Levels	values
Farm	random	3	1	2	3
Variety	fixed	4	1	2	3	4
Analysis of Varience for Yield

Source	DF	SS	MS	F	P	Click green
Farm	2	0.268861	0.134430	61.33	0.000	'hotspot' for
Variety Farm*Variety	3 6	0.038494 0.023206	0.012831 0.003868	3.32  1.76	0.098 0.190	explanation of output
Error	12	0.026305	0.002192
Total	23	0.356866
Source	Variance component	Error term	Expected Mean Square (using unrestricted model)
1 Farm	0.01653	4	(4) + 8(1)
2 Variety 3 Farm* Variety	0.00084	3 4	(4) + 2(3) + 6Q[2] (4) + 2(3)
4 Error	0.00219		(4)

So that a comparison can be made between the restricted and unrestricted models in the mixed model case, and also to make a comparison with SAS output - click to view SAS OUTPUT


 Neither the interaction F-ratio nor the variety F-ratio is significant (P = 0.190 and 0.098 respectively).  This indicates that the varieties do not differ significantly from farm to farm (the nin significant interaction) and that there is little evidence of overall significant differences between varieties.  There is, however,a significant difference within the farms (P = 0.0005).


Factorial SAS commands and Output
TWO CROSSED FACTORS WITH REPLICATES

Sas Commands and Output

For the wheat data

PROC VARCOMP METHOD = TYPE1;

CLASS FARM VARIETY;

MODEL C YIELD = VARIETY FARM VARIETY * FARM/FIXED = 1;

RUN;

Variance Components Estimation Pocedure

Dependent Variable: YIELD

Source	DF	Type I SS	Type I MS
VARIETY	3	0.03849433	0.01283144
FARM	2	0.2688610	0.13443050
FARM*VARIETY	6	0.02320567	0.00386761
Error	12	0.02630500	0.00219208
Corrected Total	23	0.35686600
Source	Expected Mean Square
VARIETY	Var(Error) + 2 Var(FARM*VARIETY) + Q(VARIETY)
FARM	Var(Error) + 2 VAR(FARM*VARIETY) + 8 VAR(FARM)
FARM*VARIETY	Var(Error) + 2 Var(FARM*VARIETY)
Error	Var(Error)
Variance ComponentVar(FARM)	Estimate
Var(FARM*VARIETY)	0.01632036
Var(FARM*VARIETY)	0.00083776
Var(Error)	0.00219208

Interpretation of the output

  There are slight differences in the statistics calculated, but the main distinction is the difference in the expected mean square column, SAS assuming an unrestricted model.  The VARCOMP procedure does not calculate the F-ratios, though it does calculate the rest of the analysis of variance table.  The ANOVA commands in SAS allows random factors and this command calculates the F-ratios.


INTERPRETATION OF THE OUTPUT

  The difference between the interpretation of the Minitab output (with the restricted model selected) and SAS output (with the unrestricted model assumed) lies in the interpretation of the farms effect.  From the SAS expected mean square column it can be seen that MSfarms should be compared with MSinteraction.  The F-ratio is 0.13443050/0.00386761 = 34.76, with 2 and 6 df.  Comparing this to F tables it can be seen that this ratio is highly significant (P , 0.001) so the conclusion is the same as that reached from the restricted model.

  In this example the interaction is not significant  and no conflict arises.  Such a conflict could arise if the interaction is significant and therefore MSinteraction > MSfarms.

  Assuming a restricted model, the random effect is compared with the error.  In the unrestricted model it is compared with the interaction .  Thus in the restricted model the overall random effect may be judged significant whereas in the unrestricted model it may be judged non-significant.

  However, even if overall it is deemed not to be a significant effect the presence of the significant interaction indicates that the effect of the random effect differs according to the level of the fixed effect, so some differences in the random effect are indicated.


Factorial 2 crossed with replicates pg 3
COMPUTING THE FACTORIAL ANOVA TABLE	page 3 of 4

Two Factors Crossed, with replicates

  The loblolly pine data are also analysed by a two-factor crossed ANOVA, but here both factors (date of sowing and method of ground clearing) are fixed effects.  This second example illustrates the interpretation of a fixed effects model, and the identification of interaction effects.

The Minitab command and output is:

ANOVA C2 = C3 | C4	[C2 - 'seedlings' C3 and C4 'date' and 'burning' classifications]

Factor	Type	Levels	Values
date	fixed	6	1	2	3	4	5	6
burning	fixed	2	1	2
Analysis of Variance for Seedling

Source	DF	SS	MS	F	P
date	5	7500086	1500017	68.88	0.000
burning	1	369252	369252	16.96	0.000
date*burning	5	686385	137277	6.30	0.000
Error	36	783925	21776
Total	47	9339648


 The F-ratios for both main effects, date and burning, and for interaction are all significant with P values less than 0.001.

  The conclusions are that there are overall differences in the mean number of seedlings germinated at different dates, and overall differences depending on whether the ground was cleared by burning or not.

  The significant interaction confuses this interpretation somewhat and is most easily sorted out by a graph.


Factorial 2 crossed with replicates pg 4
COMPUTING THE FACTORIAL ANOVA TABLE	page 4 of 4

Two Factors Crossed, with replicates

  The significant interaction confuses the interpretation somewhat and is most easily sorted out by a graph.  Cell means are plotted against date of sowing, with the different methods of ground clearing identified.

         -

         -                                  A

     1800+

         -

         -                                            A

         -                        A         B         B

         -

     1200+                        B

         -              B

         -    B         A

         -    A

         -

      600+                                                      A

         -

         -                                                      B

         -

         -

           --------+---------+---------+---------+---------+--------

                 4.0       8.0      12.0      16.0      20.0

        N* = 86

        A = burnt vs. date        B = cleared vs. date

Change Explan!!!!!

  It is clear from this figure that seedlings do better on burnt ground if sown early in March, but otherwise do worse.  Also seedlings on conventionally-cleared ground do much better than on burnt ground if sown in the middle of the month;  the differences are not so great at other times.  This graph as well as presenting the data in an easily interpreted form also suggests the possibility of using a parabola to predict the number of seedlings (see Multiple Regression).


Factorial 2 hier.nested pg 1
COMPUTING THE FACTORIAL ANOVA TABLE	page 1 of 3

Two Factors with One Hierarchically Nested within the Other

  It is vital to have replicates with hierarchically nested designs since, in the absence of replicate data, it is not possible to differentiate between the effect of the nested factor and the error.  Again before the ANOVA can be computed it must be decided which factors are random and which are fixed.

  The two-factor hierarchical analysis of variance is illustrated by using just the data relating to postgraduate students from the heterotroph data.

Water type	Student	Postgraduates	Undergraduates
Sewage	1 2	2700   2800   1700 2600   3000   3200	4900   3600   1300 3500   4100   2500
Polluted River	1 2	52      49       61 68      75       83	60     140       75 40       60       80
Clean Stream	1 2	5.9     7.6    16.0 5.6     5.9      6.3	6.4      4.2      4.1 8.3      3.1      5.2
  Here the students are a hierachical factor within water types because there are six students, two for each water type, and the students numbered do not have anything in common with each other, no distinctive feature.


Factorial 2 hier.nested pg 2
COMPUTING THE FACTORIAL ANOVA TABLE	page 2 of 3

Two Factors with One Hierarchically Nested within the Other

  It was found that the raw data did not justify the assumptions of the ANOVA and a transformation was necessary.  Thus the following results are from loge transformed data.  (The assumptions and transformation steps follow this section).

The Minitab commands and output for the Heterotroph data are:

ANOVA C2 = C3 C4(C3);	[C2 -loge'count' C3-'water type' C4 - 'student']

RANDOM C4;	[student number is random]

EMS.	[gives expected mean square]

Factor	Type	Levels	Values
Waters	fixed	3	1	2	3
Students(Waters)	random	2	1	2
Analysis of Variance for LogN

Source	DF	SS	MS	F	P
Waters	2	106.29	53.135	321.42	0.000
Students(Waters)	3	0.496	0.165	2.58	0.102
Error	12	0.768	0.064
Total	17	107.533
Source	Variance component	Error term	Expected Mean Square (using unrestricted model)
1 Waters		2	(3) + 3(2) + Q[1]
2 Students(Waters)	0.03378	3	(3) + 3(2)
3 Error	0.06398		(4)
Note C4(C3) in the command line specifies student number nested within water type.  The same convention is used in SAS.


Factorial 2 hier nested pg 3
COMPUTING THE FACTORIAL ANOVA TABLE	page 3 of 3

Two Factors with One Hierarchically Nested within the Other

  The F-ratio for students is 2.58 with 3 and 12 df, and this is not significant (P = 0.102) so no significant differences are identified between students.  The F-ratio for waters is 321.42, with 2 and 12 df.  This is highly significant (P>0.001).  So there are big differences between the mean numbers of heterotrophs from the three water types.  (This probably should not surprise a microbiologist.)

  In view of the non-significant students' effects, the extraction of a variance component for students is unlikely to be of any practical importance, but for the purposes of illustration the calculation is:

The variance for students is estimated as 0.03378 and that for error as 0.06398

Thus variation between students' accounts for

0.03378/(0.03378 + 0.06398) = 0.346 of the total variance


Factorial balanced > factors pg 1
COMPUTING THE FACTORIAL ANOVA TABLE	page 1 of 2

Balanced Experiments with more than Two Factors

  Again it must be decided whether the factors are crossed or hierarchical.  For the full heterotroph data, there are three factors, Waters, PG-UG (post or undergraduate) and Students.  Waters and PG-UG are crossed fixed factors.  Students is hierarchically nested within waters and PG-UG, and could be regarded as fixed or random, depending on the desired interpretation.  Assuming it is random the Minitab commands and output are

ANOVA  C2 = C3|C5 C4(C3 C5);	[C2 - loge'Count',  C3 - 'Waters'

RANDOM C4;	 C4 'Students'  C5 'PG-UG']

EMS.

Factor	Type	Levels	Values
Waters	fixed	3	1	2	3
PG-UG Students(Waters PG-UG)	fixed random	2 2	1 1	2 2
Analysis of Varience for log N

Source	DF	SS	MS	F	P
Waters	2	230.897	115.449	898.56	0.000
PG-UG Waters*PG-UG	1 2	0.06 0.535	0.016 0.268	0.13 2.08	0.734 0.206
Students(Waters PG-UG) Error	6 24	0.771 3.09	0.128 0.129	1.00	0.451
Total	35	235.318
Source	Variance component	Error term	Expected Mean Square (using unrestricted model)
1 Waters		4	(5) + 3(4) + Q[1,3]
2 PG-UG 3 Waters*PG-UG		4 4	(5) + 3(4) + Q[2,3] (5) + 3(4) + Q[3]
4 Students(Water PG-UG) 5 Error	-0.00021 0.12912	5	(5) + 3(4) (5)
  Note the logged data is used so that the assumptions are justified (see next step - Assumptions of the Factorial ANOVA).


Factorial balanced > 2 Factors pg 2
COMPUTING THE FACTORIAL ANOVA TABLE	page 2 of 2

Balanced Experiments with more than Two Factors

  The conclusion is straightforward, only waters effect has a significant F-ratio, and that is hugely significant (P < 0.001).  Notice the negative variance estimate for students within waters and PG-UG in the expected mean square table.  negative variances are impossible, and all that is indicated here is that, since there is no significant students' effect, the students and error mean squares are estimates of the same thing, the error variance.

  In the analysis just described, all 36 data have been treated as replicate plots, tacitly assuming that 12 equal samples of water were randomly selected from each water source and randomly divided up amongst postgraduates and undergraduates for counts to be estimated.

  If the experiment were conducted differently and two sample bottles of each water type used, one each for postgraduates and undergraduate and the replicate counts were merely subsample from each bottle, then the interpretation would have been different.

  This would then be a split-plot experiment


Factorial split-plot experiments
SPLIT PLOT EXPERIMENTS

  If two sample bottles of each water type used, one each for postgraduates and undergraduates and the replicates were sub-samples from each bottle there are then two different kinds of plot.

  The first, called the main plot, is of sample bottles.  The factors waters and PC-UG apply to these main plots.  The second, called the sub-plot, is the sub-sample taken from each bottle and plotted out.  The factor students applies to these sub-plots.

Split-plot experiments require a different model to that used for factorial ANOVA.  They are often performed unwittingly, and incorrect interpretation may follow.

  If the heterotroph counting experiment had been done this way then there would be no replicate for the waters x PG-UG combinations, and the only way of proceeding with an analysis would use the waters x PG-UG interaction as the main plot error.

For further details see Winer


Winer ref


Factorial unbalanced data
UNBALANCED EXPERIMENTS

  The calculations needed for unbalanced data are lengthy, and computer packages have to be used, but there are also difficulties of interpretation.  Several ways have been suggested for interpreting unbalanced experiments.  Four such ways, called Types I to IV, and these are calculated in the SAS GLM procedure.  Minitab (V.7) and SPSS/PC+ also allow the calculation of ANOVA tables for unbalanced data.

I don't know the commands!!

  It is particularly difficult to analyse data sets for factorial experiments where some cells have no data.  In general, care has to be exercised in the use of packages since it is easy to misunderstand the manuals and wrongly interpret the output.

Missing Data

Comparisons of means for unbalanced data


Type I to IV Interpretation of unbalanced data

Milliken and Johson ref


Factorial missing data
 MISSING DATA

  Where observations are missing in an otherwise balanced experiment an approximate procedure is to fill in the missing observations with values derived from the remaining data and analysing the experiment as if it were balanced but reducing the total and error df by the number of missing observations.  In an experiment with replicates the cell mean is often used as a surrogate for the missing observation.

Formula when there are no replicates

  Since the error df is calculated automatically, assuming no missing data, in Minitab, SPSS/PC+ or SAS, MSerror has to be calculated by hand.  Thus all F-ratios must also be recalculated by hand.

Method of unweighted means


Factorial formula for no replicates (missing data) pg 1
MISSING DATA	page 1 of 2

Formula for when there are no replicates

  In the two-factor crossed experiment with no replicates a single missing, observation, i in the (i,j)th cell, is replaced by the value calculated from:

Where Xi is the total of the b-1 observation actually obtained at the ith level of Factor UA , Xj is the total of the a-1 observations obtained at the jth level of factor B, G  is the overall total of all the data obtained and a and b are the number of levels of factors D and B respectively.


Factorial missing data formula for no replicates pg 2
MISSING DATA FORMULA EXAMPLE	page 2 of 2

Using the depleted wheat yield data (bold data only)

	Variety 1	Variety 2	Variety 3	Variety 4
Farm 1 Farm 2 Farm 3	0.327      0.280 0.532      0.526 0.269      0.277	0 500      0.510 0.599      0.637 0.308      0.268	0.442      0.360 0.516      0.499 0.241      0.328	0.471      0.460 0.638      0.455 0.305      0.301

  The total X1 of the three observations available from farm 1 is 1.24 and the total X2 of the two observations for variety 2 is 0.907.  The overall total of the 11 data is 4.648.  The formula for the missing value then gives (3 * 1.24 + 4 * 0.907 - 4.648)/2 * 3 = 0.45.  This value replaces the missing observation and the data analysed as if it were a genuine observation, but the error df are now 5 and the total df are 10.  Since the error df is calculated automatically assuming no missing data, in Minitab, SPSS/PC or SAS, MSerror has to be calculated by hand.  Thus all F-ratios must also be recalculated.


Factorial missing data method pg 1
MISSING DATA

Method of Unweighted Means

  If the discrepancies in the numbers of data in the cells of a factorial experiment are not large, an approximate procedure known as method of unweighted means can be used.  For a two-factor crossed ANOVA the harmonic mean ( [image]) of the date is calculated as follows:

  This may look complicated but can easily be broken down into a few simpler steps, using the following example


Factorial missing data method pg 2
METHOD OF UNWEIGHTED MEANS EXAMPLE

Factor B (j)

1	1  *   *    *   *	2  *  *    *	3  *   *  *   *	4  * * *    *  *	5  * * *   *  *	<<  5 = n1,5(nij) (n = number of data points)
Factor   2    A (i)	*   *    *	*   *  *   *	*   *  *   *	*   *  *   *	* * *   * *
3	*   *     *	* * *   * *	*   *    *	*    *     *	*   *  *   *
1.	In the cell (1,5) the number of data points (n) is 5 (nij = n1,5 = 5)

2.	The reciprocal of these data points in each cell is taken and are summed:

3.	This figure is then divided by ab.

4.	The reciprocal of this number is then taken, and this is the harmonic mean [image].


Factorial missing data method pg 3
METHOD OF UNWEIGHTED MEANS	page 3 of 3

Minitab commands and output

 This method can be illustrated by an analysis of the depleted data set obtained from the lobolly pine seeds

(a)	Table of cell means

(b)	Two-way crossed ANOVA for the cell means

(c)	One-way ANOVA table

The harmonic mean of the nij 's of this data set is 2.76923 and this is used to calculate:-

(d)	Complete ANOVA table (calculated by hand using the harmonic mean and the two-way table).

Source	DF	SS	MS	F
date	5	5020240	1004048	36.87
burning	1	271667	27667	9.98
date*burning	5	372943	74589	2.74
Error	23	626267	27229
Total	34	6525417


The last two lines for error and total are copied directly from the one-way ANOVA with cells as groups.


The df for date, burning and the date x burning interaction are the same as those of the date burning abd error resoectively from the two-way ANOVA on cell means.


The sums of squares for these sources of variation are each multiplied by the harmonic mean, in this case 2.76923.  Thus for dates , the sums of squares is 2.76923 x 1812865 = 5020240.


The MS's are calculated by hand by dividing the SS by the df


Finally the F-ratios are calculated by hand, in this case, with both effects assumed fixed, the denominator for all three ratios is MSerror


Rows: date    Columns: burning

1 2 3 4 5 6 ALL	1 870.0 1115.0 1410.0 1855.0 1646.7 645.0 1194.1	2 1006.7 1156.7 1223.3 1390.0 1400.0 280.0 1076.7	ALL 938.3 1140.0 1316.7 1622.5 1505.7 488.6 1133.7
CELL CONTENTS -- seedling : Mean


Analysis of Variance for meanseed

Source	DF	SS	MS	F	p
date	5	1812865	362573	13.46	0.006
burning	1	98102	98102	3.64	0.115
Error	5	134674	26935
Total	11	2045641

Two-way crossed ANOVA for the cell means


ANALYSIS OF VARIANCE  ON seedling

SOURCE	DF	SS	MS	F	p
Cells	11	5899150	536286	19.7	0.000
ERROR	23	626267	27229
TOTAL	34	6525417

One-way ANOVA for the depleted data with cells as groups


Loblolly Pine Seed Data

	Number of seedlings on
Sowing date  in march	Burnt ground		Conventionally  cleared ground
2 6 10 14 18 22	900' 760 880' 1060 1530 1390' 1970' 1820 1960 1310 830 570	810 1040 1170 910' 1160 1540 1890 2140' 1670 1480' 420 760	880 960' 1050' 1110 1140 1320 1360' 1490' 1270 1500 150 420	1100 1040 1240 1120 1270 1080 1510 1270 1380 1450 380' 270	observations marked with a ( ' )  are assumed to be missing  in some examples


Factorial Comparison of means for unbalanced data
COMPARISON OF MEANS FOR UNBALANCED EXPERIMENTS

Minitab commands and output

Using the depleted loblolly pine data

GLM C2 = C3|C4	[C2 - data, C3 - data C4 - burning]

Factor	Levels	Values
date	6	1	2        3         4         5         6
burning	2	1
Analysis of Variance for seedling

Source	DF	Seq SS	Adj SS	Adj MS	F	p
date	5	5267869	5515767	1103153	40.51	0.000
burning	1	274345	271667	271667	9.98	0.004
date*burning	5	356938	356938	71388	2.62	0.051
Error	23	626267	626267	27229
Total	34	6525417
Unusual Observations are also printed

  There is little difficulty in interpreting these data, though the interaction is only marginally significant.  An analysis of the data with no missing values gives a highly significant interaction.  Thus a loss of only 13 observations reduces the sensitivity of the statistical test.

SAS GLM PROCEDURE


Factorial SAS for Comaprison of means.  Unbalanced
COMPARISON OF MEANS FOR UNBALANCED EXPERIMENTS

SAS Commands and Output

PROC GLM;

CLASS DATE BURNING;

MODEL SEEDLING = DATE | BURNING/SS1 SS2 SS3 SS4;

RUN;

General Linear Models Procedure

Dependent Variable:SEEDLING

Source	DF	SS	MS	F-Value	Pr > F
Model	11	5899150.476	536286.407	19.7	0.0001
Error	23	626266.667	27228.986
Corrected Total	34	6525417.143
	R-Square	C.V	Root MSE	Seedling Mean
	0.904027	14.55500	165.0121	1133.71429
Source	DF	Type I SS	MS	F Value	Pr > F
DATE	5	5267868.333	1053573.667	38.69	0.0001
BURNING	1	274344.607	274344.607	10.08	0.0042
DATE*BURNING	5	356937.535	71387.507	2.62	0.0512
Source	DF	Type II SS	MS	F Value	Pr > F
DATE	5	5421607.563	1084321.513	39.82	0.0001
BURNING	1	274344.607	274344.607	10.08	0.0042
DATE*BURNING	5	356937.535	71387.507	2.62	0.0512
Source	DF	Type III SS	MS	F Value	Pr > F
DATE	5	5515766.343	1103153.269	40.51	0.0001
BURNING	1	271667.308	271667.308	9.98	0.0044
DATE*BURNING	5	356937.535	71387.507	2.62	0.0512
Source	DF	Type IV SS	MS	F Value	Pr > F
DATE	5	5515766.343	1103153.269	40.51	0.0001
BURNING	1	271667.308	271667.308	9.98	0.0044
DATE*BURNING	5	356937.535	71387.507	2.62	0.0512


In the sequential SS of the Minitab output, the terms are included in the model in the order specified by the command line.  The Type I SS given by SAS is the same as this sequential SS.


The adjusted SS of Mintab gives the SS for each term when that term is included last in the model.  It is these sums of squares that Minitab uses for the calculation of MS, F and P columns.

  The Type III SS given by SAS is the same as this adjusted SS and if the tests of the hypothesis that the effects are zero are of interest, then it is these Type III (or adjusted) SS that should be interpreted.

  However the Type III sums of squares do not, in an unbalanced experiment, add up to the total explained by the model, and if the building of a model is of interest, then the Type I (or Type II)SS should be used


The Type II SS are the extra explained when that term is included with all other terms of the same order.  Thus the SS for dates when burning (but not the interaction) is included is 5421608.  For data sets with no empty cells the Type III and Type IV SS are the same.


Factorial assumptions of fact. ANOVA
3. ASSUMPTIONS OF A FACTORIAL ANOVA

  Underpinning the analysis of the factorial experiment is a linear mathematical model.  This model contains an error term which is assumed to be independent and normally distributed. (The errors are thus assumed all to have the same variance).

  The diagnostic checks to justify these assumptions are similar to those used in the one-way ANOVA.  Diagnostic checking is an important part of the ANOVA procedure, genuine differences may not be detected if the assumptions of the model do not hold.

Homogeneity of Variance

Normality of Residuals


Linear mathematical model

xij = [image] + [image] + [image] ij


Factorial homogeneity of variance
HOMOGENEITY OF VARIANCE

  The main distinction between the multi-factorial and one-way cases is in the tests for homogeneity of variance.  These tests are done on variances calculated from the cell means, and these are not normally calculated as part of a multi-factorial ANOVA table, though the within group variances (or standard deviations) are often presented in a standard one-way analysis.

The Minitab command to calculate cell variances is:

TABLE

Alternatively a one-way ANOVA with the cells as groups can be performed.

The three main tests are:

Bartlett's [image] (best carried out on a computer)

Hartley's Fmax and Cochran's C


Factorial homogeneity of variance example
HOMOGENEITY OF VARIANCE EXAMPLE

Using the heterotroph count data:

         -

         -

         -

     1800+                                                *

         -

 SD      -

         -

         -

     1200+

         -

         -

         -                                                 *

         -

      600+                                   *

         -

         -                                           *

         -

         -

        0+ 44

           +---------+---------+---------+---------+---------+------Mean

           0       700      1400      2100      2800      3500

        N* = 19

  The cell means of each set of three counts made by a student are plotted against the cell standard deviations (using the Minitab Plot command).  There is a strong tendency for larger standard deviations chance variances) to be associated with higher estimated counts of bacteria, thus the raw data is heterogeneous.  This means that the data may need to be transformed.


factorial hartleys fmax cochrans c
HOMOGENEITY OF VARIANCE

Hartley's Fmax and Cochran's C

These tests are calculated by

where [image] is the largest group variance, [image] is the smallest group variance and [image] are all the group variances.

These values are compared with tabulated values at P = 0.05.

  In all three of these tests calculated values smaller than the table values are required for the variances to be homogeneous and this assumption of ANOVA to be upheld.

  The Fmax and C tables are designed for equal group sizes.  If the data has unequal n then the tables can still be used to give an approximate result.

  Care must be taken when interpreting the results of computer output with these tests.  The probability values often given (SPSS/PC+, Systat, Statgraphics) should be greater than P=0.05 for the variances to be significantly homogeneous.  Notice that for these tests the desirable outcome is for the P value to be large and not small as is the case for many other tests.  We can then be confident about proceeding with the one-way ANOVA assuming equal variances.


For example,some data has n = 5-15 with Fmax = 148.8 and C = 0.651.  The tabulated values for highest and lowest n are Fmax(4,15) = 4.01, Fmax(4,5) = 13.7, C(4,16) = 0.437 (the nearest value in the table, to C(4,15)) and C(4,5) = 0.590.  Thus the calculated values are larger than all possibilities from the tables and so it is clear that the variances are once again declared heterogenous.  If the calculated value for either of these test was between the two tabulated values the result would be uncertain, but if less than the smallest tabulated value (that for the largest n) then the variances would be homogeneous.


Factorial normality of residuals
NORMALITY OF RESIDUALS

The residuals from the analysis can be stored using the Minitab subcommand:

ANOVA     (etc);

RESIDUALS C7.     (empty)

The tests for normality of residuals are:

 Normal probability plots and Shapiro-Wilk test

 boxplots


Factorial normality of residuals example
NORMALITY OF RESIDUALS NORMAL PROBABILITY PLOTS

  Using the heterotroph count data and the NSCORE Minitab command the residuals from the ANOVA are plotted against the normal scores from NSCORE.

 (Linearity = normality of residuals)

   Nscore

         -                                           *

      1.5+                                       *

         -                                     **

         -                                  ** *

         -                                 3*

         -                                 4

      0.0+                                 4

         -                                 4

         -                                 4

         -                                 3

         -                            *    *

     -1.5+                       *

         -                     *

         -     *

         -

           --+---------+---------+---------+---------+---------+----resids

         -2100     -1400      -700         0       700      1400

  The curvature of the plot for the raw data is evident.  Using the correlation command in Minitab

(CORR), the correlation coefficient can be compared to table values.  A comparison of the coefficient 0.798 for the raw data against the critical table value (n = 24) 0.956 shows that since 0.798 < 0.956 the residuals are not normal - confirming the already obvious deviation from linearity of the plot. (This test is essentially thesame as the ShapiroWilkox test for normality).


factorial boxplots
NORMALITY OF RESIDUALS     BOXPLOTS

  Another straightforward graphical way to examine residuals is by using the Box-and-Whiskers plot or boxplot.  This is really a graphical representation of a letter value display.  The Minitab commands are

BOXPLOT  C1

or

GBOXPLOT  C1      (for high resolution graphics plot)

C1  contains the data to be plotted (raw or transformed)

The output will appear
  The boxplot shows the median (middle data point), and hinges (middle points between the median and maximum and minimum values), adjacent values (the end of the whiskers) and the more extreme values outside these values (outside and far outside values).  Thus residuals that are evenly distributed will show symmetrical boxplots and uneven residuals will give non-symmetrical boxplots.  As the normal distribution is even, normally distributed residuals will give an even boxplot with a central median and slightly longer whiskers than the median-hinge distance.  Boxplots are part of a useful, wider group of statistical procedures called exploratory data analysis.(see Transformations for more details)


Factorial transformations
4. TRANSFORMATIONS

  If analysis of the raw data shows that the variances are not homogeneous and the residuals are not normally distributed then none of the results of the factorial ANOVA can be relied upon.

  The usual action is to search for a suitable transformation which will modify the scale upon which the data is expressed so that the two critical assumptions are valid.  Biologists often think that transformation somehow corrupts the original data, but this is not so.  It is a perfectly valid procedure and will often make understanding the data easier as well as enabling ANOVA to be carried out.

  It is also possible that the critical assumptions might not be valid because of severe outliers within the data.  Hence this possibility should be considered whilst the data is being examined.

  After transformation the critical assumptions must be examined again to check that the transformation has worked.

HOW TO CHOOSE TRANSFORMATIONS


Factorial Choosing Transformations
HOW TO CHOOSE AN APPROPRIATE TRANSFORMATION

The approaches outlined in steps 2-5 below are given in order of ease of computation with most computer packages.

1.	Select an appropriate method according to the data structure and ease of computation
 any data; steps 2 and 5
 several groups with small numbers of replicates;  steps 2 and 3
 at least one large group;  step 4

2.	Choose a transformation based on theory.

3.	Taylor's power law.

4.	Letter value display method.

5.	Minimization of homogeneity of variance statistics

For more details see Transformation Protocol

Note

  There is no need to use all the methods suggested here to choose a transformation for unfamiliar data sets.  It is best to select one of these methods and to base the choice on the arrangement of the data and the ease with which the computations can be done with the computer package being used.


Factorial Transformation Protocol
TRANSFORMATION PROTOCOL

1.	Select an appropriate method according to the data structure and ease of computation
 any data; steps 2 and 5
 several groups with small numbers of replicates;  steps 2 and 3
 at least one large group;  step 4

2.	Choose a transformation based on theory, some examples follow
 population data;  log10 x or log10 (x+1)
 low counts in quadrats; [image] or
 percentages or proportions; the arcsine transformation,
    or

3.	Taylor's power law

(a) Calculate group means [image] and variances [image] for raw data.

(b) Plot log10 [image] against log10 [image] for raw data.

(c) Calculate slope (b) of the regression line from the above plot where

log10{bm ssq.bmp}  =  a + b log10

(d) Calculate appropriate power (p) for transformation where

(e) Use power directly to transform by [image] or choose nearest power transformation from the ladder of powers .

4.	Letter value display method

(a) Transform data according to ladder of powers .

(b) Compute midsummaries for all transformations.

(c) Calculate se from spreads/normal spreads for all transformations.

(d) The best transformation has the most even range of midsummaries and se values.

5.	Minimization of homogeneity of variance statistics.

(a) Transform all data according to ladder of powers

(b) Calculate Hartley's Fmax, Cochran's C and Bartlett's [image] statistics for all transformations.

(c) Use transformation giving the minimum value for most statistics.


Factorial Based on Thoery
CHOICES BASED ON THEORY OR EXPERIENCE

  Biologists who regularly analyse their data with correctly applied statistical procedures will soon find the best transformations for the different types of data they generate.  Thus the more time consuming methods of searching for a good transformation will rarely need to be routinely applied.

  The most commonly used are

Logarithmic Transformations

Square Root Transformations

The Arcsine Transformation


Factorial Log Transf
LOGARITHMIC TRANSFORMATIONS

  A log10 transformation is often suitable for population data.  If the data contains zero values then a small constant number should be added to all data points.  This constant should represent the smallest number that could sensibly be found, that is the limit of detection, for example log10(x+1).Many population related variables can be sensibly transformed by logarithms.  Enzyme activities of bacteria, primary production in seawater and grazing rates of zooplankton are all examples of this.  Data collected from survival experiments often naturally follow a logarithmic decline and hence can be log10 transformed.

  Using log10 the base of logarithms is easy to understand as the number of orders of magnitude; 1 million insects is [image] or 6.0 when transformed.  However, logarithms to any base can be used and some may prefer loge.  These two forms can be easily interconverted in statistical packages without one or the other transformations as log10x is equal to logex multiplied by loge10.


Factorial Sq Root Transf
SQUARE ROOT TRANSFORMATIONS

  Data collected in quadrats are distributed according to the Poisson distribution.  Haemocytometer counts of blood cells, counts of bacteria in petri dishes and plants in quadrats are examples of this type of data collection.  When the counts are low (<20) the distributions are very skewed and are best transformed by the square root of x or x+1.  When counts are high (>20) the Poisson and Normal distributions are so similar that transformation might not be necessary.


Factoeial Arcsine Transf
THE ARCSINE TRANSFORMATION

  This is a special transformation which is for proportions or percentages and can be applied almost automatically to data of this type.  The transformation for proportions.

and has been designed to make the variances homogeneous.  With percentages the values should be converted to proportions by division by 100.  It can be looked up in tables but is most conveniently calculated in degrees by computer.  The Minitab command for this conversion assumes the raw data are percentages in column 2 (C2) and the transformed data are to be stored in column 3 (C3).

LET C3 = (ASIN(SQRT(C2/100))) * 57.296

Other packages make use of similar commands.


Factorial Taylors Power Law
TAYLORS POWER LAW	page 1 of 2

  This provides a very convenient method for finding the best transformation for a set of data.  It works best when there are several groups (> 6) and is effective even when the group size is small.

  The method (Transformation Protocol , step 3) works by calculating the slope (b) of the regression line of the log10 of the group variances upon log10 group means.  The power (p) for the transformation [image] is then calculated by p = 1 - b/2.

  For the resulting power transformation to be effective the plot of log10 [image] versus log10 [image] should be linear.  That is all the points should lie close to a straight line.


Factorial Taylors Power Law pg2
TAYLORS POWER LAW	page 2 of 2

  The power transformations indicated by the Taylor's power law method will be slightly different for each new set of data collected and so do not provide continuity of statistical analysis within a research programme.  It would be better to use the calculated power as an indication of the best transformation to use.  The ladder of powers provides a logical way of doing this.

	Transformation		Addition of small  number needed for zero values
Power	Name	Formula
3 2 1 0.5 0 -0.5 -1 -2	cube square raw data square root logaritm reciprocal root reciprocal reciprocal square	x log10x 1/ 1/x 1/	no no no* yes yes yes yes

Based on a table in Velleman and Hoaglin

* Not needed but [image] or [image] is often used.

  The ladder puts the raw data and a range of equivalent transformations into a logical sequence from which a suitable transformation can be selected from the power suggested from Taylor's power law.


Reference

Velleman, P.F. and Hoaglin, D.C. (1981). Applications, Basics, and Computing of Exploratory Data Analysis. Duxberry Press, Boston


Factorial Letter Value Displays pg1
LETTER VALUE DISPLAYS	page 1 of 2

This is another technique of exploratory data analysis and Minitab will allow calculation of the associated statistics.  This is a diagrammatic representation of how the letter values are assigned and the midsummaries and spreads calculated.

  Briefly the calculations are done as follows.

  When the data is ordered the median is the middle data point and is assigned the letter value M.  The two middle data points between the median and the lowest and highest values are called hinges and assigned H.  This principle continues through eighths (E) and letter values D,C,B,A,Z etc. which do not have names.  When there is not a single central value the central two data points are averaged to calculate the number required.  The midsummaries are then calculated from averages of the pairs of letter values; hence, for example, mid H = (upper H + lower H)/2.

  For a perfectly symmetrical distribution these midsummaries will all be equal but for a distribution skewed to the left they will increase and will decrease for one skewed right.


Factorial Letter Value Displays pg2
LETTER VALUE DISPLAYS	page 2 of 2

  Spreads are calculated as the differences between the letter values and so are the magnitudes of ranges between the letter values.  Hence, for a normal distribution they will be entirely predictable.  The values of the spreads from a standard normal distribution (mean=0, standard deviation=1) are called
normal spreads Thus the relationship between a normal distribution and the data being analysed can be examined from the ratios of the calculated spreads and the normal spreads.  This ratio is an estimate of the standard deviation of the data (se) so

and will be constant if the distribution of the data is close to normal.

  As with the midsummaries the normal spreads will increase or decrease for skewed distributions.

  These ideas have been used to make a test to select the best transformation for a set of data.

  Click for details


Factorial Ladder Test
??-same for factorial??

TEST TO SELECT THE BEST TRANSFORMATION

  This test should be done on data from a single group in the factorial ANOVA and only works well for data sets with at least one large group (>ca.12).  The test is done by transforming the data from the largest group by all the transformations in the ladder of powers, taking care to add an appropriate constant (perhaps 1 or 0.5) to all the values if zeros are present.  The midsummaries and se values are then calculated and the best transformation is the one in which these values remain constant and do not increase or decrease in a regular manner.

  A Minitab macro (LADDER.MTB) can be used which calculates the midsummaries and se values.  Although this approach to selecting a transformation might not appear very precise it works remarkably well.


Spreads for letter values from a standard normal distribution

Letter value	Normal Spread
H E D C B A Z	1.349 2.301 3.068 3.726 4.308 4.836 5.320


Oneway MINIMIZATION
MINIMIZATION OF HOMOGENEITY OF VARIANCE STATISTICS

  The three statistics used to check data for homogeneity of variance (Bartlett's [image], Fmax, C) can also in principle be used to select a good transformation.  This can be done by calculating these statistics for data transformed according to the ladder of powers.  By comparing the critical values from tables at
P = 0.05 with the calculated values for the different transformations a suitable transformation can be chosen.


Factorial Minimize Table
RESULTS FROM TESTS FOR HOMOGENEITY OF VARIANCE

  Here are some example values which could have been obtained from the various tests to illustrate the method.

Transformation	Bartlett's	Hartley's Fmax	Cochran's C
Square Raw Data Square Root Log10 Reciprocal Root Reciprocal TABLE VALUE	127.9 45.9 19.09 9.80 13.99 23.5 11.07	> 7374 218.4 29.0 263.6 7075 62.0	0.964 0.697 0.587 0.620 0.482 0.309 0.532

  Hence for these data Bartlett's [image] and Fmax both reach acceptable minima at the logarithmic transformation, but C is below the critical value for both types of reciprocal transformation.

  Additionally Bartlett's [image] indicates that log10x and reciprocal are acceptable and C suggests reciprocal is satisfactory.  Clearly with different data sets these tests are likely to have different discriminatory powers.


Factorial outliers
5. TESTING FOR OUTLIERS

  Sometimes data sets might contain one or more very extreme values that may not really be part of the same data.  Such values are outliers and can occur for a variety of reasons, such as transcription or experimental errors.

  Outliers can sometimes be seen when residuals are being tested for normality.  For example probability plots of the residuals show values with high/ low residuals that might be outliers.

  Inspection of a table of residuals shows these suspect points .  These points and can also be detected easily by plotting the residuals against the order of the data.  Minitab has a time series plot for this purpose.

The most useful test for the significance of outliers is DIXON'S TEST


Factorial Mtb Outliers
MINITAB COMMANDS AND OUTPUT

The Mintab command is  TSPLOT  C10   < contains the residuals from the transformed data

An  example plot output will appear

         -

     0.60+     5

         -        8

 resids  - 1                    2

         -    4     0  3  6   0

         -  2

     0.00+      6    12 4      1 3

         -                 789

         -         9

         -               5        4

         -

    -0.60+   3

         -

         -

         -

         -

    -1.20+       7

          +---------+---------+---------+

          0        10        20        30

This confirms that the point seven in this data set is suspect.


Factorial Dixons Test pg1
TESTING FOR OUTLIERS	page 1 of 2

Dixon's Test

  The most useful test for the significance of outliers is Dixon's test which is designed for small group sizes (n<25).  However, Grubb's test is available for data sets with large group sizes (n>25).  To test the possible outliers the raw data must comply to the assumptions or transformed values must be used because Dixon's test only works with normally distributed data.The correct formula for a group size of n = 4 is

(The correct Dixon's test formula for a particular group size must be used).  The data is then ordered and values corresponding to x1, x2 and xn determined.  The outlier of interest is assigned x1 whether it is the smallest or largest value and the rest of the order maintained accordingly.


Factorial Dixons Test pg2
TESTING FOR OUTLIERS	page 2of 2

Dixon's Test

Using an illustaive example

	x	x + 1	log10(x + 1)
	0	1	0.00	x1
	13	14	1.15	x2
	42	43	1.63	x3
	61	62	1.79	x4
  Consulting a table of critical values at n = 4 shows that r10 is greater than the table values at all three probabilities quoted, so the outlier is not significantly (P > 0.1) lower than the other values in the group.  Thus the value should not be considered for exclusion from the analysis.  If an outlier is significant it should not automatically be excluded from the data set.  Firstly, field and/or laboratory notebooks should be consulted to see whether the data point is correct or whether there is a good scientific reason for exclusion.


Factorial comparisons
6.  COMPARISON OF MEANS OF MAIN EFFECTS

  A posteriori or a priori tests on means of main effects are done in multi-factorial experiments in the same way as for one-way ANOVA.

A priori tests

A posteriori tests

SAS commands and output


Factorial SAS comparisons commands
COMPARISON OF MEANS

SAS Commands and Output

  Sums of squares due to contrasts are available within the GLM procedure of SAS, and are specified by the contrast statement.  Using the wheat data to compare farms 1 against the mean of farm 2 and farm 3 and to compare farm 2 and farm 3, the commands and output are:

PROC GLM;

CLASS FARM VARIETY;

MODEL YIELD=FARM | VARIETY;

CONTRAST 'Farm 1 v Others' Farm 1 -.5 - .5;

CONTRAST 'Farm 2 v Farm 3' Farm 0 1 -1;

RUN;

Source	DF	Type I SS	Mean Square	F value	Pr > F
FARM	2	0.26886100	0.13443050	61.33	0.0001
VARIETY	3	0.03849433	0.01283144	5.85	0.0106
FARM*VARIETY	6	0.02320567	0.00386761	1.76	0.1895
Source	DF	Type III SS	Mean Square	F value	Pr > F
FARM	2	0.26886100	0.13443050	61.33	0.0001
VARIETY	3	0.03849433	0.01283144	5.85	0.0106
FARM*VARIETY	6	0.02320567	0.00386761	1.76	0.1895
Contrast	DF	Contrast SS	Mean Square	F value	Pr > F
Farm 1 v Others	1	0.00001875	0.00001875	0.01	0.9278
Farm 2 v Farm 3	1	0.26884225	0.26884225	122.64	0.0001
  (Farms are assumed fixed)

It is clear that the contest comparing farms 2 and 3 is the major component of the difference between farms, caused by the low yields obtained on farm 3.


Factorial A PRIORI TESTS1
A PRIORI TESTS	page 1 of 2

  These tests are perfomed in the same way as for one-way ANOVA comparisons, thus to illustrate the techniques used a simple one-way data set will be used-E.coli data

  The a priori tests sometimes called contrasts also include tests for similarity between means as opposed to differences between means.  For pairwise comparisons a version of the tests that work for both equal and unequal sample sizes can be used.

  These tests are very patchily represented in computer packages.  Minitab does not have any commands for these tests in release 7.1.  SPSS/PC+, Systat and Statgraphics have several tests of this type.  All three of these packages currently allow computation of the following pairwise comparisons and SPSS/PC+ and Systat allow a priori contrasts to be performed as well.


Factorial A PRIORI TESTS2
A PRIORI TESTS 	page 2 of 2

  These tests can only be done for groups which could sensibly have been tested before the experiment was done.  For the E. coli data set the researcher could sensibly have planned to compare the sites immediately above and below the tributary site 3 and 4; (see site details) and so these are candidates for an a priori pairwise comparison.  Similarly the rocky site above tributary station (site 1) could be compared a priori with the other two less rocky sites above the tributary (sites 2 and 3).  These tests are more powerful at detecting differences than the a posteriori tests and so stringent rules must be followed.  These rules are

tests can be applied regardless of significance of basic ANOVA

the maximum number of comparisons allowed is a-1 (where a = total number of groups in the ANOVA).

  Because these tests are so powerful they should be used in preference to a posteriori tests wherever validly possible.

Least significant difference (LSD)

Orthogonal decomposition


Factorial LEAST SIGNIFICANT DIFFERENCE
LEAST SIGNIFICANT DIFFERENCE (LSD)

  This pairwise comparison is valid only as an a priori test, despite being linked with a posteriori tests in the SPSS/PC+ command structure.  The LSD is the smallest difference between two means that is significant, one value is calculated for all the comparisons to be made.  Calculation of LSD is straightforward.

where the two-tailed t-statistic is quoted at a probability P and for within groups degrees of freedom (dferror) from the one-way ANOVA table; ni and nj are the number of replicates in the largest and smallest groups.  So for the E. coli data at P = 0.05

where 0.176 is the MSerror from the log10(x+1) transformed data.  As the difference between the group means at the sites above and below the tributary is 1.096 there is a significant difference between them.


Factorial ORTHOGONAL DECOMPOSITION1
ORTHOGONAL DECOMPOSITON	page 1 of 4

  These tests are to allow comparisons among groups of group means.  The method essentially performs a new one-way ANOVA on a restructured set of data.  A test sum of squares (SStest) is calculated according to the following type of formula

where g1 to g4 are the sums of the values in the four groups used in this test and n1 to n4 are all the numbers of replicates in the corresponding groups.  This formula can be extended (indicated by dots) or contracted to test more or less group means in each of the terms, or by adding more terms.


Factorial ORTHOGONAL DECOMPOSITION2
ORTHOGONAL DECOMPOSITION	page 2 of 4

  The next step is to calculate a test mean square (MStest) and F (Ftest)

where dftest is one less than the number of terms used to calculate SStest.  Ftest is then compared with a table value (FP(dftest, dferror)) and if Ftest > Ftable the result is significant.

  Thus to compare the E. coli data at site 1 with those at sites 2 and 3 (groupsums,g, ) the calculations arebmc sstest2.bmp}

From F-tables,  FP(1,18) is	4.41    for  P = 0.05

8.28    for  P = 0.01

15.4    for  P = 0.001

Hence, there is a significant difference in E. coli numbers due to rockiness above the tributary at the
 P=0.01 but not the P=0.001 level.


Factorial Orthogonal Decomposition pg3
ORTHOGONAL DECOMPOSITION 	page 3 of 4

  Particularly with larger data sets these calculations can be tedious and can easily be done by computer by re-labelling the groups appropriately.

	log10	sites	rocky	below
Using the log10(x+1) transformed E. coli data, the labels in the column headed 'rocky' compare means at sites 2and3.   Those in the column headed 'below' compare the means at sites 4, 5 and 6.   Misssing values are labelled with an asterisk	0.903 0.699 0.000 0.778 1.792 1.146 0.000 1.633 1.707 2.193 2.029 2.003 1.079 0.845 0.477 1.146 0.000 0.000 0.000 0.301 0.301 0.699 0.301 0.000	1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6	1 1 1 1 2 2 2 2 2 2 2 2 * * * * * * * * * * * *	* * * * * * * * * * * * 1 1 1 1 2 2 2 2 3 3 3 3	Click for details of Minitab commands


Factorial Ortho Mtb Commands
ORTHOGONAL DECOMPOSITION	page 4 of 4

Minitab Commands

The command is the same as for the original analysis of these daa

ONEWAY  C1  C2  C3  C4

C1  	[contains the transformed data]

C2  	[contains the labels(in this case the 'rocky' or 'below'column labels)]

[C3  and  C4 (empty)  are where Minitab stores the residuals and group means.]

  Performing a one-way ANOVA using either of these new sets of subscripts will show the SStest in the normal SSgroups position.  The Ftest however, must be calculated with the MSerror from the entire original one-way ANOVA table.

  By this method the SStest below the tributaries is 1.382 and so MStest=1.382/2=0.691 and hence
Ftest=0.691/.176=3.926.  In this test dftest is 2 because there were 3 terms compared and can be found as dfgroups in th modified ANOVA table used to calculate SStest .  As F0.05(2,18) is 3.55 there is just a significant difference between thes means

  Four means comparisons have now been done by a priori tests (above and below tributary=1 ; 'rocky'=1; 'below'=2; total=4) and this is very close to the upper limit of 5 (a-1), so it would be unwise to search artificially for a fifth a priori comparison.


 Site details, and group means and sums for the log10 (x+1) transformed E. coli data (n=4)

Site	Group Mean	Group  Sum	Site details
1 2 3 4 5 6	0.595 1.143 1.983 0.887 0.075 0.325	2.380 4.517 7.993 3.547 0.310 1.301	<<<most rocky site tributary enters between 3 and 4
Four replicate samples were taken from a river at each of six sites .  The numbers of E. coli  were estimated  in each sample taken.


Factorial a Posteriori
A POSTERIORI TESTS

  Again the one-way data set will be used to illustate these methods  E. coli data

  These tests are less stringent than the a priori tests but allow us to look for means and groups of means which appear to be different in the data after it has collected.  This feature makes them very powerful and they can be used extensively.  The basic rules for most a posteriori tests are as follows

the initial ANOVA on all the data must be significant before these tests are attempted

large numbers of comparisons are allowed, always more than (a (a - 1))/2.

  The limit set above is a lower limit and essentially allows all possible pairs to be compared.  As we are rarely interested in absolutely all pairwise comparisons, there are essentially comparisons available for some groups of means to be compared also.

The common methods used are

Tukey-Kramer

Sum of squares, simultaneous test proceedure


Factorial Tukey-Kramer
TUKEY-KRAMER METHOD

  The best a posteriori pairwise comparison is the Kramer modification of the Tukey method.  With this test a minimum significant difference (MSD) is calculated which is analogous to the LSD described in the a priori tests.

where Q is the studentized range statistic and so for the E. coli data Q0.05(6,18) is 4.495.

Hence

and so none of the group means below the tributary are different from each other.  We can also use this statistic to investigate the E. coli numbers above the tributary in more detail.  It is clear that site 1 is only individually significantly different from site 3 and not from site 2.  This does not invalidate our previous
a priori conclusion because it is perfectly possible for a priori and a posteriori tests to give conflicting results, because they have markedly different powers.


Factorial SS-STP
SUM OF SQUARES, SIMULTANEOUS TEST PROCEDURE(SS-STP)

  This is a very powerful a posteriori test for comparing groups of means.  The SStest is calculated and then compared with a critical value (SScrit) calculated as follows

SScrit = (a - 1).MSerror.FP(a-1,dferror)

If SStest > SScrit then the comparison being made is significant.  Thus for the E. coli data at P = 0.05.

SScrit = 5 x 0.176 x 2.77 = 2.438

this statistic can be used to further investigate groups of means in the E.col data.  A SStest is calculated longhand or by computer in an identical manner to that previously described for orthogonal decomposition .


Factorial Mtb Table
SUM OF SQUARES SIMULTANEOUS TEST PROCEDURE

Minitab Commands

  The table shows the method of labelling groups for calculating sums of squares by computer for the a posteriori SS-STP test used to compare multiple groups of means.  The log10 E. coli data are used.

	log10	sites	lowest
Missing values are labelled with an asterisk.   The SStest calculated from the column headed 'lowest' would compare the group with the lowest group mean (site 5) below the tributary with the other two below tributary sites (4 and 6).	0.903 0.699 0.000 0.778 1.792 1.146 0.000 1.633 1.707 2.193 2.029 2.003 1.079 0.845 0.477 1.146 0.000 0.000 0.000 0.301 0.301 0.699 0.301 0.000	1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6	* * * * * * * * * * * * 1 1 1 1 2 2 2 2 1 1 1 1
Again the  ONEWAY  command this time using the 'lowest' labels is used.

  In this case SStest was 0.751 and so the difference was not significant.  However, using the same principles a test between sites 1-4 and sites 5 and 6 was significant (SStest = 4.831).


Multiple introduction
MULTIPLE REGRESSION

Introduction

  Sometimes data are collected on two or more different predictor variables x1, x2...xk and a equation is sought to calculate y from the set of measurements of x1, x2...xk.

  The calculations required for multiple regression are usually two time-consuming to be done by hand so it is necessary to use a computer.  Fortunately, most statistical packages have regression routines and these give the user a wide choice of statistics that help in finding suitable prediction equations.

  Many of these diagnostic statistics are most easily interpreted by the use of plots so a package with built-in plotting routines should be used.

There are three main uses of multiple regression.

 finding a prediction equation for y.

 investigating the variables x that influence y

 demonstrating a physical law


Data structure for Multiple Regression
MULTIPLE REGRESSION DATA STRUCTURE

  As a first  simple illustrative example the loblolly data (used in Factorial ANOVA) will be used to show how a polynominal equation can be fitted to a set of data.

Loblolly pine data

To illustrate multiple regression analysis techniques in detail a second data set will be used.

Dipper survey data

  The aim in collecting the data was to establish the nature of the variable influencing the breeding of British dippers, thrush-sized birds in the upper reaches of rivers.  They feed on bethnic invertebrates by probing the river bed with their breaks, measurements were made of 8 variables at 22 stretches of river inhabited by dippers.

  A third data set will be used to describe the techniques of fitting different equations to sub-sets of a bivariate data set by using dummy variables and the methods of multiple regression.

Plaice survey data

  These data  are derived from surveys of the numbers of O-group plaice in Swedish and Dutch coastal waters and the calculated mortality rate assuming an exponential rate of decline in the numbers of fish.  Cautionary notes


Multiple data structure eg 1
DATE STRUCTURE EXAMPLE 1

Numbers of seedlings produced by loblolly pine seeds on six dates (in the format used by Minitab).

Row 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48	seedling 900 880 810 1100 760 960 1040 1040 880 1050 1170 1240 1060 1110 910 1120 1530 1140 1160 1270 1390 1320 1540 1080 1970 1360 1890 1510 1820 1490 2140 1270 1960 1270 1670 1380 1310 1500 1480 1450 830 150 420 380 570 420 760 270	date 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6	date-sq 1 1 1 1 1 1 1 1 4 4 4 4 4 4 4 4 9 9 9 9 9 9 9 9 16 16 16 16 16 16 16 16 25 25 25 25 25 25 25 25 36 36 36 36 36 36 36 36	The dependent variable y is the number of seedlings.  There are two predictors, the date of sowing (x1) represented as 1-6 inclusive, and the squares of these numbers (x2).  The data set for this regression therefore occupies three columns, numbers of seedlings are ordered by date. An analysis of the data (Factorial ANOVA) suggested the possibility of fitting a parabola to predict the number of seedlings with the date of sowing.


Multiple Data Structure eg 2
DATA STRUCTURE EXAMPLE 2

Data from a survey of dippers (in the format used by Minitab).

Row	Alt	WHard	Grad	Caddis	Stonefly	Mayfly	Other In	Br.Dens
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22	259 198 251 184 145 145 198 160 251 159 160 145 214 191 167 244 83 175 107 305 152 277	12.20 22.00 26.30 22.50 29.50 39.90 42.80 59.60 69.25 68.00 85.28 84.18 89.00 103.90 90.90 95.10 153.90 150.50 15.50 115.80 130.70 90.00	10.90 14.70 6.90 4.60 1.91 5.00 6.20 14.30 4.60 3.30 15.80 5.90 5.70 7.60 17.40 6.10 2.88 12.80 4.90 26.20 8.50 19.00	9 17 43 82 24 50 38 83 43 85 230 133 91 60 52 37 16 35 48 254 224 383	188 76 186 188 45 133 77 564 223 33 76 505 92 30 281 165 20 325 120 520 166 889	0 29 338 313 246 313 214 0 124 415 534 484 277 705 972 523 364 273 431 520 742 696	3 4 3 3 2 20 3 18 12 31 144 33 317 873 130 29 338 2 71 24 733 14	3.60 4.30 3.80 3.40 3.80 4.50 4.30 5.00 4.50 3.40 4.50 6.90 3.50 3.60 4.70 4.00 2.90 6.40 4.00 8.46 5.00 8.00
Data provided by Dr. S.J. Ormerod.

  Multiple regression will be used to derive an equation to predict the breeding density of dippers with some or all of the seven other variables used as predictors, possibly after suitable transformations have been made.


The first column ('row') are row numbers each of which represent separate sites on the R. Wye and its tributaries inWales.  The next three columns are environmental variables


'Alt' - Site altitude  'm' above sea level


'w hard' - water hardness, mg CaCO3/l


'grad' = gradient or slope of the river bed, m/km.


Abundances of invertebrates in riffles and marginal areas expressed as numbers from 2 x 3 min kick samples ('caddis' = caddis fly larvae)


'stonefly' = stonefly larvae


'mayfly' = mayfly larvae


'Other In' = all other invertebrates collected.


'Br. Dens'  is the number of breeding pairs of dippers per 10 km of river.


Multiple Data Structure eg 3
DATA STRUCTURE EXAMPLE 3	page 1 of 2

Data from surveys of O-group plaice

Row	M	log Dmax	Zone	x1	x2
1 2 3 4 5 6 7 8 9 10 11 12 13	0.01052 0.00889 0.00697 0.01376 0.00745 0.01013 0.00850 0.01310 0.01265 0.01082 0.02087 0.03474 0.05147	5.09987 5.25227 5.30827 5.69373 5.76832 6.21461 6.25190 6.26340 6.39693 6.68461 7.20042 7.78322 7.90101	2 2 2 2 2 2 2 1 1 1 1 1 1	1 1 1 1 1 1 1 0 0 0 0 0 0	0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 6.26340 6.39693 6.68461 7.20042 7.78322 7.90101	The column (M) is the daily mortality rate (\d) The column (log Dmax) is the natural logarithm of the maximum density (numbers/thousand square metres). The last three columns are explained when these data are analysed.
A scatter plot of the data can be done to investigate the relationship between mortality and density.


Multiple Data Structure eg 3 plot
DATA STRUCTURE EXAMPLE 3	page 2 of 2

Minitab Commands and Output

PLOT C1 C2	[C1 - 'M',  C2 - 'log Dmax']

 M       -

         -                                                     *

         -

    0.045+

         -

         -

         -                                                   *

         -

    0.030+

         -

         -

         -                                         *

         -

    0.015+                *

         -      *                  *  *   *

         -         *               2

         -         *       *

           +---------+---------+---------+---------+---------+------log Dmax

        4.80      5.40      6.00      6.60      7.20      7.80

  The scatter plot of these data suggests that one possibility to describe the relationship between these two variables is a constant mortality up to a threshold maximum density, then a straight line relationship of mortality with log density.  (A curve is another possibility)


Multiple cautionary notes
DATA STRUCTURE CAUTIONARY NOTES

  Where the regression is used for the demonstration of a physical law or for prediction, care should be exercised in generalising from the sample to populations that have not been studied, however similar to the survey sample these populations are believed to be.  It is particularly dangerous to extrapolate from the equation beyond the range of measurement of the variables.

  As in bivariate regression it is desirable that the distribution of the variables should be even.  Equations derived from data mainly concentrated at one end of the scales of measurement are usually heavily influenced by the rate extreme data, and are not therefore truly representative.

  It is preferable that the scales of measurement of the predicted variables be roughly commensurable (similar numerical orders of magnitude).  Round error in calculations can otherwise be a problem in multiple regression, and the extent to which computer packages can deal with non-commensurable variables depends on the way the numerical methods of calculations are programmed.


Multiple route through regression
ROUTE THROUGH A MULTIPLE REGRESSION

1.	Plot y against each x variable in turn

2.	Check assumptions of the multiple regression

3.	Check for evenness in the variables

4.	Transform uneven variables

5.	Compute the regression

6.	Residual analysis

7.	Influence statistics

8.	Check for multi-colinearity

9.	Reduce the number of predictor variables if one or more are not needed.

10	Calculate confidence intervals or prediction intervals for prediction equations

For more details see


Multiple protocol
MULTIPLE REGRESSION PROTOCOL

1.	Plot y against each x variable in turn.  Curvature indicates the need for transformations.  Linear plots indicate important predictors.

2.	Check assumptions of the multiple regression.

3.	Check for evenness in variables using

 boxplots

 letter value displays

 histograms

4.	Transform uneven variables

 standardizing variables

 partial residual plots

 the ladder of powers

5.	Compute the regression.  Save the studentized deleted residuals, leverages and Cook's D.  Calculate the Durbin-Watson statistic, particularly if data are in time order

6.	Analyse the residuals

7.	Check for influential data influence statistics

8.	Check for multicolinearity

9.	Investigate the possibility of reducing the number of predictor variables

10.	When an optimum equation is obtained assess its value for the purposes for which it was intended.  If it is used for prediction, calculate confidence intervals or prediction


Multiple data plotting
1. PLOTTING THE DATA

  Because there are now several potential predictor (x) variables, there are more plots to examine to decide the scale of measurement to be used and whether transformations should be applied.

  However crude the initial indicators for transformations are, it is recommended that these plots be examined at the outset of the analysis, since a more refined model is a better starting-point for an examination of the data.

  A major advantage of an initial scanning of the data is that obvious errors can be rectified from the start and anomalous data may be identified.


Multiple Data Plotting-mtb
PLOTTING THE DATA

Minitab, commands and output

For the loblolly pine seedling data

PLOT C1 C2	[C1 - 'seedling', C2 - 'date']

        -

    2100+                                   *

        -                                   2         *

seedling-                                   *

        -                                             *

        -                         2         2         2

    1400+                         *         *         2

        -               *         2         *         2

        -               5         3

        -     3         *

        -     3         *                                       *

     700+     *                                                 *

        -                                                       *

        -                                                       3

        -                                                       *

        -     *                                                 *

       0+

          ----+---------+---------+---------+---------+---------+--date

            1.0       2.0       3.0       4.0       5.0       6.0

  In these plots curvature indicates the need for transformations and linear plots indicate important predictors.

  Of course for the dipper data there are many more plots which can be examined and although a little time-consuming these plots can be very useful.  They allow an initial scrutiny of the data to check the relationships between the many variables.


Multiple Assumptions
2. CHECKING THE ASSUMPTIONS OF A MULTIPLE REGRESSION

The assumptions underpinning multiple regression are

 the x's are not random variables, and are measured without error

 the errors [image] have a mean of zero

 the errors are independent

 the errors have the same variance,

 the errors are normally distributed

  It is rarely practicable to assume that the x's are error-free and that all errors are in the measurement of y alone although models can be fitted to data with errors in both y and x.  Non-constant error variances can be allowed for, either by transforming the variables or by weighted regression.

  The assumption that the errors are normally distributed is used to establish procedures for statistical testing and the derivation of confidence intervals.

Normality of errors (residuals)

Homogeneity of Variance


Multiple normality of errors
NORMALITY OF ERRORS

  The assumption of normality can be checked by plotting the studentized residuals against the normal scores.

  The studentized residuals are assessed by the TRESIDUAL subcommand of the REGRESS command, and the normal scores by the NSCORE command in Minitab.

They can then be plotted and a straight line indicates agreement with normality, a curve indicates transformation may be needed.

  Unfortunately this check can be unreliable for small sample sizes because of a property called supernormality.

  The predictions can be written as a linear combination of the data and the central unit theorem shows that these will tend to be normally distributed even when the errors are not normal.  It n is large the first term y in the residual will dominate, and supernormality is not so much of a problem.


Multiple normality example
NORMALITY OF ERRORS

Example

  A plot of the studentized deleted residuals against thier normal scores for the dipper data indicates a transfomation is needed.

nsc     -                                                        *

        -

        -                                            *

     1.2+                                            *

        -                                            *

        -                                           **

        -                                         **

        -                                      * *

     0.0+                                    **

        -                                * *

        -                             **

        -                           2

        -                          *

    -1.2+                        *

        -                 *

        -

        -     *

          ----+---------+---------+---------+---------+---------+--Tresids

           -3.0      -2.0      -1.0       0.0       1.0       2.0

  The dependant variable, y, the breeding density of dippers, should be investigated, or the model should be further refined.


Multiple Homogeneity of Variance
HOMOGENEITY OF VARIANCE

  If the variance of the residuals is homogeneous the plots of the residuals should exhibit a random scatter.

Various patterns of residuals plots such as funnels or curves point to non constant variance.

The studentized or standardized residuals should be plotted against the

 fitted values

 x variables

 the sequence of observations if in a logical order.

These plots are examined along with the absolute magnitude of the ordinary residuals.


Multiple evenness of x & y
3. EVENNESS OF X AND Y

  For multiple regression, histograms or boxplots can be performed on the variables to check for evenness. Using the dipper survey data the Minitab commands and output is:

BOXPLOT C1	[C1 - variable to be plotted]

                              --------------------

                 -------------I     +            I------------

                              --------------------

          +---------+---------+---------+---------+---------+------Alt

         50       100       150       200       250       300

               ---------

           ----I +     I------            **  O                  O

               ---------

          +---------+---------+---------+---------+---------+------Caddis

          0        70       140       210       280       350

  Each variable in turn can be plotted and the evenness checked.  Regression analysis is far more likely to calculate an accurate model if the values of x and y are evenly spread through the range.  if a variable was very uneven with several outside and/or for outside values on one side of the boxplot only (B), then transformation might be necessary.

  Letter value displays can also be used in which the midsummaries should be equal in a perfectly symmetrical distribution.  (See one-way ANOVA Assumptions for further details).


Multiple Transformations
4. TRANSFORMATIONS

  It has not proved possible to give a rigorous strategy for determining transformations that linearize the data in the multivariate case.  Plots will sometimes indicate the need for transformations, although the precise nature of the optimum transformation is not necessarily obvious.

  There are many plots to examine in multiple regression, and the plot of y against one of the predictors may show considerable scatter even when a regression equation with two or more predictors is a satisfactory model for the data.  The detection of curvature in these plots is thus more difficult than in the bivariate case.

  Nevertheless, an examination of these plots is recommended as a first stage in multiple regression analysis, since they can often point the way to a suitable transformation.

Standardizing variables

Partial residual plots

The Ladder of Powers

Conclusions


Multiple Standardizing Variables
TRANSFORMATIONS - STANDARDIZING VARIABLES

  Transformation of a predictor variable by dividing by or subtracting a constant have no material effect on the multiple regression calculations.

 Standardization of the predictor variables, by subtracting the mean and dividing by the standard deviation, is one such transformation and this is a way of ensuring that the variables are commensurable.  Standardization of the dependent variable y is not recommended, however.  The division by the standard deviation alters and considerably complicates the results used for tests and confidence intervals.

  Any scaling or subtraction of a constant form y should be done using numbers that are not calculated from the data if the assumptions of the model are not to be compromised.


Multiple Partial Residual Plots pg 1
TRANSFORMATIONS PARTIAL RESIDUAL PLOTS	page 1 of 4

  A diagnostic method that has proved useful in the identification of appropriate transformations of the predictor variables is based on another type of residual called the partial residuals.  The idea behind these is to eliminate the (unclear) effect of all but one of the predictors in calculating residuals and then to plot these against the remaining predictor variable.

The calculation of these partial residuals is based on results from a regression package.

From the prediction equation

the prediction of y without using one of the predictors (xj) is [image] - bjxj (b = constant)

So the partial residual, which is the residual calculated from the prediction using all predictors excepting xj is:

(y - [image]) is the raw residual for the full regression equation.


Multiple Partial Residual plots pg2
TRANSFORMATIONS PARTIAL RESIDUAL PLOTS	page 2 of 4

Using the regression output to calculate the partial regression coefficient:

1.	Calculate a full regression using all predictor variables.

2.	Multiply the column of values of one predictor (xj) by the partial regression coefficient (bj) using the LET command.  [This value (bj) is in the COEF column of the regression output].

3.	Add the column of raw residuals (stored by the RESIDUALS sub-command) to this column of multiplied values, (again with LET).

4.	Plot this final column against the original predictor column values (xj)

eg

LET C8 = C2 * (bj)	<<  value from regression output

LET C9 = C8 + C7	C2 contains data from a chosen predictor variable

PLOT C9 C2	C7 contains the residuals from the regression


Multiple Partial Residual Plots pg 3
TRANSFORMATIONS - PARTIAL RESIDUAL PLOTS	page 3 of 4

Using the Dipper survey data:

         -                                 *

         -                                  *

      2.4+                                                         *

         -

 Pres    -

         -                     *

         -

      1.2+

         -         * *   *                    *

         -    * 2  *

         -   *

         -             2     *

      0.0+           * *

         -   *  **

         -

         -  *

         -

           +---------+---------+---------+---------+---------+------Stonefly

           0       160       320       480       640       800

  This is the plot for the stonefly partial residual from a regression with log (breeding density) as the dependent variable and all seven predictor variables left untransformed.

  Although there is a lot of scatter in the plot there is a distinct suggestion of curvature in the plot, indicating the need for a transformation of the stonefly numbers.  This can be compared to the plot of the partial residual for the transformed stonefly data  from the regression with all the four invertebrate numbers predictors transformed with the loge(1 + x) transformation, (1 + x since zero values are present).  This is much closer to a straight line, indicating that the correct transformation has been found.


Partial residual plot for the transformed stonefly data

         -                                                *

     1.00+

         -                                                *       *

 Presid  -

         -

         -                                    *             *

     0.75+                             *            *

         -                     *      *   * *

         -                     *

         -              *                 *       *

         -        *               *         *

     0.50+                     *

         -         *

         -                                  *

         -  *

         -

     0.25+

           --------+---------+---------+---------+---------+--------logeSt

                3.50      4.20      4.90      5.60      6.30


Multiple Partial Residual Plots pg 4
TRANSFORMATIONS PARTIAL RESIDUAL PLOTS	page 4 of 4

  The plot below is an example of a partial residual plot for a variable that is not important in the regression equation.

         -                    *

 Presalt -                                                    *

         -

         -

     0.70+                    *     *

         -                    *          *

         -            *        *       * *

         -                                         *

         -

     0.00+

         -                                         *

         -                       *                   *

         -                       *          *     *     *

         -                        *   *

    -0.70+

         -        *              *

         -

           +---------+---------+---------+---------+---------+------Alt

          50       100       150       200       250       300

??????The scatter of the altitude partial residual, is evident.  The slight negative?? slope evident in this plot is because the partial regression coefficient for altitude is negative.


Multiple Ladder of Powers pg 1
TRANSFORMATIONS The Ladder of Powers	page 1 of 2

  Boxplots and letter value displays show the distribution of the data as skewed to the left or right, or symmetrical.  Power transformations from the Ladder of Powers can be used moving along the ladder according to the shape of the distribution obtained.  Data skewed to the right (with a long right tail) will use lower powers and left skewed data (with a long left tail) will need higher powers.

  Alternatively the transformation whose normal score plot for the studentized deleted residual most closely follows a straight line can be chosen.

  The normal scores and studentized residuals are stored using the NSCORE command and TRESIDUAL sub-command of the REGRESS command in Minitab.

  A range of transformations are tried and regression equations fitted using the same predictors in each case.  The plots can then be examined.


Multiple Ladder
TRANSFORATIONS -  The Ladder of Powers

  The ladder puts the raw data and a range of equivalent transformations into a logical sequence from which a suitable transformation can be selected.

	Transformation		Addition of small  number needed for zero values
Power	Name	Formula
3 2 1 0.5 0 -0.5 -1 -2	cube square raw data square root logarythm reciprocal root reciprocal reciprocal square	x log10x 1/ 1/x 1/	no no no   no* yes yes yes yes

Based on a table in Velleman and Hoaglin

* Not needed but [image] or [image] is often used.


Reference

Velleman P.F. and Hoaglin, D.C.  (1981).  Applications, Basics, and Computing of Exploratory Data Analysis.  Duxberry Press, Boston.


Multiple Ladder of Powers pg 2
TRANSFORMATIONS - The Ladder of Powers	page 2 of 2

  Using the Dipper data a range of transformations of breeding density as dependent variable and the log transformed predictor variables.

  These and other transformation plots indicate that the y variable should not be transformed, contrary to the indications given by initial plots of the data.

  The plot that is closest to a straight line is that corresponding to untransformed y data, although there is little to choose between this and the other plots.  The ladder of powers technique is not sufficiently sensitive to identify a precise transformation to use.  It merely indicates a range of possibilities.


Multiple Conditions-transformation
TRANSFORMATIONS CONCLUSIONS

  In addition to these methods any prior knowledge of the nature of the inter-relationships between predictor variables and the dependent variable should be used in determining the form of the variables used in the multiple regression equation.  If, for example, it is known from scientific principles that y depends additively on the order of magnitude of x then the logarithm of x should be used.  (Further details based on experience or theory in one-way transformations).

  Although the techniques of transformations outlined are strongly recommended, the user of multiple regression should beware of over-refinement of the regression equation.  Only in exceptional circumstances, where a thorough (and time-consuming) analysis is justified and exhaustive use is recommended of all the techniques described.

  Extensive use of involved transformations and deletion of inconvenient outliers and lever points will often give a convincing multiple regression equation that is highly tailored to the data set used in the calculations when used on data collected outside the initial survey, the equation performs poorly.  One way of guarding against this is to deliberately collect extra data that will not be used in the calculations, but solely to justify the relationships established by the equation.


Multiple Computing the Regression
5. COMPUTING THE REGRESSION

  As a first simple example the loblolly pine data will be used to show how a polynomial equation can be fitted to a set of data.

Loblolly Data Regression

  A detailed description of the interpretation of the output from the REGRESS command is given for the Dipper data.

Dipper Survey Data Regression


Multiple computing the regression mtb pg 1
COMPUTING THE REGRESSION	page 1 of 2

Minitab Commands and Output

For the loblolly pine data:  REGRESS C1 2 C2 C3	(C1-'seedling', C2-'date', C3-'date*2')

(notice the necessity for specifying the number of predictors in this case 2)

Predictor	Coef	Stdev	t-ratio	p
Constant	-18.9	188.7	-0.10	0.921
date	913.9	123.4	7.40	0.000
date-sq	-133.10	17.26	-7.71	0.000
s =298.3	R-sq =57.1%	R-sq(adj)=55.2%

Analysis of Variance

SOURCE	DF	SS	MS	F	p
Regression	2	5335739	2667870	29.98	0.000
Error	45	4003908	88976
Total	47	9339648
SOURCE	DF	SEQ SS
date	1	44464
date-sq	1	5291275

Unusual Observations

Obs.	date	seedling	Fit	Stdev.Fit	Residual	St.Resid
31	4.00	210	1507.1	64.3	632.9	2.17R
33	5.00	1960	1223.0	58.4	737.0	2.52R

R denotes an obs. with a large st. resid.

Here the regression equation is significant. (P < 0.0005).

The data together with the fitted equation can also be plotted.


Multiple Computing the Regression mtb pg 2
COMPUTING THE REGRESSION	page 2 of 2

Minitab Commands and Output

A plot of the data with the fitted equation.

Fig 3??

In this regression the simple parabola model has explained a substantial proportion of the variation due to differences in the dates.


Multiple Computing the Regression eg 2 pg 1
COMPUTING THE REGRESSION	page 1 of 2

Minitab Commands and Output

For the Dipper data: 	(loge(1+x) used for the four invertebrate variables and loge for bird density)

REGRESS C9 7 C1-C3 C10-C13;	[C9 - 'Bird Density', C1-C3 - 'Alt,W hard,Gradient' C10-C13-'Inverts']

TRESID C14;	[C14-empty for the Studentized deleted residuals]

XPXINV M1.	[Matrix M1]

The regression equation is

log BrDn =0.440-0.000606 Alt+0.00159 W Hard+0.0152 Gradient	+0.07405 log Cadd+0.134 log Stfl

+0.0216 log Mayf	-0.0472 log OInv

Predictor	Coef	Stdev	t-ratio	p
Constant	0.4402	0.2271	1.94	0.073
Alt	-0.0006062	0.0006611	-0.92	0.375
W Hard	0.0015877	0.0009272	1.71	0.109
Gradient	0.015152	0.006688	2.27	0.040
log Cadd	0.07049	0.05128	1.37	0.191
log Stfl	0.13428	0.04657	2.88	0.012
log Mayf	0.0212	0.02133	1.01	0.328
log OInv	-0.04725	0.02464	-1.92	0.076
s =0.1347	R-sq =84.7%	R-sq(adj)=77.1%

Analysis if Variance

SOURCE	DF	SS	MS	F	p
Regression	7	1.40763	0.0109	11.08	0.000
Error	14	0.25408	0.01815
Total	21	1.66171
SOURCE	DF	SEQ SS
Alt	1	0.23906
W Hard	1	0.20079
Gradient	1	0.42897
log Cadd	1	0.15212
log Stfl	1	0.29181
log Mayf	1	0.02815
log OInv	1	0.06675

Unusual Observations

Obs.	Alt	log BrDn	Fit	Stdev.Fit	Residual	St.Resid
4	184	1.2238	1.5083	0.0637	-0.2845	-2.40R

R denotes an obs. with a large st. resid

The PREDICT SUBCOMMAND can be used to obtain Confidence and Predition Intervals (see Step 10)

Fit	Stdev.Fit	95% C.I.	95% P.I.
1.0051	0.1518	(0.6794, 1.3309)	(0.5696, 1.4406) X
1.4986	0.0287	(1.4370, 1.5603)	(1.2031, 1.7941)
1.9502	0.1046	(1.7258, 2.1746)	(1.5843, 2.3161)

X denotes a row with x values away from the center


Multiple Computing the Regression eg 2 pg 2
COMPUTING THE REGRESSION	page 2 of 2

Minitab Commands and Output

  The first four lines of output give the regression equation.  It is difficult to interpret an equation involving so many predictors and it is not possible to plot because of the large number of dimensions.  A further difficulty stems from the possibility of inter-relationships between the predictors.  The magnitude of a partial regression coefficient for a variable depends on which other predictors are included in the equation.

  The matrix referred to in Minitab as XPXINV is used in the calculation of the partial regression coefficients.  It was calculated by invoking the XPXINV subcommand.  Its usefulness to the practitioner is that when all of its elements are multiplied by the estimate of the error variance [image], the result obtained is a matrix of variances and covariances of the coefficients (the estimated variance-covariance matrix).


Multiple Residual Analysis
 6. RESIDUAL ANALYSIS

  One way of assessing the effectiveness of a multiple regression equation in predicting the dependent variable is to make a simple comparison of the observations y with the predictions [image].  The differences between them are called the raw residuals, and any large residual is an indication of an unusual value of y.

  A statistic based on the raw residuals is the sum of their squares the error sum of squares.  This is used in the analysis of variance procedure as a measure of variation not explained by the model.  Modified residuals, described later in this section are more sensitive for checking for unusual y values.

Standardized residuals and studentized deleted residuals


Multiple Standardized & Studentized residuals
STANDARDISED RESIDUALS AND STUDENTIZED DELETED RESIDUALS

  One reason for the insensitivity of the raw residuals is that their variance is not constant over the range of measurement of the predictor variables.  Allowance can be made for this non-constant variance when standardizing the raw residuals by dividing them by their standard deviations. (In Minitab these are called standardized residuals).

  Another refinement recognises the fact that an extreme outlier may influence the estimates of the partial regression coefficients and also the error variance.  This approach calculates the residual by subtracting from the observed y the prediction using the data with ith individual omitted, and standardizing by dividing by the estimated standard deviation from this depleted date set.  These are called in Minitab studentized deleted residuals and are assessed by the TRESIDUALS subcommand of the REGRESS command.

  Studentized deleted residuals are recommended for use in residual analysis (and checking the assumptions of the regression) since they are most sensitive in checking for outliers

Checks for Outliers

Checks for missing predictors

Checks on autocorrelation


Multiple checking for outliers pg 1
RESIDUAL ANALYSIS	page 1 of 3

Checking for Outliers

  A simple check for outliers is to consider all data with a studentized residual exceeding + 2 as an indication of an unusual y value.

  A more refined procedure is based on the fact that these residuals are distributed as student's t with
(n - p'- 1) df

  A simple comparison with the student's t tables is likely to lead to false declaration of outliers.  An approximation is to use the Bonferroni correction.  Here the extreme residual is compared with the 100(1 - [image] /(2n)) percentile of student's t with (n - p' - 1) df.  Minitab & SAS allow the values from the table to be calculated easily.


p' is the number of constants, including [image] o, so p' = k + 1


A false declaration of outliers is likely since with n residuals a multiple test procedure is needed to test for the significance of the maximum.


 is the significance level, n  is the number of residuals


Multiple checking for outliers pg 2
RESIDUAL ANALYSIS	page 2 of 3

Checking for Outliers - Worked Example

  When the Dipper data is analysed the 4th observation is indicated by Minitab as an outlier because the standardized residual is -2.40, thus greater than +2.

  The more refined benferrani-corrected t-test with a significance level of 0.05, uses a critical value the 100(1-0.05144) = 99.89 percentile of Student's t with (22 - 8 - 1) = 13 df.

In Minitab this can be done with the INVCDF command (see calculating confidence intervals for more details).

  The critical value obtained is 3.8018, and the studentized deleted residual for the 4th observation is smaller than this in absolute magnitude (it is -3.01), so the observation is not judged on this more conservative test for identification of an outlier.

There is probably no reason therefore to exclude observation 4 from the equation at this stage.


Multiple checking for outliers pg 3
RESIDUAL ANALYSIS	page 3 of 3

Checking for Outliers - Worked Example

  The plots for the studentized deleted residuals are illustrated by a typical example below where they are plotted against loge (1 + stonefly).

        -

        -

     2.0+                                                *

        -

Tresids -

        -        *     *      2      **                  *

        -                                *   *                   *

     0.0+                                  *       *

        -         *              *       *

        -  *                  *                            *

        -                                  *

        -

    -2.0+                                        *

        -

        -

        -                                  *

        -

          --------+---------+---------+---------+---------+--------logeSt

               3.50      4.20      4.90      5.60      6.30

  There is no pattern indicating an incorrect relationship or non-constant variance, though the suspected outlier stands out as the lowest residual.


Multiple checking for missing predictors
RESIDUAL ANALYSIS

Checking for Missing Predictors

  Plots of the studentized deleted residuals against the index number of the data i or the order in which data are collected are valuable in that any relationship may indicate a lurking variable, some missing predictor from the model.  (Index i is the order in which the data sets are entered into the data file.  Use a logical ordering where such is available).

  Plots against the predictions are also useful in this regard.  Plots against the predictor variables may indicate by curved relationships the need for transformations or the inclusion of extra terms such as quadratic, with [image] as well as x.


Multiple Autocorrelation checking
RESIDUAL ANALYSIS

Checking for Autocorrelation

  One diagnostic test that is based on the raw residuals is the Durbin-Watson test for serial autocorrelation in the residuals.

The formula for the Durbin-Watson test statistic is

(ri is the ith raw residual)

and is calculated in Minitab with the DW subcommand of the REGRESS command.  A low test statistic indicates a positive association in the error of one observation with the error of the next.  This can be indication that there is a missing variable in the model, or it may be that there is some genuine autocorrelation in the errors and the assumption of independent errors is untenable.

  If the data are in the order of increasing value of a predictor variable x, a significant value of the Durbin-Watson test statistic may result from curvature in the relationship between y and x.  A transformation of the y variable, or inclusion of powers of x as well as x as predictors, may be enough to obtain independent residuals.  So that time-dependent correlations are picked up by the Durbin-Watson statistic it is necessary for the data to be entered row by row in order of measurement in time.


Multiple Influence Stats
7. INFLUENCE STATISTICS - CHECKS ON INFLUENTIAL DATA

  There are statistics which identify observations with unusual values of the predictor variables.  The most common statistics of this type calculated by statistical packages are leverage (also called HI in Minitab, and LEVER in SPSS(PC+), Cook's D and Dfits.

  All three statistics should be calculated and compared with critical values.  If is common for several different data points to be identified as having high influence by these statistics.  Thus at least two of the statistics should be above the critical value before they are considered truly influential.  Once identified such influence points should be considered carefully for possible exclusion.

Leverages

Cook's D and Dfits

Plots based on influence statistics


Multiple leverages
INFLUENCE STATISTICS - LEVERAGES

  Leverage (hii) are based on the diagonal elements of a matrix called the hat matrix.  They measure the joint influence of the predictors by a standardized distance of the ith observation to the other n-1.

  A rough indication of an unusual observation is a leverage exceeding 3p'/n, though 2p'/n has also been suggested. (p' = the number of parameters in the model, n = number of data points in the regression).

  Minitab uses 3p'/n as the criterion for automatically displaying an unusual observation, but checking when hii exceeds 2p'/n is a prudent precaution against errors.


Multiple Cook's D and Dfits
INFLUENCE STATISTICS = COOK'S D AND DFITS

  A disadvantage of the leverages is that they do not distinguish between high leverage points that are influential in the calculation of partial regression coefficients and those that are not.  A measure that is more sensitive to such influential points of Cook's D.

  The idea behind this measure is to compare the estimates of the regression coefficients from the full data set with those when the ith observation is omitted, so that a large value indicates an influential point.  An indication of a large D value is where it exceeds F0.5(p',n-p')df, (p' = number of parameters, n = number of data points).  Another, but cruder, criterion is to check observations if D exceeds 4/n, and this is more cautious.

  Another statistic Dfits is a measure of the square root of p' times Cook's D with the sign of the studentized deleted residual attached so the two measures are essentially the same.  The criterion to check observations is if Dfits exceeds the absolute value of [image].


Multiple Influence Stats plots pg 1
PLOTS BASED ON INFLUENCE STATISTICS	page 1 of 3

  As well as the indications already mentioned for leverage and Cook's D that an observation may need checking it is useful to plot these statistics against the index number i (the order in which the data sets are entered into the data file.  Use a logical ordering where such is available).  This may indicate an unusual sequence of data that the rough checks do not identify.

  Another check for normality is to plot the square root of cook's D with the sign of the studentized deleted residual attached, against the normal scores.


Multiple Influence Stats plots pg 2
PLOTS BASED ON INFLUENCE STATISTICS	page 2 of 3

Minitab Commands and Output

For the plot of the root of Cook's D against the normal scores for the dipper data:-

LET C25=SQRT(C22) * SIGN(C21)	   [C25 - empty, C22 - Cook's D, C21 - studentized deleted residual]

NSCORE C25 C23	   [ C23 - empty]

PLOT C25 C23

        -                                                        *

        -

    0.35+

        -                                                  *

nscore  -                                      *  * *  *

        -                                 * * *

        -                            * **

    0.00+

        -                         * *

        -                     * *

        -                 *  *

        -               *

   -0.35+

        -            *

        -        *

        -

        -  *

          --------+---------+---------+---------+---------+--------ssqrtCD

              -1.40     -0.70      0.00      0.70      1.40

(Studentized deleted residuals are stored in C21 using TRESIDUAL subcommand and the Cook's D statistic are stored in C22 using the COOKD Subcommand of the REGRESS command).  It is not yet clear whether this plot is more or less sensitive than the residual plot for deleting departures form normality.


Multiple Influence stats plots pg 3
PLOTS BASED ON INFLUENCE STATISTICS	page 3 of 3

Minitab Commands and Output

  The two plots below are the index plots of the leverages hii and Cook's D.  The bunched up appearance of the latter is because Cook's D is always positive.

  No unusual leverage is identified, the largest value of Cook's D is 0.3899, and comparison with F0.5(p',n-p') = 0.964 (or 4/n) shows that this is not an unusually high value.  (p'-no.of parameters D-no.of data points).  None of the values of hii exceed 3 x 8/22 (3p'/n) because this number is greater than 1, and all hii have to be less than or equal to 1.  However, none exceed the less conservative criterion of 2 x 8/22 (2p'/n) = 0.7273 either.  So influential point are not a problem with this regression.


Multiple Multicollin
8. MULTICOLLINEARITY

  Multicollinearity is the presence of near linear relationships amongst some of the predictor variables, that is, one of the predictors is almost equal to a linear combination of other predictors.  If one of the predictors is exactly equal to the linear combination of some of the others, then it could be removed from the regression equation and is redundant as a predictor.

  Even if the linear relationship is not exact, close linear inter-relationships can cause great difficulties in interpretation though the calculation of the regression equation is still possible.  Thus it is important to recognise multi-collinearity and take appropriate steps to eliminate its effects.


Multiple Multicollinearity mtb
MULTICOLLINEARITY

  The variance inflation factor (VIF) can routinely be calculated in multiple regression packages.  Minitab calculates the VIF ( a sub-command of the REGRESS command) associated with each predictor and automatically drops from the regression any predictor where this is greater than [image] when the VIF is greater than 100, a warning of collinearity is printed, but no predictors are excluded.

  Although  VIF is valuable in identifying that variables are multicollinear, it does not help in finding the nature of the inter-relationships.  A method that is often useful in this respect is based on a Principal Components Analysis (PCA) of the standardized predictor variables.  PCA is generally used as an ordination technique, its use to investigate multi-collinearity is a somewhat special application of the technique.

Principal Components Analysis


Multiple PCA pg 1
MULTICOLLINEARITY PRINCIPAL COMPONENTS ANALYSIS	page 1 of 3

  PCA can be done using standard package such as Minitab, SAS or SPSS/PC+.  Since standardized values are used the PCA is done on the correlation matrix.

The Minitab commands and output are (for the dipper data)

PCA C1-C3  C14 - C17	[C1 - C3 raw variables C14-C17 - transformed variables]

Eigenanalysis of the Correlation Matrix

Eigenvalue	2.3765	2.2886	0.7670	0.6017	0.4624	0.3296	0.1741
Proportion	0.340	0.327	0.110	0.086	0.066	0.047	0.025
Cumulative	0.340	0.666	0.776	0.862	0.928	0.975	1.000
Variable	PC1	PC2	PC3	PC4	PC5	PC6	PC7
Alt	0.397	0.317	-0.169	0.666	0.361	-0.374	0.004
W Hard	0.242	-0459	0.421	0.358	-0.542	-0.206	0.298
Gradient	0.537	0.101	0.363	0.059	0.141	0.732	-0.102
log Cadd	0.481	-0.270	-0.323	-0.433	0.256	-0.090	0.574
log Stfl	0.503	0.252	-0.038	-0.392	-0.417	-0.340	-0.488
log Mayf	0.099	-0.479	-0.667	0.270	-0.192	0.291	-0.348
log OInv	0.033	-0.559	0.334	-0.104	0.529	-0.270	-0.460

  The raw labelled eigenvalue gives the variance of the principal components of the standardized predictors.  The sum of these is always equal to the number of predictors.  If one or more eigenvalues is small, then collinearity amongst the predictors is indicated.

  A rough test is to calculate the square root of the ratio of the largest eigenvalue to the smallest to give a number called the condition number.  For these data the condition number is 3.69 and the smallest eigenvalue is 0.17.  Only if the condition number exceeds about 10 is there any real indication of possible multicollinearity so with these data multicollinearity is not a problem.


Multiple PCA pg 2
MULTICOLLINEARITY PRINCIPAL COMPONENTS ANALYSIS	page 2 of 3

  Another data was used in the following PCA to demonstrate collinearity.  The data are four measurements made on 20 soil samples, % sand, % silt, % organic matter and pH.Soil data

PCA C1-C4	[C1-C4 contain % sand, silt & organic matter]

Eigenanalysis of the Correlation Matrix

Eigenvalue	2.0229	1.0505	0.9058	0.0209
Proportion	0.506	0.263	0.226	0.005
Cumulative	0.506	0.78	0.995	1.000
Variable	PC1	PC2	PC3	PC4
% Sand	0.684	-0.151	0.151	0.698
% Silt	-0.696	0.086	-0.052	0.711
% OrgMat	-0.215	-0.580	0.781	-0.084
pH	-0.048	-0.796	-0603	0.005

  The fourth eigenvalue is small, and the condition number for the correlation matrix is

  Clearly this component is dominated by the coefficients of 0.698 for % sand and 0.711 for % silt and a plot of these variables will be close to a straight line.  This in any regression either % sand of % silt should be used, but probably not both.


Row	%Sand	%Silt	%Organic Matter	pH
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20	77.3 82.5 66.9 47.2 65.3 83.3 81.6 47.8 48.6 61.6 58.6 69.3 61.8 67.7 57.2 67.2 59.2 80.2 82.2 69.7	13.0 10.0 20.6 33.8 20.5 10.0 12.7 36.5 37.1 25.5 26.5 22.3 30.8 25.3 31.2 22.7 31.2 13.2 11.1 20.7	1.5 1.5 2.3 2.8 1.9 2.2 2.9 2.3 2.1 1.9 2.4 4.0 2.7 4.8 2.4 3.3 2.4 2.0 2.2 3.1	6.4 6.5 7.0 5.8 6. 7.0 6.7 7.2 7.2 7.3 6.7 7.0 6.4 7.3 6.5 6.2 6.0 5.8 7.2 5.0	Data from samplesof soil


Multiple PCA pg 3
MULTICOLLINEARITY  PRINCIPAL COMPONENTS ANALYSIS	page 3 of 3

  If prediction is the aim of the regression analysis, collinearities should be removed by omitting predictor variables that are highly collinear with others.  If the regression is used to construct models to investigate the influential variables on a predictor, collinearities should be detected using the Protocol below and care exercised in interpreting the partial regression coefficients.  Different combinations of predictors should be tried, eliminating collinearities using relationships established by PCA.

Multicollinearity Investigation Protocol

1.	Perform a PCA on the correlation matrix calculated from the predictor variables.

2.	Calculate the condition number [image] where E1, and E2 are respectively the largest and smallest eigenvalues.  If this exceeds 30, a high degree of multicollinearity is present.  If it exceeds 10, multicollinearity may still need investigation.

3.	Examine the coefficients of the eigenvector (or principal component) corresponding to the smallest eigenvalue.  Ignore small coefficients.  The large coefficients indicate the variables that are collinear, and the coefficients are the multiples of the standardized variables whose sum is approximately zero.

Other approaches (biased methods of regression)

Ridge Regression

Principal component regression and latent root regression


Eigenvalues are sometimes called characteristic values or latent values


Eigenvectors are sometimes called characteristic vectors or latent vectors


Details are given by Chatterjee and Price (Chapter 8) and Gunst and Mason (chapter 10).

Ref?


Details are given by Gunst and mason (Chapter 10)

ref?


Multiple variable selection
9. SELECTION OF VARIABLES

  The selection of variables in multiple regression is based on the assessment of the variation in the dependent variable, in the sense of sum of squares explained by the regression.  The principle is to include in the regression those variables that markedly decrease the unexplained or error sum of squares.

The steps of variable selection.

1.	Identify an optimal subset of variables using a procedure that calculates the regression statistics for the best regression with a specified number of predictors.

2.	Plot Mallows Cp against q + 1 for a range of subsets of q predictors and look for one whose Cp is close to q + 1.  Ignore those subsets with most predictors included unless other subsets have a Cp much greater than q + 1.

3.	Confirm that the optimal subset is best using a stepwise procedure that both adds and removes variables to find a solution.

4.	Calculate the full regression analysis for the optimal subset and follow the multiple regression protocol.

Criteria for variable selection

Automatic selection procedures

Dummy variables in multiple regression


Multiple criteria for variable selection
CRITERIA FOR VARIABLE SELECTION

There are three related criteria for the inclusion (or exclusion) of variables in regression, all based on the sums of square.

The F-test

The coefficient of determination,

Mallow's Cp


Multiple variable selection F-test
CRITERIA FOR VARIABLE SELECTION       The F-test

  The formal F-test statistic is based on the difference between SSerror,k when all k predictor variables are included in the regression and SSerror,q when only q of the variables are included. (SSerror is the sum of squares from the error line of the ANOVA table for the relevant regression equations).

The F-test statistic is calculated by

  In practice this statistic is sometimes significant, indicating that variables should be included in the regression equation, when these variables make no practical difference to the performance of the equation.


Multiple variable selection R-sq
CRITERIA FOR VARIABLE SELECTION         The Coefficient of Determination

  A second criterion is to compare the coefficient of determination [image] from the subset of q predictors, with that obtained from all k predictors, [image].   is always bigger than [image] if the q predictors are a subset of the k.

  Unfortunately [image] is not sensitive and can remain high even when practically important predictors are removed from the regression equation.


Multiple variable selection Cp
CRITERIA FOR VARIABLE SELECTION       Mallow's Cp

  The third criterion is Mallow's Cp defined as

  If the subset of q predictors is equally as effective in predicting y as the full set of k predictors then SSerror,k/(n-k-1) is an estimate of the error variance [image] and SSerror,q is an estimate of (n-q-1) times the same error variance, so Cp is roughly equal to q+1;  where q+1 = p' the number of constants in the solved multiple regression equation.

  If Cp is much greater than Q+1 then this indicates that SSerror, Q is much greater than (n-q-1) [image] and that important variables have been omitted from the regression equation.

One approach is to try several subsets of a list of potential predictors and plot the values of Cp against Q+1, choosing a subset that gives a Cp close to the one Cp = Q+1.  A sudden jump in this plot as q decreases indicates that too many variables have been omitted from the regression equation.

  The theoretical average of Cp is q+1 if all important variables are included in the equation.  Cp can be less than q+1 in value though it is much greater if important variables are left out of the equation.  Choice of the subset with lowest Cp may lead to an unnecessarily complicated equation for practical purposes.


Multiple Automatic selection of variables
AUTOMATIC SELECTION PROCEDURES

  Various procedures have been suggested for searching through a set of potential predictors equations to decide which should be included in the final regression equation.  One is to try all possible regressions, but as the number of predictors increases it may be unpracticable to examine such a large number of regressions.

  A second, and workable, method is to decide how many predictors should be used and find the subset with the specified number of predictors that is best in the sense of smallest error sum of squares.  Algorithms are available in computer packages using this criterion, and in conjunction with plots of mallow's Cp the method can be very helpful in deciding the predictors to include.

Other methods are based on an approach in which the number of predictors is changed by one at each step.

Forward Selection

Backward Elimination

Stepwise Procedure


In forward selection, an initial selection is made of the best predictor, then the second predictors which is in conjunction with the first is included, then a third and so on. A disadvantage is that variables included earlier in the analysis may be ineffective in larger sets, because of inter-relationships with other predictors, but no mechanism is included for the omission of a variable once it is included.  It is a useful method where the predictor variables are powers of a single variable, and thus a polynomial equation in x is used as a predictor for y.


The reverse process to forward selection is backward elimination.  Here the starting point is an equation with all variables included, and the first step is to eliminate the variable that causes the smallest increase in the error sum of squares, then the variable that again adds least to prediction and so on.


A compromise is the stepwise procedure, a combination of forward and backward selection.  This starts as forward selection but includes a mechanism for the possible elimination of a variable at each step this therefore needs two criteria, an F for a variable to enter the equation and an F for a variable to be removed.  This stepwise procedure is the best readily-available, standard automatic selection method.


Multiple variables seleciton exmaple pg 1
SELECTION OF VARIABLES IN MULTIPLE REGRESSION	page 1 of 2

Using the dipper data:

Bird density (y) - untransformed dependent variable

Altitude, water hardness and gradient (x1-x3) - three intransformed predictors caddis fly, stonefly, mayfly and other invertebrates (x4-x7) - four loge(1+x) transformed predictors.

  Using the first step is to calculate Mallow's Cp for a range of prediction equations and plot these values against the number of predictors plus one.

  The command BREG (for best regression) in Minitab enables a selection of these to be obtained very rapidly.  The dependent variable and all potential predictor variables are entered on the command line.  The package then calculates the coefficient of determination [image] and Mallow's Cp for the best and next best single predictor equations, in the sense of highest [image] highest F-ratio or smallest Cp.

  The best and next best regressions are found for two predictor variables and information printed.  The process continues until all predictors entered on the command line are printed.


Multiple Variable selection exmaple pg 2
SELECTION OF VARIABLES IN MULTIPLE REGRESSION	page 2 of 2

Minitab Commands and Output

BREG C8 C1-C3 C14-C17	[C8-'B.Dens' C1-C3 & C14-C17-un/transformed variables]

							G	l	l	l	l
							r	o	o	o	o
						W	a	g	g	g	g
							d
						H	i	C	S	M	O	Minitab Commands for
					A	a	e	a	t	a	I
		Adj.			l	r	n	d	f	y	n
Vars	R-sq	R-sq	C-p	s	t	d	t	d	l	f	v	Stepwise Procedure
1	58.2	56.1	19.3	0.98492					x			Backward Elimination
1	50.4	48.0	26.3	1.0726			x					Conclusions
2	69.2	66.0	11.5	0.86708			x		x
2	69.0	65.7	11.7	0.87081		x			x
3	77.5	73.8	6.1	0.76136			x		x	x
3	75.3	71.2	8.1	0.79836		x	x		x
4	79.3	74.4	6.5	0.75220		x	x		x	x
4	78.4	73.3	7.3	0.76804		x		x	x		x
5	83.0	77.6	5.2	0.70294		x	x	x	x		x
5	81.6	75.8	6.4	0.73078		x	x		x	x	x
6	84.0	77.7	6.3	0.70283		x	x	x	x	x	x
6	83.2	76.5	7.0	0.72104	x	x	x	x	x		x
7	84.3	76.5	8.0	0.72086	x	x	x	x	x	x	x
  Unfortunately the information calculated in BREG with Minitab cannot be stored automatically, and the plot of Cp against the number of predictors plus one requires these data to be entered with the keyboard.


Multiple variable selection eg pg 3
STEPWISE PROCEDURE	page 1 of 3

  The stepwise procedure which allows both addition and removal of variables can also be performed in Minitab using the stepwise command.  The default values for the F to enter and F to remove a predictor are both equal to 4.

STEPWISE C8 C1-C3 C14-C17	[C8-'B.dens' C1-C3 & C14-C17 un/transformed variables]

STEPWISE REGRESSION OF Br. Dens ON 7 PREDICTORS, WITH N=22

STEP	1	2	3
CONSTANT	-1.0076	-0.1727	-1.7747
log Stfl	1.14	0.79	0.86
T-RATIO	5.28	3.41	4.20
Gradient		0.096	0.096
T-RATIO		2.61	2.97
log Mayf			0.232
T-RATIO			2.58
S	0.985	0.867	0.761
R-SQ	58.21	69.23	77.53
Each column of the output represents a step in the procedure.

  Notice that with the default values the steps are terminated with the 3-variable solution identified with BREG.  The three predictors were entered in turn, and none were removed at any step.  the last line of output is a subcommand inviting the user to change the default values of the F to enter or F to remove, or to force variables in or out of the regression equation.

  Forward selection and backward elimination can be done in Minitab by changing the defaults.


Multiple variable selection stepwise pg 2
STEPWISE PROCEDURE	page 2 of 3

  The F to remove (FREEMOVE) is the criterion for removal of a predictor.

STEPWISE C8 C1-C3 C14-C17;	[C8-'B.Dens' C1-C3 & C14-C17 un/untransformed]

FREEMOVE=0;

FENTER=.1.

STEPWISE REGRESSION OF Br. Dens ON 7 PREDICTORS, WITH N=22

STEP	1	2	3	4	5	6	7
CONSTANT	-1.0076	-0.1727	-1.7747	-1.8706	-1.2301	-1.1121	-0.8806
log Stfl	1.14	0.79	0.86	0.88	0.79	0.59	0.60
T-RATIO	5.28	3.41	4.20	4.34	3.75	2.44	2.42
Gradient		0.096	0.096	0.084	0.087	0.074	0.080
T-RATIO		2.61	2.97	2.50	2.66	2.26	2.24
log Mayf			0.232	0.180	0.201	0.111	0.115
T-RATIO			2.58	1.81	2.06	1.00	1.00
W Hard				0.0053	0.0087	0.0101	0.0100
T-RATIO				1.20	1.77	2.09	2.01
log OInv					-0.16	-0.25	-0.27
T-RATIO					-1.42	-2.01	-2.03
log Cadd						0.41	0.41
T-RATIO						1.52	1.49
Alt							-0.0018
T-RATIO							-0.51
S	0.985	0.867	0.761	0.752	0.731	0.703	0.721
R-SQ	58.21	69.23	77.53	79.28	81.60	84.04	84.33

The small increase in [image] achieved by the inclusion of more than three predictors indicates that the three predictor model is probably the best that can be achieved with these data. In this output it is the t-ratio that should be examined.


The F to remove (FREEMOVE) is the crterion for removal of a prdictor, if the F-ratio for removal of any predictor is less than FREEMOVE, the variable with the smallest F-ratio is removed. By setting this to zero therefore, the removal of variables is made impossible. The FENTER value is set very small so as to illustrate more steps in the forward selection routine.


Multiple variable selection stepwise pg 3
STEPWISE PROCEDURE	page 3 of 3

  The [image] value can be an unreliable guide in variable selection and since STEPWISE does not print Mallow's Cp it is the t-ratio that should be examined.

  A formal test can only be done by first calculating the degrees of freedom (n-q-1) where q is the number of predictors in the model.  At step 3 the t-ratio for the last variable introduced, loge(1+Mayfly numbers), is 2.58 and the df are (22-3-1) = 18.  The 5% critical value is 2.101 so this t-ratio is significant.

  The t-ratio for water hardness at the next step is not significant at the 5% level (nor indeed at the 10% level), so this test indicates that water hardness should not be included with the best three predictors.

  Notice that the t-ratio for loge(1+Mayfly numbers) has fallen to 1.81 at this step, but this is the t-ratio for the inclusion of this variable if water hardness is included as well, so this statistic is not relevant.


Multiple variable selection backward elimination
BACKWARD ELIMINATION

  Backward elimination is carried out in Minitab with the STEPWISE procedure with all the predictor variables entered at the first step and FENTER criterion set very large.

STEPWISE C8 C1-C3 C14-C17;	[C8-'B.Dens' C1-C3 & C14-C17-un/tansformed variables]

ENTER C1-C3 C14-C17;

FENTER=10000;

FREEMOVE=10.

STEPWISE REGRESSION OF Br. Dens ON 7 PREDICTORS, WITH N=22

STEP	1	2	3	4	5	6	7
CONSTANT	-0.8806	-1.1121	-0.7244	-1.5080	-2.3755	-1.7599	-1.0076
Alt	-0.0018
T-RATIO	-0.51
W Hard	0.0100	0.0101	0.0118	0.0137	00087	0.0114
T-RATIO	2.01	2.09	2.61	2.83	1.91	2.57
Gradient	0.080	0.074	0.065
T-RATIO	2.24	2.26	2.07
log Cadd	0.41	0.41	0.55	0.62	0.39
T-RATIO	1.49	1.52	2.43	2.55	1.65
log Stfl	0.60	0.59	0.51	0.71	0.97	1.12	1.14
T-RATIO	2.42	2.44	2.24	3.13	4.70	5.89	5.28
log Mayf	0.11	0.11
T-RATIO	1.00	1.00
log OInv	-0.27	-0.25	-0.27	-0.28
T-RATIO	-2.03	-2.01	-2.16	-2.05
S	0.721	0.703	0.703	0.768	0.834	0.871	0.985
R-SQ	84.33	84.04	82.97	78.40	73.05	68.97	58.21

Defects of the precedure


 If the default value of FREEMOVE is left at 4, the procedure stops at stepe 3 with five variables are predictors.  However when FREEMOVE is set to 10 all variables will be eliminated.  This output illustrates one of the defects of this procedure.

  The predictor variable loge(1+mayfly numbers) is eliminated at step 3 but using both BREG and STEPWISE (with default FENTER value of 4) this variable is included in the best solutions.  From step 4 onwards the solutions form the backward elimination procedure are sub-optimal.


Multiple variable selection conclusions
SELECTION OF VARIABLES IN MULTIPLE REGRESSION

Conclusions

  Forward selection often has the same defect as backward elimination, a variable included early leads to sub-optimal solutions with more predictors at subsequent steps.  The stepwise procedure to some extent ameliorates this problem, but even with this method optimal solutions are not guaranteed.  For that reason a combination of BREG and STEPWISE should be used, with careful scrutiny of Mallow's Cp.

  In the 3-variable model identified by the selection procedures all variables included are statistically significant. this equation may therefore form the basis of a useful investigation, though that does not mean it is optimal.

  The search for optimal transformation and equations is difficult, however, the solution proffered here is simple and fairly effective.  it highlights some of the important influences on dipper breeding density.  It is, however, highly tailored to this particular data set and may not be representative of other populations of dippers.


Multiple dummy variables pg 1
DUMMY VARIABLES IN MULTIPLE REGRESSION	page 1 of 4

  Qualitative factors can be included in regression models by arbitrarily ascribing numerical values to the classifications to obtain dummy or indicator variables.  In a regression study on animals, for example, the sec of each animal could be incorporated in the study using a dummy variable that takes the value 1 for females and 0 for males.  Dummy variables are also used to enable unbalanced factorial experiments to be analyzed.

  The so-called general linear model commands in packages such as Minitab and SAS calculate the necessary dummy variables automatically.


Multiple dummy variables pg 2
DUMMY VARIABLES IN MULTIPLE REGRESSION	page 2 of 4

  The use of dummy variables can be illustrated using the mortality data.  These data are the mortality rates and the logarithm of the maximum density for 0-group plaice.  Each row of data relates to a particular year class.

  The scatter plot of the data suggests that an equation that remains constant up to a threshold value for the logarithm of maximum density equal to about 6.26 and thereafter increases.  Linearly would be a sensible model for mortality.

 M       -

         -                                                     *

         -

    0.045+

         -

         -

         -                                                   *

         -

    0.030+

         -

         -

         -                                         *

         -

    0.015+                *

         -      *                  *  *   *

         -         *               2

         -         *       *

           +---------+---------+---------+---------+---------+------log Dmax

        4.80      5.40      6.00      6.60      7.20      7.80


Data from surveys of O-group plaice

Row	M	log Dmax	Zone	x1	x2
1 2 3 4 5 6 7 8 9 10 11 12 13	0.01052 0.00889 0.00697 0.01376 0.00745 0.01013 0.00850 0.01310 0.01265 0.01082 0.02087 0.03474 0.05147	5.09987 5.25227 5.30827 5.69373 5.76832 6.21461 6.25190 6.26340 6.39693 6.68461 7.20042 7.78322 7.90101	2 2 2 2 2 2 2 1 1 1 1 1 1	1 1 1 1 1 1 1 0 0 0 0 0 0	0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 6.26340 6.39693 6.68461 7.20042 7.78322 7.90101	The column (M) is the daily mortality rate (\d) The column (log Dmax) is the natural logarithm of the maximum density (numbers/thousand square metres).


Multiple dummy variables
DUMMY VARIABLES IN MULTIPLE REGRESSION	page 3 of 4

  In fact there is some justification for the use of such an equation.  The lower densities, zone 2 of the data correspond to surveys made in the Dutch Wadden Sea.  Zone 1 consists of bays in Swedish waters.  It is possible that there is a relationship between the variables for the Swedish data, but not for the Dutch.

  Variable x is simply a dummy variable constructed according to the rules on page 1 of this section.  Whilst x2 is constructed as zeros for zone 2, to represent the constant, ? value assumed in this zone and the values of logDmax for zone 1.


Multiple dummy variables
DUMMY VARIABLES IN MULTIPLE REGRESSION	page 4 of 4

  The multiple regression output for these data:

The regression equation is  M = -0.123 + 0.133x1 + 0.029x2

Predictor Constsnt x1 x2	Coef -0.12349 0.13295 0.020947	Stdev 0.02262 0.02270 0.003200	t-ratio -5.46 5.86 6.55	p 0.000 0.000 0.000
s= 0.005022    R-sq = 87.5%     R-sq(adj) = 84.9%

Analysis of Variance

SOURCE Regression Error Total	DF 2 10 12	SS 0.00175798 0.00025219 0.00201017	MS 0.0008799 0.00002522	F 34.85	p 0.000
SOURCE x1 x2	DF 1 1	SEQ SS 0.00067755 0.00108043
(Unusual Observations are also outputted)

  For the Dutch data the multiple regression equation gives as a prediction y = -0.12349 + 0.13295 = 0.00946 and this is equal to the mean of the y values of the Dutch data.  For the Swedish date the multiple regression equation gives as a prediction y = -0.112349 + 0.020947x.

  This is simply the regression equation of y on x for the Swedish data alone.  The advantage of this dummy variable method of fitting separate equations to the two parts of the data are that a comparison can be made with the simpler model in which y and x are assumed to be related by a single straight line.  Similarly, comparison could be made with some assumed curve.

The multiple regression equation can be plotted.

Plot in Minitab


Plot of the multiple regression equation for the O-group plaice mortality data

fig 24


Multiple tests and C.I's
10. TESTS AND CONFIDENCE INTERVALS IN MULTIPLE REGRESSION

  When an optimum regression equation is obtained assess its value for the purposees for which it was intended.  It is used for prediction, calculate confidence intervals or prediction intervals.

Overall tests of significance for the regression coefficients

Significance tests and confidence intervals for individual regression coefficients

Confidence intervals for predictions from the regression equation


Multiple tests of significance pg 1
OVERALL TESTS OF SIGNIFICANCE FOR THE	page 1 of 4

REGRESSION COEFFICIENTS

  As has been mentioned already the criterion for determining the optimum multiple regression equation is to make the sum of squares as small as possible.  The principle of multiple regression modelling is to include in the equation those predictors that significantly decrease the error sum of squares.

  The regression sum of squares measures the gain in explained variation obtained by suing the regression equation.  It is the difference between the error sum of square and the total sum of squares.  A formal comparison of the regression and error sum of squares is made by first dividing by the degrees of freedom to get the mean squares, then calculating the ratio of these mean squares.

Standard F tables are used to assess the significance of the ratio of these mean squares.


Multiple tests of significance pg 2
OVERALL TESTS OF SIGNIFICANCE FOR THE 	page 2 of 4

REGRESSION COEFFICIENTS

  The formal statistical hypothesis tested by the F-ratio is that all of the coefficients are equal to zero.  A large F-ratio suggests that this hypothesis should be rejected and that at least one, but not necessarily all, of the coefficients is non-zero.

  For the dipper data the F-ratio is 11.08, with 7 and 14 df, and the associated P-value is less than 0.005, indicating that the hypothesis is decisively rejected.  The initial indications are that some linear predictor equation is tenable.


Multiple tests of significance pg 3
OVERALL TESTS OF SIGNIFICANCE FOR THE 	page 3 of 4

REGRESSION COEFFICIENTS

  A related measure of the effectiveness of a multiple regression equation is the coefficient of determination or squared multiple correlation coefficient.  This is designated [image] (or R-sq in Minitab) and is calculated as the ratio of regression sum of squares to total sum of squares.  It is often expressed as a percentage my multiplication by 100.  This [image] measures the proportion of variation, as measured by sum of squares of errors, explained by the prediction equation.

  When an additional predictor variable is included in the regression equation [image] always increases to give the adjusted [image]  Since the adjusted [image] does not have a simple definition as a proportion of explained variability, it is more difficult to interpret practically than the simple


Multiple tests of significance pg 4
OVERALL TESTS OF SIGNIFICANCE FOR THE 	page 4 of 4

REGRESSION COEFFICIENTS

  The overall significance of the multiple regression equation can also be found from the multiple correlation coefficient (R = square root [image]).  The value of R can then be compared to critical table values of significance.

  As measure of the effectiveness of a multiple regression equation [image] has the defect that it is relatively insensitive, particularly where it is used in variable selection.  For the dipper data the [image] value is high, 84.7%, this value should be treated with some caution.  The overall significance of a multiple regression equation depends very much on the use to which it will be put.


Multiple significance tests and CI's pg 1
SIGNIFICANCE TESTS AND CONFIDENCE INTERVALS FOR 	page 1 of 3

INDIVIDUAL REGRESSION COEFFICIENTS

  The table underneath the analysis variance for the dipper data gives the sequential sums of squares (Seq SS) for the variables in the regression equation.

Analysis if Variance

SOURCE	DF	SS	MS	F	p
Regression	7	1.40763	0.0109	11.08	0.000
Error	14	0.25408	0.01815
Total	21	1.66171
SOURCE	DF	SEQ SS
Alt	1	0.23906
W Hard	1	0.20079
Gradient	1	0.42897
log Cadd	1	0.15212
log Stfl	1	0.29181
log Mayf	1	0.02815
log OInv	1	0.06675

  The values of the Seq SS will usually be different if the order in which the predictor variables are included in the command line is changed.  Thus the regression SS for altitude in a simple bivariate regression is 0.23906, with 1df.  The regression SS for altitude and water hardness together is 0.23906 + 0.20079 = 0.43985, with 2df.  The regression SS for water hardness alone is obtained by putting it as the first predictor variable in the REGRESS command line, (it would then be 0.16536 and with extra due to altitude the sum is still 0.43985).

  The sum of all 7 of the Seq SS is equal to the overall regression sum of squares of 1.40763 (in the ANOVA table).


Multiple significance tests and CI's
SIGNIFICANCE TESTS AND CONFIDENCE INTERVALS FOR	page 2 of 3

INDIVIDUAL REGRESSION COEFFICIENTS

  In Minitab the coefficients are tabulated together with their estimated standard deviations after the regression equation.

Predictor	Coef	Stdev	t-ratio	p
Constant	0.4402	0.2271	1.94	0.073
Alt	-0.0006062	0.0006611	-0.92	0.375
W Hard	0.0015877	0.0009272	1.71	0.109
Gradient	0.015152	0.006688	2.27	0.040
log Cadd	0.07049	0.05128	1.37	0.191
log Stfl	0.13428	0.04657	2.88	0.012
log Mayf	0.0212	0.02133	1.01	0.328
log OInv	-0.04725	0.02464	-1.92	0.076
s =0.1347	R-sq =84.7%	R-sq(adj)=77.1%

  The t-ratio also printed here is the standard student's t test statistic for the hypothesis that the constant is zero.  The associated P-value of this statistic is printed alongside.  Thus in only two cases in this analysis is the hypothesis of a zero constant rejected, the two being gradient with t  = 2.27 and P = 0.040 and loge(1 + stonefly) with t = 2.88 and P = 0.012.

  The lack of significance of the remaining t-ratios should not, however, be taken as conclusive evidence that the remaining variables should be excluded from the regression equation.  Because the sum of squares associated with any predictor variable depends on the other variables in the equation, the t-ratios will be different if the multiple regression is alone again with one or more predictors omitted.  In a regression with fewer predictors included, it may well be that another of the t-ratios is significant.


Multiple significance test and CI's pg 3
SIGNIFICANCE TESTS AND CONFIDENCE INTERVALS FOR	page 3 of 3

INDIVIDUAL REGRESSION COEFFICIENTS

  If a confidence interval for an individuals constant is required this is readily calculated from the formula.

where t is a percentile of Student's t distribution with (n - p')df and var(bj) is the estimated variance of bj ( [image] times the jth diagonal element of the XDX1NV matrix).  The square root of this variance is the standard deviation of bj, often called the standard error of the coefficient.  For a 95% confidence interval the 2-tailed critical value with P = 0.05 is needed.

  For the dipper data n = 22 and p' = 8, so the df are 14.  The 2-tailed critical value of student's t with 14df and P = 0.05 is given in tables, and is 2.145.  Thus a 95% confidence interval for the constant bo(0.4402)

is  0.4402 + 2.45 x 0.2271

or  0.4402 + 0..4871

This procedure is analogous to that used in bivariate regression.


n - no.of data points     p' - no.of predictors


Multiple CI's for predictions pg 1
CONFIDENCE INTERVALS FOR PREDICTIONS FROM 	page 1 of 5

THE REGRESSION EQUATION

  The formula for the confidence intervals for predictions is rather complicated.  Fortunately most packages print confidence intervals for predictions directly and the formula is not needed.

Where the leverages hii are calculated by the package the following formula can be used.

  To illustrate the use of this formula, a 95% confidence interval for the single prediction corresponding to the first row of data of the dipper survey data.

Row 1	Alt 259	WHard 12.20	Grad 10.90	Caddis 9	Stonefly 188	Mayfly 0	Other In 3	Br.Dens 3.60


Multiple CI's for predictions pg 2
CONFIDENCE INTERVALS FOR PREDICTIONS FROM	page 2 of 5

THE REGRESSION EQUATION

  The values of the 7 predictors used in this equation are 259, 12.2, 10.9, 2.30259, 5.24175, 0, 1.38629.  Substitution of these values into the regression equation give the predictor ( [image]) as 1.2684.  The correct t-value is the 2-tailed P = 0.05 point of student's t gives in tables and is 2.145 whilst [image] is the error mean square (MSrerror) from the output is 0.01815.

  The value of the leverage, hii for the first observation, was calculated using the HI subcommand of the Minitab REGRESS command.  the value for hii (0.557926) is put into the equation for the 95% confidence interval for the single prediction.

  Sometimes it is useful to be able to calculate confidence intervals at other sets of values of the predictor variables than those used in the regression calculations.  This can be done in Minitab:


Multiple CI's for predictions pg 3
CONFIDENCE INTERVALS FOR PREDICTIONS FROM	page 3 of 5

THE REGRESSION EQUATION

  In Minitab the sub command of the REGRESS command used is PREDICT

REGRESS etc;

PREDICT 85 12 2 3 3 0 1;

PREDICT 189.5 72.58 9.33 4.122 4.977 5.977 5.33 3.332;

PREDICT 300 150 25 6 6 6 6.

  The first set of values of predictor variables for which a prediction and confidence interval is obtained is at approximately the minimum for the data.  Notice that a warning is printed that these values are away form the centre;  some of the specified values are smaller than the minima.

Fit	Stdev.Fit	95% C.I.	95% P.I.
1.0051	0.1518	(0.6794, 1.3309)	(0.5696, 1.4406) X
1.4986	0.0287	(1.4370, 1.5603)	(1.2031, 1.7941)
1.9502	0.1046	(1.7258, 2.1746)	(1.5843, 2.3161)

X denotes a row with x values away from the center

  The second set of values is at the mean and the last set is at approximately the maxima, notice also that the standard deviation of the fit, and hence the width of the confidence interval, is much narrower for the central values than it is for the extremes.This is a general rule for confidence intervals for predictions in multiple regression.


Multiple CI's for predictions pg 4
CONFIDENCE INTERVALS FOR PREDICTIONS FROM	page 4 of 5

THE REGRESSION EQUATION

  The final column of intervals, labelled P.I. for prediction interval, are 95% interval for a further imaginary observation taken at the specified set of values of the predictor values.

  The formula for the variance of the prediction intervals differ from that of the confidence intervals for the mean predicted [image] in that the standard error is calculated slightly differently.  Thus the prediction intervals are somewhat wider, representing the scatter of further observations about the mean prediction [image].

  For a 99% interval the value of student's t with 14d.f. is 2.977, from tables, so the 99% prediction interval at the mean is


Multiple CI's for predictions pg 5
CONFIDENCE INTERVALS FOR PREDICTIONS FROM	page 5 of 5

THE REGRESSION EQUATION

  In these confidence intervals the standard deviation is multiplied by the 100(1- [image] /2) percentile of student's t with (n - p')df to give the 100(1- [image] /2)% confidence interval.  (These values are tabulated as P = [image]).  Thus for a 95% confidence interval, the multiplier is the 97.5th percentile tabulated under P = 0.05.  The confidence coefficient is at the stated level only if one such interval is required.

  If several intervals are required simultaneously, some adjustment is needed to ensure that the confidence intervals are not too narrow.  One method of ensuring that the joint confidence coefficient is at least equal to 100(1 - [image])% is to use the Bonferroni adjustment


If c confidence intervals are required simultaneously, the 100(1 - [image] /(2c)) percentile Student's t is used as the multiplier for the standard deviation.  The Bonferroni adjustment is an approximate method and can be used so long as the c different intervals are calculated for reasonably separated values of the predictor variables, and if the number of intervals c is no bigger than the number of parameters p'.

  For a 95% confidence interval [image] = 0.05, with c = 6 the adjustment requires the 100(1 - 0.05/(2 x 6)) = 99.5833 percentile of Student's t  to be used in place of the 95% percentile with 14d.f.  In Minitab:

INVCDF 0.995833;

T 14.

The percentile is calculated as 3.0687.

Confidence Interval --

1.2684 + 0.3088.