R/assrand.R

In helophilus/ColsTools: A variety of convenience tools and short-cuts

#' randomization method for test of association/independence in lieu of chi-squared
#' test/Fisher's exact test.
#'
#' \code{assrand(matr,perm=1000)}
#'
#' @param matr is the name of a flat matrix representing the contingency table.
#' @param perm is the number of randomizatins to use.
#'
#' @return \code{Chisq} is the chi-squared value from the actual data
#' @return \code{p.val.rand} is the p-value from the randomization
#'
#' @details the actual data matrix is used to calculate expected frequencies if the
#' variables were independent. A value of chi-squared is calculated based upon these
#' observed and expected contingency tables. Randomization occurs to maintain row
#' and column totals, but distribute the frequencies at random within that constraint.
#' A new chi-squared value is calculated for each of these randomized matrices and the
#' proportion at least as big as the actual value is returned as the p-value.
#'
#' @examples
#' M = matrix(1:4*10,2)
#' assrand(M)

assrand=function(matr,perm=1000){
  options(warn=-1)
  Xref=chisq.test(matr,correct=FALSE)$stat
  Rs=dim(matr)[1]
  Cs=dim(matr)[2]
  AC=NULL;BC=NULL
  for(a in 1:Rs){
    for(b in 1:Cs){
      N=matr[a,b]
      A=rep(a,N)
      B=rep(b,N)
      AC=c(AC,A)
      BC=c(BC,B)
    }}
  M=data.frame(AC,BC)
  X2=numeric(perm)
  for(c in 1:perm){
    M$BC=sample(M$BC)
    X2[c]=chisq.test(xtabs(~.,M),correct=FALSE)$stat
  }
  Pr=mean(X2>=Xref)
  cat('\n',"Randomised Association Test",'\n','\n',
      'permutations =',perm,'\n',
      ' Chi-squared =',Xref,'\n',
      ' probability =',Pr,'\n')
  options(warn=0)
  OUT=list('Chisq'=as.numeric(Xref),'p.val.rand'=Pr)
  return(invisible(OUT))
}