CODEWITHSUNDEEP
rplotquantilepoisson
no_one
no_one

Reputation: 108

How to create a Q-Q plot with Poisson as theoretical distribution

I need to create a Q-Q plot in order to check if my observed data fits a Poisson distribution.

Here is my data.frame:

df = read.table(text = 'Var1 Freq
 1975   10
 1976   12
 1977    9
 1978   14
 1979   14
 1980   11
 1981    8
 1982    7
 1983   10
 1984    8
 1985   12
 1986    9
 1987   10
 1988    9
 1989   10
 1990    9
 1991   11
 1992   12
 1993    9
 1994   10', header = TRUE)

The df$Freq column is the one that interests me as the observations represent the count of events per year.

I know I have to use the qqplot function and also the qpois one to create the theoretical quantiles, but how?

Upvotes: 1

Views: 13172

Answers (3)

paqmo
paqmo

Reputation: 3739

Also, the fitdistrplus package can do this with much less code. Compares the empirical and theoretical density and CDFs.

library('fitdistrplus')
plot(fitdist(df$Freq,"pois"))

You can get your lambda, etc. and check against other distributions as well. Not as flexible as the ggplot approach, but good for a quick check.

Upvotes: 5

no_one
no_one

Reputation: 108

here my possible answer:

#calculate Frequencies
tbl = as.data.frame(table(df$Freq))

#create theoretical poisson distr
dist = dpois(1:7, lambda = mean(tbl$Freq))
dist = dist * 20              #make values in the same scale as tbl$Freq (20 = sum(tbl$Freq))
dist = as.data.frame(dist)
dist$Var1 = tbl$Var1

#qqplot
qqplot(dist$dist, tbl$Freq, xlab = 'Theoretical Quantiles', ylab = 'Empirical Quantiles',
       main = 'Q-Q plot Poisson', xlim = c(0,5), ylim = c(0,5))
abline(0,1) #create 45° line

If you spot any mistake please let me know. Thanks

Upvotes: 2

Benjamin
Benjamin

Reputation: 17279

ggplot2 has a nice interface for doing this. Here's a QQ plot with an agreement step line in red. The QQ plot is made using stat_qq and changing the distribution argument. You'll need to provide lambda in the dparams argument.

ggplot(data = df,
       mapping = aes(sample = Freq)) + 
  stat_qq(distribution = stats::qpois,
          dparams = list(lambda = mean(df$Freq))) + 
  geom_step(data = data.frame(x = 6:16,
                              Freq = 6:16),
            mapping = aes(x = x,
                          y = Freq),
            colour = "red",
            alpha = 0.5)

Upvotes: 2

Related Questions