Estimación de efectos aleatorios y aplicación de estructura de correlación / covarianza definida por el usuario con el paquete R lme4 o nlme

Tengo el siguiente tipo de datos. He evaluado a 10 personas cada una repetidas 10 veces. Tengo una matriz de relación de 10x10 (relación entre todas las combinaciones de los individuos).

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

Este gen es de diferentes variedades de plantas, por lo que cada una puede cultivarse repetidamente y se mide el rendimiento. La matriz de covarianza es una medida de parentesco por similitud genética calculada por probabilidades ibd en experimentos separados.

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

Mi modelo es:

yld = gen + repl + error 

tanto gen como repl se consideran aleatorios y quiero obtener las estimaciones de efectos aleatorios asociados con cada gen, sin embargo, debo considerar la matriz de relación.

Si es demasiado complejo para adaptarse a modelos anidados, simplemente eliminaría las respuestas del modelo, pero idealmente lo mantendré.

yld = gen +  error 

¿Cómo puedo lograr esto usando paquetes R, quizás con nlme o lme4? Sé que ASREML puede hacerlo, pero no tengo agarre y me encanta R por ser robusto y libre.

Pruebe el kinshippaquete, que se basa en nlme. Vea este hilo en modelos mixtos r-sig para más detalles. Me había olvidado de esto cuando intentaba hacerlo para un modelo logístico. Consulte /programming/8245132 para ver un ejemplo resuelto.

Para respuestas no normales, necesitaría modificar el paquete pedigreemm, que se basa en lme4. Te acerca, pero la matriz de relación debe crearse a partir de un pedigrí. La siguiente función es una modificación de la pedigreemmfunción que toma una matriz de relación arbitraria en su lugar.

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
}

El uso es similar a, pedigreemmexcepto que usted le da la matriz de relación como relmatargumento en lugar del pedigrí como pedigreeargumento.

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

Esto no se aplica aquí, ya que tiene diez observaciones / individuo, pero para una observación / individuo necesita una línea más en esta función y un parche menor lme4para permitir solo una observación por efecto aleatorio.

Esta respuesta es una expansión potencial de la sugerencia hecha por Aaron, quien sugirió usar Pedigreem. El pedigreem puede calcular la relación de los proyectos como la siguiente sintaxis, no sé cómo podemos usar esa salida de relación de manera diferente.

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

El ajuste del modelo mixto del paquete se basa en lme4 para que la sintaxis de la función principal sea similar a la función lme4 del paquete lmer, excepto que puede colocar el objeto genealógico en él.

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

Sé que esta no es la respuesta perfecta a su pregunta, sin embargo, esto puede ayudar un poco. ¡Me alegra que hayas hecho esta pregunta, interesante para mí!

lmer()en el lme4paquete permite efectos aleatorios cruzados. Aquí, usarías algo como

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

Para una referencia completa;

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