Jak wiarygodne są przedziały ufności dla mniejszych obiektów dzięki pakietowi efektów?

EffectsPakiet zapewnia bardzo szybki i wygodny sposób kreślenia wyników liniowego modelu efektu mieszanego uzyskanego przez lme4pakiet . Te effectprzedziały ufności oblicza funkcyjne (CIS) bardzo szybko, ale jak wiarygodne są te przedziały ufności?

Na przykład:

library(lme4)
library(effects)
library(ggplot)

data(Pastes)

fm1  <- lmer(strength ~ batch + (1 | cask), Pastes)
effs <- as.data.frame(effect(c("batch"), fm1))
ggplot(effs, aes(x = batch, y = fit, ymin = lower, ymax = upper)) + 
  geom_rect(xmax = Inf, xmin = -Inf, ymin = effs[effs$batch == "A", "lower"],
        ymax = effs[effs$batch == "A", "upper"], alpha = 0.5, fill = "grey") +
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

Zgodnie z CI obliczonymi przy użyciu effectsopakowania, partia „E” nie pokrywa się z partią „A”.

Jeśli spróbuję tego samego przy użyciu confint.merModfunkcji i metody domyślnej:

a <- fixef(fm1)
b <- confint(fm1)
# Computing profile confidence intervals ...
# There were 26 warnings (use warnings() to see them)

b <- data.frame(b)
b <- b[-1:-2,]

b1 <- b[[1]]
b2 <- b[[2]]

dt <- data.frame(fit   = c(a[1],  a[1] + a[2:length(a)]), 
                 lower = c(b1[1],  b1[1] + b1[2:length(b1)]), 
                 upper = c(b2[1],  b2[1] + b2[2:length(b2)]) )
dt$batch <- LETTERS[1:nrow(dt)]

ggplot(dt, aes(x = batch, y = fit, ymin = lower, ymax = upper)) +
  geom_rect(xmax = Inf, xmin = -Inf, ymin = dt[dt$batch == "A", "lower"], 
        ymax = dt[dt$batch == "A", "upper"], alpha = 0.5, fill = "grey") + 
  geom_errorbar(width = 0.2) + geom_point() + theme_bw()

Widzę, że wszystkie elementy CI pokrywają się. Otrzymuję również ostrzeżenia wskazujące, że funkcja nie mogła obliczyć wiarygodnych CI. Ten przykład i mój aktualny zestaw danych każą mi podejrzewać, że effectspakiet korzysta ze skrótów w obliczeniach CI, które mogą nie zostać całkowicie zatwierdzone przez statystyków. Jak wiarygodne są elementy CI zwracane przez effectfunkcję z effectspakietu dla lmerobiektów?

Co próbowałem: patrząc na kod źródłowy zauważyłem, że effectfunkcja zależy od Effect.merModfunkcji, która z kolei kieruje do Effect.merfunkcji, która wygląda następująco:

effects:::Effect.mer
function (focal.predictors, mod, ...) 
{
    result <- Effect(focal.predictors, mer.to.glm(mod), ...)
    result$formula <- as.formula(formula(mod))
    result
}
<environment: namespace:effects>

mer.to.glmfunkcja wydaje się obliczać macierz wariancji-kowariancji z lmerobiektu:

effects:::mer.to.glm

function (mod) 
{
...
mod2$vcov <- as.matrix(vcov(mod))
...
mod2
}

To z kolei jest prawdopodobnie używane w Effect.defaultfunkcji do obliczania CI (mogłem źle zrozumieć tę część):

effects:::Effect.default
...
     z <- qnorm(1 - (1 - confidence.level)/2)
        V <- vcov.(mod)
        eff.vcov <- mod.matrix %*% V %*% t(mod.matrix)
        rownames(eff.vcov) <- colnames(eff.vcov) <- NULL
        var <- diag(eff.vcov)
        result$vcov <- eff.vcov
        result$se <- sqrt(var)
        result$lower <- effect - z * result$se
        result$upper <- effect + z * result$se
...

Nie wiem wystarczająco dużo o LMM, aby ocenić, czy jest to właściwe podejście, ale biorąc pod uwagę dyskusję dotyczącą obliczania przedziału ufności dla LMM, podejście to wydaje się podejrzanie proste.

Wszystkie wyniki są zasadniczo takie same ( dla tego konkretnego przykładu ). Niektóre różnice teoretyczne to:

• jak wskazuje @rvl, twoja rekonstrukcja CI bez uwzględnienia kowariancji między parametrami jest po prostu błędna (przepraszam)
• przedziały ufności dla parametrów mogą być oparte na przedziały ufności Wald (zakładając kwadratową powierzchnię log-prawdopodobieństwa): lsmeans, effects, confint(.,method="Wald"); z wyjątkiem lsmeanstych metod ignorują efekty o skończonej wielkości („stopnie swobody”), ale w tym przypadku nie robi to żadnej różnicy ( df=40jest praktycznie nie do odróżnienia od nieskończonego df)
• ... lub w przedziałach ufności profilu (metoda domyślna; ignoruje efekty o skończonych rozmiarach, ale pozwala na powierzchnie niekwadratowe)
• ... lub przy ładowaniu parametrycznym (złoty standard - zakłada, że ​​model jest poprawny [odpowiedzi są Normalne, losowe efekty są Normalnie rozmieszczone, dane są warunkowo niezależne itp.], ale poza tym nie ma wielu założeń)

Myślę, że wszystkie te podejścia są rozsądne (niektóre są bardziej przybliżone niż inne), ale w tym przypadku nie robi to żadnej różnicy, z której korzystasz. Jeśli jesteś zaniepokojony, wypróbuj kilka kontrastujących metod na swoich danych lub na symulowanych danych, które przypominają twoje, i zobacz, co się stanie ...

(PS: I nie umieszczać zbyt dużej wagi na fakt, że przedziały ufności Ai Enie pokrywają Trzeba zrobić właściwej procedury porównywania parami do przeprowadzenia wiarygodnych wniosków o różnicach między tym. Konkretnej pary szacunków. ..)

95% CI:

Kod porównawczy:

library(lme4)
fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
c0 <- confint(fm2,method="Wald")
c1 <- confint(fm2)
c2 <- confint(fm2,method="boot")
library(effects)
library(lsmeans)
c3 <- with(effect("batch",fm2),cbind(lower,upper))
c4 <- with(summary(lsmeans(fm2,spec="batch")),cbind(lower.CL,upper.CL))
tmpf <- function(method,val) {
    data.frame(method=method,
               v=LETTERS[1:10],
               setNames(as.data.frame(tail(val,10)),
                        c("lwr","upr")))
}
library(ggplot2); theme_set(theme_bw())
allCI <- rbind(tmpf("lme4_wald",c0),
      tmpf("lme4_prof",c1),
      tmpf("lme4_boot",c2),
      tmpf("effects",c3),
               tmpf("lsmeans",c4))
ggplot(allCI,aes(v,ymin=lwr,ymax=upr,colour=method))+
    geom_linerange(position=position_dodge(width=0.8))

ggsave("pastes_confint.png",width=10)

Wygląda na to, że to, co zrobiłeś w drugiej metodzie, to obliczenie przedziałów ufności dla współczynników regresji, a następnie przekształcenie ich w celu uzyskania CI dla prognoz. Ignoruje to kowariancje między współczynnikami regresji.

Spróbuj dopasować model bez przechwytywania, aby batchefekty były faktycznie przewidywaniami i confintzwróciły potrzebne interwały.

Dodatek 1

Zrobiłem dokładnie to, co zasugerowałem powyżej:

> fm2 <- lmer(strength ~ batch - 1 + (1 | cask), Pastes)
> confint(fm2)
Computing profile confidence intervals ...
           2.5 %    97.5 %
.sig01  0.000000  1.637468
.sigma  2.086385  3.007380
batchA 60.234772 64.298581
batchB 57.268105 61.331915
batchC 60.018105 64.081915
batchD 57.668105 61.731915
batchE 53.868105 57.931915
batchF 59.001439 63.065248
batchG 57.868105 61.931915
batchH 61.084772 65.148581
batchI 56.651439 60.715248
batchJ 56.551439 60.615248

Wydaje się, że te interwały pokrywają się z wynikami effects.

Dodatek 2

Inną alternatywą jest pakiet lsmeans . Uzyskuje stopnie swobody i dostosowaną macierz kowariancji z pakietu pbkrtest .

> library("lsmeans")
> lsmeans(fm1, "batch")
Loading required namespace: pbkrtest
 batch   lsmean       SE    df lower.CL upper.CL
 A     62.26667 1.125709 40.45 59.99232 64.54101
 B     59.30000 1.125709 40.45 57.02565 61.57435
 C     62.05000 1.125709 40.45 59.77565 64.32435
 D     59.70000 1.125709 40.45 57.42565 61.97435
 E     55.90000 1.125709 40.45 53.62565 58.17435
 F     61.03333 1.125709 40.45 58.75899 63.30768
 G     59.90000 1.125709 40.45 57.62565 62.17435
 H     63.11667 1.125709 40.45 60.84232 65.39101
 I     58.68333 1.125709 40.45 56.40899 60.95768
 J     58.58333 1.125709 40.45 56.30899 60.85768

Confidence level used: 0.95 

effecteffectconfint$\pm1.96\times\mbox{se}$

Wyniki effecti lsmeanssą podobne, ale z niesymetrycznym sytuacji wielu czynników, lsmeansprzez standardowe średnich niewykorzystane czynników o równych masach, podczas effectwag obserwowana częstotliwości (dostępny jako opcja lsmeans).