Как получить доверительные интервалы для нелинейной смешанной модели (кривая логистического роста) с взаимодействием

Question

Я хотел бы рассчитать доверительные интервалы, используя дельта-метод или начальную загрузку для нелинейной модели (кривая логистического роста) с двусторонним взаимодействием, но я не уверен, как это кодировать.Этот пост похож на этот https://stats.stackexchange.com/questions/231074/confidence-intervals-on-predictions-for-a-non-linear-mixed-model-nlme. Тем не менее, я не знаю, как изменить его, чтобы иметь дело с термином взаимодействия.Я использую фрейм данных ChickWeight, но я изменил его, чтобы кривые были более сигмоидальными / выглядели аналогично фрейму данных, который я использую.

До сих пор я пробовал загрузку безуспешно.

library(tidyverse)
library(broom)
library(nlme)

# Modify df so that there are more subjects (Chicks) per treatment (Diet) 
Chicks.df <- ChickWeight
Chicks.df$Chick <- factor(Chicks.df$Chick, ordered=F)
Chicks.df <- Chicks.df %>% mutate(Diet = ifelse(Diet == 3 | Diet == 4, 2, 1),
                              Diet = as.factor(Diet))

# Create df with an additional time point (25th) to make curves more sigmoidal
t.25_diet1 <- rnorm(30, mean=178, sd=1) # new observations for Diet1
t.25_diet2 <- rnorm(20, mean=215, sd=1) # new observations for Diet2
weight <- c(t.23_diet1, t.23_diet2) # bind vectors together
Diet <- c(rep(1,30), rep(2,20)) # Create remaining variables for df
Chick <- c(seq(1:30), seq(31, 50)) 
Time <- rep(25, 50)

# Bind variables together to make new df, and then combine with original df (i.e. Chicks.df) 
newdata <- data.frame(cbind(weight, Diet, Chick, Time))
newdata$Diet <- as.factor(newdata$Diet)
newdata$Chick <- as.factor(newdata$Chick)
Chicks.df <- bind_rows(Chicks.df, data.frame(newdata))

# Using nested ifelse, assign half of the individuals in each diet to new treatment
Chicks.df <- Chicks.df %>% mutate(Drug = 
                                ifelse(Chick %in% c(1:15) & 
                                         Diet == 1, 0, 
                                       ifelse(Chick %in% c(16:30) & 
                                              Diet == 1, 1, 
                                              ifelse(Chick %in% 
                                                       c(31:40) & Diet
                                                     == 2, 0, 1))))
# Check Curves
Chicks.df$Drug <- as.factor(Chicks.df$Drug)
ggplot(Chicks.df, aes(Time, weight, col=Diet, linetype=Drug)) + geom_smooth()

# Get self-starting values (use broom pkg)
 Chicks.df %>% 
     group_by(Diet, Drug) %>% 
     do(fit = nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = .)) %>% 
     tidy(fit) %>% 
     select(Diet, Drug, term, estimate) %>% 
     spread(term, estimate) 

     # Run model 
      mod <-nlme(weight ~ SSlogis(Time, Asym, xmid, scal),
      fixed=list(Asym + xmid + scal ~ Diet*Drug),
      random = Asym ~ 1 | Chick,
      start=list(fixed=c(Asym=c(209, 215, 209, 249), 
                            xmid=c(10.2, 10.5, 10.7, 9.7), 
                            scal=c(6.3, 6.3, 5.1, 5.5))),
         data=Chicks.df)

   summary(mod)

   # Build prediction df and attempt bootstrapping 
   pframe <- data.frame(Time = Chicks.df$Time, 
                    Diet = Chicks.df$Diet,
                    Drug = Chicks.df$Drug,
                    Chick = Chicks.df$Chick[1])
   pframe$weight <- predict(mod, pframe, re.form=~0) 

   myfun <- function(data, x){
            fBoot <- nlme(weight ~ SSlogis(Time, Asym, xmid, scal),
      fixed=list(Asym + xmid + scal ~ Diet*Drug),
      random = Asym ~ 1 | Chick,
      start=list(fixed=c(Asym=c(209, 215, 209, 249), 
                            xmid=c(10.2, 10.5, 10.7, 9.7), 
                            scal=c(6.3, 6.3, 5.1, 5.5))),
         data=Chicks.df)
            return((predict(fBoot, newdata=pframe)))
          }

     myboot <- boot::boot(pframe, myfun, R=100)
     myboot # no values predicted