R/hc4wmc.R

hc4wmc <-
function(x,y,nboot=599,k=2,grp=NA,con=0,SEED=TRUE,STOP=TRUE,...){
#
#   Test the hypothesis that J independent groups have identical slopes.
#   Using least squares regression
#   Data are stored in list mode or in a matrix.  In the first case
#   x[[1]] contains the data for the first group, x[[2]] the data
#   for the second group, etc. Length(x)=the number of groups = J.
#   If stored in a matrix, the columns of the matrix correspond
#   to groups.
#
#   Similarly, y[[1]] contains the data for the first group,
#   y[[2]] the data for the second groups, etc.
#
#   The argument grp can be used to analyze a subset of the groups
#   Example: grp=c(1,3,5) would compare groups 1, 3 and 5.
#
#   Missing values are allowed.
#
if(STOP)stop('Suggest ols1way. This function assumes equal n. To use anyway, set STOP=FALSE')
con<-as.matrix(con)
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop("Data must be stored in list mode or in matrix mode.")
if(is.matrix(y))y<-listm(y)
if(!is.list(y))stop("Data must be stored in list mode or in matrix mode.")
if(!is.na(sum(grp))){  # Only analyze specified groups.
xx<-list()
yy<-list()
for(i in 1:length(grp))
xx[[i]]<-x[[grp[i]]]
yy[[i]]<-y[[grp[i]]]
x<-xx
y<-yy
}
J<-length(x)
n<-length(x[[1]])
tempn<-0
slopes<-NA
covar<-NA
stemp<-NA
yhat<-numeric(J)
res<-matrix(,ncol=J, nrow=n)
for(j in 1:J){
temp<-cbind(x[[j]], y[[j]])
temp<-elimna(temp) # Remove missing values.
#n<-length(y[[j]])
tempn[j]<-length(temp)
x[[j]]<-temp[,1]
y[[j]]<-temp[,2]
tempx<-as.matrix(x[[j]])
tempy<-as.matrix(y[[j]])
#Getting yhat and residuals for wild bootstrap
yhat[j]<-mean(tempy)
res[,j]<-tempy-yhat[j]
#original Slope and SE
stemp<-lsfit(tempx, tempy)
slopes[j]<-stemp$coef[k] #Slopes for original data
covar[j]<-lsfitNci4(tempx, tempy)$cov[k,k] #original HC4 for coefficient(slope)
}
#
Jm<-J-1
#
# Determine contrast matrix
#
if(sum(con^2)==0){
ncon<-(J^2-J)/2
con<-matrix(0,J,ncon)
id<-0
for (j in 1:Jm){
jp<-j+1
for (h in jp:J){
id<-id+1
con[j,id]<-1
con[h,id]<-0-1
}}}
ncon<-ncol(con)
if(nrow(con)!=J){
stop("Something is wrong with con; the number of rows does not match the number of groups.")
}
#calculating original statistic
dif.slopes<-t(con)%*%slopes
o.se<-t(con^2)%*%covar
o.stat<-dif.slopes/sqrt(o.se) #original test statistics
#
om<-max(abs(o.stat)) #Max. absolute test statistics
#
if(SEED)set.seed(2) # set seed of random number generator so that
#             results can be duplicated.
#
data<-matrix(ifelse(rbinom(n*nboot*J,1,0.5)==1,-1,1),ncol=nboot*J) #discrete wild bootstrap sample
test<-numeric(nboot)
u<-rep(1, n)
c<-1
for (i in 1:nboot*J-J+1){
d<-data[,i:i+J-1]
ystar<-u%*%t(yhat)+res*d
ystar<-listm(ystar)
i<-i+J
test[c]<-mcslope(x,ystar, con, k)
#
c<-c+1
}
sum<-sum(test>= om)
p.val<-sum/nboot
list(p.value=p.val)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.