Math 558 Lecture #29
More on BIBDs in R The package “ibd”
In this lecture we will use the package “ibd”.
Copyright By PowCoder代写 加微信 powcoder
help(ibd) will provide a lot of information about the available functions. I will discuss some in today’s lecture. To start with we will install the package first 1
install.packages(“ibd”)
library(ibd)
Design1 <- bibd(v = 7, b = 7, r = 3, k = 3, lambda = 1)
Here the symbol v is used for the treatments. Hence we can obtain the plan for any BIBD if it exists for the given parameters2.
1it was not easy for me to install. I had to update some packages and download and install XQuartz-2.8.1.dmg
2a BIBD with given parameters may not exist. We will discuss more about that in the coming lectures
More on BIBDs in R
The design that we will obtain is given below. [,1] [,2] [,3]
Block-1 1 Block-2 5 Block-3 3 Block-4 1 Block-5 1 Block-6 2 Block-7 2
Concurrence matrix for the generated design
$NNP(Concurrence matrix)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 3 1 1 1 1 1 1 [2,] 1 3 1 1 1 1 1 [3,] 1 1 3 1 1 1 1 [4,] 1 1 1 3 1 1 1 [5,] 1 1 1 1 3 1 1 [6,] 1 1 1 1 1 3 1 [7,] 1 1 1 1 1 1 3
Concurrence matrix of a design
The matrix given in the previous slide is called the concurrence matrix of the design. This matrix gives the total number of occurrences of all the treatments in the diagonal and pair wise concurrences in the off diagonal places. We can denote this matrix by M. The general form of the matrix is given below
This a t × t matrix
r λ λ .... λ λ r λ .... λ ..... ...... ..... .... ....
M = ..... ...... ..... .... ... ..... ...... .... λ r
Incidence matrix
Another matrix that you will see in the R output is the incidence matrix given below
1 0 0 1 1 0 0 0 0 0 0 1 1 1 0 0 1 1 0 0 1
N=1 0 1 0 0 1 0 1 1 0 0 0 0 1 0 1 0 1 0 1 0
Incidence matrix
The incident matrix is given by N = (nij)t×b
1 if treatment i is in block j 0 otherwise
It can be seen that M = NNT
Note: It is the design in which a treatment occurs at most once in a block. A BIBD is always a binary design.
yij = μ + bi + τj + εij
which is identical to the model for a randomized complete block design given . However, the analysis is slightly different due to the missing observations as the blocks are incomplete.
Hereyij is the ith observation in the jth block, μi is the overall mean, τi is the effect of the ith treatment, bj is the effect of the jth block, and εij random error component. Also, εij ∼iid N(0, σ2)
Statistical Analysis
The total variability in the data is given by SST 2 y2..
SST = ∑∑yij − N ij
The total variability can be partition as
SST = SStr(adjusted) + SSB + SSE
The sum of squares for treatments are adjusted each treatment appears in a different set of r blocks. Therefore the difference between the treatment totals (unadjusted), yi. for i = 1, 2, ..t are effected by the difference in blocks as well. The block sum of squares are
1k y2 SSB= ∑y2.j−..
Statistical Analysis
y.j = y1j + y2j + ..ykj is the total in jth block. Sum of squares for blocks has b − 1 degrees of freedom. Define the quantity Qi as
Qi =yi.−k∑nijy.j i=1,2..t
Where nij = 1 if treatment i appears in block j and 0 otherwise
The adjusted treatment sum of squares is SSTr(adjusted) = k ∑ti=1 Q2i
It has t − 1 degrees of freedom. The SSE has
N − 1 − (a − 1) − (b − 1) = N − a − b + 1 degrees of freedom. The appropriate test statistic is
F = MSTr(adjusted) MSE
Balanced Incomplete Block Design for Catalyst Experiment
Treatment Block (Batch of Raw Material)
1 2 3 4 yi.
(Catalyst)
1 73 74 — 71 218 2 — 75 67 72 214 3 73 75 68 — 216 4 75 — 72 75 222
y.j 221 224 207 218 870 y.
Statistical Analysis: Catalyst example
This is a BIBD with t = 4, b = 4, k = r = 3, λ = 2, and N = 12 experimental units.
2 y2.. 8702 SST=∑∑yij−N=63.156− 12 =81.00
1k y2 SSB= ∑y2.j− ..
=1[(221)2 + (207)2 + (224)2 + (218)2] − 8702
3 12 =55.00
Statistical Analysis
Q1 = (218) − 1/3(221 + 224 + 218) = −9/3 Q2 =−7/3 Q3 =−4/3 Q4 =20/3
SSTr(adjusted) = k ∑ti=1 Q2i λt
= 3[(−9/3)2 + (−7/3)2 + (−4/3)2 + (20/3)2 ] 2×4
SSE = 81.00 − 22.75 − 55.00 = 3.25 F = 11.6, p − value = 0.107
ANOVA table
Statistical analysis in R data = taste different codes
Model1<-aov(score panelist + recipe, data = taste) > summary(Model1)
Df Sum Sq Mean Sq F value Pr(>F)
panelist 11 19.333 1.7576 2.301 0.1106
recipe 3 9.125 3.0417 3.982 0.0465 * Residuals 9 6.875 0.7639
Further ivestigation install.packages(“multcomp”) library(multcomp)
contr <- glht(Model1, linfct = mcp(recipe = "Tukey")) summary(contr, test = adjusted("none"))
Estimate B-A==0 0.750 C-A==0 1.375
D-A==0 -0.625 C-B==0 0.625 D-B==0 -1.375 D-C==0 -2.000
Std. Error 0.618 0.618 0.618 0.618 0.618 0.618
t value 1.214 2.225 -1.011 1.011 -2.225 -3.236
Pr(>|t|) 0.2558 0.0531 0.3383 0.3383 0.0531 . 0.0102 *
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com