L13-14-15 Ch9 R (sta305-class18-26Nov2019)

Today’s Class I Sample size for ANOVA I Randomized block designs I Linear model and ANOVA I Assumptions I Other Blocking Designs I Latin Square I Graeco Latin Square I hypo-Graeco Latin Square I Randomized incomplete block design

Sample size for ANOVA - Designing a study to compare more than two treatments I When H 0 is false it can be shown that: I MS Treat /σ 2 has a non-central Chi-square distribution with k - 1 degrees of freedom and non-centrality parameter δ . I MS Treat / MS E has a non-central F distribution with the numerator and denominator degrees of freedom k - 1 and N - k respectively, and non-centrality parameter I δ = ∑ k i = 1 n i ( μ i - ¯ μ ) 2 σ 2 , where n i is the number of observations in group i , ¯ μ = ∑ k i = 1 μ i / k , and σ 2 is the within group error variance . I This is dentoted by F k - 1 , N - k ( δ ) .

Blood coagulation example - sample size I μ 1 = 61, μ 2 = 66, μ 3 = 68, μ 4 = 61. I The error variance σ 2 was estimated as MS E = 5 . 6. I Assuming that the estimated values are the true values of the parameters, the non-centrality parameter of the F distribution is: I δ = 3 × ( ( 61 - 64 ) 2 + ( 66 - 64 ) 2 + ( 68 - 64 ) 2 + ( 61 - 64 ) 2 ) / 5 . 6 = 20 . 35714

Calculating power and sample size using the pwr library In the previous example δ = 20 . 35714 so f = √ 20 . 35714 = 4.5118887. library (pwr) pwr.anova.test ( k = 4 , n = 3 , f = 4.5 ) ## ## Balanced one-way analysis of variance power calculation ## ## k = 4 ## n = 3 ## f = 4.5 ## sig.level = 0.05 ## power = 1 ## ## NOTE: n is number in each group

Calculating power using simulation - R program #Simulate power of ANOVA for three groups NSIM <- 1000 # number of simulations res <- numeric (NSIM) # store p-values in res mu1 <- 2 ; mu2 <- 2.5 ;mu3 <- 2 # true mean values of treatment groups sigma1 <- 1 ; sigma2 <- 1 ; sigma3 <- 1 #variances in each group n1 <- 40 ; n2 <- 40 ; n3 <- 40 #sample size in each group for (i in 1 : NSIM) # do the calculations below N times { # generate sample of size n1 from N(mu1,sigma1^2) y1 <- rnorm ( n = n1, mean = mu1, sd = sigma1) # generate sample of size n2 from N(mu2,sigma2^2) y2 <- rnorm ( n = n2, mean = mu2, sd = sigma2) # generate sample of size n3 from N(mu3,sigma3^2) y3 <- rnorm ( n = n3, mean = mu3, sd = sigma3) y <- c (y1,y2,y3) # store all the values from the groups # generate the treatment assignment for each group trt <- as.factor ( c ( rep ( 1 ,n1), rep ( 2 ,n2), rep ( 3 ,n3))) m <- lm (y ~ trt) # calculate the ANOVA res[i] <- anova (m)[ 1 , 5 ] # p-value of F test } sum (res <= 0.05 ) / NSIM # calculate p-value ## [1] 0.642

The ANOVA identity for randomized block designs The total sum of squares can be re-expressed by adding and subtracting the treatment and block averages as: a i = 1 b j = 1 ( y ij - ¯ y ·· ) 2 = a i = 1 b j = 1 [( ¯ y i · - ¯ y ·· ) + ( ¯ y · j - ¯ y ·· ) + ( y ij - ¯ y i · - ¯ y · j + ¯ y ·· ))] 2 . After some algebra . . . SS T = ∑ a i = 1 ∑ b j = 1 ( y ij - ¯ y ·· ) 2 is equal to b a i = 1 ( ¯ y i · - ¯ y ·· ) 2 + a b j = 1 ( ¯ y · j - ¯ y ·· ) 2 + a i = 1 b j = 1 ( y ij - ¯ y i · - ¯ y · j + ¯ y ·· ) 2 So, SS T = SS Treat + SS Blocks + SS E

Linear Model for Randomized Block Design pen.model <- lm (y ~ as.factor (treatment) + as.factor (blend), data= tab0404) anova (pen.model) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) as.factor(treatment) 3 70 23.333 1.2389 0.33866 as.factor(blend) 4 264 66.000 3.5044 0.04075 * Residuals 12 226 18.833 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 Calculation of the p-value assumes that ij ∼ N ( 0 , σ 2 ) . So that MS Treat / MS E ∼ F a - 1 , ( a - 1 )( b - 1 ) , MS Blocks ∼ F b - 1 , ( a - 1 )( b - 1 ) .

Penicillin example - interpretation I There is no evidence that the four treatments produce different yields. I How could this information be used in optimizing yield in the manufacturing process? I Is one of the treatments less expensive to run? I If one of the treatments is less expensive to run then an analysis on cost rather than yield might reveal important information. I The differences between the blocks might be informative. I In particular the investigators might speculate about why blend 1 has such a different influence on yield. I Perhaps now the experimenters should study the characteristics of the different blends of corn steep liquor. (Box, Hunter, Hunter, 2005)

Other blocking designs I Latin square I Graeco-Latin squares, I Hyper-Graeco-Latin Squares, I Balanced incomplete block designs.

The Latin Square Design I There are several other types of designs that utilize the blocking principle such as The Latin Square design. I If there is more than one nuisance source that can be eliminated then a Latin Square design might be appropriate.

Latin Square Design - Automobile Emissions I An experiment to test the feasibility of reducing air pollution. I A gasoline mixture was modified by changing the amounts of certain chemicals. I This produced four different types of gasoline: A, B, C, D I These four treatments were tested with four different drivers and four different cars.

Latin Square Design - Automobile Emissions I Two blocking factors: cars and drivers. I The Latin square design was used to help eliminate possible differences between drivers I, II, III, IV and cars 1, 2, 3, 4. I Randomly allocate treatments, drivers , and cars. Driver Car 1 Car 2 Car 3 Car 4 Driver I A B D C Driver II D C A B Driver III B D C A Driver IV C A B D

Latin Square Design - Automobile Emissions I The data from the experiment. Driver Car 1 Car 2 Car 3 Car 4 Driver I A B D C 19 24 23 26 Driver II D C A B 23 24 19 30 Driver III B D C A 15 14 15 16 Driver IV C A B D 19 18 19 16

Latin Square Design - Automobile Emissions I Why not standardize the conditions and make the 16 experimental runs with a single car and single driver for the four treatments? I Could also be valid but Latin square provides a wider inductive basis.

Latin Square Design - Automobile Emissions latinsq.auto <- lm (y ~ additive + as.factor (cars) + as.factor (driver), data= tab0408) anova (latinsq.auto) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) additive 3 40 13.333 2.5 0.156490 as.factor(cars) 3 24 8.000 1.5 0.307174 as.factor(driver) 3 216 72.000 13.5 0.004466 ** Residuals 6 32 5.333 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1 SS T = SS cars + SS drivers + SS Additives + SS E

Latin Square Design - Automobile Emissions I Assumming that the residuals are independent and normally distributed and the null hypothesis that there are no treatment differences is true then the ratio of mean squares for treatments and residuals has an F 3 , 6 distribution. I This analysis assumes that treatments, cars, and drivers are additive. I If the design was replicated then this would increase the degrees of freedom for the residuals and reduce the mean square error.

General Latin Square I A Latin square for p factors of a p × p Latin square, is a square containing p rows and p columns I Each of the p 2 cells contains one of the p letters that correspond to a treatment. I Each letter occurs once and only once in each row and column. I There are many possible p × p Latin squares.

General Latin Square Which of the following is a Latin square? Col1 Col2 Col3 Row 1 B A C Row 2 A C B Row 3 C B A Col1 Col2 Col3 Row 1 A B C Row 2 C A B Row 3 B B A

Misuse of the Latin Square I Inappropriate to use Latin square to study factors that can interact. I Effects of one factor can then be mixed up with interactions of other factors. I Outliers can occur as a result of these interactions. I When interactions between factors are likely possible need to use a factorial design.

Graeco-Latin Square A Graeco-Latin square is a k × k pattern that permits study of k treatments simultaneously with three different blocking variables each at k levels. Car 1 Car 2 Car 3 Car 4 Driver I A α B β C γ D δ Driver II B δ A γ D β C α Driver III C β D α A δ B γ Driver IV D γ C δ B α A β

Graeco-Latin Square I This is a Latin square in which each Greek letter appears once and only once with each Latin letter. I Can be used to control three sources of extraneous variability (i.e. block in three different directions). Driver Car 1 Car 2 Car 3 Car 4 Driver I A α B β C γ D δ Driver II B δ A γ D β C α Driver III C β D α A δ B γ Driver IV D γ C δ B α A β

Graeco-Latin Square To generate a 3 × 3 Graeco-Latin square design, superimpose two designs using the Greek letters for the second 3 × 3 Latin square. Col1 Col2 Col3 Row 1 B A C Row 2 A C B Row 3 C B A Col1 Col2 Col3 Row 1 A B C Row 2 C A B Row 3 B C A

hyper-Graeco-Latin Square These three Latin squares can be superimposed to form a hyper-Graeco-Latin square. Can be used to control 4 nuisance factors (i.e. block 4 factors). Row Col1 Col2 Col3 Col4 Row 1 B A D C Row 2 C D A B Row 3 D B C A Row 4 A C B D Row Col1 Col2 Col3 Col4 Row 1 D A C B Row 2 A D B C Row 3 B C A D Row 4 C B D A Row Col1 Col2 Col3 Col4 Row 1 A D B C Row 2 C A D B Row 3 B C A D Row 4 D B C A

hyper-Graeco-Latin Square I A machine used for testing the wear on types of cloth. I Four pieces of cloth can be compared simultaneously on one machine. I Response is weight loss in tenths of mg when rubbed against a standard grade of emery paper for 1000 revolutions of the machine.

hyper-Graeco-Latin Square I Specimens of 4 different cloths (A, B,C,D) are compared. I The wearing qualities can be in any one of 4 positions P 1 , P 2 , P 3 , P 4 on the machine. I Each emery ( α, β, γ, δ ) paper used to cut into for quarters and each quarter used to complete a cycle C 1 , C 2 , C 3 , C 4 of 1000 revolutions. I Object was to compare treatments.

hyper-Graeco-Latin Square i) type of specimen holders 1, 2, 3, 4 ii) position on the machine P 1 , P 2 , P 3 , P 4 . iii) emory paper sheet α, β, γ, δ . iv) machine cycle C 1 , C 2 , C 3 , C 4 . The design was replicated. The first replicate is shown in the table below. P 1 P 2 P 3 P 4 C 1 A α 1 B β 2 C γ 3 D δ 4 320 297 299 313 C 2 C β 4 D α 3 A δ 2 B γ 1 266 227 260 240 C 3 D γ 2 C δ 1 B α 4 A β 3 221 240 267 252 C 4 B δ 3 A γ 4 D β 1 C α 2 301 238 243 290

hyper-Graeco-Latin Square A linear model can be fit so that the ANOVA table and parameter treatment effects can be calculated. wear.hypsq <- lm (y ~ treatment + as.factor (rep) + as.factor (position) + as.factor (cycle) + as.factor (holder) + as.factor (paper), data= tab0412) anova (wear.hypsq) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) treatment 3 1705.3 568.45 5.3908 0.021245 * as.factor(rep) 1 603.8 603.78 5.7259 0.040366 * as.factor(position) 3 2217.3 739.11 7.0093 0.009925 ** as.factor(cycle) 6 14770.4 2461.74 23.3455 5.273e-05 *** as.factor(holder) 3 109.1 36.36 0.3449 0.793790 as.factor(paper) 6 6108.9 1018.16 9.6555 0.001698 ** Residuals 9 949.0 105.45 --- Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1 1

Balanced incomplete block design I Suppose that instead of four samples to be included on each 1000 revolution cycle only three could be included, but the experimenter still wanted to compare four treatments. I The size of the block is now 3 - too small to accommodate all treatments simultaneously. I A balanced incomplete block design has the property that every pair of treatments occurs together in a block the same number of times. Cycle block 1 A B C 2 A B D 3 A C D 4 B C D Cycle block A B C D 1 x x x 2 x x x 3 x x x 4 x x x