Variance explained by similarity matrix

Question

Given a vector $Y$ for $N$ subjects and a matrix $X$ of $d$ variables for each of our $N$ subjects, we can compute matrix $M$, a $N$x$N$ random effects matrix for how similar each of our $N$ subjects are to each other, across our $d$ variables in $X$. For instance, via a correlation. I can then run a mixed effects model regression where $Y$ is my DV and where $M$ describes the similarity between the $N$ participants. My question is: How much variance in $Y$ is explained by the random effects matrix $M$? Specifically, I am looking for a worked example in R that does not use pedigrees, yielding both an estimate of the variance explained, as well as some estimate of the error.

There have been some prior related questions (one, two, three). These examples are either not fully worked, use pedigrees, or don't give an error for the estimate (perhaps this last bit may be trivial, but it would be helpful if someone could walk me through it).

For context, I'm trying to build off the approach in this paper, but I'm having difficulty interpreting their matlab code. Another paper which applies the approach I'm describing is this one ("heritability" is the variance explained by the genetic similarity between individuals). I am essentially trying to apply this approach in a general manner.

So far, based on prior questions, I've worked out something using coxme::lmekin(), but I'm not sure it's correct, particularly whether the similarity matrix I'm giving the function is being handled the way I expect, as the function was written for a pedigree. I'm probably also not calculating the variance explained by $M$ correctly at the end, but I'm having some difficulty interpreting the coxme::lmekin() documentation.

> library(clusterGeneration)
> library(MASS)
> library(Matrix)
> library(coxme)
> 
> d=201
> N = 100
> set.seed(1234)
> 
> # simulate covariance matrix
> cov_mat<-clusterGeneration::genPositiveDefMat(dim =d,covMethod = "onion")
> 
> # simulate data
> dat<-MASS::mvrnorm(N,mu = rep(0,d),Sigma = cov_mat$Sigma)
> Y = scale(dat[,1],center = T,scale = T)
> X = dat[,-1]
> 
> # subj similarity, based on data, excluding Y
> M = cor(t(X))
> # check if PD
> if(min(eigen(M)$values)<0){M=Matrix::nearPD(x = M)$mat}
> 
> # rescale correlation to a similarity matrix
> range01 <- function(x){(x-min(x))/(max(x)-min(x))}
> M.2  =range01(M)
> 
> # each subj gets their own random intercept
> rand = c(1:length(Y))
> 
> # without similarity matrix
> a1=lmekin(Y ~ 1 + (1|rand))
> a1
Linear mixed-effects kinship model fit by maximum likelihood
  Data: NULL 
  Log-likelihood = -141.3913 
  n= 100 


Model:  Y ~ 1 + (1 | rand) 
Fixed coefficients
                   Value  Std Error z p
(Intercept) 8.926083e-17 0.09949874 0 1

Random effects
 Group Variable  Std Dev    Variance  
 rand  Intercept 0.09900495 0.00980198
Residual error= 0.9900495 
> # variance explained by random effects is nearly 0
> sum(a1$vcoef$rand)/(var(a1$residuals) + sum(a1$var)+sum(a1$vcoef$rand))
[1] 0.009802
> 
> # with similarity matrix
> a2 =lmekin(Y ~ 1 + (1|rand),varlist = list(M.2))
> a2
Linear mixed-effects kinship model fit by maximum likelihood
  Data: NULL 
  Log-likelihood = -130.4529 
  n= 100 


Model:  Y ~ 1 + (1 | rand) 
Fixed coefficients
                Value Std Error    z    p
(Intercept) 0.1686349 0.4971669 0.34 0.73

Random effects
 Group Variable Std Dev   Variance 
 rand  Vmat.1   0.9976159 0.9952375
Residual error= 0.5060629 
> # variance explained by random effects is much larger
> sum(a2$vcoef$rand)/(var(a2$residuals) + sum(a2$var)+sum(a2$vcoef$rand))
[1] 0.7369908