title: Confidence Interval in mixed effect models
body: <pre><code>library(lme4)
  
  fm1 &lt;- lmer(Reaction ~ Days + (Days|Subject), data = sleepstudy)
  </code></pre>
  
  <p>To generate a 95% CI, I can use the <code>predictInterval()</code> function from the package <code>merTools</code>.</p>
  
  <pre><code>library(merTools)
  head(predictInterval(fm1, level = 0.95, seed = 123, n.sims = 100))
  # fit      upr      lwr
  # 1 255.4179 313.8781 184.1400
  # 2 273.2944 333.2005 231.3584
  # 3 291.8451 342.8701 240.8226
  # 4 311.3562 359.2908 250.4980
  # 5 330.3671 384.2520 270.7094
  # 6 353.4378 409.9307 289.4760
  </code></pre>
  
  <p>In the documentation, it says about the <code>predictInterval()</code> function</p>
  
  <blockquote>
    <p>This function provides a way to capture model uncertainty in predictions from multi-level models
         fit with lme4. By drawing a sampling distribution for the random and the fixed effects and then
         estimating the fitted value across that distribution, it is possible to generate a prediction interval for
         fitted values that includes all variation in the model except for variation in the covariance parameters,
         theta. This is a much faster alternative than bootstrapping for models fit to medium to large datasets.</p>
  </blockquote>
  
  <p>My goal is to get all the fitted values instead of the the upper and lower CI i.e. for each row, I need the 
  original n simulations from which these 95% CI are calculated. I checked the argument in the documentation and
  followed this:  </p>
  
  <pre><code>head(predictInterval(fm1, n.sims = 100, returnSims = TRUE, seed = 123, level = 0.95))
  # fit      upr      lwr
  # 1 255.4179 313.8781 184.1400
  # 2 273.2944 333.2005 231.3584
  # 3 291.8451 342.8701 240.8226
  # 4 311.3562 359.2908 250.4980
  # 5 330.3671 384.2520 270.7094
  # 6 353.4378 409.9307 289.4760
  </code></pre>
  
  <p>Instead of getting the 100 simulations, it still gives me the same output. What is it I am doing wrong here?</p>
  
  <p>A second question though I believe this is more of a StatsExchange one.     </p>
  
  <blockquote>
    <p>"By drawing a sampling distribution for the random and the fixed
    effects and then."`</p>
  </blockquote>
  
  <p>How does it draws the sampling distribution if some could explain me?        </p>