R/Scheffe.R

Scheffe <-
function(x,con=0,alpha=.05,WARN=TRUE){
#
#  Scheffe's MCP
#
#  The data are assumed to be stored in $x$ in list mode, a matrix
#  or a data frame. If in list mode,
#  length(x) is assumed to correspond to the total number of groups.
#  It is assumed all groups are independent.
#
#  con is a J by d matrix containing the contrast coefficients that are used.
#  If con is not specified, all pairwise comparisons are made.
#
#  Missing values are automatically removed.
#
#
if(WARN)print('WARNING: Suggest using lincon instead')
if(is.data.frame(x))x=as.matrix(x)
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop('Data must be stored in a matrix or in list mode.')
con<-as.matrix(con)
J<-length(x)
n=NA
xbar<-NA
for(j in 1:J){
xx<-!is.na(x[[j]])
val<-x[[j]]
x[[j]]<-val[xx]  # Remove missing values
n[j]=length(x[[j]])
xbar[j]<-mean(x[[j]])
}
N=sum(n)
df2=N-J
AOV=anova1(x)
if(sum(con^2)==0){
CC<-(J^2-J)/2
df1=J-1
psihat<-matrix(0,CC,6)
dimnames(psihat)<-list(NULL,c('Group','Group','psihat','ci.lower','ci.upper',
'p.value'))
test<-matrix(NA,CC,5)
dimnames(test)<-list(NULL,c('Group','Group','test','crit','se'))
jcom<-0
for (j in 1:J){
for (k in 1:J){
if (j < k){
jcom<-jcom+1
test[jcom,3]<-abs(xbar[j]-xbar[k])/sqrt((J-1)*AOV$MSWG*(1/n[j]+1/n[k]))
sejk<-sqrt((CC-1)*AOV$MSWG*(1/n[j]+1/n[k]))
test[jcom,5]<-sejk
psihat[jcom,1]<-j
psihat[jcom,2]<-k
test[jcom,1]<-j
test[jcom,2]<-k
psihat[jcom,3]<-(xbar[j]-xbar[k])
psihat[jcom,6]<-1-pf(test[jcom,3],df1,df2)
crit=sqrt(qf(1-alpha,df1,df2))
test[jcom,4]<-crit
psihat[jcom,4]<-(xbar[j]-xbar[k])-crit*sejk
psihat[jcom,5]<-(xbar[j]-xbar[k])+crit*sejk
}}}}
if(sum(con^2)>0){
if(nrow(con)!=length(x)){
stop('The number of groups does not match the number of contrast coefficients.')
}
psihat<-matrix(0,ncol(con),5)
dimnames(psihat)<-list(NULL,c('con.num','psihat','ci.lower','ci.upper',
'p.value'))
test<-matrix(0,ncol(con),4)
dimnames(test)<-list(NULL,c('con.num','test','crit','se'))
df1<-nrow(con)-1
df2=N-nrow(con)
crit=sqrt(qf(1-alpha,df1,df2))
for (d in 1:ncol(con)){
psihat[d,1]<-d
psihat[d,2]<-sum(con[,d]*xbar)
sejk<-sqrt(df1*AOV$MSWG*sum(con[,d]^2/n))
test[d,1]<-d
test[d,2]<-sum(con[,d]*xbar)/sejk
test[d,3]<-crit
test[d,4]<-sejk
psihat[d,3]<-psihat[d,2]-crit*sejk
psihat[d,4]<-psihat[d,2]+crit*sejk
psihat[d,5]<-(1-pf(abs(test[d,2]),df1,df2))
}
}
list(n=n,test=test,psihat=psihat)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.