[R] loops

[Heberto Ghezzo]

Hi R friends,
Reading the last issue of Stat Can J. there is an article in a single
degree of freedom test for non additivity of interactions in anova 
tables with 1 obs per cell. The authors claim it is more powerfull than 
Tukey but for the case of multiplicative interaction, which is the 
alternative studied by Tukey. B.A.Hartlaub,A.M.Dean,D.A.Wolfe. 
Rank-based test procedures for interaction in the two-way layout 
with one observation per cell. CanJ.Stat v27 p863-874 I tried to 
program it, and it works since the tables are 5 by 5 or 5 by 8 etc. 
Just for the sake of neatness how can I reprogram the core of the 
program which is:
 s <- NULL k <- 0 for ( j in
1:(nc- 1)) {
  for (jp in (j+1):nc ) {
    k <- k+1
    s[k] <- 0
    for (i in 1:(nr-1)) {
      for(ip in (i+1):nr) {
        a <- xr[i,j]+xr[ip,jp]-xr[i,jp]-xr[ip,j]
        s[k] <- s[k] + a*a
      }
    }
  }
}
cra <- max(s)

where nc=number of columns, nr number of rows and the xr's are 
the centered ranks.
the statistic being cra, there is a similar version where mean is 
replace by median in a previous step..
There must be a way in R to replace the 4 for loops i<i' and j<j' 

Peter Wolf (pwolf at wiwi.uni-bielefeld.de) suggested the function a1 
below.
This is may version of the cra.test, if somebody think of a way to 
improve it please let me know. The distribution must be simulated 
each time so speed of computation is crucial.
Thanks

cra <- function(x) {
  nc <- dim(x)[2]
  nr <- dim(x)[1]

  ca <- apply(x,2,mean)
  xc1 <- t(x)-ca
  xc <- t(xc1)
  xr1 <- apply(xc,1,rank)
  xr <- t(xr1)

  s <- NULL
  k <- 0
  for ( j in 1:(nc-1)) {
    for (jp in (j+1):nc ) {
      k <- k+1
      s[k] <- a1(xr[,j],xr[,jp]) 
     }
   }
  cra <- max(s)
}

a1 <- function(x,y){
  nr <- length(x)
  h <- x-y
  s <- 0
  for( i in 1:(nr-1)){
    ip <- (i+1):nr
    a <- h[i]-h[ip]
    s<- s+sum(a*a)
  }
  return(s)
}

cra.test <- function(x,n=1000){
# perform the cra test for interaction
#  ref: Hartlaub,Dean,Wolfe   Can.J.Stat v27,p863-874
#
# x = nc X nr matrix of the two-way table
# n = number of simulations
#
# output:
# cra : the value of the statistic
# plt : frequency less than observed in the n simulations
# ple : freq less and equal to  observed in simulations
# dist: quantiles from simulation(min, median,80%, 90%,95%,99% 
and max)
#
  nc <- dim(x)[2]
  nr <- dim(x)[1]
  tst <- cra(x)
  a <- 0
  b <- 0
  cr <- NULL
  k <- 0
  for (i in 1:n) {
    xx<-matrix(runif(nr*nc),ncol=nc)
    k <- k+1
    cr[k] <- cra(xx)
    if (tst < cr[k]) a <- a+1
    if (tst == cr[k]) b <- b+1
  }

return(list(cra=tst,plt=a/n,ple=(a+b)/n,dist=quantile(cr,probs=c(0,.5,.
8,.9,.95,.99,1)))) }


R. Heberto Ghezzo  Ph.D.
Meakins-Christie Labs
McGill University
Montreal - Canada
heberto at meakins.lan.mcgill.ca
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._