library(readr)
library(ggplot2)Statistik 4: Übung
- Lesen Sie den Datensatz steprasen_ukraine.csv in R ein. Dieser enthält Pflanzenartenzahlen (Species_richness) von 199 10 m² grossen Plots (Vegetationsaufnahmen) von Steppenrasen in der Ukraine sowie zahlreiche Umweltvariablen (Table 1)
| Name | Description |
|---|---|
| Elevation | Elevation [m a.s.l.] |
| Inclination | Inclination [°] |
| Heat_index | Heat index |
| Microrelief | Microrelief [cm] |
| Grazing_intensity | Grazing intensity (0 = 0; 1 = low; 2 = medium; 3 = high) |
| Litter | Cover litter [%] |
| Stones_and_rocks | Cover stones and rocks [%] |
| Gravel | Cover gravel [%] |
| pH | pH (H2O) |
| Conductivity | Conductivity (µS/cm) |
| CaCO3 | CaCO3 content (%) (g CaCO3 / g soil) |
| N_total | N total (mmol/gTB) |
| C_org | C org (mmol/gTB) |
| CN_ratio | C org / N org (mol/mol) |
| Temperature | Annual mean temp. °C x 10 (WorldClim, BIO_01) |
| Temperature_range | Annual temp. range °C x 10 (WorldClim, BIO_07) |
| Precipitation | Annual precipitation mm (WorldClim, BIO_12) |
- Ermitteln Sie ein minimal adäquates Modell, das den Artenreichtum in den Plots durch die Umweltvariablen erklärt.
- Bitte erklären und begründen Sie die einzelnen Schritte, die Sie unternehmen, um zu diesem Ergebnis zu kommen. Dazu erstellen Sie bitte ein mit Quarto generiertes pdf-Dokument, in das Sie Schritt für Schritt den verwendeten R-Code, die dazu gehörigen Ausgaben von R, Ihre Interpretation derselben und die sich ergebenden Schlussfolgerungen für das weitere Vorgehen dokumentieren.
- Dieser Ablauf sollte insbesondere beinhalten:
Überprüfen der Datenstruktur nach dem Einlesen: welches sind die abhängige(n) und welches die unabängige(n) Variablen, sind alle Variablen für die Analyse geeignet?
Explorative Datenanalyse
Definition eines globalen Modelles und dessen Reduktion zu einem minimal adäquaten Modell
Durchführen der Modelldiagnostik für dieses
Generieren aller Zahlen, Statistiken und Tabellen, die für eine wiss. Ergebnisdarstellung benötigt werden
Formulieren Sie abschliessend einen Methoden- und Ergebnisteil (ggf. incl. adäquaten Abbildungen) zu dieser Untersuchung in der Form einer wissenschaftlichen Arbeit (ausformulierte schriftliche Zusammenfassung, mit je einem Absatz von ca. 60-100 Worten, resp. 3-8 Sätzen für den Methoden- und Ergebnisteil). D. h. alle wichtigen Informationen sollten enthalten sein, unnötige Redundanz dagegen vermieden werden.
- Zu erstellen sind (1) Ein quarto generiertes pdf-Dokument mit begründetem Lösungsweg (Kombination aus R-Code, R Output und dessen Interpretation) und (2) ausformulierter Methoden- und Ergebnisteil (für eine wiss. Arbeit).
Lösung Übung 4
Demoscript herunterladen (.qmd)
stgl <- read_delim("datasets/stat/steprasen_ukraine.csv", delim = ";")
str(stgl)spc_tbl_ [199 × 20] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ PlotID : chr [1:199] "UA01NW" "UA01SE" "UA02NW" "UA02SE" ...
$ Species_richness : num [1:199] 44 53 48 50 53 40 46 56 30 35 ...
$ Altitude : num [1:199] 179 178 188 183 162 165 153 158 192 197 ...
$ Inclination : num [1:199] 24 17 27 33 7 33 30 32 25 18 ...
$ Heat_index : num [1:199] -0.42 -0.3 -0.51 -0.65 -0.09 -0.42 0 -0.59 0.46 0.32 ...
$ Microrelief : num [1:199] 5 2.5 2 2 3 4 16 15 5 3 ...
$ Grazing_intensity: num [1:199] 0 0 0 0 0 0 1 1 0 0 ...
$ Litter : num [1:199] 12 10 0 4 15 30 5 6 10 20 ...
$ Stones_and_rocks : num [1:199] 0 0 0 0 0 0 40 10 0 0 ...
$ Gravel : num [1:199] 0 0 0 0 0 0 0 0 0 0 ...
$ Fine_soil : num [1:199] 2 5 0 7 0 0 2 5 5 2 ...
$ pH : num [1:199] 7.32 6.91 6.72 6.44 6.1 6.23 6.79 6.43 7.19 7 ...
$ Conductivity : num [1:199] 90 115 126 90 73 76 163 119 151 69 ...
$ CaCO3 : num [1:199] 0.0754 0.1271 0.0723 0.0771 0.0829 ...
$ N_total : num [1:199] 0.14 0.17 0.24 0.26 0.29 0.2 0.34 0.29 0.18 0.2 ...
$ C_org : num [1:199] 1.54 1.97 2.99 3.22 3.77 2.5 4.59 3.67 2.16 2.38 ...
$ CN_ratio : num [1:199] 10.8 11.5 12.5 12.6 12.8 ...
$ Temperature : num [1:199] 79 79 80 80 82 82 82 82 79 83 ...
$ Temperature_range: num [1:199] 330 330 329 329 328 328 328 328 326 327 ...
$ Precipitation : num [1:199] 608 608 603 603 594 594 594 594 600 586 ...
- attr(*, "spec")=
.. cols(
.. PlotID = col_character(),
.. Species_richness = col_double(),
.. Altitude = col_double(),
.. Inclination = col_double(),
.. Heat_index = col_double(),
.. Microrelief = col_double(),
.. Grazing_intensity = col_double(),
.. Litter = col_double(),
.. Stones_and_rocks = col_double(),
.. Gravel = col_double(),
.. Fine_soil = col_double(),
.. pH = col_double(),
.. Conductivity = col_double(),
.. CaCO3 = col_double(),
.. N_total = col_double(),
.. C_org = col_double(),
.. CN_ratio = col_double(),
.. Temperature = col_double(),
.. Temperature_range = col_double(),
.. Precipitation = col_double()
.. )
- attr(*, "problems")=<externalptr> Man erkennt, dass alle Spalten bis auf die erste mit der Plot ID numerisch (num) sind. Wir können nun die Variable “PlotID” als rownames setzten
# Setze erste Spalte (colum) als rownames
library(tibble)
stgl <- stgl |>
column_to_rownames(var = "PlotID")Nun steht die abhängige Variable in Spalte 1 und die Prediktorvariablen in den Spalten 2 bis 19 stehen.
summary(stgl) Species_richness Altitude Inclination Heat_index
Min. :14.00 Min. : 73.0 Min. : 1.00 Min. :-0.94000
1st Qu.:34.00 1st Qu.:140.0 1st Qu.:12.00 1st Qu.:-0.15500
Median :40.00 Median :166.0 Median :19.00 Median : 0.01000
Mean :40.23 Mean :161.7 Mean :19.28 Mean : 0.01603
3rd Qu.:47.50 3rd Qu.:188.0 3rd Qu.:25.00 3rd Qu.: 0.21500
Max. :67.00 Max. :251.0 Max. :48.00 Max. : 0.85000
Microrelief Grazing_intensity Litter Stones_and_rocks
Min. : 0.000 Min. :0.0000 Min. : 0.00 Min. : 0.000
1st Qu.: 2.500 1st Qu.:0.0000 1st Qu.: 3.50 1st Qu.: 0.000
Median : 5.000 Median :1.0000 Median : 7.00 Median : 0.500
Mean : 7.126 Mean :0.9296 Mean :12.16 Mean : 3.994
3rd Qu.: 7.000 3rd Qu.:2.0000 3rd Qu.:13.50 3rd Qu.: 4.000
Max. :100.000 Max. :3.0000 Max. :90.00 Max. :68.000
Gravel Fine_soil pH Conductivity
Min. : 0.000 Min. : 0.00 Min. :4.890 Min. : 40.0
1st Qu.: 0.000 1st Qu.: 2.00 1st Qu.:7.240 1st Qu.:148.5
Median : 0.000 Median : 5.00 Median :7.420 Median :171.0
Mean : 2.984 Mean : 7.02 Mean :7.286 Mean :162.3
3rd Qu.: 3.000 3rd Qu.:10.00 3rd Qu.:7.545 3rd Qu.:189.5
Max. :40.000 Max. :38.00 Max. :7.790 Max. :232.0
CaCO3 N_total C_org CN_ratio
Min. : 0.0042 Min. :0.0700 Min. : 1.040 Min. : 6.04
1st Qu.: 0.4306 1st Qu.:0.2000 1st Qu.: 2.850 1st Qu.:12.24
Median : 4.6578 Median :0.2700 Median : 3.560 Median :12.95
Mean : 7.4757 Mean :0.2788 Mean : 3.689 Mean :13.48
3rd Qu.:13.0002 3rd Qu.:0.3300 3rd Qu.: 4.400 3rd Qu.:14.02
Max. :35.2992 Max. :0.9500 Max. :11.300 Max. :27.42
Temperature Temperature_range Precipitation
Min. :78.00 Min. :326.0 Min. :577.0
1st Qu.:82.00 1st Qu.:328.0 1st Qu.:583.0
Median :84.00 Median :329.0 Median :592.0
Mean :84.82 Mean :328.6 Mean :596.4
3rd Qu.:88.00 3rd Qu.:330.0 3rd Qu.:602.5
Max. :92.00 Max. :331.0 Max. :630.0 Da es sich bei Artenzahlen um Zähldaten handelt, müsste man theoretisch ein glm mit Poisson-Verteilung rechnen; bei einem Mittelwert, der hinreichend von Null verschieden ist (hier: ca. 40), ist eine Poisson-Verteilung aber praktisch nicht von einer Normalverteilung zu unterscheiden und wir können uns den Aufwand auch sparen.
cor <- cor(stgl[, 2:19])
cor[abs(cor)<0.7] <- 0
corDie Korrelationsanalyse dient dazu, zu entscheiden, ob die Prädiktorvariablen hinreichend voneinander unabhängig sind, um alle in das globale Modell hinein zu nehmen. Bei Pearson’s Korrelationskoeffizienten r, die betragsmässig grösser als 0.7 sind, würde es problematisch. Alternativ hätten wir auch den VIF (Variance Inflation Factor) als Kriterium für den möglichen Ausschluss von Variablen aus dem globalen Modell nehmen können.
Wenn man auf die matirx “cor” nun in einem separaten Fenster öffnet, sieht man, wo es problematische Korrelationen zwischen Variablenpaaren gibt. Es sind dies Altitude vs. Temperature und N.total vs. C.org. Wir müssen aus jedem dieser Paare jetzt eine Variable rauswerfen, am besten jene, die weniger gut interpretierbar ist. Ich entscheide mich dafür Temperature statt Altitude (weil das der direktere ökologische Wirkfaktor ist) und C.org statt N.total zu behalten (weil es in der Literatur mehr Daten zum Humusgehalt als zum N-Gehalt gibt, damit eine bessere Vergleichbarkeit erzielt wird). Die Aussagen, die wir für die beibehaltene Variable erzielen, stehen aber +/- auch für die entfernten Variablen.
Nun können wir das globale Modell definieren, indem wir alle verbleibenden Variablen aufnehmen, das sind 16. (Wenn das nicht eh schon so viele wären, dass es uns an die Grenze der Rechenleistung bringt, hätten wir auch noch darüber nachdenken können, einzelne quadratische Terme oder Interaktionsterme zu berücksichtigen).
global_model <- lm(Species_richness ~
Inclination +
Heat_index +
Microrelief +
Conductivity +
Grazing_intensity +
Litter +
Precipitation +
Stones_and_rocks +
Gravel +
Fine_soil +
pH +
CaCO3 +
C_org +
CN_ratio +
Temperature +
Temperature_range,
data = stgl)Nun gibt es im Prinzip zwei Möglichkeiten, vom globalen (vollen) Modell zu einem minimal adäquaten Modell zu kommen. (1) Der Ansatz der „frequentist statistic“, in dem man aus dem vollen Modell so lange schrittweise Variablen entfernt, bis nur noch signifikante Variablen verbleiben. (2) Den informationstheoretischen Ansatz, bei dem alle denkbaren Modelle berechnet und verglichen werden (also alle möglichen Kombinationen von 13,12,…, 1, 0 Parametern). Diese Lösung stelle ich im Folgenden vor:
# Multimodel inference
library("MuMIn")
options(na.action = "na.fail")
allmodels <- dredge(global_model)allmodels[1:10]Global model call: lm(formula = Species_richness ~ Inclination + Heat_index + Microrelief +
Conductivity + Grazing_intensity + Litter + Precipitation +
Stones_and_rocks + Gravel + Fine_soil + pH + CaCO3 + C_org +
CN_ratio + Temperature + Temperature_range, data = stgl)
---
Model selection table
(Int) CCO CN_rat Cnd Fin_sol Grz_int Het_ind Ltt Prc
711 46.76 0.2287 -0.5878 1.1870 -13.23 -0.10000
21191 -48.97 0.2200 -0.7969 1.0450 -13.58 -0.09218 0.09409
21127 -57.17 0.2307 -0.8809 -13.31 -0.09696 0.10710
17095 28.08 0.2021 -0.6899 1.2050 -13.15 -0.10010
8903 47.66 0.2482 -0.6332 1.1200 -13.18 -0.10530
647 49.01 0.2373 -0.6716 -12.84 -0.10690
719 44.69 0.2005 -0.5903 0.01487 1.1330 -13.17 -0.10370
21207 -49.86 0.2207 -0.8155 -0.08360 1.0940 -13.24 -0.10170 0.09454
29383 -43.85 0.2361 -0.8256 0.9964 -13.53 -0.09695 0.08950
727 47.43 0.2308 -0.5971 -0.07154 1.2290 -12.93 -0.10820
Stn_and_rck Tmp df logLik AICc delta weight
711 7 -718.604 1451.8 0.00 0.152
21191 0.5016 9 -716.461 1451.9 0.08 0.146
21127 0.5309 8 -717.740 1452.2 0.44 0.122
17095 0.2387 8 -717.875 1452.5 0.71 0.106
8903 -0.07648 8 -717.999 1452.8 0.96 0.094
647 6 -720.253 1452.9 1.15 0.086
719 8 -718.188 1453.1 1.34 0.078
21207 0.5195 10 -716.031 1453.2 1.44 0.074
29383 -0.06288 0.4808 10 -716.048 1453.3 1.47 0.073
727 8 -718.293 1453.3 1.55 0.070
Models ranked by AICc(x) Jetzt bekommen wir die besten 10 der insgesamt 65536 möglichen Modelle gelistet mit ihren Parameterschätzungen und ihrem AICc.
Das beste Modell umfasst 5 Parameter (CaCO3, CN.ratio, Grazing.intensity. Heat.index, Litter). Allerdings ist das nächstbeste Modell (mit 6 Parametern) nur wenig schlechter (delta AICc = 0.71), was sich in fast gleichen (und zudem sehr niedrigen) Akaike weights bemerkbar macht. Nach dem Verständnis des Information theoretician approach, sollte man in einer solchen Situation nicht das eine „beste“ Modell benennen, sondern eine Aussage über die Gruppe der insgesamt brauchbaren Modelle treffen. Hierzu kann man (a) Importance der Parameter über alle Modelle hinweg berechnen (= Summe der Akaike weights aller Modelle, die den betreffenden Parameter enthalten) und/oder (b) ein nach Akaike weights gemitteltes Modell berechnen.
# Importance values der Variablen
sw(allmodels) Heat_index Litter CaCO3 CN_ratio Grazing_intensity
Sum of weights: 1.00 0.91 0.81 0.77 0.63
N containing models: 32768 32768 32768 32768 32768
Temperature Precipitation Stones_and_rocks Conductivity
Sum of weights: 0.49 0.43 0.40 0.33
N containing models: 32768 32768 32768 32768
Microrelief Fine_soil Gravel pH C_org Temperature_range
Sum of weights: 0.32 0.30 0.30 0.28 0.28 0.28
N containing models: 32768 32768 32768 32768 32768 32768
Inclination
Sum of weights: 0.26
N containing models: 32768 Demnach ist Heat.index die wichtigste Variable (in 100% aller relevanten Modelle), während ferner Litter, CaCO3, CN_ratio und Grazing_intensity in mehr als 50% der relevanten Modelle enthalten sind.
# Modelaveraging (Achtung: dauert mit 13 Variablen einige Minuten)
avgmodel <- model.avg(allmodels)
summary(avgmodel)$coefficients (Intercept) CaCO3 CN_ratio Grazing_intensity Heat_index
full 16.35327 0.1814323 -0.5348808 0.7724891 -12.96429
subset 16.35327 0.2242980 -0.6981476 1.2198244 -12.96432
Litter Precipitation Temperature Stones_and_rocks Conductivity
full -0.09900428 0.03134582 0.2028271 -0.03786674 0.005393491
subset -0.10877370 0.07355793 0.4145035 -0.09379116 0.016472114
Fine_soil Gravel pH Temperature_range Microrelief
full -0.01950678 -0.02149269 -0.1777904 -0.01553740 0.01740667
subset -0.06523683 -0.07235735 -0.6363338 -0.05583764 0.05356550
C_org Inclination
full 0.0584467 0.003167463
subset 0.2093569 0.012378658Aus dem gemittelten Modell können wir die Richtung der Beziehung (positiv oder negativ) und ggf. die Effektgrössen (wie verändert sich die Artenzahl, wenn die Prädiktorvariable um eine Einheit zunimmt?) ermitteln.
# Modelldiagnostik nicht vergessen
par(mfrow = c(2, 2))
plot(global_model)
plot( lm(Species_richness ~
Heat_index + Litter + CaCO3 + CN_ratio + Grazing_intensity,
data = stgl))
Wie immer kommt am Ende die Modelldiagnostik. Wir können uns entweder das globale Modell oder das Modell mit den 5 Variablen mit importance > 50% anschauen. Das Bild sieht fast identisch aus und zeigt keinerlei problematische Abweichungen, d. h. links oben weder ein Keil, noch eine Banane, rechts oben eine nahezu perfekte Gerade.