# 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 ... – Ram Sharma Nov 18 '11 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 Nov 18 '11 at 17:21
@whuber I edited some of my typo and made my point clearer, hope helps – Ram Sharma Nov 18 '11 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. – Aaron Nov 23 '11 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 Nov 23 '11 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 http://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 Nov 18 '11 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. – Aaron Nov 18 '11 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 Nov 20 '11 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. – Aaron Nov 23 '11 at 3:41
@Aaron, do have this patch handy? I would need it often to fit a model with single individual for each class...I might want to ask as seperate question ....it might be too much in this question – Ram Sharma Nov 23 '11 at 14:51

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 Nov 24 '11 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? – Ram Sharma Nov 24 '11 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 Nov 23 '11 at 11:15