R/tkmcp.R

tkmcp <-
function(x,alpha=.05,ind.pval=T){
#
# conventional Tukey-Kramer multiple comparison procedure
# for all pairiwise comparisons.
#
#  ind.pval=T, computes p-value for each individual test
#  ind.pval=F computes p-value based on controlling the
#  familywise error rate. (The alpha level at which the
#  Tukey-Kramer test would reject.)
#
if(is.matrix(x))x<-listm(x)
J<-length(x)
A<-0
B<-0
C<-0
N<-0
for(j in 1:J){
N<-N+length(x[[j]])
A<-A+sum(x[[j]]^2)
B<-B+sum(x[[j]])
C<-C+(sum(x[[j]]))^2/length(x[[j]])
}
SST<-A-B^2/N
SSBG<-C-B^2/N
SSWG<-A-C
nu1<-length(x)-1
nu2<-N-length(x)
MSBG<-SSBG/nu1
MSWG<-SSWG/nu2
numcom<-length(x)*(length(x)-1)/2
output<-matrix(nrow=numcom,ncol=7)
dimnames(output)<-list(NULL,c("Group","Group","t.test","est.difference",
"ci.lower","ci.upper","p.value"))
ic<-0
for (j in 1:J){
for (k in 1:J){
if (j < k){
ic<-ic+1
output[ic,1]<-j
output[ic,2]<-k
dif<-mean(x[[j]])-mean(x[[k]])
output[ic,3]<-abs(dif)/sqrt(MSWG*(1/length(x[[j]])+1/length(x[[k]]))/2)
output[ic,4]<-dif
crit<-qtukey(1-alpha,length(x),nu2)
output[ic,5]<-dif-crit*sqrt(MSWG*(1/length(x[[j]])+1/length(x[[k]]))/2)
output[ic,6]<-dif+crit*sqrt(MSWG*(1/length(x[[j]])+1/length(x[[k]]))/2)
if(!ind.pval)output[ic,7]<-1-ptukey(output[ic,3],length(x),nu2)
if(ind.pval)output[ic,7]<-2*(1-pt(output[ic,3],nu2))
}}}
output
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.