title: How to specify all levels for R{stats} predict() function in non linear mixed model after model fit with package medrc
body: <p>I have 3 trials (trial: e1, e2, e3), 2 products/trial (products: A, B), 5 rates/product (.1,1,10,100,1000), total of 6 curves (curve: c1,...,c6).
  After fitting a non linear mixed model, I want to plot all the curves and the curves resulting from the model in the same chart.
  This is the reference (package medrc in github): <a href="https://doseresponse.github.io/medrc/articles/medrc.html" rel="nofollow noreferrer">https://doseresponse.github.io/medrc/articles/medrc.html</a></p>
  
  <p>This is the code to generate the non-linear mixed model.</p>
  
  <pre><code>#packages
  library(drc)
  library(medrc)
  library(dplyr)
  library(tidyr)
  
  #my data
  trial &lt;- c("e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1",
             "e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1","e1",
             "e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2",
             "e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2","e2",
             "e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3",
             "e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3","e3")
  
  curve &lt;- c("c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1","c1",
             "c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2","c2",
             "c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3","c3",
             "c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4","c4",
             "c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5","c5",
             "c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6","c6")
  
  rates &lt;- c(.1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
             .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
             .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
             .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
             .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000,
             .1,.1,.1,1,1,1,10,10,10,100,100,100,1000,1000,1000)
  
  product &lt;- c("A","A","A","A","A","A","A","A","A","A","A","A","A","A","A",
               "B","B","B","B","B","B","B","B","B","B","B","B","B","B","B",
               "A","A","A","A","A","A","A","A","A","A","A","A","A","A","A",
               "B","B","B","B","B","B","B","B","B","B","B","B","B","B","B",
               "A","A","A","A","A","A","A","A","A","A","A","A","A","A","A",
               "B","B","B","B","B","B","B","B","B","B","B","B","B","B","B")
  
  resp &lt;- c(.295,.3232,.3015,.155,.1501,.1483,.0511,.036,.0445,.0021,.0022,.0035,.0015,.0025,.0009,         
        .312,.3373,.2994,.265,.2501,.2482,.1022,.103,.1142,.0220,.0198,.0159,.0036,.0099,.0100,
        .289,.3122,.3093,.141,.1612,.1398,.0722,.022,.0581,.0019,.0015,.0011,.0018,.0009,.0014,
        .325,.3451,.2952,.267,.2412,.2398,.1125,.109,.1019,.0554,.0547,.0118,.0029,.0075,.0078,
        .294,.2452,.2991,.121,.1925,.1485,.0871,.025,.0658,.0019,.0019,.0010,.0025,.0019,.0008,
        .285,.3412,.3069,.124,.1861,.1958,.1276,.132,.1985,.0325,.0201,.0225,.0031,.0089,.0094)
  
  
  data.test &lt;- data.frame(trial,curve,rates,product,resp) #my data frame
  
  #my model
  m1 &lt;- medrm(resp ~ rates, 
              curveid=b + c + d + e ~ product, 
              data = data.test, 
              fct=LL.4(), 
              random = c + d ~ 1|trial,
              start=NULL)
  </code></pre>
  
  <p>To make the plot:</p>
  
  <pre><code>#plotting
  pdata &lt;- data.test%&gt;%
    group_by(curve, product) %&gt;%
    expand(rates=exp(seq(-3, 10, length=50)))
  #pdata$resp_ind &lt;- predict(m1, newdata=pdata)
  pdata$resp &lt;- predict(m1, newdata=pdata, level=0)
  
  ggplot(data.test, aes(x=log(rates), y=resp, 
                        colour=product, group=curve, shape=product)) +
    geom_point() +
    geom_line(data=pdata) +
    #geom_line(data=pdata, aes(y=resp_ind), linetype=2) +
    theme_bw() +
    scale_x_continuous("DOSE", 
                       breaks=log(c(.1, 1, 10, 100, 1000)), 
                       labels=c(.1, 1, 10, 100, 1000))
  </code></pre>
  
  <p>Note that two rows of code are commented.
  <strong>When extracting the predict data per curve, I cannot specify the levels, i.e. the curves that gave the random components. What am I missing?</strong></p>
  
  <pre><code>pdata$resp_ind &lt;- predict(m1, newdata=pdata)
  </code></pre>
  
  <p>is resulting the error:</p>
  
  <pre><code>Error in predict.nlme(object$fit, newdata = newdata, level = level) : 
  cannot evaluate groups for desired levels on 'newdata'
  </code></pre>
  
  <p>So I cannot use this row of code to plot each curve line</p>
  
  <pre><code>geom_line(data=pdata, aes(y=resp_ind), linetype=2) +
  </code></pre>
  
  <p>These are similar questions, but I'm getting the average trend with the code:</p>
  
  <pre><code>pdata$resp &lt;- predict(m1, newdata=pdata, level=0)
  </code></pre>
  
  <p><strong>I wanted to specify the levels to get all curves.</strong>
  <a href="https://stackoverflow.com/questions/28684645/r-lme-cannot-evaluate-groups-for-desired-levels-on-newdata">R: lme, cannot evaluate groups for desired levels on &#39;newdata&#39;</a>
  <a href="https://stats.stackexchange.com/questions/58031/prediction-on-mixed-effect-models-what-to-do-with-random-effects">https://stats.stackexchange.com/questions/58031/prediction-on-mixed-effect-models-what-to-do-with-random-effects</a></p>