Bayesian and Classical Analysis of Multi-Stratum Response Surface Designs

1. Bayesian and Classical Analysis of Multi-Stratum Response Surface Designs Steven Gilmour Queen Mary, University of London Peter Goos Universiteit Antwerpen

2. Outline • Split-plot and other multi-stratum designs • “State-of-the-art” analysis of data • REML/generalized least squares • Problems • Estimation of variance components • Degrees of freedom • Three possible solutions • Fix value(s) of variance-component(s) • Use randomization-based estimation • Bayesian analysis

3. Multi-stratum designs • Randomization of treatments to experimental units is restricted in such a way that particular sets of units must receive the same level of one or more treatment factors • Includes classical orthogonal split-plot, split-split-plot, criss-cross, etc. designs (regular factorial treatment sets) • Also includes nonorthogonal designs with similar structures (irregular factorial or response surface treatment sets) • Are (nested) block designs with at least one main effect totally confounded with block effects • Often necessary when some factors are hard to change

4. Multi-stratumdesigns • I refer to the runs as units and the groups of units defined by the randomization restrictions as blocks, superblocks, … • Randomization is performed by randomly relabelling …, superblocks, blocks and units • Implies random effects for …, superblocks, blocks, units (error) in derived linear model • Fixed treatment effects can be modelled using usual polynomial response surface model

5. Freeze-dried coffee experiment • Response: amount of retained volatile compounds in freeze-dried coffee • Treatment factors: • Pressure in drying chamber (dial-controlled) • Heating temperature, Initial solids content, Slab thickness, Freezing rate (all easy to change) • 5 runs during each of 6 days • Randomization restricted so that all runs in a day have the same pressure

6. Freeze-dried coffee experiment

7. Model • Generalized least squares (GLS) estimation • Variance-covariance matrix Model and analysis

8. Model • Generalized least squares (GLS) estimation • Variance-covariance matrix Model and analysis ? ?

9. Variance component estimation • REML: REsidual Maximum Likelihood • Yields the same answers as ANOVA in orthogonal designs (e.g. standard split-plots) • Applicable when designs are not orthogonal (e.g. nonorthogonal split-plots) • State of the art in many disciplines • Available in many statistical software packages

10. Analysis • Different implementations: • Variance components allowed to be negative or not • Various methods for obtaining effective degrees of freedom • Estimates generally consistent with each other, given different implementations • Effective degrees of freedom can be inconsistent with each other • All methods can give surprising results

11. Freeze-dried coffee experiment • Estimates of whole-plot error variance • SAS: 0 • GenStat: 0 • R: 0.0051 • Degrees of freedom for testing linear effect of pressure (full model) • SAS proc mixed with Kenward & Roger: 9 df • SAS proc mixed with containment method: 6 df • R lme default: 3 df

12. Freeze-dried coffee experiment Simplified model:

13. Freeze-dried coffee experiment • Data are treated as if they come from a completely randomized experiment • OLS estimates are obtained • Degrees of freedom for testing linear effect of pressure are too optimistic • Upper bound for full second-order model: 3 df • Because of nonorthogonality: less than 3 df • SAS proc mixed with Kenward & Roger: 9 df • SAS proc mixed with containment method: 6 df • R lme default: 3 df

14. Artificial example

15. Artificial example

16. Solution II: Randomization-based analysis • Even nonorthogonal multi-stratum designs have simple orthogonal block structures (if each block/superblock/... is the same size) [Nelder, 1965] • Ignoring treatment structure, randomization-based analysis gives minimum variance unbiased estimators of variance components (pure error) • Only assumption is that treatment and unit effects are additive

17. Randomization-based analysis • Proposed analysis: • Use discrete treatments defined by combinations of factor levels (ignoring treatment model) • Anova gives correct estimates of variance components with correct degrees of freedom • Use these estimates to fit treatment model using GLS • Base inferences on these estimates • “Extra sums of squares” represent lack of fit • Not clear that GLS is best, but is same as with REML

18. Freeze-dried coffee experiment

19. Freeze-dried coffee experiment • There are 27 treatments, 3 replicated twice • 0 residual degrees of freedom for blocks • 3 residual degrees of freedom for runs • Blocks variance component cannot be estimated, unit variance badly estimated • Full polynomial model can be fitted, but no global inference is possible • Weak inference is possible for all individual parameters except main effects of pressure • This design is too small for a frequentist analysis

20. Solution III: Bayesian approach • Advantages: • Takes into account uncertainty in prior beliefs • Prior beliefs can be contradicted by the data • No problems determining the appropriate degrees of freedom for hypothesis tests • WinBUGS software is free

21. The Bayesian approach • Requires a user-specified (joint) distribution for all model parameters (, 2, 2) • Posterior marginal distributions can be used for inference about parameters • Results: • Similar to REML/GLS if data contain enough information • Similar to prior distribution if data don’t contain enough information

22. The Bayesian approach • Noninformative priors for r: N(0,) • Weakly informative priors for r: • Linear and interaction effects: N(0,25) • Quadratic effects: N(0,100)

23. The Bayesian approach • Variance components

24. weakly informative highly informative not informative The Bayesian approach • Variance components

25. Results: linear effect of pressure

26. Linear effect of temperature

27. Interaction of slab thickness and freezing rate

28. Summary of results • Prior information on  has little impact • Prior information on 2 not important at all • Some results strongly depend on prior information about 2 • Hard-to-change factor coefficients • Sub-plot factor interaction coefficients that are not nearly orthogonal to whole plots • Results for other coefficients insensitive to the choice of the prior for 2

29. Discussion • REML/GLS analysis can be misleading as it often leads to an analysis that ignores the multi-stratum nature of the design • Likelihood methods have good asymptotic properties, i.e. large numbers of units in each stratum, so should not be expected to work in small experiments • Problem is due to a lack of information in the blocks stratum • We should honestly admit that there is no information and/or provide prior information

30. Discussion • Randomization based analysis should always be done (in every experiment!) as a first step • Makes very few assumptions, so is much more robust than any other analysis • Provides a “reality check” • Might make extra assumptions unnecessary • Bayesian analysis can help • Prior information is taken into account • Prior information can be overruled • Depends heavily on prior assumptions, but these are clearly and honestly expressed

31. References • Multi-stratum response surface designs Luzia A. Trinca and Steven G. Gilmour Technometrics, 2001 • A split-unit response surface design for improving aroma retention in freeze dried coffee Steven G. Gilmour, J. Mauricio Pardo, Luzia A. Trinca, K. Niranjan and Don Mottram Proceedings of the 6th European Conference on Food-Industry and Statistics, 2000 • Analysis of data from unbalanced multi-stratum designs Steven G. Gilmour and Peter Goos Submitted