3
$\begingroup$

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[1],x[2],sum(x[3:5]))
exptd2<-c(exptd[1],exptd[2],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.

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Maarten Punt Jul 19 '17 at 20:39
  • $\begingroup$ Did you really intend what you have here: mu<--fitted$estimate ?? That flips the sign. $\endgroup$ – Glen_b -Reinstate Monica Jul 19 '17 at 23:36
  • $\begingroup$ 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. $\endgroup$ – Glen_b -Reinstate Monica Jul 20 '17 at 0:25
  • $\begingroup$ As it currently stands this seems to be a purely "how do I write code for this in R" sort of question. $\endgroup$ – Glen_b -Reinstate Monica Jul 20 '17 at 8:57
  • $\begingroup$ 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. $\endgroup$ – John Jul 21 '17 at 7:04
1
$\begingroup$

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')}
  pacman::p_load(MASS)
  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[1] "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[1] "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

$\endgroup$
  • $\begingroup$ Thanks Hemant. Solves the problem but sadly I think my students will struggle with this :( $\endgroup$ – John Jul 21 '17 at 7:08
  • $\begingroup$ @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. $\endgroup$ – Hemant Rupani Jul 21 '17 at 8:09
  • $\begingroup$ @John or maybe you find a better option (Y) $\endgroup$ – Hemant Rupani Jul 21 '17 at 8:10
  • $\begingroup$ Good point - I'll just need to persuade the IFoA that this is a good idea. Thanks so much for your help. $\endgroup$ – John Jul 24 '17 at 7:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.