ggpredict for mclogit model predicts high CIs are low or 0, low CI high

[szojka]

The issue is that confidence intervals predicted by ggpredict are incorrect. My multinomial model is fit with mclogit and has nested random effects, so below I've attached a reproducible example of the problem with simulated data of the same structure.

#################

# load packages
library(tidyverse)
library(mclogit)
library(ggeffects)
set.seed(222)

#################

#1. create fake data with nested random effects and dependent response variables

# random effects have nested structure: site > grid > block 
site <- c(rep("a",times = 150), 
          rep("b",times = 150),
          rep("c",times = 150),
          rep("d",times = 150),
          rep("e",times = 150),
          rep("f",times = 150)) 
grid <- rep(c(1:18), times = 50) # 3 grids within each site 6*3=18
grid <- sort(grid) # keep reps near one another
block <- rep(c(1:450), times = 2) # 25 blocks per grid 18*25 = 450. Block links treatments 
block <- sort(block) # keep reps near one another

# treatments are repeated within blocks
treatment <- c(rep(c("W","Z"), times= 450)) # 2 treatments within each block 450*2 = 900

# x values stay the same per treatment within block
x_seq <- seq(from = -4,to = 2, length.out = 450)
x_seq <- sample(x_seq) # randomize order of x
x_seq <- rep(x_seq, each = 2)# repeat each number in this sequence and  keep same numbers next to one another

# response variables are dependent
resp1 <- c(sample((rep(c(0,1), length.out = 900))))
resp2 <- c(sample((rep(c(0,1), length.out = 900))))
resp3 <- c(sample((rep(c(0,1), length.out = 900))))
resp4 <- c(sample((rep(c(0,1), length.out = 900))))
identical(resp1, resp2) # should be false

# bring into dataframe
dat <- data.frame(site=site, grid=grid, block=block, 
                  treatment=treatment, 
                  x = x_seq,
                  resp1=resp1, resp2=resp2, resp3=resp3, resp4=resp4) 

# make response variables into proportions. 
dat <- dat %>%
  mutate(sum = rowSums(.[6:9]))

dat$resp1 <- dat$resp1/dat$sum
dat$resp2 <- dat$resp2/dat$sum
dat$resp3 <- dat$resp3/dat$sum
dat$resp4 <- dat$resp4/dat$sum

dat[is.na(dat)] <- 0

# make random effects factors

dat$site <- as.factor(dat$site)
dat$grid <- as.factor(dat$grid)
dat$block <- as.factor(dat$block)
dat$treatment <- as.factor(dat$treatment)

#select response variables (proportion data)

y_mat <-  as.matrix(dat[, 6:9]) # make response matrix

#################

#2. fit multinomial model

mfit <- mblogit(y_mat ~ x*treatment, random = c(~ 1|site/grid, ~1|block), data = dat,
                control = mmclogit.control(epsilon = 1e-08,
                                           maxit = 50, trace=TRUE,
                                           trace.inner=FALSE,
                                           avoid.increase = FALSE,
                                           break.on.increase = FALSE,
                                           break.on.infinite = FALSE,
                                           break.on.negative = FALSE))

#################

#3. visualize problem using ggeffects

vis_mfit <-
  ggpredict(mfit,
            terms = c("x[all]", "treatment"),
            type = "fe",
            ci.lvl = 0.95,
            interval = "confidence"); plot(vis_mfit)

head(vis_mfit)

# PROBLEM: 95% confidence intervals do not make sense -  high values are low, low values are high.