Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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.

share|improve this question
    
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
add comment

3 Answers

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.

share|improve this answer
    
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
show 3 more comments

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 !

share|improve this answer
    
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
add comment

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

share|improve this answer
1  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.