title: How to get standard deviations using predict?
body: <p>I am trying to plot a three-way interaction. First, i use <code>expand.grid</code> to create different combinations of my variables. Then i use <code>predict</code> to get the predicted values. It works well, however, i would like to get the standard deviation or confidence interval. Once i include the option <code>se.fit = T</code> in <code>predict</code>, i get back some results that i cannot interpret as two combinations get one predicted value and the other two get a range of four outcomes. I am not sure if the <code>se.fit = T</code> option is wrong in this context or whether the issue is in my interpretation of the results. Could someone help please?</p>
  <p>This is the code that works without the <code>se.fit</code> option:</p>
  <pre><code>set.seed(20170925) 
  dat &lt;- data.frame(time=gl(n = 2,k = 5000),
                    y= rnorm(10000, mean=1000, sd=400),
                    a=factor(sample(1:2, 10000, replace=TRUE)),
                    b=factor(sample(1:2, 10000, replace=TRUE)),
                    c=factor(rep(rep(1:3,times=c(2500,2500,5000)),2)))
  
  m = lm(y ~ a*b*c, data=dat)
  pred &lt;- expand.grid(a=factor(1:2),b=factor(1:2), c=factor(1, levels=c(1,2,3)))
  pred$y &lt;- predict(m,pred)
  </code></pre>
  <p>To add the <code>se.fit</code> option, we run the same code but change the last line to:</p>
  <pre><code>pred$y &lt;- predict(m,pred, se.fit = T)
  </code></pre>