У меня есть многоуровневая отрицательная биномиальная модель, подходящая для brms (library(brms))
fit1 <- brm(TOTAL_VIOLATIONS ~ LN_POP + Source_binary + Source_purchased + (1|TYPE_consolidated) + (1|COUNTY), data = Data, family = negbinomial())
Вот как выглядят данные:
> dput(droplevels(Data[1:20, c(3, 9, 20, 21, 22, 23)]))
structure(list(COUNTY = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L), .Label = c("ALAMEDA",
"ALPINE", "AMADOR"), class = "factor"), TYPE_consolidated = structure(c(9L,
6L, 3L, 2L, 5L, 7L, 1L, 1L, 4L, 12L, 1L, 1L, 1L, 1L, 8L, 10L,
6L, 5L, 11L, 2L), .Label = c("City", "County Water District",
"CSA", "CSD", "IOU", "MHP", "MUD", "MWC", "PA", "Private", "PUD",
"Special Act District"), class = "factor"), TOTAL_VIOLATIONS = c(0L,
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2L, 0L, 0L,
8L, 0L, 0L), LN_POP = c(3.91202300542815, 6.29710931993394, 6.21260609575152,
12.7367008965923, 10.9852927228879, 14.1374128813017, 11.9290007521904,
11.1991321074213, 11.2374881189349, 12.332000202128, 10.2255710517052,
6.10924758276437, 6.62007320653036, 6.21460809842219, 3.91202300542815,
3.2188758248682, 4.24849524204936, 7.88231491898027, 8.96839619119826,
4.91265488573605), Source_binary = structure(c(1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L
), .Label = c("GW", "SW"), class = "factor"), Source_purchased = structure(c(1L,
1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L, 2L), .Label = c("No", "Purchased"), class = "factor")), row.names = c(NA,
20L), class = "data.frame")
Я хочу изучить значение включения первого случайного эффекта (TYPE_consolidated) по сравнению с той же моделью только со вторым случайным перехватом (COUNTY), но я не могу, хоть убей, выяснить, как получить icc(), чтобы сообщить об этом информация по группам с использованием by_group. Вывод будет таким же, как с этим аргументом, так и без него (см. Ниже). У меня есть ощущение, что это потому, что это объект brms, поскольку, согласно справке, вывод для этих моделей отличается, но мне еще предстоит придумать другой способ получить это. Кто-нибудь знает способ получить коэффициент дисперсии для отдельных случайных эффектов (или, альтернативно, может увидеть, что я делаю не так с by_group)? Если нет, существует ли стандартный способ сравнения I CC между вложенными моделями? Если это так, возможно, я мог бы просто рассчитать это для двух разных версий моей модели и сделать это вместо этого?
> performance::icc(fit1, by_group = TRUE)
# Random Effect Variances and ICC
Conditioned on: all random effects
## Variance Ratio (comparable to ICC)
Ratio: 0.94 CI 95%: [0.80 0.99]
## Variances of Posterior Predicted Distribution
Conditioned on fixed effects: 7.39 CI 95%: [ 2.50 20.03]
Conditioned on rand. effects: 117.57 CI 95%: [59.15 331.86]
## Difference in Variances
Difference: 109.10 CI 95%: [50.23 320.86]
> performance::icc(fit1)
# Random Effect Variances and ICC
Conditioned on: all random effects
## Variance Ratio (comparable to ICC)
Ratio: 0.94 CI 95%: [0.79 0.99]
## Variances of Posterior Predicted Distribution
Conditioned on fixed effects: 7.42 CI 95%: [ 2.48 20.19]
Conditioned on rand. effects: 117.90 CI 95%: [59.53 349.90]
## Difference in Variances
Difference: 109.71 CI 95%: [51.20 340.30]```