Eu tenho um conjunto de dados de medidas repetidas desequilibrado para analisar e li que a maioria dos pacotes estatísticos lida com isso com ANOVA (ou seja, soma dos quadrados do tipo III) está errada. Portanto, eu gostaria de usar um modelo de efeitos mistos para analisar esses dados. Eu li muito sobre modelos mistos R, mas ainda sou muito novo Re modelos de efeito misto e não estou muito confiante de que estou fazendo as coisas corretamente. Observe que ainda não posso me divorciar completamente dos métodos "tradicionais" e ainda preciso de valores- $p$ e testes post hoc.

Gostaria de saber se a abordagem a seguir faz sentido ou se estou fazendo algo terrivelmente errado. Aqui está o meu código:

# load packages
library(lme4)
library(languageR)
library(LMERConvenienceFunctions)
library(coda)
library(pbkrtest)

# import data
my.data <- read.csv("data.csv")

# create separate data frames for each DV & remove NAs
region.data <- na.omit(data.frame(time=my.data$time, subject=my.data$subject, dv=my.data$dv1))

# output summary of data
data.summary <- summary(region.data)

# fit model
# "time" is a factor with three levels ("t1", "t2", "t3")
region.lmer <- lmer(dv ~ time + (1|subject), data=region.data)

# check model assumptions
mcp.fnc(region.lmer)

# remove outliers (over 2.5 standard deviations)
rm.outliers <- romr.fnc(region.lmer, region.data, trim=2.5)
region.data <- rm.outliers$data
region.lmer <- update(region.lmer)

# re-check model assumptions
mcp.fnc(region.lmer)

# compare model to null model
region.lmer.null <- lmer(dv ~ 1 + (1|subject), data=region.data)
region.krtest <- KRmodcomp(region.lmer, region.lmer.null)

# output lmer summary
region.lmer.summary <- summary(region.lmer)

# run post hoc tests
t1.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t2") + (1|subject), data=region.data)
t2.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

region.lmer <- lmer(dv ~ relevel(time,ref="t3") + (1|subject), data=region.data)
t3.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE)

# Get mcmc mean and 50/95% HPD confidence intervals for graphs
# repeated three times and stored in a matrix (not shown here for brevity)
as.numeric(t1.pvals$fixed$MCMCmean)
as.numeric(t1.pvals$fixed$HPD95lower)
as.numeric(t1.pvals$fixed$HPD95upper)
HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)
    HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)

Algumas perguntas específicas que tenho:

1. Essa é uma maneira válida de analisar modelos de efeitos mistos? Caso contrário, o que devo fazer em seu lugar.
2. Os gráficos de crítica gerados pelo mcp.fnc são bons o suficiente para verificar as suposições do modelo ou devo tomar medidas adicionais.
3. Entendo que, para modelos mistos serem válidos, os dados precisam respeitar suposições de normalidade e homoscedasticidade. Como julgar o que é "aproximadamente normal" e o que não é olhando os gráficos de críticas gerados pelo mcp.fnc? Eu só preciso ter uma idéia disso, ou é uma maneira prescrita de fazer as coisas? Quão robustos são os modelos mistos em relação a essas suposições?
4. Preciso avaliar as diferenças entre os três momentos para ~ 20 características (biomarcadores) dos sujeitos da minha amostra. A montagem e o teste de modelos separados para cada um são aceitáveis, desde que eu relate todos os testes realizados (significativos ou não) ou precise de alguma forma de correção para comparações múltiplas.

Para ser um pouco mais preciso em relação ao experimento, aqui estão mais alguns detalhes. Seguimos um número de participantes longitudinalmente, enquanto eles eram submetidos a um tratamento. Medimos vários biomarcadores antes do início do tratamento e em dois momentos. O que eu gostaria de ver é se há diferença nesses biomarcadores entre os três momentos.

Estou baseando a maior parte do que estou fazendo aqui neste tutorial , mas fiz algumas alterações com base nas minhas necessidades e nas coisas que li. As alterações que eu fiz são:

1. relevel o fator "time" para obter comparações t1-t2, t2-t3 e t1-t3 com pvals.fnc (do pacote languageR)
2. compare meu modelo misto com o modelo nulo usando um teste F aproximado com base na abordagem de Kenward-Roger (usando o pacote pbkrtest) em vez de um teste de razão de verossimilhança (porque eu li que Kenward-Roger é mais bem visto no momento)
3. Use o pacote LMERConvenienceFunctions para verificar suposições e remover discrepâncias (porque li que modelos mistos são muito sensíveis a discrepâncias)

Suas perguntas são um pouco "grandes", então começarei com alguns comentários e dicas gerais.

Algumas leituras em segundo plano e pacotes úteis

Você provavelmente deve dar uma olhada em algumas das introduções do tutorial sobre o uso de modelos mistos, eu recomendaria começar com Baayen et al (2008) - Baayen é o autor de languageR. Barr et al (2013) discutem alguns problemas com a estrutura de efeitos aleatórios, e Ben Bolker é um dos desenvolvedores atuais do lme4. Baayen et al é especialmente bom para suas perguntas, porque elas gastam muito tempo discutindo o uso de testes de inicialização / permutação (o que está por trás mcp.fnc).

Literatura

• Baayen, RH, D. Davidson e D. Bates (2008). Modelagem de efeitos mistos com efeitos aleatórios cruzados para assuntos e itens . Journal of Memory and Language, 59 (390-412).
• Barr, Dale J, R. Levy, C. Scheepers e HJ Tily (2013). Estrutura de efeitos aleatórios para teste de hipótese confirmatória: Mantenha-o no máximo . Journal of Memory and Language, 68: 255– 278.
• Bolker, Benjamin M, ME Brooks, CJ Clark, SW Geange, JR Poulsen, MHH Stevens e J.-SS White (2009). Modelos mistos lineares generalizados: um guia prático para ecologia e evolução . Trends Ecol Evol, 24 (3): 127–35.

Florian Jaeger também tem um monte de coisas sobre o lado prático de modelos mistos no blog de seu laboratório .

Há também vários pacotes R úteis para visualizar e testar modelos mistos, como lmerTeste effects. O effectspacote é especialmente bom porque permite plotar os intervalos de regressão linear e de confiança acontecendo nos bastidores.

Qualidade do ajuste e comparação de modelos

$p$lmerTest

anova()mer$\chi^2$$\chi^2$$p$-valores para comparar diretamente dois modelos. A desvantagem disso é que não está claro imediatamente o quão bom é o seu ajuste.

$t$summary()$|t|>2$fixef()

Você também deve garantir que nenhum dos seus efeitos fixos esteja fortemente correlacionado - a multicolinearidade viola as suposições do modelo. Florian Jaeger escreveu um pouco sobre isso, bem como algumas soluções possíveis.

Comparações múltiplas

Abordarei sua pergunta # 4 um pouco mais diretamente. A resposta é basicamente a mesma que as boas práticas da ANOVA tradicional, infelizmente este parece ser um ponto em que há muita incerteza para a maioria dos pesquisadores. (É o mesmo que a ANOVA tradicional porque os modelos lineares de efeitos mistos e ANOVA são baseados no modelo linear geral, é apenas que os modelos mistos têm um termo extra para os efeitos aleatórios.) Se você está assumindo que as janelas de tempo produzem um efeito diferença e deseja comparar os efeitos do tempo, inclua o tempo no seu modelo. Aliás, isso também fornecerá uma maneira conveniente para você julgar se o tempo fez alguma diferença: existe um efeito principal (fixo) para o tempo? A desvantagem de seguir esse caminho é que seu modelo ficará MUITO mais complexo e o único "super" um modelo com tempo como paramater provavelmente levará mais tempo para computar do que três modelos menores sem tempo como paramater. De fato, o exemplo clássico de tutorial para modelos mistos ésleepstudy que usa o tempo como um parâmetro.

$t$foreachlme4$\chi^2$

Se suas características forem a variável dependente, será necessário calcular modelos diferentes de qualquer maneira e, em seguida, você poderá usar o AIC e o BIC para comparar os resultados.