Estimation des effets aléatoires et application d'une structure de corrélation / covariance définie par l'utilisateur avec le package R lme4 ou nlme

J'ai le type de données suivant. J'ai évalué 10 individus chacun répété 10 fois. J'ai une matrice de relations 10x10 (relation entre toutes les combinaisons des individus).

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

Cette génération est constituée de différentes variétés de plantes, chacune peut donc être cultivée à plusieurs reprises et le rendement est mesuré. La matrice de covariance est une mesure de parenté par similitude génétique calculée par les probabilités ibd dans des expériences séparées.

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

Mon modèle est:

yld = gen + repl + error 

gen et repl sont considérés comme aléatoires et je veux obtenir les estimations de l'effet aléatoire associées à chaque génération, mais je dois prendre en compte la matrice de relation.

S'il est trop complexe pour s'adapter à des modèles imbriqués, je supprimerais simplement repl du modèle, mais idéalement je le garderai.

yld = gen +  error 

Comment puis-je y parvenir en utilisant des packages R, peut-être avec nlme ou lme4? Je sais qu'ASREML peut le faire mais je n'ai pas de prise et j'aime R pour être robuste aussi bien que gratuit.

Essayez le kinshippackage, qui est basé sur nlme. Voir ce fil sur les modèles r-sig-mixed-pour plus de détails. J'avais oublié cela alors que j'essayais de le faire pour un modèle logistique. Voir /programming/8245132 pour un exemple élaboré.

Pour les réponses non normales, vous devez modifier le package pedigreemm, qui est basé sur lme4. Cela vous rapproche, mais la matrice des relations doit être créée à partir d'un pedigree. La fonction ci-dessous est une modification de la pedigreemmfonction qui prend à la place une matrice de relation arbitraire.

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
}

L'utilisation est similaire, pedigreemmsauf que vous lui donnez la matrice de relation comme relmatargument au lieu du pedigree comme pedigreeargument.

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

Cela ne s'applique pas ici car vous avez dix observations / individu, mais pour une observation / individu, vous avez besoin d'une ligne de plus dans cette fonction et d'un patch mineur lme4pour ne permettre qu'une seule observation par effet aléatoire.

Cette réponse est une extension potentielle de la suggestion faite par Aaron, qui a suggéré d'utiliser Pedigreem. Le pedigreem peut calculer la relation à partir des projets comme la syntaxe suivante, je ne sais pas comment nous pouvons utiliser une telle sortie de relation de manière différente.

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

L'ajustement de modèle mixte du package est basé sur lme4 car la syntaxe de la fonction principale est similaire à la fonction de package lme4, fonction lmer, sauf que vous pouvez y insérer l'objet pedigree.

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

Je sais que ce n'est pas la réponse parfaite à votre question, mais cela peut aider un peu. je suis content que vous ayez posé cette question qui m'intéresse!

lmer()dans l' lme4emballage permet des effets aléatoires croisés. Ici, vous utiliseriez quelque chose comme

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

Pour une référence complète;

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