Die Gleitmethode nach Mitscherlich und was sie mit Geostatistik zu tun hat

Mitscherlich's slide method and what it has to do with geostatistics

Journal für Kulturpflanzen, 72 (10-11). S. 527–540, 2020, ISSN 1867-0911, DOI: 10.5073/JfK.2020.10-11.03, Verlag Eugen Ulmer KG, Stuttgart

Kontaktanschrift

Prof. Dr. Hans-Peter Piepho, Universität Hohenheim, Institut für Kulturpflanzenwissenschaften, Fruwirthstr. 23, 70599 Stuttgart, E-Mail: piepho@uni-hohenheim.de

Zur Veröffentlichung angenommen

27. August 2020

Dies ist ein Open-Access-Artikel, der unter den Bedingungen der Creative Commons Namensnennung 4.0 International Lizenz (CC BY 4.0) zur Verfügung gestellt wird (https://creativecommons.org/licenses/by/4.0/deed.de). This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/deed.en).

Einleitung

Eine der ältesten Anlagen für Feldversuche ist die Mitscherlich-Anlage (1925; siehe auch hier: Möller-Arnold & Feichtinger, 1929; Wermke, 1959; von Lochow & Schuster, 1961). Bei dieser Anlage werden die Versuchsglieder (Behandlungen) in systematischer Reihenfolge nebeneinander geprüft. In Abb. 1 findet sich ein Beispiel mit fünf Behandlungen (A-E) in fünf Wiederholungen.

Abb. 1. Eine systematische Versuchsanlage (Mitscherlich-Anlage) mit fünf Behandlungen und fünf Wiederho­lungen für die Auswertung mit der Gleitmethode (von Lochow & Schuster, 1961, S. 5, Abb. 2).

Für die Ausgleichsrechnung mit der Gleitmethode nach Mitscherlich wird nun ein gleitendes Fenster der Größe eines vollständigen Blocks über die Anlage gelegt. So bilden die ersten fünf Parzellen in Abb. 1 einen vollständigen Block. Wird dieses Fenster um eine Parzelle nach rechts verschoben, ergibt sich wieder ein vollständiger Block. Der Mittelwert der jeweils fünf Prüfglieder des gleitenden Blocks liefert ein Maß für den räumlichen Trend der Ertragsfähigkeit von links nach rechts. Die beobach­teten Ertragswerte können mit diesem Mittelwert korrigiert werden, typischerweise durch Differenzbildung (von Lochow & Schuster, 1961, S. 96ff.) oder Divi­sion (von Boguslawski, 1942). Solche korrigierten Werte werden für alle möglichen Platzierungen des gleitenden Blocks berechnet und dann je Behandlung gemittelt, um bereinigte Behandlungseffekte zu berechnen. Wird zu diesen Schätzwerten ein Gesamtmittelwert hinzuaddiert, ergeben sich korrigierte Behandlungsmittelwerte. Aus den einzelnen korrigierten Werten werden darü­ber hinaus Standardabweichungen und Standardfehler sowie Grenzdifferenzen für die Mittelwerte berech­net (von Lochow & Schuster, 1961, S.99).

Die Gleitmethode von Mitscherlich stellt ein Verfahren für Anlagen dar, bei denen die Parzellen wie in Abb. 1 in einer einzigen Reihe angeordnet sind. Es wurde von von Boguslawski (1942) für zweidimensionale Anlagen unter der Bezeichnung „Ertragsflächenmethode“ erweitert. Hierbei liegt der Fokus auf zwei Behandlungsfaktoren, von denen der eine systematisch in den Zeilen und der andere in den Spalten variiert wird. Das gleitende Fenster umfasst jeweils alle Behandlungen einmal, wobei die Zahl der Zeilen des Fensters der Zahl der Stufen des einen Faktors entspricht und die Zahl der Spalten der Zahl der Stufen des anderen Faktors. Die vorgeschlagene Analyse lehnt sich an die Gleitmethode von Mitscherlich an, wobei von Boguslawski für die Trendkorrektur Quo­tien­ten aus Parzellenwert und Mittelwert des gleitenden Blocks verwendet.

Die Gleitmethode von Mitscherlich wie auch die Ertrags­flächenmethode von von Boguslawski bringen es mit sich, dass die Behandlungen systematisch angeordnet werden müssen und keinerlei Randomisation erfolgt, damit die gleitenden Blöcke jeweils vollständig sind. Damit stehen diese Methoden im Gegensatz zu der von R. A. Fisher (1925) eingeführten Methode zur Anlage von Feldversuchen, welche eine zufällige Zuordnung der Behandlungen zu den Parzellen erfordert. Da Fishers Prinzipien der Versuchsplanung sich inzwischen nicht nur in der Landwirtschaft als Standard etabliert haben (John & Williams, 1995; Thomas, 2006; Bailey, 2008), mag die Befassung mit den systematischen Versuchsan­lagen nach Mitscherlich und von Boguslawski aus heutiger Sicht als von eher historischem Interesse erscheinen (Thomas, 2006). Es gibt allerdings noch aktuell laufende Dauerversuche, die nach diesen systematischen Plänen angelegt wurden und auch noch heute ausgewertet werden (Macholdt et al., 2019, 2020). Daher stellt sich die Frage nach der angemessenen Auswertung solcher Versuche weiterhin.

Ziel der vorliegenden Arbeit ist es, ein neues Licht auf die von Mitscherlich vorgeschlagene Auswertungsmethode zu werfen und eine Brücke zu moderneren geostatistischen Verfahren der Modellierung von Daten aus Feldversuchen zu schlagen. Nach einem einführenden Beispiel mit einem eindimensionalen Versuchsplan wird ein zweites Beispiel betrachtet, bei dem der Versuch in Zeilen und Spalten angelegt ist.

Die Grundidee der Gleitmethode von Mitscherlich ist, dass gleitende Blöcke es erlauben, die beobachteten Parzellendaten um Bodentrends zu bereinigen. Falls die Korrektur mit Differenzen erfolgt (von Lochow & Schuster, 1961, S. 96–97), so kann diese für einen gleitenden Block geschrieben werden als

wobei y_i der Parzellenwert für die i-te Behandlung ist (i = 1, 2,..., v) und

der arithmetische Mittelwert über alle v Beobachtungen des Blocks, wobei v gleich der Zahl der Behandlungen ist. Bei dieser Formulierung, welche lediglich einen einzigen gleitenden Block betrachtet, ist eine Indizierung für die Blöcke der Einfachheit halber weggelassen worden. Es wird außerdem vorausgesetzt, dass der gleitende Block jeweils ein vollständiger ist, dass also jede Behandlung genau einmal vorkommt.

Die Formulierung der Mitscherlich-Methode basiert ausschließlich auf Beobachtungswerten und deren direkter Verrechnung. Es gibt keinerlei statistisches Modell, welches von Mitscherlich explizit formuliert worden wäre, um sein Verfahren zu begründen oder herzuleiten. Aller­dings impliziert die Betrachtung von vollständigen Blöcken zur Korrektur nach Gleichung (1), dass Blockmittelwerte es erlauben, einen räumlichen Trend zu erfas­sen. Der Blockmittelwert, um den die Parzellen­werte dabei korrigiert werden, kann als Trendwert aufgefasst werden. Um diesen Gedanken zu konkretisieren, ist es hilfreich, ein lineares Modell für die Parzellenwerte zu formulieren. Hierzu gehen wir zunächst von der Situa­tion aus, dass auf allen Parzellen dieselbe Behandlung geprüft wird. Dies führt zu einem sog. Nullmodell (Nelder, 1965), also einem Modell, welches bei Abwesenheit von Behandlungseffekten gilt. Für einen betrachteten Block können wir das Mittel über die Parzellenwerte als Trendwert für den gesamten Block betrachten. Für diesen Trendwert kann das Modell + b verwendet werden, wobei ein Gesamtmittel aller Parzellenwerte über den ganzen Versuch ist und b die Abweichung des Trendwertes für den betrachteten Block vom Gesamtmittel, also der Blockeffekt. Da es bei der Gleitmethode nicht nur einen gleitenden Block gibt, ist es sinnvoll, einen Index für den Block einzuführen. Dies führt dann zu dem Modell + b_j für den im j-ten gleitenden Block für jede Parzelle erwarteten Beobachtungswert, wobei der Index von j = 1 bis zur Zahl der Blöcke läuft. Allerdings variieren die Beobachtungswerte innerhalb eines Blocks von Parzelle zu Parzelle. Um diese Schwankung abzubilden, führen wir einen Parzelleneffekt e ein, der auch als Parzellenfehler oder Residualeffekt bezeichnet wird. Das Nullmodell für den Beobachtungswert in der i-ten Behandlung im j-ten Block lautet dann

Bis hierher wurde angenommen, dass es keine Behandlungseffekte gibt. Wenn den Effekt der i-ten Behandlung bezeichnet, dann lautet das um diesen Effekt erweiterte Modell

Nachdem nun ein lineares Modell für die Beobachtungswerte formuliert ist, können wir die Korrektur in Gl. (1) untersuchen, indem wir dieses Modell anstelle der Beobachtungswerte dort einsetzen, wobei der Index für den Blockeffekt hinzugefügt wird:

Offensichtlich bereinigt die Korrektur die Beobachtungswerte um die Trendwerte + b_j. Dies ist die Kernidee des Verfahrens von Mitscherlich, auch wenn er für dessen Begründung kein lineares Modell verwendet hat.

Das Modell in Gleichung (4) ist unschwer zu erkennen als dasselbe Modell, welches auch für eine randomisierte vollständige Blockanlage verwendet wird (Fisher, 1925; John & Williams, 1995). Allerdings besteht ein wesentlicher Unterschied darin, dass eine Blockanlage disjunkte, also räumlich eindeutig getrennte Blöcke hat, während die Gleitmethode überlappende Blöcke betrachtet. Dies bedeutet, dass bei der Gleitmethode jede Parzelle nicht nur einem Block zuzuordnen ist, sondern v verschiedenen Blöcken. Diese Tatsache ist in Gleichung (4) so nicht abgebildet, weil wir bisher nur einen einzigen gleitenden Block betrachtet haben, weswegen dieses vereinfachte Modell verwendet werden konnte. Hierbei sind die Effekte der weiteren Blöcke, denen eine Parzelle im Gleitverfahren zugeordnet wird, durch den Residualeffekt e_ij subsumiert. Aber wegen der überlappenden Blöcke ist das Modell für die gesamte Anlage in dieser Form noch nicht vollständig und muss im Folgenden noch erweitert werden. Hierzu muss die mathematische Notation entsprechend modifiziert werden. Bisher wurde der Beobachtungswert nach dem Block j und der Behandlung i indiziert. Diese Indizierung wird jetzt ersetzt durch einen einzigen Index k für die Parzelle. Es ist dann y_k der Beobachtungswert auf der k-ten Parzelle. Dabei läuft der Index k von 1 bis N, wobei N die Zahl der Parzellen des Versuchs bezeichnet. Ein lineares Modell, welches alle Effekte derjenigen Blöcke umfasst, denen die k-te Par­zelle zugeordnet ist, kann wie folgt geschrieben werden:

Hierin bezeichnet M_k die Menge aller Indices j für diejenigen gleitenden Blöcke, welche der k-ten Parzelle zugeord­net sind, i(k) ist der Index für die Behandlung, welche der k-ten Parzelle zugeordnet ist, der dazugehörige Behandlungseffekt, und e_k ist der Residualeffekt, also die Abweichung des k-ten Parzellenwertes y_k vom erwarteten Wert

Da die Gleitmethode vollständige Blöcke verwendet, welche jeweils alle v Behandlungen umfassen, sind jedem gleitenden Block v Parzellen zugeordnet. Umgekehrt bedeu­tet dies auch, dass jeder Parzelle v verschiedene Blockeffekte zugeordnet sind, dass jede Menge M_k also v Blockindices umfasst. Dabei haben unmittelbar benachbarte Parzellen, die nur einen „Parzellenschritt“ von einander entfernt sind, jeweils v–1 Blockeffekte oder Trendwerte b_j gemeinsam. Parzellen, die zwei Parzellenschritte voneinander entfernt sind, haben noch v–2 Blockeffekte gemeinsam, usw. Parzellen, die v Parzellenschritte oder mehr voneinander entfernt sind, haben keinen Block­effekt mehr gemeinsam. Generell ist die Zahl der übereinstimmenden Blockeffekte für zwei Parzellen k und k' gegeben durch

also gleich dem Maximum der Werte 0 und der Zahl der Behandlungen v abzüglich Parzellenschritte d_{k, k'} = |k – k'|. Je mehr übereinstimmende Blockeffekte zwei Parzellen haben, um so ähnlicher sind sich die Beobachtungswerte gemäß dem Modell (6). Daher bildet dieses Modell die Erwartung ab, dass Parzellen um so ähnlichere Ertragsfähigkeit haben, je näher sie sich räumlich sind.

Der Vollständigkeit halber sei angemerkt, dass wir oben zwar die Annahme getroffen haben, dass jede Parzelle v Blocks zugeordnet ist, dass dies aber physisch offen­sichtlich nicht für die Randparzellen zutrifft. Jedoch ist es ohne weiteres möglich, von einer hypothetischen Erweiterung der Anlage nach links und rechts auszugehen dergestalt, dass auch die jeweils am Rande liegenden v–1 Parzellen jeweils v (zum Teil hypothetischen) Blöcken zugeordnet werden können. Damit lässt sich das aus unserem Ansatz hergeleitete Modell tatsächlich auf alle Parzellen anwenden.

Parzellen mit übereinstimmenden Blockeffekten haben außerdem eine Kovarianz. So haben zwei benachbarte Parzellen wie oben bemerkt v–1 Blockeffekte gemeinsam und daher eine Kovarianz von . Bei zwei Parzellenschritten beträgt die Kovarianz noch usw. Im allgemeinen ist die Kovarianz gleich

wobei h_{k, k'} in Gleichung (7) angegeben ist. Je höher die Kovarianz, desto ähnlicher sind die beiden Parzellen k und k' in den Beobachtungswerten. Nach v Parzellenschritten von der k-ten zur k'-ten Parzelle ist die Kova­rianz gleich Null. Bis zu diesem Abstand, der Reichweite (engl. range) des Modells, fällt die Kovarianz linear mit dem Abstand ab, um dann beim Wert Null stehen zu bleiben (Abb. 2).

Abb. 2. Kovarianzfunktionen für drei geostatistische Modelle für das erste Beispiel (Parameter­schätzwerte aus Tab. 2). AR(1) = autoregressives Modell erster Ordnung; LIN = lineares geostati­stisches Modell (rho frei: ρ geschätzt mit der REML Methode; rho = 0.2: Fixierung ρ = 1/5).

Dieses Muster einer mit dem Abstand sinkenden Kovarianz ist charakteristisch für eine Vielzahl von Modellen der Geostatistik, welche sich lediglich in der genauen funktionalen Form dieses Abfallens der Kovarianz unterscheiden. Die Stärke des Abfallens hängt dabei von den Modellparametern ab. Je schneller der Abfall erfolgt, desto geringer ist die Reichweite der räumlichen Korrelation. Die meisten Modelle nehmen einen nichtlinearen Verlauf der Kovarianzfunktion an (Schabenberger & Gotway, 2005; für eine anwendungsorientierte Einführung siehe auch Slaets et al., 2020). Das hier beschriebene Modell ist dagegen ein Spezialfall des linearen Kova­rianzmodells der Geostatistik (Williams, 1986). Die Reichweite dieses linearen geostatistischen Modells ist gegeben durch die Größe der gleitenden Blöcke und somit durch die Zahl der Behandlungen fixiert. Es besteht allerdings kein Grund zu der Annahme, dass dies der beste anzunehmende Wert für die Reichweite ist. Besser ist es, andere Werte für die Reichweite zuzulassen und den optimalen Wert aus den Daten zu schätzen. Dies kann mit Standardverfahren wie der Residual Maximum Likelihood (REML) Methode erfolgen (Schabenberger & Gotway, 2005). Es ist außerdem möglich Modelle anzupassen, welche anstatt eines linearen einen nichtlinearen Abfall der Kovarianz mit der räumlichen Distanz annehmen. Diese Betrachtungen stellen die Gleitmethode von Mitscherlich in einem gemeinsamen Zusammenhang mit geostatistischen Methoden, die heute routinemäßig für die Analyse von Feldversuchen eingesetzt werden (Gilmour et al., 1997; Piepho et al., 2008; Rodríguez-Álvarez et al., 2018). Wegen seiner Linearität besteht eine besonders enge Beziehung zu Verfahren der sog. Nächstnachbaranalyse, deren Entwicklung parallel mit der Erfindung der Varianzanalyse in den 30er Jahren des letzten Jahrhunderts vor allem im angelsächsischen Bereich ihren Anfang nahm und die auch heute noch eine große Bedeu­tung hat (Piepho et al., 2008).

Ein lineares Kovarianzmodell wird allerdings nicht zum selben Analyseergebnis führen wie die Gleitmethode von Mitscherlich, selbst wenn die Kovarianz-Struktur der Gleichungen (8) und (9) angenommen wird. Es ist nicht ohne weiteres klar, welches Analyseverfahren die besseren Ergebnisse liefert. Für den geostatistischen Ansatz spricht allerdings, dass ihm ein statistisches Modell zugrunde gelegt werden kann, während es für die Gleitmethode von Mitscherlich in der ursprünglich vorgeschlagenen Form kein Modell gibt. Auf Basis eines wohldefinierten statistischen Modells kann auch eine statistische Analyse mit etablierten Verfahren durchgeführt werden, während es für die Gleitmethode keine solche Basis gibt. Zwar wird bei von Lochow & Schuster (1961) vorgeschlagen, wie Standardfehler und Grenzdifferenzen für Behandlungsmittelwerte berechnet werden können, es besteht aber wegen des fehlenden Modells keine Basis für die Ableitung einer weitergehenden inferenz-statistischen Analyse, etwa mit einer Varianzanalyse oder multiplen Mittelwertvergleichen. Diese Lücke wird durch den hier gemachten Brückenschlag zur Geo­statistik geschlossen.

Die Kovarianz nach dem linearen geostatistischen Modell ist in der MIXED Prozedur wie folgt parametrisiert:

Hierbei ist eine räumliche Varianzkomponente, d_k,k' = |k – k'| ist die Distanz zwischen den Parzellen k und k', und ist ein Parameter, der den Abfall der Kovarianz mit der Distanz bestimmt. Die Funktion 1(.) ist die Indikator-Funktion, welche den Wert 1 annimmt, falls der Ausdruck im Argument wahr ist und sonst den Wert 0. Dieser Faktor stellt hier sicher, dass die Kovarianz Null wird, wenn d_k,k' > 1 wird. Ein Koeffizientenvergleich mit Gleichung (9) führt zu den Entsprechungen und  = v^–1. Hierbei stellt die Spezifikation  = v^–1 sicher, dass für Parzellen mit d_k,k' = |k – k'| = v, also für Parzellen, die v Parzellenschritte voneinander entfernt sind, d_k,k' = 1 ist und somit die Kovarianz gerade auf Null abgefallen ist. Alternativ kann der Reichweiten-Parameter mit der REML-Methode auch aus den Daten geschätzt werden.

Die Varianzparameterschätzwerte nach dem linearen geostatistischen Modell sind in Tab. 2 angegeben. Zum Vergleich haben wir auch das AR(1)-Modell angepasst (; Piepho et al., 2015b), ein Modell mit unabhängigen Fehlern sowie ein Modell für eine Bockanlage mit zufälligen Effekten für sechs disjunkte Blöcke. Das AR(1)-Modell ist das am häufigsten verwendete nichtlineare geostatistische Modell für Feldversuche (Gilmour et al., 1997), weswegen es hier betrachtet wird. Bei diesem Modell sinkt die Kovarianz exponentiell mit dem räumlichen Abstand (Abb. 2). Gängige Maße für die Güte der Modellanpassung bei Verwendung der Likelihood-Methode sind die Devianz (minus der zweifache Wert der residualen logarithmierten Likelihood; je kleiner der Wert, umso besser die Anpassung) und das Akaikesche Informationskriterium (AIC = Devianz + 2p, wobei p = Zahl der Varianzparameter; je kleiner der Wert, umso besser das Modell; Wolfinger, 1996). Die Devianzen wie auch die AIC-Werte in Tab. 2 zeigen, dass die geostatistischen Modelle für diesen Datensatz vergleichsweise gut passen. Die Kovarianzfunktionen nach den geostatistischen Modellen aus Tab. 2 sind in Abb. 2 dargestellt. Es ist erkennbar, dass nach allen diesen Modellen die Kovarianz mit der Distanz sinkt. Dabei ist die Übereinstimmung des AR(1)-Modells mit dem linearen Modell relativ gut, wenn der Parameter des letzteren frei geschätzt wird (rote, gestrichelte Linie), während mit der Fixierung bei  = 1/5 die Übereinstimmung weniger gut ausfällt. In Tab. 3 finden sich die Behandlungsmittel­werte nach verschiedenen Verfahren und Tab. 4 zeigt die paarweisen Differenzen der Behandlungsmittelwerte nach dem linearen geostatistischen Modell bei Schätzung des Reichweitenparameters und bei dessen Fixierung beim Wert  = 1/5. In Tab. 3 sind zum Vergleich außer­dem die arithmetischen Mittelwerte sowie die in von Lochow & Schuster (1961) mit der Gleitmethode berech­neten Mittelwerte aufgeführt. Für das lineare geostatistische Modell finden wir mit fixiertem Reichweiten-Parameter eine etwas höhere und damit schlechtere Devianz als bei freier Schätzung mit  = 0.398 (Tab. 2). Die geschätzte Reichweite (1/ Parzellenschritte; bei dieser Distanz ist die Kovarianz auf Null gesunken) ist etwa halb so weit wie für das Modell mit fixiertem ; nach etwa nach 2.5 Parzellenschritten ist die Kovarianz bereits auf Null abgesunken (Abb. 1), was ein ähnlich schneller Abfall ist wie beim AR(1)-Modell, bei dem die Kovarianz sich allerdings nur asymptotisch der Null nähert. Wegen der geringen Reichweite bei Schätzung des Parameters liegen die adjustieren Mittelwerte relativ nahe bei den arithmetischen Mittelwerten (Tab. 3). Die paarweisen Differenzen der Behandlungsmittelwerte haben bemerkenswerterweise für beide linearen geostatistischen Modelle relativ ähnliche Varianzen (Tab. 4). Die kompletten SAS Anweisungen sind im Elektronischen Anhang verfügbar, sowohl für das erste Beispiel als auch für das später folgende zweite Beispiel.

Die Nähe des Mitscherlichschen Ansatzes zu geostatistischen Modellen bietet allerdings die Möglichkeit, die genannten Restriktionen auf ihre Implikationen für die Effizienz des Designs hin zu untersuchen. Eine wichtige Forderung beim Design von Versuchen ist, dass alle paarweisen Behandlungsvergleiche möglichst mit derselben Genauigkeit durchgeführt werden können. Unter einem geostatistischen Modell ist die Varianz einer Differenz zweier Parzellenwerte um so kleiner je näher die Parzellen für die beiden Behandlungen beieinanderliegen. Bei einer systematischen Anordnung wie in Abb. 1 ist jede Behandlung immer von denselben anderen Behandlungen umgeben. Die Behandlung C in Abb. 1 findet sich zum Beispiel immer in der direkten Nachbarschaft der Behandlungen B und D. Dagegen ist die Behandlung C nie in direkter Nachbarschaft zu den Behandlungen A und E. Daher könnte erwartet werden, dass Vergleiche von C mit B und D genauer sind als die Vergleiche von C mit den Behandlungen A und E, was gegen die die systematische Anordnung sprechen könnte, abgesehen von der fehlenden Randomisierung. Zwar zeigte die Analyse für das erste Beispiel, dass die Standardfehler einer Differenz bei dem angenommenen linearen geostatistischen Modell relativ homogen sind, aber diese folgen trotzdem dem erwarteten Muster, da die Vergleiche A-B, B-C, C-D und D-E die kleineren Standardfehler haben (Tab. 4). Die Ergebnisse für Vergleiche der Behandlungsmittelwerte nach dem AR(1)-Modell sind relativ ähnlich denen nach dem linearen geostatistischen Modell mit geschätztem Wert für den Parameter . Die Standardfehler einer Differenz für das Modell einer Blockanlage und das Modell mit unabhängigen Fehlern betragen 0.1242 und 0.1468. Diese Werte liegen höher als die durchschnittlichen Standardfehler einer Differenz nach den geostatistischen Model­len in Tab. 4.

Interessant ist es in diesem Zusammenhang zu fragen, was unter dem oben abgeleiteten linearen geostatistischen Modell das optimale Design ist. Es ist mit heute gängigen numerischen Verfahren möglich, für jedes vorgegebene Modell bezüglich der Varianzen und Kovarianzen zwischen den Parzellen ein jeweils optimales Design zu suchen. Dazu wird jeweils ein Zielkriterium definiert, welches in der computer-gestützten Suche optimiert wird. Zwei gängige Kriterien sind die D-Effizienz und die A-Effizienz (John & Williams, 1995). Beide beruhen auf der sog. Informations-Matrix der Behandlungsmittelwerte. Die D-Effizienz entspricht der Determinante dieser Matrix, und das Ziel ist diese zu maximieren. Die A-Effizienz betrachtet die durchschnittliche paarweise Varianz einer Mittelwertdifferenz und das Ziel ist diese zu minimieren.

Wenn wir nun für das im vorangegangenen Abschnitt analysierte Beispiel die geschätzte Varianzstruktur für das lineare geostatistische Modell mit fixiertem Reichweiteparameter = 1/5 zugrunde legen und das optimale Design für fünf Behandlungen mit jeweils sechs Wiederholungen suchen, so finden wir für beide Kriterien die systematische Anordnung wie in Abb. 1 als das optimale Design! Dabei haben wir im Fall der D-Optimalität mit der SAS Prozedur OPTEX gerechnet (Piepho, 2015) und im Fall der A-Optimalität mit einem selbstgeschriebenen Programm in der Programmiersprache Julia (Vo-Thanh & Piepho, 2020). Dieses Resultat ist für uns, wie vermutlich auf für die Leser und Leserinnen auch, auf den ersten Blick sehr überraschend gewesen. Es bedeutet in gewisser Weise eine nachträgliche zusätzliche Rechtfertigung für den systematischen Ansatz der Mitscherlichschen Methode, denn legen wir das besagte lineare geostatistische Modell zugrunde, das wir hier hergeleitet haben, passen Auswertung und Design gewissermaßen optimal zusammen. Allerdings muss einschränkend angemerkt werden, dass das genannte Modell nicht notwendigerweise das für einen Datensatz am besten passende ist, und andere geostatistische Modelle führen zu anderen optimalen Designs (Cullis et al., 2006; Williams & Piepho, 2013). Und die weiterhin wichtigste Methode zur Design-Generierung, nämlich für die randomisierte vollständige Blockanlage mit disjunkten Blöcken, legt kein geostatistisches Modell zugrunde, sondern ein lineares Modell mit Blockeffekten (John & Williams, 1995). Es gibt auch keinen Grund davon auszugehen, dass unter den vielen möglichen geostatistischen Modellen (Schaben­berger & Gotway, 2005) ausgerechnet das lineare Modell, und dann auch noch mit einem Reichweiteparameterwert von genau  = 1/v am besten an die Daten passt. Daher möchten wir aus diesem interessanten und in der Tat überraschenden Ergebnis nicht den Schluss ziehen, dass ein systematisches Design vorzuziehen sei. Unsere Empfehlung ist, trotz dieses Befundes, so weit als möglich den Fisherschen Prinzipien zu folgen (Rundfeldt, 1953; John & Williams, 1995; Bailey, 2008). Im vorliegenden Fall würde dies weiterhin natürlich zur Empfehlung einer randomisierten vollständigen Blockanlage führen. Wir werden aber die Design-Frage im Folgenden wieder aufgreifen, wenn wir die Erweiterung der Mischerlichschen Gleitmethode auf den zweidimensionalen Fall betrachtet haben.

von Boguslawski (1942) hat die Mitscherlich-Methode auf zwei Dimensionen erweitert, wobei er von einer zweifaktoriellen Behandlungsstruktur ausgeht. In seinem Beispiel werden v₁ Stufen des einen Faktors systematisch in den Zeilen eines Blocks angeordnet und die v₂ Stufen des anderen Faktors systematisch in den Spalten eines Blocks. Der gleitende Block hat somit v₁ Zeilen und v₂ Spalten. Die Gesamtzahl der Behandlungen, und daher auch die Gesamtzahl der Parzellen je Block haben die Größe v = v₁ × v₂. Ein bei von Boguslawski (1942) angegebenes Beispiel für zwei Faktoren mit jeweils vier Stufen und sechs Wiederholungen ist in Abb. 3 reproduziert.

Abb. 3. Anlageschema zur Er­tragsflächenmethode von von Boguslawski (1942, Abb. 1). Faktor 1 hat die vier Stufen I-IV (v₁ = 4), Faktor 2 hat die vier Stufen 1–4 (v₂ = 4).

Betrachten wir nun eine zweidimensionale Erweiterung des eindimensionalen Ansatzes für die Mitscherlich-Methode. Hierbei können alle denkbaren gleitenden Blöcke der Dimension v₁ × v₂ berücksichtigt werden. Diese finden wir, wenn oben links in der Anlage beginnend der gleitende Block schrittweise von links nach rechts jeweils um eine Spalte weitergerückt wird und darü­ber hinaus eine schrittweise Verschiebung Zeile für Zeile von oben nach unten erfolgt. Sodann kann für jedes Parzellenpaar in der Anlage abgezählt werden, wie viele gleitende Blöcke beide Parzellen gemeinsam haben. Wird dann wie im eindimensionalen Fall angenommen, dass die Blöcke unabhängige Effekte mit Varianz haben, ergibt sich die folgende Kovarianz-Struktur:

wobei y_rc der Beobachtungswert in der r-ten Zeile und c-ten Spalte ist. Diese Kovarianzstruktur ist in der Geosta­tistik keine übliche, aber sie kann trotzdem mit der REML Methode angepasst werden.

Abb. 4. Beispiel eines Mohn­versuchs (1940) zur Er­tragsflächenmethode von von Boguslawski (1942, Abb. 2) (zweites Beispiel). Die Behand­lungen entsprechen der Struktur in Abb. 3. Kor­nerträge bei 92 % Troc­kensubstanz in kg/Teilstück; Teilstückgrö­ße 4 m × 4 m = 16 qm.

Abb. 5. Anlageschema für ein Design in 8 Zeilen und 12 Spalten mit 16 Be­handlungen. D-effizien­tes Design unter Annahme des linearen geostatistischen Mo­dells (13), angepasst für die Daten aus dem zwei­ten Beispiel (Tab. 5), er­zeugt mit der SAS Prozedur OPTEX (Pie­pho, 2015).

Abb. 6. Anlageschema für ein Design in 8 Zeilen und 12 Spalten mit 16 Be­handlungen. A-effi­­zientes Design unter Annahme des linearen geostatistischen Modells (13), angepasst für die Daten aus dem zweiten Beispiel (Tab. 5), erzeugt mit einem in Julia selbst geschriebe­nen Programm (Vo-Thanh & Piepho, 2020).

Es sei auch hier betont, dass die Designs in Abb. 5 und Abb. 6 ein spezielles Kovarianzmodell voraussetzen und dass es in der Praxis nicht zu erwarten ist, dass genau dieses Modell gut zu den Daten passen wird. Daher wäre dieses Design nicht unbedingt zu empfehlen. Es gibt natür­lich viele Alternativen, von denen eine in Abb. 7 gezeigt ist. Dieses randomisierte Design folgt einem auflösbaren Zeilen-Spalten-Plan mit sechs vollständigen Wiederholungen. Es ist außerdem in den Zeilen latinisiert, was bedeutet, dass keine Behandlung mehr als einmal vorkommt. In diesem Fall ist es sogar so, dass jede Zeile einen vollständigen Block bildet. Außerdem bilden jeweils zwei benachbarte Spalten, also die Spalten 1 und 2, 3 und 4, usw., einen vollständigen Block. Dies wird als t-Latinisierung bezeichnet. Desweiteren ist dieses Design hinsichtlich der räumlichen Nachbarschaft von Behandlungen optimiert. Für die Designoptimierung wird ein linea­res Modell zugrundegelegt, welches neben Behandlungseffekte nur Effekte für die verschiedenen disjunkten Blockeinheiten umfasst, aber keine geostatistische Kovarianzstruktur wie bei den Designs in Abb. 5 und Abb. 6. Nähere Details zu dieser Art von Design finden sich in Piepho et al. (2015a) und Piepho et al. (2020).

Abb. 7. Anlageschema für ein A-optimales Design in 8 Zeilen und 12 Spalten mit 16 Behand­lungen. Auflösbarer Plan mit 6 Wiederholungen, Lati­nisierung in den Zeilen und t-Latinisierung in jeweils t = 2 Spalten, er­zeugt mitCycDesigN 7.0 (VSNi; https://www.vsni.co.uk/soft­ware/cycdesign).

Eine Betrachtung der an die Daten von zwei Beispielen angepassten Kovarianzstrukturen im Rahmen einer Erzeu­gung optimaler Versuchsanlagen führte zu dem überraschenden Befund, dass das systematische Versuchsdesign von Mitscherlich und mit Abstrichen auch das von von Boguslawski optimal unter diesem angenommenen Modell sind. Insofern führt das vorgeschlagene geostatistische Modell zu einer Kohärenz von Versuchsplan und Auswertung dieser Anlagen. Allerdings ist nicht davon auszugehen, dass dieses Modell notwendiger­weise immer optimal an gegebene Daten passt, so dass für zukünf­tige Versuche nicht der Schluss gezogen werden kann, systematische Versuchsdesigns seien grundsätzlich zu empfehlen. Vielmehr ist es aus unserer Sicht empfehlenswert, einen randomisierten Versuchsplan zu verwenden (Rundfeldt, 1953). Einer der großen Vorteile eines randomisierten Versuchsplanes ist, dass die Randomisationsstruktur ein auf jeden Fall zulässiges Analysemodell induziert, welches als Baseline-Modell bezeichnet werden kann (Piepho & Williams, 2010). Ein solches Modell lässt darüber hinaus verschiedene Auswertungsmöglichkeiten offen, einschließlich der Einbeziehung einer geostatistischen Kovarianzstruktur. Likelihood-basierte Kriterien wie die Devianz und das AIC können dann herangezogen werden, um zu entscheiden, ob eine geostatistische Modellierung sinnvoll ist und wenn ja, welche Kovarianzstruktur am besten zu den Daten passt. Einschränkend muss hinzugefügt werden, dass bei Verwendung desselben Datensatzes zur Modellselektion die Signifikanztests einer anschließenden Analyse mit dem ausgewählten Modell das nominelle Signifikanzniveau oft nicht einhalten (Richter et al., 2015).

Für bereits vorliegende Versuche, die nach einem systematischen Design angelegt wurden, liefert das hier vorgeschlagene Modell ebenfalls eine mögliche Auswertungsoption, aber nicht die einzige. Wegen der fehlenden Randomisation kann jedoch kein Baseline-Modell formuliert werden, welches mit einer Randomisation begründet werden könnte. Man könnte versucht sein, das hier angegebene lineare geostatistische Modell als Baseline-Modell zu verstehen. Diesem Eindruck möchten wir aus verschiedenen Gründen entgegentreten. Zum einen könnte durch so eine Betrachtungsweise der falsche Eindruck entstehen, eine systematische Versuchsanlage sei auf jeden Fall vertretbar, weil es ja ein solches Modell gibt. Zum anderen muss berücksichtigt werden, dass ein Paradigmenwechsel erfolgt in dem Moment, wo eine geostatistische Analyse erwogen wird (Piepho et al., 2013). Wie oben betont, induziert eine Randomisation in der Regel ein ganz bestimmtes Analysemodell, welches immer verwendet werden kann. In dem Moment aber, wo eine geostatistische Modellierung erfolgt, wird diese Basis verlassen. Stattdessen wird die Annahme getroffen, die Daten seien nach einem noch zu identifizierenden geostatistischen Modell erzeugt worden. Die große Herausforderung ist dann, das jeweils beste geostatistische Modell zu identifizieren, und zwar aus einer potentiell sehr großen Klasse von Modellen, wobei dann immer die Unsicherheit bleibt, ob tatsächlich das beste Modell gefunden wurde. All diese Komplikationen treten dann nicht auf, wenn ein Versuch randomisiert wurde, weil dann immer das Baseline-Modell zur Verfügung steht. Dieses kann zwar geostatistisch erweitert werden, muss es aber nicht. Bei systematischen Versuchsanlagen besteht dagegen wegen des Fehlens eine Baseline-Modells die Notwendigkeit, verschiedene Modelle auszuprobieren, einschließlich geostatistischer Optionen, und das jeweils beste zu identifizieren. Dabei bleibt immer die Unsicherheit, ob auch das beste Modell gefunden werden konnte und ob es wegen der systematischen Anlage nicht zu verzerrten Schätzungen von Behandlungseffekten kommt.