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

Answer 1

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 ...