library("sf")
library("terra")
library("dplyr")
library("readr")
library("ggplot2")
library("PerformanceAnalytics")
library("pastecs")
library("lme4")
library("bbmle")
library("MuMIn")
library("MASS")BiEc5_N Variablenselektion MM
Libraries laden
Variablenselektion
→ Vorgehen analog Coppes u. a. (2017)
Aufgabe 1
Mit dem folgenden Code kann eine simple Korrelationsmatrix aufgebaut werden, vergl. Aufgabe 5 vorangehende Woche
DF_mod <- read_delim("datasets/fallstudie_n/Aufgabe4_Datensatz_Habitatnutzung_Modelle_231027_moodle.csv", delim = ";")
DF_mod_day <- DF_mod |>
filter(time_of_day == "day")
round(cor(DF_mod_day[, 6:14], method = "kendall"), 2)
## slope topo_pos dist_road_all dist_road_only dist_sett remoteness
## slope 1.00 0.07 0.14 0.17 0.10 0.59
## topo_pos 0.07 1.00 0.05 0.08 0.04 0.00
## dist_road_all 0.14 0.05 1.00 0.85 -0.05 0.35
## dist_road_only 0.17 0.08 0.85 1.00 -0.06 0.34
## dist_sett 0.10 0.04 -0.05 -0.06 1.00 0.06
## remoteness 0.59 0.00 0.35 0.34 0.06 1.00
## forest_prop 0.17 0.00 -0.09 -0.09 0.40 0.12
## us_2014 0.22 -0.04 -0.06 -0.03 0.09 0.13
## os_2014 0.34 0.04 -0.06 -0.04 0.23 0.24
## forest_prop us_2014 os_2014
## slope 0.17 0.22 0.34
## topo_pos 0.00 -0.04 0.04
## dist_road_all -0.09 -0.06 -0.06
## dist_road_only -0.09 -0.03 -0.04
## dist_sett 0.40 0.09 0.23
## remoteness 0.12 0.13 0.24
## forest_prop 1.00 0.31 0.52
## us_2014 0.31 1.00 0.42
## os_2014 0.52 0.42 1.00
# hier kann die Schwelle für die Korrelation gesetzt werden, 0.7 ist liberal /
# 0.5 konservativ
cor <- round(cor(DF_mod_day[, 6:14], method = "kendall"), 2)
cor[abs(cor) < 0.7] <- 0
cor
## slope topo_pos dist_road_all dist_road_only dist_sett remoteness
## slope 1 0 0.00 0.00 0 0
## topo_pos 0 1 0.00 0.00 0 0
## dist_road_all 0 0 1.00 0.85 0 0
## dist_road_only 0 0 0.85 1.00 0 0
## dist_sett 0 0 0.00 0.00 1 0
## remoteness 0 0 0.00 0.00 0 1
## forest_prop 0 0 0.00 0.00 0 0
## us_2014 0 0 0.00 0.00 0 0
## os_2014 0 0 0.00 0.00 0 0
## forest_prop us_2014 os_2014
## slope 0 0 0
## topo_pos 0 0 0
## dist_road_all 0 0 0
## dist_road_only 0 0 0
## dist_sett 0 0 0
## remoteness 0 0 0
## forest_prop 1 0 0
## us_2014 0 1 0
## os_2014 0 0 1Aufgabe 2
Skalieren der Variablen, damit ihr Einfluss vergleichbar wird (Befehl scale(); Problem verschiedene Skalen der Variablen (bspw. Neigung in Grad, Distanz in Metern)); Umwandeln der Reh-ID in einen Faktor, damit dieser als Random Factor ins Model eingespiesen werden kann.
DF_mod_day <- DF_mod_day |>
mutate(
slope_scaled = scale(slope),
topo_pos_scaled = scale(topo_pos),
us_scaled = scale(us_2014),
os_scaled = scale(os_2014),
forest_prop_scaled = scale(forest_prop),
dist_road_all_scaled = scale(dist_road_all),
dist_road_only_scaled = scale(dist_road_only),
dist_sett_scaled = scale(dist_sett),
remoteness_scaled = scale(remoteness),
id = as.factor(id)
)Aufgabe 3
Selektion der Variablen in einem univariaten Model
Ein erstes GLMM (Generalized Linear Mixed Effects Modell) aufbauen: Funktion und Modelformel
wichtige Seite auf der man viele Hilfestellungen zu GLMM’s finden kann.
# wir werden das package lme4 mit der Funktion glmer verwenden
# die Hilfe von glmer aufrufen: ?glmer
# glmer(formula, data = , family = binomial)
# 1) formula:
# Abhängige Variable ~ Erklärende Variable + Random Factor
# In unseren Modellen kontrollieren wir für individuelle Unterschiede bei den Rehen
# indem wir einen Random Factor definieren => (1 | id)
# 2) data:
# euer Datensatz
# 3) family:
# hier binomial
# warum binomial? Verteilung Daten der Abhängigen Variable Präsenz/Absenz
ggplot(DF_mod_day, aes(pres_abs)) +
geom_histogram()
# --> Binäre Verteilung => Binomiale Verteilung mit n = 1
# und wie schaut die Verteilung der Daten der Abhängigen Variable Nutzungsintensität
# (nmb, werden wir in diesem Kurs aber nicht genauer anschauen) aus?Aufgabe 4
Mit der GLMM Formel bauen wir in einem ersten Schritt eine univariate Variablenselektion auf.
Als abhängige Variable verwenden wir die Präsenz/Absenz der Rehe in den Kreisen
# Die erklärende Variable in m1 ist die erste Variable der korrelierenden Variablen
# Die erklärende Variable in m2 ist die zweite Variable der korrelierenden Variablen
m1 <- glmer(Abhaengige_Variable ~ Erklaerende_Variable + (1 | id),
data = DF_mod_day,
family = binomial
)
m2 <- glmer(Abhaengige_Variable ~ Erklaerende_Variable + (1 | id),
data = DF_mod_day,
family = binomial
)
# mit dieser Funktion können die Modellergebnisse inspiziert werden
summary(m1)
# Mit dieser Funktion kann der Informationgehalt der beiden Modelle gegeneinander
# abgeschätzt werden
bbmle::AICtab(m1, m2)
# tieferer AIC -> besser (AIC = Akaike information criterion)
# ==> dieses Vorgehen muss nun für alle korrelierten Variablen für jeden Teildatensatz
# (Tag/Nacht) durchgeführt werden, um nur noch nicht (R < 0.7) korrelierte Variablen
# in das Modell einfliessen zu lassenCode
m1 <- glmer(pres_abs ~ dist_road_all_scaled + (1 | id), data = DF_mod_day, family = binomial)
m2 <- glmer(pres_abs ~ dist_road_only_scaled + (1 | id), data = DF_mod_day, family = binomial)
summary(m1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: pres_abs ~ dist_road_all_scaled + (1 | id)
## Data: DF_mod_day
##
## AIC BIC logLik deviance df.resid
## 5104.9 5123.7 -2549.4 5098.9 3961
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1046 -0.7564 -0.6277 1.0812 1.9219
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 0.1996 0.4467
## Number of obs: 3964, groups: id, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.31120 0.13418 -2.319 0.0204 *
## dist_road_all_scaled 0.37745 0.03848 9.808 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## dst_rd_ll_s -0.004
bbmle::AICtab(m1, m2)
## dAIC df
## m2 0.0 3
## m1 6.5 3
# tieferer AIC -> besser (AIC = Akaike information criterion) -> als deltaAIC
# ausgewiesen besser == Distanz zu Strassen
# ==> dieses Vorgehen muss nun für alle korrelierten Variablen für jeden
# Teildatensatz (geringe Störung/starke Störung) durchgeführt werden, um nur
# noch nicht (R < 0.7) korrelierte Variablen in das Modell einfliessen zu
# lassenAufgabe 5
Selektion der Variablen in einem multivariaten Model
Mit folgendem Code kann eine automatisierte Variablenselektion (dredge-Funktion) und ein Modelaveraging aufgebaut werden (siehe auch Stats-Skript von J.Dengler & Team)
# hier wird die Formel für die dredge-Funktion vorbereitet (die Variablen V1-V8
# sind jene welche nach der univariaten Variablenselektion noch übrig bleiben)
f <- pres_abs ~
V1 +
V2 +
V3 +
V4 +
V5 +
V6 +
V7 +
V8
# in diesem Befehl kommt der Random-Factor (das Reh) hinzu und es wird eine Formel
# daraus gemacht
f_dredge <- paste(c(f, "+ (1 | id)"), collapse = " ") |> as.formula()
# Das Modell mit dieser Formel ausführen
m <- glmer(f_dredge, data = DF_mod_day, family = binomial, na.action = "na.fail")
# Das Modell in die dredge-Funktion einfügen (siehe auch unbedingt ?dredge)
all_m <- dredge(m)
# Importance values der einzelnen Variablen (Gibt an, wie bedeutsam eine bestimmte
# Variable ist, wenn man viele verschiedene Modelle vergleicht (multimodel inference))
sw(all_m)
# Schlussendlich wird ein Modelaverage durchgeführt (Schwellenwert für das delta-AIC = 2)
avgmodel <- model.avg(all_m, rank = "AICc", subset = delta < 2)
summary(avgmodel)
# ==> für den Nachtdatensatz muss der gleiche Prozess der Variablenselektion
# durchgespielt werden.Code
# hier wird die Formel für die dredge-Funktion vorbereitet (die Variablen V1-V8
# sind jene welche nach der univariaten Variablenselektion noch übrig bleiben)
f <- pres_abs ~
slope_scaled +
topo_pos_scaled +
us_scaled +
os_scaled +
forest_prop_scaled +
dist_road_only_scaled +
dist_sett_scaled +
remoteness_scaled
# inn diesem Befehl kommt der Random-Factor (das Reh) hinzu und es wird eine Formel
# daraus gemacht
f_dredge <- paste(c(f, "+ (1 | id)"), collapse = " ") |> as.formula()
# Das Modell mit dieser Formel ausführen
m <- glmer(f_dredge, data = DF_mod_day, family = binomial, na.action = "na.fail")
# Das Modell in die dredge-Funktion einfügen (siehe auch unbedingt ?dredge)
all_m <- dredge(m)
# Importance values der einzelnen Variablen (Gibt an, wie bedeutsam eine bestimmte
# Variable ist, wenn man viele verschiedene Modelle vergleicht (multimodel inference))
sw(all_m)
## dist_road_only_scaled forest_prop_scaled us_scaled
## Sum of weights: 1.00 1.00 1.00
## N containing models: 128 128 128
## remoteness_scaled dist_sett_scaled slope_scaled os_scaled
## Sum of weights: 1.00 0.42 0.32 0.29
## N containing models: 128 128 128 128
## topo_pos_scaled
## Sum of weights: 0.27
## N containing models: 128
# Schlussendlich wird ein Modelaverage durchgeführt (Schwellenwert für das delta-AIC = 2)
avgmodel <- model.avg(all_m, rank = "AICc", subset = delta < 2)
summary(avgmodel)
##
## Call:
## model.avg(object = get.models(object = all_m, subset = delta <
## 2), rank = "AICc")
##
## Component model call:
## glmer(formula = pres_abs ~ <4 unique rhs>, data = DF_mod_day, family =
## binomial, na.action = na.fail)
##
## Component models:
## df logLik AICc delta weight
## 1357 6 -2261.70 4535.41 0.00 0.39
## 12357 7 -2261.02 4536.07 0.65 0.28
## 13567 7 -2261.52 4537.06 1.65 0.17
## 13457 7 -2261.60 4537.23 1.82 0.16
##
## Term codes:
## dist_road_only_scaled dist_sett_scaled forest_prop_scaled
## 1 2 3
## os_scaled remoteness_scaled slope_scaled
## 4 5 6
## us_scaled
## 7
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) -0.437730 0.148042 0.148087 2.956 0.00312 **
## dist_road_only_scaled 0.510519 0.052598 0.052614 9.703 < 2e-16 ***
## forest_prop_scaled 0.844263 0.058809 0.058826 14.352 < 2e-16 ***
## remoteness_scaled -0.255199 0.056758 0.056775 4.495 7e-06 ***
## us_scaled 0.382072 0.040946 0.040959 9.328 < 2e-16 ***
## dist_sett_scaled -0.020550 0.046737 0.046744 0.440 0.66021
## slope_scaled -0.006037 0.027844 0.027851 0.217 0.82839
## os_scaled 0.004232 0.026310 0.026317 0.161 0.87225
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) -0.43773 0.14804 0.14809 2.956 0.00312 **
## dist_road_only_scaled 0.51052 0.05260 0.05261 9.703 < 2e-16 ***
## forest_prop_scaled 0.84426 0.05881 0.05883 14.352 < 2e-16 ***
## remoteness_scaled -0.25520 0.05676 0.05678 4.495 7e-06 ***
## us_scaled 0.38207 0.04095 0.04096 9.328 < 2e-16 ***
## dist_sett_scaled -0.07295 0.06270 0.06272 1.163 0.24479
## slope_scaled -0.03532 0.05918 0.05919 0.597 0.55069
## os_scaled 0.02690 0.06157 0.06159 0.437 0.66226
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
#| eval: false
#| error: true
#| echo: false
# hier wird die Formel für die dredge-Funktion vorbereitet (die Variablen V1-V8
# sind jene welche nach der univariaten Variablenselektion noch übrig bleiben)
f <- pres_abs ~
slope_scaled +
topo_pos_scaled +
us_scaled +
os_scaled +
forest_prop_scaled +
dist_road_only_scaled +
dist_road_all_scaled +
dist_sett_scaled +
remoteness_scaled
# inn diesem Befehl kommt der Random-Factor (das Reh) hinzu und es wird eine Formel
# daraus gemacht
f_dredge <- paste(c(f, "+ (1 | id)"), collapse = " ") |> as.formula()
# Das Modell mit dieser Formel ausführen
m <- glmer(f_dredge, data = DF_mod_day, family = binomial, na.action = "na.fail")
car::vif(m)
## slope_scaled topo_pos_scaled us_scaled
## 2.228320 1.098040 1.239462
## os_scaled forest_prop_scaled dist_road_only_scaled
## 2.613227 2.737806 5.816134
## dist_road_all_scaled dist_sett_scaled remoteness_scaled
## 6.032018 1.434214 2.554936
mean(car::vif(m))
## [1] 2.861573