g04ca:: Analysis of Variance (NAG Toolbox)

An experiment consists of a collection of units, or plots, to which a number of treatments are applied. In a factorial experiment the effects of several different sets of conditions are compared, e.g., three different temperatures,

T_{1}

T_{2}

and

T_{3}

, and two different pressures,

P_{1}

and

P_{2}

. The conditions are known as factors and the different values the conditions take are known as levels. In a factorial experiment the experimental treatments are the combinations of all the different levels of all factors, e.g.,

\begin{array}{l} T_{1} P_{1}, & T_{2} P_{1}, & T_{3} P_{1} \\ T_{1} P_{2}, & T_{2} P_{2}, & T_{3} P_{2} \end{array}

The effect of a factor averaged over all other factors is known as a main effect, and the effect of a combination of some of the factors averaged over all other factors is known as an interaction. This can be represented by a linear model. In the above example if the response was

y_{i j k}

for the

k

th replicate of the

i

th level of

T

and the

j

th level of

P

the linear model would be

y_{i j k} = μ + t_{i} + p_{j} + γ_{i j} + e_{i j k}

where

μ

is the overall mean,

t_{i}

is the main effect of

T

p_{j}

is the main effect of

P

γ_{i j}

is the

T \times P

interaction and

e_{i j k}

is the random error term. In order to find unique estimates constraints are placed on the arguments estimates. For the example here these are:

\begin{array}{l} \sum_{i = 1}^{3} {\hat{t}}_{i} = 0, \\ \sum_{j = 1}^{2} {\hat{p}}_{j} = 0, \\ \sum_{i = 1}^{3} {\hat{γ}}_{i j} = 0, & for ​ j = 1, 2 ​ and \\ \sum_{j = 1}^{2} {\hat{γ}}_{i j} = 0, & for ​ i = 1, 2, 3, \end{array}

where

\hat{​ ​}

denotes the estimate.

If there is variation in the experimental conditions (e.g., in an experiment on the production of a material different batches of raw material may be used, or the experiment may be carried out on different days), then plots that are similar are grouped together into blocks. For a balanced complete factorial experiment all the treatment combinations occur the same number of times in each block.

nag_anova_factorial (g04ca) computes the analysis of variance (ANOVA) table by sequentially computing the totals and means for an effect from the residuals computed when previous effects have been removed. The effect sum of squares is the sum of squared totals divided by the number of observations per total. The means are then subtracted from the residuals to compute a new set of residuals. At the same time the means for the original data are computed. When all effects are removed the residual sum of squares is computed from the residuals. Given the sums of squares an ANOVA table is then computed along with standard errors for the difference in treatment means.

The data for nag_anova_factorial (g04ca) has to be in standard order given by the order of the factors. Let there be

k

factors,

f_{1}, f_{2}, \dots, f_{k}

in that order with levels

l_{1}, l_{2}, \dots, l_{k}

respectively. Standard order requires the levels of factor

f_{1}

are in order

1, 2, \dots, l_{1}

and within each level of

f_{1}

the levels of

f_{2}

are in order

1, 2, \dots, l_{2}

and so on.

For an experiment with blocks the data is for block

1

then for block

2

, etc. Within each block the data must be arranged so that the levels of factor

f_{1}

are in order

1, 2, \dots, l_{1}

and within each level of

f_{1}

the levels of

f_{2}

are in order

1, 2, \dots, l_{2}

and so on. Any within block replication of treatment combinations must occur within the levels of

f_{k}

The ANOVA table is given in the following order. For a complete factorial experiment the first row is for blocks, if present, then the main effects of the factors in their order, e.g.,

f_{1}

followed by

f_{2}

, etc. These are then followed by all the two factor interactions then all the three factor interactions, etc., the last two rows being for the residual and total sums of squares. The interactions are arranged in lexical order for the given factor order. For example, for the three factor interactions for a five factor experiment the

10

interactions would be in the following order:

\begin{array}{l} f_{1} f_{2} f_{3} \\ f_{1} f_{2} f_{4} \\ f_{1} f_{2} f_{5} \\ f_{1} f_{3} f_{4} \\ f_{1} f_{3} f_{5} \\ f_{1} f_{4} f_{5} \\ f_{2} f_{3} f_{4} \\ f_{2} f_{3} f_{5} \\ f_{2} f_{4} f_{5} \\ f_{3} f_{4} f_{5} \end{array}

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Let there be

k

factors with levels

l_{i}

for

i = 1, 2, \dots, k

, and let

t

be the maximum number of terms in an interaction then the number of rows in the ANOVA tables is:

\sum_{i = 1}^{t} (\begin{array}{l} k \\ i \end{array}) + 3 .

The number of treatment means is:

\sum_{i = 1}^{t} \prod_{j \in S_{i}} l_{j},

where

S_{i}

is the set of all combinations of the

k

factors

i

at a time.

To estimate missing values the Healy and Westmacott procedure or its derivatives may be used, see John and Quenouille (1977). This is an iterative procedure in which estimates of the missing values are adjusted by subtracting the corresponding values of the residuals. The new estimates are then used in the analysis of variance. This process is repeated until convergence. A suitable initial value may be the grand mean. When using this procedure irdf should be set to the number of missing values plus one to obtain the correct degrees of freedom for the residual sum of squares.

For analysis of covariance the residuals are obtained from an analysis of variance of both the response variable and the covariates. The residuals from the response variable are then regressed on the residuals from the covariates using, say, nag_correg_linregs_noconst (g02cb) or nag_correg_linregm_fit (g02da). The coefficients obtained from the regression can be examined for significance and used to produce an adjusted dependent variable using the original response variable and covariate. An approximate adjusted analysis of variance table can then be produced by using the adjusted dependent variable. In this case irdf should be set to one plus the number of fitted covariates.

For designs such as Latin squares one more of the blocking factors has to be removed in a preliminary analysis before the final analysis. This preliminary analysis can be performed using nag_anova_random (g04bb) or a prior call to nag_anova_factorial (g04ca) if the data is reordered between calls. The residuals from the preliminary analysis are then input to nag_anova_factorial (g04ca). In these cases irdf should be set to the difference between n and the residual degrees of freedom from preliminary analysis. Care should be taken when using this approach as there is no check on the orthogonality of the two analyses.

Example

The data, given by John and Quenouille (1977), is for the yield of turnips for a factorial experiment with two factors, the amount of phosphate with

6

levels and the amount of liming with

3

levels. The design was replicated in

3

blocks. The data is input and the analysis of variance computed. The analysis of variance table and tables of means with their standard errors are printed.

function g04ca_example


fprintf('g04ca example results\n\n');

y = [274; 361; 253; 325; 317; 339; 326; 402; 336; 379;
     345; 361; 352; 334; 318; 339; 393; 358; 350; 340;
     203; 397; 356; 298; 382; 376; 355; 418; 387; 379;
     432; 339; 293; 322; 417; 342;  82; 297; 133; 306;
     352; 361; 220; 333; 270; 388; 379; 274; 336; 307;
     266; 389; 333; 353];

% Sizes
nblock = int64(3);
inter  = int64(2);

% Levels for the 2 factors
lfac   = [int64(6); 3];
irdf   = int64(0);

% Get routine to calculate maxt and mterm
ws = warning;
warning('Off');
warn_state = nag_issue_warnings();
nag_issue_warnings(true);

mterm  = int64(1);
maxt   = int64(0);
[table, itotal, tmean, e, imean, semean, bmean, r, ifail] = ...
  g04ca( ...
         y, lfac, nblock, inter, irdf, mterm, maxt);

nag_issue_warnings(warn_state);
warning(ws);

mterm = itotal;
maxt = imean(1);

% Calculate ANOVA table
[table, itotal, tmean, e, imean, semean, bmean, r, ifail] = ...
  g04ca( ...
         y, lfac, nblock, inter, irdf, mterm, maxt);

% Display results
fprintf('ANOVA table\n\n');
fprintf(' Source      df         SS          MS          F        Prob\n\n');
fmt6 = '%s%2d%5.0f%12.0f%12.0f%12.3f%11.4f\n';
fmt5 = '%s%5.0f%12.0f%12.0f%12.3f%11.4f\n';
fmt3 = '%s%5.0f%12.0f%12.0f\n';
fmt2 = '%s%5.0f%12.0f\n';
k = 0;
if nblock>1
  k = k + 1;
  fprintf(fmt5, 'Blocks    ', table(1,1:5));
end
ntreat = itotal - 2 - k;
for i = 1:ntreat
  fprintf(fmt6, 'Effect  ', i, table(k+i,1:5));
end
fprintf(fmt3, 'Residual  ', table(itotal-1,1:3));
fprintf(fmt2, 'Total     ', table(itotal,1:2));
fprintf('\nTreatment Means and Standard Errors\n\n');
k = 1;
for i = 1:ntreat
  l = imean(i);
  fprintf('Effect %5d\n\n', i);
  for j = k:6:l
    fprintf('%10.2f', tmean(j:min(j+5,l)));
    fprintf('\n');
  end
  fprintf('\nSE of difference in means  = %10.2f\n\n', semean(i));
  k = l + 1;
end

g04ca example results

ANOVA table

 Source      df         SS          MS          F        Prob

Blocks        2       30119       15059       7.685     0.0018
Effect   1    5       73008       14602       7.451     0.0001
Effect   2    2       21596       10798       5.510     0.0085
Effect   3   10       31192        3119       1.592     0.1513
Residual     34       66628        1960
Total        53      222543

Treatment Means and Standard Errors

Effect     1

    254.78    339.00    333.33    367.78    330.78    360.67

SE of difference in means  =      20.87

Effect     2

    334.28    353.78    305.11

SE of difference in means  =      14.76

Effect     3

    235.33    332.67    196.33    342.67    341.67    332.67
    309.33    370.33    320.33    395.00    370.33    338.00
    373.33    326.67    292.33    350.00    381.00    351.00

SE of difference in means  =      36.14

On entry,	$n < 4$ ,
or	$nfac < 1$ ,
or	$nblock < 0$ ,
or	$inter < 0$ ,
or	$inter > nfac$ ,
or	$irdf < 0$ .

On entry,	$lfac (i) \leq 1$ , for some $i = 1, 2, \dots, nfac$ ,
or	the value of maxt is too small,
or	the value of mterm is too small,
or	$nblock > 1$ and n is not a multiple of nblock,
or	the number of plots per block is not a multiple of the number of treatment combinations.

NAG Toolbox: nag_anova_factorial (g04ca)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example