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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._