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

[Viechtbauer Wolfgang (STAT)]

Interesting, thanks for the feedback. I am replying to you and cc-ing r-sig-mixed-models, as this is getting a bit off-topic.

The degree of sparseness of course depends on the application. For phylogenies, both the matrix and its inverse are equally sparse (due to the split at the root node). For pedigrees, this may be different, but I don't know enough about this. The types of applications I have been working on are phylogenetic meta-analyses, so here using sparse matrices works just fine (and in fact, we have been doing comparisons between metafor and MCMCglmm and things match up pretty nicely).

At any rate, one can use rma.mv() just fine to fit the model that Matthew was asking about. Maybe it is not the most efficient way of doing this, but I don't see how it is wrong.

Best,
Wolfgang

> -----Original Message-----
> From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
> Sent: Thursday, February 26, 2015 18:18
> To: Viechtbauer Wolfgang (STAT)
> Cc: r-sig-mixed-models at r-project.org; Matthew Keller
> Subject: RE: [R-sig-ME] estimating variance components for arbitrarily
> defined var/covar matrices
> 
> Hi Wolfgang,
> 
> The matrix itself is not necessarily sparse (although usually is) its
> the inverse that is sparse.  The pattern of sparsity is also important
> because the equations can be permuted so they are more easily solved
> (if you use something like minimum degree ordering it will probably do
> a good job automatically). In addition, the nice structure is lost if
> individuals connecting two relatives are omitted. In this situation
> (the usual case) the random effect vector is usually augmented with
> these `missing' individuals and so the matrix is of greater dimension
> than the number of observations, and the corresponding columns of the
> design matrix are set to zero. I'm not sure if metafor handles these
> sorts of issues? The issues are not difficult to solve, but often
> general-purpose packages don't accommodate them because for most
> problems its not clear why you would ever need to worry about them.
> 
> Cheers,
> 
> Jarrod
> 
> Quoting "Viechtbauer Wolfgang (STAT)"
> <wolfgang.viechtbauer at maastrichtuniversity.nl> on Thu, 26 Feb 2015
> 17:50:36 +0100:
> 
> > Good point. Using sparse matrix methods is one way of handling that.
> > Set 'sparse=TRUE' for rma.mv() to do so.
> >
> > Best,
> > Wolfgang
> >
> >> -----Original Message-----
> >> From: Jarrod Hadfield [mailto:j.hadfield at ed.ac.uk]
> >> Sent: Thursday, February 26, 2015 17:42
> >> To: Viechtbauer Wolfgang (STAT)
> >> Cc: r-sig-mixed-models at r-project.org; Matthew Keller
> >> Subject: Re: [R-sig-ME] estimating variance components for arbitrarily
> >> defined var/covar matrices
> >>
> >> Hi,
> >>
> >> If the `genetic' matrix is pedigree derived (?) the matrix is a
> >> correlation matrix (in the absence of inbreeding) or proportional to
> >> the covariance matrix (with inbreeding). However, the matrix has a lot
> >> of structure that can be exploited by *not* treating it as just
> >> another correlation matrix. Specifically, the inverse of the matrix is
> >> all zeros, except for the diagonal elements and those that correspond
> >> to mates and/or parent-offspring. ASReml, WOMBAT, MCMCglmm, pedigreem
> >> (?) + others all exploit this structure leading to faster more robust
> >> results. For small problems this might not be an issue.
> >>
> >> Cheers,
> >>
> >> Jarrod
> >>
> >> Quoting "Viechtbauer Wolfgang (STAT)"
> >> <wolfgang.viechtbauer at maastrichtuniversity.nl> on Thu, 26 Feb 2015
> >> 17:24:39 +0100:
> >>
> >> > 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!