GLM-Poisson-NegBin-crabs
========================  STAT 818R, 10/8/26

Now we look at "satell" count as response variable

   install.packages("CatDataAnalysis")  
   crabs = table_4.13

   fit2 = glm(satell ~ color + spine + width + weight, data=crabs, family=poisson)
    anova(fit2)
Analysis of Deviance Table,      Model: poisson, link: log

Response: satell,    Terms added sequentially (first to last)

       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                     172     632.79              
color   1   22.131       171     610.66 2.547e-06 ***
spine   1    0.089       170     610.57  0.765346    
width   1   50.409       169     560.16 1.248e-12 ***
weight  1    8.310       168     551.85  0.003943 ** 

  table(crabs$color)
 2  3  4  5 
12 95 44 22 

### Now for comparison, same fit with "color" made into 4-level factor

> anova(update(fit2, data=crabsF))

Response: satell   ,  Terms added sequentially (first to last)

       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                     172     632.79              
color   3   23.653       169     609.14 2.952e-05 ***
spine   1    0.100       168     609.04   0.75207    
width   1   49.798       167     559.24 1.704e-12 ***
weight  1    8.639       166     550.60   0.00329 ** 

#-----------------------------------------

   fit2S = step(update(fit2, satell ~ color), scope=list(lower = satell ~ 1,
               upper = satell ~ (color + spine + width + weight)^2))
   round(summary(fit2S)$coef,5)
             Estimate Std. Error  z value Pr(>|z|)
(Intercept)   3.02227    1.05050  2.87700  0.00401
color        -1.02443    0.31760 -3.22556  0.00126
weight       -0.00055    0.00040 -1.36607  0.17192
color:weight  0.00034    0.00012  2.76855  0.00563


   library(glmtoolbox)
   anova(glmgee(satell ~ color+weight+color:weight, data=crabs, family=poisson))

  Wald test 

Model 1 :  satell ~ 1 
Model 2 :  satell ~ color 
Model 3 :  satell ~ color + weight 
Model 4 :  satell ~ color + weight + color:weight 

         Chi    df   Pr(>Chi)    
1 vs 2  4.9271   1    0.02644 *  
2 vs 3 23.5877   1  1.193e-06 ***
3 vs 4  1.9861   1    0.15875       ### Look at this: the EE analysis 
                                  ## does not find color:weight  significant !!!

##=================================================================

   fit2Q = update(fit2S, family="quasipoisson")
   summary(fit2Q)$disp
[1] 3.261803              ### recover this another way !)

 with(crabs, sum( (satell - fit2S$fitted)^2/(fit2S$fitted) )/(173-4) ) ## = 3.2618
 with(crabs, sum( (satell - fit2$fitted)^2/(fit2$fitted) )/(173-5) ) ## = 3.24178


### So here we have legitimate OVERdispersion. Try NegBin next

   library(MASS)
   help(glm.nb)
   fit2NB = glm.nb(satell ~ (color + width)^2, data=crabs)
> summary(fit2NB)

Call:
glm.nb(formula = satell ~ (color + width)^2, data = crabs, init.theta = 0.9438649434, 
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.67444    5.72376   0.467    0.640
color       -1.85309    1.66585  -1.112    0.266
width       -0.04110    0.21657  -0.190    0.849
color:width  0.06401    0.06352   1.008    0.314

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

    Null deviance: 218.04  on 172  degrees of freedom
Residual deviance: 196.81  on 169  degrees of freedom
AIC: 757.88

Number of Fisher Scoring iterations: 1


              Theta:  0.944       ### we covered in class that this "theta" is the
          Std. Err.:  0.172       ### shape parameter for the random Gamma exp(beta0)

 2 x log-likelihood:  -747.882     ### compare with logLik(fit2S)

> c(logLik(fit2S), logLik(fit2NB))
[1] -450.1405 -373.9409

  with(crabs, sum( (satell - fit2NB$fitted)^2/(fit2NB$fitted + 
                               fit2NB$fitted^2/fit2NB$theta) )/(173-5) ) 
[1] 0.9553971

> anova(fit2NB)            ### This was a likelihood-based fit, not EE-based !!

Analysis of Deviance Table,   Model: Negative Binomial(0.9439), link: log

Response: satell,    Terms added sequentially (first to last)


            Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                          172     218.04              
color        1   5.4884       171     212.56 0.0191432 *  
width        1  14.6908       170     197.87 0.0001267 ***
color:width  1   1.0582       169     196.81 0.3036188    
---

Warning message:
In anova.negbin(fit2NB) : tests made without re-estimating 'theta'  ### NOTE THIS !!
                 ### This diminishes the relevance of the p-values for color and width
               ### deviances separately, but does NOT invalidate the test for color:width !!