R/olsLmcp.R

olsLmcp <-
function(x,y,xout=TRUE,outfun=outpro,ISO=FALSE,STAND=TRUE,alpha=.05,pr=TRUE,BLO=FALSE,HC3=FALSE,...){
#
#  All pairwise comparison of regression models among J independent groups
#  That is, for groups j and k, all j<k, test H_0: all corresponding
#  parameters are equal
#
#  Strategy: Use HC4  estimate of standard errors followed by
#  Johansen type test statistic.
#
#  Hochberg to control FWE
#
#  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
#
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=5,nrow=tot)
x=lapply(x,as.matrix)
xx=list()
yy=list()
iall=0
ivp=c(1,tot)-tot
i=0
for(j in 1:J){
for(k in 1:J){
if(j < k){
i=i+1
xx[[1]]=x[[j]]
xx[[2]]=x[[k]]
yy[[1]]=y[[j]]
yy[[2]]=y[[k]]
if(!ISO)all=ols1way(xx,yy,HC3=HC3)
if(ISO)all=ols1wayISO(xx,yy,HC3=HC3)
temp=all$adjusted.p.value
if(is.null(temp))temp=all$p.value
outp[i,1]=j
outp[i,2]=k
outp[i,3]=temp
}}
temp2<-order(0-outp[,3])
icc=c(1:tot)
icc[temp2]=dvec
outp[,4]=icc
}
flag=(outp[,3]<=outp[,4])
outp[,5]=rep(0,tot)
outp[flag,5]=1
dimnames(outp)=list(NULL,c('Group','Group','p.value','p.crit','sig'))
list(n=nv,n.keep=nv.keep,output=outp)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.