R/rananova.R

rananova <-
function(x,tr=.2,grp=NA){
#
#  A heteroscedastic one-way random effects ANOVA for trimmed means.
#
#  The data are assumed to be stored in a matrix on in list mode.
#  If in list mode,
#  Length(x) is assumed to correspond to the total number of groups.
#  If the data are stored in a matrix, groups correspond to columns.
#  By default, the null hypothesis is that all group 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.
#
if(is.matrix(x))x<-listm(x)
if(is.na(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
#if(pr)print("The number of groups to be compared is")
#print(J)
h<-1
xbar<-1
ybar<-1
wvar<-1
ell<-0
for(j in 1:J){
ell[j]<-length(x[[grp[j]]])/(length(x[[grp[j]]])+1)
h[j]<-length(x[[grp[j]]])-2*floor(tr*length(x[[grp[j]]]))
   # h is the number of observations in the jth group after trimming.
ybar[j]<-winmean(x[[grp[j]]],tr)
xbar[j]<-mean(x[[grp[j]]],tr)
wvar[j]<-winvar(x[[grp[j]]],tr)
}
q<-NA
bsst<-var(xbar)
for (j in 1:J)q[j]<-(length(x[[grp[j]]]-1)-1)*wvar[j]/(h[j]*(h[j]-1))
wssw<-mean(q)
D<-bsst/wssw
g<-q/J
nu1<-((J-1)*sum(q))^2/((sum(q))^2+(J-2)*J*sum(q^2))
nu2<-(sum(J*q))^2/sum((J*q)^2/(h-1))
sig<-1-pf(D,nu1,nu2)
# Next, estimate the Winsorized intraclass correlation
sighat<-mean(ell*(ybar-(sum(ell*ybar)/sum(ell)))^2)
rho<-sighat/(sighat+winmean(wvar,tr))
list(teststat=D,df=c(nu1,nu2),p.value=sig,rho=rho,num.groups=J)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.