Quel est le nombre minimum recommandé de groupes pour un facteur d'effets aléatoires?

26

J'utilise un modèle mixte dans R( lme4) pour analyser certaines données de mesures répétées. J'ai une variable de réponse (teneur en fibres des matières fécales) et 3 effets fixes (masse corporelle, etc.). Mon étude ne compte que 6 participants, avec 16 mesures répétées pour chacun (même si deux n'ont que 12 répétitions). Les sujets sont des lézards qui ont reçu différentes combinaisons d'aliments dans différents «traitements».

Ma question est: puis-je utiliser l'ID de sujet comme un effet aléatoire?

Je sais que c'est la ligne de conduite habituelle dans les modèles longitudinaux à effets mixtes, pour tenir compte de la nature échantillonnée au hasard des sujets et du fait que les observations au sein des sujets seront plus étroitement corrélées que celles entre les sujets. Mais, traiter l'ID du sujet comme un effet aléatoire implique d'estimer une moyenne et une variance pour cette variable.

  • Comme je n'ai que 6 sujets (6 niveaux de ce facteur), est-ce suffisant pour obtenir une caractérisation précise de la moyenne et de la variance?

  • Le fait d'avoir plusieurs mesures répétées pour chaque sujet aide-t-il à cet égard (je ne vois pas en quoi cela est important)?

  • Enfin, si je ne peux pas utiliser l'ID de sujet comme effet aléatoire, l'inclure comme effet fixe me permettra-t-il de contrôler le fait que j'ai des mesures répétées?

Edit: Je voudrais juste préciser que lorsque je dis "puis-je" utiliser l'ID du sujet comme un effet aléatoire, je veux dire "est-ce une bonne idée de". Je sais que je peux adapter le modèle avec un facteur avec seulement 2 niveaux, mais ce serait sûrement in défendable? Je demande à quel moment devient-il raisonnable de penser à traiter les sujets comme des effets aléatoires? Il semble que la littérature indique que 5-6 niveaux est une limite inférieure. Il me semble que les estimations de la moyenne et de la variance de l'effet aléatoire ne seraient pas très précises tant qu'il n'y aurait pas plus de 15 facteurs.

Chris
la source

Réponses:

21

Réponse courte: Oui, vous pouvez utiliser ID comme effet aléatoire avec 6 niveaux.

Réponse légèrement plus longue: La FAQ GLMM de @ BenBolker dit (entre autres) ce qui suit sous le titre " Dois-je traiter le facteur xxx comme fixe ou aléatoire? ":

Un point particulièrement pertinent pour l’estimation de modèle mixte «moderne» (plutôt que l’estimation de la méthode des moments «classique») est que, pour des raisons pratiques, il doit y avoir un nombre raisonnable de niveaux d’effets aléatoires (par exemple, des blocs) - plus de 5 ou 6 au minimum.

Vous êtes donc à la limite inférieure, mais à droite.

Henrik
la source
12

Pour tenter de déterminer le nombre minimum de groupes pour un modèle à plusieurs niveaux, j'ai examiné le livre Data Analysis Using Regression and Mulitilevel / Hierarchical models de Gelman et Hill (2007).

Ils semblent aborder ce sujet dans le chapitre 11, section 5 (page 247) où ils écrivent que lorsqu'il y a moins de 5 groupes, les modèles à plusieurs niveaux ajoutent généralement peu par rapport aux modèles classiques. Cependant, ils semblent écrire qu'il y a peu de risques à appliquer un modèle à plusieurs niveaux.

Les mêmes auteurs semblent revenir sur ce sujet dans le chapitre 12, section 9 (pages 275-276). Là, ils écrivent que des conseils sur le nombre minimum de groupes pour un modèle à plusieurs niveaux sont erronés. Là encore, ils disent que les modèles à plusieurs niveaux ajoutent souvent peu par rapport aux modèles classiques lorsque le nombre de groupes est petit. Cependant, ils écrivent également que les modèles multiniveaux ne devraient pas faire pire que la régression sans regroupement (où le non-regroupement semble signifier que les indicateurs de groupe sont utilisés dans la régression classique).

Aux pages 275-276, les auteurs ont une sous-section spécifique pour le cas d'un ou deux groupes (par exemple, hommes contre femmes). Ici, ils écrivent qu'ils expriment généralement le modèle sous une forme classique. Cependant, ils indiquent que la modélisation à plusieurs niveaux peut être utile même avec seulement un ou deux groupes. Ils écrivent qu'avec un ou deux groupes, la modélisation à plusieurs niveaux se réduit à une régression classique.

Mon impression est que la régression classique est une extrémité d'un continuum de modèles, c'est-à-dire un cas particulier d'un modèle à plusieurs niveaux.

Sur la base de ce qui précède, j'ai l'impression que la régression classique et la modélisation multiniveaux renverront des estimations presque identiques lorsqu'il n'y a que deux groupes et que l'utilisation de modèles multiniveaux avec seulement un, deux, trois, quatre, cinq ou six groupes est acceptable.

J'essaierai de modifier cette réponse à l'avenir avec du Rcode et un petit ensemble de données comparant les estimations obtenues avec les deux approches lors de l'utilisation de deux groupes.

Mark Miller
la source
10

Pour ce que ça vaut, j'ai fait une petite étude de simulation pour examiner la stabilité de l'estimation de la variance pour un LMM relativement simple (en utilisant l' sleepstudyensemble de données disponible via lme4). La première méthode génère toutes les combinaisons ngroupsde sujets possibles pour le nombre de sujets et réajuste le modèle pour chaque combinaison possible. Le second prend plusieurs sous-ensembles aléatoires de sujets.

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18 
# generate all possible combinations
combos <- lapply(X = ngroups, 
                 FUN = function(x) combn(unique(sleepstudy$Subject), x)) 

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

La ligne noire en pointillés est l'estimation ponctuelle d'origine de la variance, et les facettes représentent différents nombres de sujets ( s3étant des groupes de trois sujets, s4quatre, etc.). entrez la description de l'image ici

Et la manière alternative:

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

entrez la description de l'image ici

Il apparaît (pour cet exemple, de toute façon) que la variance ne se stabilise pas vraiment jusqu'à ce qu'il y ait au moins 14 sujets, sinon plus tard.

alexforrence
la source
1
+1. Bien sûr, plus le nombre de sujets est petit, plus la variance de l'estimateur de variance est grande. Mais je ne pense pas que ce soit ce qui importe ici. La question est, quel nombre de sujets permet d'obtenir des résultats sensibles? Si nous définissons un résultat "non sensible" comme l'obtention d'une variance nulle, alors dans votre simulation, cela se produit assez souvent avec n = 5 ou moins. En partant de n = 6 ou n = 7, vous n'obtenez presque jamais d'estimation exacte de la variance, c'est-à-dire que le modèle converge vers une solution non dégénérée. Ma conclusion serait que n = 6 est limite acceptable.
amibe dit Réintégrer Monica
1
BTW c'est en ligne avec rpubs.com/bbolker/4187 .
amibe dit Réintégrer Monica
8

Angrist et Pischke "Mostly Harmless Econometrics" a une section intitulée "Moins de 42 grappes", dans laquelle ils disent en plaisantant à moitié,

Par conséquent, suivant le ... dicton selon lequel la réponse à la vie, à l'univers et à tout est 42, nous pensons que la question est: combien de grappes suffisent pour une inférence fiable en utilisant l'ajustement de grappe standard [semblable à l'estimateur de variance dans GEE]?

La façon dont mon professeur d'économétrie répondait à des questions comme la vôtre est la suivante: «L'Amérique est un pays libre, vous pouvez faire ce que vous voulez. Mais si vous voulez que votre article soit publié, vous devez être en mesure de défendre ce que vous avez fait. " En d'autres termes, vous pourrez probablement exécuter le code R ou Stata ou HLM ou Mplus ou SAS PROC GLIMMIX avec 6 sujets (et basculer vers ces packages alternatifs si celui de votre choix ne le fait pas), mais vous aurez probablement moment très difficile pour défendre cette approche et justifier des tests asymptotiques.

Je crois que par défaut, inclure une variable comme une pente aléatoire implique également l'inclure comme un effet fixe, et vous devez sauter à travers beaucoup de cercles de syntaxe si vous voulez seulement avoir cela comme un effet aléatoire avec la moyenne de zéro. C'est un choix judicieux que les développeurs de logiciels ont fait pour vous.

StasK
la source
1
Je suppose que la réponse à la question est, dans une certaine mesure, "combien de temps dure un morceau de ficelle". Mais, je ne ferais pas beaucoup confiance à l'estimation d'une moyenne ou d'une variance à partir d'un échantillon inférieur à 15-20, donc la même règle empirique ne s'appliquerait-elle pas aux niveaux d'un effet aléatoire. Je n'ai jamais vu personne inclure l'ID de sujet comme effet fixe et aléatoire dans les études longitudinales - est-ce une pratique courante?
Chris
En plus d'un petit nombre de sujets dans un modèle mixte, leurs effets aléatoires ne sont pas observés, vous devez donc les taquiner hors des données, et sans doute vous avez besoin de relativement plus de données pour le faire de manière fiable que pour simplement estimer la moyenne et la variance quand tout est observé. Ainsi 42 vs 15-20 :). Je pense que je voulais dire des pentes aléatoires, comme vous avez raison dans les ID de sujet qui ne représentent que des effets aléatoires, sinon ils ne seront pas identifiés. Les économistes ne croient d'ailleurs pas aux effets aléatoires et publient presque exclusivement ce qu'ils appellent des «effets fixes», c'est-à-dire des estimations intra-sujet.
StasK
2
+1 @StasK pour une très bonne réponse à une question très difficile à traiter. Je pense cependant qu'il y a une teinte de sarcasme inutile et vous pourriez envisager de modifier votre réponse afin d'être un peu plus respectueux du PO.
Michael R. Chernick
@Michael, vous avez probablement raison de dire que c'est une réponse de mauvaise humeur, et probablement inutilement. Le PO a accepté la réponse qu'il voulait entendre, alors il a obtenu une résolution à ce sujet. Une réponse plus sérieuse indiquerait soit une bonne preuve de simulation, soit une analyse asymptotique d'ordre supérieur, mais malheureusement je n'ai pas connaissance de telles références.
StasK
3
Pour ce que ça vaut, je pense que le nombre magique "42" ne concerne pas le moment où les effets aléatoires sont justifiés, mais quand on peut s'en tirer sans se soucier des corrections de taille finie (par exemple en pensant aux degrés de liberté du dénominateur effectif / corrections Kenward-Roger / d'autres approches similaires).
Ben Bolker
7

Vous pouvez également utiliser un modèle mixte bayésien - dans ce cas, l'incertitude dans l'estimation des effets aléatoires est entièrement prise en compte dans le calcul des intervalles crédibles de prédiction à 95%. Le nouveau package brmset la nouvelle fonction R brm, par exemple, permettent une transition très facile d'un lme4modèle mixte fréquentiste à un modèle bayésien, car sa syntaxe est presque identique.

Tom Wenseleers
la source
4

Je n'utiliserais pas un modèle à effets aléatoires avec seulement 6 niveaux. Les modèles utilisant un effet aléatoire à 6 niveaux peuvent parfois être exécutés à l'aide de nombreux programmes statistiques et donnent parfois des estimations non biaisées, mais:

  1. Je pense qu'il y a un consensus arbitraire dans la communauté statistique que 10-20 est le nombre minimum. Si vous souhaitez publier votre recherche, il vous sera conseillé de rechercher une revue sans revue statistique (ou de pouvoir justifier votre décision en utilisant un langage assez sophistiqué).
  2. Avec si peu de grappes, la variance entre grappes est susceptible d'être mal estimée. Une mauvaise estimation de la variance entre grappes se traduit généralement par une mauvaise estimation de l'erreur-type des coefficients d'intérêt. (les modèles à effets aléatoires reposent sur le nombre de clusters allant théoriquement à l'infini).
  3. Souvent, les modèles ne convergent tout simplement pas. Avez-vous essayé d'exécuter votre modèle? Je serais surpris avec seulement 12-16 mesures par sujet que les modèles convergent. Quand j'ai réussi à faire converger ce type de modèle, j'ai eu des centaines de mesures par cluster.

Ce problème est abordé dans la plupart des manuels standard dans le domaine et vous les avez en quelque sorte abordés dans votre question. Je ne pense pas que je vous donne de nouvelles informations.

Charles
la source
Ce vote a-t-il été voté pour une raison liée à son contenu technique?
N Brouwer
Avec quel type de données travaillez-vous? Je ne sais pas pourquoi vous êtes surpris d'entendre que le modèle convergera avec 12 à 16 mesures par individu. Je ne peux pas commenter le biais dans les modèles résultants, mais je n'ai jamais eu de problèmes de convergence dans lme4les modèles mixtes et je les exécute souvent sur des tailles d'échantillons similaires à l'OP (je travaille également avec des ensembles de données de biologie).
RTbecard
1

Cela fait longtemps que la question d'origine n'a pas été posée, mais j'ai pensé ajouter quelques points pertinents à la sélection des modèles.

1 - Tant que le modèle est identifié (c'est-à-dire que vous avez des degrés de liberté dans l'espace des paramètres), vous devriez pouvoir ESSAYER pour s'adapter au modèle. Selon la méthode d'optimisation, le modèle peut converger ou non. En tout cas, je n'essaierais pas d'inclure plus de 1 ou 2 effets aléatoires et certainement pas plus d'une interaction croisée. Dans le cas spécifique du problème présenté ici, si nous soupçonnons une interaction entre les caractéristiques spécifiques du lézard (par exemple, l'âge, la taille, etc.) et les caractéristiques du traitement / mesure, la taille du groupe 6 peut ne pas être suffisante pour faire des estimations suffisamment précises.

2 - Comme quelques réponses le mentionnent, la convergence peut être un problème. Cependant, mon expérience est que, alors que les données des sciences sociales ont un énorme problème de convergence en raison de problèmes de mesure, les sciences de la vie et en particulier les mesures répétées biochimiques ont des erreurs standard beaucoup plus petites. Tout dépend du processus de génération des données. Dans les données sociales et économiques, nous devons travailler à différents niveaux d'abstraction. Dans les erreurs de mesure des données biologiques et chimiques et très certainement astronomiques, c'est moins un problème.

m_e_s
la source