Solved – Estimating random effects and applying user defined correlation/covariance structure with R lme4 or nlme package

mixed modelr

I have following type of data. I have evaluated 10 individuals each repeated 10 times. I have 10×10 relation matrix (relationship between all combination of the individuals).

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

This gen is different varieties of plant, so each can be repeatedly grown and yield is measured. The covariance matrix is relatedness measure by genetic similarity calculated by ibd probabilities in seperate experiments. 

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)

> covmat                   
10 x 10 Matrix of class "dgeMatrix"                    
      1    2    3    4    5    6    7    8    9   10
1  1.00 0.08 0.06 0.03 0.09 0.09 0.10 0.08 0.07 0.10
2  0.08 1.00 0.08 0.09 0.04 0.12 0.08 0.08 0.11 0.09
3  0.06 0.08 1.00 0.10 0.05 0.09 0.09 0.07 0.04 0.13
4  0.03 0.09 0.10 1.00 0.02 0.11 0.09 0.06 0.04 0.12
5  0.09 0.04 0.05 0.02 1.00 0.06 0.07 0.05 0.02 0.08
6  0.09 0.12 0.09 0.11 0.06 1.00 0.12 0.08 0.07 0.14
7  0.10 0.08 0.09 0.09 0.07 0.12 1.00 0.08 0.03 0.15
8  0.08 0.08 0.07 0.06 0.05 0.08 0.08 1.00 0.06 0.09
9  0.07 0.11 0.04 0.04 0.02 0.07 0.03 0.06 1.00 0.03
10 0.10 0.09 0.13 0.12 0.08 0.14 0.15 0.09 0.03 1.00

My model is: 

yld = gen + repl + error

both gen and repl are considered random and I want to get the random effect estimates associated with each gen, however I need to consider the relationship matrix.

If it is too complex to fit nested models, I would just remove repl from the model, but ideally I will keep it. 

yld = gen +  error

How can I achieve this using R packages, perhaps with nlme or lme4? I know that ASREML can do it but I do not have hold and I love R for being robust as well as free.

Best Answer

Try the kinship package, which is based on nlme. See this thread on r-sig-mixed-models for details. I'd forgotten about this as I was trying to do it for a logistic model. See https://stackoverflow.com/questions/8245132 for a worked-out example.

For non-normal responses, you'd need to modify the pedigreemm package, which is based on lme4. It gets you close, but the relationship matrix has to be created from a pedigree. The below function is a modification of the pedigreemm function which takes an arbitrary relationship matrix instead.

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
}

Usage is similar to pedigreemm except you give it the relationship matrix as the relmat argument instead of the pedigree as the pedigree argument.

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

This doesn't apply here as you have ten observations/individual, but for one observation/individual you need one more line in this function and a minor patch to lme4 to allow for only one observation per random effect.

Related Question