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

[Anthony Ives]

Steve,

Just a quick follow-up. Even though your work (and communityPGLMM) does not allow transformations in the covariance matrices themselves (like the OU), it does in effect perform Pagel’s lambda in the form sigma2*V + I ~= lambda*V + (1-lambda)*I.

Cheers, Tony


Anthony Ives
Department of Zoology
459 Birge Hall (4th floor, E end of bldg)
UW-Madison
Madison, WI 53706
608-262-1519

> On Feb 26, 2015, at 1:09 PM, Steve Walker <steve.walker at utoronto.ca> wrote:
> 
> Thanks Tony.  For those interested, you can check out the development of my work on phylogenetic models in lme4 here:
> 
> https://github.com/stevencarlislewalker/lme4ord
> 
> The section in the README file about phylogenetic models gives an example and briefly describes the kinds of models that can be fitted.
> 
> Regarding speed, the glmerc (with a 'c' for known covariance over grouping factor levels) function can fit models similar to the one in the README but with a 500 tip phylogeny in about 20 seconds on my pretty standard macbook pro.  lme4 uses sparse matrices for all random effect-related matrix computations, so most of the speed comes from this infrastructure.  However, it may be possible to speed things up more by exploiting Kronecker-product-type structure, but this would require some substantial additions to lme4, I think.
> 
> I should also note that the models fitted by glmerc do _not_ allow for fitting parameters that scale branch lengths (e.g. Ornstein-Uhlenbeck).  The fundamental challenge with these models in lme4 is that lme4 is based on a Cholesky factor parameterization of the random effects covariance matrix.  To parameterize the covariance matrix on the covariance or branch length scale, one would need to compute a Cholesky factor every time the deviance function is evaluated.  This could get costly.  The recent methods of Ho and Ané (2014) in the phylolm package might help here.  However, in their current form these methods only apply to quadratic terms involving the inverse covariance matrix and to the log determinant of the covariance matrix.  These computations are certainly important in mixed modelling in general, but unfortunately not _so_ much for lme4, because of the Cholesky parameterization.  However, I bet that their ideas could be modified to be applicable to lme4.  In particular I bet one could compute a series of small Cholesky decompositions at the tips and then iteratively rank-one update them by traversing the tree to the root.  My guess is that this procedure would scale linearly, as do the other methods of Ho and Ané (2014).
> 
> Cheers,
> Steve
> 
> On 2015-02-26 7:04 AM, Anthony R Ives wrote:
>> 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!
>>>> 
>>>> [[alternative HTML version deleted]]
>>>> 
>>>> _______________________________________________
>>>> R-sig-mixed-models at r-project.org mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>> 
>>>> 
>>> 
>>> 
>>> 
>>> -- The University of Edinburgh is a charitable body, registered in
>>> Scotland, with registration number SC005336.
>>> 
>>> _______________________________________________
>>> R-sig-mixed-models at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> 
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> 
>> 
>