1 Das GLM Lasso penalty Modell:

Der Lasso-Penalty wird verwendet, um die Anzahl der verwendeten Prädiktoren (Einflussgrößen) in einem Modell zu reduzieren und gleichzeitig die Modellkomplexität zu kontrollieren. Dabei wird ein Strafterm hinzugefügt, der die Summe der absoluten Werte der Koeffizienten mit einem Tuning-Parameter (Lambda) multipliziert. Durch die Optimierung des Strafterms werden einige der Koeffizienten auf Null reduziert, wodurch das Modell automatisch eine Variablenselektion durchführt und irrelevante Prädiktoren ausschließt. Dies hat den Vorteil, dass Modelle erzeugt werden, die nur die wichtigsten Prädiktoren nutzen. Das erleichtert die Interpretation der Ergebnisse, insbesondere wenn die Anzahl der Prädiktoren groß ist im Vergleich zur Anzahl der Beobachtungen.

Im Zusammenhang mit dem GLM mit Lasso Penalty existieren drei verschiedene Arten von Regressionsmodellen, die je nach Art der abhängigen Variablen verwendet werden können:

1.Lineare Regression: family = “gaussian”

  min(β0,β)∈Rp+112N∑i=1N(yi−β0−xTiβ)2+λ[(1−α)∥β∥22/2+α∥β∥1],

Dies ist das Standard-Argument für die Funktion glmnet. Dieses wird verwendet, wenn die abhängige Variable kontinuierlich ist und einer “Gaußschen” (normalen) Verteilung folgt. Das Ziel der linearen Regression ist es, einen linearen Zusammenhang zwischen den Variablen herzustellen, der die beobachteten Daten möglichst gut widerspiegelt und die Abweichung zwischen den beobachteten und vorhergesagten Werten minimiert. Die Summe der quadrierten Residuen ist ein Maß für die Güte der Anpassung des Modells an die Daten.

Residuen sind die Differenzen zwischen den beobachteten Werten einer abhängigen Variablen und den vorhergesagten Werten, die durch das lineare Regressionsmodell berechnet werden. Sie stellen die nicht erklärte Varianz der abhängigen Variablen dar. Der Strafterm in der Ziel-Funktion der Lasso-Regression und Ridge-Regression ist eine Kombination aus L2 (Ridge) und L1 (Lasso) Straftermen.

Die L2-Strafterm (Ridge) verwendet die quadrierten Koeffizienten, während der L1-Strafterm (Lasso) die Beträge der Koeffizienten verwendet. Durch die Kombination beider Strafterme wird eine Regularisierung erreicht, die dazu beiträgt, die Koeffizienten zu schrumpfen und den Einfluss einzelner Prädiktoren zu reduzieren.

Der L2-Strafterm (Ridge) hilft, Overfitting zu reduzieren, indem er die Koeffizienten schrumpft und die Modellkomplexität kontrolliert

Der L1-Strafterm (Lasso) hat eine sparsifizierende Wirkung, da er einige der Koeffizienten auf Null reduziert, was zu einer Variablenselektion führt.

2.Lineare Regression: family = “mgaussian” (Multi-Response)

  min(β0,β)∈R(p+1)×K12N∑i=1N∥yi−β0−βTxi∥2F+λ[(1−α)∥β∥2F/2+α∑j=1p∥βj∥2].

Die Multi-Response “Gaußsche” family wird verwendet, wenn es mehrere korrelierte abhängige Variablen gibt. Dies ist nützlich im Kontext von Multi-Task-Learning-Problemen. Die abhängige Variable wird als Matrix von quantitativen Variablen dargestellt. Das Ziel der linearen Regression ist ähnlich wie bei der Gaußschen Familie, es wird jedoch auf die Matrix der abhängigen Variablen und Koeffizienten angewendet.

Es besteht darin, einen linearen Zusammenhang zwischen den unabhängigen Variablen und der Matrix der abhängigen Variablen herzustellen, der die beobachteten Daten möglichst gut widerspiegelt und die Abweichung zwischen den beobachteten und vorhergesagten Werten minimiert.

Der Unterschied zur “linear regression”normalen linear regression” liegt darin, dass die Koeffizientenmatrix in der Multi-Response Regression verwendet wird, um den Zusammenhang zwischen den unabhängigen Variablen, der Matrix und der abhängigen Variablen zu modellieren.

3.Logistische Regression: family = “binomial”

   min(β0,β)∈Rp+1−[1N∑i=1Nyi⋅(β0+xTiβ)−log(1+e(β0+xTiβ))]+λ[(1−α)∥β∥22/2+α∥β∥1].

Die logistische Regression wird häufig verwendet, wenn die abhängige Variable binär ist, d.h. sie nimmt nur zwei Werte an, wie zum Beispiel Erfolg/Misserfolg oder Ja/Nein. Die Verteilung der abhängigen Variable wird als binomial angenommen, da sie entweder zu einer der beiden Kategorien gehört. Ähnlich wie bei der linearen Regression werden auch in der logistischen Regression Strafterme verwendet, um die Modellkomplexität zu kontrollieren und Überanpassungen zu vermeiden. Der Strafterm kombiniert L2 (Ridge) und L1 (Lasso) Strafterme, um eine angemessene Schrumpfung der Koeffizienten zu erreichen.

Diese verschiedenen Regressionsmodelle ermöglichen verschiedene Arten von abhängigen Variablen und Verteilungsannahmen und bieten somit Flexibilität bei der Modellierung verschiedener Datentypen.

2 Laden aller benötigten R Packete

# Laden aller erforderlichen Pakete
library("terra")
library("geodata")
library("raster")
library("glmnet")
library("readr")
library("caret")
library("spdep")
library("sf")
library("lwgeom")

3 Einlesen der Datensätze und erstellen einer geographischen Karte von Pakistan

# Datensatz einlesen
d <- read_csv("PakistanLadakh.csv")

# Karte von Pakistan erstellen
pak <- geodata::gadm(country = 'PAK', level = 0, path = getwd())
plot(pak)

4 Geographische Karte von Pakistan mit Ländergrenzen und Geographischen Höhen

# Höhendaten einladen und auf Pakistan-Grenzen zuschneiden
elevation <- geodata::worldclim_global(var = "elev", res = 10, path = getwd())
cropped_elevation <- crop(elevation, ext(pak))

# Maske basierend auf Pakistan-Grenzen erstellen
clipped_elevation <- mask(cropped_elevation, pak)

# Karte mit Höhendaten plotten
plot(clipped_elevation$wc2.1_10m_elev)

# Karte von Pakistan mit Grenzen plotten
plot(pak, add = TRUE)

5 Visualiesierung der Pakistan Karte mit Verbreitung aller Arten

Dies dient als Überblick für die Verbreitung aller Schmertterlingsarten

#Karte plotten 
plot(clipped_elevation$wc2.1_10m_elev)

# Karte von Pakistan mit Grenzen plotten
plot(pak, add = TRUE)
# Schmetterlingsarten innerhalb von Pakistan filtern
butterflies_pakistan <- subset(d, x >= xmin(pak) & x <= xmax(pak) & y >= ymin(pak) & y <= ymax(pak))

# Schmetterlingsdaten auf Karte plotten
points(butterflies_pakistan$x, butterflies_pakistan$y, col = 'black')

6 Visualiesierung der Karte von Pakistan mit zufällig generierten Background Points

#Karte  plotten 
plot(clipped_elevation$wc2.1_10m_elev)

# Karte von Pakistan mit Grenzen plotten
plot(pak, add = TRUE)
# Zufällige Hintergrundpunkte generieren und auf Karte plotten
background_points <- sf::st_sample(st_as_sf(pak), size = 10000)
background_points <- st_coordinates(background_points)
background_points <- as.data.frame(background_points)
points(background_points$X, background_points$Y, col ="red")

7 Datenvorbereitung:

Filtern der Höhendaten, Beobachtungen und Aufteilen der Datensätze in Test und Trainingsdaten für die Modellierung

# Datenvorbereitung
elevation_train <- extract(clipped_elevation$wc2.1_10m_elev, d[,c("x", "y")])
d$elevation <- elevation_train$wc2.1_10m_elev
d <- na.omit(d)

# Filtern der Klassen mit ausreichender Anzahl an Beobachtungen
class_counts <- table(d$species)
filtered_classes <- names(class_counts[class_counts > 14])
filtered_data <- subset(d, species %in% filtered_classes)
X_filtered <- as.data.frame(filtered_data[, c("x", "y", "elevation")])
Y_filtered <- factor(filtered_data$species)

# Daten aufteilen in Trainings- und Testdaten
set.seed(123)
train_indices <- createDataPartition(Y_filtered, p = 0.7, list = FALSE)
X_train_filtered <- X_filtered[train_indices, ]
Y_train_filtered <- Y_filtered[train_indices]
X_test_filtered <- X_filtered[-train_indices, ]
Y_test_filtered <- Y_filtered[-train_indices]

8 GLM Lasso Modell generieren und darstellen

Als “2norm” werden alle Daten als ein Plot dargestellt, aber auch möglich als “coef”. Bei “coef” wird für jede Species eine seperate Grafik generiert

# GLM Lasso-Modell```{r}
glm_lasso <- glmnet(X_train_filtered, factor(Y_train_filtered), family = "multinomial")

# Plot des Modells
plot(glm_lasso, xvar = "lambda", type.coef = "2norm")

9 Ermitteln eines optimalen Lambda wertes

# Ermitteln des optimalen Lambda-Wertes
cv_fit <- cv.glmnet(as.matrix(X_train_filtered), factor(Y_train_filtered), family = "multinomial")
lambda_min <- cv_fit$lambda.min
plot(cv_fit)

# Vorhersagen für Testdaten mit lambda.min
Y_hat <- predict(cv_fit, newx = as.matrix(X_test_filtered), s = "lambda.min", type = "class")

10 Evaluierung der Modelgenauigkeit

Der errechnete Wert muss * 100 genommen werden um die % für die Genauigkeit zu bekommen

# Evaluierung der Modellgenauigkeit
accuracy <- sum(Y_hat == Y_test_filtered) / length(Y_test_filtered)
print(accuracy)    
## [1] 0.02744425

11 Generierte Vorhersagen

Der Befehl “print(background_points_with_predicted_species)” erzeugt eine Tabelle alle Arten mit den jeweiligen Koordinaten und die dazugehörige Höhe für die jeweils wahrscheinlichste Schmetterlingsart

# Vorhersagen für alle Hintergrundpunkte
elevation_background <- extract(clipped_elevation$wc2.1_10m_elev, background_points[,c("X", "Y")])
background_points$elevation <- elevation_background$wc2.1_10m_elev
predicted_butterflies <- predict(cv_fit, newx = as.matrix(background_points), s = "lambda.min", type = "class")
background_points_with_predicted_species <- background_points
background_points_with_predicted_species$species <- predicted_butterflies

12 Resultate als csv. Datei

Ergebnisse für alle Schmetterlingsarten werden als csv. Datei abgespeichert

write.csv(background_points_with_predicted_species, "result.csv")

13 Beispiel

Ein Beispiel anhand zwei konkreten Schmetterlingsarten( Colotis_vestalis und Colias_erate) mit Vorhersage -> Die Arten können nach Wunsch geändert werden

# Vorhersage für einzelne Schmetterlingsarten plotten
plot(clipped_elevation$wc2.1_10m_elev)
plot(pak, add = TRUE)

selection <- which(background_points_with_predicted_species$species == "Colotis_vestalis")
test <- background_points_with_predicted_species[selection,]
points(test$X, test$Y)
text(test$X[1], test$Y[1], "Colotis_vestalis", col = "red")

selection <- which(background_points_with_predicted_species$species == "Colias_erate")
test <- background_points_with_predicted_species[selection,]
points(test$X, test$Y, col = 6)
text(test$X[1], test$Y[1], "Colias_erate", col = "red")

14 Vorhersage auf der Pakistan Karte mit allen relevanten Arten

# Vorhersage für alle Schmetterlingsarten plotten
i <- 1
plot(clipped_elevation$wc2.1_10m_elev)
plot(pak, add = TRUE)
legend_colors <- vector("integer", length = length(unique(predicted_butterflies)))

for (specie in unique(predicted_butterflies)) {
  selection <- which(background_points_with_predicted_species$species == specie)
  test <- background_points_with_predicted_species[selection, ]
  points(test$X, test$Y, col = i)
  legend_colors[i] <- i
  i <- i + 1
}

15 Erstellen einer Legende für die Schmetterlingsarten

plot(NULL, type = “n”) num_cols <- 3 legend(“center”, legend = as.vector(unique(predicted_butterflies)), fill = legend_colors, title = “Schmetterlingsarten”, cex = 0.5, pt.cex = 1, ncol = num_cols)

16 Literatur