Entwicklung eines Individual-Tree Basal Area Increment Modells mit einem linearen Mixed-Effects-Ansatz

Zusammenfassung

Modelle mit gemischten Effekten sind flexible und nützliche Werkzeuge zur Analyse von Daten mit einer hierarchischen stochastischen Struktur in der Forstwirtschaft und könnten auch verwendet werden, um die Leistung von Waldwachstumsmodellen deutlich zu verbessern. Hier wird ein Protokoll vorgestellt, das Informationen zu linearen Mixed-Effekt-Modellen synthetisiert.

Zusammenfassung

Hier haben wir ein individuelles Baummodell von 5-Jahres-Basalflächenschritten entwickelt, das auf einem Datensatz basiert, einschließlich 21898 Picea Asperata-Bäumen aus 779 Musterparzellen in der Provinz Xinjiang, Nordwestchina. Um hohe Korrelationen zwischen Beobachtungen derselben Stichprobeneinheit zu verhindern, haben wir das Modell mit einem linearen Mixed-Effekt-Ansatz mit zufälligem Ploteffekt entwickelt, um die stochastische Variabilität zu berücksichtigen. Verschiedene Variablen auf Baum- und Standebene, wie Z. B. Indizes für Baumgröße, Wettbewerb und Standortbedingung, wurden als feste Effekte zur Erklärung der Restvariabilität einbezogen. Darüber hinaus wurden Heteroskedastizität und Autokorrelation durch die Einführung von Varianzfunktionen und Autokorrelationsstrukturen beschrieben. Das optimale modellte lineare Mixed-Effekte wurde durch mehrere Anpassungsstatistiken bestimmt: das Informationskriterium von Akaike, das Bayessche Informationskriterium, die Logarithmuswahrscheinlichkeit und ein Wahrscheinlichkeitsverhältnistest. Die Ergebnisse zeigten, dass signifikante Variablen der Einzelbaum-Basalflächenzunahme die inverse Transformation des Durchmessers in Brusthöhe, die Basalfläche von Bäumen größer als der Themenbaum, die Anzahl der Bäume pro Hektar und die Höhe waren. Darüber hinaus wurden Fehler in der Varianzstruktur am erfolgreichsten durch die exponentielle Funktion modelliert, und die Autokorrelation wurde durch die autoregressive Struktur erster Ordnung (AR(1)) signifikant korrigiert. Die Leistung des linearen Mixed-Effekt-Modells wurde im Vergleich zum Modell mit der Regression der gewöhnlichen kleinsten Quadrate deutlich verbessert.

Einleitung

Im Vergleich zur monogetagten Monokultur hat die walduneinheitliche Waldbewirtschaftung mit mehreren Zielen in jüngster Zeit erhöhte Aufmerksamkeit erhalten^1,2,3. Die Vorhersage verschiedener Bewirtschaftungsalternativen ist für die Formulierung robuster Waldbewirtschaftungsstrategienerforderlich, insbesondere für komplexe, ungleich gealterte Mischarten4 . Waldwachstums- und Ertragsmodelle wurden ausgiebig zur Vorhersage der Baum- oder Standentwicklung und -ernte im Rahmen verschiedener Bewirtschaftungsschemata^5,6,7verwendet. Waldwachstums- und Ertragsmodelle werden in Einzelbaummodelle, Größenklassenmodelle und Vollwertmodell6^,7,8klassifiziert. Leider eignen sich Modelle der Größenklasse und Ganzstandsmodelle nicht für ungleich gealterte Mischartenwälder, die eine detailliertere Beschreibung erfordern, um den Entscheidungsprozess der Waldbewirtschaftung zu unterstützen. Aus diesem Grund haben individuelle Baumwachstum und Ertragsmodelle in den letzten Jahrzehnten erhöhte Aufmerksamkeit erhalten, weil sie in der Lage sind, Vorhersagen für Waldbestände mit einer Vielzahl von Artenzusammensetzungen, Strukturen und Managementstrategien^9,10,11zu machen.

Die Regression der gewöhnlichen kleinsten Quadrate (OLS) ist die am häufigsten verwendete Methode für die Entwicklung von Einzelbaumwachstumsmodellen^12,13,14,15. Die Datensätze für einzelne Baumwachstumsmodelle, die wiederholt über einen festen Zeitraum auf derselben Stichprobeneinheit (d. h. Stichprobendiagramm oder Baum) gesammelt wurden, weisen eine hierarchische stochastische Struktur auf, mit einem Mangel an Unabhängigkeit und hoher räumlicher und zeitlicher Korrelation zwischen den Beobachtungen^10,16. Die hierarchische stochastische Struktur verstößt gegen die grundannahmen der OLS-Regression: nämlich unabhängige Residuen und normal verteilte Daten mit gleichen Varianzen. Daher führt die Verwendung der OLS-Regression unweigerlich zu verzerrten Schätzungen des Standardfehlers von Parameterschätzungen für diese Daten^13,14.

Modelle mit gemischten Effekten bieten ein leistungsstarkes Werkzeug für die Analyse von Daten mit komplexen Strukturen, z. B. Daten mit wiederholten Messeinheiten, Längsschnittdaten und mehrstufigen Daten. Modelle mit gemischten Effekten bestehen aus festen Komponenten, die der gesamten Grundgesamtheit gemeinsam sind, und zufälligen Komponenten, die für jede Stichprobenstufe spezifisch sind. Darüber hinaus berücksichtigen Modelle mit gemischten Effekten Heteroskedastizität und Autokorrelation in Raum und Zeit, indem nicht diagonale Varianz-Kovarianzstrukturmatrizen^17,18,19definiert werden. Aus diesem Grund wurden Mixed-Effekt-Modelle in der Forstwirtschaft ausgiebig eingesetzt, wie z.B. in den Durchmesser-Höhen-Modellen^20,21, Kronenmodelle^22,23, selbstverdünnende Modelle^24,25und Wachstumsmodelle^26,27.

Dabei ging es vor allem darum, ein individuelles Grundflächen-Inkrementmodell mit einem linearen Mixed-Effekt-Ansatz zu entwickeln. Wir hoffen, dass der Ansatz der gemischten Auswirkungen umfassend angewandt werden kann.

Protokoll

1. Datenaufbereitung

1. Vorbereiten von Modellierungsdaten, die Einzelbauminformationen (Art und Durchmesser bei Brusthöhe bei 1,3 m) und Diagramminformationen (Neigung, Aspekt und Höhe) enthalten. In dieser Studie wurden die Daten aus dem 8. (2009) und 9. (2014) Chinesischen National Forest Inventory in der Provinz Xinjiang, Nordwestchina, gewonnen, das 21.898 Beobachtungen von 779 Stichprobenparzellen enthält. Diese Stichprobenparzellen sind quadratisch mit einer Größe von 1 Mu (chinesische Einheit mit einer Fläche von 0,067 ha) und systematisch über ein Raster von 4 km x 8 km angeordnet.
   ANMERKUNG: Die Daten für die Modellierung (Basalfläche) erfordern mindestens eine Wachstumsperiode (d. h. zwei Beobachtungen).
2. Teilen Sie die Daten zufällig in zwei Datasets auf, wobei 80 % der Daten aus den Stichprobendiagrammen für die Modellanpassung verwendet werden (Modellentwicklungsdatensatz), der aus 17.145 Beobachtungen aus 623 Stichprobendiagrammen und 20 % für die Modellvalidierung (Modellvalidierungs-Dataset) besteht, die aus 4.753 Beobachtungen aus 156 Stichprobendiagrammen besteht. Beschreibende Statistiken für die verwendeten Schlüsselvariablen sind in Tabelle 1enthalten.
   HINWEIS: Dieser Schritt der Modellierungsprozedur kann weggelassen werden, und alle Daten werden für die Modellentwicklung verwendet.

Variablen	Anpassen von Daten					Validierungsdaten
	Min	Max	Bedeuten	Sd.		Min	Max	Bedeuten	Sd.
DBH1 (cm)	5	124.8	19.9	13.2		5	101.5	19.5	13.4
QMD (cm)	6.7	82.3	22.5	8.5		9.2	73.3	21.8	9.2
ID (cm)	0.1	14.4	1.1	1		0.1	16.9	1	1.1
BAL (m3)	0	5.2	1.7	0.9		0	5.4	1.7	1
NT (Bäume/ha)	14.9	3642	1072	673.7		14.9	3418	1205	829.3
BA (m2/ha)	0.1	77.5	34.2	13.9		0.1	80.6	34.5	15.3
EL (m)	2	3302	2189	340.3		1441	3380	2256	308.3

Tabelle 1. Zusammenfassungsstatistiken für Anpassungs- und Validierungsdaten. DBH₁: Anfangsdurchmesser bei Brusthöhe bei 1,3 m (DBH), DBH₂: DBH gemessen nach 5 Jahren Wachstum, QMD: quadratischer Mittlerer Durchmesser, ID: Durchmesserzuwachs für 5 Jahre (DBH₂ – DBH₁), BAL: die Grundfläche von Bäumen größer als der Themenbaum (der Subjektbaum: der Baum, der die Wettbewerbsindizes berechnet wurde), NT: die Anzahl der Bäume pro Hektar, BA: Basalfläche pro Hektar, EL: Höhe, S.D.: Standardabweichung.

2. Grundlegende Modellentwicklung

1. Konsultieren Sie Verweise, um Variablen zu identifizieren, die sich auf die Einzelnenbaum-Basalflächeninkremente auswirken.
2. Wählen Sie Variablen basierend auf den Daten aus und berechnen Sie diese. Im Allgemeinen wird das Einzelbaum-Basalflächeninkrement durch drei Variablengruppen beeinflusst: Baumgröße, Wettbewerb und Standortbedingung^27,28,29,30.
   1. Betrachten Sie Baumgrößeneffekte wie DBH₁, Quadrat von DBH₁ ( ), die inverse Transformation von DBH₁ (1/DBH₁), und den gemeinsamen Logarithmus von DBH₁ (logDBH₁) oder Kombinationen von ihnen.
   2. Betrachten Sie Wettbewerbseffekte wie ein- und zweiseitige Wettbewerbsindizes, um das Wettbewerbsniveau eines Baumes sowie seine soziale Position innerhalb des Standes umfassender zu quantifizieren. Zu den einseitigen Wettbewerben gehören BAL und der relative Dichteindex (RD=DBH₁/QMD); Zu den beidseitigen Wettbewerben gehören NT und BA.
      HINWEIS: Die entfernungsabhängigen Wettbewerbsindizes sollten berücksichtigt werden, wenn Daten verfügbar sind.
   3. Betrachten Sie Standorteffekte wie Aspekt (ASP), Neigung (SL) und EL. SL und ASP sollten mithilfe der Stage-Transformation³¹einbezogen werden.
3. Wählen Sie log( - +1) ( bezeichnet Quadrat von DBH₂) als abhängige Variable.
4. Entwickeln Sie das Basismodell mit der schrittweisen Regressionsmethode. Stellen Sie sicher, dass das Modell biologisch vertretbar ist und erhebliche Unterschiede zwischen unabhängigen Variablen aufweist. Verwenden Sie den Varianzinflationsfaktor (VIF), um die Multikollinearität zu überprüfen.
5. Lassen Sie die unabhängigen Variablen mit p < 0.05 und VIF < 5 im Basismodell.
6. Geben Sie die Basismodellergebnisse und das Restdiagramm aus. Das hier produzierte Basismodell dient als Grundlage für die Weiterentwicklung eines Mixed-Effekt-Modells.

3. Lineare Mixed-Effekt-Modellentwicklung mit dem Paket "nlme" in R-Software

1. Lesen Sie das Modellentwicklungs-Dataset, und laden Sie das Paket "nlme".
   >model.development.dataset=read.csv("E:/DATA/JoVE/modelingdata.csv",
   header=TRUE)
   >Bibliothek(nlme)
2. Wählen Sie Stichprobendiagramme als Zufallseffekte aus, um das Modell für gemischte Effekte zu entwickeln.
3. Passen Sie alle möglichen Kombinationen von Zufallseffekten mit der Methode der maximalen Wahrscheinlichkeit (ML) an und geben Sie die Ergebnisse aus.
   >Modell<-lme(Y1/DBH₁+BAL+NT+EL,data=model.development.dataset,
   method="ML", zufällig = 1| PLOT)
   >Zusammenfassung(Modell)
   1. Set zufällig = 1 ist der Abfang zu zufälligen Parametern. Ändern Sie die zufälligen Anweisungen, bis alle Kombinationen angepasst sind. Um beispielsweise 1/DBH₁ und BAL als zufällige Parameter festzulegen, lautet der Code wie folgt: zufällig = 1/DBH₁+BAL-1. Darüber hinaus können die Codes bei der Anpassung Fehler melden, die auf die Nichtkonvergenz des angepassten Modells zurückzuführen sind.
4. Wählen Sie das beste Modell nach dem Informationskriterium (AIC) von Akaike, dem Bayesschen Informationskriterium (BIC), der Logarithmuswahrscheinlichkeit (Loglik) und dem Wahrscheinlichkeitsverhältnistest (LRT).
   >anova(Modell.1, Modell.6)
   >anova(Modell.6, Modell.23)
   >anova(Modell.23, Modell.30)
5. Bestimmen Sie die Struktur von R_i. Adressieren Sie die Heteroskedastizität und Autokorrelation von R_i³². Das R_i ist wie folgt geschrieben:
   (1)
   wobei σ² ein unbekannter Skalierungsfaktor ist, der der Modell-Restvarianz entspricht, ist G_i eine diagonale Matrix, die Heteroskedastizität beschreibt, und Γ_i ist eine Matrix, die die Autokorrelation beschreibt.
   1. Beobachten Sie, ob die Residuen heteroskedasitizität aus dem Restdiagramm haben. Wenn heteroskedastistizistiz ist (die Residuen haben ein klares Muster oder einen klaren Trend), führen Sie drei häufig verwendete Varianzfunktionen ein – die Konstante-Plus-Leistungsfunktion, die Leistungsfunktion und die exponentielle Funktion –, um die Fehlervarianzstruktur zu modellieren.
      >Modell.30.1<-lme(Y'1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",random='1/DBH₁+BAL+NT| Plot
      weights=varConstPower(form='angepasst(.)))
      >Zusammenfassung(Modell.30.1)
      >Modell.30.2<-lme(Y'1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",random='1/DBH₁+BAL+NT| Plot
      gewichte=varPower(form='angepasst(.)))
      >Zusammenfassung(Modell.30.2)
      >Modell.30.3<-lme(Y'1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",random='1/DBH₁+BAL+NT| Plot
      weights=varExp(form='angepasst(.)))
      >Zusammenfassung(Modell.30.3)
   2. Bestimmen Sie die beste Varianzfunktion für das Modell gemäß AIC, BIC, Loglik und LRT.
      >anova(Modell.30, Modell.30.1)
      >anova (Modell.30, Modell.30.2)
      >anova(Modell.30, Modell.30.3)
   3. Führen Sie drei häufig verwendete Autokorrelationsstrukturen ein – die zusammengesetzte Symmetriestruktur (CS), die autoregressive Struktur erster Ordnung [AR(1)] und eine Kombination aus autoregressiven und gleitenden Durchschnittsstrukturen erster Ordnung [ARMA(1,1)]–, um die Autokorrelation zu berücksichtigen.
      >Modell.30.3.1<-lme(Y'1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",
      zufällig = 1 /DBH₁+BAL+NT| PLOT, weights=varExp(form='fitted(.)), corr= corCompSymm())
      >Zusammenfassung(Modell.30.3.1)
      >Modell.30.3.2<-lme(Y'1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",
      zufällig = 1 /DBH₁+BAL+NT| PLOT,weights=varExp(form='angepasst(.)), corr=corAR1())
      >Zusammenfassung(Modell.30.3.2)
      >Modell.30.3.3<-lme(Y1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="ML",
      zufällig = 1 /DBH₁+BAL+NT| PLOT,weights=varExp(form='angepasst(.)), corr=corARMA(q=1,p=1))
      >Zusammenfassung(Modell.30.3.3)
   4. Bestimmen Sie die beste Autokorrelationsstruktur gemäß AIC, BIC, Loglik und LRT.
      >anova(Modell.30.3, Modell.30.3.2)
      ANMERKUNG: Die G_i und i können nicht definiert werden, wenn es keine Heteroskedastizität und Autokorrelation gibt.
   5. Geben Sie die Endergebnisse des Mixed-Effekt-Modells mit der REML-Methode (Restricted Maximum Likelihood) aus.
      >Mixed.model<-lme(Y'1/DBH₁+BAL+NT+EL,data=model.development.dataset, method="REML",random='1/DBH₁+BAL+NT| Plot
      weights=varExp(form=' fitted(.)), corr=corAR1()))
      >Zusammenfassung(Mixed.model)

4. Bias-Korrektur

1. Transformieren Sie die vorhergesagten Werte des Basalflächeninkrements mithilfe des endgültigen Modells auf einer logarithmischen Skala in den ursprünglichen Maßstab. Eine solche lineare Rücktransformation des vorhergesagten Wertes aus einem logtransformierten Modell erzeugt jedoch eine zugeordnete Log-Transformationsverzerrung. Um die Log-Bias zu behandeln, wurde ein Korrekturfaktor abgeleitet und in die Vorhersagegleichung integriert, die das tatsächliche vorhergesagte Grundflächeninkrement für einen bestimmten Baum schätzt [Gleichung (2)]:
   (2)
   wobei der prognostizierte logarithmische Wert des Basalflächeninkrements aus dem Modell vorhergesagt wird, während der vorhergesagte rücktransformierte Wert des Basalflächeninkrements für 5 Jahre nach Korrektur für Log-Transformations-Bias ist. ist die Varianz von Zufälligeffekten im Diagramm und σ² ist Restvarianz.
2. Konvertieren Sie das Grundflächeninkrement ( ) in das Durchmesserinkrement.

5. Modellvorhersage und -auswertung

1. Bereiten Sie das in Abschnitt 1.2 erstellte Modellvalidierungs-Dataset für die Vorhersage vor.
2. Verwenden Sie das lineare Modell für gemischte Effekte, um das Inkrement der einzelnen Baumbasalflächen vorherzusagen. Die zufälligen Komponenten wurden mit dem folgenden besten linearen unvoreingenommenen Prädiktor berechnet:
   (3)
   wobei ein Vektor für die zufälligen Komponenten ist; ist die Varianz-Kovarianz-Matrix für die Variabilität zwischen den Plots; ist die Konstruktionsmatrix für die zufälligen Komponenten, die auf die komplementären Beobachtungen wirken; ist der Restvektor, dessen Komponenten durch die Differenz zwischen den Basalflächenschritten und den vorhergesagten Inkrementen mit dem Modell der festen Effekte angegeben werden.
3. Bewerten und vergleichen Sie die Vorhersagefähigkeit des Basismodells und des linearen Mixed-Effekt-Modells anhand der folgenden drei statistischen Indikatoren^23,33.
   (4)
   (5)
   (6)
   wobei obj_i die Basalflächeninkremente ist, est_i die vorhergesagten Basalflächenschritte, das Mittel der Beobachtungen und N die Anzahl der Beobachtungen.

Ergebnisse

Das grundlegende Grundflächeninkrementmodell für P. asperata wurde als Gleichung (7) ausgedrückt. Die Parameterschätzungen, die entsprechenden Standardfehler und die Nichtanpassungsstatistiken sind in Tabelle 2dargestellt. Das Restdiagramm ist in Abbildung 1dargestellt. Es wurde eine ausgeprägte Heteroskedastizität der Residuen beobachtet.
(7)

...

Diskussion

Eine entscheidende Frage für die Entwicklung von Modellen mit gemischten Effekten ist die Bestimmung, welche Parameter als Zufallseffekte behandelt werden können und welche als feste Effekte^34,35betrachtet werden sollten. Es wurden zwei Methoden vorgeschlagen. Der häufigste Ansatz besteht darin, alle Parameter als Zufallseffekte zu behandeln und dann das beste Modell von AIC, BIC, Loglik und LRT auswählen zu lassen. Dies war die Methode unserer Studie

Materialien

Name	Company	Catalog Number	Comments
Computer	acer
Microsoft Office 2013
R x64 3.5.1