NAG Library Routine Document
G04BBF
1 Purpose
G04BBF computes the analysis of variance and treatment means and standard errors for a randomized block or completely randomized design.
2 Specification
SUBROUTINE G04BBF ( 
N, Y, IBLOCK, NT, IT, GMEAN, BMEAN, TMEAN, TABL, LDTABL, C, LDC, IREP, R, EF, TOL, IRDF, WK, IFAIL) 
INTEGER 
N, IBLOCK, NT, IT(*), LDTABL, LDC, IREP(NT), IRDF, IFAIL 
REAL (KIND=nag_wp) 
Y(N), GMEAN, BMEAN(abs(IBLOCK)), TMEAN(NT), TABL(LDTABL,5), C(LDC,NT), R(N), EF(NT), TOL, WK(3*NT) 

3 Description
In a completely randomized design, experimental material is divided into a number of units, or plots, to which a treatment can be applied. In a randomized block design the units are grouped into blocks so that the variation within blocks is less than the variation between blocks. If every treatment is applied to one plot in each block it is a complete block design. If there are fewer plots per block than treatments then the design will be an incomplete block design and may be balanced or partially balanced.
For a completely randomized design, with
$t$ treatments and
${n}_{t}$ plots per treatment, the linear model is
where
${y}_{ij}$ is the
$i$th observation for the
$j$th treatment,
$\mu $ is the overall mean,
${\tau}_{j}$ is the effect of the
$j$th treatment and
${e}_{ij}$ is the random error term. For a randomized block design, with
$t$ treatments and
$b$ blocks of
$k$ plots, the linear model is
where
${\beta}_{i}$ is the effect of the
$i$th block and the
$ij\left(l\right)$ notation indicates that the
$l$th treatment is applied to the
$j$th plot in the
$i$th block.
The completely randomized design gives rise to a oneway analysis of variance. The treatments do not have to be equally replicated, i.e., do not have to occur the same number of times. First the overall mean, $\hat{\mu}$, is computed and subtracted from the observations to give ${y}_{ij}^{\prime}={y}_{ij}\hat{\mu}$. The estimated treatment effects, ${\hat{\tau}}_{j}$ are then computed as the treatment means of the mean adjusted observations, ${y}_{ij}^{\prime}$, and the treatment sum of squares can be computed from the sum of squares of the treatment totals of the ${y}_{ij}^{\prime}$ divided by the number of observations per treatment total, ${n}_{j}$. The final residuals are computed as ${r}_{ij}={y}_{ij}^{\prime}{\hat{\tau}}_{j}$, and, from the residuals, the residual sum of squares is calculated.
For the randomized block design the mean is computed and subtracted from the observations to give ${y}_{ij\left(l\right)}^{\prime}={y}_{ij\left(l\right)}\hat{\mu}$. The estimated block effects, ignoring treatment effects, ${\hat{\beta}}_{i}$, are then computed using the block means of the ${y}_{ij\left(l\right)}^{\prime}$ and the unadjusted sum of squares computed as the sum of squared block totals for the ${y}_{ij\left(l\right)}^{\prime}$ divided by number of plots per block, $k$. The block adjusted observations are then computed as ${y}_{ij\left(l\right)}^{\prime \prime}={{y}_{i}^{\prime}j}_{\left(l\right)}={\hat{\beta}}_{i}$. In the case of the complete block design, with the same replication for each treatment within each block, the blocks and treatments are orthogonal, and so the treatment effects are estimated as the treatment means of the block adjusted observations, ${y}_{ij\left(l\right)}^{\prime \prime}$. The treatment sum of squares is computed as the sum of squared treatment totals of the ${y}_{ij\left(l\right)}^{\prime \prime}$ divided by the number of replicates to the treatments, $r=bk/t$. Finally the residuals, and hence the residual sum of squares, are given by ${r}_{ij\left(l\right)}={y}_{ij\left(l\right)}^{\prime \prime}{\hat{\tau}}_{l}$.
For a design without the same replication for each treatment within each block the treatments and the blocks will not be orthogonal, so the treatments adjusted for blocks need to be computed. The adjusted treatment effects are found as the solution to the equations
where
$q$ is the vector of the treatment totals for block adjusted observations,
${y}_{ij\left(l\right)}^{\prime \prime}$,
$R$ is a diagonal matrix with
${R}_{ll}$ equal to the number of times the
$l$th treatment is replicated, and
${\mathbf{N}}$ is the
$t$ by
$b$ incidence matrix, with
${N}_{lj}$ equal to the number of times treatment
$l$ occurs in block
$j$. The solution to the equations can be written as
where
$\Omega $ is a generalized inverse of
$\left(RN{N}^{\mathrm{T}}/k\right)$. The solution is found from the eigenvalue decomposition of
$\left(RN{N}^{\mathrm{T}}/k\right)$. The residuals are first calculated by subtracting the estimated treatment effects from the block adjusted observations to give
${r}_{ij\left(l\right)}^{\prime}={y}_{ij\left(l\right)}^{\prime \prime}{\hat{\tau}}_{l}$. However, since only the unadjusted block effects have been removed and blocks and treatments are not orthogonal, the block means of the
${r}_{ij\left(l\right)}^{\prime}$ have to be subtracted to give the correct residuals,
${r}_{ij\left(l\right)}$ and residual sum of squares.
The mean squares are computed as the sum of squares divided by the degrees of freedom. The degrees of freedom for the unadjusted blocks is
$b1$, for the completely randomized and the complete block designs the degrees of freedom for the treatments is
$t1$. In the general case the degrees of freedom for treatments is the rank of the matrix
$\Omega $. The
$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 the completely randomized or the complete block designs, are given by:
where
${s}^{2}$ is the residual mean square and
${n}_{j}={n}_{j*}=b$ in the complete block design. 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.
In the complete block design all the information on the treatment effects is given by the within block analysis. In other designs there may be a loss of information due to the nonorthogonality of treatments and blocks. The efficiency of the within block analysis in these cases is given by the (canonical) efficiency factors, these are the nonzero eigenvalues of the matrix
$\left(RN{N}^{\mathrm{T}}/k\right)$, 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 some treatments can only be compared using a between block analysis.
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 Parameters
 1: $\mathrm{N}$ – INTEGERInput

On entry: the number of observations.
Constraint:
${\mathbf{N}}\ge 2$ and if
$\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\ge 2$,
N must be a multiple of
$\mathrm{abs}\left({\mathbf{IBLOCK}}\right)$.
 2: $\mathrm{Y}\left({\mathbf{N}}\right)$ – REAL (KIND=nag_wp) arrayInput

On entry: the observations in the order as described by
IBLOCK and
NT.
 3: $\mathrm{IBLOCK}$ – INTEGERInput

On entry: indicates the block structure.
 $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\le 1$
 There are no blocks, i.e., it is a completely randomized design.
 ${\mathbf{IBLOCK}}\ge 2$
 There are IBLOCK blocks and the data should be input by blocks, i.e., Y must contain the observations for block $1$ followed by the observations for block $2$, etc.
 ${\mathbf{IBLOCK}}\le 2$
 There are $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)$ blocks and the data is input in parallel with respect to blocks, i.e., ${\mathbf{Y}}\left(1\right)$ must contain the first observation for block $1$, ${\mathbf{Y}}\left(2\right)$ must contain the first observation for block $2\cdots {\mathbf{Y}}\left(\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\right)$ must contain the first observation for block $\mathrm{abs}\left({\mathbf{IBLOCK}}\right),{\mathbf{Y}}\left(\mathrm{abs}\left({\mathbf{IBLOCK}}+1\right)\right)$ must contain the second observation for block $1$, etc.
Constraint:
${\mathbf{IBLOCK}}=1$, $2$ or $2$.
 4: $\mathrm{NT}$ – INTEGERInput

On entry: the number of treatments. If only blocks are required in the analysis then set ${\mathbf{NT}}=1$.
Constraints:
 if $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\ge 2$, ${\mathbf{NT}}\ge 1$;
 otherwise ${\mathbf{NT}}\ge 2$.
 5: $\mathrm{IT}\left(*\right)$ – INTEGER arrayInput

Note: the dimension of the array
IT
must be at least
${\mathbf{N}}$ if
${\mathbf{NT}}\ge 2$, and at least
$1$ otherwise.
On entry:
${\mathbf{IT}}\left(\mathit{i}\right)$ indicates which of the
NT treatments plot
$\mathit{i}$ received, for
$\mathit{i}=1,2,\dots ,{\mathbf{N}}$.
If
${\mathbf{NT}}=1$,
IT is not referenced.
Constraint:
$1\le {\mathbf{IT}}\left(\mathit{i}\right)\le {\mathbf{NT}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{N}}$.
 6: $\mathrm{GMEAN}$ – REAL (KIND=nag_wp)Output

On exit: the grand mean, $\hat{\mu}$.
 7: $\mathrm{BMEAN}\left(\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\right)$ – REAL (KIND=nag_wp) arrayOutput

On exit: if $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\ge 2$,
${\mathbf{BMEAN}}\left(\mathit{j}\right)$ contains the mean for the $\mathit{j}$th block, ${\hat{\beta}}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,b$.
 8: $\mathrm{TMEAN}\left({\mathbf{NT}}\right)$ – REAL (KIND=nag_wp) arrayOutput

On exit: if ${\mathbf{NT}}\ge 2$,
${\mathbf{TMEAN}}\left(\mathit{l}\right)$ contains the (adjusted) mean for the $\mathit{l}$th treatment, ${\hat{\mu}}^{*}+{\hat{\tau}}_{\mathit{l}}$, for $\mathit{l}=1,2,\dots ,t$, where ${\hat{\mu}}^{*}$ is the mean of the treatment adjusted observations, ${y}_{ij\left(l\right)}{\hat{\tau}}_{l}$.
 9: $\mathrm{TABL}\left({\mathbf{LDTABL}},5\right)$ – REAL (KIND=nag_wp) arrayOutput

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 Blocks, row 2 for Treatments, row 3 for Residual and row 4 for Total. Mean squares are computed for all but the Total row; $F$statistics and significance are computed for Treatments and Blocks, if present. Any unfilled cells are set to zero.
 10: $\mathrm{LDTABL}$ – INTEGERInput

On entry: the first dimension of the array
TABL as declared in the (sub)program from which G04BBF is called.
Constraint:
${\mathbf{LDTABL}}\ge 4$.
 11: $\mathrm{C}\left({\mathbf{LDC}},{\mathbf{NT}}\right)$ – REAL (KIND=nag_wp) arrayOutput

On exit: if
${\mathbf{NT}}\ge 2$, 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.,
${\mathbf{C}}\left(\mathit{i},\mathit{j}\right)$ contains the covariance of treatment
$\mathit{i}$ and
$\mathit{j}$ if
$\mathit{j}\ge \mathit{i}$ and the standard error of the difference between treatment
$\mathit{i}$ and
$\mathit{j}$ if
$\mathit{j}<\mathit{i}$, for
$\mathit{i}=1,2,\dots ,t$ and
$\mathit{j}=1,2,\dots ,t$.
 12: $\mathrm{LDC}$ – INTEGERInput

On entry: the first dimension of the array
C as declared in the (sub)program from which G04BBF is called.
Constraint:
${\mathbf{LDC}}\ge {\mathbf{NT}}$.
 13: $\mathrm{IREP}\left({\mathbf{NT}}\right)$ – INTEGER arrayOutput

On exit: if ${\mathbf{NT}}\ge 2$, the treatment replications,
${R}_{\mathit{l}\mathit{l}}$, for $\mathit{l}=1,2,\dots ,{\mathbf{NT}}$.
 14: $\mathrm{R}\left({\mathbf{N}}\right)$ – REAL (KIND=nag_wp) arrayOutput

On exit: the residuals,
${r}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{\mathbf{N}}$.
 15: $\mathrm{EF}\left({\mathbf{NT}}\right)$ – REAL (KIND=nag_wp) arrayOutput

On exit: if ${\mathbf{NT}}\ge 2$, the canonical efficiency factors.
 16: $\mathrm{TOL}$ – REAL (KIND=nag_wp)Input

On entry: the tolerance value used to check for zero eigenvalues of the matrix $\Omega $. If ${\mathbf{TOL}}=0.0$ a default value of ${10}^{5}$ is used.
Constraint:
${\mathbf{TOL}}\ge 0.0$.
 17: $\mathrm{IRDF}$ – INTEGERInput

On entry: an adjustment to the degrees of freedom for the residual and total.
 ${\mathbf{IRDF}}\ge 1$
 The degrees of freedom for the total is set to ${\mathbf{N}}{\mathbf{IRDF}}$ and the residual degrees of freedom adjusted accordingly.
 ${\mathbf{IRDF}}=0$
 The total degrees of freedom for the total is set to ${\mathbf{N}}1$, as usual.
Constraint:
${\mathbf{IRDF}}\ge 0$.
 18: $\mathrm{WK}\left(3\times {\mathbf{NT}}\right)$ – REAL (KIND=nag_wp) arrayWorkspace

 19: $\mathrm{IFAIL}$ – INTEGERInput/Output

On entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, because for this routine the values of the output parameters may be useful even if
${\mathbf{IFAIL}}\ne {\mathbf{0}}$ on exit, the recommended value is
$1$.
When the value $\mathbf{1}\text{ or}1$ is used it is essential to test the value of IFAIL on exit.
On exit:
${\mathbf{IFAIL}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6 Error Indicators and Warnings
If on entry
${\mathbf{IFAIL}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Note: G04BBF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$

On entry,  ${\mathbf{N}}<2$, 
or  ${\mathbf{NT}}\le 0$, 
or  ${\mathbf{NT}}=1$ and $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\le 1$, 
or  ${\mathbf{LDTABL}}<4$, 
or  ${\mathbf{LDC}}<{\mathbf{NT}}$, 
or  ${\mathbf{TOL}}<0.0$, 
or  ${\mathbf{IRDF}}<0$. 
 ${\mathbf{IFAIL}}=2$

On entry,  $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)\ge 2$ and N is not a multiple of $\mathrm{abs}\left({\mathbf{IBLOCK}}\right)$. 
 ${\mathbf{IFAIL}}=3$

On entry,  ${\mathbf{IT}}\left(i\right)<1$ or ${\mathbf{IT}}\left(i\right)>{\mathbf{NT}}$ for some $i$ when ${\mathbf{NT}}\ge 2$, 
or  no value of ${\mathbf{IT}}=j$ for some $j=1,2,\dots ,{\mathbf{NT}}$, when ${\mathbf{NT}}\ge 2$. 
 ${\mathbf{IFAIL}}=4$

On entry,  the values of Y are constant. 
 ${\mathbf{IFAIL}}=5$

A computed standard error is zero due to rounding errors, or the eigenvalue computation failed to converge. Both are unlikely error exits.
 ${\mathbf{IFAIL}}=6$

The treatments are totally confounded with blocks, so 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.
 ${\mathbf{IFAIL}}=7$

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$ or
${s}^{2}$.
 ${\mathbf{IFAIL}}=8$

The design is disconnected; the standard errors may not be valid. The design may be nested.
 ${\mathbf{IFAIL}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.8 in the Essential Introduction for further information.
 ${\mathbf{IFAIL}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.7 in the Essential Introduction for further information.
 ${\mathbf{IFAIL}}=999$
Dynamic memory allocation failed.
See
Section 3.6 in the Essential Introduction for further information.
7 Accuracy
The algorithm used by G04BBF, described in
Section 3, achieves greater accuracy than the traditional algorithms based on the subtraction of sums of squares.
8 Parallelism and Performance
G04BBF is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
G04BBF 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 implementationspecific information.
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
$\hat{\mu}$. 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 designs such as Graeco–Latin squares one or more of the blocking factors has to be removed in a preliminary analysis before the final analysis using calls to G04BBF or
G04BCF. The residuals from the preliminary analysis are then input to G04BBF. 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.
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 blocks only.
10 Example
The data, given by
John and Quenouille (1977), are for a balanced incomplete block design with
$10$ blocks and
$6$ treatments and with
$3$ plots per block. The observations are the degree of pain experienced and the treatments are penicillin of different potency. The data is input and the analysis of variance table and treatment means are printed.
10.1 Program Text
Program Text (g04bbfe.f90)
10.2 Program Data
Program Data (g04bbfe.d)
10.3 Program Results
Program Results (g04bbfe.r)