title: How to obtain confidence intervals for non-linear mixed model (logistic growth curve) with an interaction
body: <p>I'd like to calculate confidence intervals using the delta method, or bootstrapping, for a non-linear model (logistic growth curve) with a two-way interaction, but am unsure of how to code this. This post is similar to this one 
  <a href="https://stats.stackexchange.com/questions/231074/confidence-intervals-on-predictions-for-a-non-linear-mixed-model-nlme">https://stats.stackexchange.com/questions/231074/confidence-intervals-on-predictions-for-a-non-linear-mixed-model-nlme</a>. However, I do not know how to modify it to deal with the interaction term. I'm using the ChickWeight dataframe, but I've modified it so that the curves are more sigmoidal/looks similar to a dataframe I'm using. </p>
  
  <p>So far, I've tried bootstrapping with no success. </p>
  
  <pre><code>library(tidyverse)
  library(broom)
  library(nlme)
  
  # Modify df so that there are more subjects (Chicks) per treatment (Diet) 
  Chicks.df &lt;- ChickWeight
  Chicks.df$Chick &lt;- factor(Chicks.df$Chick, ordered=F)
  Chicks.df &lt;- Chicks.df %&gt;% 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 &lt;- rnorm(30, mean=178, sd=1) # new observations for Diet1
  t.25_diet2 &lt;- rnorm(20, mean=215, sd=1) # new observations for Diet2
  weight &lt;- c(t.23_diet1, t.23_diet2) # bind vectors together
  Diet &lt;- c(rep(1,30), rep(2,20)) # Create remaining variables for df
  Chick &lt;- c(seq(1:30), seq(31, 50)) 
  Time &lt;- rep(25, 50)
  
  # Bind variables together to make new df, and then combine with original df (i.e. Chicks.df) 
  newdata &lt;- data.frame(cbind(weight, Diet, Chick, Time))
  newdata$Diet &lt;- as.factor(newdata$Diet)
  newdata$Chick &lt;- as.factor(newdata$Chick)
  Chicks.df &lt;- bind_rows(Chicks.df, data.frame(newdata))
  
  # Using nested ifelse, assign half of the individuals in each diet to new treatment
  Chicks.df &lt;- Chicks.df %&gt;% mutate(Drug = 
                                  ifelse(Chick %in% c(1:15) &amp; 
                                           Diet == 1, 0, 
                                         ifelse(Chick %in% c(16:30) &amp; 
                                                Diet == 1, 1, 
                                                ifelse(Chick %in% 
                                                         c(31:40) &amp; Diet
                                                       == 2, 0, 1))))
  # Check Curves
  Chicks.df$Drug &lt;- 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 %&gt;% 
       group_by(Diet, Drug) %&gt;% 
       do(fit = nls(weight ~ SSlogis(Time, Asym, xmid, scal), data = .)) %&gt;% 
       tidy(fit) %&gt;% 
       select(Diet, Drug, term, estimate) %&gt;% 
       spread(term, estimate) 
  
       # Run model 
        mod &lt;-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 &lt;- data.frame(Time = Chicks.df$Time, 
                      Diet = Chicks.df$Diet,
                      Drug = Chicks.df$Drug,
                      Chick = Chicks.df$Chick[1])
     pframe$weight &lt;- predict(mod, pframe, re.form=~0) 
  
     myfun &lt;- function(data, x){
              fBoot &lt;- 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 &lt;- boot::boot(pframe, myfun, R=100)
       myboot # no values predicted
  </code></pre>