Trying to figure out an easy way to do Pearson’s Chi-Squared Test in R when verifying goodness of fit. I thought the function chisq.test() would help!

Problems?

- There’s no way to tell chisq.test the correct degrees of freedom (“df”) – it cannot figure this out itself. My stats text tells me to reduce df by 1 for each value derived from the observations. If I have 10 buckets of counts
**AND**I have calculated the mean and sd from the data, df = 7**NOT**9 as assumed by chisq.test(). - In the test distribution, again according to my text, no expected count should be less than 5. All the chisq.test does is tell you the results are unreliable – it doesn’t fix them leaving the pre-processing up to you.

So here’s how to use R’s chisq.test()

First, I use a little function to set up some likely buckets for the observations “obs”:

Ian.Chisq.Buckets <- function(obs){

obs.hist <- hist(obs,breaks=seq(min(obs), max(obs), length=(1 + length(obs) / 5))) return(obs.hist)

}

This returns an object of class “histogram” which can be passed to the next step once the distribution to be tested for goodness of fit is set up. For example: test.cdf <- pnorm(obs.hist$breaks, mean(obs), sd(obs)) generates a cdf that corresponds to the buckets in the histogram of the observations.

Passing the histogram of the observations and the cdf to the next function accumulates the buckets to ensure a minimum count of 5 in every bucket for the expected distribution, and accumulates the buckets in parallel for the observed distribution. Then the chisq.test() is run, setting rescale.p=T, because counts not probabilities are being passed.

Ian.Chisq.Test <- function(obs.hist, cdf){

if (!(class(obs.hist) == “histogram” | class(obs.hist) == “numeric”)) stop(“problem with obs.hist”)

if (length(obs.counts) != length(cdf) – 1) stop(“observations and cum dist function mismatch”)

if (any(obs.counts < 0) | any(cdf < 0) | any(cdf > 1)) stop(“Counts or cdf messed up”)

if (class(obs.hist) == “histogram”){

obs.counts <- obs.hist$counts

}else{

obs.counts <- obs.hist

}

cdf <- sort(cdf)

cdf[length(cdf)] <- 1

exp.counts <- length(obs) * (cdf[-1] - cdf[-length(cdf)])

mid.point <- round(length(exp.counts)/2, 0)

for (i in 1:mid.point){

if (exp.counts[i] < 5){

exp.counts[i+1] <- exp.counts[i+1] + exp.counts[i]

exp.counts[i] <- 0

obs.counts[i+1] <- obs.counts[i+1] + obs.counts[i]

obs.counts[i] <- 0

}

}

for (i in length(exp.counts):mid.point){

if (exp.counts[i] < 5){

exp.counts[i-1] <- exp.counts[i-1] + exp.counts[i]

exp.counts[i] <- 0

obs.counts[i-1] <- obs.counts[i-1] + obs.counts[i]

obs.counts[i] <- 0

}

}

obs.counts <- obs.counts[exp.counts != 0]

exp.counts <- exp.counts[exp.counts != 0]

if (length(exp.counts) < 2){

stop(“insufficient expected observations”)

}else{

mess.str <- "If 'df' is incorrect, use Ian.Chisq.P(chisq.results, df) to correctly calculate P-Value"

return(c(chisq.test(obs.counts, p=exp.counts, rescale.p=T), “Message”=mess.str))

}

}

The chisq.test() results are returned with a reminder to think about df.

If df is actually different, then pass the whole set of test results plus the correct df to the last function:

Ian.Chisq.P <- function(chisq.results, df){

p.stat <- pchisq(chisq.results$statistic, df, lower.tail=F)

names(p.stat) <- "P-Statistic"

return(p.stat)

}

The correct P-Statistic is returned!

## Share

If you found this post informative, please share it!