GLMM

Bonjour,

J'ai un petit souci avec la fonction glmer mais pour plus de compréhension, voici le cheminement.

Mon jeu de données se compose de variables réponses qui sont des comptes de loirs dans des nichoirs. Les comptes ont été transformés via racine carré.
Mes variables explicatives sont des différentes mesures de l'habitat. Cependant, vu que le nombre de variables étaient bien trop nombreuses, nous avons regroupés ces variables en différentes groupes caractéristiques (sol, climat, ...) et effectués des acp. Donc mon fichier contient les variables réponses (sqrt transformées) et les variables explicatives (coordonnées de variance des axes retenus).

Ensuite, je construit mon glmer avec mes variables et quelques interactions (uniquement celle qui sont plausibles). Nous avons ajoutés la taille des nichoirs (petit ou grand) comme variable aléatoire. Je lance un dregde() pour sélectionner le meilleur modèle via AICc (avec REML=FALSE). Et je reconstruit mon modèle avec les variables sélectionnées et REML=TRUE. Voilà ce que j'obtiens :

Code:

> model2<- glmer(sqrtCapt ~ Canopy + NestHeight + NestSpecies + WoodDiversity + Woodland + Canopy:Woodland + NestHeight:NestSpecies + (1|BoxSize), data=data, family=gaussian)
> summary(model2)
Linear mixed model fit by REML
Formula: sqrtCapt ~ Canopy + NestHeight + NestSpecies + WoodDiversity +      Woodland + Canopy:Woodland + NestHeight:NestSpecies + (1 |      BoxSize)
   Data: data
   AIC   BIC logLik deviance REMLdev
 166.5 185.6 -73.23    134.7   146.5
Random effects:
 Groups   Name        Variance Std.Dev.
 BoxSize  (Intercept) 0.62503  0.79059
 Residual             0.91787  0.95805
Number of obs: 50, groups: BoxSize, 2

Fixed effects:
                       Estimate Std. Error t value
(Intercept)              2.0896     0.5930   3.524
Canopy                  -0.3745     0.2122  -1.765
NestHeight               0.4713     0.2124   2.218
NestSpecies              0.0895     0.1719   0.521
WoodDiversity            0.4332     0.1628   2.661
Woodland                -0.5307     0.2541  -2.088
Canopy:Woodland         -0.5113     0.2530  -2.021
NestHeight:NestSpecies   0.5411     0.1764   3.068

Correlation of Fixed Effects:
            (Intr) Canopy NstHgh NstSpc WdDvrs Wodlnd Cnpy:W
Canopy       0.026                                          
NestHeight   0.041 -0.326                                    
NestSpecies -0.096 -0.138 -0.046                            
WoodDivrsty  0.037  0.182  0.123  0.252                      
Woodland     0.039 -0.571  0.743 -0.074 -0.022              
Canpy:Wdlnd -0.244 -0.107 -0.167  0.395 -0.153 -0.161        
NstHght:NsS  0.047 -0.379 -0.105  0.163 -0.245 -0.035 -0.193

Comment vous pouvez le voir, j'ai des t value et non des p-value. Et là, je suis coincé.
Je voudrais connaître les effets de la variable aléatoire (BoxSize) et des autre variables sur ma variable réponse.

J'ai par la suite trouvé la fonction pvals.fnc() du package languageR qui calcule des p-value et autres utilisant Markov Chain Monte Carlo test (MCMC).

Dans la litterature, il serait préférable d'utiliser MCMC pour calculer les p-value qu'un test de Wald après une GLMM.

Donc voici la sortie R de la fonction

Code:

> pvals.fnc(model2)
$fixed
                       Estimate MCMCmean HPD95lower HPD95upper  pMCMC Pr(>|t|)
(Intercept)              2.0896   2.0942     0.0539     4.4158 0.0564   0.0010
Canopy                  -0.3745  -0.3744    -0.7960     0.0609 0.0854   0.0849
NestHeight               0.4713   0.4676     0.0335     0.8941 0.0338   0.0320
NestSpecies              0.0895   0.0864    -0.2778     0.4159 0.6204   0.6054
WoodDiversity            0.4332   0.4346     0.1104     0.7664 0.0104   0.0110
Woodland                -0.5307  -0.5344    -1.0523    -0.0222 0.0448   0.0429
Canopy:Woodland         -0.5113  -0.5148    -1.0321    -0.0085 0.0508   0.0497
NestHeight:NestSpecies   0.5411   0.5394     0.1841     0.8943 0.0054   0.0038

$random
    Groups        Name Std.Dev. MCMCmedian MCMCmean HPD95lower HPD95upper
1  BoxSize (Intercept)   0.7906     0.9346   1.1895     0.1466     3.0004
2 Residual               0.9581     0.9704   0.9817     0.7731     1.2037

D'où les questions suivantes:

Q1: La méthode et les fonctions/packages vous semblent-t-ils correctes ou est-ce que j'ai manqué quelque chose ?

Q2: J'obtiens des p-value pour les effets fixes mais compte tenue que ces p-value sont recalculés, est-ce que je peux les utiliser pour déterminer les effets de ces variables sur mon nombre de capture, et ce , pour une publie ?

Q3: Question un peu plus naïve, comment je fais pour déterminer si l'effet de ma variable BoxSize est significatif ou non ?
J'ai entendu parlé du LR test pour calculer l'effet de cette variable mais je ne sais pas quelle fonction utiliser.

Merci
Cordialement

Mélanie

Bonjour,

Q1 : a priori ta démarche semble correcte d'après ce que j'ai pu voir dans les aides de la fonction etc. Si j'ai bien compris il existe également une alternative, regarde l'aide de ?pvalues du package lme4 :

for fixed effects, F tests via Kenward-Roger approximation using ‘KRmodcomp’ from the ‘pbkrtest’ package (MC)

Cela dit : 1 - je pense que cette aide est antérieur à la fonction p.vals.fnc que tu utilises, et 2 - il s'agit d'un F test... En tout cas, l'utilisation de MCMC pour le calcul des p.values dans les modèles mixte étant préconisé par Bates, la question est bien de savoir si ce sont ces valeurs que p.vals.fnc() calcule, ce qui n'est pas limpide dans leur aide. Je pense tout de même qu'il l'aurait précisé dans le cas contraire étant donné que les auteurs prétendent remplacer la fonction mcmcsamp implémenté à l'origine dans lme4 (et qui faisait bien du Monte Carlo) et que en en-tête de l'aide on a :

Compute p-values and MCMC confidence intervals for mixed models

Q2 : J'imagine que oui, les p.values de base sont critiquées car elles sont basées sur un df fixe (voir le lien que je t'ai donné sur le forum des utilisateurs de R). Le calcul a posteriori que tu fais corrige cela, sensément, et elles sont du coup plus adaptées. Le mieux est peut être de poser directement la question à Bates.

Q3 : Une liste des possibilité est donnée dans l'aide de la fonction pvalues citée plus haut, notamment :

for random effects, simulation tests via the ‘RLRsim’ package (MC,*)

Je tiens à te rappeler que je suis néophyte en matière de modèles mixtes, mais étant particulièrement intéressé par ces modèles et confronté aux mêmes difficultés, je cherche à alimenter la discussion et à comprendre, donc n'hésite pas à partager tes avancées sur la voie de la sagesse du modèle mixte

Bonne journée!

Aussi, ce lien risque de t'intéresser si tu ne l'a pas déjà parcouru : https://stat.ethz.ch/pipermail/r-help/2006-August/110736.html

Bonjour,

La vraie question est plutôt quel est l'intérêt de calculer ces p-values étant donné que tu as fait une procédure de sélection de modèle (même si passer par du data dredging est pas ce qu'il y a de plus scientifique) et que tu as la valeur des paramètres ?
Le seul paramètre qui pose question selon son IC via MCMC c'est NestSpecies mais la valeur du paramètre (l'effet) est 5 à 6 fois moindre que les autres facteurs retenus. Donc certes selon les données, le meilleur modèle comprend ce paramètre mais il est quantitativement peu important au regard des autres modèles. Il serait plus intéressant d'aller regarder les poids d'Akaike pour chacune des variables pour voir celles qui sont retenues le plus souvent dans la procédure de dredging.

Pour ta question 3, si je la comprends bien c'est est ce que la structure aléatoire apporte quelque chose au modèle ou pas ?
Pour cela, tu peux passer par un gls qui ne contient pas la structure aléatoire en utilisant une estimation au MV et faire de même sur le modèle mixte. Tu peux alors comparer ces 2 modèles par AICc et tu auras une indication sur l'intérêt de la structure aléatoire. Par contre, j'ai un petit doute s'il faut faire cette comparaison à partir du modèle complet ou si tu peux partir sur le modèle déjà sélectionné via dredging. Le livre de Zuur et al Mixed Effects Models and Extensions in Ecology with R détaille cette procédure.

HTH

Nik

Merci beaucoup pour vos réponses,

En attendant de continuer grâce à vos remarques, j'aimerais d'abord revenir sur un point qu'a soulevé Nick et qui est correct sur ma sélection de modèle.

Je suis passé par du data dredging et j'ai pris le meilleur modèle, ce qui n'est pas forcément la bonne solution.

Q1: Quelle autre méthode de sélection de variables me conseillerez-vous ?

Pour ma part, j'aimerais bien selectionner un sous-modèle sur la base de l'AICc comme le fait la fontion step(), mais cette fonction ne fonctionne malheureusement pas pour les GLMM.

Je suis donc passé par la fonction dredge(). Si il y a une autre façon de faire, je suis preneuse.

Un autre procédé, sans doute plus correct pour selectionner le meilleur modèle, passerait par un "model-averaging approach" (décrit notamment dans l'article de Grueber et al. 2010. Multimodel inference in ecology and evolution: challenges and solutions).

Pour faire simple, on prend le modèle global, on effectue une standardisation du modèle, on lance la fonction dredge() et on fait une compile des meilleurs modèles.

J'ai vu cette procédure dans quelques autres articles mais il y a un point sur lequel je reste indecise.
Pour sélectionner les meilleurs modèles, il y aurait 4 façons de faire, soit on prend tous les modèles avec un delta inférieur à 2, 6 ou 10 soit on prend tous les modèles jusqu'à ce que la somme du poids des modèles atteigne 95%.

D'où cette question:

Q2: Comment savoir quel mode de sélection choisir ? De quoi cela dépend-t-il ?
Pour ma part, le dredge() compare quelques 300 modèles.

Merci
Mélanie

Quelle autre méthode de sélection de variables me conseillerez-vous ?

Le fait de partir sur du data dredging indique que tu pars à l'aveuglette. Or une sélection de modèle n'est pas un remède miracle à l'absence de connaissance sur un sujet. On ne peut pas prendre tous les paramètres possibles et imaginables (et surtout disponibles) et attendre que le résultat soit cohérent écologiquement parlant. Plus tu ajoutes des paramètres plus tu as de chance qu'un modèle sorte du processus de sélection alors qu'il est éminemment faux.
Le processus de sélection de modèle doit donc correspondre avant tout à un processus de sélection d'hypothèse. Ce qui suppose une construction des modèles candidats selon des hypothèses raisonnables. Donc il faut passer par une définition des modèles quasiment à la main.
Dans ton cas, le fait de passer par des sorties d'ACP ne va pas te simplifier la tâche car c'est très compliqué à interpréter et donc difficile d'émettre des hypothèses à priori. Par contre on peut voir dans tes effets fixe, qu'il existe au moins une corrélation assez forte donc ce genre de chose ne devrait pas être entré dans un modèle ni dans un processus de sélection de modèle.

Pour ce qui est du critère de choix pour le modèle averaging, c'est vrai qu'il n'y a pas de consensus là dessus car finalement c'est un champs de développement assez récent. Mais pour ma part, je suis assez réfractaire à l'idée des 95% car cela peut amener des modèles dont le poids est très faible (donc très peu supportés par les données) et on est tout de même pas sur la qualification de paramètres d'une variable aléatoire alors ressortir l'histoire de l'intervalle des 95% me parait douteuse. Ensuite, pour ce qui est de la valeur de delta AIC, il faut tout de même relever qu'on est sur une échelle log donc une différence de qq unités d'AIC peut déjà représenter une différence très importante entre les modèles (distance relative avec le "vrai" modèle). La valeur de 2 me semble donc correct voire 4 mais au delà je pense qu'on commence à rentrer dans des échelles de tolérance assez larges...

Nik

Merci Nick pour tes éclaircissements !

Le fait de passer par des sorties d'acp est un choix calculé dans notre cas (au moins sur le papier, je ne suis pas encore passé à l'interprétation). Pour faire simple, chaque axe n'est en fait expliqué que par une ou deux variables dont certaines sont des facteurs. C'est donc un choix délibéré d'utiliser les sorties d'acp pour limiter la sur-dispersion des modèles via ces variables facteurs.

Mélanie

Pour alimenter la discussion et partager mes avancées :

Une solution pour éviter de passer par du data dredging et surtout par la méthode du model-averaging est d'utiliser le package lmerTest.

Non tu tournes en rond. Tu repars sur du test d'hypothèse basé sur une distribution F. Le package lmertest retranscrit les fonctions utilisées sous SAS et sérieusement critiquées par Douglas Bates.
Le model averaging n'est pas à éviter, c'est une méthode théoriquement très intéressante et bien pensée. Le data dredging n'est pas forcément faux non plus mais ce n'est pas une démarche scientifique, en tout cas pas dans une science qui repose beaucoup sur des hypothèses à tester.

Nik

Bonjour,

Je continue à partager mes avancées comme l'a demandé Alexandre.

J'ai choisit d'utiliser la méthode du model-averaging avec le package MuMIn pour voir les effets de mes variables.

J'ai aussi abandonner les modèles mixtes et je me suis redirigé vers les GLM. En effet, Benjamin Bolker a évoqué le fait qu'il est plus raisonnable d'utiliser une variable aléatoire qui contient au moins 5-6 niveaux, sinon les GLMM ne marcheraient pas proprement. Dans mon cas, je n'ai que 2 niveaux.

Merci encore Nik pour tes conseils et éclaircissements

Mélanie

Bonjour,

Merci pour ce partage, je n'ai plus répondu depuis un moment mais soit sûre que je garde un œil sur cette discussion. Les modèles mixtes m'intéressent dans le cadre de données longitudinales. Je n'ai pas encore eu le temps de m'y coller et me suis rabattu sur les GLM en simplifiant mes données le temps de mieux comprendre ces méthodes et ne pas faire n'importe quoi...

Aussi, je ne suis pas sûr de comprendre ce qu'est "une variable à 5-6 niveaux". As-tu un lien vers l'article (ou la source en tout cas) de Benjamin Bolker?

Cordialement.
Alex

L'article de Ben Bolker que tu trouveras sur internet sans pb:

Bolker et al. 2008. Generalized linear mixed models: a practical guide for ecology and evolution. Trends in Ecology and Evolution. 24(3):127-135.

Bon courage
Mélanie