PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_anova_rowcol (g04bc)
Purpose
nag_anova_rowcol (g04bc) computes the analysis of variance for a general row and column design together with the treatment means and standard errors.
Syntax
[
gmean,
tmean,
tabl,
c,
irep,
rpmean,
rmean,
cmean,
r,
ef,
ifail] = g04bc(
nrep,
nrow,
ncol,
y,
nt,
it,
tol,
irdf)
[
gmean,
tmean,
tabl,
c,
irep,
rpmean,
rmean,
cmean,
r,
ef,
ifail] = nag_anova_rowcol(
nrep,
nrow,
ncol,
y,
nt,
it,
tol,
irdf)
Description
In a row and column design the experimental material can be characterised by a twoway 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 nonorthogonal 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$t$ treatments in
r$r$ rows and
c$c$ columns and
b$b$ replicates or squares with
n = brc$n=brc$ observations the linear model is:
for
i = 1,2 … ,b$i=1,2\dots ,b$,
j = 1,2, … ,r$j=1,2,\dots ,r$,
k = 1,2 … ,c$k=1,2\dots ,c$ and
l = 1,2, … ,t$l=1,2,\dots ,t$, where
β_{i}${\beta}_{i}$ is the effect of the
i$i$th replicate,
ρ_{j}${\rho}_{j}$ is the effect of the
j$j$th row,
γ_{k}${\gamma}_{k}$ is the effect of the
k$k$th column and the
ijk(l)$ijk\left(l\right)$ notation indicates that the
l$l$th treatment is applied to the unit in row
j$j$, column
k$k$ of replicate
i$i$.
To compute the analysis of variance for a row and column design the mean is computed and subtracted from the observations to give, y_{ijk(l)}^{ ′ } = y_{ijk(l)} − μ̂${y}_{ijk\left(l\right)}^{\prime}={y}_{ijk\left(l\right)}\hat{\mu}$. Since the replicates, rows and columns are orthogonal the estimated effects, ignoring treatment effects, β̂_{i}${\hat{\beta}}_{i}$, ρ̂_{j}${\hat{\rho}}_{j}$, γ̂_{k}${\hat{\gamma}}_{k}$, can be computed using the appropriate means of the y_{ijk(l)}^{ ′ }${y}_{ijk\left(l\right)}^{\prime}$, and the unadjusted sum of squares computed as the appropriate sum of squared totals for the y_{ijk(l)}^{ ′ }${y}_{ijk\left(l\right)}^{\prime}$ 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 y_{ijk(l)}^{ ′ }${y}_{ijk\left(l\right)}^{\prime}$ to give y_{ijk(l)}^{ ′ ′ }${y}_{ijk\left(l\right)}^{\prime \prime}$.
In the case of a Latin square design the treatments are orthogonal to replicates, rows and columns and so the treatment effects, τ̂_{l}${\hat{\tau}}_{l}$, can be estimated as the treatment means of the adjusted observations, y_{ijk(l)}^{ ′ ′ }${y}_{ijk\left(l\right)}^{\prime \prime}$. The treatment sum of squares is computed as the sum of squared treatment totals of the y_{ij(l)}^{ ′ ′ }${y}_{ij\left(l\right)}^{\prime \prime}$ divided by the number of times each treatment is replicated. Finally the residuals, and hence the residual sum of squares, are given by r_{ij(l)} = y_{ij(l)}^{ ′ ′ } − τ̂_{l}${r}_{ij\left(l\right)}={y}_{ij\left(l\right)}^{\prime \prime}{\hat{\tau}}_{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:
where
q$q$ is the vector of the treatment totals of the observations adjusted for replicates, rows and columns,
y_{ijk(l)}^{ ′ ′ }${y}_{ijk\left(l\right)}^{\prime \prime}$,
R$R$ is a diagonal matrix with
R_{ll}${R}_{ll}$ equal to the number of times the
l$l$th treatment is replicated, and
N_{b}${N}_{b}$ is the
t$t$ by
b$b$ incidence matrix, with
N_{l,i}${N}_{l,i}$ equal to the number of times treatment
l$l$ occurs in replicate
i$i$, with
N_{r}${N}_{r}$ and
N_{c}${N}_{c}$ being similarly defined for rows and columns. The solution to the equations can be written as:
where,
Ω$\Omega $ is a generalized inverse of
A$A$. The solution is found from the eigenvalue decomposition of
A$A$. The residuals are first calculated by subtracting the estimated adjusted treatment effects from the adjusted observations to give
r_{ij(l)}^{ ′ } = y_{ij(l)}^{ ′ ′ } − τ̂_{l}${r}_{ij\left(l\right)}^{\prime}={y}_{ij\left(l\right)}^{\prime \prime}{\hat{\tau}}_{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
r_{ij(l)}^{ ′ }${r}_{ij\left(l\right)}^{\prime}$ have to be subtracted to give the correct residuals,
r_{ij(l)}${r}_{ij\left(l\right)}$ 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$b1$,
r − 1$r1$ and
c − 1$c1$ respectively and for the Latin square designs the degrees of freedom for the treatments is
t − 1$t1$. In the general case the degrees of freedom for treatments is the rank of the matrix
Ω$\Omega $. The
F$F$statistic given by the ratio of the treatment mean square to the residual mean square tests the hypothesis:
The standard errors for the difference in treatment effects, or treatment means, for Latin square designs, are given by:
where
s^{2}${s}^{2}$ is the residual mean square. In the general case the variances of the treatment effects are given by:
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 rowcolumn 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 nonorthogonality 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$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.
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
Parameters
Compulsory Input Parameters
 1:
nrep – int64int32nag_int scalar
b$b$, the number of replicates.
Constraint:
nrep ≥ 1${\mathbf{nrep}}\ge 1$.
 2:
nrow – int64int32nag_int scalar
r$r$, the number of rows per replicate.
Constraint:
nrow ≥ 2${\mathbf{nrow}}\ge 2$.
 3:
ncol – int64int32nag_int scalar
c$c$, the number of columns per replicate.
Constraint:
ncol ≥ 2${\mathbf{ncol}}\ge 2$.
 4:
y(nrep × nrow × ncol${\mathbf{nrep}}\times {\mathbf{nrow}}\times {\mathbf{ncol}}$) – double array
The
n = brc$n=brc$ observations ordered by columns within rows within replicates. That is
y(rc(i − 1) + r(j − 1) + k)${\mathbf{y}}(rc(\mathit{i}1)+r(\mathit{j}1)+\mathit{k})$ contains the observation from the
k$\mathit{k}$th column of the
j$\mathit{j}$th row of the
i$\mathit{i}$th replicate, for
i = 1,2, … ,b$\mathit{i}=1,2,\dots ,b$,
j = 1,2, … ,r$\mathit{j}=1,2,\dots ,r$ and
k = 1,2, … ,c$\mathit{k}=1,2,\dots ,c$.
 5:
nt – int64int32nag_int scalar
The number of treatments. If only replicates, rows and columns are required in the analysis then set
nt = 1${\mathbf{nt}}=1$.
Constraint:
nt ≥ 1${\mathbf{nt}}\ge 1$.
 6:
it( : $:$) – int64int32nag_int array

Note: the dimension of the array
it
must be at least
nrep × nrow × ncol${\mathbf{nrep}}\times {\mathbf{nrow}}\times {\mathbf{ncol}}$ if
nt > 1${\mathbf{nt}}>1$, and at least
1$1$ otherwise.
If
nt > 1${\mathbf{nt}}>1$,
it(i)${\mathbf{it}}\left(\mathit{i}\right)$ indicates which of the
nt treatments unit
i$\mathit{i}$ received, for
i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$.
If
nt = 1${\mathbf{nt}}=1$,
it is not referenced.
Constraint:
if
nt ≥ 2${\mathbf{nt}}\ge 2$,
1 ≤ it(i) ≤ nt$1\le {\mathbf{it}}\left(\mathit{i}\right)\le {\mathbf{nt}}$, for
i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$.
 7:
tol – double scalar
The tolerance value used to check for zero eigenvalues of the matrix
Ω$\Omega $. If
tol = 0.0${\mathbf{tol}}=0.0$ a default value of
0.00001$0.00001$ is used.
Constraint:
tol ≥ 0.0${\mathbf{tol}}\ge 0.0$.
 8:
irdf – int64int32nag_int scalar
An adjustment to the degrees of freedom for the residual and total.
 irdf ≥ 1${\mathbf{irdf}}\ge 1$
 The degrees of freedom for the total is set to n − irdf$n{\mathbf{irdf}}$ and the residual degrees of freedom adjusted accordingly.
 irdf = 0${\mathbf{irdf}}=0$
 the total degrees of freedom for the total is set to n − 1$n1$, as usual.
Constraint:
irdf ≥ 0${\mathbf{irdf}}\ge 0$.
Optional Input Parameters
None.
Input Parameters Omitted from the MATLAB Interface
 ldtabl ldc wk
Output Parameters
 1:
gmean – double scalar
The grand mean, μ̂$\hat{\mu}$.
 2:
tmean(nt) – double array
If
nt ≥ 2${\mathbf{nt}}\ge 2$,
tmean(l)${\mathbf{tmean}}\left(\mathit{l}\right)$ contains the (adjusted) mean for the
l$\mathit{l}$th treatment,
μ̂^{*} + τ̂_{l}${\hat{\mu}}^{*}+{\hat{\tau}}_{\mathit{l}}$, for
l = 1,2, … ,t$\mathit{l}=1,2,\dots ,t$, where
μ̂^{*}${\hat{\mu}}^{*}$ is the mean of the treatment adjusted observations
y_{ijk(l)} − τ̂_{l}${y}_{ijk\left(l\right)}{\hat{\tau}}_{l}$. Otherwise
tmean is not referenced.
 3:
tabl(ldtabl,5$5$) – double array
ldtabl ≥ 6$\mathit{ldtabl}\ge 6$.
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$F$statistic and column 5 the significance level of the
F$F$statistic. Row 1 is for replicates, row 2 for rows, row 3 for columns, row 4 for treatments (if
nt > 1${\mathbf{nt}}>1$), row 5 for residual and row 6 for total. Mean squares are computed for all but the total row,
F$F$statistics and significance are computed for treatments, replicates, rows and columns. Any unfilled cells are set to zero.
 4:
c(ldc,nt) – double array
ldc ≥ nt$\mathit{ldc}\ge {\mathbf{nt}}$.
The upper triangular part of
c contains the variancecovariance matrix of the treatment effects, the strictly lower triangular part contains the standard errors of the difference between two treatment effects (means),
i.e.,
c(i,j)${\mathbf{c}}\left(\mathit{i},\mathit{j}\right)$ contains the covariance of treatment
i$\mathit{i}$ and
j$\mathit{j}$ if
j ≥ i$\mathit{j}\ge \mathit{i}$ and the standard error of the difference between treatment
i$\mathit{i}$ and
j$\mathit{j}$ if
j < i$\mathit{j}<\mathit{i}$, for
i = 1,2, … ,t$\mathit{i}=1,2,\dots ,t$ and
j = 1,2, … ,t$\mathit{j}=1,2,\dots ,t$.
 5:
irep(nt) – int64int32nag_int array
If
nt > 1${\mathbf{nt}}>1$, the treatment replications,
R_{ll}${R}_{\mathit{l}\mathit{l}}$, for
l = 1,2, … ,nt$\mathit{l}=1,2,\dots ,{\mathbf{nt}}$. Otherwise
irep is not referenced.
 6:
rpmean(nrep) – double array
If
nrep > 1${\mathbf{nrep}}>1$,
rpmean(i)${\mathbf{rpmean}}\left(\mathit{i}\right)$ contains the mean for the
i$\mathit{i}$th replicate,
μ̂ + β̂_{i}$\hat{\mu}+{\hat{\beta}}_{\mathit{i}}$, for
i = 1,2, … ,b$\mathit{i}=1,2,\dots ,b$. Otherwise
rpmean is not referenced.
 7:
rmean(nrep × nrow${\mathbf{nrep}}\times {\mathbf{nrow}}$) – double array
rmean(j)${\mathbf{rmean}}\left(\mathit{j}\right)$ contains the mean for the
j$\mathit{j}$th row,
μ̂ + ρ̂_{i}$\hat{\mu}+{\hat{\rho}}_{i}$, for
j = 1,2, … ,r$\mathit{j}=1,2,\dots ,r$.
 8:
cmean(nrep × ncol${\mathbf{nrep}}\times {\mathbf{ncol}}$) – double array
cmean(k)${\mathbf{cmean}}\left(\mathit{k}\right)$ contains the mean for the
k$\mathit{k}$th column,
μ̂ + γ̂_{k}$\hat{\mu}+{\hat{\gamma}}_{\mathit{k}}$, for
k = 1,2, … ,c$\mathit{k}=1,2,\dots ,c$.
 9:
r(nrep × nrow × ncol${\mathbf{nrep}}\times {\mathbf{nrow}}\times {\mathbf{ncol}}$) – double array
The residuals,
r_{i}${r}_{\mathit{i}}$, for i = 1,2, … ,n$\mathit{i}=1,2,\dots ,n$.
 10:
ef(nt) – double array
If
nt ≥ 2${\mathbf{nt}}\ge 2$, the canonical efficiency factors. Otherwise
ef is not referenced.
 11:
ifail – int64int32nag_int scalar
ifail = 0${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see
[Error Indicators and Warnings]).
Error Indicators and Warnings
Note: nag_anova_rowcol (g04bc) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:
Cases prefixed with W are classified as warnings and
do not generate an error of type NAG:error_n. See nag_issue_warnings.
 ifail = 1${\mathbf{ifail}}=1$
On entry,  nrep < 1${\mathbf{nrep}}<1$, 
or  nrow < 2${\mathbf{nrow}}<2$, 
or  ncol < 2${\mathbf{ncol}}<2$, 
or  nt < 1${\mathbf{nt}}<1$, 
or  ldtabl < 6$\mathit{ldtabl}<6$, 
or  ldc < nt$\mathit{ldc}<{\mathbf{nt}}$, 
or  tol < 0.0${\mathbf{tol}}<0.0$, 
or  irdf < 0${\mathbf{irdf}}<0$. 
 ifail = 2${\mathbf{ifail}}=2$
On entry,  it(i) < 1${\mathbf{it}}\left(i\right)<1$ or it(i) > nt${\mathbf{it}}\left(i\right)>{\mathbf{nt}}$ for some i$i$ when nt ≥ 2${\mathbf{nt}}\ge 2$, 
or  no value of it = j${\mathbf{it}}=j$ for some j = 1,2, … ,nt$j=1,2,\dots ,{\mathbf{nt}}$, when nt ≥ 2${\mathbf{nt}}\ge 2$. 
 ifail = 3${\mathbf{ifail}}=3$
On entry,  the values of y are constant. 
 ifail = 4${\mathbf{ifail}}=4$
A computed standard error is zero due to rounding errors, or the eigenvalue computation failed to converge. Both are unlikely error exits.
 W ifail = 5${\mathbf{ifail}}=5$
The treatments are totally confounded with replicates, rows and columns, so the treatment sum of squares and degrees of freedom are zero. The analysis of variance table is not computed, except for replicate, row, column and total sums of squares and degrees of freedom.
 W ifail = 6${\mathbf{ifail}}=6$
The residual degrees of freedom or the residual sum of squares are 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$s$ or
s^{2}${s}^{2}$.
 W ifail = 7${\mathbf{ifail}}=7$
The design is disconnected, the standard errors may not be valid. The design may have a nested structure.
Accuracy
The algorithm used in
nag_anova_rowcol (g04bc), described in
Section [Description], achieves greater accuracy than the traditional algorithms based on the subtraction of sums of squares.
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,
nag_correg_linregs_noconst (g02cb) or
nag_correg_linregm_fit (g02da). The results from those functions 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.
Example
Open in the MATLAB editor:
nag_anova_rowcol_example
function nag_anova_rowcol_example
nrep = int64(1);
nrow = int64(5);
ncol = int64(5);
y = [6.67;
7.15;
8.29;
8.95;
9.62;
5.4;
4.77;
5.4;
7.54;
6.93;
7.32;
8.53;
8.5;
9.99;
9.68;
4.92;
5;
7.29;
7.85;
7.08;
4.88;
6.16;
7.83;
5.38;
8.51];
nt = int64(5);
it = [int64(5);4;1;3;2;2;5;4;1;3;3;2;5;4;1;1;3;2;5;4;4;1;3;2;5];
tol = 1e05;
irdf = int64(0);
[gmean, tmean, table, c, irep, rpmean, rmean, cmean, r, ef, ifail] = ...
nag_anova_rowcol(nrep, nrow, ncol, y, nt, it, tol, irdf)
gmean =
7.1856
tmean =
7.3180
7.2440
7.2060
6.9000
7.2600
table =
0 0 0 0 0
4.0000 29.4231 7.3558 9.0266 0.0013
4.0000 22.9950 5.7487 7.0545 0.0037
4.0000 0.5423 0.1356 0.1664 0.9514
12.0000 9.7788 0.8149 0 0
24.0000 62.7392 0 0 0
c =
0.1304 0.0326 0.0326 0.0326 0.0326
0.5709 0.1304 0.0326 0.0326 0.0326
0.5709 0.5709 0.1304 0.0326 0.0326
0.5709 0.5709 0.5709 0.1304 0.0326
0.5709 0.5709 0.5709 0.5709 0.1304
irep =
5
5
5
5
5
rpmean =
0
rmean =
8.1360
6.0080
8.8040
6.4280
6.5520
cmean =
5.8380
6.3220
7.4620
7.9420
8.3640
r =
0.1928
0.1632
0.2548
0.0372
0.2472
0.6812
0.4488
0.5988
0.6432
0.2768
0.1568
0.5312
0.6548
0.7152
0.4348
0.2928
0.5848
0.5272
0.5912
0.2408
0.0388
0.3392
0.9812
1.9868
0.7052
ef =
1
1
1
1
1
ifail =
0
Open in the MATLAB editor:
g04bc_example
function g04bc_example
nrep = int64(1);
nrow = int64(5);
ncol = int64(5);
y = [6.67;
7.15;
8.29;
8.95;
9.62;
5.4;
4.77;
5.4;
7.54;
6.93;
7.32;
8.53;
8.5;
9.99;
9.68;
4.92;
5;
7.29;
7.85;
7.08;
4.88;
6.16;
7.83;
5.38;
8.51];
nt = int64(5);
it = [int64(5);4;1;3;2;2;5;4;1;3;3;2;5;4;1;1;3;2;5;4;4;1;3;2;5];
tol = 1e05;
irdf = int64(0);
[gmean, tmean, table, c, irep, rpmean, rmean, cmean, r, ef, ifail] = ...
g04bc(nrep, nrow, ncol, y, nt, it, tol, irdf)
gmean =
7.1856
tmean =
7.3180
7.2440
7.2060
6.9000
7.2600
table =
0 0 0 0 0
4.0000 29.4231 7.3558 9.0266 0.0013
4.0000 22.9950 5.7487 7.0545 0.0037
4.0000 0.5423 0.1356 0.1664 0.9514
12.0000 9.7788 0.8149 0 0
24.0000 62.7392 0 0 0
c =
0.1304 0.0326 0.0326 0.0326 0.0326
0.5709 0.1304 0.0326 0.0326 0.0326
0.5709 0.5709 0.1304 0.0326 0.0326
0.5709 0.5709 0.5709 0.1304 0.0326
0.5709 0.5709 0.5709 0.5709 0.1304
irep =
5
5
5
5
5
rpmean =
0
rmean =
8.1360
6.0080
8.8040
6.4280
6.5520
cmean =
5.8380
6.3220
7.4620
7.9420
8.3640
r =
0.1928
0.1632
0.2548
0.0372
0.2472
0.6812
0.4488
0.5988
0.6432
0.2768
0.1568
0.5312
0.6548
0.7152
0.4348
0.2928
0.5848
0.5272
0.5912
0.2408
0.0388
0.3392
0.9812
1.9868
0.7052
ef =
1
1
1
1
1
ifail =
0
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013