1
$\begingroup$

I have the following DF:

> str(df)
'data.frame':   80 obs. of  3 variables:
 $ Transect: Factor w/ 3 levels "A","B","E": 1 1 1 1 1 1 1 1 2 2 ...
 $ BI      : Factor w/ 2 levels "B","I": 2 2 2 2 2 2 2 2 2 2 ...
 $ Nfor    : int  51 5 15 5 24 17 20 9 3 2 ...

I'm trying to run some linear models with Nfor as DV, BI as a fixed effect and Transect as a random effect. When I use glmer with Poisson family, I get differing intercepts (which is what I need):

> model_2 <- glmer(Nfor ~ BI + (1|Transect),  family=poisson, data=DF)
> coef(model_2)
$Transect
  (Intercept)        BII
A    4.189409 -0.7751179
B    3.690552 -0.7751179
E    4.024364 -0.7751179

However that model has overdispersion, so I tried Quasi-Poisson (using MASS::glmmPQL()):

model_1 <- glmmPQL(Nfor ~ BI, random = ~1|Transecto, 
                       family=quasipoisson, data=DF)
coef(model_1)

> coef(model_1)
(Intercept)         BII 
  3.8809629  -0.4663443

As you can see, the intercept does not change between the levels of Transect, which I guess means that the random effects are not being specified correctly. I tried using glmmADMB::glmmadmb() (with "nbinom" family) and the same thing happens. I then ran models using using lmer and lme, just to see if that would make any difference, but it doesn't. Does anybody have a clue about this?

EDIT:

Summary of both models:

> model_2 <- glmer(Nfor ~ BI + (1|Transect), family = poisson,
+ data = df)
> summary(model_2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: Nfor ~ BI + (1 | Transect)
   Data: df

     AIC      BIC   logLik deviance df.resid 
  1880.7   1886.3   -937.3   1874.7       45 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.8111 -3.9006 -2.2600  0.9005 24.5002 

Random effects:
 Groups   Name        Variance Std.Dev.
 Transect (Intercept) 0.04471  0.2115  
Number of obs: 48, groups:  Transect, 3

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  3.96734    0.12530   31.66   <2e-16 ***
BII         -0.77512    0.04945  -15.68   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
    (Intr)
BII -0.124


> model_1 <- glmmPQL(Nfor ~ BI, random = ~1|Transect, family = quasipoisson,
+ data = df)
iteration 1
> summary(model_1)
Linear mixed-effects model fit by maximum likelihood
 Data: df 
  AIC BIC logLik
   NA  NA     NA

Random effects:
 Formula: ~1 | Transect
         (Intercept) Residual
StdDev: 3.723954e-05 6.752193

Variance function:
 Structure: fixed weights
 Formula: ~invwt 
Fixed effects: Nfor ~ BI 
                Value Std.Error DF   t-value p-value
(Intercept)  3.988984 0.1915937 44 20.820019   0.000
BII         -0.775121 0.3411702 44 -2.271947   0.028
 Correlation: 
    (Intr)
BII -0.562

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.0278538 -0.6084323 -0.3377733  0.2642417  4.2524930 

Number of Observations: 48
Number of Groups: 3 
$\endgroup$
  • 1
    $\begingroup$ You might be better asking on an R-specific site. You could also revise the documentation for the two different extractor functions called coef which you are using before you do that. $\endgroup$ – mdewey Dec 14 '16 at 18:28
  • $\begingroup$ Thank you @mdewey. I first thought about posting on StackOverflow, but then I thought that the problem might be related to the statistics itself, and not the code or functions, that's why I tried here first. I also looked at the documentation for coef but I did not find anything useful there (although maybe that's just because I'm the one that didn't understand it). $\endgroup$ – Lfppfs Dec 14 '16 at 18:35
  • $\begingroup$ I just examined one of my recent quasipoisson models built with glmmPQL and I definitely have different intercepts for each random effect. In eyeballing what you've presented above, it looks like you only have one level to your random effect in the glmmPQL model. However, I can't confirm this is the case. Could you recreate the problem in a dummy dataset? And subsequently share that dataset? That might help us to troubleshoot your problem. $\endgroup$ – Reid Dec 14 '16 at 18:35
  • $\begingroup$ @Reid there are three levels in Transect. By dummy dataset you mean transforming Transect into a dummy variable? If that's so, I did that but nothing changes. Here, I just put a .txt on GitHub (it's "MyDF.txt"): dataset $\endgroup$ – Lfppfs Dec 14 '16 at 19:27
  • $\begingroup$ Another option is to allow for overdispersion by including an observation specific random effect in glmer. $\endgroup$ – HStamper Dec 14 '16 at 20:04
1
$\begingroup$

Summary

glmmPQL uses the penalized quasilikelihood estimation method, which can produce biased estimates with discrete dependent variables and a low sample size. Low sample size is usually defined as less than 5 per cell, but you're pushing that with 8 per cell.

My guess is that you're looking at a combination a low sample size and a low effect size.

Exercise 1

1) Try running your model in the glmmPQL module, but change the 'quasipoisson' to a regular 'poisson.' You'll notice that nothing really changes.

model_1b <- glmmPQL(Nfor ~ BI, random = ~1|Transect, 
               family=poisson, data=DF)

coef(model_1b)

| (Intercept)| BII | | -----------|----------| | 3.988984 |-0.7751208| | 3.988984 |-0.7751208| | 3.988984 |-0.7751208|

That's because the glmmPQL module is still using a quasilikelihood estimation method, which runs into trouble with low sample size.

Exercise 2:

1) Download a modified version of your file here where Transect 1 has consistently larger numbers.

2) Plug these new data into your model and examine the coefficients.

model_1c <- glmmPQL(Nfor ~ BI, random = ~1|Transect, 
               family=quasipoisson, data=altered_df)

coef(model_1c)

| (Intercept) | BII | | ----------- |--------- | | 5.551994 |-0.7827377| | 3.801897 |-0.7827377| | 4.081167 |-0.7827377|

In summary I wouldn't worry too much about this. It looks like the Transect variable is not really contributing a lot to the model.

You can try leaving the Transect variable out, then subsequently test whether it changes the model significantly. If it has an effect, then keep it included. If it doesn't have an effect, then you'll have enough evidence to exclude this variable.

$\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.