Генерация доверительных интервалов предельного прогноза из объекта glmer с помощью predInterval () из merTools

Сопровождающий пакета для merTools здесь. То, как мы реализовали эту функцию, непросто, но это возможно.

Вам нужно добавить шаг, чтобы добавить средний случайный эффект к вашему data.frame. В большинстве случаев средний случайный эффект должен быть равен 0 или достаточно близок, чтобы приблизительно соответствовать тому, что вы ищете. Для этого вам нужно просто немного изменить код и использовать функцию REquantile из merTools:

medEff = REquantile(mod, quantile = 0.5, 
                    groupFctr = "ID", 
                    term = "(Intercept)")

PredDat <- data.frame(ELEVATION = seq(-1.97056, 2.52926, length.out = 1000),
                      SLOPE = 0, ET = 0, HLI = 0, ID = medEff)

Preds <- predictInterval(mod, newdata = PredDat, type = "probability")

Это дает прогноз, но включает неопределенность в отношении случайного эффекта, включая средний случайный эффект, равный 0. В приведенном выше примере это приводит к размыванию эффекта переменной ELEVATION по наблюдениям, поскольку средний случайный эффект оценивается не очень точно. Итак, это может быть не то, что вам нужно.

Кроме того, если у вас есть более сложная спецификация случайных эффектов с наклонами и пересечениями, тогда этот подход становится сложнее, потому что средний эффект для точки пересечения может быть равен 0, но не для наклона.

Если вы действительно хотите просто зафиксировать дисперсию в прогнозе, основанную только на фиксированных эффектах и ​​их неопределенности - что-то с тех пор, как я узнал о создании пакета, довольно распространено - есть способы сделать это за пределами merTools. Это не самый элегантный вариант, но именно это происходит под капотом predictInterval, чтобы получить изменчивость прогнозов фиксированного эффекта:

PredDat <- data.frame(Intercept = 1, 
           ELEVATION = seq(-1.97056, 2.52926,length.out = 1000), 
           SLOPE = 0, ET = 0, HLI = 0)

 fe.tmp <- fixef(mod)
 vcov.tmp <- as.matrix(vcov(mod))
 n.sims <- 1000
 sigmahat <- rep(1, n.sims)

 # Make n.sims draws for each element of the fixed effects

 betaSim <- abind::abind(lapply(1:n.sims,
  function(x) mvtnorm::rmvnorm(n = 1, mean = fe.tmp, 
       sigma = sigmahat[x]*vcov.tmp, method = "chol")), along=1)
# Calculate n.sims predictions for each row in PredDat
fixed <- as.matrix(PredDat) %*% t(betaSim)
# For each row (observation) in PredDat calculate the median, upr and lwr 
Preds <- data.frame(fit = apply(fixed, 1, median), 
                upr = apply(fixed, 1, quantile, 0.9), 
                lwr = apply(fixed, 1, quantile, 0.1))
# Calculate the probability from the linear predictor
Preds <- apply(Preds, 2, invlogit)

У вас должно получиться что-то вроде этого:

head(Preds)

     fit       upr       lwr
1 0.1860053 0.2482220 0.1427370
2 0.1860058 0.2482226 0.1427373
3 0.1860062 0.2482231 0.1427377
4 0.1860066 0.2482237 0.1427380
5 0.1860071 0.2482242 0.1427384
6 0.1860075 0.2482248 0.1427388

Это не включает какую-либо неопределенность на уровне наблюдения, связанную с вариациями факторов группирования или самой модели.