Qual é o número mínimo recomendado de grupos para um fator de efeitos aleatórios?

Estou usando um modelo misto em R( lme4) para analisar alguns dados de medidas repetidas. Eu tenho uma variável de resposta (teor de fibras nas fezes) e 3 efeitos fixos (massa corporal, etc.). Meu estudo tem apenas 6 participantes, com 16 medidas repetidas para cada um (embora dois tenham apenas 12 repetições). Os sujeitos são lagartos que receberam diferentes combinações de alimentos em diferentes 'tratamentos'.

Minha pergunta é: posso usar o ID do assunto como um efeito aleatório?

Sei que esse é o curso de ação usual em modelos longitudinais de efeitos mistos, para levar em consideração a natureza amostrada aleatoriamente dos sujeitos e o fato de que as observações dentro dos sujeitos serão mais intimamente correlacionadas do que as entre os sujeitos. Mas, tratar o ID do sujeito como um efeito aleatório envolve estimar uma média e uma variação para essa variável.

• Como eu tenho apenas 6 sujeitos (6 níveis desse fator), isso é suficiente para obter uma caracterização precisa da média e variância?

• O fato de eu ter várias medidas repetidas para cada assunto ajuda nesse sentido (não vejo como isso importa)?

• Por fim, se eu não puder usar o ID do assunto como um efeito aleatório, incluí-lo como um efeito fixo permitirá que eu controle o fato de ter repetido medidas?

Edit: Eu gostaria de esclarecer que quando digo "posso" usar a ID do assunto como um efeito aleatório, quero dizer "é uma boa idéia". Eu sei que posso ajustar o modelo com um fator com apenas 2 níveis, mas certamente isso seria in defensável? Estou perguntando em que momento é sensato pensar em tratar assuntos como efeitos aleatórios? Parece que a literatura recomenda que 5-6 níveis sejam um limite inferior. Parece-me que as estimativas da média e variância do efeito aleatório não seriam muito precisas até que houvesse mais de 15 níveis de fatores.

Resposta curta: Sim, você pode usar o ID como efeito aleatório com 6 níveis.

Resposta um pouco mais longa: A FAQ do GLMM do @ BenBolker diz (entre outras coisas) o seguinte sob o título " Devo tratar o fator xxx como fixo ou aleatório? ":

Um ponto de particular relevância para a estimativa do modelo misto "moderno" (em vez da estimativa do método dos momentos "clássica") é que, para fins práticos, deve haver um número razoável de níveis de efeitos aleatórios (por exemplo, blocos) - mais de 5 ou 6 no mínimo.

Então você está no limite inferior, mas no lado direito.

Na tentativa de descobrir o número mínimo de grupos para um modelo multinível, observei o livro Análise de Dados Usando Regressão e Modelos Mulitilevel / Hierárquicos de Gelman e Hill (2007).

Eles parecem abordar esse tópico no Capítulo 11, Seção 5 (p 247), onde escrevem que, quando há <5 grupos, os modelos multiníveis geralmente adicionam pouco sobre os modelos clássicos. No entanto, eles parecem escrever que há pouco risco de aplicar um modelo multinível.

Os mesmos autores parecem retornar a esse tópico no capítulo 12, seção 9 (páginas 275-276). Lá eles escrevem que o conselho sobre o número mínimo de grupos para um modelo multinível é equivocado. Lá, eles dizem novamente que os modelos multiníveis geralmente acrescentam pouco aos modelos clássicos quando o número de grupos é pequeno. No entanto, eles também escrevem que os modelos multiníveis não devem fazer pior do que a regressão sem pool (onde sem pool parece significar que os indicadores de grupo são usados ​​na regressão clássica).

Nas páginas 275-276, os autores têm uma subseção específica para o caso de um ou dois grupos (por exemplo, homem versus mulher). Aqui eles escrevem que geralmente expressam o modelo na forma clássica. No entanto, eles afirmam que a modelagem multinível pode ser útil mesmo com apenas um ou dois grupos. Eles escrevem que, com um ou dois grupos, a modelagem multinível se reduz à regressão clássica.

Minha impressão disso é que a regressão clássica é um final de um continuum de modelos, isto é, um caso especial de um modelo multinível.

Com base no exposto, minha impressão é que a regressão clássica e a modelagem multinível retornarão estimativas quase idênticas quando houver apenas dois grupos, e que o uso de modelos multiníveis com apenas um, dois, três, quatro, cinco ou seis está bem.

Tentarei modificar essa resposta no futuro com Rcódigo e um pequeno conjunto de dados comparando estimativas obtidas com as duas abordagens ao usar dois grupos.

Pelo que vale a pena, fiz um pouco de um estudo de simulação para examinar a estabilidade da estimativa de variância para um LMM relativamente simples (usando o sleepstudyconjunto de dados disponível por meio lme4). A primeira maneira gera todas as combinações possíveis de assuntos para o ngroupsnúmero de assuntos e reajusta o modelo para cada combinação possível. O segundo leva vários subconjuntos aleatórios de assuntos.

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)

A linha preta pontilhada é a estimativa original do ponto da variação, e as facetas representam diferentes números de sujeitos ( s3sendo grupos de três sujeitos, s4sendo quatro, etc.).

E a maneira alternativa:

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)

Parece (para este exemplo, de qualquer maneira) que a variação não se estabiliza realmente até que haja pelo menos 14 assuntos, se não mais tarde.

A "Econometria Principalmente Inofensiva" de Angrist e Pischke tem uma seção intitulada "Menos de 42 aglomerados", na qual eles dizem sem brincadeira:

Portanto, seguindo o ditado de que a resposta para a vida, o universo e tudo é 42, acreditamos que a pergunta é: quantos clusters são suficientes para inferência confiável usando o ajuste de cluster padrão [semelhante ao estimador de variância no GEE]?

A maneira como meu instrutor de econometria costumava responder perguntas como a sua é: "Os Estados Unidos são um país livre, você pode fazer o que quiser. Mas, se você quiser publicar seu trabalho, precisará defender o que fez. " Em outras palavras, você provavelmente poderá executar o código R ou Stata ou HLM ou Mplus ou SAS PROC GLIMMIX com 6 indivíduos (e alternar para esses pacotes alternativos se o de sua escolha não executar isso), mas provavelmente terá muito difícil defender essa abordagem e justificar testes assintóticos.

Acredito que, por padrão, incluir uma variável como uma inclinação aleatória implica incluir isso também como um efeito fixo, e você precisa passar por muitos obstáculos na sintaxe se quiser ter isso como um efeito aleatório com a média de zero. Essa é uma escolha sensata que os desenvolvedores de software fizeram para você.

Você também pode usar um modelo misto bayesiano - nesse caso, a incerteza na estimativa dos efeitos aleatórios é totalmente resolvida no cálculo dos intervalos credíveis de previsão de 95%. O novo pacote brmse função R brm, por exemplo, permite uma transição muito fácil de um lme4modelo misto freqüentista para um modelo bayesiano, pois possui uma sintaxe quase idêntica.

Eu não usaria um modelo de efeitos aleatórios com apenas 6 níveis. Os modelos que usam um efeito aleatório de 6 níveis podem, em algum momento, ser executados usando muitos programas estatísticos e, às vezes, fornecem estimativas imparciais, mas:

1. Eu acho que existe um consenso arbitrário na comunidade estatística de que 10-20 é o número mínimo. Se você deseja publicar sua pesquisa, será aconselhado a procurar uma revista sem revisão estatística (ou poderá justificar sua decisão usando uma linguagem bastante sofisticada).
2. Com tão poucos clusters, é provável que a variação entre os cluster seja pouco estimada. A estimativa ruim da variação entre os agrupamentos geralmente se traduz em uma estimativa ruim do erro padrão dos coeficientes de interesse. (os modelos de efeitos aleatórios dependem do número de clusters teoricamente indo para o infinito).
3. Muitas vezes, os modelos simplesmente não convergem. Você já tentou executar seu modelo? Eu ficaria surpreso com apenas 12 a 16 medidas por assunto que os modelos convergem. Quando consegui convergir esse tipo de modelo, tive centenas de medições por cluster.

Esse problema foi solucionado na maioria dos livros-padrão do campo e você os abordou na sua pergunta. Acho que não estou lhe dando nenhuma informação nova.

Já faz muito tempo desde a pergunta original, mas pensei em acrescentar alguns pontos pertinentes à seleção de modelos.

1 - Desde que o modelo seja identificado (ou seja, você tenha graus de liberdade no espaço de parâmetros), você poderá TENTAR ajustar-se ao modelo. Dependendo do método de otimização, o modelo pode ou não convergir. De qualquer forma, eu não tentaria incluir mais de 1 ou 2 efeitos aleatórios e, definitivamente, não mais que 1 interação de nível cruzado. No caso específico do problema apresentado aqui, se suspeitarmos de uma interação entre características específicas do lagarto (por exemplo, idade, tamanho, etc.) e características do tratamento / medida, o tamanho do grupo 6 pode não ser suficiente para fazer estimativas suficientemente precisas.

2 - Como mencionam algumas respostas, a convergência pode ser um problema. No entanto, minha experiência é que, embora os dados das ciências sociais tenham um enorme problema de convergência devido a problemas de medição, as ciências da vida e especialmente as medidas repetidas bioquímicas têm erros padrão muito menores. Tudo depende do processo de geração de dados. Em dados sociais e econômicos, temos que trabalhar em vários níveis de abstração. Em dados biológicos e químicos e, certamente, astronômicos, o erro de medição é um problema menor.