Variance component based parameter estimation of incomplete block designs

Introduction

Variance component models are also suited for analysis of incomplete block designs, besides complete block designs. This post aims to demonstrate exactly that. Using a dataset generated from alpha lattice design, I show how the design can be properly modeled and fit using OLS regression having various fixed model components. This system of model fitting is analogous to classical ANOVA based technique of estimating parameters.

The dataset

The data is same as introduced earlier in “Design and analysis of balanced randomized complete block (RCBD) experiment” post. It comprises of plant height trait for the maize recorded in multiple replication units, only difference from the earlier post being that the arrangement of observation plots into blocks which are themselves nested within replication will be accounted for. This is the objective of this post: to account for effects of incomplete blocks which was incorporated in the design.

An overview of how the data looks like is shown in Table 1.

Model formulation

model_pht <- lm(formula(paste("Plant_height", "~ Rep + Entry + Block%in%Rep")), 
            data = ihyb_multiple)
ANOVA_pht <- anova(model_pht) # model anova

mean_pht # overall mean

## [1] 267

Fstat <- data.frame("Fit Statistics" = c(AIC = AIC(model_pht), BIC = BIC(model_pht)))
Fstat # fit statistics

##     Fit.Statistics
## AIC           1178
## BIC           1419

E <- (ntr - 1) * (nrep - 1)/((ntr - 1) * (nrep - 1) + nrep * (s - 1)) 
# ntr = number of treatments, 
# nrep = number of replications
# s = number of blocks per replication 
E # efficiency factor

## [1] 0.784

CV <- sqrt(Ee) * 100/mean_pht # where, Ee = deviance(model_pht)/df_resid_pht 
CV # model CV 

## [1] 3.93

Alternatively, analysis and results from analysis could be obtained via agricolae function PBIB.test().

# vc_model_alpha <- agricolae::PBIB.test(block = block, trt = trt, replication = replication, 
#                      k = sblock, y = y, method = "VC", console = TRUE)
# vc_model_alpha$groups