title: Get Residual Variance-Covariance Matrix in lme4
body: <p>I am fitting a linear mixed effects model using <code>lme4</code>:</p>
  
  <pre><code>library(lme4)
  data(Orthodont)
  dent &lt;- Orthodont
  d.test &lt;- lmer(distance ~ age + (1|Subject), data=dent)
  </code></pre>
  
  <p>If we say generically <code>Y = X * B + Z * d + e</code> is the form of a linear mixed effects model, then I am trying to get <code>Var(Y) = Z * Var(d) * Z^t + Var(e)</code> from the results of the model.</p>
  
  <p>Is the following formulation the right way to do this?</p>
  
  <pre><code>k &lt;- table(dent$Subject)[1]
  vars &lt;- VarCorr(d.test)
  v &lt;- as.data.frame(vars)
  sigma &lt;- attr(vars, "sc")
  s.tech &lt;- diag(v$vcov[1], nrow=k)
  icc &lt;- v$vcov[1]/sum(v$vcov)
  s.tech[upper.tri(s.tech)] &lt;- icc
  s.tech[lower.tri(s.tech)] &lt;- icc
  sI &lt;- diag(sigma^2, nrow=length(dent$age))
  var.b &lt;- kronecker(diag(1, nrow=length(dent$age)/k), s.tech)
  var.y &lt;- sI + var.b
  </code></pre>
  
  <p>I think this is a simple question, but I can't find anywhere code for doing this, so I'm asking if I'm doing it right.</p>