Syntaxe pour lmer?

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Syntaxe pour lmer?

Messagepar Karine Poitrineau » 04 Juil 2006, 23:00

Bonjour à tous,
J'essaie de comprendre comment utiliser la fonction lmer pour réaliser une régression logistique.
Mes données ont cette allure:

Code : Tout sélectionner

        souche  repetition  encapsul= Réponse= variable binaire
1059    A45             bt           0 
1060    A45             bt           0 
1061    A45             bt           0 
1062    A45             bt           0 
1063    A45             bt           0 
1064    SF1             bu           1 
1065    SF1             bu           1 
 ...
 ...
1078   SF1              bv           0
1079   SF1              bv           0
1080   SF1              bv           0
1081   SF1              bv           0
1082   SF1              bv           0

J'ai une réponse binaire, et je cherche à déterminer si elle diffère entre mes souches (n=30) sachant que pour chaque souche on a des répétition ( généralement 4 répétitions, pour lesquelle on obtient un nombre variable de réponses 0 ou 1).

J'ai donc tenté de modéliser ceci, mais je ne suis pas sûre de moi...
Est-ce que la formule:
lmer(encapsul ~souche+ (1|souche:repetition), family="binomial")->mix3 est correcte par rapport à ce que je souhaite tester?

Si j'utilise ensuite une anova, j'obtiens:

Code : Tout sélectionner

> anova(mix3)
Analysis of Variance Table
       Df  Sum Sq  Mean Sq   Denom F value Pr(>F)
souche 29    0.14 0.004704 1660.00  0.0072      1


Si je ne me trompe pas, la conclusion de tout ceci est que je n'ai pas de différence significative entre mes souches?
Malheureusement, si je m'en reporte à cette discussion: http://finzi.psych.upenn.edu/R/Rhelp02a ... 61884.html je ne sais pas quel crédit je dois accorder au résultat de l'anova.

Le sommaire de la régression donnant ceci:

Code : Tout sélectionner

> summary(mix3)
Generalized linear mixed model fit using PQL
Formula: encapsul ~ souche + (1 | souche:repetition)
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 1065.961 1239.801 -500.9807 1001.961
Random effects:
           Groups              Name          Variance          Std.Dev.
souche:repetition       (Intercept)            5.1336            2.2657
# of obs: 1690, groups: souche:repetition, 111

Estimated scale (compare to 1)  0.7750342

Fixed effects:
             Estimate Std. Error  z value Pr(>|z|) 
(Intercept) -3.464559   1.380345 -2.50992  0.01208 *
soucheA21    2.566020   1.883634  1.36227  0.17311 
soucheA22    1.680399   1.958827  0.85786  0.39097 
soucheA28    1.473629   1.888869  0.78016  0.43529 
soucheA3     0.983925   1.922013  0.51192  0.60870 
soucheA32    1.780925   1.816156  0.98060  0.32679 
soucheA34    1.169993   2.067607  0.56587  0.57148 
soucheA41    1.900044   1.829047  1.03882  0.29889 
soucheA45    0.254297   2.024243  0.12563  0.90003 
soucheA5     2.605590   1.870204  1.39321  0.16356 
soucheAc11   1.968235   1.910949  1.02998  0.30302 
soucheAc13  -0.891319   2.061694 -0.43232  0.66551 
soucheAc14   3.994075   1.929978  2.06949  0.03850 *
soucheAc16   2.593285   1.917688  1.35230  0.17628 
soucheAc21   2.074794   1.979327  1.04823  0.29453 
soucheAc3    0.673147   1.951983  0.34485  0.73021 
soucheAc9    3.825062   1.952796  1.95876  0.05014 .
soucheAcX    1.004453   1.879326  0.53448  0.59301 
soucheAcY    1.935509   1.824665  1.06075  0.28880 
soucheSF1    0.873992   1.867653  0.46796  0.63981 
soucheSF11   4.128431   1.940136  2.12791  0.03334 *
soucheSF12   1.930819   1.882131  1.02587  0.30495 
soucheSF2    3.537742   2.023653  1.74820  0.08043 .
soucheSF3    0.923809   2.023648  0.45651  0.64803 
soucheSF4    1.606224   2.107019  0.76232  0.44587 
soucheSF5   -1.453207   2.172547 -0.66890  0.50356 
soucheSF6    2.309963   1.976730  1.16858  0.24257 
soucheSF7    0.080365   1.902902  0.04223  0.96631 
soucheSF8    3.558574   1.916392  1.85691  0.06332 .
soucheSF9    1.357004   1.983106  0.68428  0.49380 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Bref, je ne sais pas si j'ai bien procédé, je ne sais pas quoi conclure, ni si je peux conlure quelque chose, ou s'il y a d'autres façons de procéder...

Merci d'avance pour tout conseil ou remarque...

K.i espère ne pas se couvrir de ridicule...

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 05 Juil 2006, 09:50

Bonjour,

J'ai une réponse binaire, et je cherche à déterminer si elle diffère entre mes souches (n=30) sachant que pour chaque souche on a des répétition ( généralement 4 répétitions, pour lesquelle on obtient un nombre variable de réponses 0 ou 1).

J'ai donc tenté de modéliser ceci, mais je ne suis pas sûre de moi... Est-ce que la formule:
lmer(encapsul ~souche+ (1|souche:repetition), family="binomial")->mix3 est correcte par rapport à ce que je souhaite tester?


je ne comprends pas vraiment à quoi correspond la variable répétition ?

les valeurs de cette variables varient-elles pour une souche donnée ?
Ou s’agit-il juste de simples répétitions (réplicats) ? dans ce cas tu n’as pas de raison de faire intervenir un effet aléatoire.

Une autre petite chose me chagrine … c’est le nombre d'observation.
Tu as environs 120 observations (environs 4 replicats*30 souches) pour 30 paramètres estimés , je pense que c'est trop peu. En plus dans les modèles mixtes, il faut estimer aussi les paramètres associés aux effets aléatoires ...
Dix observations pour estimer un paramètre dans un modèle, c’est vraiment un minimum (voir Anderson et Burnham, 1998).

Tu peux peut-être regrouper certaines souches entre elles ?
Ou augmenter ton nombre de réplicat ?


désolé de ne pas avoir été plus constructif
en espérant t'avoir aidé un peu :)


Pierre


Référence :
Burnham K.P. & Anderson D.R. (1998) Model selection and inference. Springer, Berlin.

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 05 Juil 2006, 13:40

Bonjour,
je ne comprends pas vraiment à quoi correspond la variable répétition ?

Pour chacune de mes souches j'avais 4 (ou parfois 3) tubes dans lesquels j'ai mesuré ma variable pour plusieurs individus.
Ne maîtrisant pas la syntaxe de la fonction, j'ai essayé de me baser sur des exemples trouvés dans le site de R, et il est fort possible que je n'aie pas utilisé le bon code...

Code : Tout sélectionner

Une autre petite chose me chagrine … c’est le nombre d'observation.
Tu as environs 120 observations (environs 4 replicats*30 souches) pour 30 paramètres estimés , je pense que c'est trop peu. En plus dans les modèles mixtes, il faut estimer aussi les paramètres associés aux effets aléatoires ...

J'ai 1690 observations (=1690 individus pour lesquels j'ai un effet 0 ou 1), réparties dans 30 souches (avec des réplicats dans mes souches). J'espère que cette fois j'ai été plus claire.

Tu peux peut-être regrouper certaines souches entre elles ?
Ou augmenter ton nombre de réplicat ?

C'était une manip déjà très lourde et que je ne peux malheureusement pas refaire...
J'avais à l'époque procédé en travaillant sur les proportions d'individus présentant la réponse 0 ou 1 pour chacun des réplicats (proportions normalisées par arcsin-transformation), ce qui me permettait de faire une ANOVA, mais vu que mes pourcentages étaient très faible, la normalisation était mauvaise... et c'est pourquoi je tente de vérifier si j'obtiens un résultat identique en utilisant une régression logistique.

En tout cas merci beaucoup pour vous être penché sur mon problème.

K.

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 05 Juil 2006, 15:00

rebonjour,

# of obs: 1690, groups: souche:repetition, 111


oups ... désolé ... j'ai mal lu le listing :(

J'ai 1690 observations (=1690 individus pour lesquels j'ai un effet 0 ou 1), réparties dans 30 souches (avec des réplicats dans mes souches). J'espère que cette fois j'ai été plus claire.


ok ... disons que je n'ai rien dit :)

donc si je reprends les données, on a :
- effet souche -> fixe (S)
- effet tube -> aléatoire (r)
- individus


personnellement, j'aurais tendance à utiliser plutôt ce modèle ...
logit(yij) = mu + S(i) + r(j) + eij
avec eij correspondant à l'erreur.

Code : Tout sélectionner

lmer(encapsul ~ souche+ (1|repetition), family="binomial")


il ne me semble pas exister de lien logique de hiérarchisation
de la variabilité entre la souche et les tubes (à vérifier).



en espérant t'avoir aidé :)

@+++

Pierre

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 05 Juil 2006, 16:49

Bonjour,
Merci beaucoup pour l'aide!
lmer(encapsul ~ souche+ (1|repetition), family="binomial")


il ne me semble pas exister de lien logique de hiérarchisation
de la variabilité entre la souche et les tubes (à vérifier).


OK... mais avec ce modèle, j'ai en fait quasiment le même résultat:

Code : Tout sélectionner

> summary( lmer(encapsul ~ souche+ (1|repetition), family="binomial"))
Generalized linear mixed model fit using PQL
Formula: encapsul ~ souche + (1 | repetition)
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 1065.961 1239.801 -500.9807 1001.961
Random effects:
     Groups        Name    Variance    Std.Dev.
 repetition (Intercept)      5.1336      2.2657
# of obs: 1690, groups: repetition, 111

Estimated scale (compare to 1)  0.7750342

Fixed effects:
             Estimate Std. Error  z value Pr(>|z|) 
(Intercept) -3.464559   1.380345 -2.50992  0.01208 *
soucheA21    2.566020   1.883634  1.36227  0.17311 
soucheA22    1.680399   1.958827  0.85786  0.39097 
soucheA28    1.473629   1.888869  0.78016  0.43529 
soucheA3     0.983925   1.922013  0.51192  0.60870 
soucheA32    1.780925   1.816156  0.98060  0.32679 
soucheA34    1.169993   2.067607  0.56587  0.57148 
soucheA41    1.900044   1.829047  1.03882  0.29889 
soucheA45    0.254297   2.024243  0.12563  0.90003 
soucheA5     2.605590   1.870204  1.39321  0.16356 
soucheAc11   1.968235   1.910949  1.02998  0.30302 
soucheAc13  -0.891319   2.061694 -0.43232  0.66551 
soucheAc14   3.994075   1.929978  2.06949  0.03850 *
soucheAc16   2.593285   1.917688  1.35230  0.17628 
soucheAc21   2.074794   1.979327  1.04823  0.29453 
soucheAc3    0.673147   1.951983  0.34485  0.73021 
soucheAc9    3.825062   1.952796  1.95876  0.05014 .
soucheAcX    1.004453   1.879326  0.53448  0.59301 
soucheAcY    1.935509   1.824665  1.06075  0.28880 
soucheSF1    0.873992   1.867653  0.46796  0.63981 
soucheSF11   4.128431   1.940136  2.12791  0.03334 *
soucheSF12   1.930819   1.882131  1.02587  0.30495 
soucheSF2    3.537742   2.023653  1.74820  0.08043 .
soucheSF3    0.923809   2.023648  0.45651  0.64803 
soucheSF4    1.606224   2.107019  0.76232  0.44587 
soucheSF5   -1.453207   2.172547 -0.66890  0.50356 
soucheSF6    2.309963   1.976730  1.16858  0.24257 
soucheSF7    0.080365   1.902902  0.04223  0.96631 
soucheSF8    3.558574   1.916392  1.85691  0.06332 .
soucheSF9    1.357004   1.983106  0.68428  0.49380 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) sochA21 schA22 schA28 sochA3 schA32 schA34 schA41 schA45 sochA5 schA11 schA13 schA14 schA16 schAc21 schAc3 schAc9 schAcX schAcY schSF1
soucheA21  -0.733   ...           


Et:

Code : Tout sélectionner

> anova(
+ lmer(encapsul ~ souche+ (1|repetition), family="binomial"))
Analysis of Variance Table
       Df  Sum Sq  Mean Sq   Denom F value Pr(>F)
souche 29    0.14 0.004704 1660.00  0.0072      1


Et je ne sais toujours pas quoi conclure de ce résultat (est-ce que je dois faire confiance à la p-value de l'ANOVA?)... sinon, comment faire?

K.

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 06 Juil 2006, 08:22

bonjour,


Et je ne sais toujours pas quoi conclure de ce résultat (est-ce que je dois faire confiance à la p-value de l'ANOVA?)... sinon, comment faire?


si on a confiance en la fonction, on garde ces résultats. La discussion concernant le bug date de septembre 2005, il a peut-être été corrigé (à vérifier).
sinon, voici quelques propositions :

- regarder les résultats obtenus avec le modèle sans effet aléatoire
- regarder les autres paramétres de diagnostiques du modèle
- utiliser la fonction glmmPQL de la library MASS pour comparer.
- prévoir des simulations
- jeter un coup d'oeil du coté des modèles marginaux avec GEE.


en espérant t'avoir aidé :)

@+++

Pierre

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 06 Juil 2006, 13:03

Bonjour,

Je viens de tester la fonction glmmPQL... j'espère ne pas m'être trompée sur la syntaxe, qui est différente.
Si j'ai bien posé ma fonction, on n'a effectivement pas le même résultat!

Code : Tout sélectionner

>   summary(aov(glmmPQL(encapsul ~ souche, random = ~ 1|repetition, family = binomial)))
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
              Df  Sum Sq Mean Sq F value    Pr(>F)   
souche        29  20.573   0.709  6.7783 < 2.2e-16 ***
Residuals   1660 173.737   0.105                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Le sommaire présente aussi des différences:

Code : Tout sélectionner

 >   summary(glmmPQL(encapsul ~ souche, random = ~ 1 | repetition, family = binomial)
+ )
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
Linear mixed-effects model fit by maximum likelihood
 Data: NULL
      AIC    BIC   logLik
  9555.66 9729.5 -4745.83

Random effects:
 Formula: ~1 | repetition
        (Intercept)  Residual
StdDev:     1.75592 0.7750255

Variance function:
 Structure: fixed weights
 Formula: ~invwt
Fixed effects: encapsul ~ souche
                Value Std.Error   DF   t-value p-value
(Intercept) -3.464541  1.079332 1579 -3.209894  0.0014
soucheA21    2.565997  1.472867   81  1.742178  0.0853
soucheA22    1.680383  1.531694   81  1.097075  0.2759
soucheA28    1.473622  1.476975   81  0.997730  0.3214
soucheA3     0.983933  1.502869   81  0.654703  0.5145
soucheA32    1.780907  1.420129   81  1.254046  0.2134
soucheA34    1.170021  1.616662   81  0.723726  0.4713
soucheA41    1.900027  1.430209   81  1.328495  0.1877
soucheA45    0.254302  1.582817   81  0.160664  0.8728
soucheA5     2.605578  1.462377   81  1.781741  0.0785
soucheAc11   1.968212  1.494232   81  1.317206  0.1915
soucheAc13  -0.891285  1.612020   81 -0.552899  0.5819
soucheAc14   3.994030  1.509096   81  2.646637  0.0098
soucheAc16   2.593229  1.499494   81  1.729403  0.0875
soucheAc21   2.074797  1.547697   81  1.340571  0.1838
soucheAc3    0.673158  1.526300   81  0.441039  0.6604
soucheAc9    3.825027  1.526924   81  2.505053  0.0142
soucheAcX    1.004448  1.469513   81  0.683524  0.4962
soucheAcY    1.935492  1.426781   81  1.356545  0.1787
soucheSF1    0.873989  1.460383   81  0.598466  0.5512
soucheSF11   4.128393  1.517023   81  2.721377  0.0080
soucheSF12   1.930794  1.471699   81  1.311949  0.1932
soucheSF2    3.537704  1.582360   81  2.235714  0.0281
soucheSF3    0.923807  1.582358   81  0.583817  0.5610
soucheSF4    1.606194  1.647539   81  0.974905  0.3325
soucheSF5   -1.453187  1.698720   81 -0.855460  0.3948
soucheSF6    2.309944  1.545696   81  1.494436  0.1389
soucheSF7    0.080378  1.487921   81  0.054020  0.9571
soucheSF8    3.558551  1.498481   81  2.374772  0.0199
soucheSF9    1.357008  1.550651   81  0.875122  0.3841


C'est marrant comme un problème qui me semblait au départ simple pose en fait des problèmes statistiques que je n'aurais même pas soupçonnés...

K.

Pierre Bady
Messages : 405
Enregistré le : 02 Mai 2006, 07:46

Messagepar Pierre Bady » 07 Juil 2006, 10:31

bonjour,

C'est marrant comme un problème qui me semblait au départ simple pose en fait des problèmes statistiques que je n'aurais même pas soupçonnés...


ce type de situation est toujours délicate. Qui croire ? Que croire ?
dans ton cas, il plane quand même un doute sur la fonction lmer :
Malheureusement, si je m'en reporte à cette discussion: http://finzi.psych.upenn.edu/R/Rhelp02a ... 61884.html je ne sais pas quel crédit je dois accorder au résultat de l'anova.


lorsque, j'ai un doute sur une fonction, j'ai tendance à utiliser d'autres fonctions (il y a qlq doublons dans R ... et heureusement). ça permet (il me semble) d'obtenir des résultats plus fiables et moins dépendant du mode de programmation de la procédure utilisée .

et si vraiment le doute persiste ... ben ... j'utilise aussi d'autres logiciels de stat. je ne voudrais pas blasphémer sur ce forum, mais la justesse de mes analyses passe avant l'utilisation de R. :)


en espérant t'avoir aidé un peu :)

@+++

Pierre

PS: Au passage, il y a aussi la commande glmmML (disponible dans la librairie du même nom) qui te permet de réaliser des glmm avec des effets aléatoires sur l’intercept.

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 10 Juil 2006, 15:04

Bonjour,
Et merci beaucoup de m'avoir aidé à débroussailler mon problème!
J'ai donc aussi essayé la fonction glmmML, qui me donne des résultats encore un peu différents:


Code : Tout sélectionner

> glmmML(encapsul ~ souche, cluster= repetition, family="binomial")

Call:  glmmML(formula = encapsul ~ souche, family = "binomial", cluster = repetition)


               coef se(coef)       z Pr(>|z|)
(Intercept) -4.4310    1.523 -2.9096 0.003620
soucheA21    3.8825    1.727  2.2477 0.024600
soucheA22    2.4549    1.591  1.5430 0.123000
soucheA28    2.7057    1.678  1.6128 0.107000
soucheA3     0.9443    1.827  0.5167 0.605000
soucheA32    2.2575    1.537  1.4686 0.142000
soucheA34    0.7112    1.702  0.4180 0.676000
soucheA41    3.0590    1.609  1.9013 0.057300
soucheA45    1.8133    1.602  1.1315 0.258000
soucheA5     3.4715    1.917  1.8104 0.070200
soucheAc11   2.6736    1.699  1.5738 0.116000
soucheAc13   1.0879    1.611  0.6755 0.499000
soucheAc14   5.3103    1.885  2.8170 0.004850
soucheAc16   4.5930    1.894  2.4254 0.015300
soucheAc21   2.3090    1.547  1.4926 0.136000
soucheAc3    0.5420    1.685  0.3216 0.748000
soucheAc9    5.2100    2.026  2.5714 0.010100
soucheAcX    1.8544    1.638  1.1323 0.258000
soucheAcY    1.6621    1.546  1.0749 0.282000
soucheSF1    1.8311    1.635  1.1198 0.263000
soucheSF11   5.8861    2.234  2.6350 0.008410
soucheSF12   6.0159    1.801  3.3395 0.000839
soucheSF2    3.9444    1.657  2.3800 0.017300
soucheSF3    1.6741    1.591  1.0524 0.293000
soucheSF4    1.8698    2.124  0.8805 0.379000
soucheSF5   -1.5838    2.184 -0.7254 0.468000
soucheSF6    2.2766    1.768  1.2874 0.198000
soucheSF7    2.3205    1.635  1.4194 0.156000
soucheSF8    3.3138    1.664  1.9912 0.046500
soucheSF9    0.5698    1.569  0.3631 0.717000

Standard deviation in mixing distribution:  2.515
Std. Error:                                 0.4577

Je n'ai pas trouvé de test associé à la fonction pour obtenir un p-valeur globale...

Je vais encore tester ceci avec quelquesautres fonctions qui apparemment peuvent être utilisées:
http://www.biostat.umn.edu/~nali/teachi ... 08GLMM.pdf
Fitting GLMM: Software
• SAS:
– PROC NLMIXED uses Gauss-Hermite or Adaptive Quadrature methods.
– PROC GLIMMIX uses linearization idea to fit penalized/marginal quasilikelihood
models with Gaussian random eects.
• R
– glmmPQL (MASS): penalized quasi-likelihood, allows the use of an
additional correlation structure.
– GLMM (lme4 a newer reimplementation of nlme): similar to (same
as?) glmmPQL.
– glmmML (glmmML): maximum likelihood using Gauss-Hermite quadrature,
only can fit models with a random intercept. (work in progress)
– glmm (repeated): Gauss-Hermite quadrature, models with random
intercept.
– gnlmix (repeated): non-linear regression with mixed random eects
for the location parameters. Non-Gaussian mixing distributions are
allowed.
– glmm (GLMMGibbs): Gibbs sampling.
• BUGS (OpenBUGS/WinBUGS)


Bon, bref, il me semble que les résultats concordent pour dire que l'on a des souches qui ont des caractéristiques significativement différentes des autres, mais ça ne sera pas évident pour moi à présenter dans un article tout ça...

K.

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 03 Aoû 2006, 16:41

Bonjour,

Comme la méthode m'intéresse, que j'en aurai peut-être besoin d'autre fois, et que je n'aime pas ne pas comprendre certaines choses, j'ai encore des interrogations, qui touchent à l'interprétation des résultats des tests.

Avec la fonction lmer ce résultat est disponible:

Code : Tout sélectionner

> anova( lmer(encapsul ~ souche+ (1|repetition), family="binomial") )
Analysis of Variance Table
       Df  Sum Sq  Mean Sq   Denom F value Pr(>F)
souche 29    0.14 0.004704 1660.00  0.0072      1


Avec la fonction glmmPQ je peux obtenir:

Code : Tout sélectionner

> anova(glmmPQL(encapsul ~ souche, random = ~ 1|repetition, family = binomial))
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
            numDF denDF  F-value p-value
(Intercept)     1  1579 78.52857  <.0001
souche         29    81  1.47565  0.0885


Ici j'ai donc un résultat sur l'intercept et sur l'effet "souche". Ce test n'est donc pas équivalent à celui réalisé avec la fonction lmer?

Je peux aussi obtenir ce résultat:

Code : Tout sélectionner

> summary(aov(glmmPQL(encapsul ~ souche, random = ~ 1|repetition, family = binomial)))
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
iteration 6
              Df  Sum Sq Mean Sq F value    Pr(>F)   
souche        29  20.573   0.709  6.7783 < 2.2e-16 ***
Residuals   1660 173.737   0.105                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Donc, pour le test glmmPQL, ça signifierait que l'ajustement par le modèle (globalement) est bon, mais que l'effet de la souche n'est pas significatif???


Là je me sens un peu stupide, car jusqu'à maintenant sur des tests plus classiques, j'avais souvent utilisé alternativement aov et anova en considérant que c'était équivalent (et obtenu des résultats similaires).

Donc, au risque de paraître complètement stupide, je viens chercher ici quelques explications, s'il y en a.

Merci d'avance.
K

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 04 Aoû 2006, 13:29

Karine,

glmmPQL est un "bricolage" qui utilise la fonction lme pour réaliser l'essentiel des calculs. Son auteur (B Ripley) n'a pas écrit de méthode particulière pour les modèles ajustés avec glmmPQL. Je serais circonspect sur les résultats obtenus avec anova. Note d'ailleurs qu'avec les objets de classe lme, il faut utiliser l'argument type = "m" pour avoir un test de Wald correct ==> à essayer dans ton exemple pour voir si ça change qque chose. Pour plus de détails sur glmmPQL, voir la v4 du bouquin MASS (Venables et Ripley).

Je ferais donc plutôt confiance aux résultats de lmer qu'à ceux de glmmPQL. Je vois que tu as aussi utilisé le package glmmML. Une nouvelle version est sortie il y a qques jours avec une nouvelle méthode d'estimation: quadrature de Gauss-Hermite. C'est le top actuel (hors MCMC), prévu mais pas encore disponible avec lmer. Les résultats devraient être moins biaisés pour l'effet aléatoire. Tu devrais vérifier que sa variance est la plus petite avec PQL (glmmPQL ou lmer avec method = "PQL"), la plus grande avec Gauss-Hermite et intermédiaire avec Laplace.

Voir aussi les archives de R-Help pour de nombreuses discussions sur les tests de Wald avec les modèles GLMM, notamment les interventions de Douglas Bates (auteur de lmer), et les pbs de calcul des ddl. Avec une structure aléatoire aussi simple que celle que tu utilises, je ne pense pas que ça soit un vrai pb.

Il est toutefois préférable d'utiliser un test du rapport des vraisemblances pour évaluer la signification stat d'un effet fixe... ce qui pose un autre pb avec les GLMM, dans la mesure où la vraisemblance n'est pas calculée de manière exacte. Il faudrait utiliser l'option method = "Laplace" qui fournit des estimateurs plus précis qu'avec PQL.

Dans ton exemple, cela donnerait:

Code : Tout sélectionner

fm1 <- lmer(encapsul ~ souche+ (1|repetition), family="binomial", method = "Laplace")
fm2 <- lmer(encapsul ~ 1 + (1|repetition), family="binomial", method = "Laplace")
anova(fm1, fm2)


Le package glmmML propose également une méthode bootstrap pour évaluer la signification de l'effet aléatoire et des effets fixes. Je ne l'ai pas encore essayée.

Tu peux aussi faire le test de Wald "à la main" en utilisant la fonction wald.test du package aod, ce qui te permettrait d'évaluer la sensibilité du résultat du test aux ddl. Voir l'aide de cette fonction et utiliser les fonctions vcov et fixef pour extraire la matrice de var-cov et les coef des effets fixes.

Amicalement,

Renaud

Karine Poitrineau
Messages : 18
Enregistré le : 30 Mai 2006, 19:37

Messagepar Karine Poitrineau » 08 Aoû 2006, 15:36

Bonjour Renaud,
Merci pour tes explications et conseils.

Je progresse doucement dans ma compréhension des modèles logistiques mixtes... :D

Mais voilà, R me fait encore quelques misères, il refuse de me calculer les paramètres de ce modèle:

Code : Tout sélectionner

> fm2 <- lmer(encapsul ~ 1 + (1|repetition), family="binomial", method = "Laplace")
Message d'avis :
optim or nlminb returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
 in: LMEopt(x = mer, value = cv)


Ca ne me peine pas tant que ça, vu que j'ai réussi à me faire à l'idée de présenter essentiellement une bête ANOVA hiérarchique sur des données arcsin-transformées et juste présenter le résultat d'un modèle logistique mixte qui ne sera de toute façon pas meilleur (et au vu de mes données, bon, finalement je me dis que c'était prévisible: pour l'ANOVA hiérarchique, j'ai une p-value de 0.9 pour l'effet de mes lignées :roll: -> avec lmer le F de l'anova donne p= 0.933). Mais au moins je saurai un peu mieux me débrouiller pour des futures analyses (si je continue à travailler sur des problèmes de sex-ratio, je risque de retomber sur ce type d'analyses).

K

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 09 Aoû 2006, 16:44

Bonjour Karine,

C'est sans doute un problème dans tes données, comme une très faible variance de l'effet aléatoire, par exemple. Essaie avec les données transformées (arcsin) et lme pour voir si ça plante aussi.

Amitiés,

Renaud


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité