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

Question

Que dispone de los siguientes tipos de datos. He evaluado de 10 individuos cada uno repite 10 veces. He 10x10 relación de la matriz de relación entre todos los combinación 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 la planta, por lo que cada uno puede ser repetidamente crecido y la rentabilidad se mide. La matriz de covarianza es la relación, medida por la similitud genética calculada por la eii probabilidades 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 

ambos gen y repl se consideran aleatorios y quiero conseguir el efecto aleatorio de las estimaciones asociadas con cada generación, sin embargo, tengo la necesidad de considerar la relación de la matriz.

Si es demasiado complejo como para adaptarse a modelos anidados, me acaba de quitar repl de la modelo, pero lo ideal será mantenerlo.

yld = gen +  error 

¿Cómo puedo lograr esto mediante R paquetes, tal vez con nlme o lme4? Sé que ASREML puede hacerlo, pero no me han presionado y me encanta R por ser robusto así como libre.

Answer 1

Pruebe el kinship paquete, que se basa en nlme. Ver este hilo en r-sig-mixto-modelos para más detalles. Me había olvidado de este tema como yo estaba tratando de hacer un modelo logístico. Ver http://stackoverflow.com/questions/8245132 para un resueltos de ejemplo.

Para los no-respuestas normales, sería necesario modificar el pedigreemm paquete, que se basa en lme4. Se pone de cerrar, pero la relación de la matriz tiene que ser creado a partir de un árbol genealógico. La siguiente función es una modificación de la pedigreemm función que toma una relación arbitraria de la matriz 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 pedigreemm , salvo que le den la relación de la matriz como de la relmat argumento en lugar de la pedigrí como pedigree argumento.

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

Esto no se aplica aquí tienes diez observaciones/persona, pero para una observación o la persona de que se necesita una línea más en esta función y un menor de edad parche lme4 para permitir sólo una observación por efecto aleatorio.

Answer 2

Esta respuesta es la expansión del potencial de la sugerencia hecha por Aarón, que ha sugerido el uso de Pedigreem. El pedigreem puede calcular la relación de los proyectos de la siguiente sintaxis, no me doy cuenta de cómo podemos usar este tipo de relación de salida 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 modelo mixto de ajuste del paquete se basa en lme4 para la sintaxis de la función principal es similar al paquete lme4 función de lmer función, excepto que usted puede poner el árbol genealógico de los objetos en ella.

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 alegro de que haya hecho esta pregunta, interesante para mí !

Answer 3

lmer() en la lme4 paquete permite cruzado de efectos aleatorios. Aquí, tendría que usar algo como

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

Para una referencia completa;

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