R/Association and GoF.R

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

#' randomization test for goodness of fit
#'
#' \code{gofRand(exp, obs, perm = 1000)}
#'
#' @details produces a p-value from proportion of random data sets which
#' result in a chi-squared value at least as big as that calculated from the
#' actual data. Similar to the \code{simulate.p.value} option in \code{chisq.test}
#'
#' @param obs a vector of observed frequencies
#' @param exp a vector of expected frequencies, proportions, percentages or relative values
#' @param perm number of randomizations
#' @return list of output variables;
#' @return \code{randomisation p-value} = calculated p-value from randomization
#' @return \code{original chi-squared} = chi-squared calculated from actual data
#' @return \code{chi-squared distrib. p} = p-value from chi-squared distribution;
#' for comparison
#' @return \code{number of permutations}
#' @examples
#' Exp = c(9,3,3,1)
#' Obs = c(80,35,25,10)
#' gofRand(Obs,Exp)


gofRand=function(obs,exp,perm=1000){
  options(warn=-1)
  exp=exp/sum(exp)
  L=length(obs)
  S=sum(obs)
  M=1:L
  CHI=numeric(perm)
  Chiref=chisq.test(obs,p=exp)$stat
  Chip=chisq.test(obs,p=exp)$p.val
  for(a in 1:perm){
    Rand=sample(rep(M,round(100*S*exp)),S,replace=TRUE)
    Rand=factor(Rand)
    levels(Rand)=1:L
    Robs=table(Rand)[]
    CHI[a]=chisq.test(Robs,p=exp)$stat
  }
  pval=mean(CHI>Chiref)
  options(warn=0)
  options(digits=4)
  cat("\n","G-o-F Randomization test","\n",
      "\n","  original chi-squared =",Chiref,"\n",
      "chi-squared distrib. p =",Chip,"\n",
      "\n","randomisation p-value =",pval,"\n",
      "\n","number of permutations =",perm,"\n","\n")
  OUT=list('p.value.rand'=pval,'Chisq'=as.numeric(Chiref),'Chidist.p.val'=Chip)
  return(invisible(OUT))
}