knitr::opts_chunk$set(echo = FALSE, results = 'asis') devtools::load_all()
r6obj_docstat <- rmddochelper::R6ClassDocuStatus$new() r6obj_docstat$set_current_status(psVersion = "0.0.901", psStatus = "Initialisation", psProject = "ZL_HS_2016") r6obj_docstat$include_doc_stat(psTitle = "## Document Status")
r6ob_abbrtable <- rmddochelper::R6ClassTableAbbrev$new() ### # include table of abbreviations only, if there are any if (!r6ob_abbrtable$is_empty_abbr()) r6ob_abbrtable$include_abbr_table(psAbbrTitle = "## Abbreviations")
In der Vorlesung wurden die linearen gemischten Modelle zur Voraussage von Zuchtwerten vorgestellt. Damit der Unterschied zwischen linearen Modellen mit nur fixen Effekten (Regressionsmodelle) und den linearen gemischten Modellen besser verständlich wird, wollen in dieser Übung nochmals anschauen, wie multiple lineare Regressionen an Daten angepasst werden.
Gegeben sei der folgende Datensatz.
nNrRecords <- 9 dfMlrData <- data.frame(Tochter = c(1:nNrRecords), Herde = c("1","1","2","2","2","3","3","3","3"), Vater = c("C","A","B","A","C","C","C","A","B"), Leistung = c(110,100,110,100,100,110,110,100,100)) ### # write data to outfile sMlrOutfile <- "mlrdata.csv" if (!file.exists(sMlrOutfile)) write.csv2(dfMlrData, file = sMlrOutfile) ### # create the table knitr::kable(dfMlrData)
Die Zahlen in der Kolonne Leistung der oben gezeigten Tabelle sollen als beobachtete phänotypische Leistungen aufgefasst werden. Im anzupassenden Regressionsmodell sollen Herde
und Vater
als fixe Effekte verwendet werden.
a) Stellen Sie das Regressionsmodell für die oben aufgelisteten Daten.
b) Berechnen Sie die Schätzungen für die fixen Effekte der Herde
und des Vaters
.
R
lautet lm()
lm()
muss über ein Dataframe
gemacht werden. Dies kann entweder über die direkte Angabe der Daten in der Funktion data.frame()
gemacht werden, oder die Daten können auch mit read.csv2()
von der Webseite (https://charlotte-ngs.github.io/LBGHS2016/w10/mlrdata.csv) eingelesen werden. Die direkte Angabe des Dataframes sieht wie folgt aus:nNrRecords <- 9 dfMlrData <- data.frame(Tochter = c(1:nNrRecords), Herde = c("1","1","2","2","2","3","3","3","3"), Vater = c("C","A","B","A","C","C","C","A","B"), Leistung = c(110,100,110,100,100,110,110,100,100))
Leistung ~ -1 + Herde + Vater
summary()
angezeigt werden.a) Das Regressionsmodell lautet
$$y = Xb + e$$
wobei der Beobachtungsvektor $y$
vecYTask1 <- dfMlrData$Leistung cat("$$y = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecYTask1, ncol = 1), pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
Die Inzidenzmatrix $X$
matXTask1 <- matrix(data = c(1,0,0,0,0,1, 1,0,0,1,0,0, 0,1,0,0,1,0, 0,1,0,1,0,0, 0,1,0,0,0,1, 0,0,1,0,0,1, 0,0,1,0,0,1, 0,0,1,1,0,0, 0,0,1,0,1,0), nrow = nNrRecords, byrow = TRUE) cat("$$X = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matXTask1, pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
Der Vektor $b$ der fixen Effekte
vecFixEffects <- c("Herde_1", "Herde_2", "Herde_3", "Vater_A", "Vater_B", "Vater_C") cat("$$b = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecFixEffects, ncol = 1), pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
b) Angenommen, die Daten sind in einem Dataframe namens dfMlrData
gespeichert. Dann wird das Regressionsmodell mit folgenden Anweisungen angepasst.
lmMlrData <- lm(Leistung ~ -1 + Herde + Vater, data = dfMlrData) summary(lmMlrData)
Anstelle des Regressionsmodells in der Aufgabe 1 verwenden wir jetzt ein gemischtes lineares Modell. In diesem Modell werden die Herden als fixe Effekte und die Väter als zufällige Effekte modelliert. In der Vorlesung werden wir diese Art von Modellen noch genauer unter dem Namen "Vatermodell" anschauen. Die Covarianzmatrix $G$ der zufälligen Vatereffekte setzen wir dabei zu
$$G = 0.1*I$$
wobei $I$ die Einheitsmatrix bezeichnet.
a) Stellen Sie die Mischmodellgleichungen für die Daten in Aufgabe 1 auf. b) Berechnen Sie die Schätzer für die fixen Effekte der Herde und machen Sie eine Voraussage für die zufälligen Vatereffekte
a) Die Mischmodellgleichungen lauten
$$\left[ \begin{array}{rr} X^TX & X^TZ\ Z^TX & Z^TZ + G^{-1} \end{array} \right] \left[ \begin{array}{c} \hat{Herde}_1\ \hat{Herde}_2\ \hat{Herde}_3\ \hat{Vater}_A\ \hat{Vater}_B\ \hat{Vater}_C \end{array} \right] = \left[ \begin{array}{c} X^Ty\ Z^Ty \end{array} \right] $$
Abgekürzt können wir die Mischmodellgleichungen schreiben als
$$M * \hat{s} = r$$
Die einzelnen Komponenten lauten
matXTask2 <- matXTask1[,1:3] matXtXTask2 <- crossprod(matXTask2) cat("$$X^TX = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matXtXTask2, pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
matZTask2 <- matXTask1[,4:6] matXtZTask2 <- crossprod(matXTask2,matZTask2) cat("$$X^TZ = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matXtZTask2, pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
matZtZTask2 <- crossprod(matZTask2) cat("$$Z^TZ = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matZtZTask2, pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
$$G^{-1} = 10*I$$
vecXtYTask2 <- crossprod(matXTask2,vecYTask1) vecZtYTask2 <- crossprod(matZTask2,vecYTask1) cat("$$X^Ty = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecXtYTask2, ncol = 1), pnDigits = 0), collapse = "\n")) cat("\\right]\n$$") cat("$$Z^Ty = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecZtYTask2, ncol = 1), pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
Stellen wir alle Komponenten zusammen erhalten wir
matCoeffTask2 <- cbind(rbind(matXtXTask2,t(matXtZTask2)),rbind(matXtZTask2,matZtZTask2 + 10*diag(1,nrow = nrow(matZtZTask2)))) vecRhsTask2 <- c(vecXtYTask2,vecZtYTask2) cat("$$\\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matCoeffTask2, pnDigits = 0), collapse = "\n")) cat("\\right]\n") vecFixEffectsHat <- c("\\hat{Herde}_1", "\\hat{Herde}_2", "\\hat{Herde}_3", "\\hat{Vater}_A", "\\hat{Vater}_B", "\\hat{Vater}_C") cat("\\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecFixEffectsHat, ncol = 1)), collapse = "\n")) cat("\\right]\n") cat(" = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecRhsTask2, ncol = 1), pnDigits = 0), collapse = "\n")) cat("\\right]\n$$")
b) Die Lösungen berechnen wir als
$$\hat{s} = M^{-1} * r$$
vecSolTask2 <- solve(matCoeffTask2,vecRhsTask2) cat("$$\\hat{s} = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = as.matrix(vecSolTask2, ncol = 1), pnDigits = 4), collapse = "\n")) cat("\\right]\n$$")
r6ob_abbrtable$writeToTsvFile()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.