# Using chisq.test in R to measure goodness of fit of a fitted distribution

I have a frequency table of claims made by policyholders:

Claims       0,     1,    2,    3,   4
Frequency  2962,  382,   47,   25,   4


And I am going to carry out a chi-squared goodness of fit test to see if it conforms to a Poisson distribution (there are probably far better methods - but I'm teaching basic stats - so go with the flow please).

I converted the frequency table into a vector as follows:

n<-c(0,1,2,3,4)
x<-c(2962,382,47,25,4)
data <- rep(n,x)


I then fitted a Poisson(mu) as follows:

library(MASS)
fitted <-fitdistr(data,"Poisson")
mu<-fitted$estimate  I obtained the expected probabilities under this model using: exptd<-c(dpois(0,mu),dpois(1,mu),dpois(2,mu),dpois(3,mu),1-ppois(3,mu))  Using round(sum(x)*exptd,3) I noticed that the expected frequency of the last two groups were both <5 and so I combined the last three groups together: x2<-c(x,x,sum(x[3:5])) exptd2<-c(exptd,exptd,sum(exptd[3:5]))  I now carried out my chi-squared goodness of fit test: chisq.test(x2,p=exptd2)  I have two issues: 1. The chi-square test has the wrong degrees of freedom as I estimated a parameter on the data. 2. This seems dreadfully long-winded - there must be a quicker way of doing this test. Unfortunately despite searching I can't find examples of this online (which makes me think that this method isn't used despite what our textbooks say) so I need your help. Many thanks in advance for the help. • The lsr package from Daniel Navarro that comes with the book Learning Statistics with R has a nice wrapper function for the chi-square test. Input are a vector of observed frequencies and probability vector. Output is a more verbose version of the chisq.test(). That should reduce your problem by a few steps Jul 19, 2017 at 20:39 • Did you really intend what you have here: mu<--fitted$estimate ?? That flips the sign. Jul 19, 2017 at 23:36
• vcd::goodfit will do it with one function call, or see fitdistr::gofstat. The more direct alternative (and perhaps better from a teaching point of view) would be to call pchisq to get the p-value with the correct degrees of freedom. Jul 20, 2017 at 0:25
• As it currently stands this seems to be a purely "how do I write code for this in R" sort of question. Jul 20, 2017 at 8:57
• Only if the chisq.test doesn't do what I want it to do. Which given the replies seems to be the case. :( So I shall start exploring the suggested alternatives.
– John
Jul 21, 2017 at 7:04

User Defined Function can be used for such requirements

I tried to make a naive function for mentioned purpose. you can expand this.

chi2test=function(data,model,dens,df=length(unique(data))-length(params)-1){
if(!require(pacman)) {install.packages('pacman')}
params=fitdistr(data,model)\$estimate
exptd=dens(unique(data),params)
exptd[length(exptd)]=1-sum(exptd[-length(exptd)])
normalise=sum((table(data)-length(data)*exptd)^2/(length(data)*exptd))
cat('pvalue','=',pchisq(normalise,df,lower.tail = F),';','df','=',df)
if(any(length(data)*exptd<5)) {print("Warning: there are some expected frequencies that are less than 5")}
}


data should be a vector class.

model is the model defined under 'densfunc' under 'fitdistr' under package MASS... should not be write as object, but under quotation marks.

dens is the density function from default stats package.

df is set default, but you can also add by df=n. where n is positive Integer.

Data

data
0    1    2    3    4
2962  382   47   25    4

data2
0    1    2
2962  382   76


Now, trials

chi2test(data,'Poisson',dpois)


Ouput:- pvalue = 6.845853e-91 ; df = 3 "Warning: there are some expected frequencies that are less than 5"

 chi2test(data2,'Poisson',dpois)


Output:- pvalue = 5.72538e-13 ; df = 1

chi2test(data2,'Poisson',dpois,df=2)


Output:- pvalue = 5.267489e-12 ; df = 2

chi2test(data,'geometric',dgeom)


Output:- pvalue = 2.281883e-09 ; df = 3 "Warning: there are some expected frequencies that are less than 5"

chi2test(data2,'geometric',dgeom)


Output:- pvalue = 0.05292066 ; df = 1

chi2test(data2,'geometric',dgeom,df=2)


Output:- pvalue = 0.1536295 ; df = 2

• Thanks Hemant. Solves the problem but sadly I think my students will struggle with this :(
– John
Jul 21, 2017 at 7:08
• @John Hi... But you can put time one time for Function or even a package, you'll anyway teach your students to install some packages. Jul 21, 2017 at 8:09
• @John or maybe you find a better option (Y) Jul 21, 2017 at 8:10
• Good point - I'll just need to persuade the IFoA that this is a good idea. Thanks so much for your help.
– John
Jul 24, 2017 at 7:35