Porównanie modelu mieszanego (podmiot jako efekt losowy) z prostym modelem liniowym (przedmiot jako efekt stały)

Kończę analizę dużej grupy danych. Chciałbym wziąć model liniowy zastosowany w pierwszej części pracy i ponownie go dopasować za pomocą liniowego modelu mieszanego (LME). LME byłby bardzo podobny, z tym wyjątkiem, że jedna ze zmiennych zastosowanych w modelu byłaby zastosowana jako efekt losowy. Te dane pochodzą z wielu obserwacji (> 1000) w małej grupie podmiotów (~ 10) i wiem, że modelowanie efektu podmiotu lepiej jest wykonać jako efekt losowy (jest to zmienna, którą chcę przesunąć). Kod R wyglądałby następująco:

my_modelB <- lm(formula = A ~ B + C + D)    
lme_model <- lme(fixed=A ~ B + C, random=~1|D, data=my_data, method='REML')

Wszystko działa dobrze, a wyniki są bardzo podobne. Byłoby miło, gdybym mógł użyć czegoś takiego jak RLRsim lub AIC / BIC, aby porównać te dwa modele i zdecydować, który jest najbardziej odpowiedni. Moi koledzy nie chcą zgłaszać LME, ponieważ nie ma łatwo dostępnego sposobu wyboru, który jest „lepszy”, mimo że uważam, że LME jest bardziej odpowiednim modelem. Jakieś sugestie?

Ma to na celu dodanie odpowiedzi @ ocram, ponieważ jest zbyt długi, aby opublikować komentarz. Traktowałbym A ~ B + Cjak twój model zerowy, abyś mógł ocenić istotność statystyczną Dprzechwytywania losowego na poziomie w konfiguracji modelu zagnieżdżonego. Jak wskazał ocram, warunki regularności są naruszane, gdy , a statystyka testu wskaźnika prawdopodobieństwa (LRT) niekoniecznie będzie asymptotycznie rozłożona χ 2 . Rozwiązaniem, którego nauczyłem, było ładowanie LRT (którego rozkład prawdopodobnie nie będzie wynosił χ 2 ) parametrycznie i obliczanie wartości p ładowania w następujący sposób:$H_0: \sigma^2 = 0$$\chi^2$$\chi^2$

library(lme4)
my_modelB <- lm(formula = A ~ B + C)
lme_model <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
lrt.observed <- as.numeric(2*(logLik(lme_model) - logLik(my_modelB)))
nsim <- 999
lrt.sim <- numeric(nsim)
for (i in 1:nsim) {
    y <- unlist(simulate(mymodlB))
    nullmod <- lm(y ~ B + C)
    altmod <- lmer(y ~ B + C + (1|D), data=my_data, REML=F)
    lrt.sim[i] <- as.numeric(2*(logLik(altmod) - logLik(nullmod)))
}
mean(lrt.sim > lrt.observed) #pvalue

Odsetek LRT z ładowaniem początkowym bardziej ekstremalny niż obserwowany LRT jest wartością p.

$0$$H_0:variance = 0$$H_1: variance > 0$...

EDYTOWAĆ

Aby uniknąć nieporozumień: Wspomniany powyżej test jest czasem wykorzystywany do podjęcia decyzji, czy efekt losowy jest znaczący ... ale nie do podjęcia decyzji, czy należy go przekształcić w efekt stały.