Как я могу получить матрицу дисперсий, связанных со случайными эффектами (LMER) - PullRequest
Новая
       92

Как я могу получить матрицу дисперсий, связанных со случайными эффектами (LMER)

0 голосов
/ 31 октября 2019

Я пытался получить матрицу, сообщающую значения (оценочной) дисперсии, связанной как со случайными факторами, так и с остатками, из линейной модели смешанных эффектов. Это позволяет получить CI (конфиденциальный интервал) этих случайных факторов для последующего сравнения различных наборов данных.

Схема выборки является двусторонней и состоит из: коэффициента площадки (случайный, 3 уровня) и коэффициента площади (случайный), вложенный в сайт, 4 уровня). Были отобраны 4 копии для каждого Района. Набор данных: 

>head(dataset)
          Variable  Site Area Replica
          90        A    1    a
          65        A    1    b
          65        A    1    c
          55        A    1    d
          55        A    2    a
          50        A    2    b
42 more rows

Линейная модель смешанных эффектов: 

lm=(lmer(Variable~(1|Site)+(1|Site:Area),dataset))`

>summary(lm)
...
...
    Min      1Q  Median      3Q     Max
-1.7826 -0.1301 -0.0461 -0.0461  4.5427
Random effects:
 Groups    Name        Variance Std.Dev.
 Site:Area (Intercept)  2.951   1.718
 Site      (Intercept) 17.536   4.188
 Residual              33.160   5.758
Number of obs: 48, groups:  Site:Area, 12; Site, 3
etc...

Я правильно извлек дисперсии как остаточных, так и двух случайных эффектов с помощью:

>Var_Residual<- attr(VarCorr(lm), "sc")^2

>Var_Area<-summary(lm)$varcor[[1]][[1]]

>Var_Site<-summary(lm)$varcor[[2]][[1]]

и, наконец, я загрузился для B = 10000, чтобы получить 10000 значений каждого из этих отклонений, используя формулу:

Bootstrap_residual<-matrix(sample(Var_Residual <- attr(VarCorr(lm),"sc")^2,size=1*B,,replace=TRUE),nrow=1,ncol=B)

и т. Д. Для Area и Site тоже.

Но когда я проверил загрузочные значения, я увидел ненормальное расхождение (форма была прямоугольной)) и, кроме того, мои наблюдаемые значения (2.95, 17.53, 33.16) всегда были одними из более высоких (это для многих других наборов данных и переменных, к которым я применил тот же конвейер)

>mean(attr(VarCorr(lm), "sc")^2<=Bootstrap_residual) [1] 0.005

Итак, мой вопрос: какой правильный сценарий для получения матриц (или одной матрицы), сообщающих о значении дисперсии (связанной как с остаточным, так и со случайным коэффициентом) после начальной загрузки моего набора данных?

Спасибо вамзаранее

Ferrante

...