R/msmed.R

msmed <-
function(x,y=NA,con=0,alpha=.05){
#
# Test a set of linear contrasts using Medians
#
#  The data are assumed to be stored in $x$ in a matrix or in list mode.
#  Length(x) is assumed to correspond to the total number of groups, J
#  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(!is.na(y[1])){
xx<-list()
xx[[1]]<-x
xx[[2]]<-y
if(is.matrix(x) || is.list(x))stop("When y is speficied, x should not have list mode or be a matrix")
x<-xx
}
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)
h<-vector("numeric",J)
w<-vector("numeric",J)
xbar<-vector("numeric",J)
for(j in 1:J){
xx<-!is.na(x[[j]])
val<-x[[j]]
if(sum(duplicated(val)>0)){
print(paste("Warning: Group",j, "has tied values. Might want to used medpb"))
}
x[[j]]<-val[xx]  # Remove missing values
xbar[j]<-median(x[[j]])
w[j]<-msmedse(x[[j]])^2 # Squared standard error.
}
if(sum(con^2!=0))CC<-ncol(con)
if(sum(con^2)==0){
CC<-(J^2-J)/2
psihat<-matrix(0,CC,5)
dimnames(psihat)<-list(NULL,c("Group","Group","psihat","ci.lower","ci.upper"))
test<-matrix(NA,CC,6)
dimnames(test)<-list(NULL,c("Group","Group","test","crit","se","p.value"))
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(w[j]+w[k])
test[jcom,6]<-2*(1-pt(test[jcom,3],999))
sejk<-sqrt(w[j]+w[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])
crit<-NA
if(CC==1)crit<-qnorm(1-alpha/2)
if(CC>1){
if(alpha==.05)crit<-smmcrit(500,CC)
if(alpha==.01)crit<-smmcrit01(500,CC)
if(is.na(crit))warning("Can only be used with alpha=.05 or .01")
}
test[jcom,4]<-crit
psihat[jcom,4]<-psihat[jcom,3]-crit*test[jcom,5]
psihat[jcom,5]<-psihat[jcom,3]+crit*test[jcom,5]
}}}}
if(sum(con^2)>0){
if(nrow(con)!=length(x))warning("The number of groups does not match the number of contrast coefficients.")
psihat<-matrix(0,ncol(con),4)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper"))
test<-matrix(0,ncol(con),5)
dimnames(test)<-list(NULL,c("con.num","test","crit","se","p.value"))
for (d in 1:ncol(con)){
psihat[d,1]<-d
psihat[d,2]<-sum(con[,d]*xbar)
sejk<-sqrt(sum(con[,d]^2*w))
test[d,1]<-d
test[d,2]<-sum(con[,d]*xbar)/sejk
test[d,5]<-2*(1-pt(abs(test[d,2]),999))
crit<-NA
if(CC==1)crit<-qnorm(1-alpha/2)
if(alpha==.05)crit<-smmcrit(500,ncol(con))
if(alpha==.01)crit<-smmcrit01(500,ncol(con))
test[d,3]<-crit
test[d,4]<-sejk
psihat[d,3]<-psihat[d,2]-crit*sejk
psihat[d,4]<-psihat[d,2]+crit*sejk
}}
list(test=test,psihat=psihat)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.