[R-sig-ME] estimating variance components for arbitrarily defined var/covar matrices

[Viechtbauer Wolfgang (STAT)]

And just to throw another package in the mix, you can also use the metafor package for this. While this doesn't quite sound like a meta-analysis, in the end, meta-analysis models are just mixed-effects models. And in a phylogenetic meta-analysis, we add random effects with known correlations matrices to the model.

The syntax would be:

id.e <- 1:nrow(dat)
id.r <- 1:nrow(dat)
rma.mv(y ~ <fixed effects>, V = 0, random = list(~ 1 | id.r, ~ 1 | id.e), R = list(id.r = GRM), data=dat)

where 'GRM' is the n x n matrix of similarities. The V = 0 part seems a bit strange, but in meta-analytic models, we usually don't estimate the error variance and instead have known sampling variances (or even a known variance-covariance matrix of the sampling errors). Here, we don't, so we just set that part to 0 (you'll get a warning that 'V appears to be not positive definite.' but you can safely ignore this). This will fit the model that you specified below.

Besides the estimates of the fixed effects, the results will include two variance components, one for id.r (this is VG) and one for id.e (this is sigma^2_error). The default is REML estimation (method="ML" is you want MLEs).

By default, rma.mv() tries to rescale the matrix supplied via the R argument into a correlation matrix. Not sure if your matrix of similarities is really a correlation matrix or not. In case you don't want the function to mess with the matrix, set 'Rscale=FALSE' and it won't touch it.

Best,
Wolfgang

--   
Wolfgang Viechtbauer, Ph.D., Statistician   
Department of Psychiatry and Psychology   
School for Mental Health and Neuroscience   
Faculty of Health, Medicine, and Life Sciences   
Maastricht University, P.O. Box 616 (VIJV1)   
6200 MD Maastricht, The Netherlands   
+31 (43) 388-4170 | http://www.wvbauer.com   

> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-
> project.org] On Behalf Of Anthony R Ives
> Sent: Thursday, February 26, 2015 13:04
> To: Jarrod Hadfield; Matthew Keller
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] estimating variance components for arbitrarily
> defined var/covar matrices
> 
> Matthew,
> 
> You should be able to do this in communityPGLMM in {pez}. Also, Steve
> Walker is currently working on a way to do this in lmer/glmer.
> 
> Cheers, Tony
> 
> 
> On 02/26/15, Jarrod Hadfield  wrote:
> > Hi Matthew,
> >
> > Both MCMCglmm and asreml-r fit these models in R.
> >
> > Cheers,
> >
> > Jarrod
> >
> > Quoting Matthew Keller <mckellercran at gmail.com> on Wed, 25 Feb 2015
> 16:42:32 -0700:
> >
> > >Hi all,
> > >
> > >This is a typical problem in genetics and I'm trying to figure out
> whether
> > >there's any way to solve it using lmer or similar, and if not, why it
> isn't
> > >possible.
> > >
> > >Often in genetics, we have an n-by-n matrix (n=sample size) of genetic
> > >relationships, where the diagonal is how related you are to yourself
> (~1,
> > >depending on inbreeding) and off-diagonals each pairwise relationship.
> I'd
> > >like to be able to use lmer or some other function in R to estimate
> the
> > >variance attributable to this genetic relationship matrix. Thus:
> > >y = b0 + b*X + g*Z + error
> > >where y is a vector of observations, b is a vector of fixed covariate
> > >effects and g is a vector of random genetic effects. X and Z are
> incidence
> > >matrices for b & g respectively, and we assume g ~ N(0, VG). The
> variance
> > >of y is therefore
> > >var(y) = Z*Z' * VG + I*var(e)
> > >
> > >Z*Z' is the observed n-by-n genetic relationship matrix. Given an
> observed
> > >Z*Z' genetic relationship matrix, is there a way to estimate VG?
> > >
> > >I guess this boils down to, if we have an observed n-by-n matrix of
> > >similarities, can we use mixed models in R to get the variance in y
> that is
> > >explained by that similarity?
> > >
> > >Thanks in advance!