med1way <-
function(x,grp=NA,alpha=.05,crit=NA,iter=1000,SEED=TRUE,pr=T){
#
# A heteroscedastic one-way ANOVA for medians.
#
# If
# crit=NA, an appropriate critical value is determined
# for the alpha value specified and
# a p-value is returned.
#
# If a value for crit is specified, it is used as the critical
# value, but no p-value is reported. Specifying a value for
# crit reduces execution time.
# With crit=NA, the critical value is a function of the sample
# sizes and is determined by calling the function med1way.crit.
#
# The data are assumed to be stored in $x$ in list mode.
# Length(x) is assumed to correspond to the total number of groups.
# By default, the null hypothesis is that all groups have a common mean.
# To compare a subset of the groups, use grp to indicate which
# groups are to be compared. For example, if you type the
# command grp<-c(1,3,4), and then execute this function, groups
# 1, 3, and 4 will be compared with the remaining groups ignored.
#
# Missing values are automatically removed.
#
if(is.data.frame(x))x=as.matrix(x)
if(is.matrix(x)){
y<-list()
for(j in 1:ncol(x))y[[j]]<-x[,j]
x<-y
}
if(is.na(sum(grp[1])))grp<-c(1:length(x))
if(!is.list(x))stop("Data are not stored in a matrix or in list mode.")
J<-length(grp) # The number of groups to be compared
n<-vector("numeric",J)
w<-vector("numeric",J)
xbar<-vector("numeric",J)
for(j in 1:J){
xx<-!is.na(x[[j]])
val<-x[[j]]
x[[j]]<-val[xx] # Remove missing values
w[j]<-1/msmedse(x[[grp[j]]])^2
xbar[j]<-median(x[[grp[j]]])
n[j]<-length(x[[grp[j]]])
}
pval<-NA
u<-sum(w)
xtil<-sum(w*xbar)/u
TEST<-sum(w*(xbar-xtil)^2)/(J-1)
if(is.na(crit)){
temp<-med1way.crit(n,alpha,SEED=SEED,iter=iter,TEST=TEST)
crit.val<-temp$crit.val
}
if(!is.na(crit))crit.val<-crit
list(TEST=TEST,crit.val=crit.val,p.value=temp$p.value)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.