lme4qtl compared to lmekin?

[m66k]

I am really happy that I recently found your package lme4qtl. This seems to be exactly the package I have hoped for! So far, I have used the function lmekin from package coxme to model data with phylogenetic covariance matrix to account for species independence, but that function is not well developed (e.g not able to plot, predict). My only concern is that lme4qtl gives a bit different model parameters than lmekin. Most notably, the estimates and standard errors of the variables directly connected to species identity, such as mean body mass of species, are 2-3 times higher than in lmekin.

For example, it seems a bit odd that the model output from lme4qtl indicates a significant interaction between mass and habitat but the confidence bands for mass are so wide that this makes the interaction hard to believe.
Do you know what could be the reason for the large difference between lme4qtl and lmekin? Is one of them wrong or which package should I prefer for using with my type of data?

I will send my data sample, phylogenetic tree, and code by e-mail.

[variani]

Hi,

I added scripts related to your problem in this commit. Basically, I just considered a single random effect, and it seems the outputs are now more consistent.

I am not an experienced user of neither pylogenetic data or visreg. I also tried lmekin a little. My initial guess is that you've might misused lmekin (its argument varlist) when specifying two random effects. It happened to me in my experiments:

lmeqtl_2 <- relmatLmer(APTT ~ AGE + (1|ID) + (1|HHID), dat, relmat = list(id = dkin)),
lmekin_2 <- lmekin(APTT ~ AGE + (1|ID) + (1|HHID), dat, varlist=list(coxmeMlist(dkin), coxmeFull))

I don't remember "why", by the syntax with coxmeFull solved my problem that time.

Here results for your analysis go:

library(coxme)
m_lmekin <- lmekin(variable ~ mass*habitat + (1|species), data=dataX,
                      varlist=phylo_all, method = "REML")
print(m_lmekin)

#  Log-likelihood = -111.105 
#  n= 1868 
#
# Model:  variable ~ mass * habitat + (1 | species) 
#
# Fixed coefficients
#                   Value  Std Error      z       p
# (Intercept)    1.0926508 0.22466203   4.86 1.2e-06
# mass           0.1483706 0.05843213   2.54 1.1e-02
# habitatU      -0.2663249 0.01286789 -20.70 0.0e+00
# mass:habitatU -0.1065064 0.02695353  -3.95 7.8e-05

library(lme4qtl)
m_lme4qtl <- relmatLmer(variable ~ mass*habitat + (1|species), data=dataX, 
                     relmat = list(species = phylo_all), REML = T)

summary(m_lme4qtl)[["logLik"]]

# 'log Lik.' -107.7412 (df=6)

summary(m_lme4qtl)[["coefficients"]]

#               Estimate Std. Error    t value
# (Intercept)    0.9814953 0.23427874   4.189434
# mass           0.3642359 0.15830427   2.300860
# habitatU      -0.2655202 0.01284164 -20.676496
# mass:habitatU -0.1103905 0.02686246  -4.109473

The estimates are more consistent, and the differences are due to the fact that lme4qtl did slightly a better job based on logLik numbers.

It would be nice if you can explore further on the correct syntax of using lmekin for two-random-effects models. Is the agreement between lme4qtl and lmekin for one-random-effect model OK for you?

You also might be interested in help("pvalues", package = "lme4").

Best,
Andrey

[variani]

Hi,

It seems I fixed the issue, which actually was a part of lmekin functionality. One needs to take a special care of the order in two variables: id grouping variable (species) and rownames/colnames of the covariance matrix (phylo_all). Note that lme4qtl is robust to this ordering issue (I implemented this kind of control someday, when I faced with a similar inconsistent models' behavior).

The code needed to make lmekin working properly:

species_vals <- levels(dataX$species)
phylo_all_ordered <- phylo_all[species_vals, species_vals]

You can also take a look at the R file I created in the last commit 02-check-ids.R.

I guess your original example should work now.

Best,
Andrey

[m66k]

Hi!

Oh, that's interesting. Thank you so much for helping me figure out what caused the difference!

All the best!