Logistic Regression — Species Distribution Modelling

Penulis

Afiliasi

Dr. Arif K Wijayanto

Laboratory of Environmental Analysis and Geospatial Modeling, Dept of Forest Resources Conservation and Ecotourism, IPB University

Diterbitkan

8 Juni 2026

Petunjuk: Jalankan setiap chunk secara individual dengan tombol ▶ Run, atau jalankan semua sekaligus melalui Run › Run All. Output muncul di bawah setiap chunk.

1 Install dan Load Package

options(repos = c(CRAN = "https://cran.rstudio.com/"))

pkgs     <- c("car", "ResourceSelection", "DescTools", "corrplot", "pROC")
new_pkgs <- pkgs[!pkgs %in% installed.packages()[, "Package"]]
if (length(new_pkgs)) install.packages(new_pkgs)

library(car)
library(ResourceSelection)
library(DescTools)
library(ggplot2)
library(corrplot)
library(pROC)

2 Baca Data

# Ganti dengan lokasi file Anda
data <- read.csv("sample.csv")

head(data)

summary(data)

3 Buat Model Awal (Semua Prediktor)

Membangun model regresi logistik penuh dengan semua variabel prediktor.

model_all <- glm(
  PA ~ `euclid.jalan` + `euclid.desa` + `euclid.sungai` + `euclid.hutan` +
       `euclid.situs` + Slope + Elevasi + NDVI,
  data   = data,
  family = binomial
)

summary(model_all)

4 Cek Multikolinearitas — VIF

Nilai VIF > 10 (atau > 5 secara konservatif) mengindikasikan multikolinearitas yang bermasalah.

cat("=== VIF ===\n")
print(vif(model_all))

5 Korelasi Antar Prediktor

Matriks korelasi Pearson untuk melihat hubungan linear antar prediktor.

vars <- data[, c("euclid.jalan", "euclid.desa", "euclid.sungai", "euclid.hutan",
                 "Slope", "Elevasi", "NDVI")]

correlation_matrix <- cor(vars, method = "pearson")
print(correlation_matrix)

corrplot(
  correlation_matrix,
  method     = "number",
  type       = "upper",
  tl.col     = "black",
  tl.srt     = 45,
  tl.cex     = 1.2,
  cl.pos     = "n",
  number.cex = 1.2
)

6 Perbandingan Model Alternatif

Jika Slope dan Elevasi berkorelasi tinggi, kita bandingkan dua model yang masing-masing hanya mengandung satu dari keduanya. Model dengan AIC lebih rendah adalah yang lebih baik.

model_slope <- glm(
  PA ~ `euclid.jalan` + `euclid.desa` + `euclid.sungai` + `euclid.hutan` +
       Slope + NDVI,
  data   = data,
  family = binomial
)

summary(model_slope)

model_elevasi <- glm(
  PA ~ `euclid.jalan` + `euclid.desa` + `euclid.sungai` + `euclid.hutan` +
       Elevasi + NDVI,
  data   = data,
  family = binomial
)

summary(model_elevasi)

cat("=== Perbandingan AIC ===\n")
print(AIC(model_slope, model_elevasi))

7 Evaluasi Model Terbaik

Menggunakan model_elevasi sebagai model terbaik. Ganti dengan model_slope jika AIC-nya lebih rendah.

best_model <- model_elevasi  # Ganti sesuai model yang dipilih

7.1 Kurva ROC dan AUC

Kurva ROC mengukur kemampuan diskriminasi model. Nilai AUC > 0.7 menunjukkan performa yang dapat diterima; AUC > 0.9 sangat baik.

roc_obj <- roc(data$PA, fitted(best_model))
plot(roc_obj, col = "steelblue", main = paste("AUC =", round(auc(roc_obj), 3)))

7.2 Hosmer–Lemeshow Goodness-of-Fit

p-value > 0.05 menunjukkan model fit yang baik (tidak ada perbedaan signifikan antara nilai observasi dan prediksi).

cat("=== Hosmer-Lemeshow ===\n")
print(hoslem.test(data$PA, fitted(best_model), g = 10))

7.3 Pseudo-R² (McFadden & Nagelkerke)

Mengukur kekuatan penjelas model. Nagelkerke R² > 0.2 umumnya dianggap dapat diterima.

cat("=== Nagelkerke R² ===\n")
print(PseudoR2(best_model, c("McFadden", "Nagel")))

7.4 Signifikansi Koefisien

Koefisien estimasi, standard error, z-value, dan p-value untuk setiap prediktor.

cat("=== Signifikansi tiap variabel ===\n")
print(summary(best_model)$coefficients)

8 Informasi Sesi

sessionInfo()

R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: Asia/Jakarta
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] pROC_1.19.0.1           corrplot_0.95           ggplot2_4.0.1          
[4] DescTools_0.99.60       ResourceSelection_0.3-6 car_3.1-3              
[7] carData_3.0-5          

loaded via a namespace (and not attached):
 [1] generics_0.1.4     class_7.3-23       lattice_0.22-7     hms_1.1.4         
 [5] digest_0.6.39      magrittr_2.0.4     evaluate_1.0.5     grid_4.5.2        
 [9] RColorBrewer_1.1-3 mvtnorm_1.3-7      fastmap_1.2.0      cellranger_1.1.0  
[13] jsonlite_2.0.0     Matrix_1.7-4       e1071_1.7-16       Formula_1.2-5     
[17] httr_1.4.7         scales_1.4.0       abind_1.4-8        cli_3.6.5         
[21] rlang_1.1.6        expm_1.0-0         withr_3.0.2        yaml_2.3.11       
[25] rootSolve_1.8.2.4  tools_4.5.2        tzdb_0.5.0         lmom_3.3          
[29] gld_2.6.8          dplyr_1.1.4        Exact_3.3          forcats_1.0.1     
[33] boot_1.3-32        vctrs_0.6.5        R6_2.6.1           proxy_0.4-27      
[37] lifecycle_1.0.4    fs_1.6.6           htmlwidgets_1.6.4  MASS_7.3-65       
[41] pkgconfig_2.0.3    pillar_1.11.1      gtable_0.3.6       data.table_1.17.8 
[45] glue_1.8.0         Rcpp_1.1.0         tidyselect_1.2.1   haven_2.5.5       
[49] xfun_0.54          tibble_3.3.0       rstudioapi_0.18.0  knitr_1.50        
[53] farver_2.1.2       htmltools_0.5.9    rmarkdown_2.30     readr_2.1.6       
[57] compiler_4.5.2     S7_0.2.1           readxl_1.4.5

Dokumen ini dibuat dengan Quarto dan dapat direproduksi dengan menjalankan ulang seluruh chunk menggunakan data yang sesuai.

--- title: "Logistic Regression — Species Distribution Modelling" author: - name: "Dr. Arif K Wijayanto" affiliation: "Laboratory of Environmental Analysis and Geospatial Modeling, Dept of Forest Resources Conservation and Ecotourism, IPB University" date: today lang: id format: html: theme: cosmo toc: true toc-depth: 3 toc-title: "Daftar Isi" toc-location: left number-sections: true number-depth: 3 code-fold: false code-tools: true code-copy: true code-line-numbers: true highlight-style: github df-print: kable fig-align: center fig-cap-location: bottom tbl-cap-location: top embed-resources: true smooth-scroll: true page-layout: full execute: echo: true eval: false warning: false message: false --- > **Petunjuk:** Jalankan setiap chunk secara individual dengan tombol **▶ Run**, atau jalankan semua sekaligus melalui **Run › Run All**. Output muncul di bawah setiap chunk. --- # Install dan Load Package ```{r} #| label: packages #| eval: true options(repos = c(CRAN = "https://cran.rstudio.com/")) pkgs <- c("car", "ResourceSelection", "DescTools", "corrplot", "pROC") new_pkgs <- pkgs[!pkgs %in% installed.packages()[, "Package"]] if (length(new_pkgs)) install.packages(new_pkgs) library(car) library(ResourceSelection) library(DescTools) library(ggplot2) library(corrplot) library(pROC) ``` --- # Baca Data ```{r} #| label: load-data # Ganti dengan lokasi file Anda data <- read.csv("sample.csv") head(data) ``` ```{r} #| label: data-summary summary(data) ``` --- # Buat Model Awal (Semua Prediktor) Membangun model regresi logistik penuh dengan semua variabel prediktor. ```{r} #| label: model-all model_all <- glm( PA ~ `euclid.jalan` + `euclid.desa` + `euclid.sungai` + `euclid.hutan` + `euclid.situs` + Slope + Elevasi + NDVI, data = data, family = binomial ) summary(model_all) ``` --- # Cek Multikolinearitas — VIF Nilai VIF **> 10** (atau > 5 secara konservatif) mengindikasikan multikolinearitas yang bermasalah. ```{r} #| label: vif cat("=== VIF ===\n") print(vif(model_all)) ``` --- # Korelasi Antar Prediktor Matriks korelasi Pearson untuk melihat hubungan linear antar prediktor. ```{r} #| label: correlation-matrix vars <- data[, c("euclid.jalan", "euclid.desa", "euclid.sungai", "euclid.hutan", "Slope", "Elevasi", "NDVI")] correlation_matrix <- cor(vars, method = "pearson") print(correlation_matrix) ``` ```{r} #| label: fig-corrplot #| fig-cap: "Matriks Korelasi Pearson antar variabel prediktor" #| fig-width: 8 #| fig-height: 6 corrplot( correlation_matrix, method = "number", type = "upper", tl.col = "black", tl.srt = 45, tl.cex = 1.2, cl.pos = "n", number.cex = 1.2 ) ``` --- # Perbandingan Model Alternatif Jika `Slope` dan `Elevasi` berkorelasi tinggi, kita bandingkan dua model yang masing-masing hanya mengandung satu dari keduanya. Model dengan **AIC lebih rendah** adalah yang lebih baik. ```{r} #| label: model-slope model_slope <- glm( PA ~ `euclid.jalan` + `euclid.desa` + `euclid.sungai` + `euclid.hutan` + Slope + NDVI, data = data, family = binomial ) summary(model_slope) ``` ```{r} #| label: model-elevasi model_elevasi <- glm( PA ~ `euclid.jalan` + `euclid.desa` + `euclid.sungai` + `euclid.hutan` + Elevasi + NDVI, data = data, family = binomial ) summary(model_elevasi) ``` ```{r} #| label: aic-comparison cat("=== Perbandingan AIC ===\n") print(AIC(model_slope, model_elevasi)) ``` --- # Evaluasi Model Terbaik Menggunakan `model_elevasi` sebagai model terbaik. Ganti dengan `model_slope` jika AIC-nya lebih rendah. ```{r} #| label: best-model best_model <- model_elevasi # Ganti sesuai model yang dipilih ``` ## Kurva ROC dan AUC Kurva ROC mengukur kemampuan diskriminasi model. Nilai **AUC > 0.7** menunjukkan performa yang dapat diterima; **AUC > 0.9** sangat baik. ```{r} #| label: fig-roc #| fig-cap: "Kurva ROC model terbaik" #| fig-width: 6 #| fig-height: 6 roc_obj <- roc(data$PA, fitted(best_model)) plot(roc_obj, col = "steelblue", main = paste("AUC =", round(auc(roc_obj), 3))) ``` ## Hosmer–Lemeshow Goodness-of-Fit **p-value > 0.05** menunjukkan model fit yang baik (tidak ada perbedaan signifikan antara nilai observasi dan prediksi). ```{r} #| label: hosmer-lemeshow cat("=== Hosmer-Lemeshow ===\n") print(hoslem.test(data$PA, fitted(best_model), g = 10)) ``` ## Pseudo-R² (McFadden & Nagelkerke) Mengukur kekuatan penjelas model. Nagelkerke R² **> 0.2** umumnya dianggap dapat diterima. ```{r} #| label: pseudo-r2 cat("=== Nagelkerke R² ===\n") print(PseudoR2(best_model, c("McFadden", "Nagel"))) ``` ## Signifikansi Koefisien Koefisien estimasi, standard error, z-value, dan p-value untuk setiap prediktor. ```{r} #| label: coefficients cat("=== Signifikansi tiap variabel ===\n") print(summary(best_model)$coefficients) ``` --- # Informasi Sesi ```{r} #| label: session-info #| eval: true sessionInfo() ``` --- *Dokumen ini dibuat dengan [Quarto](https://quarto.org/) dan dapat direproduksi dengan menjalankan ulang seluruh chunk menggunakan data yang sesuai.*