Modèle linéaire mixte - variables aléatoires.

Bonjour,

Je travaille sur des données longitudinales et on vient de me conseiller l'utilisation de modèle linéaire mixte. J'ai cherché des informations à droites à gauche et expérimenté des tutoriaux mais j'ai quelques difficultés encore.

Mes données sont les mesures d'un paramètre respiratoires chez des personnes atteintes d'une maladie. Chaque individus ne présentent pas le même nombre d'examen au cours de leurs suivi, ni le même intervalle temps entre deux examens.

Je cherche à savoir si le paramètre en question évolue en fonction du temps, dans mon échantillon.

Mes variables :

variable réponse, quantitative, « Pcv »
variables à effet fixe :   - « forme » qualitatif : la forme de la maladie, codé 1 / 2 /3 / d.
                               - « sexe », qualitatif : codé 0/ 1
                               -« age », quantitatif : l'âge du patient lors de l'examen.
variables aléatoires : - « id », qualitatif : identité du patient.

library(lme4)
model1<- lmer(Pcvassis ~ age+sexe+forme+(1|id), data=fin)
model2 <- lmer(Pcvassis ~ age+sexe+forme+(1+age|id), data=fin)

J'ai quelques difficultés à définir la forme de la partie variable entre le model1 et model2 (ou autre?).
Avec le model2 j'essaye d'introduire une pente aléatoire mais je ne sais pas si cela suit la bonne formulation.
Qu'en pensez vous ?

Une fois le modèle choisi je pensais le comparer avec un modèle emboîté ne comportant pas la variable « age » à partir d'une anova.

Je vous remercie d'avance pour votre aide,
Cordialement.

Benjamin

1) Sur le modèle fixe, tu supposes aucune interaction entre les variables explicatives, donc notamment que la forme de la maladie et le sexe sont indépendants ou que l'effet âge ne dépend pas de la forme de la maladie. En es-tu sûr ou utilises-tu un modèle simplifié car tu considères n'avoir pas assez de données pour un modèle plus complexe / puissant ? Il existe des approches pour chercher le modèle le plus adapté.

2)

Avec le model2 j'essaye d'introduire une pente aléatoire mais je ne sais pas si cela suit la bonne formulation.

Si, la formule est exact. Cependant, tu n'as pas besoin d'écrire (1 + age)|id, le '1+' est implicite dans ce cas.

3) On peut utiliser la fonction "anova" pour comparer deux modèles fixes emboités (dans ce cas, il faut utiliser comme méthode de calcul du modèle ML et non REML, mais anova réécrit dans lme4 fait elle-même le recalcul si nécessaire), par contre "anova" ne peut pas servir pour comparer deux modèles aléatoires emboîtés. En effet, "anova" calcule la significativité de la différence à partir des degrés de liberté or personne ne sait définir ce que sont les degrés de libertés d'un modèle aléatoire. La procédure de comparaison est plus complexe. Voici une procédure que j'ai trouvée pour comparer les effets aléatoires, elle est basée sur une méthode rééchantillonnage (http://www.stat.wisc.edu/~ane/st572/notes/index.html ) :

Merci pour ta réponse.

1) Hmm en effet.
Il y a des arguments pour que la variables sexe n'ait pas tant d'importance que ça sur ma variable prédite.

Par contre pour l'interaction forme*age c'est l'inconnu.
Pour 2 des 4 catégories de la variable « forme » j'ai des effectifs relativement faibles (moins de 5 sur 85 individus ;  soit < 20 / 324 examens).


3) Donc si je résume la fonction « anova » me permet de comparer les parties fixes de modèles mixtes emboîtés,
tandis que si je cherche à comparer deux modèles mixtes emboîtés dont la partie aléatoires diffère il faut que je passe par d'autres étapes.


Si oui je peux vérifier la présence d'interaction forme*age et l’intérêt de la variable age:

modelsimple <- lmer(Pcvassis ~ age+sexe+forme+(age|id), data=fin)
modelintera <- lmer(Pcvassis ~ age+sexe+forme+ age*forme +(age|id), data=fin)
modeltest <- lmer(Pcvassis ~ sexe+forme+(age|id), data=fin)

anova(modelsimple, modelintera, modeltest)

Sans pour autant savoir si la partie aléatoire du modèle est pertinente.
(En attendant d'avoir intégré la procédure que tu m'as envoyé en lien).

En ce qui concerne la vérification de la pertinence de l'interaction et l'intérêt de la variable âge :
1) Dans anova, il faut que tu listes tes modèles dans l'ordre, ou l'ordre inverse, d'emboîtement. Tu dois donc écrire :

Code:

anova( modelintera, modelsimple, modeltest)
# ou
anova( modeltest, modelsimple, modelintera)

De plus, formellement cette comparaison se fait à partir des estimations par ML tandis que le calcul du modèle utilise REML.

2) Il existe une méthode plus simple pour chercher le modèle 'optimal' par la fonction stepAIC. Le code serait, en reprenant tes modèles :

Code:

modelintera <- lmer( Pcvassis ~ age + sexe + forme + age*forme + (age|id), data=fin, REML=FALSE)
modeloptL <- stepAIC( modelintera , scope=list( lower=~ sexe + forme))
# Formule du modèle optimal
formula( modeloptL)
# Estimation du modèle
est.modeleoptL <- update( modeleoptL, REML=TRUE)
 

on peut aussi écrire :

Code:

modeltest <- lmer( Pcvassis ~ sexe + forme + (age|id), data=fin, REML=FALSE)
modeloptU <- stepAIC( modeltest , scope=list( lower=~ age + sexe + forme + age*forme))
# Formule du modèle optimal
formula( modeloptU)
# Estimation du modèle
modeleoptU <- update( modeleoptU, REML=TRUE)
 

Il arrive en effet qu'avec stepAIC, les deux modèles soient différents. Cela provient du fait que l'AIC est dérivé du maximum de vraisemblance corrigé du nombre de degré de liberté. Si c'est le cas, ce que je conseille toujours est de comparer alors les deux modèles par anova, et s'ils ne sont pas significativement différents de prendre le plus simple (principe du rasoir d'Occam). S'ils sont significativement différents mais que la différence entre les deux AIC est faible par rapport à la valeur absolue des AIC (par exemple, inférieure au %), je conseille malgré tout d'utiliser le modèle le plus simple. En effet, l'AIC à une interprétation simple :
exp( (AIC1 - AIC2) / 2), avec AIC2 > AIC1
[*]est la vraisemblance relative entre les deux modèles et peut être interprété comme étant proportionnel à la probabilité que le second modèle minimise la perte estimée d'informations par rapport au premier. De ce fait, si les AIC sont grands en valeur absolue et
- positifs, de toute façon l'ajustement n'est pas très bon ;
- négatifs, l'ajustement est très bon et la différence peut être un artefact car se situer dans le bruit.

[*]Plus l'AIC est faible, meilleur est le modèle.

1) En utilisant anova:

Code:

Data: fin
Models:
modeltest: Pcvassis ~ sexe + forme + (age | id)
modelsimple: Pcvassis ~ age + sexe + forme + (age | id)
modelintera: Pcvassis ~ age + sexe + forme + age * forme + (age | id)
            Df     AIC     BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
modeltest    9 -776.67 -743.04 397.33  -794.67                          
modelsimple 10 -774.92 -737.55 397.46  -794.92 0.2482      1     0.6184  
modelintera 13 -777.86 -729.28 401.93  -803.86 8.9400      3     0.0301 *

J'ai donc une différence significative entre modèles mais une variation absolue qui est négligeable.

L'AIC est négatif.
Y a t'il des références de grandeur / des règles de décisions pour les AIC ?
Par exemple: un AIC à -775, ça te donne une information ou c'est uniquement sa comparaison à ceux des autres modèles qui a un sens?

2) Message d'erreur en utilisant ton code, avec la library(MASS).

Code:

> modeloptL <- stepAIC( modelintera , scope=list( lower=~ sexe + forme))
Error in `$<-`(`*tmp*`, formula, value = Terms) :
  no method for assigning subsets of this S4 class

En tout cas, un grand merci pour ton aide.

Pour la première question, une grande valeur négative signifie un bon ajustement. La différence relative entre les deux AIC est de l'ordre de 0.15 %. De plus, c'est modeltest qui a le meilleur BIC, indice qui n'a pas vraiment d'interprétation probabiliste contrairement à l'AIC mais qui est interprété comme un indice de simplicité du modèle compte-tenu des degrés de liberté du modèle, de sa déviance et du nombre de paramètres à estimer. Il prévient du sur-ajustement un peu à la manière du R2 ajusté du modèle linéaire classique à effets fixes.


Pour ta seconde question, au lieu d'utiliser lmer, utilise lme du package nlme. A ce moment, l'argument pour dire que tu utilise ML au leiu de REML est method="ML". Puis dans le update, c'est method="REML". lmer est une réécriture avec extensions de lme par la même équipe. Et pour être sûr de la pertinence du modèle je partirai plutôt de la formule fixe :

Code:

Pcvassis ~ (age + sexe + forme)^2

que de celle que tu utilises dans modelintra.

Sinon, globalement au vu de tes résultats, je ne serais pas choqué personnellement que tu partes sur modeltest mais il faudra discuter ensuite (compromis puissance / robustesse, ...) car cliniquement il est quand même peu satisfaisant à moins que tes groupes soient bien équilibrés en âge et que les âges soient resserrés.

Deux réponses différentes en fonction de la formule avec laquelle je démarre le stepAIC.
Au final des AIC très proches, une différence relative < 0.0%.
Par contre, si j'ai bien compris, le BIC du modèle "simple" devrait être + favorable.

Code:

Start:  AIC=-773.28
Pcvassis ~ (age + sexe + forme)^2

             Df     AIC
- sexe:forme  3 -778.01
- age:sexe    1 -773.37
<none>          -773.28
- age:forme   3 -769.11

Step:  AIC=-778.01
Pcvassis ~ age + sexe + forme + age:sexe + age:forme

            Df     AIC
<none>         -778.01
- age:sexe   1 -777.86
- age:forme  3 -773.93

Code:

Start:  AIC=-774.92
Pcvassis ~ age + sexe + forme

       Df     AIC
- age   1 -776.67
<none>    -774.92

Step:  AIC=-776.67
Pcvassis ~ sexe + forme

Hmmm, probablement non.
C'est une maladie rare (tant mieux) donc j'ai peu d'individus au total, et une répartition de la variable "forme" qui n'est pas équilibrée.
C'est une maladie grave donc j'ai une forte concentration de personnes qui sont suivis sur un intervalle 20-40 ans, beaucoup moins par la suite.

Existe t'il un moyen de différencier ces modèles autrement que par leur AIC ?
Ou mettre en rapport avec une mesure de dispersion de ma variable age ?

L'algorithme utilisé dans le modèle mixte est peu sensible au déséquilibre des groupes contrairement à celui du modèle linéaire à effets fixes. Donc, le fait que la variable forme ne soit pas équilibrée n'est pas un problème en soi, sauf si un des niveaux ne contient que quelques individus (par exemple 2 ou 3). Dans ce cas, tu as quand même intérêt à faire des regroupements même si tu perds en précision.

Maintenant, si c'est une maladie rare et que tu as peu d'individus, tu peux toujours arguer que tu as choisi un compromis entre robustesse et qualité d'ajustement et puissance de l'analyse et que tu n'avais pas assez de sujets pour détecter des interactions de type age:sexe ou age:forme et donc qu'elles peuvent être considérées comme négligeables. Tu peux alors choisir le modèle ~ age + forme

Super. Tes explications m'ont grandement aidé.
Merci.