### #Ryan's Q in R #from Day and Quinn Ecological Monographs 1989 #pages 460-461 in appendix 1 #and a few extra pieces from me # Jarrett Byrnes - byrnes@msi.ucsb.edu #last updated 11/20/08 ### #TODO # - calculate sec for unequal variances # - make this a general framework for more tests in Day & Quinn tests (perhaps a library) # - pretty printing of results (groupings a,b,c, etc...) # - integrate with plotting? # - correction for ANCOVA ### #some useful stuff for calculating number of combinations ## #n = number of unique entries, k=length of combination ncombos<-function(n,k) factorial(n)/(factorial(n-k)*factorial(k)) #how many treatments were there, given the number of combinations bsolve.ncombos<-function(x) (1+sqrt(1+8*x))/2 ## #dif levels ## order.diff<-function(x){ ret<-c() for (i in 1:x) ret<-c(ret,rep(i,x-i)) return(ret) } ### #Kramer's modification for unequal sample sizes #If sample sizes are equal, will return the same value #as if uncorrected ## sec.unequal.n<-function(object, n1, n2){ msd<-sum(object$residuals^2)/object$df.residual sec<-sqrt(msd*(1/n1 + 1/n2)) return(sec) } ### #calculate Ryan's Q value #a= chosen sig level #b = adj sig level #nmeans = # of means in group to be tested (e.g. 4, then 3, then 2...) #m = number of means in the experiment ### Q.ryans<-function(nmeans, m, df, a=0.05){ b<-1-((1-a)^(nmeans/m)) Q<-qtukey(b, nmeans, df, lower.tail=F) return(Q) } ### #calculate a ryan's Q table #the output isn't pretty, but, it does tell you which #differences are significantly different from 0 at the specified alpha value ### ryans_q<-function(obj, a=0.05){ #get a table of means and sample sizes mm<-model.tables(obj, "means") tabs<-mm$tables[-1] reps<-mm$n #DF of the standard error. #Again, assumes equal sample size and variance df<-obj$df.residual #our ryan's Q object - a list of matrices rq.obj<-list() for (trt in names(tabs)){ tab<-tabs[[trt]] rep<-reps[[trt]] #gakked from TukeyHSD code levs <- if(length(d <- dim(tab)) > 1) { dn <- dimnames(tab) apply(do.call("expand.grid", dn), 1, paste, collapse=":") } else names(tab) means<-as.vector(tab) n<-as.vector(rep) if(length(n) CV) issig<-1 tmp.mat[i,4]<-issig } rq.obj[[trt]]<-tmp.mat } return(rq.obj) } ###DEMO using data from Day & Quinn #data in http://homes.msi.ucsb.edu/~byrnes/r_files/dq_barnacle.csv #barnacles<-read.csv("dq_barnacle.csv") #banova<-aov(Barnacles ~ Treatment, data=barnacles) # #ryans_q(banova) #TukeyHSD(banova) ###