R/reg1mcp.R

reg1mcp <-
function(x,y,regfun=tsreg,SEED=TRUE,nboot=100,xout=FALSE,outfun=outpro,STAND=TRUE,alpha=.05,pr=TRUE,MC=FALSE,...){
#
#  Perform all pairwise comparisons of intercepts among J independent groups
#  Do the same of the first slope, followed by the 2nd slope, etc.
#  
#  Control FWE via Hochberg's methods for each set of 
#  (J^2-J)/2 parameters. That is, control FWE for the intercepts
#  Do the same for the first slope, etc. 
#
#  #  x and y are assumed to have list mode having 
#  length J equal to the number of groups
#  For example, x[[1]] and y[[1]] contain the data for group 1.
#
#   xout=T will eliminate leverage points using the function outfun, 
#   which defaults to the projection method.
#
#  OUTPUT: 
#   n is sample size after missing values are removed
#   nv.keep is sample size after leverage points are removed.
#   output contains all pairwise comparisons  
#   For each parameter, FWE is controlled using Hochberg's method
#   So by default, for the intercepts, 
#   all pairwise comparisons are performed with FWE=.05
#   For the first slope, all pairwise comparisons are performed 
#   with FWE=.05, etc.
#
if(SEED)set.seed(2)
if(!is.list(x))stop('Argument x should have list mode')
if(!is.list(y))stop('Argument y should have list mode')
J=length(x) # number of groups
x=lapply(x,as.matrix)
pchk=lapply(x,ncol)
temp=matl(pchk)
if(var(as.vector(temp))!=0)stop('Something is wrong. Number of covariates differs among the groups being compared')
nv=NULL
p=ncol(x[[1]])
p1=p+1
for(j in 1:J){
xy=elimna(cbind(x[[j]],y[[j]]))
x[[j]]=xy[,1:p]
y[[j]]=xy[,p1]
x[[j]]=as.matrix(x[[j]])
nv=c(nv,nrow(x[[j]]))
}
nv.keep=nv
critrad=NULL
if(xout){
temp=lapply(x,outfun,plotit=FALSE,STAND=STAND,...)
for(j in 1:J){
x[[j]]=x[[j]][temp[[j]]$keep,]
y[[j]]=y[[j]][temp[[j]]$keep]
nv.keep[j]=length(y[[j]])
}}
tot=(J^2-J)/2
dvec<-alpha/c(1:tot)
outl=list()
nr=tot*p1
outp=matrix(NA,ncol=7,nrow=nr)
x=lapply(x,as.matrix)
rlab=rep('Intercept',tot)
xx=list()
yy=list()
iall=0
ivp=c(1,tot)-tot
for(ip in 1:p){
#iv=ip-1
rlab=c(rlab,rep(paste('slope',ip),tot))
}
i=0
sk=1+tot*p
st=seq(1,sk,tot)
st=st-1
for(j in 1:J){
for(k in 1:J){
if(j < k){
i=i+1
st=st+1
xx[[1]]=x[[j]][,1:p]
xx[[2]]=x[[k]][,1:p]
yy[[1]]=y[[j]]
yy[[2]]=y[[k]]
if(!MC)temp=reg2ci(xx[[1]],yy[[1]],xx[[2]],yy[[2]],regfun=regfun)$output
if(MC)temp=reg2ci(xx[[1]],yy[[1]],xx[[2]],yy[[2]],regfun=regfun)$output
iall=iall+1
outp[iall,1]=j
outp[iall,2]=k
outp[st,3]=temp[,4]
outp[st,5]=temp[,2]
outp[st,6]=temp[,3]
}}}
for(i in 1:p1){
ivp=ivp+tot
temp2<-order(0-outp[ivp[1]:ivp[2],3])
icc=c(ivp[1]:ivp[2])
icc[temp2]=dvec
outp[ivp[1]:ivp[2],4]=icc
}
flag=(outp[,3]<=outp[,4])                                   
outp[,7]=rep(0,nr)
outp[flag,7]=1   
v=outp[1:tot,1]
vall=rep(v,p1)
outp[,1]=vall
v=outp[1:tot,2]
vall=rep(v,p1)
outp[,2]=vall
dimnames(outp)=list(rlab,c('Group','Group','p.value','p.crit','ci.low','ci.hi','Sig'))
list(n=nv,n.keep=nv.keep,output=outp)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.