R/t1wayF.R

t1wayF <-
function(x,fac,tr=.2,nboot=100,SEED=TRUE){
#
# Same a t1way, but now the data are assumed to be
# stored in a matrix or data frame where one of the columns contain
# the data to be analyzed and another column contains the group
# identification.
#
# For example, if dat is a data frame, with column 1 containing
# the outcome measures of interest, and column 2 is a factor variable
# indicating to  which group a value in column 1 belongs, then
#  t1wayF(dat[,1],dat[,2])
#  will test the hypothesis that all J groups have identical
#  trimmed means.
#
#  Missing values are automatically removed.
#
library(MASS)
if(SEED)set.seed(2)
x=fac2list(x,fac)
J<-length(x)
h<-vector("numeric",J)
w<-vector("numeric",J)
xbar<-vector("numeric",J)
pts=NULL
nval=0
for(j in 1:J)x[[j]]=elimna(x[[j]])
for(j in 1:J){
val<-x[[j]]
val<-elimna(val)
nval[j]=length(val)
pts=c(pts,val)
x[[j]]<-val # missing values have been removed
h[j]<-length(x[[j]])-2*floor(tr*length(x[[j]]))
   # h is the number of observations in the jth group after trimming.
w[j]<-h[j]*(h[j]-1)/((length(x[[j]])-1)*winvar(x[[j]],tr))
xbar[j]<-mean(x[[j]],tr)
}
u<-sum(w)
xtil<-sum(w*xbar)/u
A<-sum(w*(xbar-xtil)^2)/(J-1)
B<-2*(J-2)*sum((1-w/u)^2/(h-1))/(J^2-1)
TEST<-A/(B+1)
nu1<-J-1
nu2<-1./(3*sum((1-w/u)^2/(h-1))/(J^2-1))
sig<-1-pf(TEST,nu1,nu2)
#
# Determine explanatory effect size
#
chkn=var(nval)
if(chkn==0){
top=var(xbar)
cterm=NULL
if(tr==0)cterm=1
if(tr==0.2)cterm=.642
if(is.null(cterm))cterm=area(dnormvar,qnorm(tr),qnorm(1-tr))+
2*(qnorm(tr)^2)*tr
bot=winvar(pts,tr=tr)/cterm
e.pow=top/bot
}
if(chkn!=0){
vals=0
N=min(nval)
xdat=list()
for(i in 1:nboot){
for(j in 1:J){
xdat[[j]]=sample(x[[j]],N)
vals[i]=t1way.effect(xdat,tr=tr)$Var.Explained
}}
e.pow=mean(vals,na.rm=TRUE)
}
list(TEST=TEST,nu1=nu1,nu2=nu2,p.value=sig,Var.Explained=e.pow,
Effect.Size=sqrt(e.pow))
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.