title: Using clustered covariance matrix in predict.lm()
body: <p>I am analyzing a dataset in which data is clustered in several groups (towns in regions). The dataset looks like: </p>
  
  <pre><code>R&gt; df &lt;- data.frame(x = rnorm(10), 
                       y = 3*rnorm(x), 
                       groups = factor(sample(c('0','1'), 10, TRUE)))
  R&gt; head(df)
          x     y groups
  1 -0.8959  1.54      1
  2 -0.1008 -2.73      1
  3  0.4406  0.44      0
  4  0.0683  1.62      1
  5 -0.0037 -0.20      1
  6 -0.8966 -2.34      0
  </code></pre>
  
  <p>I want my lm() estimates to account for intraclass correlation in groups and for that purpose I am using a function <code>cl()</code> that takes an <code>lm()</code> and returns the robust clustered covariance matrix (original <a href="http://people.su.se/~ma/clustering.pdf" rel="noreferrer">here</a>):</p>
  
  <pre><code>cl  &lt;- function(fm, cluster) {
    library(sandwich)
    M &lt;- length(unique(cluster))   
    N &lt;- length(cluster)              
    K &lt;- fm$rank                   
    dfc &lt;- (M/(M-1))*((N-1)/(N-K-1))
    uj  &lt;- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum));
    vcovCL &lt;- dfc * sandwich(fm, meat = crossprod(uj)/N)
    return(vcovCL)
  }
  </code></pre>
  
  <p>Now, </p>
  
  <pre><code>output &lt;- lm(y ~ x, data = df)
  clcov &lt;- cl(output, df$groups)
  coeftest(output, clcov, nrow(df) - 1)
  </code></pre>
  
  <p>gives me the estimates I need. The problem now is that I want to use the model for prediction, and I need the standard error of the prediction to be calculated with the new covariance matrix <code>clcov</code>. That is, I need</p>
  
  <pre><code>predict(output, se.fit = TRUE)
  </code></pre>
  
  <p>but using <code>clcov</code> instead of <code>vcov(output)</code>. Something like a <code>vcov() &lt;-</code> would be perfect. </p>
  
  <p>Of course, I could write my own function to do predictions, but I am just wondering whether there is a more practical method that allows me to use methods for signature <code>lm</code> (like arm::sim).</p>