Como posso testar se um efeito aleatório é significativo?

34

Estou tentando entender quando usar um efeito aleatório e quando é desnecessário. Foi-me dito uma regra prática: se você tem 4 ou mais grupos / indivíduos, eu (15 alces individuais). Alguns desses alces foram experimentados 2 ou 3 vezes, para um total de 29 tentativas. Quero saber se eles se comportam de maneira diferente quando estão em paisagens de maior risco do que não. Então, pensei em definir o indivíduo como um efeito aleatório. No entanto, agora me disseram que não há necessidade de incluir o indivíduo como um efeito aleatório, porque não há muita variação em sua resposta. O que não consigo descobrir é como testar se realmente há algo sendo levado em consideração ao definir o indivíduo como um efeito aleatório. Talvez uma pergunta inicial seja: Que teste / diagnóstico posso fazer para descobrir se o indivíduo é uma boa variável explicativa e deve ser um efeito fixo - gráficos qq? histogramas? gráficos de dispersão? E o que eu procuraria nesses padrões.

Eu executei o modelo com o indivíduo como um efeito aleatório e sem, mas depois li http://glmm.wikidot.com/faq onde eles afirmam:

não compare modelos lmer com os ajustes lm correspondentes ou glmer / glm; as probabilidades de log não são proporcionais (ou seja, incluem diferentes termos aditivos)

E aqui presumo que isso significa que você não pode comparar entre um modelo com efeito aleatório ou sem. Mas eu realmente não saberia o que devo comparar entre eles de qualquer maneira.

No meu modelo com o efeito Aleatório, eu também estava tentando analisar o resultado para ver que tipo de evidência ou significado o ER tem

lmer(Velocity ~ D.CPC.min + FD.CPC + (1|ID), REML = FALSE, family = gaussian, data = tv)

Linear mixed model fit by maximum likelihood 
Formula: Velocity ~ D.CPC.min + FD.CPC + (1 | ID) 
   Data: tv 
    AIC    BIC logLik deviance REMLdev
 -13.92 -7.087  11.96   -23.92   15.39
Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.00000  0.00000 
 Residual             0.02566  0.16019 
Number of obs: 29, groups: ID, 15

Fixed effects:
              Estimate Std. Error t value
(Intercept)  3.287e-01  5.070e-02   6.483
D.CPC.min   -1.539e-03  3.546e-04  -4.341
FD.CPC       1.153e-04  1.789e-05   6.446

Correlation of Fixed Effects:
          (Intr) D.CPC.
D.CPC.min -0.010       
FD.CPC    -0.724 -0.437

Você vê que minha variação e DP do ID individual como efeito aleatório = 0. Como isso é possível? O que significa 0? Isso esta certo? Então, meu amigo que disse "como não há variação usando o ID como efeito aleatório é desnecessário" está correto? Então, então eu o usaria como um efeito fixo? Mas o fato de haver tão pouca variação significa que isso não vai nos dizer muito?

Kerry
fonte
Para obter uma variação exata de 0 de um efeito aleatório, consulte stats.stackexchange.com/questions/115090 .
ameba diz Restabelecer Monica

Respostas:

21

A estimativa, IDvariância de = 0, indica que o nível de variabilidade entre os grupos não é suficiente para garantir a incorporação de efeitos aleatórios no modelo; ie seu modelo é degenerado.

Como você se identifica corretamente: provavelmente, sim; IDcomo efeito aleatório é desnecessário. Poucas coisas vêm à mente para testar esta suposição:

  1. Você pode comparar (usando REML = Fsempre) o AIC (ou seu IC favorito em geral) entre um modelo com e sem efeitos aleatórios e ver como isso acontece.
  2. Você examinaria a anova()saída dos dois modelos.
  3. Você pode fazer uma inicialização paramétrica usando a distribuição posterior definida pelo seu modelo original.

Lembre-se de que as opções 1 e 2 têm um problema: você está procurando algo que esteja nos limites do espaço dos parâmetros para que, na verdade, eles não sejam tecnicamente sólidos. Dito isto, acho que você não terá idéias erradas deles e muitas pessoas os usam (por exemplo, Douglas Bates, um dos desenvolvedores do lme4, os usa em seu livro, mas afirma claramente essa ressalva sobre os valores dos parâmetros que estão sendo testados no limite do conjunto de valores possíveis). A escolha 3 é a mais entediante das 3, mas na verdade dá a melhor idéia sobre o que está acontecendo. Algumas pessoas também são tentadas a usar o bootstrap não paramétrico, mas acho que, considerando que você está fazendo suposições paramétricas para começar, também pode usá-las.

usεr11852 diz Reinstate Monic
fonte
6
O pacote RLRsim é uma maneira realmente conveniente de testar efeitos aleatórios usando testes de razão de verossimilhança com base em simulação.
Atrichornis
@atrichornis: +1. Pacote interessante; Eu não sabia disso. Acabei de dar uma olhada no código, bem direto, posso dizer. Eu gostaria que eles o incorporassem (ou algo parecido) lme4especialmente agora que mcmcsamp()está quebrado e as pessoas ficam com apenas suas próprias implementações de inicialização ad-hoc para obter alguns valores p decentes etc.
usεr11852 diz Reinstate Monic
É verdade que os modelos mistos não são diretos em R. Muitas aproximações e soluções alternativas ... Embora eu colete SAS etc., apenas encobrir algumas das mesmas incertezas? Ben Bolker é co-autor de ambos os pacotes, ele pode ter suas razões para não incorporá-lo. Provavelmente hora!
Atrichornis
4
A inicialização no limite do espaço de parâmetro tem seu próprio conjunto de questões e problemas que levam à inconsistência . O bootstrap não é uma panacéia e não deve ser jogado na bolsa levemente, assumindo que resolverá tudo.
StasK
2
Dê uma olhada, o argumento é muito sutil. Tanto quanto me lembro, tudo se resume ao fato de você estar fazendo o bootstrap a partir de uma distribuição diferente da nula; e dadas as distribuições não padronizadas obtidas no limite, as condições de regularidade são violadas e a distribuição de autoinicialização não converge para o destino. Eu acho que o bootstrap não paramétrico ainda pode ser construído aqui, retirando-se as médias de resíduos do grupo. No entanto, com a violação da independência das observações entre os grupos, outra camada de complicações pode surgir.
StasK
3

Não tenho certeza de que a abordagem que vou sugerir seja razoável; portanto, aqueles que sabem mais sobre esse tópico me corrigem se eu estiver errado.

Minha proposta é criar uma coluna adicional em seus dados que tenha um valor constante de 1:

IDconst <- factor(rep(1, each = length(tv$Velocity)))

Em seguida, você pode criar um modelo que use essa coluna como seu efeito aleatório:

fm1 <- lmer(Velocity ~ D.CPC.min + FD.CPC + (1|IDconst), 
  REML = FALSE, family = gaussian, data = tv)

Nesse ponto, você pode comparar (AIC) seu modelo original com o efeito aleatório ID(vamos chamá-lo fm0) com o novo modelo que não leva em consideração, IDpois IDconsté o mesmo para todos os seus dados.

anova(fm0,fm1)

Atualizar

O usuário11852 estava pedindo um exemplo, porque, na sua opinião, a abordagem acima nem será executada. Pelo contrário, posso mostrar que a abordagem funciona (pelo menos com o lme4_0.999999-0que estou usando atualmente).

set.seed(101)
dataset <- expand.grid(id = factor(seq_len(10)), fac1 = factor(c("A", "B"),
  levels = c("A", "B")), trial = seq_len(10))
dataset$value <- rnorm(nrow(dataset), sd = 0.5) +
      with(dataset, rnorm(length(levels(id)), sd = 0.5)[id] +
      ifelse(fac1 == "B", 1.0, 0)) + rnorm(1,.5)
    dataset$idconst <- factor(rep(1, each = length(dataset$value)))

library(lme4)
fm0 <- lmer(value~fac1+(1|id), data = dataset)
fm1 <- lmer(value~fac1+(1|idconst), data = dataset)

anova(fm1,fm0)

Saída:

  Data: dataset
  Models:
  fm1: value ~ fac1 + (1 | idconst)
  fm0: value ~ fac1 + (1 | id)

      Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
  fm1  4 370.72 383.92 -181.36                      
  fm0  4 309.79 322.98 -150.89 60.936      0  < 2.2e-16 ***

De acordo com este último teste, devemos manter o efeito aleatório, uma vez que o fm0modelo tem a AIC mais baixa e a BIC.

Atualização 2

A propósito, essa mesma abordagem é proposta por NW Galwey em 'Introdução à modelagem mista: além da regressão e análise de variância' nas páginas 213-214.

VLC
fonte
Você já testou sua ideia? Por favor, prove que estou errado, mas acho que sua ideia nem será executada. Se o IDconstmesmo for para todos os seus dados, você não terá nenhum agrupamento. Você precisa de um fator de agrupamento para ter pelo menos um nível amostrado e a maneira como você configura o modelo não possui nenhum. Talvez eu acredite na lógica de usar um "agrupamento aleatório", mas esse é um jogo diferente. Teste sua abordagem com alguns dados fictícios. Eu acredito firmemente que a configuração proposta lmer()não será executada. ( lme4_0.99999911-1
Usava
@ user11852 Por favor, veja minha atualização e deixe-nos saber se essa abordagem também funciona lme4_0.99999911-1.
VLC
Z
3
Sim, eu fiz o que você sugere; não vai funcionar / computar. Error in lFormula(formula = value ~ fac1 + (1 | idconst), data = dataset) : grouping factors must have at least 1 sampled level. E como eu disse, conceitualmente está errado. Não é uma questão de enganar o software para fornecer alguns números, é uma questão de se o que você diz é razoável. Você não tem um segundo modelo misto para comparar se, nesse modelo, o efeito aleatório é por construção uma constante. Você também pode excluí-lo e tentar um modelo linear.
usεr11852 diz Reinstate Monic
11
Atualização de concerto definindo uma variável aleatória de grupo único em lme4. Isso pode ser feito se você definir a opção: control=lmerControl(check.nlev.gtr.1="ignore"). Ben Bolker menciona aqui: github.com/lme4/lme4/issues/411 .
Robin Beaumont
1

Eu gostaria de responder à pergunta mais 'inicial'.

Se você suspeitar de algum tipo de heterogeneidade na variação entre a variável dependente devido a alguns fatores, prossiga e plote os dados usando gráficos de dispersão e caixa. Alguns padrões comuns para verificar, eu coloquei esta lista abaixo de várias fontes na web.

Padrões de heterocedasticidade

Além disso, plote sua variável dependente por fator / grupos de tratamento para ver se há variação constante. Caso contrário, você pode explorar efeitos aleatórios ou regressões ponderadas. Por exemplo. Este gráfico abaixo é um exemplo de variação em forma de funil em meus grupos de tratamento. Então, eu escolho fazer efeitos aleatórios e testar os efeitos na inclinação e nas interceptações.

Boxplot para verificar a heterocedasticidade

A partir daqui, as respostas acima atendem à sua pergunta principal. Também existem testes que verificam a heterocedasticidade, um deles está aqui - https://dergipark.org.tr/download/article-file/94971 . Mas não tenho certeza se existem testes para detectar heterocedasticidade no nível do grupo.

Uma corrida
fonte
Por favor, use apenas o campo "Sua resposta" para fornecer respostas à pergunta do OP. O CV é um site estrito de perguntas e respostas, não um fórum de discussão. A última parte em negrito da sua postagem é uma nova pergunta, não uma resposta a essa pergunta. Se você tiver uma nova pergunta, clique no cinza ASK QUESTIONna parte superior e faça a pergunta lá. Como você é novo aqui, convém fazer um tour , que contém informações para novos usuários.
gung - Restabelece Monica