R/ancmg1.R

ancmg1 <-
function(x,y,pool=TRUE,jcen=1,fr=1,depfun=fdepth,nmin=8,op=3,tr=.2,
SEED=TRUE,pr=TRUE,pts=NA,con=0,nboot=NA,alpha=.05,bhop=FALSE){
#
# ANCOVA
# for two or more groups based on trimmed means or medians
# Single  covariate is assumed.
#
# FOR TWO OR MORE COVARIATES, USE ANCMG
#
# for two groups and one covariate, also consider ancova and ancGpar
#
# op=1 use omnibus test for trimmed means, with trimming given by tr
# op=2 use omnibus test for medians.
#       (Not recommended when there are tied values, use op=4)
# op=3 multiple comparisons using trimming and percentile bootstrap.
#     This method seems best for general use.
# op=4 multiple comparisons using medians and percentile bootstrap
#
# y is matrix with J columns, so have J groups.
# or y can have list mode with length J
#
# x is a matrix with Jp columns, so first p columns
# correspond to the p covariates in the first group, etc.
#
# Or,
# x can have list mode with length J
#
# nmin is the minimum sample size allowed for any group
# when testing hypotheses.
# If a design point results in a sample size <nmin,
# that point is eliminated.
#
#  pool=T means pool the data when determining the center of the
#  design points and the measure of scatter when applying the smooth
#
#  pool=F, does not pool but rather use the data from group
#  jcen to determine center and the measure of scatter
#
#  The output includes a matrix of sample sizes. The ith row
#  corresponds to the ith point used to compare groups.
#  The jth column indicates the number of points (the sample size)
#  that was found for the jth group. That is, how many points
#  in the jth group were found that are close to the design point
#  under consideration.
#
if(SEED)set.seed(2) # set the seed so that MVE always gives same result
if(pr){
if(op==1)print("Trimmed means are to be compared. For medians, use op=2")
if(op==2)print("Medians are to be compared. For trimmed means, use op=1")
if(op==3)print("20% trimmed means are compared. For medians, use op=4")
if(op==4)print("medians are compared. For 20% trimmed means, use op=3")
}
output<-NULL
conout=NULL
nval<-NA
if(is.matrix(y))J<-ncol(y)
if(is.matrix(x))pval=ncol(x)
if(is.list(y))J<-length(y)
if(is.list(x))pval<-ncol(as.matrix(x[[1]]))
if(pval>1)stop("More than one covariate. Use ancmg")
if(J==1)stop("Only have one group stored in y")
if(is.matrix(x)){
if(ncol(x)%%J!=0)stop("Number of columns of x should be a multiple of ncol(y)")
}
if(is.matrix(x)){
xcen<-x[,jcen]
}
if(is.list(x))xcen<-x[[jcen]]
if(is.na(pts[1])){
if(pool){
if(is.matrix(x))xval<-stackit(x,1)
if(is.list(x))xval<-stacklist(x)
temp<-idealf(xval)
pts<-temp$ql
pts[2]<-median(xval)
pts[3]<-temp$qu
}
if(!pool){
temp<-idealf(xcen)
pts<-temp$ql
pts[2]<-median(xval)
pts[3]<-temp$qu
}}
nval<-matrix(NA,ncol=J,nrow=length(pts))
for(j in 1:J){
for(i in 1:length(pts)){
if(is.matrix(x)  && is.matrix(y)){
nval[i,j]<-length(y[near(x[,j],pts[i],fr=fr)])
}
if(is.matrix(x)  && is.list(y)){
tempy<-y[[j]]
nval[i,j]<-length(tempy[near(x[,j],pts[i],fr=fr)])
}
if(is.list(x)  && is.matrix(y)){
xm<-as.matrix(x[[j]])
nval[i,j]<-length(y[near(xm,pts[i],fr=fr),j])
}
if(is.list(x)  && is.list(y)){
tempy<-y[[j]]
xm<-as.matrix(x[[j]])
nval[i,j]<-length(tempy[near(xm,pts[i],fr=fr)])
}
#
}}
flag<-rep(TRUE,length(pts))
for(i in 1:length(pts)){
if(min(nval[i,])<nmin)flag[i]<-F
}
nflag<-F
if(sum(flag)==0){
print("Warning: No design points found with large enough sample size")
nflag<-T
}
if(!nflag){
pts<-pts[flag] # eliminate points for which the sample size is too small
nval<-nval[flag,]
if(!is.matrix(pts))pts<-t(as.matrix(pts))
output<-matrix(NA,nrow=length(pts),ncol=3)
dimnames(output)<-list(NULL,c("point","test.stat","p-value"))
if(op==3 || op==4)output<-list()
}
for(i in 1:length(pts)){
if(op==1 || op==2)output[i,1]<-i
icl<-0-pval+1
icu<-0
yval<-list()
for(j in 1:J){
if(is.matrix(x)  && is.matrix(y)){
yval[[j]]<-y[near(x[,j],pts[i],fr=fr),j]
}
if(is.matrix(x)  && is.list(y)){
tempy<-y[[j]]
yval[[j]]<-tempy[near(x[,j],pts[i],fr=fr)]
}
if(is.list(x)  && is.matrix(y)){
yval[[j]]<-y[near3d(x[[j]],pts[i],fr=fr),j]
}
if(is.list(x)  && is.list(y)){
tempy<-y[[j]]
yval[[j]]<-tempy[near(x[[j]],pts[i],fr=fr)]
}
#
}
if(op==1)temp<-t1way(yval,tr=tr)
if(op==2)temp<-med1way(yval,SEED=SEED,pr=FALSE)

if(op==1 || op==2){
output[i,2]<-temp$TEST
output[i,3]<-temp$p.value
}
if(op==3){
output[[i]]<-linconpb(yval,alpha=alpha,SEED=SEED,con=con,bhop=bhop,est=tmean,nboot=nboot)
}
if(op==4){output[[i]]<-medpb(yval,alpha=alpha,SEED=SEED,con=con,bhop=bhop,nboot=nboot)
}
}
if(op==1 || op==2)tempout=output
if(nflag)output<-NULL
if(op==3 || op==4){
conout=list()
tempout=list()
for(j in 1:length(output))tempout[[j]]=output[[j]]$output
for(j in 1:length(output))conout[[j]]=output[[j]]$con
}
#list(points.chosen=pts,sample.sizes=nval,output=tempout,contrast.coef=conout[[1]])
list(points.chosen=pts,sample.sizes=nval,point=tempout,contrast.coef=conout[[1]])
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.