Alpha lattice design in R studio

1. Background

The use of incomplete block design is used among researchers to control random variation in evaluation trails of many genotypes.  It has been predominately used by plant breeders for the evaluation of various genotypes.  It could even be useful where there is no factorial type relationship among treatments and the number of treatment is sufficiently large and soil heterogeneity is very high in test site.

1.       In breeding programs, the number of treatments may be large
2.      As blocks get larger, the conditions within blocks become more heterogeneous
3.      Other factors may limit the number of experimental units in a block

•     Size of a growth chamber, greenhouse, or available field area
•      Number of samples that can be processed at one time

Breeding lines are tested in micro varietal trails at research stations and substations and then evaluated at outstations . These trails determine the suitability, adaptability and disease responses of those selected varieties and are usually conducted in randomized completely block design (RCBD). However, using RCBD in such condition is risky because in large number of  treatments, maintaining homogeneity within blocks is very difficult. In such condition, varietal effects are estimated with less precision as the residual (error) mean square will be larger than actually it should be. To get rid of such condition and increase efficiency of overall research, incomplete block designs (IBD) are used as alternative.  IBD have better chances of adhering the homogeneous condition within blocks and allows recovering information through inter block analysis to improve the efficiency of research.  In this regard, Patterson and Williams developed generalized lattice designs called as alpha designs which are more concerned about the condition of research area and the number of experimental units available. 

Reduction of experimental error is one of the basic principle of any experimental design. And it is achieved through blocking. In RCBD and LSD, full set of treatments are kept in each block. When the number of varieties are too large, it is impossible to fit them all in one block of homogeneous experimental units. For this IBD divides each complete blocks into smaller blocks. 

2. Types of Incomplete block designs

1. Balanced incomplete block design

    Each treatment occurs together in the same block with every other treatment an equal number of times. 

    •t = number of treatments

     •k = number of units per block (block size)

     •b = total number of blocks in the experiment

     •r = number of replicates of each treatment

     •l= number of times that treatments occur together in the same bloc

All pairs of treatments are compared with the same level of precision even though differences between blocks may be large

For balance, r = l(t-1)/(k-1)

•t = number of treatments

•k = number of units per block (block size)

•b = total number of blocks in the experiment

•r = number of replicates of each treatment

•l = number of times that treatments occur together in the same block

•N = total number of experimental units

Lamda must be an integer; N = b*k = r*t

 Example:   t = 10 treatments with block size k = 4

             •r = 6, b = 15, l = 2  

                 N = b*k = r*t = 6*10 = 60

In plant breeding, the minimum number of replications required for balance is often too large to be practical. 

2. Partially balanced incomplete block designs 

Different treatment pairs occur in the same blocks an unequal number of times or some treatment pairs never occur together in the same block

•      Mean comparisons have differing levels of precision
•     Greater precision when treatments occur in the same block
•     Statistical analysis more complex

Common for plant breeding trials due to large number of entries

3. Resolvable incomplete block designs 

Blocks are grouped so that each group of blocks constitute one complete replication of the treatment

•“block” = incomplete block = “sub-block”

• replication = “super block”

Trials can be managed in the field on a rep-by-rep basis

Field operations can be conducted in stages (planting, weeding, data collection, harvest)

Complete replicates can be lost without losing the whole experiment

If you have two or more complete replications, you can analyze as an RCBD if the blocking turns out to be ineffective. 

3. Lattice designs are resolvable 

Lattice designs are a well-known type of resolvable incomplete block design. 

Here S*R=B

t = number of treatments =15 

k = number of units per block (block size)=  5 

b = total number of blocks in the experiment = 6

r = number of replicates of each treatment =2

s = number of blocks in each complete replication = 3

Type of lattice design 

Square lattice designs

•         Number of treatments must be a perfect square (t = k2)
•     Blocks per replicate (s) and plots per block (k) are equal (s = k) and are the square root of the number of treatments (t)
•      For complete balance, number of replicates (r) = k+1

Rectangular lattice designs

•        t = s*(s-1) and k = s-1
•     Example: 4 x 5 lattice has 4 plots per block, 5 blocks per replicate, and 20 treatments

Alpha lattices

•        t = s*k
•     more flexibility in choice of s and k

The use of alpha lattice design allows the adjustment of treatment means for block effects.  This in turn brings benefit from the small incomplete blocks which help varietal comparisons under more homogenous conditions. The alpha lattice design also provides effective control within replicate variability. 

Information about alpha lattice design: 

•        Used when there is large number of treatments and the area available is small.
•         There are no check varieties for estimation error. Error checks within the same varieties.
•         They are used to reduce the effect of within complete block variation.
•         They increases precision and provide repeatability in trails. 
•         The design permits the removal of incomplete block effects from the plot residuals.
•         Maximizes the use of comparison between genotypes in the same incomplete block. 

ANOVA in alpha lattice design

4. Alpha lattice design in R studio 

File used: Alpha_lattice.xlsx (Output are in blue color)
https://drive.google.com/file/d/1ctoRzEaveGQXygTBITH6hhzaMGC5MCWZ/view?usp=sharing
A. Import data set to R studio
B.Attach the data set
> attach(Alpha_lattice)
C. Explain block, replication  and treatment as factor
> data$Treatment=as.factor(Treatment)
> data$Replication=as.factor(Replication)
> data$Block=as.factor(Block)
D. See the structure of your dataset
> str(data)
tibble [50 x 4] (S3: tbl_df/tbl/data.frame)
 $ Treatment  : Factor w/ 25 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ Replication: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ Block      : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
 $ Yield      : num [1:50] 6 7 5 8 6 16 12 12 13 8 ...
E. Load Package: Agricolae
> library(agricolae)
F. Analysis
There are two ways of doing analysis. One is through PBIB test and another is through making a model and doing analysis. I will show the both ways.
1st way
> model<-aov(Yield~Replication+Treatment+Replication:Block, data=data)
> summary(model)
                                    Df            Sum Sq         Mean Sq         F value         Pr(>F)   
Replication                   1               212.2             212.18            15.539      0.00117 **
Treatment                   24                559.3             23.30               1.707       0.13579   
Replication:Block        8                501.8             62.73               4.594     0.00463 **
Residuals                   16                218.5             13.66                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Note: No Significant difference is observed for treatment.

> out<-HSD.test(Yield, Treatment, DFerror = 16, MSerror = 13.66, alpha = 0.05,console=TRUE)
Note: You can use LSD.test (for least significant difference) and duncan.test (for Duncan's Multiple Range Test) instead of HSD.test( 

Study: Yield ~ Treatment

HSD Test for Yield 

Mean Square Error:  13.66 

Treatment,  means

   Yield        std r Min Max
1   15.0    12.7279221 2   6  24
2   14.0     9.8994949 2   7  21
3   10.5    7.7781746 2   5  16
4   12.5    6.3639610 2   8  17
5   10.5    6.3639610 2   6  15
6   14.5    2.1213203 2  13  16
7   11.5    0.7071068 2  11  12
8    8.0    5.6568542 2   4  12
9   11.5  2.1213203 2  10  13
10  11.5  4.9497475 2   8  15
11  20.5  4.9497475 2  17  24
12  10.5  4.9497475 2   7  14
13   9.5  3.5355339 2   7  12
14  19.5 14.8492424 2   9  30
15  18.0  5.6568542 2  14  22
16  14.5  4.9497475 2  11  18
17  13.5  3.5355339 2  11  16
18  12.5  0.7071068 2  12  13
19  11.0  2.8284271 2   9  13
20  15.0  1.4142136 2  14  16
21  11.0  4.2426407 2   8  14
22  19.0  5.6568542 2  15  23
23  11.5  0.7071068 2  11  12
24  18.5  6.3639610 2  14  23
25  16.5  3.5355339 2  14  19

Alpha: 0.05 ; DF Error: 16 
Critical Value of Studentized Range: 6.125725 

Minimun Significant Difference: 16.00913 

Treatments with the same letter are not significantly different.

   Yield groups
11  20.5      a
14  19.5      a
22  19.0      a
24  18.5      a
15  18.0      a
25  16.5      a
1   15.0      a
20  15.0      a
6   14.5      a
16  14.5      a
2   14.0      a
17  13.5      a
4   12.5      a
18  12.5      a
7   11.5      a
9   11.5      a
10  11.5      a
23  11.5      a
19  11.0      a
21  11.0      a
3   10.5      a
5   10.5      a
12  10.5      a
13   9.5      a
8    8.0      a

2nd way: 

> modelLattice<-PBIB.test(Block,Treatment,Replication,Yield,k=5,console = TRUE,method =c("VC"),test = "tukey",alpha = 0.05,group = TRUE)

ANALYSIS PBIB:  Yield 

Class level information
Block : 10 
Treatment : 25

Number of observations:  50 

Estimation Method:  Variances component model 

    Fit Statistics
AIC       285.6274
BIC       352.5482

Analysis of Variance Table

Response: Yield
                                        Df         Sum Sq        Mean Sq       F value         Pr(>F)   
Replication                      1            212.18         212.180        15.5386     0.001166 **
Treatment.unadj            24            559.28           23.303          1.7066     0.135789   
Block/Replication           8            501.84           62.730          4.5939     0.004629 **
Residual                        16            218.48           13.655                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Coefficient of variation: 27.1 %
Yield Means: 13.62 

Parameters PBIB
                   .
Treatment         25
Block size         5
Block/Replication  5
Replication        2

Efficiency factor 0.75 

Comparison test tukey 

Treatments with the same letter are not significantly different.

    Yield.adj groups
11  53.564079      a
1   52.376785      a
3   43.564079     ab
13  38.251373    abc
15  27.344196    abc
21  24.656902   abcd
5   24.156902   abcd
6   22.047235   bcde
2   21.547235   bcde
23  20.844196  bcdef
14  16.984255  bcdef
16  16.296961  bcdef
4   14.296961 bcdefg
12  13.734529  cdefg
8   11.234529 cdefgh
18   9.984255 cdefgh
25   6.437019 cdefgh
22   2.827353 cdefgh
24  -3.422922  defgh
10  -4.672647  defgh
20  -6.922922   efgh
7  -10.782314    fgh
17 -14.532589     gh
9  -16.532589      h
19 -22.782863      h

How to plot the data
>bar.group(out$groups,ylim=c(0,20),xlab="treatment", ylab="Yield")

> bar.err(out$means,ylim=c(0,20),xlab="treatment", ylab="Yield")

> boxplot(Yield~Treatment,data = data,xlab = "Treatment",ylab = "yield")

Note: This analysis is for one year only. But for combined analysis (including, multiple site and year) you can use following command.

 Multiple sites:

>model.(name of y-variable)<-aov(name of y-variable~Rep+Entry+Loc+Entry*Loc+Rep:Block, data=name of excel data)

Multiple years:

>model.(name of y-variable)<-aov(name of y-variable~Rep+Entry+Year+Entry*Year+Rep:Block, data=name of excel data)

Multiple sites and multiple years:

>model.(name of y-variable)<-aov(name of y variable~Rep+Entry+Loc+Year+Entry*Loc+Entry*Year++Loc*Year+Entry*Loc*Year+ Rep:Block, data=name of excel data)