GLMMs

 

Hi Mr Bolker,

 

I would be very grateful to you if you could help with GLMMs.

 

I have the next model (R 11.1):

*y < cbind (fruits, flowers-fruits) *# none cell is = 0

*model <- lmer (y ~ plantheight + pH + ffertilization * fclipping + (1
|flocality), family=binomial, data= fruits), *where

y= proportion of flowers setting fruits

plantheight, pH: covariates

ffertilization, fclipping: fixed factors, two levels each one

flocality: random factor, six levels, four observations per locality
(balanced design)

Generalized linear mixed model fit by the Laplace approximation

Formula: y ~ plantheight + pH + ffertilization * fclipping + (1 |
flocality)

   Data: frutos

   AIC   BIC logLik deviance

 101.1 109.4 -43.57    87.13

Random effects:

 Groups     Name        Variance Std.Dev.

 flocality(Intercept) 0.53584  0.73201

Number of obs: 24, groups: flocality, 6

 

Fixed effects:

                                               Estimate   Std. Error
    z value  Pr(>|z|)   

(Intercept)                              -0.465952   1.127642  -0.413 
  0.6795   

plantheight                           0.016867   0.007252   2.326  
  0.0200 * 

pH                                         -0.118768   0.248236 
-0.478     0.6323   

ffertilization2                        0.150188   0.100795   1.490  
  0.1362   

fclipping2                             0.712214   0.096554   7.376
    1.63e-13 ***

ffertilization2:fclipping2    -0.669289   0.140147  -4.776     1.79e-06 ***

The point it is that data exhibit large overdispersion, so I refit the
model using quasi-binomial.

Generalized linear mixed model fit by the Laplace approximation

Formula: y ~ plantheight + pH + ffertilization * fclipping + (1 |
flocality)

   Data: frutos

   AIC   BIC logLik deviance

 103.1 112.6 -43.57    87.13

Random effects:

 Groups     Name        Variance  Std.Dev.

 flocality (Intercept) 0.0016303 0.040378

 Residual               0.0030426 0.055160

Number of obs: 24, groups: flocality, 6

 

Fixed effects:

                                                          Estimate
               Std. Error                  t value

(Intercept)                                     -0.465952  
        0.062200                     -7.49

plantheight                                 0.016867  
            0.000400                     42.17

pH                                                -0.118768  
            0.013693                       -8.67

ffertilization2                               0.150188  
        0.005560                    27.01

fclipping2                                     0.712214  
         0.005326                     133.73

ffertilization2:fclipping2      -0.669289            0.007731
                    -86.58

 

All factors seem to have significant effect in the model. I am aware of
the ?statistics limitations? to get pvalues in this case so I decided to
use MCMC:

*markov1=mcmcsamp (model,5000)*

Error en .local(object, n, verbose, ...) : Update not yet written

*HPDinterval(markov1)*

Error en HPDinterval(markov1) :

Error in evaluating the argument 'object' in selecting a method for
function 'HPDinterval': Error: objeto 'markov1' no encontrado

 

?Is it possible that update is not written?

 

Thus, I thought that may be could be better to try with a newer R
version (R 13.0), but I was very surprised when the next error appeared:

 

*model <- lmer (y ~ plantheight + pH + ffertilization * fclipping + (1

|flocality), family=quasibinomial, data= fruits)*

Error en glmer(formula = y ~ alturaplanta + pH + ffertilizaci?n *  :

  "quasi" families cannot be used in glmer

 

First, I wasn?t using ?glmer?, but ?lmer?; and second, I was able to
obtain ?t? values for the same model with R 11.1. Why not now???????

 

 

 

An alternative is to compare models in order to obtain significance for
the fixed factors.

 

*model <- lmer (y ~ plantheight + pH + ffertilization * fclipping + (1
|flocality), family=quasibinomial, data= fruits)*

*model2<- lmer (y ~ plantheight + pH + ffertilization + (1 |flocality),
family=quasibinomial, data= fruits)*

*anova (model, model2) *# to get significance of the factor fclipping

 

Data: fruits

Models:

model2: y ~ alturaplanta + pH + ffertilizaci?n + (1 | flocalidad)

model: y ~ alturaplanta + pH + ffertilizaci?n * fherbivor?a + (1 |
flocalidad)

       Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)   

model2  6 153.79 160.85 -70.893                            

model   8 103.13 112.56 -43.566 54.654      2  1.355e-12 ***

---

Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

 

The problem here is that F (instead of Chi) values are supposed to be
given when you use quasi-.

 

*> anova (model, model2, test = "F")*

Data: fruits

Models:

model2: y ~ alturaplanta + pH + ffertilizaci?n + (1 | flocalidad)

model: y ~ alturaplanta + pH + ffertilizaci?n * fherbivor?a + (1 |
flocalidad)

       Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)   

model2  6 153.79 160.85 -70.893                            

model   8 103.13 112.56 -43.566 54.654      2  1.355e-12 ***

---

Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

 

It is exactly the same result as above. So I afraid that the model
comparison is not taking into account the overdispersion.

 

 

Thank you very much,

 

Regards

J

 

*________________________________________*

Joaquina Pato

 

?rea de Ecolog?a

Departamento de Biolog?a de Organismos y Sistemas
Universidad de Oviedo

&

Unidad Mixta de Investigaci?n en Biodiversidad

www.unioviedo.es/icab <http://www.unioviedo.es/icab>

 

 

Campus del Cristo

C/ Catedr?tico Rodrigo Ur?a s/n
E-33071 Oviedo (Asturias)
Espa?a


Tfno: +34 985 104831
Fax: +34 985 104866
e-mail: patojoaquina at uniovi.es <mailto:patojoaquina at uniovi.es>