NAG FL Interface
g04bcf (rowcol)

1 Purpose

g04bcf computes the analysis of variance for a general row and column design together with the treatment means and standard errors.

2 Specification

Fortran Interface
Subroutine g04bcf ( nrep, nrow, ncol, y, nt, it, gmean, tmean, tabl, ldtabl, c, ldc, irep, rpmean, rmean, cmean, r, ef, tol, irdf, wk, ifail)
Integer, Intent (In) :: nrep, nrow, ncol, nt, it(*), ldtabl, ldc, irdf
Integer, Intent (Inout) :: ifail
Integer, Intent (Out) :: irep(nt)
Real (Kind=nag_wp), Intent (In) :: y(nrep*nrow*ncol), tol
Real (Kind=nag_wp), Intent (Inout) :: tabl(ldtabl,5), c(ldc,nt)
Real (Kind=nag_wp), Intent (Out) :: gmean, tmean(nt), rpmean(nrep), rmean(nrep*nrow), cmean(nrep*ncol), r(nrep*nrow*ncol), ef(nt), wk(3*nt)
C Header Interface
#include <nag.h>
void  g04bcf_ (const Integer *nrep, const Integer *nrow, const Integer *ncol, const double y[], const Integer *nt, const Integer it[], double *gmean, double tmean[], double tabl[], const Integer *ldtabl, double c[], const Integer *ldc, Integer irep[], double rpmean[], double rmean[], double cmean[], double r[], double ef[], const double *tol, const Integer *irdf, double wk[], Integer *ifail)
The routine may be called by the names g04bcf or nagf_anova_rowcol.

3 Description

In a row and column design the experimental material can be characterised by a two-way classification, nominally called rows and columns. Each experimental unit can be considered as being located in a particular row and column. It is assumed that all rows are of the same length and all columns are of the same length. Sets of equal numbers of rows/columns can be grouped together to form replicates, sometimes known as squares or rectangles, as appropriate.
If for a replicate, the number of rows, the number of columns and the number of treatments are equal and every treatment occurs once in each row and each column then the design is a Latin square. If this is not the case the treatments will be non-orthogonal to rows and columns. For example in the case of a lattice square each treatment occurs only once in each square.
For a row and column design, with t treatments in r rows and c columns and b replicates or squares with n=brc observations the linear model is:
y ijkl = μ+ βi+ ρj+ γk+ τl+ eijk  
for i=1,2,b, j=1,2,,r, k=1,2,c and l=1,2,,t, where βi is the effect of the ith replicate, ρj is the effect of the jth row, γk is the effect of the kth column and the ijkl notation indicates that the lth treatment is applied to the unit in row j, column k of replicate i.
To compute the analysis of variance for a row and column design the mean is computed and subtracted from the observations to give, yijkl=yijkl-μ^. Since the replicates, rows and columns are orthogonal the estimated effects, ignoring treatment effects, β^i, ρ^j, γ^k, can be computed using the appropriate means of the yijkl, and the unadjusted sum of squares computed as the appropriate sum of squared totals for the yijkl divided by number of units per total. The observations adjusted for replicates, rows and columns can then be computed by subtracting the estimated effects from yijkl to give yijkl .
In the case of a Latin square design the treatments are orthogonal to replicates, rows and columns and so the treatment effects, τ^l, can be estimated as the treatment means of the adjusted observations, yijkl . The treatment sum of squares is computed as the sum of squared treatment totals of the yijl divided by the number of times each treatment is replicated. Finally the residuals, and hence the residual sum of squares, are given by rijl=yijl -τ^l.
For a design which is not orthogonal, for example a lattice square or an incomplete Latin square, the treatment effects adjusted for replicates, rows and columns need to be computed. The adjusted treatment effects are found as the solution to the equations:
Aτ^=R-NbNbT/rc-NrNrT/bc-NcNcT/brτ^=q  
where q is the vector of the treatment totals of the observations adjusted for replicates, rows and columns, yijkl , R is a diagonal matrix with Rll equal to the number of times the lth treatment is replicated, and Nb is the t by b incidence matrix, with Nl,i equal to the number of times treatment l occurs in replicate i, with Nr and Nc being similarly defined for rows and columns. The solution to the equations can be written as:
τ^=Ωq  
where, Ω is a generalized inverse of A. The solution is found from the eigenvalue decomposition of A. The residuals are first calculated by subtracting the estimated adjusted treatment effects from the adjusted observations to give rijl=yijl -τ^l. However, since only the unadjusted replicate, row and column effects have been removed and they are not orthogonal to treatments, the replicate, row and column means of the rijl have to be subtracted to give the correct residuals, rijl and residual sum of squares.
Given the sums of squares, the mean squares are computed as the sums of squares divided by the degrees of freedom. The degrees of freedom for the unadjusted replicates, rows and columns are b-1, r-1 and c-1 respectively and for the Latin square designs the degrees of freedom for the treatments is t-1. In the general case the degrees of freedom for treatments is the rank of the matrix Ω. The F-statistic given by the ratio of the treatment mean square to the residual mean square tests the hypothesis:
H0:τ1=τ2==τt=0.  
The standard errors for the difference in treatment effects, or treatment means, for Latin square designs, are given by:
seτ^j-τ^j*=2s2/bt  
where s2 is the residual mean square. In the general case the variances of the treatment effects are given by:
Varτ^=Ωs2  
from which the appropriate standard errors of the difference between treatment effects or the difference between adjusted means can be calculated.
The analysis of a row-column design can be considered as consisting of different strata: the replicate stratum, the rows within replicate and the columns within replicate strata and the units stratum. In the Latin square design all the information on the treatment effects is given at the units stratum. In other designs there may be a loss of information due to the non-orthogonality of treatments and replicates, rows and columns and information on treatments may be available in higher strata. The efficiency of the estimation at the units stratum is given by the (canonical) efficiency factors, these are the nonzero eigenvalues of the matrix, A, divided by the number of replicates in the case of equal replication, or by the mean of the number of replicates in the unequally replicated case, see John (1987). If more than one eigenvalue is zero then the design is said to be disconnected and information on some treatment comparisons can only be obtained from higher strata.

4 References

Cochran W G and Cox G M (1957) Experimental Designs Wiley
Davis O L (1978) The Design and Analysis of Industrial Experiments Longman
John J A (1987) Cyclic Designs Chapman and Hall
John J A and Quenouille M H (1977) Experiments: Design and Analysis Griffin
Searle S R (1971) Linear Models Wiley

5 Arguments

1: nrep Integer Input
On entry: b, the number of replicates.
Constraint: nrep1.
2: nrow Integer Input
On entry: r, the number of rows per replicate.
Constraint: nrow2.
3: ncol Integer Input
On entry: c, the number of columns per replicate.
Constraint: ncol2.
4: ynrep×nrow×ncol Real (Kind=nag_wp) array Input
On entry: the n=brc observations ordered by columns within rows within replicates. That is y rci-1 +rj-1 +k contains the observation from the kth column of the jth row of the ith replicate, for i=1,2,,b, j=1,2,,r and k=1,2,,c.
5: nt Integer Input
On entry: the number of treatments. If only replicates, rows and columns are required in the analysis then set nt=1.
Constraint: nt1.
6: it* Integer array Input
Note: the dimension of the array it must be at least nrep×nrow×ncol if nt>1, and at least 1 otherwise.
On entry: if nt>1, iti indicates which of the nt treatments unit i received, for i=1,2,,n.
If nt=1, it is not referenced.
Constraint: if nt2, 1itint, for i=1,2,,n.
7: gmean Real (Kind=nag_wp) Output
On exit: the grand mean, μ^.
8: tmeannt Real (Kind=nag_wp) array Output
On exit: if nt2, tmeanl contains the (adjusted) mean for the lth treatment, μ^*+τ^l, for l=1,2,,t, where μ^* is the mean of the treatment adjusted observations yijkl-τ^l. Otherwise tmean is not referenced.
9: tablldtabl5 Real (Kind=nag_wp) array Output
On exit: the analysis of variance table. Column 1 contains the degrees of freedom, column 2 the sum of squares, and where appropriate, column 3 the mean squares, column 4 the F-statistic and column 5 the significance level of the F-statistic. Row 1 is for replicates, row 2 for rows, row 3 for columns, row 4 for treatments (if nt>1), row 5 for residual and row 6 for total. Mean squares are computed for all but the total row, F-statistics and significance are computed for treatments, replicates, rows and columns. Any unfilled cells are set to zero.
10: ldtabl Integer Input
On entry: the first dimension of the array tabl as declared in the (sub)program from which g04bcf is called.
Constraint: ldtabl6.
11: cldcnt Real (Kind=nag_wp) array Output
On exit: the upper triangular part of c contains the variance-covariance matrix of the treatment effects, the strictly lower triangular part contains the standard errors of the difference between two treatment effects (means), i.e., cij contains the covariance of treatment i and j if ji and the standard error of the difference between treatment i and j if j<i, for i=1,2,,t and j=1,2,,t.
12: ldc Integer Input
On entry: the first dimension of the array c as declared in the (sub)program from which g04bcf is called.
Constraint: ldcnt.
13: irepnt Integer array Output
On exit: if nt>1, the treatment replications, Rll, for l=1,2,,nt. Otherwise irep is not referenced.
14: rpmeannrep Real (Kind=nag_wp) array Output
On exit: if nrep>1, rpmeani contains the mean for the ith replicate, μ^+β^i, for i=1,2,,b. Otherwise rpmean is not referenced.
15: rmeannrep×nrow Real (Kind=nag_wp) array Output
On exit: rmeanj contains the mean for the jth row, μ^+ρ^i, for j=1,2,,r.
16: cmeannrep×ncol Real (Kind=nag_wp) array Output
On exit: cmeank contains the mean for the kth column, μ^+γ^k, for k=1,2,,c.
17: rnrep×nrow×ncol Real (Kind=nag_wp) array Output
On exit: the residuals, ri, for i=1,2,,n.
18: efnt Real (Kind=nag_wp) array Output
On exit: if nt2, the canonical efficiency factors. Otherwise ef is not referenced.
19: tol Real (Kind=nag_wp) Input
On entry: the tolerance value used to check for zero eigenvalues of the matrix Ω. If tol=0.0 a default value of 0.00001 is used.
Constraint: tol0.0.
20: irdf Integer Input
On entry: an adjustment to the degrees of freedom for the residual and total.
irdf1
The degrees of freedom for the total is set to n-irdf and the residual degrees of freedom adjusted accordingly.
irdf=0
The total degrees of freedom for the total is set to n-1, as usual.
Constraint: irdf0.
21: wk3×nt Real (Kind=nag_wp) array Workspace
22: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value -1 is recommended since useful values can be provided in some output arguments even when ifail0 on exit. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
Note: in some cases g04bcf may return useful information.
ifail=1
On entry, irdf=value.
Constraint: irdf0.
On entry, ldc=value and nt=value.
Constraint: ldcnt.
On entry, ldtabl=value.
Constraint: ldtabl6.
On entry, ncol=value.
Constraint: ncol2.
On entry, nrep=value.
Constraint: nrep1.
On entry, nrow=value.
Constraint: nrow2.
On entry, nt=value.
Constraint: nt1.
On entry, tol=value.
Constraint: tol0.0.
ifail=2
On entry, at least one treatment is not present. Treatment value does not appear in it.
On entry, i=value, iti=value and nt=value.
Constraint: 1itint.
ifail=3
On entry, the values in y are constant.
ifail=4
A computed standard error is zero. This is due to rounding errors and is an unlikely error exit.
The computation of the eigenvalues has failed to converge. This is an unlikely error exit.
ifail=5
The treatments are totally confounded with blocks. The treatment sum of squares and degrees of freedom are zero. The analysis of variance table is not computed, except for block and total sums of squares and degrees of freedom.
ifail=6
The residual degrees of freedom is zero. Columns 3,4 and 5 of the analysis of variance table will not be computed and the matrix of standard errors and covariances, c, will not be scaled by s or s2.
The residual mean square is zero. Columns 3,4 and 5 of the analysis of variance table will not be computed and the matrix of standard errors and covariances, c, will not be scaled by s or s2.
ifail=7
The design is disconnected. The standard errors may not valid. The design may be nested.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

The algorithm used in g04bcf, described in Section 3, achieves greater accuracy than the traditional algorithms based on the subtraction of sums of squares.

8 Parallelism and Performance

g04bcf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
g04bcf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

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, g02cbf or g02daf. The results from those routines can be used to test for the significance of the covariates. To test the significance of the treatment effects after fitting the covariate, the residual sum of squares from the regression should be compared with the residual sum of squares obtained from the equivalent regression but using the residuals from fitting replicates, rows and columns only.

10 Example

The data for a 5×5 Latin square is input and the ANOVA and treatment means computed and printed. Since the design is orthogonal only one standard error need be printed

10.1 Program Text

Program Text (g04bcfe.f90)

10.2 Program Data

Program Data (g04bcfe.d)

10.3 Program Results

Program Results (g04bcfe.r)