Hi,
I am trying to run a multinomial model with a car structure. I am using car to account for spatial correlation between columns and rows for a field trial. If I do not have car(), I usually have no problems to use compare_levels() with the multinomial family. For this model in particular, I did not figure out how to make it work for getting contrasts for my explanatory variable. Contrast for the response variables seems to work with a few modifications. I have a sample code that reproduce the error.
df2<-structure(list(rep = c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L,
6L, 6L), spatial_id = c("6_4", "6_2", "5_10", "5_9", "4_6", "4_5",
"3_1", "3_15", "2_15", "2_13", "1_3", "1_1"), row = structure(c(6L,
6L, 5L, 5L, 4L, 4L, 3L, 3L, 2L, 2L, 1L, 1L), levels = c("1",
"2", "3", "4", "5", "6"), class = "factor"), column = c(1, 2,
2, 1, 1, 2, 2, 1, 1, 2, 2, 1), total = c(48, 42, 53, 29, 74,
39, 140, 55, 79, 55, 28, 92), var1 = c("x1", "x2", "x1", "x2",
"x1", "x2", "x1", "x2", "x1", "x2", "x1", "x2"), a = c(4L, 5L,
3L, 1L, 9L, 5L, 10L, 6L, 8L, 5L, 2L, 10L), b = c(19L, 18L, 24L,
13L, 27L, 12L, 61L, 22L, 33L, 30L, 13L, 42L), c = c(25L, 19L,
26L, 15L, 38L, 22L, 69L, 27L, 38L, 20L, 13L, 40L), y = structure(c(4L,
5L, 3L, 1L, 9L, 5L, 10L, 6L, 8L, 5L, 2L, 10L, 19L, 18L, 24L,
13L, 27L, 12L, 61L, 22L, 33L, 30L, 13L, 42L, 25L, 19L, 26L, 15L,
38L, 22L, 69L, 27L, 38L, 20L, 13L, 40L), dim = c(12L, 3L), dimnames = list(
NULL, c("a", "b", "c")))), row.names = c(NA, -12L), class = c("tbl_df",
"tbl", "data.frame"))
Grid <- expand.grid(unique(df2$row), unique(df2$column))
K <- nrow(Grid)
# set up distance and neighbourhood matrices
distance <- as.matrix(dist(Grid))
W <- array(0, c(K, K))
W[distance == 1] <- 1
rownames(W) <- df2$spatial_id
colnames(W) <- df2$spatial_id
brm.ly2 <- brm(y|trials(total) ~ var1+ car(W,gr=spatial_id),
family= brms::multinomial() ,
data = df2,
data2 = list(W = W),
core=4)
# Point estimates
est_<-df2%>%
add_epred_draws(brm.ly2)%>%group_by(var1,.category)%>%mean_qi(.epred)
# I would usually use data.frame(var1=unique(df2$var1)%>%add_epred_draws(brm.ly2)%>%compare_levels(variable =.epred, by =var1)%>%mean_qi()
#For the response it worked like that:
df2%>%
add_epred_draws(brm.ly2)%>%compare_levels(variable =.epred, by =.category)%>%group_by(var1,.category)%>%
mean_qi(.epred)
#For my explanatory it is not working:
df2%>%
add_epred_draws(brm.ly2)%>%compare_levels(variable =.epred, by =var1)%>%group_by(var1,.category)%>%
mean_qi(.epred)
#> Error in combn(levels(x), 2, simplify = FALSE) : n < m
# So I tried:
options("marginaleffects_posterior_center" = "mean")
c_ly3<-avg_comparisons(brm.ly2,variables = "var1") #the estimates do not make sense
est_[est_$var1 == "x2" & est_$.category == "a", ]$.epred-est_[est_$var1 == "x1" & est_$.category == "a", ]$.epred
est_[est_$var1 == "x2" & est_$.category == "b", ]$.epred-est_[est_$var1 == "x1" & est_$.category == "b", ]$.epred
- Operating System: windows 11 enterprise
- brms Version: 2.22.0
I really appreciate your help.