Schätzen zufälliger Effekte und Anwenden einer benutzerdefinierten Korrelations- / Kovarianzstruktur mit dem Paket R lme4 oder nlme

Ich habe folgende Art von Daten. Ich habe 10 Personen bewertet, die jeweils 10 Mal wiederholt wurden. Ich habe eine 10x10 Beziehungsmatrix (Beziehung zwischen allen Kombinationen der Individuen).

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

Bei diesem Gen handelt es sich um verschiedene Pflanzensorten, so dass jede wiederholt angebaut werden kann und der Ertrag gemessen wird. Die Kovarianzmatrix ist ein Maß für die Verwandtschaft durch genetische Ähnlichkeit, berechnet durch ibd-Wahrscheinlichkeiten in separaten Experimenten.

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

Mein Modell ist:

yld = gen + repl + error 

Sowohl gen als auch repl werden als zufällig betrachtet und ich möchte die zufälligen Effektschätzungen erhalten, die jedem Gen zugeordnet sind, jedoch muss ich die Beziehungsmatrix berücksichtigen.

Wenn es zu komplex ist, um verschachtelte Modelle anzupassen, würde ich repl einfach aus dem Modell entfernen, aber im Idealfall werde ich es behalten.

yld = gen +  error 

Wie kann ich dies mit R-Paketen erreichen, vielleicht mit nlme oder lme4? Ich weiß, dass ASREML das kann, aber ich habe keinen Halt und ich liebe R, weil es sowohl robust als auch frei ist.

Probieren Sie das kinshipPaket aus, auf dem basiert nlme. Weitere Informationen finden Sie in diesem Thread zu R-Sig-Mixed-Modellen. Ich hatte das vergessen, als ich versuchte, es für ein logistisches Modell zu tun. Ein ausgearbeitetes Beispiel finden Sie unter /programming/8245132 .

Für nicht normale Antworten müssen Sie das pedigreemm-Paket ändern, das auf lme4 basiert. Es bringt Sie in die Nähe, aber die Beziehungsmatrix muss aus einem Stammbaum erstellt werden. Die folgende Funktion ist eine Modifikation der pedigreemmFunktion, die stattdessen eine beliebige Beziehungsmatrix verwendet.

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
}

Die Verwendung ähnelt der, pedigreemmaußer dass Sie die Beziehungsmatrix als relmatArgument anstelle des Stammbaums als pedigreeArgument angeben.

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

Dies gilt hier nicht, da Sie zehn Beobachtungen / Einzelperson haben, aber für eine Beobachtung / Einzelperson benötigen Sie eine weitere Zeile in dieser Funktion und einen kleinen Patch lme4, um nur eine Beobachtung pro Zufallseffekt zu ermöglichen.

Diese Antwort ist eine mögliche Erweiterung des Vorschlags von Aaron, der vorgeschlagen hat, Pedigreem zu verwenden. Das Pedigreem kann die Beziehung aus den Projekten wie folgt berechnen: Ich weiß nicht, wie wir eine solche Beziehungsausgabe auf unterschiedliche Weise verwenden können.

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

Die gemischte Modellanpassung des Pakets basiert auf lme4, da die Syntax für die Hauptfunktion der lmer4-Paketfunktion lmer ähnelt, außer dass Sie das Stammbaumobjekt darin einfügen können.

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

Ich weiß, dass dies keine perfekte Antwort auf Ihre Frage ist, aber dies kann ein wenig helfen. Ich bin froh, dass Sie diese für mich interessante Frage gestellt haben!

lmer()in der lme4Packung erlaubt gekreuzte zufällige Effekte. Hier würden Sie so etwas wie verwenden

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

Für eine vollständige Referenz;

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