Загрузочные доверительные интервалы для прогнозируемых данных перекрываются, когда фиксированные эффекты не совпадают

Примечание. Это может быть более подходящим для перекрестной проверки, и я могу переместить его туда, если потребуется, но я решил сначала попробовать здесь, поскольку это может быть связано с r.

Я сравниваю оценки параметров кривой роста между двумя группами, используя нелинейные смешанные модели с nlme. Я использую boot_nlme() из nlraa для начальной загрузки доверительных интервалов для оценок параметров модели и прогнозируемых данных.

Ниже приведен код, который я использую, но с penguin.data из FlexParamCurve. Он не имеет такой же проблемы, но предоставляет код в удобном формате. Если этого недостаточно, я могу редактировать свои собственные данные.

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

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

Любая помощь/понимание будет оценено по достоинству.

library(FlexParamCurve) #loads dataset and  package 'nlme'
library(ggplot2)
library(nlraa)
library(boot)
library(dplyr)
library(car)



VarFunc.Auto<-varPower(form=~fitted(.))

##creating model

richards.func <- function(age, A, Ti, k, d){
  A * (1 + (d - 1) * exp(( (-k) * (age - Ti)) / ( d ^ ( d / (1 - d))))) ^ (1 / (1 - d))
}

ggplot(penguin.data, aes (x = ckage, y = weight, color = year)) +
  geom_smooth()

#fixed asymptote
mod <- nlme(weight ~ SSlogis(ckage, Asym, R0, lrc),
            data = penguin.data,
            fixed= list(Asym ~ year,
                        R0 ~ year,
                        lrc ~ year),
            random = Asym ~ 1,
            start = c(Asym = 1000, 0,
                      R0 = 21, 0,
                      lrc = 1, 0),
            control = list(maxIter = 100),
            na.action = na.pass)


summary(mod)

peng.mdiff <- function(x, ckage = seq(0, 80, length.out = 500)){
  ndat <- expand.grid(ckage = ckage, 
                      year = c('2000','2002'), 
                     nest = NA,
                     bandid = NA,
                      stringsAsFactors = TRUE)
  
  prd <- predict(x, newdata = ndat, level = 0)

}
set.seed(123)
#this takes a few minutes to run on my computer
system.time(peng.bt <- boot_nlme(mod, peng.mdiff, R = 500, cores = 3))
set.seed(123)
system.time(peng.bt.ci <- confint(peng.bt, level = 0.95))


prd.df <- as.data.frame(peng.bt.ci)

prd.df <- prd.df %>%
  rename_at('2.5 %', ~'lower.ci') %>%
  rename_at('97.5 %', ~'upper.ci')

ndat1 <- with(penguin.data,
              expand.grid(
                ckage=seq(0, 80 , length.out=500),
                year = c('2000','2002'), 
                nest = NA,
                bandid = NA
              ))

newX <- ndat1 %>%
  mutate(prd = predict_nlme(mod, newdata = ndat1, level = 0))

comb.df1 <- cbind(prd.df, ndat1, 'prd' = newX$prd)

ci.plot <- ggplot() +
  geom_ribbon(comb.df1, mapping = aes(x = ckage, ymin = lower.ci, ymax = upper.ci, fill = year ), alpha = 0.50) +
  geom_line(data = newX, aes(x = ckage, y = prd, group = year), color = 'black', alpha = 0.75)

ci.plot