[R-sig-ME] GLMMs
Ben Bolker
bbolker at gmail.com
Mon Jun 13 22:26:45 CEST 2011
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
Dear Joaquina,
I'm trying to cut down on off-list answers. Please send this to
r-sig-mixed-models at r-project.org ...
sincerely
Ben Bolker
On 06/13/2011 09:36 AM, Joaquina wrote:
>
>
> 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>
>
>
>
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/
iEYEARECAAYFAk32coUACgkQc5UpGjwzenMDmQCeKsB7kGTA3KfkHS6RveKnTXSL
enoAnjTTBzONLbBNUJSdMdRXBXyLw0V8
=Y/CG
-----END PGP SIGNATURE-----
More information about the R-sig-mixed-models
mailing list