Math 3083 - 1 Coal Cleaning Example: 2^3 Experiment Feb. 22, 2016 with Replications via Yates Method Data from a VPI study of a filtering system for coal. A coagulant was added to coal sludge which was washed in a recirculation system. Three factors were considered. The response is percent solids by weight in the underflow, which is a measure of how clean the coal has become. The three factors at HI/LOW levels were A percent solids circulated initially in the underflow B flow rate of the polymer C pH of the tank Two replicates were taken for each experimental condition. Example from Walpole, Myers & Myers, Probability and Statistics for Engineers and Scientists, 6th ed., Prentice Hall 1998. ============================OUTPUT FROM R==================================== R version 2.13.1 (2011-07-08) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i386-apple-darwin9.8.0/i386 (32-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. [R.app GUI 1.41 (5874) i386-apple-darwin9.8.0] [Workspace restored from /Users/andrejstreibergs/.RData] [History restored from /Users/andrejstreibergs/.Rapp.history] %============ENTER DATA AND DEFINE FACTORS==================== > x=scan() 1: 4.65 5.81 21.42 21.35 12.66 12.56 18.27 16.62 9: 7.93 7.88 13.18 12.87 6.51 6.26 18.23 17.83 17: Read 16 items > matrix(x,byrow=T,ncol=2) [,1] [,2] [1,] 4.65 5.81 [2,] 21.42 21.35 [3,] 12.66 12.56 [4,] 18.27 16.62 [5,] 7.93 7.88 [6,] 13.18 12.87 [7,] 6.51 6.26 [8,] 18.23 17.83 > repl=factor(rep(1:2,times=8));repl [1] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 Levels: 1 2 > A=factor(rep(c(-1,-1,1,1),times=4));A [1] -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 Levels: -1 1 > B=factor(rep(c(-1,-1,-1,-1,1,1,1,1),times=2));B [1] -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 1 1 1 1 Levels: -1 1 > C=factor(c(rep(-1,times=8),rep(1,times=8)));C [1] -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 1 Levels: -1 1 %=========COMPUTE CELL SUMS============================ > xm=xtabs(x~A+B+C) > xm , , C = -1 B A -1 1 -1 10.46 25.22 1 42.77 34.89 , , C = 1 B A -1 1 -1 15.81 12.77 1 26.05 36.06 > xx=as.vector(xm) > cbind(matrix(c(x),byrow=T,ncol=2,dimnames=list( c("1","a","b","ab","c","ac","bc","abc"),c("Rep 1","Rep 2"))),xx) Rep 1 Rep 2 xx 1 4.65 5.81 10.46 a 21.42 21.35 42.77 b 12.66 12.56 25.22 ab 18.27 16.62 34.89 c 7.93 7.88 15.81 ac 13.18 12.87 26.05 bc 6.51 6.26 12.77 abc 18.23 17.83 36.06 %=========DO YATES METHOD TO COMPUTE EFFECTS CONTRASTS (X3)=========== > Yates=function(z){c(z[2]+z[1],z[4]+z[3],z[6]+z[5],z[8]+z[7], z[2]-z[1], z[4]-z[3],z[6]-z[5],z[8]-z[7])} > x1=Yates(xx);x2=Yates(x1);x3=Yates(x2) %==========COMPUTE SUM-SQUARES FROM EFFECT CONTRASTS====================== > x3sq=x3*x3/16 %==============TABLE OF YATES STEPS AND SUM-SQUARES======================= > cbind(xx,x1,x2,x3,x3sq) xx x1 x2 x3 x3sq [1,] 10.46 53.23 113.34 204.03 2.601765e+03 [2,] 42.77 60.11 90.69 75.51 3.563600e+02 [3,] 25.22 41.86 41.98 13.85 1.198891e+01 [4,] 34.89 48.83 33.53 -9.59 5.748006e+00 [5,] 15.81 32.31 6.88 -22.65 3.206391e+01 [6,] 26.05 9.67 6.97 -8.45 4.462656e+00 [7,] 12.77 10.24 -22.64 0.09 5.062500e-04 [8,] 36.06 23.29 13.05 35.69 7.961101e+01 %=================RUN ANOVA. COMPARE SUM-SQUARES========================= > a1=aov(x~A*B*C) > print(a1) Call: aov(formula = x ~ A * B * C) Terms: A B C A:B A:C Sum of Squares 356.3600 11.9889 32.0639 5.7480 4.4627 Deg. of Freedom 1 1 1 1 1 B:C A:B:C Residuals Sum of Squares 0.0005 79.6110 2.2020 Deg. of Freedom 1 1 8 %================COMPUTE EFFECTS. COMPARE WITH MODEL TABLES================ > x3/16 [1] 12.751875 4.719375 0.865625 -0.599375 -1.415625 -0.528125 [7] 0.005625 2.230625 > model.tables(a1) Tables of effects A A -1 1 -4.719 4.719 B B -1 1 -0.8656 0.8656 C C -1 1 1.4156 -1.4156 A:B B A -1 1 -1 -0.5994 0.5994 1 0.5994 -0.5994 A:C C A -1 1 -1 -0.5281 0.5281 1 0.5281 -0.5281 B:C C B -1 1 -1 0.005625 -0.005625 1 -0.005625 0.005625 A:B:C , , C = -1 B A -1 1 -1 -2.2306 2.2306 1 2.2306 -2.2306 , , C = 1 B A -1 1 -1 2.2306 -2.2306 1 -2.2306 2.2306 %=============COMPUTE ANOVA TABLE "BY HAND" FROM COMPUTED EFFECTS============ > df=c(1,1,1,1,1,1,1,8,15);SST=sum(x*x-sum(x)^2/256);SST [1] 492.437 > SSE=SST-sum(x3sq[2:8]);SSE [1] 2.20205 > MS=c(x3sq[-1],SSE/8) > F=c(MS[1:7]/MS[8],-1,-1) > PV=c(pf(F[1:7],1,8,lower.tail=FALSE),-1,-1);PV [1] 3.899232e-10 1.694476e-04 1.826448e-03 4.788462e-06 [5] 3.806531e-03 9.668436e-01 1.450796e-07 -1.000000e+00 [9] -1.000000e+00 > matrix(c(df,x3sq[-1],SSE,SST,MS,-1,F,PV),ncol=5, dimnames= list(c("a","b","ab","c","ac","bc","abc","Error","Total"), c("DF","SS","MS","F","P-Value"))) DF SS MS F P-Value a 1 356.36000625 356.36000625 1.294648e+03 3.899232e-10 b 1 11.98890625 11.98890625 4.355544e+01 1.694476e-04 ab 1 5.74800625 5.74800625 2.088238e+01 1.826448e-03 c 1 32.06390625 32.06390625 1.164875e+02 4.788462e-06 ac 1 4.46265625 4.46265625 1.621273e+01 3.806531e-03 bc 1 0.00050625 0.00050625 1.839195e-03 9.668436e-01 abc 1 79.61100625 79.61100625 2.892251e+02 1.450796e-07 Error 8 2.20205000 0.27525625 -1.000000e+00 -1.000000e+00 Total 15 492.43704375 -1.00000000 -1.000000e+00 -1.000000e+00 %==================COMPARE TO================================================= > anova(a1) Analysis of Variance Table Response: x Df Sum Sq Mean Sq F value Pr(>F) A 1 356.36 356.36 1294.6482 3.899e-10 *** B 1 11.99 11.99 43.5554 0.0001694 *** C 1 32.06 32.06 116.4875 4.788e-06 *** A:B 1 5.75 5.75 20.8824 0.0018264 ** A:C 1 4.46 4.46 16.2127 0.0038065 ** B:C 1 0.00 0.00 0.0018 0.9668436 A:B:C 1 79.61 79.61 289.2251 1.451e-07 *** Residuals 8 2.20 0.28 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 %=============CHECK IF NORMALITY OF ERRORS IS REASONABLE============== > shapiro.test(rstandard(a1)) Shapiro-Wilk normality test data: rstandard(a1) W = 0.9267, p-value = 0.2158 > plot(a1)