Kovarianzstrukturen
library(data.table) # bessere Datenmanipulation
library(nlme) # gemischtes Modell package 1
library(asreml) # gemischtes Modell package 2 - benötigt R version 3.2.3. und LizenzDatensatz
In diesem Experiment wurde in einer randomisierten vollständigen Blockanlage (RCBD) mit 5 Wiederholungen der Blattflächenindex von 5 Sorghumsorten verglichen. Allerdings wurde der Blattflächenindex mehrfach, nämlich in 5 aufeinanderfolgenden Wochen, gemessen, sodass Messwiederholungen vorliegen.
Dieses Beispiel basiert auf Example 4 des
agriTutorialpackages und der dazugehörigen Veröffentlichung
Piepho, H. P., & Edmondson, R. N. (2018). A tutorial on the statistical analysis of factorial experiments with qualitative and quantitative treatment factor levels. Journal of Agronomy and Crop Science.
Messwiederholungen werfen eine Neuerung gegenüber den vorangegangenen Beispielen auf: zum ersten Mal ist die kleinste Randomisationseinheit (=Parzelle) nicht gleichzeitig die Beobachtungseinheit, da wir mehrere Beobachtungen pro Parzelle haben. Der wichtige Punkt hier ist, dass der Faktor Woche nicht randomisiert werden kann. Statt der üblichen Annahme von unabhängigen Messwerten, sind Messwerte derselben Parzelle von aufeinanderfolgenden Wochen wahrscheinlich miteinander korreliert. Um dies zu modellieren, soll im ersten Schritt vorerst das Modell für die Analyse einer einzelnen Woche aufgestellt werden.
print(repmes, nrows=10)## week gen rep plot y row col
## 1: 1 A rep1 plot1 5.00 1 1
## 2: 1 D rep1 plot2 5.86 2 1
## 3: 1 B rep1 plot3 5.82 3 1
## 4: 1 C rep1 plot4 5.65 4 1
## 5: 1 C rep2 plot5 5.39 1 2
## ---
## 96: 5 B rep4 plot16 3.99 4 4
## 97: 5 D rep5 plot17 2.96 1 5
## 98: 5 C rep5 plot18 3.00 2 5
## 99: 5 B rep5 plot19 2.95 3 5
## 100: 5 A rep5 plot20 2.16 4 5str(repmes, width=40, strict.width="cut")## Classes 'data.table' and 'data.frame': 100 obs. of 7 variables:
## $ week: Factor w/ 5 levels "1","2",""..
## $ gen : Factor w/ 4 levels "A","B",""..
## $ rep : Factor w/ 5 levels "rep1","r"..
## $ plot: Factor w/ 20 levels "plot1","..
## $ y : num 5 5.86 5.82 5.65 5.39 4...
## $ row : num 1 2 3 4 1 2 3 4 1 2 ...
## $ col : num 1 1 1 1 2 2 2 2 3 3 ...
## - attr(*, ".internal.selfref")=<exter..Kovarianzstrukturen
nlme & asreml-R
Wir werden hier die Modelle mit der gls() Funktion aus dem nlme package anpassen. Zu jedem Modell wird auch das Gegenstück in asreml() Syntax gezeigt, allerdings nicht ausgeführt. Das asreml-R version 3.0 package funktioniert nur mit älteren R-Version 3.2.3, die man hier herunterladen kann. Außerdem muss man eine Lizenz besitzen. An der Universität Hohenheim gibt es diese im ILIAS und man muss auch nach der Installation mit dem Hohenheimer VPN verbunden sein, damit das package funktioniert. Der Code dieses Beispiels ist auch hier auf meinem GitHub verfügbar.
Mehr Infos zu Varianzstrukturen gibt es beispielsweise in Prof. Piephos Skript zu gemischten Modellen Kapitel 6.
ID: unabhängige, homogene Varianzen
nlme
gls.id <- gls(y ~ week + week*gen + week*rep,
data=repmes)
gls.id$sigma^2 # ID.var## [1] 0.02317583asreml-R v3
asr.id <- asreml(data = repmes,
fixed = y ~ week + week:gen + week:rep)
summary(asr.id)$varcomp[,c(2,3,5)] # ID.varDIAG: unabhängige, heterogene Varianzen
nlme
gls.dg <- gls(y ~ week + week*gen + week*rep,
weights = varIdent(form = ~ 1|week),
data=repmes)
gls.dg$sigma^2 * c(1, coef(gls.dg$modelStruct$varStruct, unc=F))^2 # DG.vars## 2 3 4 5
## 0.01974885 0.02056508 0.02230259 0.03269346 0.02056925asreml-R v3
asr.dg <- asreml(data = repmes,
fixed = y ~ week + week:gen + week:rep,
rcov = ~ at(week):plot)
summary(asr.dg)$varcomp[,c(2,3,5)] # DG.varsAR1: first order autoregressive
nlme
gls.ar <- gls(y ~ week + week*gen + week*rep,
corr = corExp(form = ~ week|plot),
data=repmes)
gls.ar$sigma^2 # AR1.var ## [1] 0.02317583as.numeric(exp(-1/coef(gls.ar$modelStruct$corStruct, unconstrained=F))) # AR1.cor## [1] 0.7776763asreml-R v3
asr.ar <- asreml(data = repmes,
fixed = y ~ week + week:gen + week:rep,
rcov = ~ ar1(week):plot)
summary(asr.ar)$varcomp[,c(2,3,5)] # AR var + corCS: compound symmetry
nlme
gls.cs <- gls(y ~ week + week*gen + week*rep,
corr = corCompSymm(form = ~ week|plot),
data=repmes)
gls.cs$sigma^2 # CS.var ## [1] 0.02317584as.numeric(coef(gls.cs$modelStruct$corStruct, unconstrained=F)) # CS.cor## [1] 0.7007553asreml-R v3
asr.cs <- asreml(data = repmes,
fixed = y ~ week + week:gen + week:rep,
rcov = ~ cor(week):plot)
summary(asr.cs)$varcomp[,c(2,3,5)] # CS var + corUN: unstrukturiert
nlme
gls.un <- gls(y ~ week + week*gen + week*rep,
corr = corSymm(form = ~ 1|plot),
weights = varIdent(form = ~ 1|week),
data=repmes)
c(gls.un$sigma^2, coef(gls.un$modelStruct$varStruct, unconstrained=T)) # UN.vars## [1] 0.01974916 0.02023961 0.06079338 0.25202993 0.02034066gls.un$modelStruct$corStruct # UN.cors## Correlation structure of class corSymm representing
## Correlation:
## 1 2 3 4
## 2 0.708
## 3 0.644 0.634
## 4 0.748 0.744 0.934
## 5 0.534 0.549 0.718 0.776asreml-R v3
asr.un <- asreml(data = repmes,
fixed = y ~ week + week:gen + week:rep,
rcov = ~ us(week):plot)
summary(asr.un)$varcomp[, c(2,3,5)] # UN vars + corsKomplexere Kovarianzsstrukturen
In asreml, nicht aber in nlme, ist es möglich komplexere Kovarianzsstrukturen zu erzeugen. So können bis zu 3 Kovarianzstrukturmatrizen mittels Kronecker-Produkt miteinander kombiniert werden. Siehe dazu Seite 11 der “asreml-R.pdf” package Dokumentation, sowie die Beispielcodes hier.
Bei Fragen kannst du mir gerne schreiben!
schmidtpaul@hotmail.de
