Wahl des statistischen Verfahrens

[Verena_St.]

Hallo,

ich habe im Rahmen eines Projekts Daten (Treibhausgasemissionen, sowie Bodenparameter (Wassergehalt, pH)) auf 2 verschiedenen Feldern (unterschiedliche Bewirtschaftungsform) erhoben.
Jedes Feld wurde in 4 Zonen untergliedert, mit 10 Messpunkten pro Zone.
Die Zonen ergaben sich nach der Entfernung zur Futterstelle, bzw. eine Zone bildet der Wald. Anzunehmen ist, dass in der Nähe der Futterstelle bzw. im Wald höhere Verweildauer der Tiere (und damit höhere THG Emissionen), als auf den weiter entfernten Flächen.

Mich würde nun interessieren, wie z.B. die THG mit dem pH/dem Wassergehalt/ (v.a.) mit der Bewirtschaftungsform/dem Abschnitt auf dem Feld zusammenhängen.

Ich habe nun im ersten Schritt die Mittelwerte, min/max pro Feld, ... verglichen und Grafiken in R erstellt.

Meine Frage wäre nun, welche weiteren statistischen Analysen/Schritte (in R) sich hierzu anbieten würden, um die Zusammenhänge zwischen Emissionen und Bodenparametern, Managementform bzw. räumlicher Lage der Messpunkte zu erkennen?

Frage des Versuchs ist allgemein: Welche Auswirkung die Managementform auf die THG-Emissionen hat?

Ich hoffe meine Frage ist nachvollziehbar beschrieben und wäre sehr dankbar über jeden Input.

Mit lieben Grüssen und vielen Dank im Voraus
Verena

[bigben]

Hallo Verena,

Deine Daten gehören innerhalb der Bewirtschaftungsform und da innerhalb der Zonen jeweils in Gruppen zusammen. Deshalb dürfte es sich anbieten, sie in einem sog. Hierarchischen Modell, beispielsweise einem mixed effects Modell, zu untersuchen.
Das kommt in der Ausbildung halt viel später als Minimum und Maximum. Was ist denn Dein Ausbildungsstand? Weißt Du, was multiple lineare Regression ist? Was fixed und was random effects sind?

Kannst Du das mit den Zonen nochmal erklären? Ist Zone 1 auf Feld 1 ähnlich wie Zone 1 auf Feld 2 oder ist die Zoneneinteilung auf beiden Feldern unabhängig?

LG,
Bernhard

[Verena_St.]

Lieber Bernhard,

vielen Dank für die rasche Rückmeldung und die konkreten Vorschläge - das hilft bereits sehr!

Zu den Zonen:
Genau, die Zonen/strata sind auf beiden Feldern ident (nur die Flächen in den einzelnen Zonen sind unterschiedlich groß, da beide Felder unterschiedlich groß sind).
Die Unterteilung sieht auf beiden Feldern wie folgt aus:
Zone 1: 100 m um die Futterstelle
Zone 2: 100 - 500 m um die Futterstelle
Zone 3: > 500 m von der Futterstelle entfernt
Zone 4: Waldflächen auf dem Feld.

Ich bin im Master (nun Abschlussarbeit) und hatte im Bachelor und Master jeweils eine Statistik-VO. Lineare Regression (musste ich zwar noch kurz nachlesen), haben wir tatsächlich gemacht.
Für mein Beispiel würde mir die multiple lineare Regression den Einfluss der Bodenparameter auf meine THG zeigen, wenn ich das so richtig verstanden habe?

Fixed und random effects kommen mir nicht wirklich bekannt vor. Habe eine Beschreibung und ein Beispiel kurz nachgelesen. Grundsätzlich ist die Beschreibung nachvollziehbar, aber die Anwendung auf mein Beispiel und das Ziel des mixed effect Modells fallen mir noch etwas schwer.

D.h. ich werde nun einmal die multiple lineare Regression pro Feld durchführen und noch das mixed effect Modell genauer nachlesen.

Liebe Grüße
Verena

[bigben]

...ich werde nun einmal die multiple lineare Regression pro Feld durchführen und noch das mixed effect Modell genauer nachlesen.

Hallo Verena,

das lineare Modell muss nicht unbedingt das beste sein, aber es ist ein gutes Standardverfahren, dass man gerne mal als erstes ausprobiert. Wenn Du nur das im Studium gelernt hast, umso mehr. Mixed effects musst Du erstmal nicht nachlesen, ich denke, für Deine Zwecke kommst Du ohne aus.

Du rechnest am besten nicht eine lineare Regression pro Feld sondern eine mit allen Daten und dummycodiertem Feld.

Schau Dir dieses Beispiel mal an:

Code: Alles auswählen

boden <- data.frame(
  feld = rep(c("A", "B"), each = 12),
  zone = rep(c("nah", "mittel", "fern"), 8),
  pH = c(8.8, 7.5, 6, 8.5, 7.7, 6.7, 5.8, 6.5, 5.9, 7.3, 6.4, 8, 7.5, 
         7.6, 7.2, 6.9, 6.3, 6.7, 5.6, 8.4, 6.3, 8, 7.1, 8.7),
  thg = c(2.17, 2.73, 1.97, 2.16, 2.74, 1.97, 2.22, 2.67, 1.96, 2.17, 
          2.7, 2.03, 4.24, 4.74, 4.04, 4.23, 4.73, 3.98, 4.19, 4.69, 3.98, 
          4.18, 4.7, 3.98)
  )
print(daten)

modell_1 <- lm(thg ~ feld + zone + pH, data = boden)
summary(modell_1)  |>  summary()

qqnorm(modell_1$residuals, main = "QQ-Plot der Residuen")
qqline(modell_1$residuals)

plot(modell_1$fitted.values, modell_1$residuals,
      xlab = "Fitted Values | Vorhergesagt",
      ylab = "Residuals | Residuen")

plot(1:5, modell_1$coefficients, xaxt = "n", pch = 15)
axis(1, at = 1:5, labels = names(modell_1$coefficients))

LG,
Bernhard

[Verena_St.]

Liebe Bernhard,

vielen Dank für die weiteren Ergänzungen und Erläuterungen!
Ich habe nun versucht es mit meinen Daten durchzuführen, folgende Fragen habe ich dazu bitte noch:

1.) Wie erhält man z.b. bei pH die Werte in einer Zeile? Ich habe Sie aus meiner Excel Datei kopiert, dann werden Sie allerdings als Spalte in R eingefügt.
2.) Worauf bezieht sich "Daten" in "print(daten)"?
3.) Wofür wird der Teil "|>summary()" gebraucht in summary(modell_1) |> summary()
4.) Folgenden Befehl verstehe ich leider auch nicht (Was ich daraus herauslesen kann? Woraus sich 1:5 ergibt? :
plot(1:5, modell_1$coefficients, xaxt = "n", pch = 15)
axis(1, at = 1:5, labels = names(modell_1$coefficients))

Vielen Dank nochmals und lG
Verena

[bigben]

Hallo Verena,

2.) Worauf bezieht sich "Daten" in "print(daten)"?

Das war ein Fehler. Ich wollte dort print(boden) schreiben, da ich meine ausgedachten Daten davor ja boden genannt hatte. Abends um zehn kommen solche Fehler schon mal vor.

1.) Wie erhält man z.b. bei pH die Werte in einer Zeile? Ich habe Sie aus meiner Excel Datei kopiert, dann werden Sie allerdings als Spalte in R eingefügt.

Am besten gar nicht. Du solltest die Daten in Excel als ordentliche Tabelle speichern, wobei jede Spalte für einen Messwert steht und jede Zeile für eine Messung. Dann speicherst Du die Daten in Excel als csv-Datei und liest sie dann mit read.csv2 ein, oder Du speicherst sie als Excel-Tabelle und liest sie mit Hilfe von Zusatzpaketen wie rio oder openxlsx oder readxl als R Dataframe ein. Du kannst sie in RStudio auch über den Menüpunkt "File" und dort "Import Dataset" einlesen. Es macht auf jeden Fall enorm viel Sinn, im Internet ein wenig hin und her zu lesen, wie man Excel-Daten in R einliest. Copypaste aus Excel sollte dabei nicht der Normalfall sein.

Google und Youtube sind deine Freunde, z.B. http://www.sthda.com/english/wiki/readi ... lsx-into-r oder https://www.youtube.com/watch?v=46yqTMYW7U0 oder https://www.youtube.com/watch?v=uEYEBIkQAs8 oder https://www.youtube.com/watch?v=fPjtqSjYDMs oder https://www.youtube.com/watch?v=JYVWufSQ4OI oder https://www.spsanderson.com/steveondata ... 024-06-14/ oder...

3.) Wofür wird der Teil "|>summary()" gebraucht in summary(modell_1) |> summary()

Auch das ein Tippfehler. Die Funktion summary beschreibt das Ergebnis der Regression. Diese beiden Zeilen bedeuten das gleiche, das was ich geschrieben habe war Müll (ich wurde beim coden von einem Telefonat gestört):

Code: Alles auswählen

summary(modell_1)
modell_1 |> summary()

4.) Folgenden Befehl verstehe ich leider auch nicht (Was ich daraus herauslesen kann? Woraus sich 1:5 ergibt? :
plot(1:5, modell_1$coefficients, xaxt = "n", pch = 15)
axis(1, at = 1:5, labels = names(modell_1$coefficients))

So wie ich das Modell aufgestellt habe erzeugt es fünf Koeffizienten. Um die grafisch darzustellen habe ich auf der x-Achse die Werte 1 bis 5 und auf der y-Achse die Koeffizienten dargestellt. 1:5 war nur dazu da, die x-Werte 1 bis 5 zu erzeugen.
Du musst das nicht nachmachen. Einfacher, schöner und besser kann man so ein Diagramm beispielsweise mit dem Zusatzpaket dotwhisker machen. Wenn man das von CRAN herunterlädt, dann kann man daraus die Funktion dwplot machen, um ein Modell oder auch mehrere im Vergleich darzustellen. Bei Interesse findest Du mehr dazu unter https://cran.r-project.org/web/packages ... nette.html

Eine Überarbeitete Version von meinem Code, diesmal ausgeschlafen, könnte dann beispielsweise so aussehen:

Code: Alles auswählen

install.packages("dotwhisker")   # diese Zeile braucht man nur beim ersten Mal

boden <- data.frame(
  feld = rep(c("A", "B"), each = 12),
  zone = rep(c("nah", "mittel", "fern"), 8),
  pH = c(8.8, 7.5, 6, 8.5, 7.7, 6.7, 5.8, 6.5, 5.9, 7.3, 6.4, 8, 7.5, 
         7.6, 7.2, 6.9, 6.3, 6.7, 5.6, 8.4, 6.3, 8, 7.1, 8.7),
  thg = c(2.17, 2.73, 1.97, 2.16, 2.74, 1.97, 2.22, 2.67, 1.96, 2.17, 
          2.7, 2.03, 4.24, 4.74, 4.04, 3.23, 4.73, 3.98, 4.19, 4.69, 3.98, 
          4.18, 4.7, 3.98)
)
print(boden)

modell_1 <- lm(thg ~ feld + zone + pH, data = boden)
summary(modell_1)  

# Visualisierung der gefundenen Koeffizienten. Wenn hier Werte zu fehlen
# scheinen, dann nochmal die Grundlagen der Dummycodierung studieren.
dotwhisker::dwplot(modell_1)

# Die folgende Abbildung kann dabei helfen zu beurteilen, ob eine multiple
# lineare Regression in dieser Form zu den Daten passt, oder nicht. 
# Mehr dazu unter 
# http://www.sthda.com/english/articles/39-regression-model-diagnostics/161-linear-regression-assumptions-and-diagnostics-in-r-essentials/#diagnostic-plots
par(mfrow = c(2, 2))
plot(modell_1)

LG,
Bernhard

[Verena_St.]

Perfekt, vielen Dank für die Hilfe! Ich habe nun die einfache und multiple lineare Regression mit H2O/pH/Zone/Feld (letztere zwei dummykodiert) durchgeführt und die Grafiken erstellt.
So wie ich das verstehe, mit dem Ergebnis, dass lediglich die Zone signifikanten Einfluss hat.

Besten Dank für den Input und die Unterstützung!
LG:)
Verena

[bigben]

So wie ich das verstehe, mit dem Ergebnis, dass lediglich die Zone signifikanten Einfluss hat.

Erstmal herzlichen Glückwunsch zur gelungenen Analyse. Wenn Du unsicher bist ("So wie ich das verstehe") bist Du gerne herzlich eingeladen, die summary und die diagnostischen Plots hier hochzuladen.

VG,
Bernhard

[Verena_St.]

Ja, vielen Dank, das nehme ich dann gerne in Anspruch:)

Folgendes Ergebnis habe ich bekommen:

Code: Alles auswählen

lm(formula = N2O ~ H2O + pH + DRX + Db + Dc + Dd, data = R_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.24821 -0.10432 -0.04552  0.04249  0.44983 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.91055    1.10944   0.821    0.428
H2O          0.01860    0.02144   0.867    0.403
pH          -0.12366    0.17966  -0.688    0.504
DRX         -0.02317    0.13738  -0.169    0.869
Db           0.10212    0.17302   0.590    0.566
Dc          -0.17819    0.15023  -1.186    0.259
Dd          -0.18844    0.22669  -0.831    0.422

Residual standard error: 0.2163 on 12 degrees of freedom
  (43 observations deleted due to missingness)
Multiple R-squared:  0.2021,	Adjusted R-squared:  -0.1969 
F-statistic: 0.5065 on 6 and 12 DF,  p-value: 0.7924

DRX...dummykodiertes Feld Nummer 2
Dc/Dd/Db...dummykodiertes Zonen

Da mein p-Wert >0,05 ist, hätte ich gesagt, dass der Einfluss aller unabhängigen Variablen nicht signifikant ist, d.h. die Emissionen nicht wirklich durch pH/H2O usw. erklärt werden können (Was für mich kein tolles Ergebnis ist?). Ich habe zuvor eine einfache Regression für die einzelnen Parameter (pH, H2O,...)durchgeführt, hier war dann lediglich der p-Wert bei den Zonen signifikant (zone a als Referenz). Somit hätte ich gesagt, das heißt, dass durch die Zonen die N2O Emissionen erklärt werden können (in Zone a höhere).


Zwischenzeitlich habe ich ergänzend dazu einem Beispiel zu gstat gefolgt (https://cran.r-project.org/web/packages ... /gstat.pdf) und für meine Daten umgesetzt (THG auf Karte dargestellt mit Bubble plot).

Nun würde ich gerne ein Variogramm (Kap. 5) erstellen (als Vorarbeit zu kriging). Toll wäre nämlich eine Darstellung, in der man den Verlauf der THG über die beiden Felder sieht. Allerdings ist mir hier noch ein absolutes Rätsel, welche Spalten/Daten ich hierzu (für variogram-Befehl) in meiner Excel brauche.
Hast du hierzu evt. auch noch einen Tip?

In weiterer Folge würde ich dann vielleicht noch ein Mixed-effect Modell probieren.

Ich hoffe, mein Thema ist halbwegs nachvollziehbar beschrieben!
Vielen Dank nochmals im Voraus und lG
Verena

[bigben]

Code: Alles auswählen

43 observations deleted due to missingness

Du hast viel zu wenige Beobachtungen für so viele Prädiktoren! Wieso sind denn die Daten so unvollständig, dass da 2/3 fehlen?
Oben hattest Du 2x4x10 = 80 Beobachtungen in Aussicht gestellt, jetzt sind es keine zwanzig mehr und unvollständige Beobachtungen kann eine lineare Regression nicht verarbeiten..

Code: Alles auswählen

summary(R_data)

Könnte helfen herauszufinden, ob da immer die gleiche Variable fehlt. Dann könnte man die vielleicht aus dem Modell nehmen, wenn es nicht gerade N20 ist.

LG,
Bernhard