Reproducing results from previous answer is not working due to using new version of lme4

Question

I have tried to reproduce the results from the answers for this question “Estimating random effects and applying user defined correlation/covariance structure with R lme4 or nlme package “ https://stats.stackexchange.com/questions/18563/estimating-random-effects-and-applying-user-defined-correlation-covariance-struc

Aaron Rendahl's codes

 library(pedigreemm)
 relmatmm <- function (formula, data, family = NULL, REML = TRUE, relmat = list(), 
control = list(), start = NULL, verbose = FALSE, subset, 
weights, na.action, offset, contrasts = NULL, model = TRUE, 
x = TRUE, ...) 
{
mc <- match.call()
lmerc <- mc
lmerc[[1]] <- as.name("lmer")
lmerc$relmat <- NULL
if (!length(relmat)) 
    return(eval.parent(lmerc))
stopifnot(is.list(relmat), length(names(relmat)) == length(relmat))
lmerc$doFit <- FALSE
lmf <- eval(lmerc, parent.frame())
relfac <- relmat
relnms <- names(relmat)
stopifnot(all(relnms %in% names(lmf$FL$fl)))
asgn <- attr(lmf$FL$fl, "assign")
for (i in seq_along(relmat)) {
    tn <- which(match(relnms[i], names(lmf$FL$fl)) == asgn)
    if (length(tn) > 1) 
        stop("a relationship matrix must be associated with only one random effects term")
    Zt <- lmf$FL$trms[[tn]]$Zt
    relmat[[i]] <- Matrix(relmat[[i]][rownames(Zt), rownames(Zt)], 
        sparse = TRUE)
    relfac[[i]] <- chol(relmat[[i]])
    lmf$FL$trms[[tn]]$Zt <- lmf$FL$trms[[tn]]$A <- relfac[[i]] %*% Zt
}
ans <- do.call(if (!is.null(lmf$glmFit)) 
    lme4:::glmer_finalize
else lme4:::lmer_finalize, lmf)
ans <- new("pedigreemm", relfac = relfac, ans)
ans@call <- match.call()
ans
}

the original example

 set.seed(1234)
 mydata <- data.frame (gen = factor(rep(1:10, each = 10)),
                  repl = factor(rep(1:10, 10)),
                  yld = rnorm(10, 5, 0.5))
library(lme4)
covmat <- round(nearPD(matrix(runif(100, 0, 0.2), nrow = 10))$mat, 2)
diag(covmat) <- diag(covmat)/10+1
rownames(covmat) <- colnames(covmat) <- levels(mydata$gen)
m <- relmatmm(yld ~ (1|gen) + (1|repl), relmat=list(gen=covmat), data=mydata)

here is the error message

Error in lmf$FL : $ operator not defined for this S4 class
In addition: Warning message:
In checkArgs("lmer", doFit = FALSE) : extra argument(s) ‘doFit’ disregarded

I will appreciate any help ? Thanks

this is using internal structures from lme4 that have changed in the new version of lme4. I don't have time to look this over right now: you may need to ask the original author (Aaron Rendahl?) to update the code to work with the new version of lme4 (see help("modular",package="lme4") for a start. It might also be helpful to contact the maintainer of the pedigreemm package to see if they can help. — Ben Bolker
PS this should actually be easier to do with the new version of lme4 than it was before due to the improved modularity (see previous comment), but it still means the code needs to be redone ... — Ben Bolker

Ben Bolker Ben Bolker · Accepted Answer · 2013-10-15T13:12:07

This is a re-implementation of the previous code -- I have done some slight modifications, and I have not tested it in any way -- test yourself and/or use at your own risk.

First create a slightly more modularized function that constructs the deviance function and fits the model:

doFit <- function(lmod,lmm=TRUE) {
    ## see ?modular
    if (lmm) {
        devfun <- do.call(mkLmerDevfun, lmod)
        opt <- optimizeLmer(devfun)
        mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
    } else {
        devfun <- do.call(mkGlmerDevfun, lmod)
        opt <- optimizeGlmer(devfun)
        devfun <- updateGlmerDevfun(devfun, lmod$reTrms)
        opt <- optimizeGlmer(devfun, stage=2)
        mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)
    }
}

Now create a function to construct the object that doFit needs and modify it:

relmatmm <- function (formula, ..., lmm=TRUE, relmat = list()) {
    ff <- if (lmm) lFormula(formula, ...) else glFormula(formula, ...)
    stopifnot(is.list(relmat), length(names(relmat)) == length(relmat))
    relnms <- names(relmat)
    relfac <- relmat
    flist <- ff$reTrms[["flist"]]   ## list of factors
     ## random-effects design matrix components
    Ztlist <- ff$reTrms[["Ztlist"]]
    stopifnot(all(relnms %in% names(flist)))
    asgn <- attr(flist, "assign")
    for (i in seq_along(relmat)) {
        tn <- which(match(relnms[i], names(flist)) == asgn)
        if (length(tn) > 1) 
            stop("a relationship matrix must be",
                 " associated with only one random effects term")
        zn <- rownames(Ztlist[[i]])
        relmat[[i]] <- Matrix(relmat[[i]][zn,zn],sparse = TRUE)
        relfac[[i]] <- chol(relmat[[i]])
        Ztlist[[i]] <-  relfac[[i]] %*% Ztlist[[i]]
    }
    ff$reTrms[["Ztlist"]] <- Ztlist
    ff$reTrms[["Zt"]] <- do.call(rBind,Ztlist)
    fit <- doFit(ff,lmm)
}

Example

set.seed(1234)
mydata <- data.frame (gen = factor(rep(1:10, each = 10)),
                  repl = factor(rep(1:10, 10)),
                  yld = rnorm(10, 5, 0.5))

library(lme4)
covmat <- round(nearPD(matrix(runif(100, 0, 0.2), nrow = 10))$mat, 2)
diag(covmat) <- diag(covmat)/10+1
rownames(covmat) <- colnames(covmat) <- levels(mydata$gen)

m <- relmatmm(yld ~ (1|gen) + (1|repl), relmat=list(gen=covmat), 
     data=mydata)

This runs -- I don't know if the output is correct. It also doesn't make the resulting object into a pedigreemm object ...

Reproducing results from previous answer is not working due to using new version of lme4

the original example

here is the error message

1 Answers