I have 10 different data sets in R and for the coefficients I have combined the estimates using Rubin's rules. How would I combine the variance covariance matrix? Should I find the mean for each value in my 10 matrices or be extracting the actual random effects estimates and be combining them and finding the overall variances? Likewise for the covariance values?
I think the guts of this package give you the answer. That is, provided that the authors of that package (and the referenced paper, Nassiri et al 2018) know what they're doing, this code applies Rubin's rules at the level of the covariance matrix:
# estimate covariance matrix of imputed values and combine estimates
cov_imps <- lapply(data_imps, stats::cov)
cov_comb <- Reduce("+", cov_imps) / length(cov_imps)
In other words, they're computing the mean of the full covariance matrix (this is equivalent to taking elementwise means). I would not compute the covariance matrix of the combined BLUPs/conditional modes (for reasons).
Technical note: Since a linear combination of positive-definite matrices with positive coefficients is positive definite (the linked answer refers to positive semi-definiteness, but I think the argument should extend to positive definiteness), the average of covariance matrices should also be a sensible covariance matrix.
Nassiri, V., Lovik, A., Molenberghs, G., & Verbeke, G. (2018). On using multiple imputation for exploratory factor analysis of incomplete data. Behavioral Research Methods 50, 501–517. doi: 10.3758/s13428-017-1013-4