R/ranki.R

ranki <-
function(J,K,x,grp=c(1:p),alpha=.05,p=J*K){
#
#  Compute a confidence interval for all interaction terms
#  in J by K (two-way) anova using the modified Patel-Hoel method.
#
#  This method is not recommended if there are tied observations among the
#  pooled data.
#
#  All JK groups are independent.
#
#  The R variable x is assumed to contain the raw
#  data stored in list mode.
#  If grp is unspecified, it is assumed x[[1]] contains the data
#  for the first level of both factors: level 1,1.
#  x[[2]] is assumed to contain the data for level 1 of the
#  first factor and level 2 of the second factor: level 1,2
#  x[[j+1]] is the data for level 2,1, etc.
#  If the data are in wrong order, grp can be used to rearrange the
#  groups. For example, for a two by two design, grp<-c(2,4,3,1)
#  indicates that the second group corresponds to level 1,1;
#  group 4 corresponds to level 1,2; group 3 is level 2,1;
#  and group 1 is level 2,2.
#
#  It is assumed that the input variable x has length JK, the total number of
#  groups being tested. If not, a warning message is printed.
#
#  The estimated standard error is based on Sen's Jackknife as used by
#  Mee (1990).
#
if(!is.list(x))stop("Data are not stored in list mode")
if(p!=length(x)){
print("Warning: The number of groups in your data is not equal to JK")
}
jtk<-J*K
tl<-0
com<-x[[1]]
for(i in 1:jtk)tl<-tl+length(x[[i]])
for(i in 2:jtk)com<-c(com,x[[i]])
tiex<-sum(abs(c(1:tl)-sort(rank(com))))
if(tiex > 0)
print("Tied values detected. Interchanging columns might give different results. That is, comparing rows based on P(X<Y) is not necessarily the same as comparing rows based on P(X>Y)")
ck<-(K^2-K)/2
cj<-(J^2-J)/2
tc<-ck*cj
if(tc>28){
print("Warning: The number of contrasts exceeds 28.")
print("The critical value being used is based on 28 contrasts")
tc<-28
}
idmat<-matrix(NA,nrow=tc,ncol=8)
dimnames(idmat)<-list(NULL,c("row","row","col","col","ci.lower","ci.upper","estimate","test.stat"))
crit<-smmcrit(300,tc)
if(alpha != .05){
crit<-smmcrit01(300,tc)
if(alpha != .01){print("Warning: Only alpha = .05 and .01 are allowed,")
print("alpha = .01 is being assumed.")
}
}
phatsqse<-0
phat<-0
allit<-0
jcount<-0-K
it<-0
for(j in 1:J){
for(jj in 1:J){
if(j < jj){
for(k in 1:K){
for(kk in 1:K){
if(k < kk){
it<-it+1
idmat[it,1:4]<-c(j,jj,k,kk)
}}}}}
jcount<-jcount+K
for(k in 1:K){
for(kk in 1:K){
if(k < kk){
allit<-allit+1
xx<-x[[grp[k+jcount]]]
yy<-x[[grp[kk+jcount]]]
temp<-rankisub(xx,yy)
phat[allit]<-temp$phat
phatsqse[allit]<-temp$sqse
}}}}
#
# Compute the contrast matrix. Each row contains a 1, -1 and the rest 0
#  That is, all pairwise comparisons among K groups.
#
con<-matrix(0,cj,J)
id<-0
Jm<-J-1
for (j in 1:Jm){
jp<-j+1
for (k in jp:J){
id<-id+1
con[id,j]<-1
con[id,k]<-0-1
}}
IK<-diag(ck)
B<-kron(con,IK)
ntest<-ck*(J^2-J)/2
test<-0
civecl<-0
civecu<-0
for (itest in 1:ntest){
temp1<-sum(B[itest,]*phat)
idmat[itest,7]<-temp1
idmat[itest,8]<-temp1/sqrt(sum(B[itest,]^2*phatsqse))
idmat[itest,5]<-temp1-crit*sqrt(sum(B[itest,]^2*phatsqse))
idmat[itest,6]<-temp1+crit*sqrt(sum(B[itest,]^2*phatsqse))
}
nsig<-sum((abs(idmat[,8])>crit))
list(phat=phat,ci=idmat,crit=crit,nsig=nsig)
}
musto101/wilcox_R documentation built on May 23, 2019, 10:52 a.m.