Version:	3.2.30
Date:	2024-12-01
Title:	Functions Useful in the Design and ANOVA of Experiments
Depends:	R (≥ 3.5.0), ggplot2
Imports:	ggpubr, graphics, methods, plyr, stats, tryCatchLog
Suggests:	testthat, vdiffr, R.rsp
VignetteBuilder:	R.rsp
Description:	The content falls into the following groupings: (i) Data, (ii) Factor manipulation functions, (iii) Design functions, (iv) ANOVA functions, (v) Matrix functions, (vi) Projector and canonical efficiency functions, and (vii) Miscellaneous functions. There is a vignette describing how to use the design functions for randomizing and assessing designs available as a vignette called 'DesignNotes'. The ANOVA functions facilitate the extraction of information when the 'Error' function has been used in the call to 'aov'. The package 'dae' can also be installed from http://chris.brien.name/rpackages/.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	http://chris.brien.name
BugReports:	https://github.com/briencj/dae/issues
RoxygenNote:	5.0.1
NeedsCompilation:	no
Packaged:	2024-12-01 10:50:34 UTC; briencj
Author:	Chris Brien [aut, cre]
Maintainer:	Chris Brien <chris.brien@adelaide.edu.au>
Repository:	CRAN
Date/Publication:	2024-12-01 11:20:03 UTC

Functions Useful in the Design and ANOVA of Experiments

Description

The content falls into the following groupings: (i) Data, (ii) Factor manipulation functions, (iii) Design functions, (iv) ANOVA functions, (v) Matrix functions, (vi) Projector and canonical efficiency functions, and (vii) Miscellaneous functions. There is a vignette describing how to use the design functions for randomizing and assessing designs available as a vignette called 'DesignNotes'. The ANOVA functions facilitate the extraction of information when the 'Error' function has been used in the call to 'aov'. The package 'dae' can also be installed from <http://chris.brien.name/rpackages/>.

Version: 3.2.30

Date: 2024-12-01

Index

Author(s)

Chris Brien [aut, cre] (<https://orcid.org/0000-0003-0581-1817>)

Maintainer: Chris Brien <chris.brien@adelaide.edu.au>

Randomly generated set of values indexed by three factors

Description

This data set has randomly generated values of the response variable MOE (Measure Of Effectiveness) which is indexed by the two-level factors A, B and C.

Usage

data(ABC.Interact.dat)

Format

A data.frame containing 8 observations of 4 variables.

Source

Generated by Chris Brien

Data for a balanced incomplete block experiment

Description

The data set comes from Joshi (1987) and is the data from an experiment to investigate six varieties of wheat that employs a balanced incomplete block design (BIBD) with ten blocks, each consisting of three plots. For more details see the vignette accessed via vignette("DesignNotes", package="dae").

Usage

data(BIBDWheat.dat)

Format

A data.frame containing 30 observations of 4 variables.

Source

Joshi, D. D. (1987) Linear Estimation and Design of Experiments. Wiley Eastern, New Delhi.

A design for one of the growth cabinets in an experiment with 50 lines and 4 harvests

Description

The systematic design for a lattice square for 49 treatments consisting of four 7 x 7 squares. For more details see the vignette daeDesignNotes.pdf.

Usage

data(Cabinet1.des)

Format

A data.frame containing 160 observations of 15 variables.

Data for an experiment with rows and columns from Williams (2002)

Description

Williams (2002, p.144) provides an example of a resolved, Latinized, row-column design with four rectangles (blocks) each of six rows by ten columns. The experiment investigated differences between 60 provenances of a species of Casuarina tree, these provenances coming from 18 countries; the trees were inoculated prior to planting at two different times, time of inoculation being assigned to the four replicates so that each occurred in two replicates. At 30 months, diameter at breast height (Dbh) was measured. For more details see the vignette accessed via vignette("DesignNotes", package="dae").

Usage

data(Casuarina.dat)

Format

A data.frame containing 240 observations of 7 variables.

Source

Williams, E. R., Matheson, A. C. and Harwood, C. E. (2002) Experimental design and analysis for tree improvement. 2nd edition. CSIRO, Melbourne, Australia.

Systematic, main-unit design for an experiment to be run in a greenhouse

Description

In this main-unit design, there are 24 lanes by 11 Positions, the lanes being blocked into 6 Zones of 4 lanes. The design for the main units is for assigning 75 wheat lines, of which 73 are Nested Association Mapping (NAM) wheat lines and the other two are two check lines, Scout and Gladius. A row and column design was generated with DiGGer (Coombes, 2009). For more details see the vignette accessed via vignette("DesignNotes", package="dae").

Usage

data(Exp249.munit.des)

Format

A data.frame containing 264 observations of 3 variables.

Source

Coombes, N. E. (2009) Digger: design search tool in R. URL: http://nswdpibiom.org/austatgen/software/, (accessed June 3, 2015).

Data for a 2^4 factorial experiment

Description

The data set come from an unreplicated 2^4 factorial experiment to investigate a chemical process. The response variable is the Conversion percentage (Conv) and this is indexed by the 4 two-level factors Catal, Temp, Press and Conc, with levels “-” and “+”. The data is aranged in Yates order. Also included is the 16-level factor Runs which gives the order in which the combinations of the two-level factors were run.

Usage

data(Fac4Proc.dat)

Format

A data.frame containing 16 observations of 6 variables.

Source

Table 10.6 of Box, Hunter and Hunter (1978) Statistics for Experimenters. New York, Wiley.

A Lattice square design for 49 treatments

Description

The systematic design for a lattice square for 49 treatments consisting of four 7 x 7 squares. For more details see the vignette accessed via vignette("DesignNotes", package="dae").

Usage

data(LatticeSquare_t49.des)

Format

A data.frame containing 196 observations of 4 variables.

Source

Cochran and Cox (1957) Experimental Designs. 2nd edn Wiley, New York.

The design and data from McIntyre's (1955) two-phase experiment

Description

McIntyre (1955) reports an investigation of the effect of four light intensities on the synthesis of tobacco mosaic virus in leaves of tobacco Nicotiana tabacum var. Hickory Pryor. It is a two-phase experiment: the first phase is a treatment phase, in which the four light treatments are randomized to the tobacco leaves, and the second phase is an assay phase, in which the tobacco leaves are randomized to the half-leaves of assay plants. For more details see the vignette accessed via vignette("DesignNotes", package="dae").

Usage

data(McIntyreTMV.dat)

Format

A data.frame containing 196 observations of 4 variables.

Source

McIntyre, G. A. (1955) Design and Analysis of Two Phase Experiments. Biometrics, 11, 324–334.

Data for an experiment to investigate nitrogen response of 3 oats varieties

Description

Yates (1937) describes a split-plot experiment that investigates the effects of three varieties of oats and four levels of Nitrogen fertilizer. The varieties are assigned to the main plots using a randomized complete block design with 6 blocks and the nitrogen levels are randomly assigned to the subplots in each main plot.

The columns in the data frame are: Blocks, Wplots, Subplots, Variety, Nitrogen, xNitrogen, Yield. The column xNitrogen is a numeric version of the factor Nitrogen. The response variable is Yield.

Usage

data(Oats.dat)

Format

A data.frame containing 72 observations of 7 variables.

Author(s)

Chris Brien

Source

Yates, F. (1937). The Design and Analysis of Factorial Experiments. Imperial Bureau of Soil Science, Technical Communication, 35, 1-95.

Data for an experiment to investigate the effects of grazing patterns on pasture composition

Description

The response variable is the percentage area covered by the principal grass (Main.Grass). The design for the experiment is a split-unit design. The main units are arranged in 3 Rows x 3 Columns. Each main unit is split into 2 SubRows x 2 SubColumns.

The factor Period, with levels 3, 9 and 18 days, is assigned to the main units using a 3 x 3 Latin square. The two-level factors Spring and Summer are assigned to split-units using a criss-cross design within each main unit. The levels of each of Spring and Summer are two different grazing patterns in its season.

Usage

data(SPLGrass.dat)

Format

A data.frame containing 36 observations of 8 variables.

Source

Example 14.1 from Mead, R. (1990). The Design of Experiments: Statistical Principles for Practical Application. Cambridge, Cambridge University Press.

Data for the three-phase sensory evaluation experiment in Brien, C.J. and Payne, R.W. (1999)

Description

The data is from an experiment involved two phases. In the field phase a viticultural experiment was conducted to investigate the differences between 4 types of trellising and 2 methods of pruning. The design was a split-plot design in which the trellis types were assigned to the main plots using two adjacent Youden squares of 3 rows and 4 columns. Each main plot was split into two subplots (or halfplots) and the methods of pruning assigned at random independently to the two halfplots in each main plot. The produce of each halfplot was made into a wine so that there were 24 wines altogether.

The second phase was an evaluation phase in which the produce from the halplots was evaluated by 6 judges all of whom took part in 24 sittings. In the first 12 sittings the judges evaluated the wines made from the halfplots of one square; the final 12 sittings were to evaluate the wines from the other square. At each sitting, each judge assessed two glasses of wine from each of the halplots of one of the main plots. The main plots allocated to the judges at each sitting were determined as follows. For the allocation of rows, each occasion was subdivided into 3 intervals of 4 consecutive sittings. During each interval, each judge examined plots from one particular row, these being determined using two 3x3 Latin squares for each occasion, one for judges 1-3 and the other for judges 4-6. At each sitting judges 1-3 examined wines from one particular column and judges 4-6 examined wines from another column. The columns were randomized to the 2 sets of judges x 3 intervals x 4 sittings using duplicates of a balanced incomplete block design for v=4 and k=2 that were latinized. This balanced incomplete block design consists of three sets of 2 blocks, each set containing the 4 "treatments". For each interval, a different set of 2 blocks was taken and each block assigned to two sittings, but with the columns within the block placed in reverse order in one sitting compared to the other sitting. Thus, in each interval, a judge would evaluate a wine from each of the 4 columns.

The data.frame contains the following factors, in the order give: Occasion, Judges, Interval, Sittings, Position, Squares, Rows, Columns, Halfplot, Trellis, Method. They are followed by the simulated response variable Score.

The scores are ordered so that the factors Occasion, Judges, Interval, Sittings and Position are in standard order; the remaining factors are in randomized order.

See also the vignette accessed via vignette("DesignNotes", package="dae").

Usage

data(Sensory3Phase.dat)
data(Sensory3PhaseShort.dat)

Format

A data.frame containing 576 observations of 12 variables. There are two versions, one with shorter factor names than the other.

References

Brien, C.J. and Payne, R.W. (1999) Tiers, structure formulae and the analysis of complicated experiments. The Statistician, 48, 41-52.

Calculates the design matrix for fitting the random component of a natural cubic smoothing spline

Description

Calculates the design matrix, Z, of the random effects for a natural cubic smoothing spline as described by Verbyla et al., (1999). An initial design matrix, \bold{\Delta} \bold{\Delta}^{-1} \bold{\Delta}, based on the knot points is computed. It can then be post multiplied by the power of the tri-diagonal matrix \bold{G}_s that is proportional to the variance matrix of the random spline effects. If the power is set to 0.5 then the random spline effects based on the resulting Z matrix will be independent with variance \sigma_s^2.

Usage

Zncsspline(knot.points, Gpower = 0, print = FALSE)

Arguments

Value

A matrix containing the design matrix.

Author(s)

Chris Brien

References

Verbyla, A. P., Cullis, B. R., Kenward, M. G., and Welham, S. J. (1999). The analysis of designed experiments and longitudinal data by using smoothing splines (with discussion). Journal of the Royal Statistical Society, Series C (Applied Statistics), 48, 269-311.

See Also

mat.ncssvar.

Examples

Z <- Zncsspline(knot.points = 1:10, Gpower = 0.5)

Coerces a pstructure.object to a data.frame.

Description

Coerces a pstructure.object, which is of class pstructure, to a data.frame. One can choose whether or not to include the marginality matrix in the data.frame. The aliasing component is excluded.

Usage

## S3 method for class 'pstructure'
as.data.frame(x, row.names = NULL, optional = FALSE, ..., 
              omit.marginality = FALSE)

Arguments

Value

A data.frame with as many rows as there are non-aliased terms in the pstructure.object. The columns are df, terms, sources and, if omit.marginality is FALSE, the columns of the generated levels with columns of the marginality matrix that is stored in the marginality component of the object.

Author(s)

Chris Brien

See Also

as.data.frame.

Examples

## Generate a data.frame with 4 factors, each with three levels, in standard order
ABCD.lay <- fac.gen(list(A = 3, B = 3, C = 3, D = 3))

## create a pstructure object based on the formula ((A*B)/C)*D
ABCD.struct <- pstructure.formula(~ ((A*B)/C)*D, data =ABCD.lay)

## print the object either using the Method function or the generic function 
ABCS.dat <- as.data.frame.pstructure(ABCD.struct)
as.data.frame(ABCD.struct)

Convert a factor to a numeric vector, possibly centering or scaling the values

Description

Converts a factor to a numeric vector with approximately the numeric values of its levels. Hence, the levels of the factor must be numeric values, stored as characters. It uses the method described in factor. Use as.numeric to convert a factor to a numeric vector with integers 1, 2, ... corresponding to the positions in the list of levels. The numeric values can be centred and/or scaled. You can also use fac.recast to recode the levels to numeric values. If a numeric is supplied and both center and scale are FALSE, it is left unchanged; otherwise, it will be centred and scaled according to the settings of center and scale.

Usage

as.numfac(factor, center = FALSE, scale = FALSE)

Arguments

Details

The value of center specifies how the centring is done. If it is TRUE, the mean of the unique values of the factor are subtracted, after the factor is converted to a numeric. If center is numeric, it is subtracted from factor's converted numeric values. If center is FALSE no scaling is done.

The value of scale specifies how scaling is performed, after any centring is done. If scale is TRUE, the unique converted values of the factor are divided by (i) the standard deviaton when the values have been centred and (ii) the root-mean-square error if they have not been centred; the root-mean-square is given by \sqrt{\Sigma(x^2)/(n-1)}, where x contains the unique converted factor values and n is the number of unique values. If scale is FALSE no scaling is done.

Value

A numeric vector. An NA will be stored for any value of the factor whose level is not a number.

Author(s)

Chris Brien

See Also

as.numeric, fac.recast in package dae, factor, scale.

Examples

## set up a factor and convert it to a numeric vector
a <- factor(rep(1:6, 4))
x <- as.numfac(a)
x <- as.numfac(a, center = TRUE)
x <- as.numfac(a, center = TRUE, scale = TRUE)

This function plots a block boundary on a plot produced by designPlot.

Description

This function plots a block boundary on a plot produced by designPlot. It allows control of the starting unit, through rstart and cstart, and the number of rows (nrows) and columns (ncolumns) from the starting unit that the blocks to be plotted are to cover.

Usage

blockboundaryPlot(blockdefinition = NULL, blocksequence = FALSE, 
                  rstart= 0, cstart = 0, nrows, ncolumns, 
                  blocklinecolour = 1, blocklinewidth = 2)

Arguments

Value

no values are returned, but modifications are made to the currently active plot.

Author(s)

Chris Brien

See Also

designPlot, par, DiGGer

Examples

## Not run: 
    SPL.Lines.mat <- matrix(as.numfac(Lines), ncol=16, byrow=T)
    colnames(SPL.Lines.mat) <- 1:16
    rownames(SPL.Lines.mat) <- 1:10
    SPL.Lines.mat <- SPL.Lines.mat[10:1, 1:16]
    designPlot(SPL.Lines.mat, labels=1:10, new=TRUE,
               rtitle="Rows",ctitle="Columns", 
               chardivisor=3, rcellpropn = 1, ccellpropn=1,
               plotcellboundary = TRUE)
    #Plot Mainplot boundaries
    blockboundaryPlot(blockdefinition = cbind(4,16), rstart = 1, 
                      blocklinewidth = 3, blockcolour = "green", 
                      nrows = 9, ncolumns = 16)
    blockboundaryPlot(blockdefinition = cbind(1,4), 
                      blocklinewidth = 3, blockcolour = "green", 
                      nrows = 1, ncolumns = 16)
    blockboundaryPlot(blockdefinition = cbind(1,4), rstart= 9, nrows = 10, ncolumns = 16, 
                      blocklinewidth = 3, blockcolour = "green")
    #Plot all 4 block boundaries            
    blockboundaryPlot(blockdefinition = cbind(8,5,5,4), blocksequence=T, 
                      cstart = 1, rstart= 1, nrows = 9, ncolumns = 15, 
                      blocklinewidth = 3,blockcolour = "blue")
    blockboundaryPlot(blockdefinition = cbind(10,16), blocklinewidth=3, blockcolour="blue", 
                      nrows=10, ncolumns=16)
    #Plot border and internal block boundaries only
    blockboundaryPlot(blockdefinition = cbind(8,14), cstart = 1, rstart= 1, 
                      nrows = 9, ncolumns =  15,
                      blocklinewidth = 3, blockcolour = "blue")
    blockboundaryPlot(blockdefinition = cbind(10,16), 
                      blocklinewidth = 3, blockcolour = "blue", 
                      nrows = 10, ncolumns = 16)
## End(Not run)

Check the degrees of freedom in an object of class projector

Description

Check the degrees of freedom in an object of class "projector".

Usage

correct.degfree(object)

Arguments

Details

The degrees of freedom of the projector are obtained as its number of nonzero eigenvalues. An eigenvalue is regarded as zero if it is less than daeTolerance, which is initially set to.Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

TRUE or FALSE depending on whether the correct degrees of freedom have been stored in the object of class "projector".

Author(s)

Chris Brien

See Also

degfree, projector in package dae.

projector for further information about this class.

Examples

## set up a 2 x 2 mean operator that takes the mean of a vector of 2 values
m <- matrix(rep(0.5,4), nrow=2)

## create a projector based on the matrix m
proj.m <- new("projector", data=m)

## add its degrees of freedom
degfree(proj.m) <- 1
    
## check degrees of freedom are correct
correct.degfree(proj.m)

Deprecated Functions in Package dae

Description

These functions have been renamed and deprecated in dae.

Usage

Ameasures(...)
blockboundary.plot(...)
design.plot(...)
proj2.decomp(...)
proj2.ops(...)
projs.canon(...)
projs.structure(...)

Arguments

Author(s)

Chris Brien

The intermittent, randomly-presented, startup tips.

Description

The intermittent, randomly-presented, startup tips.

Startup tips

Need help? Enter help(package = 'dae') and click on 'User guides, package vignettes and other docs'.

Need help? The manual is in the doc subdirectory of the package's install directory.

Find out what has changed in dae: enter news(package = 'dae').

Need help to produce randomized designs? Enter vignette('DesignNotes', package = 'dae').

Need help to do the canonical analysis of a design? Enter vignette('DesignNotes', package = 'dae'). Use suppressPackageStartupMessages() to eliminate all package startup messages.

To see all the intermittent, randomly-presented, startup tips enter ?daeTips.

For versions between CRAN releases (and more) go to http://chris.brien.name/rpackages.

Author(s)

Chris Brien

Examines the relationship between the eigenvectors for two decompositions

Description

Two decompositions produced by proj2.eigen are compared by computing all pairs of crossproduct sums of eigenvectors from the two decompositions. It is most useful when the calls to proj2.eigen have the same Q1.

Usage

decomp.relate(decomp1, decomp2)

Arguments

Details

Each element of the r1 x r2 matrix is the sum of crossproducts of a pair of eigenvectors, one from each of the two decompositions. A sum is regarded as zero if it is less than daeTolerance, which is initially set to .Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

A matrix that is r1 x r2 where r1 and r2 are the numbers of efficiencies of decomp1 and decomp2, respectively. The rownames and columnnames of the matrix are the values of the efficiency factors from decomp1 and decomp2, respectively.

Author(s)

Chris Brien

See Also

proj2.eigen, proj2.combine in package dae, eigen.

Examples

## PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 
## 2nd edn Wiley, New York
PBIBD2.unit <- list(Block = 6, Unit = 4)
PBIBD2.nest <- list(Unit = "Block")
trt <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
PBIBD2.lay <- designRandomize(allocated = trt, 
                              recipient = PBIBD2.unit, 
                              nested.recipients = PBIBD2.nest)

##obtain sets of projectors
unit.struct <- pstructure(~ Block/Unit, data = PBIBD2.lay)
trt.struct <- pstructure(~ trt, data = PBIBD2.lay)

## obtain intra- and inter-block decompositions
decomp.inter <- proj2.eigen(unit.struct$Q[["Block"]], trt.struct$Q[["trt"]])
decomp.intra <- proj2.eigen(unit.struct$Q[["Unit[Block]"]], trt.struct$Q[["trt"]])

## check that intra- and inter-block decompositions are orthogonal
decomp.relate(decomp.intra, decomp.inter) 

Degrees of freedom extraction and replacement

Description

Extracts the degrees of freedom from or replaces them in an object of class "projector".

Usage

degfree(object)

degfree(object) <- value

Arguments

Details

There is no checking of the correctness of the degrees of freedom, either already stored or as a supplied integer value. This can be done using correct.degfree.

When the degrees of freedom of the projector are to be calculated, they are obtained as the number of nonzero eigenvalues. An eigenvalue is regarded as zero if it is less than daeTolerance, which is initially set to .Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

An object of class "projector" that consists of a square, summetric, idempotent matrix and degrees of freedom (rank) of the matrix.

Author(s)

Chris Brien

See Also

correct.degfree, projector in package dae.

projector for further information about this class.

Examples

## set up a 2 x 2 mean operator that takes the mean of a vector of 2 values
m <- matrix(rep(0.5,4), nrow=2)

## coerce to a projector
proj.m <- projector(m)

## extract its degrees of freedom
degfree(proj.m)

## create a projector based on the matrix m
proj.m <- new("projector", data=m)

## add its degrees of freedom and print the projector
degfree(proj.m) <- proj.m
print(proj.m)

Calculates the average variance of pairwise differences from the variance matrix for predictions

Description

Calculates the average variance of pairwise differences between, or of elementary contrasts of, predictions using the variance matrix for the predictions. The weighted average variance of pairwise differences can be computed from a vector of replications, as described by Williams and Piepho (2015). It is possible to compute either A-optimality measure for different subgroups of the predictions. If groups are specified then the A-optimality measures are calculated for the differences between predictions within each group and for those between predictions from different groups. If groupsizes are specified, but groups are not, the predictions will be sequentially broken into groups of the size specified by the elements of groupsizes. The groups can be named.

Usage

designAmeasures(Vpred, replications = NULL, groupsizes = NULL, groups = NULL)

Arguments

Details

The variance matrix of pairwise differences is calculated as v_{ii} + v_{jj} - 2 v_{ij}, where v_{ij} is the element from the ith row and jth column of Vpred. if replication is not NULL then weights are computed as r_{i} * r_{j} / \mathrm{mean}(\mathbf{r}), where \mathbf{r} is the replication vector and r_{i} and r_{j} are elements of \mathbf{r}. The (i,j) element of the variance matrix of pairwise differences is multiplied by the (i,j)th weight. Then the mean of the variances of the pairwise differences is computed for the nominated groups.

Value

A matrix containing the within and between group A-optimality measures.

Author(s)

Chris Brien

References

Smith, A. B., D. G. Butler, C. R. Cavanagh and B. R. Cullis (2015). Multi-phase variety trials using both composite and individual replicate samples: a model-based design approach. Journal of Agricultural Science, 153, 1017-1029.

Williams, E. R., and Piepho, H.-P. (2015). Optimality and contrasts in block designs with unequal treatment replication. Australian & New Zealand Journal of Statistics, 57, 203-209.

See Also

mat.Vpred, designAnatomy.

Examples

## Reduced example from Smith et al. (2015) 
## Generate two-phase design
mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3))
field.lay <- fac.gen(list(Frep = 2, Fplot = 4))
field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), 
                            levels = c("Y","W","G","M","D","E"))
start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),])
rownames(start.design) <- NULL

## Set up matrices
n <- nrow(start.design)
W <- model.matrix(~ -1+ Variety, start.design)
ng <- ncol(W)
Gg<- diag(1, ng)
Vu <- with(start.design, fac.vcmat(Mrep, 0.3) + 
                         fac.vcmat(fac.combine(list(Mrep, Mday)), 0.2) + 
                         fac.vcmat(Frep, 0.1) + 
                         fac.vcmat(fac.combine(list(Frep, Fplot)), 0.2))
R <- diag(1, n)
  
## Calculate the variance matrix of the predicted random Variety effects
Vp <- mat.Vpred(W = W, Gg = Gg, Vu = Vu, R = R)
  
## Calculate A-optimality measure
designAmeasures(Vp)
designAmeasures(Vp, groups=list(fldUndup = c(1:4), fldDup = c(5,6)))
grpsizes <- c(4,2)
names(grpsizes) <- c("fldUndup", "fldDup")
designAmeasures(Vp, groupsizes = grpsizes)
designAmeasures(Vp, groupsizes = c(4))
designAmeasures(Vp, groups=list(c(1,4),c(5,6)))

## Calculate the variance matrix of the predicted fixed Variety effects, elminating the grand mean
Vp.reduc <- mat.Vpred(W = W, Gg = 0, Vu = Vu, R = R, 
                      eliminate = projector(matrix(1, nrow = n, ncol = n)/n))
## Calculate A-optimality measure
designAmeasures(Vp.reduc)

Given the layout for a design, obtain its anatomy via the canonical analysis of its projectors to show the confounding and aliasing inherent in the design.

Description

Computes the canonical efficiency factors for the joint decomposition of two or more structures or sets of mutually orthogonally projectors (Brien and Bailey, 2009; Brien, 2017; Brien, 2019), orthogonalizing projectors in a set to those earlier in the set of projectors with which they are partially aliased. The results can be summarized in the form of a decomposition table that shows the confounding between sources from different sets. For examples of the function's use also see the vignette accessed via vignette("DesignNotes", package="dae") and for a discussion of its use see Brien, Sermarini and Demetro (2023).

Usage

designAnatomy(formulae, data, keep.order = TRUE, grandMean = FALSE, 
              orthogonalize = "hybrid", labels = "sources", 
              marginality = NULL, check.marginality = TRUE, 
              which.criteria = c("aefficiency","eefficiency","order"), 
              aliasing.print = FALSE, 
              omit.projectors = c("pcanon", "combined"), ...)

Arguments

Details

For each formula supplied in formulae, the set of projectors is obtained using pstructure; there is one projector for each term in a formula. Then projs.2canon is used to perform an analysis of the canonical relationships between two sets of projectors for the first two formulae. If there are further formulae, the relationships between its projectors and the already established decomposition is obtained using projs.2canon. The core of the analysis is the determination of eigenvalues of the products of pairs of projectors using the results of James and Wilkinson (1971). However, if the order of balance between two projection matrices is 10 or more or the James and Wilkinson (1971) methods fails to produce an idempotent matrix, equation 5.3 of Payne and Tobias (1992) is used to obtain the projection matrices for their joint decompostion.

The hybrid method is recommended for general use. However, of the three methods, eigenmethods is least likely to fail, but it does not establish the marginality between the terms. It is often needed when there is nonorthogonality between terms, such as when there are several linear covariates. It can also be more efficeint in these circumstances.

The process can be computationally expensive, particularly for a large data set (500 or more observations) and/or when many terms are to be orthogonalized.

If the error Matrix is not idempotent should occur then, especially if there are many terms, one might try using set.daeTolerance to reduce the tolerance used in determining if values are either the same or are zero; it may be necessary to lower the tolerance to as low as 0.001. Also, setting orthogonalize to eigenmethods is worth a try.

Value

A pcanon.object.

Author(s)

Chris Brien

References

Brien, C. J. (2017) Multiphase experiments in practice: A look back. Australian & New Zealand Journal of Statistics, 59, 327-352.

Brien, C. J. (2019) Multiphase experiments with at least one later laboratory phase . II. Northogonal designs. Australian & New Zealand Journal of Statistics, 61, 234-268.

Brien, C. J. and R. A. Bailey (2009). Decomposition tables for multitiered experiments. I. A chain of randomizations. The Annals of Statistics, 36, 4184-4213.

Brien, C. J., Sermarini, R. A., & Demetrio, C. G. B. (2023). Exposing the confounding in experimental designs to understand and evaluate them, and formulating linear mixed models for analyzing the data from a designed experiment. Biometrical Journal, accepted for publication.

James, A. T. and Wilkinson, G. N. (1971) Factorization of the residual operator and canonical decomposition of nonorthogonal factors in the analysis of variance. Biometrika, 58, 279-294.

Payne, R. W. and R. D. Tobias (1992). General balance, combination of information and the analysis of covariance. Scandinavian Journal of Statistics, 19, 3-23.

See Also

designRandomize, designLatinSqrSys, designPlot,
pcanon.object, p2canon.object, summary.pcanon, efficiencies.pcanon, pstructure ,
projs.2canon, proj2.efficiency, proj2.combine, proj2.eigen, efficiency.criteria, in package dae,
eigen.

projector for further information about this class.

Examples

## PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 
## 2nd edn Wiley, New York
PBIBD2.unit <- list(Block = 6, Unit = 4)
PBIBD2.nest <- list(Unit = "Block")
trt <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
PBIBD2.lay <- designRandomize(allocated = trt, 
                              recipient = PBIBD2.unit, 
                              nested.recipients = PBIBD2.nest)

##obtain combined decomposition and summarize
unit.trt.canon <- designAnatomy(formulae = list(unit=~ Block/Unit, trt=~ trt),
                                data = PBIBD2.lay)
summary(unit.trt.canon, which.criteria = c("aeff","eeff","order"))
summary(unit.trt.canon, which.criteria = c("aeff","eeff","order"), labels.swap = TRUE)

## Three-phase sensory example from Brien and Payne (1999)
## Not run: 
data(Sensory3Phase.dat)
Eval.Field.Treat.canon <- designAnatomy(formulae = list(
                              eval= ~ ((Occasions/Intervals/Sittings)*Judges)/Positions, 
                              field= ~ (Rows*(Squares/Columns))/Halfplots,
                              treats= ~ Trellis*Method),
                                        data = Sensory3Phase.dat)
summary(Eval.Field.Treat.canon, which.criteria =c("aefficiency", "order"))

## End(Not run)

Adds block boundaries to a plot produced by designGGPlot.

Description

This function adds block boundaries to a plot produced by designGGPlot. It allows control of the starting unit, through originrow and origincolumn, and the number of rows (nrows) and columns (ncolumns) from the starting unit that the blocks to be plotted are to cover.

Usage

designBlocksGGPlot(ggplot.obj, blockdefinition = NULL, blocksequence = FALSE, 
                   originrow= 0, origincolumn = 0, nrows, ncolumns, 
                   blocklinecolour = "blue", blocklinesize = 2, 
                   facetstrips.placement = "inside", 
                   printPlot = TRUE)

Arguments

Value

An object of class "ggplot", formed by adding to the input ggplot.obj and which can be plotted using print.

Author(s)

Chris Brien

Source

Brien, C.J., Harch, B.D., Correll, R.L., and Bailey, R.A. (2011) Multiphase experiments with at least one later laboratory phase. I. Orthogonal designs. Journal of Agricultural, Biological, and Environmental Statistics, 16:422-450.

See Also

designGGPlot, par, DiGGer

Examples

## Construct a randomized layout for the split-unit design described by 
## Brien et al. (2011, Section 5)
split.sys <- cbind(fac.gen(list(Months = 4, Athletes = 3, Tests = 3)),
                   fac.gen(list(Intensities = LETTERS[1:3], Surfaces = 3), 
                           times = 4))
split.lay <- designRandomize(allocated = split.sys[c("Intensities", "Surfaces")],
                             recipient = split.sys[c("Months", "Athletes", "Tests")], 
                             nested.recipients = list(Athletes = "Months", 
                                                      Tests = c("Months", "Athletes")),
                             seed = 2598)
## Plot the design
cell.colours <- c("lightblue","lightcoral","lightgoldenrod","lightgreen","lightgrey",
                  "lightpink","lightsalmon","lightcyan","lightyellow","lightseagreen")

split.lay <- within(split.lay, 
                    Treatments <- fac.combine(list(Intensities, Surfaces), 
                                              combine.levels = TRUE))
plt <- designGGPlot(split.lay, labels = "Treatments", 
                    row.factors = "Tests", column.factors = c("Months", "Athletes"),
                    colour.values = cell.colours[1:9], label.size = 6, 
                    blockdefinition = rbind(c(3,1)), blocklinecolour = "darkgreen",
                    printPlot = FALSE)
#Add Month boundaries
designBlocksGGPlot(plt, nrows = 3, ncolumns = 3, blockdefinition = rbind(c(3,3)))



#### A layout for a growth cabinet experiment that allows for edge effects
data(Cabinet1.des)
plt <- designGGPlot(Cabinet1.des, labels = "Combinations", cellalpha = 0.75,
                    title = "Lines and Harvests allocation for Cabinet 1", 
                    printPlot = FALSE)

## Plot Mainplot boundaries
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(4,16), originrow= 1 , 
                          blocklinecolour = "green", nrows = 9, ncolumns = 16, 
                          printPlot = FALSE)
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(1,4), 
                          blocklinecolour = "green", nrows = 1, ncolumns = 16, 
                          printPlot = FALSE)
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(1,4), originrow= 9, 
                          blocklinecolour = "green", nrows = 10, ncolumns = 16, 
                          printPlot = FALSE)
## Plot all 4 block boundaries            
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(8,5,5,4), blocksequence = TRUE, 
                          origincolumn = 1, originrow= 1, 
                          blocklinecolour = "blue", nrows = 9, ncolumns = 15, 
                          printPlot = FALSE)
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(10,16), 
                          blocklinecolour = "blue", nrows = 10, ncolumns = 16, 
                          printPlot = FALSE)
## Plot border and internal block boundaries only
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(8,14), origincolumn = 1, originrow= 1, 
                          blocklinecolour = "blue", nrows = 9, ncolumns = 15, 
                          printPlot = FALSE)
plt <- designBlocksGGPlot(plt, blockdefinition = cbind(10,16), 
                          blocklinecolour = "blue", nrows = 10, ncolumns = 16)

Plots labels on two-way grids of coloured cells using ggplot2 to represent an experimental design

Description

Plots the labels in a grid of cells specified by row.factors and column.factors. The cells can be coloured by the values of the column specified by column.name and can be divided into facets by specifying multiple row and or column factors.

Usage

designGGPlot(design, labels = NULL, label.size = NULL, 
             row.factors = "Rows", column.factors = "Columns", 
             scales.free = "free", facetstrips.switch = NULL, 
             facetstrips.placement = "inside", 
             cellfillcolour.column = NULL, colour.values = NULL, 
             cellalpha = 1, celllinetype = "solid", celllinesize = 0.5, 
             celllinecolour = "black", cellheight = 1, cellwidth = 1,
             reverse.x = FALSE, reverse.y = TRUE, x.axis.position = "top", 
             xlab, ylab, title, labeller = label_both, 
             title.size = 15, axis.text.size = 15, 
             blocksequence = FALSE, blockdefinition = NULL, 
             blocklinecolour = "blue", blocklinesize = 2, 
             printPlot = TRUE, ggplotFuncs = NULL, ...)

Arguments

Value

An object of class "ggplot", which can be plotted using print.

Author(s)

Chris Brien

See Also

designBlocksGGPlot, fac.combine in package dae, designPlot.

Examples

#### Plot a randomized complete block design
Treatments <- factor(rep(1:6, times = 5))
RCBD.lay <- designRandomize(allocated = Treatments,
                            recipient = list(Blocks = 5, Units = 6),
                            nested.recipients = list(Units = "Blocks"),
                            seed = 74111)
designGGPlot(RCBD.lay, labels = "Treatments", label.size = 5, 
             row.factors = "Blocks", column.factors = "Units", 
             blockdefinition = cbind(1,5))
             
## Plot without labels
designGGPlot(RCBD.lay, cellfillcolour.column = "Treatments", 
             row.factors = "Blocks", column.factors = "Units", 
             colour.values = c("lightblue","lightcoral","lightgoldenrod",
                               "lightgreen","lightgrey", "lightpink"), 
             blockdefinition = cbind(1,6))

             
#### Plot a lattice square design
data(LatticeSquare_t49.des)
designGGPlot(LatticeSquare_t49.des, labels = "Lines", label.size = 5, 
             row.factors = c("Intervals", "Runs"), column.factors = "Times", 
             blockdefinition = cbind(7,7))

Generate a systematic plan for a Latin Square design

Description

Generates a systematic plan for a Latin Square design using the method of cycling the integers 1 to the number of treatments. The start of the cycle for each row, or the first column, can be specified as a vector of integers.

Usage

designLatinSqrSys(order, start = NULL)

Arguments

Value

A numeric containing order x order integers between 1 and order such that, when the numeric is considered as a square matrix of size order, each integer occurs once and only once in each row and column of the matrix.

See Also

designRandomize, designPlot, designAnatomy in package dae.

Examples

   matrix(designLatinSqrSys(5, start = c(seq(1, 5, 2), seq(2, 5, 2))), nrow=5)
   designLatinSqrSys(3)

A graphical representation of an experimental design using labels stored in a matrix.

Description

This function uses labels, usually derived from treatment and blocking factors from an experimental design and stored in a matrix, to build a graphical representation of the matrix, highlighting the position of certain labels . It is a modified version of the function supplied with DiGGer. It includes more control over the labelling of the rows and columns of the design and allows for more flexible plotting of designs with unequal block size.

Usage

designPlot(designMatrix, labels = NULL, altlabels = NULL, plotlabels = TRUE, 
           rtitle = NULL, ctitle = NULL, 
           rlabelsreverse = FALSE, clabelsreverse = FALSE, 
           font = 1, chardivisor = 2, rchardivisor = 1, cchardivisor = 1, 
           cellfillcolour = NA, plotcellboundary = TRUE, 
           rcellpropn = 1, ccellpropn = 1, 
           blocksequence = FALSE, blockdefinition = NULL, 
           blocklinecolour = 1, blocklinewidth = 2, 
           rotate = FALSE, new = TRUE, ...)

Arguments

Value

no values are returned, but a plot is produced.

Author(s)

Chris Brien

References

Coombes, N. E. (2009). DiGGer design search tool in R. http://nswdpibiom.org/austatgen/software/

See Also

blockboundaryPlot, designPlotlabels, designLatinSqrSys, designRandomize, designAnatomy in package dae.
Also, par, polygon, DiGGer

Examples

## Not run: 
  designPlot(des.mat, labels=1:4, cellfillcolour="lightblue", new=TRUE, 
             plotcellboundary = TRUE, chardivisor=3, 
             rtitle="Lanes", ctitle="Positions", 
             rcellpropn = 1, ccellpropn=1)
  designPlot(des.mat, labels=5:87, plotlabels=TRUE, cellfillcolour="grey", new=FALSE,
             plotcellboundary = TRUE, chardivisor=3)
  designPlot(des.mat, labels=88:434, plotlabels=TRUE, cellfillcolour="lightgreen", 
             new=FALSE, plotcellboundary = TRUE, chardivisor=3, 
             blocksequence=TRUE, blockdefinition=cbind(4,10,12), 
             blocklinewidth=3, blockcolour="blue")
## End(Not run)

Plots labels on a two-way grid using ggplot2

Description

Plots the labels in a grid specified by grid.xand grid.y. The labels can be coloured by the values of the column specified by column.name.

Usage

designPlotlabels(data, labels, grid.x = "Columns", grid.y = "Rows", 
                 colour.column=NULL, colour.values=NULL, 
                 reverse.x = FALSE, reverse.y = TRUE, 
                 xlab, ylab, title, printPlot = TRUE, ggplotFuncs = NULL, ...)

Arguments

Value

An object of class "ggplot", which can be plotted using print.

Author(s)

Chris Brien

See Also

fac.combine in package dae, designPlot.

Examples

Treatments <- factor(rep(1:6, times = 5))
RCBD.lay <- designRandomize(allocated = Treatments,
                            recipient = list(Blocks = 5, Units = 6),
                            nested.recipients = list(Units = "Blocks"),
                            seed = 74111)
designPlotlabels(RCBD.lay, labels = "Treatments", 
                 grid.x = "Units", grid.y = "Blocks",
                 colour.column = "Treatments", size = 5)

Randomize allocated to recipient factors to produce a layout for an experiment

Description

A systematic design is specified by a set of allocated factors that have been assigned to a set of recipient factors. In textbook designs the allocated factors are the treatment factors and the recipient factors are the factors indexing the units. To obtain a randomized layout for a systematic design it is necessary to provide (i) the systematic arrangement of the allocated factors, (ii) a list of the recipient factors or a data.frame with their values, and (iii) the nesting of the recipient factors for the design being randomized. Given this information, the allocated factors will be randomized to the recipient factors, taking into account the nesting between the recipient factors for the design. However, allocated factors that have different values associated with those recipient factors that are in the except vector will remain unchanged from the systematic design.

Also, if allocated is NULL then a random permutation of the recipient factors is produced that is consistent with their nesting as specified by nested.recipients.

For examples of its use also see the vignette accessed via vignette("DesignNotes", package="dae") and for a discussion of its use see Brien, Sermarini and Demetro (2023).

Usage

designRandomize(allocated = NULL, recipient, nested.recipients = NULL, 
                except = NULL, seed = NULL, unit.permutation = FALSE, ...)

Arguments

Details

A systematic design is specified by the matching of the supplied allocated and recipient factors. If recipient is a list then fac.gen is used to generate a data.frame with the combinations of the levels of the recipient factors in standard order. Although, the data.frames are not combined at this stage, the systematic design is the combination, by columns, of the values of the allocated factors with the values of recipient factors in the recipient data.frame.

The method of randomization described by Bailey (1981) is used to randomize the allocated factors to the recipient factors. That is, a permutation of the recipient factors is obtained that respects the nesting for the design, but does not permute any of the factors in the except vector. A permutation is generated for all combinations of the recipient factors, except that a nested factor, specifed using the nested.recipients argument, cannot occur in a combination without its nesting factor(s). These permutations are combined into a single, units permutation that is applied to the recipient factors. Then the data.frame containing the permuted recipient factors and that containng the unpermuted allocated factors are combined columnwise, as in cbind. To produce the randomized layout, the rows of the combined data.frame are reordered so that its recipient factors are in either standard order or, if a data.frame was suppled to recipient, the same order as for the supplied data.frame.

The .Units and .Permutation vectors enable one to swap between this combined, units permutation and the randomized layout. The ith value in .Permutation gives the unit to which unit i was assigned in the randomization.

Value

A data.frame with the values for the recipient and allocated factors that specify the layout for the experiment and, if unit.permutation is TRUE, the values for .Units and .Permutation vectors.

Author(s)

Chris Brien

References

Bailey, R.A. (1981) A unified approach to design of experiments. Journal of the Royal Statistical Society, Series A, 144, 214–223.

See Also

fac.gen, designLatinSqrSys, designPlot, designAnatomy in package dae.

Examples

## Generate a randomized layout for a 4 x 4 Latin square
## (the nested.recipients argument is not needed here as none of the 
## factors are nested)
## Firstly, generate a systematic layout
LS.sys <- cbind(fac.gen(list(row = c("I","II","III","IV"), 
                             col = c(0,2,4,6))),
                treat = factor(designLatinSqrSys(4), label = LETTERS[1:4]))
## obtain randomized layout
LS.lay <- designRandomize(allocated = LS.sys["treat"], 
                          recipient = LS.sys[c("row","col")], 
                          seed = 7197132, unit.permutation = TRUE) 
LS.lay[LS.lay$.Permutation,]

## Generate a randomized layout for a replicated randomized complete 
## block design, with the block factors arranged in standard order for 
## rep then plot and then block
## Firstly, generate a systematic order such that levels of the 
## treatment factor coincide with plot
RCBD.sys <- cbind(fac.gen(list(rep = 2, plot=1:3, block = c("I","II"))),
                  tr = factor(rep(1:3, each=2, times=2)))
## obtain randomized layout
RCBD.lay <- designRandomize(allocated = RCBD.sys["tr"], 
                            recipient = RCBD.sys[c("rep", "block", "plot")], 
                            nested.recipients = list(plot = c("block","rep"), 
                                                     block="rep"), 
                            seed = 9719532, 
                            unit.permutation = TRUE)
#sort into the original standard order
RCBD.perm <- RCBD.lay[RCBD.lay$.Permutation,]
#resort into randomized order
RCBD.lay <- RCBD.perm[order(RCBD.perm$.Units),]

## Generate a layout for a split-unit experiment in which: 
## - the main-unit factor is A with 4 levels arranged in 
##   a randomized complete block design with 2 blocks;
## - the split-unit factor is B with 3 levels.
## Firstly, generate a systematic layout
SPL.sys <- cbind(fac.gen(list(block = 2, main.unit = 4, split.unit = 3)),
                 fac.gen(list(A = 4, B = 3), times = 2))
## obtain randomized layout
SPL.lay <- designRandomize(allocated = SPL.sys[c("A","B")], 
                           recipient = SPL.sys[c("block", "main.unit", "split.unit")], 
                           nested.recipients = list(main.unit = "block", 
                                                    split.unit = c("block", "main.unit")), 
                           seed=155251978)

## Generate a permutation of Seedlings within Species
seed.permute <- designRandomize(recipient = list(Species = 3, Seedlings = 4),
                                nested.recipients = list(Seedlings = "Species"),
                                seed = 75724, except = "Species", 
                                unit.permutation = TRUE)

Given the layout for a design and three structure formulae, obtain the anatomies for the (i) two-phase, (ii) first-phase, (iii) cross-phase, treatments, and (iv) combined-units designs.

Description

Uses designAnatomy to obtain the four species of designs, described by Brien (2019), that are associated with a standard two-phase design: the anatomies for the (i) two-phase, (ii) first-phase, (iii) cross-phase, treatments, and (iv) combined-units designs. (The names of the last two designs in Brien (2019) were cross-phase and second-phase designs.) For the standard two-phase design, the first-phase design is the design that allocates first-phase treatments to first-phase units. The cross-phase, treatments design allocates the first-phase treatments to the second-phase units and the combined-units design allocates the the first-phase units to the second-phase units. The two-phase design combines the other three species of designs. However, it is not mandatory that the three formula correspond to second-phase units, first-phase units and first-phase treatments, respectively, as is implied above; this is just the correspondence for a standard two-phase design. The essential requirement is that three structure formulae are supplied. For example, if there are both first- and second-phase treatments in a two-phase design, the third formula might involve the treatment factors from both phases. In this case, the default anatomy titles when printing occurs will not be correct, but can be modifed using the titles argument.

Usage

designTwophaseAnatomies(formulae, data, which.designs = "all", 
                        printAnatomies = TRUE, titles,
                        orthogonalize = "hybrid", 
                        marginality = NULL, 
                        which.criteria = c("aefficiency", "eefficiency", 
                                           "order"), ...)

Arguments

Details

To produce the anatomies, designAnatomy is called. The two-phase anatomy is based on the three formulae supplied in formulae, the first-phase anatomy uses the second and third formulae, the cross-phase, treatments anatomy derives from the first and third formulae and the combined-units anatomy is obtained with the first and second formulae.

Value

A list containing the components twophase, first, cross and combined.Each contains the pcanon.object for one of the four designs produced by designTwophaseAnatomies, unless it is NULL because the design was omitted from the which.designs argument. The returned list has an attribute titles, being a character vector of length four and containing the titles used in printing the anatomies.

Author(s)

Chris Brien

References

Brien, C. J. (2017) Multiphase experiments in practice: A look back. Australian & New Zealand Journal of Statistics, 59, 327-352.

Brien, C. J. (2019) Multiphase experiments with at least one later laboratory phase . II. Northogonal designs. Australian & New Zealand Journal of Statistics61, 234-268.

See Also

designAnatomy, pcanon.object, p2canon.object, summary.pcanon, efficiencies.pcanon,
pstructure , projs.2canon, proj2.efficiency, proj2.combine, proj2.eigen,
efficiency.criteria, in package dae, eigen.

projector for further information about this class.

Examples

  #'## Microarray example from Jarrett & Ruggiero (2008) - see Brien (2019)
  jr.lay <- fac.gen(list(Set = 7, Dye = 2, Array = 3))
  jr.lay <- within(jr.lay, 
                   { 
                     Block <- factor(rep(1:7, each=6))
                     Plant <- factor(rep(c(1,2,3,2,3,1), times=7))
                     Sample <- factor(c(rep(c(2,1,2,2,1,1, 1,2,1,1,2,2), times=3), 
                                        2,1,2,2,1,1))
                     Treat <- factor(c(1,2,4,2,4,1, 2,3,5,3,5,2, 3,4,6,4,6,3, 
                                       4,5,7,5,7,4, 5,6,1,6,1,5, 6,7,2,7,2,6, 
                                       7,1,3,1,3,7),
                                     labels=c("A","B","C","D","E","F","G"))
                   })
  
  jr.anat <- designTwophaseAnatomies(formulae = list(array = ~ (Set:Array)*Dye,
                                                     plot = ~ Block/Plant/Sample,
                                                     trt = ~ Treat),
                                     which.designs = c("first","cross"), 
                                     data = jr.lay)  

## Three-phase sensory example from Brien and Payne (1999)
## Not run: 
data(Sensory3Phase.dat)
Sensory.canon <- designTwophaseAnatomies(formulae = list(
                              eval= ~ ((Occasions/Intervals/Sittings)*Judges)/Positions, 
                              field= ~ (Rows*(Squares/Columns))/Halfplots,
                              treats= ~ Trellis*Method),
                                        data = Sensory3Phase.dat)

## End(Not run)

Computes the detectable difference for an experiment

Description

Computes the delta that is detectable for specified replication, power, alpha.

Usage

detect.diff(rm=5, df.num=1, df.denom=10, sigma=1, alpha=0.05, power=0.8, 
            tol = 0.001, print=FALSE)

Arguments

Value

A single numeric value containing the computed detectable difference.

Author(s)

Chris Brien

See Also

power.exp, no.reps in package dae.

Examples

## Compute the detectable difference for a randomized complete block design 
## with four treatments given power is 0.8 and alpha is 0.05. 
rm <- 5
detect.diff(rm = rm, df.num = 3, df.denom = 3 * (rm - 1),sigma = sqrt(20))

Extracts the canonical efficiency factors from a pcanon.object or a p2canon.object.

Description

Produces a list containing the canonical efficiency factors for the joint decomposition of two or more sets of projectors (Brien and Bailey, 2009) obtained using designAnatomy or projs.2canon.

Usage

## S3 method for class 'pcanon'
efficiencies(object, which = "adjusted", ...)
## S3 method for class 'p2canon'
efficiencies(object, which = "adjusted", ...)

Arguments

Value

For a pcanon.object, a list with a component for each component of object, except for the last component – for more information about the components see pcanon.object .

For a p2canon object, a list with a component for each element of the Q1 argument from projs.2canon. Each component is list, each its components corresponding to an element of the Q2 argument from projs.2canon

Author(s)

Chris Brien

References

Brien, C. J. and R. A. Bailey (2009). Decomposition tables for multitiered experiments. I. A chain of randomizations. The Annals of Statistics, 36, 4184 - 4213.

See Also

designAnatomy, summary.pcanon, proj2.efficiency, proj2.combine, proj2.eigen,
pstructure in package dae, eigen.

projector for further information about this class.

Examples

## PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 
## 2nd edn Wiley, New York
PBIBD2.unit <- list(Block = 6, Unit = 4)
PBIBD2.nest <- list(Unit = "Block")
trt <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
PBIBD2.lay <- designRandomize(allocated = trt, 
                              recipient = PBIBD2.unit, 
                              nested.recipients = PBIBD2.nest)

##obtain combined decomposition using designAnatomy and get the efficiencies
unit.trt.canon <- designAnatomy(list(unit=~ Block/Unit, trt=~ trt), data = PBIBD2.lay)
efficiencies.pcanon(unit.trt.canon)

##obtain the projectors for each formula using pstructure
unit.struct <- pstructure(~ Block/Unit, data = PBIBD2.lay)
trt.struct <- pstructure(~ trt, data = PBIBD2.lay)

##obtain combined decomposition projs.2canon and get the efficiencies
unit.trt.p2canon <- projs.2canon(unit.struct$Q, trt.struct$Q)
efficiencies.p2canon(unit.trt.p2canon)

Computes efficiency criteria from a set of efficiency factors

Description

Computes efficiency criteria from a set of efficiency factors.

Usage

efficiency.criteria(efficiencies)

Arguments

Details

The aefficiency criterion is the harmonic mean of the nonzero efficiency factors. The mefficiency criterion is the mean of the nonzero efficiency factors. The eefficiency criterion is the minimum of the nonzero efficiency factors. The sefficiency criterion is the variance of the nonzero efficiency factors. The xefficiency is the maximum of the efficiency factors. The order is the order of balance and is the number of unique nonzero efficiency factors. The dforthog is the number of efficiency factors that are equal to one.

Value

A list whose components are aefficiency, mefficiency, sefficiency, eefficiency, xefficiency, order and dforthog.

Author(s)

Chris Brien

See Also

proj2.efficiency, proj2.eigen, proj2.combine in package dae, eigen.

projector for further information about this class.

Examples

## PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 
## 2nd edn Wiley, New York
PBIBD2.unit <- list(Block = 6, Unit = 4)
PBIBD2.nest <- list(Unit = "Block")
trt <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
PBIBD2.lay <- designRandomize(allocated = trt, 
                              recipient = PBIBD2.unit, 
                              nested.recipients = PBIBD2.nest)

## obtain sets of projectors
unit.struct <- pstructure(~ Block/Unit, data = PBIBD2.lay)
trt.struct <- pstructure(~ trt, data = PBIBD2.lay)

## save intrablock efficiencies
eff.inter <- proj2.efficiency(unit.struct$Q[["Unit[Block]"]], trt.struct$Q[["trt"]])

## calculate efficiency criteria
efficiency.criteria(eff.inter)

Extract the elements of an array specified by the subscripts

Description

Elements of the array x corresponding to the rows of the two dimensional object subscripts are extracted. The number of columns of subscripts corresponds to the number of dimensions of x. The effect of supplying less columns in subscripts than the number of dimensions in x is the same as for "[".

Usage

elements(x, subscripts)

Arguments

Value

A vector containing the extracted elements and whose length equals the number of rows in the subscripts object.

Author(s)

Chris Brien

See Also

Extract

Examples

## Form a table of the means for all combinations of Row and Line.
## Then obtain the values corresponding to the combinations in the data frame x,
## excluding Row 3.
x <- fac.gen(list(Row = 2, Line = 4), each =2)
x$y <- rnorm(16)
RowLine.tab <- tapply(x$y, list(x$Row, x$Line), mean)
xs <- elements(RowLine.tab, subscripts=x[x$"Line" != 3, c("Row", "Line")])

Expands the values in table to a vector

Description

Expands the values in table to a vector according to the index.factors that are considered to index the table, either in standard or Yates order. The order of the values in the vector is determined by the order of the values of the index.factors.

Usage

extab(table, index.factors, order="standard")

Arguments

Value

A vector of length equal to the factors in index.factor whose values are taken from table.

Author(s)

Chris Brien

Examples

## generate a small completely randomized design with the two-level 
## factors A and B 
n <- 12
CRD.unit <- list(Unit = n)
CRD.treat <- fac.gen(list(A = 2, B = 2), each = 3)
CRD.lay <- designRandomize(allocated = CRD.treat, recipient = CRD.unit, 
                           seed = 956)

## set up a 2 x 2 table of A x B effects	
AB.tab <- c(12, -12, -12, 12)

## add a unit-length vector of expanded effects to CRD.lay
attach(CRD.lay)
CRD.lay$AB.effects <- extab(table=AB.tab, index.factors=list(A, B))

forms the ar1 correlation matrix for a (generalized) factor

Description

Form the correlation matrix for a (generalized) factor where the correlation between the levels follows an autocorrelation of order 1 (ar1) pattern.

Usage

fac.ar1mat(factor, rho)

Arguments

Details

The method is: a) form an n x n matrix of all pairwise differences in the numeric values corresponding to the observed levels of the factor by taking the difference between the following two n x n matrices are equal: 1) each row contains the numeric values corresponding to the observed levels of the factor, and 2) each column contains the numeric values corresponding to the observed levels of the factor, b) replace each element of the pairwise difference matrix with rho raised to the absolute value of the difference.

Value

An n x n matrix, where n is the length of the factor.

Author(s)

Chris Brien

See Also

fac.vcmat, fac.meanop, fac.sumop in package dae.

Examples

## set up a two-level factor and a three-level factor, both of length 12
A <- factor(rep(1:2, each=6))
B <- factor(rep(1:3, each=2, times=2))

## create a 12 x 12 ar1 matrix corrresponding to B
ar1.B <- fac.ar1mat(B, 0.6)

Combines several factors into one

Description

Combines several factors into one whose levels are the combinations of the used levels of the individual factors.

Usage

fac.combine(factors, order="standard", combine.levels=FALSE, sep=",", ...)

Arguments

Value

A factor whose levels are formed form the observed combinations of the levels of the individual factors.

Author(s)

Chris Brien

See Also

fac.uncombine, fac.split, fac.divide in package dae.

Examples

## set up two factors
A <- factor(rep(1:2, each=6))
B <- factor(rep(1:3, each=2, times=2))

## obtain six-level factor corresponding to the combinations of A and B
AB <- fac.combine(list(A,B))

Divides a factor into several separate factors

Description

Takes a factor and divides it into several separate factors as if the levels in the original combined.factor are numbered from one to its number of levels and correspond to the numbering of the levels combinations of the new factors when these are arranged in standard or Yates order.

Usage

fac.divide(combined.factor, factor.names, order="standard")

Arguments

Value

A data.frame whose columns consist of the factors listed in factor.names and whose values have been computed from the combined factor. All the factors will be of the same length.

Note

A single factor name may be supplied in the list in which case a data.frame is produced that contains the single factor computed from the numeric vector. This may be useful when calling this function from others.

Author(s)

Chris Brien

See Also

fac.split, fac.uncombine, fac.combine in package dae.

Examples

## generate a small completely randomized design for 6 treatments 
n <- 12
CRD.unit <- list(Unit = n)
treat <- factor(rep(1:4, each = 3))
CRD.lay <- designRandomize(allocated = treat, recipient = CRD.unit, seed=956)

## divide the treatments into two two-level factors A and B
CRD.facs <- fac.divide(CRD.lay$treat, factor.names = list(A = 2, B = 2))

Generate all combinations of several factors and, optionally, replicate them

Description

Generate all combinations of several factors and, optionally, replicate them.

Usage

fac.gen(generate, each=1, times=1, order="standard")

Arguments

Details

The levels of each factor are generated in a hierarchical pattern, such as standard order, where the levels of one factor are held constant while those of the adjacent factor are cycled through the complete set once. If a number is supplied instead of a name, the pattern is generated as if a factor with that number of levels had been supplied in the same position as the number. However, no levels are stored for this unamed factor.

Value

A data.frame of factors whose generated levels are those supplied in the generate list. The number of rows in the data.frame will equal the product of the numbers of levels of the supplied factors and the values of the each and times arguments.

Warning

Avoid using factor names F and T as these might be confused with FALSE and TRUE.

Author(s)

Chris Brien

See Also

fac.genfactors , fac.combine in package dae

Examples

## generate a 2^3 factorial experiment with levels - and +, and 
## in Yates order
mp <- c("-", "+")
fnames <- list(Catal = mp, Temp = mp, Press = mp, Conc = mp)
Fac4Proc.Treats <- fac.gen(generate = fnames, order="yates")

## Generate the factors A, B and D. The basic pattern has 4 repetitions
## of the levels of D for each A and B combination and 3 repetitions of 
## the pattern of the B and D combinations for each level of A. This basic 
## pattern has each combination repeated twice, and the whole of this 
## is repeated twice. It generates 864 A, B and D combinations.
gen <- list(A = 3, 3, B = c(0,100,200), 4, D = c("0","1"))
fac.gen(gen, times=2, each=2)

Generate all combinations of the levels of the supplied factors, without replication

Description

Generate all combinations of the levels of the supplied factors, without replication. This function extracts the levels from the supplied factors and uses them to generate the new factors. On the other hand, the levels must supplied in the generate argument of the function fac.gen.

Usage

fac.genfactors(factors, ...)

Arguments

Details

The levels of each factor are generated in standard order, unless order is supplied to fac.gen via the ‘...’ argument. The levels of the new factors will be in the same order as in the supplied factors.

Value

A data.frame whose columns correspond to factors in the factors list. The values in a column are the generated levels of the factor. The number of rows in the data.frame will equal the product of the numbers of levels of the supplied factors.

Author(s)

Chris Brien

See Also

fac.gen in package dae

Examples

## generate a treatments key for the Variety and Nitrogen treatments factors in Oats.dat
data(Oats.dat)
trts.key <- fac.genfactors(factors = Oats.dat[c("Variety", "Nitrogen")])
trts.key$Treatment <- factor(1:nrow(trts.key))

Match, for each combination of a set of columns in x, the row that has the same combination in table

Description

Match, for each combination of a set of columns in x, the rows that has the same combination in table. The argument multiples.allow controls what happens when there are multple matches in table of a combination in x.

Usage

fac.match(x, table, col.names, nomatch = NA_integer_, multiples.allow = FALSE)

Arguments

Value

A vector of length equal to x that gives the rows in table that match the combinations of col.names in x. The order of the rows is the same as the order of the combintions in x. The value returned if a combination is unmatched is specified in the nomatch argument.

Author(s)

Chris Brien

See Also

match

Examples

## Not run: 
#A single unmatched combination
kdata <- data.frame(Expt="D197-5", 
                    Row=8, 
                    Column=20, stringsAsFactors=FALSE)
index <- fac.match(kdata, D197.dat, c("Expt", "Row", "Column"))

# A matched and an unmatched combination
kdata <- data.frame(Expt=c("D197-5", "D197-4"), 
                    Row=c(8, 10), 
                    Column=c(20, 8), stringsAsFactors=FALSE)
index <- fac.match(kdata, D197.dat, c("Expt", "Row", "Column"))

## End(Not run)

computes the projection matrix that produces means

Description

Computes the symmetric projection matrix that produces the means corresponding to a (generalized) factor.

Usage

fac.meanop(factor)

Arguments

Details

The design matrix X for a (generalized) factor is formed with a column for each level of the (generalized) factor, this column being its indicator variable. The projection matrix is formed as X %*% (1/diag(r) %*% t(X), where r is the vector of levels replications.

A generalized factor is a factor formed from the combinations of the levels of several original factors. Generalized factors can be formed using fac.combine.

Value

A projector containing the symmetric, projection matrix and its degrees of freedom.

Author(s)

Chris Brien

See Also

fac.combine, projector, degfree, correct.degfree, fac.sumop in package dae.

projector for further information about this class.

Examples

## set up a two-level factor and a three-level factor, both of length 12
A <- factor(rep(1:2, each=6))
B <- factor(rep(1:3, each=2, times=2))

## create a generalized factor whose levels are the combinations of A and B
AB <- fac.combine(list(A,B))

## obtain the operator that computes the AB means from a vector of length 12
M.AB <- fac.meanop(AB)

Creates several factors, one for each level of the nesting factor and each of whose values are either generated within those of a level of the nesting factor or using the values of the nested factor within the levels of the nesting factor.

Description

Creates several factors, one for each level of nesting.fac and each of whose values are either (i) generated within those of the level of nesting.fac or (ii) using the values of nested.fac within the levels of the nesting.fac. For (i), all elements having the same level of nesting.fac are numbered from 1 to the number of different elements having that level. For (ii), the values of nested.fac for a level of nesting.fac are copied. In both cases, for the values of nested.fac not equal to the level of the values of nested.fac for which a nested factor is being created, the levels are set to outlevel and labelled using outlabel. A factor is not created for a level of nesting.fac with label equal to outlabel. The names of the factors are equal to the levels of nesting.fac; optionally fac.prefix is added to the beginning of the names of the factors. The function is used to split up a nested term into separate terms for each level of nesting.fac.

Usage

fac.multinested(nesting.fac, nested.fac = NULL, fac.prefix = NULL, 
                nested.levs = NA, nested.labs = NA, 
                outlevel = 0, outlabel = "rest", ...)

Arguments

Value

A data.frame containing a factor for each level of nesting.fac.

Note

The levels of nesting.fac do not have to be equally replicated.

Author(s)

Chris Brien

See Also

fac.gen, fac.nested in package dae, factor.

Examples

  lay <- data.frame(A = factor(rep(c(1:3), c(3,6,4)), labels = letters[1:3]))
  lay$B <-fac.nested(lay$A)

  #Add factors for B within each level of A
  lay2 <- cbind(lay, fac.multinested(lay$A))
  canon2 <- designAnatomy(list(~A/(a+b+c)), data = lay2)
  summary(canon2)

  #Add factors for B within each level of A, but with levels and outlabel given
  lay2 <- cbind(lay, fac.multinested(lay$A, nested.levs = seq(10,60,10), outlabel = "other"))
  
  canon2 <- designAnatomy(list(~A/(a+b+c)), data = lay2)
  summary(canon2)
  

  #Replicate the combinations of A and B three times and index them with the factor sample
  lay3 <- rbind(lay,lay,lay)
  lay3$sample <- with(lay3, fac.nested(fac.combine(list(A,B))))
  
  #Add factors for B within each level of A
  lay4 <- cbind(lay3, fac.multinested(nesting.fac = lay$A, nested.fac = lay$B))
  
  canon4 <- designAnatomy(list(~(A/(a+b+c))/sample), data = lay4)
  summary(canon4)
  

  #Add factors for sample within each combination of A and B
  lay5 <- with(lay4, cbind(lay4, 
                           fac.multinested(nesting.fac = a, fac.prefix = "a"),
                           fac.multinested(nesting.fac = b, fac.prefix = "b"),
                           fac.multinested(nesting.fac = c, fac.prefix = "c")))
  
  canon5 <- designAnatomy(list(~A/(a/(a1+a2+a3)+b/(b1+b2+b3+b4+b5+b6)+c/(c1+c2+c3))), data = lay5)
  summary(canon5)

  #Add factors for sample within each level of A
  lay6 <- cbind(lay4, 
                fac.multinested(nesting.fac = lay4$A, nested.fac = lay$sample, fac.prefix = "samp"))
  canon6 <- designAnatomy(list(~A/(a/sampa+b/sampb+c/sampc)), data = lay6)
  summary(canon6)

creates a factor, the nested factor, whose values are generated within those of the factor nesting.fac

Description

Creates a nested factor whose levels are generated within those of the factor nesting.fac. All elements of nesting.fac having the same level are numbered from 1 to the number of different elements having that level.

Usage

fac.nested(nesting.fac, nested.levs=NA, nested.labs=NA, ...)

Arguments

Value

A factor that is a character vector with class attribute "factor" and a levels attribute which determines what character strings may be included in the vector. It has a different level for of the values of the nesting.fac with the same level.

Note

The levels of nesting.fac do not have to be equally replicated.

Author(s)

Chris Brien

See Also

fac.gen, fac.multinested in package dae, factor.

Examples

## set up factor A
A <- factor(c(1, 1, 1, 2, 2))

## create nested factor
B <- fac.nested(A)

Recasts a factor by modifying the values in the factor vector and/or the levels attribute, possibly combining some levels into a single level.

Description

A factor is comprised of a vector of values and a levels attribute. This function can modify these separately or jointly. The newlevels argument recasts both the values of a factor vector and the levels attribute, using each value in the newlevels vector to replace the corresponding value in both factor vector and the levels attribute. The factor, possibly with the new levels, can have its levels attribute reordered and/or new labels associated with the levels using the levels.order and newlabels arguments.

Usage

fac.recast(factor, newlevels = NULL, levels.order = NULL, newlabels = NULL, ...)

Arguments

Value

A factor.

Author(s)

Chris Brien

See Also

fac.uselogical, as.numfac and mpone in package dae, factor, relevel.

Examples

## set up a factor with labels
Treats <- factor(rep(1:4, 4), labels=letters[1:4])
 
## recast to reduce the levels: "a" and "d" to 1 and "b" and "c" to 2, i.e. from 4 to 2 levels
A <- fac.recast(Treats, newlevels = c(1,2,2,1), labels = letters[1:2])
A <- fac.recast(Treats, newlevels = letters[c(1,2,2,1)])

#reduce the levels from 4 to 2, with re-ordering the levels vector without changing the values 
#of the new recast factor vector
A <- fac.recast(Treats, newlevels = letters[c(1,2,2,1)], levels.order = letters[2:1])  

#reassign the values in the factor vector without re-ordering the levels attribute
A <- fac.recast(Treats, newlevels = letters[4:1])  

#reassign the values in the factor vector, with re-ordering the levels attribute
A <- fac.recast(Treats, newlabels = letters[4:1])

#reorder the levels attribute with changing the values in the factor vector
A <- fac.recast(Treats, levels.order = letters[4:1])  

#reorder the values in the factor vector without changing the levels attribute
A <- fac.recast(Treats, newlevels = 4:1, newlabels = levels(Treats))

Recodes factor levels using values in a vector. The values in the vector do not have to be unique.

Description

Recodes the levels and values of a factor using each value in the newlevels vector to replace the corresponding value in the vector of levels of the factor.

This function has been superseded by fac.recast, which has extended functionality. Calls to fac.recast that use only the factor and newlevels argument will produce the same results as a call to fa.recode. fac.recode may be deprecated in future versions of dae and is being retained for now to maintain backwards compatibility.

Usage

fac.recode(factor, newlevels, ...)

Arguments

Value

A factor.

Author(s)

Chris Brien

See Also

fac.recast, fac.uselogical, as.numfac and mpone in package dae, factor, relevel.

Examples

## set up a factor with labels
Treats <- factor(rep(1:4, 4), labels=c("A","B","C","D"))
 
## recode "A" and "D" to 1 and "B" and "C" to 2
B <- fac.recode(Treats, c(1,2,2,1), labels = c("a","b"))

Splits a factor whose levels consist of several delimited strings into several factors

Description

Splits a factor, whose levels consist of strings delimited by a separator character, into several factors. It uses the function strsplit, with fixed = TRUE to split the levels.

Usage

fac.split(combined.factor, factor.names, sep=",", ...)

Arguments

Value

A data.frame containing the new factors.

Author(s)

Chris Brien

See Also

fac.divide, fac.uncombine, fac.combine in package dae and strsplit.

Examples

## Form a combined factor to split
data(Oats.dat)
tmp <- within(Oats.dat, Trts <- fac.combine(list(Variety, Nitrogen), combine.levels = TRUE))

##Variety levels sorted into alphabetical order
trts.dat <- fac.split(combined.factor = tmp$Trts, 
                      factor.names = list(Variety = NULL, Nitrogen = NULL))

##Variety levels order from Oats.dat retained
trts.dat <- fac.split(combined.factor = tmp$Trts, 
                      factor.names = list(Variety = levels(tmp$Variety), Nitrogen = NULL))

computes the summation matrix that produces sums corresponding to a (generalized) factor

Description

Computes the matrix that produces the sums corresponding to a (generalized) factor.

Usage

fac.sumop(factor)

Arguments

Details

The design matrix X for a (generalized) factor is formed with a column for each level of the (generalized) factor, this column being its indicator variable. The summation matrix is formed as X %*% t(X).

A generalized factor is a factor formed from the combinations of the levels of several original factors. Generalized factors can be formed using fac.combine.

Value

A symmetric matrix.

Author(s)

Chris Brien

See Also

fac.combine, fac.meanop in package dae.

Examples

## set up a two-level factoir and a three-level factor, both of length 12
A <- factor(rep(1:2, each=6))
B <- factor(rep(1:3, each=2, times=2))

## create a generlaized factor whose levels are the combinations of A and B
AB <- fac.combine(list(A,B))

## obtain the operator that computes the AB means from a vector of length 12
S.AB <- fac.sumop(AB)

Cleaves a single factor, each of whose levels has delimited strings, into several factors using the separated strings.

Description

Cleaves a single factor into several factors whose levels, the levels of the original factor consisting of several delimited strings that can be separated to form the levels of the new.factors. That is, it reverses the process of combining factors that fac.combine performs.

Usage

fac.uncombine(factor, new.factors, sep=",", ...)

Arguments

Value

A data.frame whose columns consist of the factors listed in new.factors and whose values have been computed from the values of the combined factor.

Author(s)

Chris Brien

See Also

fac.split, fac.combine, fac.divide in package dae and strsplit.

Examples

## set up two factors and combine them
facs <- fac.gen(list(A = letters[1:3], B = 1:2), each = 4)
facs$AB <- with(facs, fac.combine(list(A, B), combine.levels = TRUE))

## now reverse the proces and uncombine the two factors
new.facs <- fac.uncombine(factor = facs$AB, 
                          new.factors = list(A = letters[1:3], B = NULL), 
                          sep = ",")
new.facs <- fac.uncombine(factor = facs$AB, 
                          new.factors = list(A = NULL, B = NULL), 
                          sep = ",")

Forms a two-level factor from a logical object.

Description

Forms a two-level factor from a logical object. It can be used to recode a factor when the resulting factor is to have only two levels.

Usage

fac.uselogical(x, levels = c(TRUE, FALSE), labels = c("yes", "no"), ...)

Arguments

Value

A factor.

Author(s)

Chris Brien

See Also

fac.recast, as.numfac and mpone in package dae, factor, relevel.

Examples

## set up a factor with labels
Treats <- factor(rep(1:4, 4), labels=c("A","B","C","D"))
 
## recode "A" and "D" to "a" and "B" and "C" to "b"
B <- fac.uselogical(Treats %in% c("A", "D"), labels = c("a","b"))
B <- fac.uselogical(Treats %in% c("A", "D"), labels = c(-1,1))

## suppose level A in factor a is a control treatment
## set up a factor Control to discriminate between control and treated
Control <- fac.uselogical(Treats == "A")

forms the variance matrix for the variance component of a (generalized) factor

Description

Form the variance matrix for a (generalized) factor whose effects for its different levels are independently and identically distributed, with their variance given by the variance component; elements of the matrix will equal either zero or sigma2 and displays compound symmetry.

Usage

fac.vcmat(factor, sigma2)

Arguments

Details

The method is: a) form the n x n summation or relationship matrix whose elements are equal to zero except for those elements whose corresponding elements in the following two n x n matrices are equal: 1) each row contains the numeric values corresponding to the observed levels of the factor, and 2) each column contains the numeric values corresponding to the observed levels of the factor, b) multiply the summation matrix by sigma2.

Value

An n x n matrix, where n is the length of the factor.

Author(s)

Chris Brien

See Also

fac.ar1mat, fac.meanop, fac.sumop in package dae.

Examples

## set up a two-level factor and a three-level factor, both of length 12
A <- factor(rep(1:2, each=6))
B <- factor(rep(1:3, each=2, times=2))

## create a 12 x 12 ar1 matrix corrresponding to B
vc.B <- fac.vcmat(B, 2)

Extract the fitted values for a fitted model from an aovlist object

Description

Extracts the fitted values as the sum of the effects for all the fitted terms in the model, stopping at error.term if this is specified. It is a method for the generic function fitted.

Usage

## S3 method for class 'aovlist'
fitted(object, error.term=NULL, ...)

Arguments

Value

A numeric vector of fitted values.

Note

Fitted values will be the sum of effects for terms from the model, but only for terms external to any Error function. If you want effects for terms in the Error function to be included, put them both inside and outside the Error function so they are occur twice.

Author(s)

Chris Brien

See Also

fitted.errors, resid.errors, tukey.1df in package dae.

Examples

## set up data frame for randomized complete block design in Table 4.4 from 
## Box, Hunter and Hunter (2005) Statistics for Experimenters. 2nd edn 
## New York, Wiley.
RCBDPen.dat <- fac.gen(list(Blend=5, Flask=4))
RCBDPen.dat$Treat <- factor(rep(c("A","B","C","D"), times=5))
RCBDPen.dat$Yield <- c(89,88,97,94,84,77,92,79,81,87,87,
                       85,87,92,89,84,79,81,80,88)

## perform the analysis of variance
RCBDPen.aov <- aov(Yield ~ Blend + Treat + Error(Blend/Flask), RCBDPen.dat)
summary(RCBDPen.aov)

## two equivalent ways of extracting the fitted values
fit  <- fitted.aovlist(RCBDPen.aov)
fit <- fitted(RCBDPen.aov, error.term = "Blend:Flask")

Extract the fitted values for a fitted model

Description

An alias for the generic function fitted. When it is available, the method fitted.aovlist extracts the fitted values, which is provided in the dae package to cover aovlist objects.

Usage

## S3 method for class 'errors'
fitted(object, error.term=NULL, ...)

Arguments

Value

A numeric vector of fitted values.

Warning

See fitted.aovlist for specific information about fitted values when an Error function is used in the call to the aov function.

Author(s)

Chris Brien

See Also

fitted.aovlist, resid.errors, tukey.1df in package dae.

Examples

## set up data frame for randomized complete block design in Table 4.4 from 
## Box, Hunter and Hunter (2005) Statistics for Experimenters. 2nd edn 
## New York, Wiley.
RCBDPen.dat <- fac.gen(list(Blend=5, Flask=4))
RCBDPen.dat$Treat <- factor(rep(c("A","B","C","D"), times=5))
RCBDPen.dat$Yield <- c(89,88,97,94,84,77,92,79,81,87,87,
                       85,87,92,89,84,79,81,80,88)

## perform the analysis of variance
RCBDPen.aov <- aov(Yield ~ Blend + Treat + Error(Blend/Flask), RCBDPen.dat)
summary(RCBDPen.aov)

## three equivalent ways of extracting the fitted values
fit  <- fitted.aovlist(RCBDPen.aov)
fit <- fitted(RCBDPen.aov, error.term = "Blend:Flask")
fit <- fitted.errors(RCBDPen.aov, error.term = "Blend:Flask")

Gets the value of daeRNGkind for the package dae from the daeEnv environment

Description

A function that gets the character value of daeRNGkind from the daeEnv environment. The value specifies the name of the Random Number Generator to use in dae.

Usage

get.daeRNGkind()

Value

The character value of daeRNGkind.

Author(s)

Chris Brien

See Also

set.daeRNGkind.

Examples

## get daeRNGkind.
get.daeRNGkind()

Gets the value of daeTolerance for the package dae

Description

A function that gets the vector of values such that, in dae functions, values less than it are considered to be zero.

Usage

get.daeTolerance()

Value

The vector of two values for daeTolerance, one named element.tol that is used for elements of matrices and a second named element.eigen that is used for eigenvalues and quantities based on them, such as efficiency factors.

Author(s)

Chris Brien

See Also

set.daeTolerance.

Examples

## get daeTolerance.
get.daeTolerance()

Calcuates the harmonic mean.

Description

A function to calcuate the harmonic mean of a set of nonzero numbers.

Usage

harmonic.mean(x)

Arguments

Details

All the elements of x are tested as being less than daeTolerance, which is initially set to .Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

A numeric. Returns Inf if x contains a value close to zero

Examples

y <- c(seq(0.1,1,0.2))
harmonic.mean(y)

Plots an interaction plot for three factors

Description

Plots a function (the mean by default) of the response for the combinations of the three factors specified as the x.factor (plotted on the x axis of each plot), the groups.factor (plotted as separate lines in each plot) and the trace.factor (its levels are plotted in different plots). Interaction plots for more than three factors can be produced by using fac.combine to combine all but two of them into a single factor that is specified as the trace.factor.

Usage

interaction.ABC.plot(response, x.factor, groups.factor, 
       trace.factor,data, fun="mean", title="A:B:C Interaction Plot", 
       xlab, ylab, key.title, lwd=4, columns=2, ggplotFuncs = NULL, ...)

Arguments

Value

An object of class "ggplot", which can be plotted using print.

Author(s)

Chris Brien

See Also

fac.combine in package dae, interaction.plot.

Examples

## Not run: 
## plot for Example 14.1 from Mead, R. (1990). The Design of Experiments: 
## Statistical Principles for Practical Application. Cambridge, 
## Cambridge University Press.  
## use ?SPLGrass.dat for details
data(SPLGrass.dat)
interaction.ABC.plot(Main.Grass, x.factor=Period,
                     groups.factor=Spring, trace.factor=Summer,
                     data=SPLGrass.dat,
                     title="Effect of Period, Spring and Summer on Main Grass")

## plot for generated data
## use ?ABC.Interact.dat for data set details
data(ABC.Interact.dat)
## Add standard errors for plotting 
## - here data contains a single value for each combintion of A, B and C
## - need to supply name for data twice 
ABC.Interact.dat$se <- rep(c(0.5,1), each=4)
interaction.ABC.plot(MOE, A, B, C, data=ABC.Interact.dat,
                     ggplotFunc=list(geom_errorbar(data=ABC.Interact.dat, 
                                                   aes(ymax=MOE+se, ymin=MOE-se), 
                                                   width=0.2)))

## End(Not run)

Tests whether all elements are approximately zero

Description

A single-line function that tests whether all elements are zero (approximately).

Usage

is.allzero(x)

Arguments

Details

The mean of the absolute values of the elements of x is tested to determine if it is less than daeTolerance, which is initially set to .Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

A logical.

Author(s)

Chris Brien

Examples

## create a vector of 9 zeroes and a one
y <- c(rep(0,9), 1)

## check that vector is only zeroes is FALSE 
is.allzero(y)

Tests whether an object is a valid object of class projector

Description

Tests whether an object is a valid object of class "projector".

Usage

is.projector(object)

Arguments

Details

The function is.projector tests whether the object consists of a matrix that is square, symmetric and idempotent. In checking symmetry and idempotency, the equality of the matrix with either its transpose or square is tested. In this, a difference in elements is considered to be zero if it is less than daeTolerance, which is initially set to .Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

TRUE or FALSE depending on whether the object is a valid object of class "projector".

Warning

The degrees of freedom are not checked. correct.degfree can be used to check them.

Author(s)

Chris Brien

See Also

projector, correct.degfree in package dae.

projector for further information about this class.

Examples

## set up a 2 x 2 mean operator that takes the mean of a vector of 2 values
m <- matrix(rep(0.5,4), nrow=2)

## create an object of class projector
proj.m <- projector(m)

## check that it is a valid projector
is.projector(proj.m)

Extracts the marginality matrix (matrices) from a pstructure.object or a pcanon.object.

Description

Produces (i) a marginality matrix for the formula in a call to pstructure.formula or (ii) a list containing the marginlity matrices, one for each formula in the formulae argument of a call to designAnatomy.

A marginality matrix for a set of terms is a square matrix with a row and a column for each ternon-aliased term. Its elements are zeroes and ones, the entry in the ith row and jth column indicates whether or not the ith term is marginal to the jth term i.e. the column space of the ith term is a subspace of that for the jth term and so the source for the jth term will be orthogonal to that for the ith term.

Usage

## S3 method for class 'pstructure'
marginality(object, ...)
## S3 method for class 'pcanon'
marginality(object, ...)

Arguments

Value

If object is a pstructure.object then a matrix containing the marginality matrix for the terms obtained from the formuula in the call to pstructure.formula.

If object is a pcanon.object then a list with a component for each formula, each component having a marginality matrix that corresponds to one of the formulae in the call to designAnatomy. The components of the list will have the same names as the componeents of the formulae list and so will be unnamed if the components of the latter list are unnamed.

Author(s)

Chris Brien

See Also

pstructure.formula, designAnatomy, summary.pcanon, proj2.efficiency, proj2.combine, proj2.eigen,
pstructure in package dae, eigen.

projector for further information about this class.

Examples

## PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 
## 2nd edn Wiley, New York
PBIBD2.unit <- list(Block = 6, Unit = 4)
PBIBD2.nest <- list(Unit = "Block")
trt <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
PBIBD2.lay <- designRandomize(allocated = trt, 
                              recipient = PBIBD2.unit, 
                              nested.recipients = PBIBD2.nest)

##obtain pstructure.object and extract marginality matrix
unit.struct <- pstructure(~ Block/Unit, data = PBIBD2.lay)
unit.marg <- marginality(unit.struct)

##obtain combined decomposition and extract marginality matrices
unit.trt.canon <- designAnatomy(list(unit=~ Block/Unit, trt=~ trt), data = PBIBD2.lay)
marg <- marginality(unit.trt.canon)

Forms a unit matrix

Description

Form the unit or identity matrix of order order.

Usage

mat.I(order)

Arguments

Value

A square matrix whose diagonal elements are one and its off-diagonal are zero.

Author(s)

Chris Brien

See Also

mat.J, mat.ar1

Examples

    col.I <- mat.I(order=4)

Forms a square matrix of ones

Description

Form the square matrix of ones of order order.

Usage

mat.J(order)

Arguments

Value

A square matrix all of whose elements are one.

Author(s)

Chris Brien

See Also

mat.I, mat.ar1

Examples

    col.J <- mat.J(order=4)

Calculates the variances of a set of predicted effects from a mixed model

Description

For n observations, w effects to be predicted, f nuiscance fixed effects and r nuisance random effects, the variances of a set of predicted effects is calculated using the incidence matrix for the effects to be predicted and, optionally, a variance matrix of the effects, an incidence matrix for the nuisance fixed factors and covariates, the variance matrix of the nuisance random effects in the mixed model and the residual variance matrix.

This function has been superseded by mat.Vpredicts, which allows the use of both matrices and formulae.

Usage

mat.Vpred(W, Gg = 0, X = matrix(1, nrow = nrow(W), ncol = 1), Vu = 0, R, eliminate)

Arguments

Details

Firstly the information matrix is calculated as
A <- t(W) %*% Vinv %*% W + ginv(Gg) - A%*%ginv(t(X)%*%Vinv%*%X)%*%t(A), where Vinv <- ginv(Vu + R), A = t(W) %*% Vinv %*% X and ginv(B) is the unique Moore-Penrose inverse of B formed using the eigendecomposition of B.

If eliminate is set and the effects to be predicted are fixed then the reduced information matrix is calculated as A <- (I - eliminate) Vinv (I - eliminate).

Finally, the variance of the predicted effects is calculated: Vpred <- ginv(A).

Value

A w x w matrix containing the variances and covariances of the predicted effects.

Author(s)

Chris Brien

References

Smith, A. B., D. G. Butler, C. R. Cavanagh and B. R. Cullis (2015). Multi-phase variety trials using both composite and individual replicate samples: a model-based design approach. Journal of Agricultural Science, 153, 1017-1029.

See Also

designAmeasures, mat.Vpredicts.

Examples

## Reduced example from Smith et al. (2015)
## Generate two-phase design
mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3))
field.lay <- fac.gen(list(Frep = 2, Fplot = 4))
field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), 
                            levels = c("Y","W","G","M","D","E"))
start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),])
rownames(start.design) <- NULL

## Set up matrices
n <- nrow(start.design)
W <- model.matrix(~ -1+ Variety, start.design)
ng <- ncol(W)
Gg<- diag(1, ng)
Vu <- with(start.design, fac.vcmat(Mrep, 0.3) + 
                         fac.vcmat(fac.combine(list(Mrep, Mday)), 0.2) + 
                         fac.vcmat(Frep, 0.1) + 
                         fac.vcmat(fac.combine(list(Frep, Fplot)), 0.2))
R <- diag(1, n)
  
## Calculate the variance matrix of the predicted random Variety effects
Vp <- mat.Vpred(W = W, Gg = Gg, Vu = Vu, R = R)
designAmeasures(Vp)

## Calculate the variance matrix of the predicted fixed Variety effects, 
## elminating the grand mean
Vp.reduc <- mat.Vpred(W = W, Gg = 0, Vu = Vu, R = R, 
                      eliminate = projector(matrix(1, nrow = n, ncol = n)/n))
designAmeasures(Vp.reduc)

Calculates the variances of a set of predicted effects from a mixed model, based on supplied matrices or formulae.

Description

For n observations, w effects to be predicted, f nuiscance fixed effects, r nuisance random effects and n residuals, the variances of a set of predicted effects is calculated using the incidence matrix for the effects to be predicted and, optionally, a variance matrix of these effects, an incidence matrix for the nuisance fixed factors and covariates, the variance matrix of the nuisance random effects and the residual variance matrix. The matrices can be supplied directly or using formulae and a matrix specifying the variances of the nuisance random effects. The difference between mat.Vpredicts and mat.Vpred is that the former has different names for equivalent arguments and the latter does not allow for the use of formulae.

Usage

mat.Vpredicts(target, Gt = 0, fixed = ~ 1, random, G, R, design, 
              eliminate, keep.order = TRUE, result = "variance.matrix")

Arguments

Details

The mixed model for which the predictions are to be obtained is of the form \bold{Y} = \bold{X\beta} + \bold{Ww} + \bold{Zu} + \bold{e}, where \bold{W} is the incidence matrix for the target predicted effects \bold{w}, \bold{X} is the is incidence matrix for the fixed nuisance effects \bold{\beta}, \bold{Z} is the is incidence matrix for the random nuisance effects \bold{u}, \bold{e} are the residuals; the \bold{u} are assumed to have variance matrix \bold{G} so that their contribution to the variance matrix for \bold{Y} is \bold{Vu} = \bold{ZGZ}^T and \bold{e} is assumed to have variance matrix \bold{R}. If the target effects are random then the variance matrix for \bold{w} is \bold{G}_t so that their contribution to the variance matrix for \bold{Y} is \bold{WG}_t\bold{W}^T.

As described in Hooks et al. (2009, Equation 19), the information matrix is calculated as
A <- t(W) %*% Vinv %*% W + ginv(Gg) - A%*%ginv(t(X)%*%Vinv%*%X)%*%t(A), where Vinv <- ginv(Vu + R), A = t(W) %*% Vinv %*% X and ginv(B) is the unique Moore-Penrose inverse of B formed using the eigendecomposition of B.

Then, if eliminate is set and the effects to be predicted are fixed then the reduced information matrix is calculated as A <- (I - eliminate) Vinv (I - eliminate).

Finally, if result is set to variance.matrix, the variance of the predicted effects is calculated: Vpred <- ginv(A) and returned; otherwise the information matrix A is returned. The rank of the matrix to be returned is obtain via a singular value decomposition of the information matrix, it being the number of nonzero eigenvalues. An eigenvalue is regarded as zero if it is less than daeTolerance, which is initially set to.Machine$double.eps ^ 0.5 (about 1.5E-08). The function set.daeTolerance can be used to change daeTolerance.

Value

A w x w matrix containing the variances and covariances of the predicted effects or the information matrix for the effects, depending on the setting of result. The matrix has its rank as an attribute.

Author(s)

Chris Brien

References

Hooks, T., Marx, D., Kachman, S., and Pedersen, J. (2009). Optimality criteria for models with random effects. Revista Colombiana de Estadistica, 32, 17-31.

Smith, A. B., D. G. Butler, C. R. Cavanagh and B. R. Cullis (2015). Multi-phase variety trials using both composite and individual replicate samples: a model-based design approach. Journal of Agricultural Science, 153, 1017-1029.

See Also

designAmeasures, mat.random, mat.Vpred.

Examples

## Reduced example from Smith et al. (2015)
## Generate two-phase design
mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3))
field.lay <- fac.gen(list(Frep = 2, Fplot = 4))
field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), 
                            levels = c("Y","W","G","M","D","E"))
start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),])
rownames(start.design) <- NULL

## Set gammas
terms <- c("Variety", "Frep", "Frep:Fplot", "Mrep", "Mrep:Mday", "Mrep:Mday:Mord")
gammas <- c(1, 0.1, 0.2, 0.3, 0.2, 1)
names(gammas) <- terms

## Specify matrices to calculate the variance matrix of the predicted fixed Variety effects 
W <- model.matrix(~ -1 + Variety, start.design)
Vu <- with(start.design, fac.vcmat(Mrep, gammas["Mrep"]) + 
                         fac.vcmat(fac.combine(list(Mrep,Mday)), gammas["Mrep:Mday"]) + 
                         fac.vcmat(Frep, gammas["Frep"]) + 
                         fac.vcmat(fac.combine(list(Frep,Fplot)), gammas["Frep:Fplot"]))
R <- diag(1, nrow(start.design))
  
## Calculate variance matrix
Vp <- mat.Vpredicts(target = W, random=Vu, R=R, design = start.design)

## Calculate the variance matrix of the predicted random Variety effects using formulae
Vp <- mat.Vpredicts(target = ~ -1 + Variety, Gt = 1, 
                    fixed = ~ 1, 
                    random = ~ -1 + Mrep/Mday + Frep/Fplot, 
                    G = as.list(gammas[c(4,5,2,3)]), 
                    R = R, design = start.design)
designAmeasures(Vp)

## Calculate the variance matrix of the predicted fixed Variety effects, 
## elminating the grand mean
n <- nrow(start.design)
Vp.reduc <- mat.Vpredicts(target = ~ -1 + Variety, 
                          random = ~ -1 + Mrep/Mday + Frep/Fplot, 
                          G = as.list(gammas[c(4,5,2,3)]), 
                          eliminate = projector(matrix(1, nrow = n, ncol = n)/n), 
                          design = start.design)
designAmeasures(Vp.reduc)

Forms an ar1 correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the ar1 pattern. The matrix is banded and has diagonal elements equal to one and the off-diagonal element in the ith row and jth column equal to \rho^k where k = |i- j|.

Usage

mat.ar1(rho, order)

Arguments

Value

A banded correlation matrix whose elements follow an ar1 pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.exp, mat.gau, mat.banded, mat.ar2, mat.ar3, mat.sar2, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.ar1(rho=0.4, order=4)

Forms an ar2 correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the ar2 pattern. The resulting matrix is banded.

Usage

mat.ar2(ARparameters, order)

Arguments

Details

The correlations in the correlation matrix, corr say, are calculated from the autoregressive parameters, ARparameters. The values in

• the diagonal (k = 1) of corr are one;

• the first subdiagonal band (k = 2) of corr are equal to
  ARparameters[1]/(1-ARparameters[2]);

• in subsequent disgonal bands, (k = 3:order), of corr are
  ARparameters[1]*corr[k-1] + ARparameters[2]*corr[k-2].

Value

A banded correlation matrix whose elements follow an ar2 pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.exp, mat.gau, mat.banded, mat.ar1, mat.ar3, mat.sar2, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.ar2(ARparameters = c(0.4, 0.2), order = 4)

Forms an ar3 correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the ar3 pattern. The resulting matrix is banded.

Usage

mat.ar3(ARparameters, order)

Arguments

Details

The correlations in the correlation matrix, corr say, are calculated from the autoregressive parameters, ARparameters.
Let omega = 1 - ARparameters[2] - ARparameters[3] * (ARparameters[1] + ARparameters[3]). Then the values in

• the diagonal of corr (k = 1) are one;

• the first subdiagonal band (k = 2) of corr are equal to
  (ARparameters[1] + ARparameters[2]*ARparameters[3]) / omega;

• the second subdiagonal band (k = 3) of corr are equal to
  (ARparameters[1] * (ARparameters[1] + ARparameters[3]) +
ARparameters[2] * (1 - ARparameters[2])) / omega;

• the subsequent subdiagonal bands, (k = 4:order), of corr are equal to
  ARparameters[1]*corr[k-1] + ARparameters[2]*corr[k-2] + ARparameters[3]*corr[k-3].

Value

A banded correlation matrix whose elements follow an ar3 pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.banded, mat.exp, mat.gau, mat.ar1, mat.ar2, mat.sar2, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.ar3(ARparameters = c(0.4, 0.2, 0.1), order = 4)

Forms an arma correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the arma pattern. The resulting matrix is banded.

Usage

mat.arma(ARparameter, MAparameter, order)

Arguments

Details

The correlations in the correlation matrix, corr say, are calculated from the correlation parameters, ARparameters. The values in

• the diagonal (k = 1) of corr are one;

• the first subdiagonal band (k = 2) of corr are equal to
  ARparameters[1]/(1-ARparameters[2]);

• in subsequent disgonal bands, (k = 3:order), of corr are
  ARparameters[1]*corr[k-1] + ARparameters[2]*corr[k-2].

Value

A banded correlation matrix whose elements follow an arma pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.exp, mat.gau, mat.banded, mat.ar1, mat.ar3, mat.sar2, mat.ma1, mat.ma2

Examples

    corr <- mat.arma(ARparameter = 0.4, MAparameter = -0.2, order = 4)

Form a banded matrix from a vector of values

Description

Takes the first value in x and places it down the diagonal of the matrix. Takes the second value in x and places it down the first subdiagonal, both below and above the diagonal of the matrix. The third value is placed in the second subdiagonal and so on, until the bands for which there are elements in x have been filled. All other elements in the matrix will be zero.

Usage

mat.banded(x, nrow, ncol)

Arguments

Value

An nrow \times ncol matrix.

Author(s)

Chris Brien

See Also

mat.cor, mat.corg, mat.ar1, mat.ar2, mat.ar3, mat.sar2, mat.exp, mat.gau, mat.ma1, mat.ma2, mat.arma mat.I, mat.J

Examples

      m <- mat.banded(c(1,0.6,0.5), 5,5)
      m <- mat.banded(c(1,0.6,0.5), 3,4)
      m <- mat.banded(c(1,0.6,0.5), 4,3)

Forms a correlation matrix in which all correlations have the same value.

Description

Form the correlation matrix of order order in which all correlations have the same value.

Usage

mat.cor(rho, order)

Arguments

Value

A correlation matrix.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.corg, mat.banded, mat.exp, mat.gau, mat.ar1, mat.ar2, mat.sar2, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.cor(rho = 0.4, order = 3)

Forms a general correlation matrix

Description

Form the correlation matrix of order order for which all correlations potentially differ.

Usage

mat.corg(rhos, order, byrow = FALSE)

Arguments

Value

A correlation matrix.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.banded, mat.exp, mat.gau, mat.ar1, mat.ar2, mat.sar2, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.corg(rhos = c(0.4, 0.2, 0.1), order = 3)

Forms the direct product of two matrices

Description

Form the direct product of the m \times n matrix A and the p \times q matrix B. It is also called the Kroneker product and the right direct product. It is defined to be the result of replacing each element of A, a_{ij}, with a_{ij}\bold{B}. The result matrix is mp \times nq.

The method employed uses the rep function to form two mp \times nq matrices: (i) the direct product of A and J, and (ii) the direct product of J and B, where each J is a matrix of ones whose dimensions are those required to produce an mp \times nq matrix. Then the elementwise product of these two matrices is taken to yield the result.

Usage

mat.dirprod(A, B)

Arguments

Value

An mp \times nq matrix.

Author(s)

Chris Brien

See Also

matmult, mat.dirprod

Examples

    col.I <- mat.I(order=4)
    row.I <- mat.I(order=28)
    V <- mat.dirprod(col.I, row.I)

Forms the direct sum of a list of matrices

Description

The direct sum is the partitioned matrices whose diagonal submatrices are the matrices from which the direct sum is to be formed and whose off-diagonal submatrices are conformable matrices of zeroes. The resulting matrix is m \times n, where m is the sum of the numbers of rows of the contributing matrices and n is the sum of their numbers of columns.

Usage

mat.dirsum(matrices)

Arguments

Value

An m \times n matrix.

Author(s)

Chris Brien

See Also

mat.dirprod, matmult

Examples

       m1 <- matrix(1:4, nrow=2)
       m2 <- matrix(11:16, nrow=3)
       m3 <- diag(1, nrow=2, ncol=2)
       dsum <- mat.dirsum(list(m1, m2, m3))

Forms an exponential correlation matrix

Description

Form the correlation matrix of order equal to the length of coordinates. The matrix has diagonal elements equal to one and the off-diagonal element in the ith row and jth column equal to \rho^k where k = |coordinate[i]- coordinate[j]|.

Usage

mat.exp(rho, coordinates)

Arguments

Value

A correlation matrix whose elements depend on the power of the absolute distance apart.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.banded, mat.ar1, mat.ar2, mat.ar3, mat.sar2, mat.ma1, mat.ma2, mat.arma, mat.gau

Examples

    corr <- mat.exp(coordinates=c(3:6, 9:12, 15:18), rho=0.1)

Forms an exponential correlation matrix

Description

Form the correlation matrix of order equal to the length of coordinates. The matrix has diagonal elements equal to one and the off-diagonal element in the ith row and jth column equal to \rho^k where k = (coordinate[i]- coordinate[j])^2.

Usage

mat.gau(rho, coordinates)

Arguments

Value

A correlation matrix whose elements depend on the power of the absolute distance apart.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.banded, mat.ar1, mat.ar2, mat.ar3, mat.sar2, mat.ma1, mat.ma2, mat.arma, mat.exp

Examples

    corr <- mat.gau(coordinates=c(3:6, 9:12, 15:18), rho=0.1)

Computes the generalized inverse of a matrix

Description

Computes the Moore-Penrose generalized inverse of a matrix.

Usage

mat.ginv(x, tol = .Machine$double.eps ^ 0.5)

Arguments

Value

A matrix. An NA is returned if svd fails during the compution of the generalized inverse.

Author(s)

Chris Brien

Examples

## PBIBD(2) from p. 379 of Cochran and Cox (1957) Experimental Designs. 
## 2nd edn Wiley, New York
PBIBD2.unit <- list(Block = 6, Unit = 4)
PBIBD2.nest <- list(Unit = "Block")
trt <- factor(c(1,4,2,5, 2,5,3,6, 3,6,1,4, 4,1,5,2, 5,2,6,3, 6,3,4,1))
PBIBD2.lay <- designRandomize(allocated = trt, 
                              recipient = PBIBD2.unit, 
                              nested.recipients = PBIBD2.nest)
## Compute the projector for a linear trend across Blocks
PBIBD2.lay <- within(PBIBD2.lay,
                     {
                       cBlock <- as.numfac(Block)
                       cBlock <- cBlock - mean(unique(cBlock))
                     })
X <- model.matrix(~ cBlock, data = PBIBD2.lay)
Q.cB <- projector((X %*% mat.ginv(t(X) %*% X) %*% t(X)))

Forms an ma1 correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the ma1 pattern. The matrix is banded and has diagonal elements equal to one and subdiagonal element equal to
-MAparameter / (1 + MAparameter*MAparameter).

Usage

mat.ma1(MAparameter, order)

Arguments

Value

A banded correlation matrix whose elements follow an ma1 pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.exp, mat.gau, mat.banded, mat.ar2, mat.ar3, mat.sar2, mat.ma2, mat.arma

Examples

    corr <- mat.ma1(MAparameter=0.4, order=4)

Forms an ma2 correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the ma2 pattern. The resulting matrix is banded.

Usage

mat.ma2(MAparameters, order)

Arguments

Details

The correlations in the correlation matrix, corr say, are calculated from the moving average parameters, MAparameters. The values in

• the diagonal (k = 1) of corr are one;

• the first subdiagonal band (k = 2) of corr are equal to
  -MAparameters[1]*(1 - MAparameters[2]) / div;

• the second subdiagonal bande (k = 3) of corr are equal to -MAparameters[2] / div;

• in subsequent disgonal bands, (k = 4:order), of corr are zero,

where div = 1 + MMAparameters[1]*MAparameter[1] + MAparameters[2]*MAparameters[2].

Value

A banded correlation matrix whose elements follow an ma2 pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.exp, mat.gau, mat.banded, mat.ar1, mat.ar3, mat.sar2, mat.ma1, mat.arma

Examples

    corr <- mat.ma2(MAparameters = c(0.4, -0.2), order = 4)

Calculates the variance matrix of the random effects for a natural cubic smoothing spline

Description

Calculates the variance matrix of the random effects for a natural cubic smoothing spline. It is the tri-diagonal matrix \bold{G}_s given by Verbyla et al., (1999) multiplied by the variance component for the random spline effects.

Usage

mat.ncssvar(sigma2s = 1, knot.points, print = FALSE)

Arguments

Value

A matrix containing the variances and covariances of the random spline effects.

Author(s)

Chris Brien

References

Verbyla, A. P., Cullis, B. R., Kenward, M. G., and Welham, S. J. (1999). The analysis of designed experiments and longitudinal data by using smoothing splines (with discussion). Journal of the Royal Statistical Society, Series C (Applied Statistics), 48, 269-311.

See Also

Zncsspline.

Examples

Gs <- mat.ncssvar(knot.points = 1:10)

Calculates the variance matrix for the random effects from a mixed model, based on a supplied formula or a matrix.

Description

For n observations, compute the variance matrix of the random effects. The matrix can be specified using a formula for the random effects and a list of values of the variance components for the terms specified in the random formula. If a matrix specifying the variances of the nuisance random effects is supplied then it is returned as the value from the function.

Usage

mat.random(random, G, design, keep.order = TRUE)

Arguments

Details

If \bold{Z}_i is the is incidence matrix for the random nuisance effects in \bold{u}_i for a term in random and \bold{u}_i has variance matrix \bold{G}_i so that the contribution of the random effectst to the variance matrix for \bold{Y} is \bold{V}_u = \Sigma (\bold{Z}_i\bold{G}_i(\bold{Z}_i)^T).

Value

A n x n matrix containing the variance matrix for the random effects.

Author(s)

Chris Brien

See Also

mat.Vpredicts.

Examples

## Reduced example from Smith et al. (2015)
## Generate two-phase design
mill.fac <- fac.gen(list(Mrep = 2, Mday = 2, Mord = 3))
field.lay <- fac.gen(list(Frep = 2, Fplot = 4))
field.lay$Variety <- factor(c("D","E","Y","W","G","D","E","M"), 
                            levels = c("Y","W","G","M","D","E"))
start.design <- cbind(mill.fac, field.lay[c(3,4,5,8,1,7,3,4,5,8,6,2),])
rownames(start.design) <- NULL

## Set gammas
terms <- c("Variety", "Frep", "Frep:Fplot", "Mrep", "Mrep:Mday", "Mrep:Mday:Mord")
gammas <- c(1, 0.1, 0.2, 0.3, 0.2, 1)
names(gammas) <- terms

## Specify matrices to calculate the variance matrix of the predicted fixed Variety effects 
Vu <- with(start.design, fac.vcmat(Mrep, gammas["Mrep"]) + 
                         fac.vcmat(fac.combine(list(Mrep,Mday)), gammas["Mrep:Mday"]) + 
                         fac.vcmat(Frep, gammas["Frep"]) + 
                         fac.vcmat(fac.combine(list(Frep,Fplot)), gammas["Frep:Fplot"]))

## Calculate the variance matrix of the predicted random Variety effects using formulae
Vu <- mat.random(random = ~ -1 + Mrep/Mday + Frep/Fplot, 
                 G = as.list(gammas[c(4,5,2,3)]), 
                 design = start.design)

Forms an sar correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the sar pattern. The resulting matrix is banded.

Usage

mat.sar(SARparameter, order)

Arguments

Details

The values of the correlations in the correlation matrix, corr say, are calculated from the SARparameter, gamma as follows. The values in

• the diagonal of corr (k = 1) are one;

• the first subdiagonal band (k = 2) of corr are equal to gamma/(1 + (gamma * gamma / 4));

• the subsequent subdiagonal bands, (k = 3:order), of corr are equal to
  gamma * corr[k-1] - (gamma * gamma/4) * corr[k-2].

Value

A banded correlation matrix whose elements follow an sar pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.banded, mat.exp, mat.gau, mat.ar1, mat.ar2, mat.ar3, mat.sar2, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.sar(SARparameter = -0.4, order = 4)

Forms an sar2 correlation matrix

Description

Form the correlation matrix of order order whose correlations follow the sar2 pattern, a pattern used in crop competition models. The resulting matrix is banded and is a constrained AR3 matrix.

Usage

mat.sar2(gamma, order, print = NULL)

Arguments

Details

The values of the AR3 parameters, phi, are calculated from the gammas as follows:
phi[1] = gamma[1] + 2 * gamma[2]; phi[2] = -gamma[2] * (2*gamma[2] + gamma[1]);
phi[3] = gamma[1] * gamma[2] * gamma[2].

Then the correlations in the correlation matrix, corr say, are calculated from the correlation parameters, phi. Let omega = 1 - phi[2] - phi[3] * (phi[1] + phi[3]). Then the values in

• the diagonal of corr (k = 1) are one;

• the first subdiagonal band (k = 2) of corr are equal to (phi[1] + phi[2]*phi[3]) / omega;

• the second subdiagonal band (k = 3) of corr are equal to
  (phi[1] * (phi[1] + phi[3]) + phi[2] * (1 - phi[2])) / omega;

• the subsequent subdiagonal bands, (k = 4:order), of corr are equal to
  phi[1]*corr[k-1] + phi[2]*corr[k-2] + phi[3]*corr[k-3].

Value

A banded correlation matrix whose elements follow an sar2 pattern.

Author(s)

Chris Brien

See Also

mat.I, mat.J, mat.cor, mat.corg, mat.banded, mat.exp, mat.gau, mat.ar1, mat.ar2, mat.ar3, mat.sar, mat.ma1, mat.ma2, mat.arma

Examples

    corr <- mat.sar2(gamma = c(-0.4, 0.2), order = 4)
    corr <- mat.sar2(gamma = c(-0.4, 0.2), order = 4, print = "ar3")

computes the projection matrix that produces means

Description

Replaced by fac.meanop.

Converts the first two levels of a factor into the numeric values -1 and +1

Description

Converts the first two levels of a factor into the numeric values -1 and +1.

Usage

mpone(factor)

Arguments

Value

A numeric vector.

Warning

If the factor has more than two levels they will be coerced to numeric values.

Author(s)

Chris Brien

See Also

mpone in package dae, factor, relevel.

Examples

## generate all combinations of two two-level factors
mp <- c("-", "+")
Frf3.trt <- fac.gen(list(A = mp, B = mp))

## add factor C, whose levels are the products of the levles of A and B
Frf3.trt$C <- factor(mpone(Frf3.trt$A)*mpone(Frf3.trt$B), labels = mp)

Computes the number of replicates for an experiment

Description

Computes the number of pure replicates required in an experiment to achieve a specified power.

Usage

no.reps(multiple=1., df.num=1.,
        df.denom=expression((df.num + 1.) * (r - 1.)), delta=1.,
        sigma=1., alpha=0.05, power=0.8, tol=0.1, print=FALSE)

Arguments

Value

A list containing nreps, a single numeric value containing the computed number of pure replicates, and power, a single numeric value containing the power for the computed number of pure replicates.

Author(s)

Chris Brien

See Also

power.exp, detect.diff in package dae.

Examples

## Compute the number of replicates (blocks) required for a randomized 
## complete block design with four treatments. 
no.reps(multiple = 1, df.num = 3,
        df.denom = expression(df.num * (r - 1)), delta = 5,
	        sigma = sqrt(20), print = TRUE)

Description of a p2canon object

Description

An object of class p2canon that contains information derived from two formulae using projs.2canon.

Value

A list of class p2canon. It has two components: decomp and aliasing. The decomp component iscomposed as follows:

• It has a component for each component of Q1.

• Each of the components for Q1 is a list; each of these lists has one component for each of Q2 and a component Pres.

• Each of the Q2 components is a list of three components: pairwise, adjusted and Qproj. These components are based on an eigenalysis of the relationship between the projectors for the parent Q1 and Q2 components.

  1. Each pairwise component is based on the nonzero canonical efficiency factors for the joint decomposition of the two parent projectors (see proj2.eigen).

  2. An adjusted component is based on the nonzero canonical efficiency factors for the joint decomposition of the Q1 component and the Q2 component, the latter adjusted for all Q2 projectors that have occured previously in the list.

  3. The Qproj component is the adjusted projector for the parent Q2 component.

• The pairwise and adjusted components have the following components: efficiencies, aefficiency, mefficiency, sefficiency, eefficiency, xefficiency, order and dforthog – for details see efficiency.criteria.

The aliasing component is a data.frame decribing the aliasing between terms corresponding to two Q2 projectors when estimated in subspaces corresponding to a Q1 projector.

Author(s)

Chris Brien

See Also

projs.2canon, designAnatomy, pcanon.object.

Description of a pcanon object

Description

An object of class pcanon that contains information derived from several formulae using designAnatomy.

Value

A list of class pcanon that has four components: (i) Q, (ii) terms, (iii) sources, (iv) marginality, and (v) aliasing. Each component is a list with as many components as there are formulae in the formulae list supplied to designAnatomy.

The Q list is made up of the following components:

1. The first component is the joint decomposition of two structures derived from the first two formulae, being the p2canon.object produced by projs.2canon.

2. Then there is a component for each further formulae; it contains the p2canon.object obtained by applying projs.2canon to the structure for a formula and the already established joint decomposition of the structures for the previous formulae in the formulae.

3. The last component contains the the list of the projectors that give the combined canonical decomposition derived from all of the formulae.

The terms, sources, marginalty and aliasing lists have a component for each formula in the formulae argument to designAnatomy, Each component of the terms and sources lists has a character vector containing the terms or sources derived from its formula. For the marginality component, each component is the marginality matrix for the terms derived from its formula. For the aliasing component, each component is the aliasing data.frame for the source derived from its formula. The components of these four lists are produced by pstructure.formula and are copied from the pstructure.object for the formula. The names of the components of these four lists will be the names of the components in the formulae list.

The object has the attribute labels, which is set to "terms" or "sources" according to which of these were used to label the projectors when the object was created.

Author(s)

Chris Brien

See Also

designAnatomy, p2canon.object.

Takes a list of projectors and constructs a pstructure.object that includes projectors, each of which has been orthogonalized to all projectors preceding it in the list.

Description

Constructs a pstructure.object that includes a set of mutually orthogonal projectors, one for each of the projectors in the list. These specify a structure, or an orthogonal decomposition of the data space. This function externalizes the process previously performed within pstructure.formula to orthogonalize projectors. There are three methods available for carrying out orthogonalization: differencing, eigenmethods or the default hybrid method.

It is possible to use this function to find out what sources are associated with the terms in a model and to determine the marginality between terms in the model. The marginality matrix can be saved.

Usage

## S3 method for class 'list'
porthogonalize(projectors, formula = NULL, keep.order = TRUE, 
               grandMean = FALSE, orthogonalize = "hybrid", labels = "sources", 
               marginality = NULL, check.marginality = TRUE, 
               omit.projectors = FALSE, 
               which.criteria = c("aefficiency","eefficiency","order"), 
               aliasing.print = TRUE, ...)

Arguments

Details

It is envisaged that the projectors in the list supplied to the projectors argument correspond to the terms in a linear model. One way to generate them is to obtain the design matrix X for a term and then calculate its projector as \mathbf{X(X'X)^-X'}, There are three methods available for orhtogonalizing the supplied projectors: differencing, eigenmethods or the default hybrid method.

Differencing relies on comparing the factors involved in two terms, one previous to the other, to identify whether to subtract the orthogonalized projector for the previous term from the primary projector of the other. It does so if factors/variables for the previous term are a subset of the factors/variablesfor for the other term. This relies on ensuring that all projectors whose factors/variables are a subset of the current projector occur before it in the expanded formula. It is checked that the set of matrices are mutually orthogonal. If they are not then a warning is given. It may happen that differencing does not produce a projector, in which case eigenmethods must be used.

Eigenmethods forces each projector to be orthogonal to all terms previous to it in the expanded formula. It uses equation 4.10 of James and Wilkinson (1971), which involves calculating the canonical efficiency factors for pairs of primary projectors. It produces a table of efficiency criteria for partially aliased terms. Again, the order of terms is crucial. This method has the disadvantage that the marginality of terms is not determined and so sources names are set to be the same as the term names, unless a marginality matrix is supplied.

The hybrid method is the most general and uses the relationships between the projection operators for the terms in the formula to decide which projectors to subtract and which to orthogonalize using eigenmethods. If \mathbf{Q}_i and \mathbf{Q}_j are two projectors for two different terms, with i < j, then

1. if \mathbf{Q}_j\mathbf{Q}_i \neq \mathbf{0} then have to orthogonalize \mathbf{Q}_j to \mathbf{Q}_i.

2. if \mathbf{Q}_j\mathbf{Q}_i = \mathbf{Q}_j then, if \mathbf{Q}_i = \mathbf{Q}_j, they are equal and \mathbf{Q}_j will be removed from the list of terms; otherwise they are marginal and \mathbf{Q}_i is subtracted from \mathbf{Q}_j.

3. if have to orthogonalize and \mathbf{Q}_j\mathbf{Q}_i = \mathbf{Q}_i then \mathbf{Q}_j is aliased with previous terms and will be removed from the list of terms; otherwise \mathbf{Q}_i is partially aliased with \mathbf{Q}_j and \mathbf{Q}_j is orthogonalized to \mathbf{Q}_i using eigenmethods.