1

我试图从这个问题的答案中重现结果“Estimating random effects and apply user defined correlation/covariance structure with R lme4 or nlme package” https://stats.stackexchange.com/questions/18563/estimating-random-效果和应用用户定义的相关协方差结构

亚伦伦达尔的代码 
 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
}

原始示例

 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)

这是错误消息

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

我将不胜感激任何帮助？谢谢

4

1 回答 1

2

这是对先前代码的重新实现——我做了一些轻微的修改，并且我没有以任何方式对其进行测试——测试你自己和/或使用你自己承担风险。

首先创建一个稍微模块化的函数来构造偏差函数并拟合模型：

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)
    }
}

现在创建一个函数来构造doFit需要的对象并对其进行修改：

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)
}

例子

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)

这运行 - 我不知道输出是否正确。它也不会将生成的对象变成pedigreemm对象...

于 2013-10-15T13:12:07.383 回答