R/msmedsub.R

msmedsub <-
function(n,con=0,alpha=.05,se.fun=bpmedse,iter=1000,SEED=TRUE){
#
# Determine a Studentized critical value, assuming normality
# and homoscedasticity, for the function msmedv2
#
# Goal: 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.
#
if(SEED)set.seed(2)
con<-as.matrix(con)
J<-length(n)
h<-vector("numeric",J)
w<-vector("numeric",J)
xbar<-vector("numeric",J)
x<-list()
test<-NA
testmax<-NA
for (it in 1:iter){
for(j in 1:J){
x[[j]]<-rnorm(n[j])
xbar[j]<-median(x[[j]])
 w[j]<-se.fun(x[[j]])^2
}
if(sum(con^2!=0))CC<-ncol(con)
if(sum(con^2)==0){
CC<-(J^2-J)/2
jcom<-0
for (j in 1:J){
for (k in 1:J){
if (j < k){
jcom<-jcom+1
test[jcom]<-abs(xbar[j]-xbar[k])/sqrt(w[j]+w[k])
}}}}
if(sum(con^2)>0){
for (d in 1:ncol(con)){
sejk<-sqrt(sum(con[,d]^2*w))
test[d]<-sum(con[,d]*xbar)/sejk
}}
testmax[it]<-max(abs(test))
}
testmax<-sort(testmax)
testmax
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.