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

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.

• Aaron, thank you for your thoughts, hope will get more robust suggestion on this ... Commented Nov 18, 2011 at 10:46
• The example is extremely confusing because it strongly suggests a different kind of dataset altogether; it contradicts the question. Please either delete this example or provide a realistic one.
– whuber
Commented Nov 18, 2011 at 17:21
• @whuber I edited some of my typo and made my point clearer, hope helps Commented Nov 18, 2011 at 19:05
• @RamSharma: I took the liberty to make a sample positive definite covariance matrix, and made a few minor edits; feel free to edit back if I've changed something important. Commented Nov 23, 2011 at 3:43
• I think we should migrate this to stackoverflow, to get more views. I do not how to do it, can somebody help ?
– John
Commented Nov 23, 2011 at 12:04

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.

• How about: lme(yld ~ 1, data = mydata, random = ~ gen + repl, correlation = covmat)# the formula is giving and error and I am not sure that if the correlation structure applies to replication or residual term, what do you think ?
– John
Commented Nov 18, 2011 at 11:17
• @John: Crossed random effects are tricky with nlme and something more complicated is needed that gen + repl; also, I think the correlation needs to call one of nlme's covariance/correlations functions with covmat as a parameter. The notation for this is really tricky so to say more I'd need Pinhiero/Bates on hand, which I don't today. Commented Nov 18, 2011 at 15:38
• then if there is no repl effect, do you think lme(yld ~ 1, data = mydata, random = ~ gen, correlation = covmat) is appropriate or equivalent to asreml code: asreml(yld ~ 1, random = ~ gen, ginverse = list (gen = inv.covmat), data = mydata), where inv.covmat is lower triangle melted matrix (see asreml-r documentation)
– John
Commented Nov 20, 2011 at 1:10
• The right notation would be probably be something like lme(yld ~ 1, data = mydata, random = ~ 1 | gen, correlation = corSymm(value=covmatX, form= ~ gen, fixed=TRUE)), where covmatX is a modified version of covmat to make it however corSymm wants it. I'm not quite sure the form term is right either. Commented Nov 23, 2011 at 3:41
• I had a link to my version of lme4 here but removed it as it's not well documented and I decided I'd rather not have it too visible, especially as I think the kinship package will meet most needs. If anyone wants my version of lme4 please contact me by email; I shouldn't be too hard to find. Commented Nov 23, 2011 at 20:15

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 !

• Suggested by RamSharma: Solution using kinship: require(kinship); model1 <- lmekin(yld ~1, random = ~ 1|gen, varlist =list(covmat), data=mydata).
– chl
Commented Nov 24, 2011 at 12:49
• further edits to my suggestion: model1 <- lmekin(yld ~1, random = ~ 1|gen, varlist =list(covmat), data=mydata), still there is a problem, sorry for premature posting. Can somebody fix it? Commented Nov 24, 2011 at 12:52

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

• yes, fitting crossed random effect is not issue alone however fitting user defined co-variance structure is main issue and the question primarily seeks answer for that.
– John
Commented Nov 23, 2011 at 11:15

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


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