Residual Variance-Covariance matrix

Hi,

I am using the following the code to extract the residual
variance-covariance matrix cov(Yi/Xi):

I first fit the model with the name lmemod_lme4. I have an unbalanced
cluster dataset.

Then, I extracted components with the following:

var.d <- crossprod(getME(lmemod_lme4,"Lambdat"))
Zt <- getME(lmemod_lme4,"Zt")
vr <- sigma(lmemod_lme4)^2

Then, I combine them with the following:

var.b <- t(Zt) %*% var.d %*% Zt
sI <- vr * Diagonal(nrow(Nameofdataset))
var.y <- var.b + sI

I have 2799 observations in my dataset. MY var.y matrix has dimension
2799 times 2799. Is there now a way to extract the

cov(Yi/Xi) for each observation? Also, I get a very large value for
the determinant of var.y.

I would be grateful to get further guidance on this.

Thank you very much.

Regards