1
$\begingroup$

This is my first post on this forum so let me know if I am missing anything. I have count data which is over-dispersed, which I am treating as the outcome variable. I am looking to model the relationship between this and a set of predictive variables. I have done a poisson, quasi-poisson and now a negative binomial. My data file is named LR_tab7. This is my code for the three glm (bc sorry for the weird data names):

> LR_tab7 <- read.csv("~/ongoing/24) Toine Run/LR_tab7.csv")            
> library(MASS)

###quaispoisson

> pois<-glm(lang~Q8_early_lang+Group+Age+sex.2,family=quasipoisson(identity), data=LR_tab7) >summary(pois)                                                                            

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8729  -1.0403  -0.8129   0.6000   2.7736  

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   0.47234    0.14498   3.258 0.001181 ** 
Q8_early_lang[T.two_to_2.5]   0.27025    0.07606   3.553 0.000409 ***
Q8_early_lang[T.greater_2.5]  0.56247    0.14042   4.006 6.90e-05 ***
Q8_early_lang[T.greater_3]    1.27225    0.31448   4.046 5.86e-05 ***
Group[T.PR]                   0.15055    0.07768   1.938 0.053049 .  
Group[T.AF]                   0.25314    0.07901   3.204 0.001423 ** 
Age                          -0.01997    0.01132  -1.765 0.078020 .  
sex.2[T.M]                    0.02245    0.05970   0.376 0.706965    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasipoisson family taken to be 1.081649)

    Null deviance: 770.83  on 649  degrees of freedom
Residual deviance: 691.00  on 642  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 8
###poisson

> poiswithoutquasi<-glm(lang~Q8_early_lang+Group+Age+sex.2,family=poisson(identity),
data=LR_tab7)                                        
>summary(poiswithoutquasi)
   Coefficients:
                                 Estimate Std. Error z value Pr(>|z|)    
    (Intercept)                   0.47234    0.13940   3.388 0.000703 ***
    Q8_early_lang[T.two_to_2.5]   0.27025    0.07313   3.695 0.000220 ***
    Q8_early_lang[T.greater_2.5]  0.56247    0.13501   4.166 3.10e-05 ***
    Q8_early_lang[T.greater_3]    1.27225    0.30238   4.207 2.58e-05 ***
    Group[T.PR]                   0.15055    0.07469   2.016 0.043834 *  
    Group[T.AF]                   0.25314    0.07597   3.332 0.000862 ***
    Age                          -0.01997    0.01088  -1.836 0.066393 .  
    sex.2[T.M]                    0.02245    0.05740   0.391 0.695682    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    (Dispersion parameter for poisson family taken to be 1)

        Null deviance: 770.83  on 649  degrees of freedom
    Residual deviance: 691.00  on 642  degrees of freedom
    AIC: 1332.4

    Number of Fisher Scoring iterations: 8

###negative binomial

> NB<-glm.nb(lang~Q8_early_lang+Group+Age+sex.2, data=LR_tab7)          
> > summary(NB)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8692  -1.0287  -0.8478   0.6020   2.7092  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                  -0.783814   0.293258  -2.673 0.007522 ** 
Q8_early_lang[T.two_to_2.5]   0.452437   0.117263   3.858 0.000114 ***
Q8_early_lang[T.greater_2.5]  0.790672   0.146478   5.398 6.74e-08 ***
Q8_early_lang[T.greater_3]    1.272113   0.187517   6.784 1.17e-11 ***
Group[T.PR]                   0.220222   0.192779   1.142 0.253306    
Group[T.AF]                   0.436158   0.186818   2.335 0.019561 *  
Age                          -0.029722   0.021942  -1.355 0.175557    
sex.2[T.M]                    0.005163   0.104520   0.049 0.960603    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(239.2127) family taken to be 1)

    Null deviance: 768.89  on 649  degrees of freedom
Residual deviance: 690.69  on 642  degrees of freedom
AIC: 1335.8

Number of Fisher Scoring iterations: 1


              Theta:  239 
          Std. Err.:  5570 
Warning while fitting theta: alternation limit reached 

 2 x log-likelihood:  -1317.8

I am trying to use the likelihood ratio to show the negative binomial is a better fit model. I couldn't use the quasipoisson in the likelihood ratio test because it kept coming up as NA. This is my code for the regular poisson and negative binomial lrtest:

> X2<-2*(logLik(NB)-logLik(poiswithoutquasi))                           
> X2                                                                    
> 'log Lik.' -1.383039(df=9)                                            
> pchisq(X2, df = 1, lower.tail=FALSE)                                  
> 'log Lik.' 1 (df=9)

Could someone help me interpret what this means and what I am doing wrong. I have very little statistical or R background and I am stuck on how to validate your model.

sample is 650, average count is 0.6 (there are options 0-3 on a scale)

$\endgroup$
1
  • $\begingroup$ What's your sample size, what's the average count? $\endgroup$ – Aksakal Oct 14 '16 at 18:27
1
$\begingroup$

First of all your negative binomial model was not fitted properly as the warning says, the overdispersion is quite large and the standard error is huge. You'd better fix this first, since it does not convince me that the negative binomial actually provides a better fit.

Secondly you cannot use a regular likelihood ratio test to show that the negative binomial fits better than the poisson (or more formally: prove that the overdispersion parameter is significantly different from zero). This is because the overdispersion parameter is restricted to be positive, and your null hypothesis

$$\theta=0$$

is thus on the edge of the parameter space. As a result the null distribution of the test statistic is an even mixture of a point mass at zero and the $\chi_{(1)}$ distribution (see e.g. Loeys et al.,2012). If you're not an expert at the subject I can recommend the odTest() function in the pscl package.

$\endgroup$

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.