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

Question

I have following type of data. I have evaluated 10 individuals each repeated 10 times. I have 10x10 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.

Answer 1

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.

Answer 2

This answer is potential expansion of the suggestion made by Aaron, who has suggested to use Pedigreem. The pedigreem can compute relationship from the projects as following syntax, I am unaware how we can use such relation output from different way.

# just example from the manual to create pedigree structure and relation matrix 
  # (although you have already the matrix in place) 
p1 <- new("pedigree",
sire = as.integer(c(NA,NA,1, 1,4,5)),
dam = as.integer(c(NA,NA,2,NA,3,2)),
label = as.character(1:6))
p1
(dtc <- as(p1, "sparseMatrix")) # T-inverse in Mrode’s notation
solve(dtc)
inbreeding(p1) 

The mixed model fit of the package is based on lme4 for the syntax for the main function is similar to lme4 package function lmer function except you can put the pedigree object in it.

pedigreemm(formula, data, family = NULL, REML = TRUE, pedigree = list(),
 control = list(),
start = NULL, verbose = FALSE, subset, weights, na.action, 
  offset, contrasts = NULL, model = TRUE, x = TRUE, ...)

I know this is not perfect answer to your question, however this can help a little bit. i am glad you asked this question, interesting to me !

Answer 3

lmer() in the lme4 package permits crossed random effects. Here, you'd use something like

y ~ (1|gen) + (1|repl)

For a full reference;

http://www.stat.wisc.edu/~bates/PotsdamGLMM/LMMD.pdf

Answer 4

Your title says "with lme4 or nlme package", but your text says

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.

This approach is not based on these two packages, but it is open source and very flexible. GBLUP with arbitrary covariance structure is a special case of RKHS regression aka Kernel Ridge Regression. The package BGLR estimates the variance components in a Bayesian Framework. An alternative is the package KRMM that seems to solve the same model but using Expectation Maximization instead of a Bayesian approach (Gibbs sampling). But I didn't test that.

This excerpt from the BGLR extended documentation computes

y ~ a + g + e

where a is a random effect with a pedigree-derived covariance structure, g is a random effect using a marker-derived covariance structure (you can use another genetic distance instead of the definition shown here) and e is the residual. For your problem, you can of course just omit a (=list(K=A, ...). The genomic relationship matrices (G and A in this example) must relate 1-to-1 to the genotype order in y, so if a genotype occurs multiple times in y, it must do so in the matrices as well.

Box 4a: Fitting a Pedigree + Markers regression using Gaussian Processes

#1# Loading and preparing the input data
library(BGLR);
data(wheat);Y<-wheat.Y; X<-wheat.X; A<-wheat.A;
y<-Y[,1]

#2# Computing the genomic relationship matrix
X<-scale(X,center=TRUE,scale=TRUE)
G<-tcrossprod(X)/ncol(X)

#3# Computing the eigen-value decomposition of G
EVD <-eigen(G)

#3# Setting the linear predictor
ETA<-list(list(K=A, model='RKHS'),
          list(V=EVD$vectors,d=EVD$values, model='RKHS'))

#4# Fitting the model
fm<-BGLR(y=y,ETA=ETA, nIter=12000, burnIn=2000,saveAt='PGBLUP_') 
save(fm,file='fmPG_BLUP.rda')

See also these examples of different ways to compute GBLUP.

This documentation page shows an example including fixed effects (and other methods, such as BayesB, just use those models you need):

pheno=mice.pheno

fm=BGLR(y=pheno$Obesity.BMI,
        ETA=list(
          fixed=list(~factor(GENDER)+factor(Litter),data=pheno,model='FIXED'),
          cage=list(~factor(cage),data=pheno,model='BRR'),
          ped=list(K=A,model='RKHS'),
          mrk=list(X=X,model='BayesB')
     )
)