phylocov

phylocov.R

library(ape)
cat("(((Strix_aluco:4.2,Asio_otus:4.2):3.1,",
    "Athene_noctua:7.3):6.3,Tyto_alba:13.5);",
    file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
set.seed(101)
dd <- data.frame(y = rnorm(4),
                 obs = factor(1:4),
                 g = factor(1))
library(lme4)
## set up structures, ignoring the fact that we have as many RE as observations
## and that we have a dummy grouping variable
m1 <- lFormula(resp ~ 1 + (obs|g), dd,
               control = lmerControl(check.nobs.vs.nRE = "ignore",
                                     check.nobs.vs.nlev = "ignore",
                                     check.nlev.gtr.1 = "ignore"))
devfun <- do.call(mkLmerDevfun, m1)

## Cholesky factor of phylogenetic correlation matrix
phylocor <- chol(cov2cor(vcv(tree.owls)))
## need to extract *Lower* triangle in column-major order for
##  Lind to work right ...
phylovals <- t(phylocor)[lower.tri(phylocor, diag = TRUE)]

## practice inserting phylovals into Lambda
Lambdat <- m1$reTrms$Lambdat
Lambdat@x <- phylovals[m1$reTrms$Lind]
all.equal(unname(phylocor), unname(as.matrix(Lambdat)))

optfun <- function(r) {
    devfun(r*phylovals)
}
## this will fit residual variance + RE variance scaled relative to residual var
optim(optfun, par = 1, method = "Brent", lower = 0, upper = 1000)
## singular fit in this case