R/qdmcpdif.R

qdmcpdif <-
function(x, con = 0,alpha = 0.05){
#
# MCP with medians on difference scores
# FWE controlled with Rom's method
#
if(is.data.frame(x))x=as.matrix(x)
if(!is.matrix(x))x<-matl(x)
if(!is.matrix(x))stop("Data must be stored in a matrix or in list mode.")
con<-as.matrix(con)
J<-ncol(x)
xbar<-vector("numeric",J)
x<-elimna(x)  # Remove missing values
nval<-nrow(x)
h1<-nrow(x)
df<-h1-1
if(sum(con^2!=0))CC<-ncol(con)
if(sum(con^2)==0)CC<-(J^2-J)/2
ncon<-CC
if(alpha==.05){
dvec<-c(.05,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncon > 10){
avec<-.05/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.01,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(ncon > 10){
avec<-.01/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha != .05 && alpha != .01)dvec<-alpha/c(1:ncon)
if(sum(con^2)==0){
psihat<-matrix(0,CC,5)
dimnames(psihat)<-list(NULL,c("Group","Group","psihat","ci.lower","ci.upper"))
test<-matrix(NA,CC,5)
dimnames(test)<-list(NULL,c("Group","Group","p-value","p.crit","se"))
temp1<-0
jcom<-0
for (j in 1:J){
for (k in 1:J){
if (j < k){
jcom<-jcom+1
dv<-x[,j]-x[,k]
test[jcom,5]<-msmedse(dv)
temp<-sintv2(dv,alpha=alpha/CC)
temp1[jcom]<-temp$p.value
test[jcom,3]<-temp$p.value
psihat[jcom,1]<-j
psihat[jcom,2]<-k
test[jcom,1]<-j
test[jcom,2]<-k
psihat[jcom,3]<-median(dv)
psihat[jcom,4]<-temp$ci.low
psihat[jcom,5]<-temp$ci.up
}}}
temp2<-order(0-temp1)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2]>=zvec)
if(sum(sigvec)<ncon){
dd<-ncon-sum(sigvec) #number that are sig.
ddd<-sum(sigvec)+1
zvec[ddd:ncon]<-dvec[ddd]
}
test[temp2,4]<-zvec
}
if(sum(con^2)>0){
if(nrow(con)!=ncol(x))print("WARNING: The number of groups does not match the number of contrast coefficients.")
ncon<-ncol(con)
psihat<-matrix(0,ncol(con),4)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper"))
test<-matrix(0,ncol(con),4)
dimnames(test)<-list(NULL,c("con.num","sig","crit.sig","se"))
temp1<-NA
for (d in 1:ncol(con)){
psihat[d,1]<-d
for(j in 1:J){
if(j==1)dval<-con[j,d]*x[,j]
if(j>1)dval<-dval+con[j,d]*x[,j]
}
temp3<-sintv2(dval)
temp1[d]<-temp3$p.value
test[d,1]<-d
test[d,4]<-msmedse(dval)
psihat[d,2]<-median(dval)
psihat[d,3]<-temp3$ci.low
psihat[d,4]<-temp3$ci.up
}
test[,2]<-temp1
temp2<-order(0-temp1)
zvec<-dvec[1:ncon]
print(c(ncon,zvec))
sigvec<-(test[temp2,2]>=zvec)
if(sum(sigvec)<ncon){
dd<-ncon-sum(sigvec) #number that are sig.
ddd<-sum(sigvec)+1
zvec[ddd:ncon]<-dvec[ddd]
}
test[temp2,3]<-zvec
}
if(sum(con^2)==0)num.sig<-sum(psihat[,4]>0)+ sum(psihat[,5]<0)
if(sum(con^2)>0)num.sig<-sum(psihat[,3]>0)+ sum(psihat[,4]<0)
list(test=test,psihat=psihat,con=con,num.sig=num.sig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.