19.2 18.7 12.5 20.3 19.9 18.7 14.3 14.3 22.5 24.3 21.3 20.2 8.7 17.6 17.6 16.5 17.6 11.4 18.4 20.2 17.3 19.3 9.5 15.9 18.4 22.4 16.1 16.5 19.0 19.1 tab2.16=matrix(scan(),byrow=T, ncol=5) N=length(tab2.16) n=dim(tab2.16)[1] m=dim(tab2.16)[2] syms = matrix(rep(c(2,6,15,5,1), 6), byrow=T, ncol=5) ids = matrix(rep(c(1:5), 6), byrow=T, ncol=5) plot(c(tab2.16)) matplot(t(tab2.16), pch=t(syms)) ma=apply(tab2.16, 2, mean) xbar = paste(sep=".", "Xb", (1:m)) SEx = paste(sep=".", "SE", (1:m)) maxes = apply(tab2.16, 2, max) mins = apply(tab2.16, 2, min) ###################### PLOT the data figure 2:39 in Davis plot(range(1:6), range(tab2.16), type='n', xlab="sample", ylab="CaCO3, %", ann=FALSE, axes=FALSE) axis(2) axis(1, at=1:m) box() points(ids, tab2.16, pch=syms) abline(h=ma, lty=2, col=grey(.8)) title( xlab="sample", ylab="CaCO3 %", main="ANOVA") text((1:m)+.1, ma, labels=xbar, pos=4) arrows(5.7, min(ma), 5.7, max(ma), code=3, length=.1, col='blue') arrows(6.1, min(tab2.16), 6.1, max(tab2.16), code=3, length=.1) arrows((1:m)-.1, mins, (1:m)-.1, maxes, code=3, length=.1, col='red') text(6.1, mean(range(tab2.16)), "SST", srt=90) text(5.7, mean(range(ma)), "SSA", srt=90, col='blue') text((1:m)-.1, (mins+maxes)/2, labels=SEx, srt=90, col='red') ###################### corterm=(sum(tab2.16)^2)/length(tab2.16) sst=sum(tab2.16^2)-corterm a=apply(tab2.16, 2, sum) if(FALSE){ ###a = c(115.4 106.2 72.9 113.7 119.5) ### ### ja = sum(tab2.16)/N ### ma=apply(tab2.16, 2, mean) ### sa = sweep(tab2.16, 2, apply(tab2.16, 2, mean), "-") ### ra=apply(sa^2, 2, sum)/(6-1) ### ja = sum(tab2.16)/N ma=tab2.16-ja a=apply(ma^2, 2, sum) } ssa=sum( (a^2)/n)-corterm ssw=sst-ssa MSa=ssa/(m-1) MSw=ssw/(N-m) Ft=MSa/MSw qf(0.95, m-1, N-m) ######################################## ######################################## OnewayAOV<-function(tabl, sig) { N=length(tabl) n=dim(tabl)[1] m=dim(tabl)[2] corterm=(sum(tabl)^2)/N sst=sum(tabl^2)-corterm a=apply(tabl, 2, sum) ssa=sum( (a^2)/n)-corterm ssw=sst-ssa doft = N-1 dofa = (m-1) dofw = (N-m) MSa=ssa/dofa MSw=ssw/dofw Ft=MSa/MSw Ftable=qf(sig, m-1, N-m) print("*********************************************", quote=FALSE); print("Source SS DoF MeanSQ TEST", quote=FALSE); print(paste(sep=" ", "Among ", format(ssa) , format(dofa), format(MSa), format(Ft) ) , quote=FALSE) print(paste(sep=" ", "Within ", format(ssw), format(dofw), format(MSw) ) , quote=FALSE) print(paste(sep=" ", "Total ", format(sst), format(doft) ) , quote=FALSE) print("*********************************************", quote=FALSE); print(paste(sep=" ","F=", format(Ftable), " at ",sig," level, with ", dofa, " and ", dofw ,"degrees of freedom"), quote=FALSE) print("THUS:", quote=FALSE) if(Ft>Ftable) { print("reject Null hypothesis", quote=FALSE) } else { print("failure to reject Null hypothesis", quote=FALSE) } list(MSa=MSa, MSw=MSw, SSt=sst) } ################################ ################################ OnewayAOV(tab2.16, 0.99) TwowayAOV<-function(tabl) { N=length(tabl) n=dim(tabl)[1] m=dim(tabl)[2] corterm=(sum(tabl)^2)/N sst=sum(tabl^2)-corterm a=apply(tabl, 2, sum) b=apply(tabl, 1, sum) ssa=sum( (a^2)/n)-corterm ssb=sum( (b^2)/m)-corterm sse=sst-(ssa+ssb) MSa=ssa/(m-1) MSb=ssw/(n-1) MSe=sse/( (m-1)*(n-1)) Fta=MSa/Mse Ftb=MSb/Mse ### etc.....compare F with F-tables } ###################################### PEN = c(89, 88, 97, 94, 84, 77, 92, 79, 81, 87, 87, 85, 87, 92, 89, 84, 79, 81, 80, 88) pen = matrix(PEN, byrow=TRUE, ncol=4) N=length(pen) n=dim(pen)[1] m=dim(pen)[2] ga = sum(pen)/N GA = matrix(rep(ga, length=N), ncol=4) blox = matrix(rep(apply(pen-GA, 1, mean), m), byrow=FALSE, ncol=4) treatz = matrix(rep(apply(pen-GA, 2, mean), n), byrow=TRUE, ncol=4) res = pen-(GA+blox+treatz) c(sum(pen^2), sum(GA^2) , sum(blox^2), sum(treatz^2) ,sum(res^2)) 148480 147920 264 70 226 dof = c(20, 1, 4, 3, 12) 89 88 97 94 84 77 92 79 81 87 87 85 87 92 89 84 79 81 80 88 86 86 86 86 86 86 86 86 86 86 86 86 86 86 86 86 86 86 86 86 6 6 6 6 -3 -3 -3 -3 -1 -1 -1 -1 2 2 2 2 -4 -4 -4 -4 -2 -1 3 0 -2 -1 3 0 -2 -1 3 0 -2 -1 3 0 -2 -1 3 0 -1 -3 2 2 3 -5 6 -4 -2 3 -1 0 1 5 -2 -4 -1 0 -5 6 ########################################### ########################################### ########################################### d1 = rnorm(n= 15, m=92.1, sd=4) d2 = rnorm(n= 15, m=92.3, sd=4) d3 = rnorm(n= 15, m=93.1, sd=4) d4 = rnorm(n= 15, m=90.5, sd=4) boxplot(list(d1, d2, d3, d4)) D = cbind(d1,d2,d3,d4) save(D, file="D:/LEES/CLASSES/GEOL_520 Data_Analysis/LAB_NOTES", ascii=TRUE) write.table(formatC(D, width=4, digits=4), file="D:/LEES/CLASSES/GEOL_520 Data_Analysis/LAB_NOTES/volcdates.txt", row.names = FALSE, col.names = FALSE, quote=FALSE) X = matrix(scan(file="D:/LEES/CLASSES/GEOL_520 Data_Analysis/LAB_NOTES/volcdates.txt"), byrow=TRUE, ncol=4) boxplot(X[,1], X[,2],X[,3],X[,4]) OnewayAOV(D, 0.95) OnewayAOV(X, 0.95) OnewayAOV(X, 0.90)