library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
library(MASS)
library(readr)
## Warning: package 'readr' was built under R version 4.3.3
data_karhutla <- read_csv("D:/download/data_tubes_statling - Sheet1.csv")
## Rows: 14 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Provinsi
## dbl (4): Jumlah_kejadian, Suhu, Luas_Lahan, Kepadatan_Penduduk
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Tampilkan data
print(colnames(data_karhutla))
## [1] "Provinsi" "Jumlah_kejadian" "Suhu"
## [4] "Luas_Lahan" "Kepadatan_Penduduk"
print(head(data_karhutla))
## # A tibble: 6 × 5
## Provinsi Jumlah_kejadian Suhu Luas_Lahan Kepadatan_Penduduk
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 BARITO SELATAN 19 26.5 10495. 135966
## 2 BARITO TIMUR 1 26.5 1245. 116341
## 3 BARITO UTARA 2 26.1 832. 159735
## 4 GUNUNG MAS 2 25.9 1006. 130241
## 5 KAPUAS 17 26.1 38062. 415210
## 6 KATINGAN 13 27.2 8946. 177106
# Standarisasi Variabel X
data_karhutla$X1_Std <- scale(data_karhutla$Suhu)
data_karhutla$X2_Std <- scale(data_karhutla$Luas_Lahan)
data_karhutla$X3_Std <- scale(data_karhutla$Kepadatan_Penduduk)
#-------------------------------------------------------------------------------
#Uji Multikolinearitas (VIF)
# Hitung VIF
model_vif <- lm(Jumlah_kejadian ~ X1_Std + X2_Std + X3_Std, data = data_karhutla)
vif_result <- vif(model_vif)
# Tampilkan hasil VIF
print("=== Hasil Uji Multikolinearitas (VIF) ===")
## [1] "=== Hasil Uji Multikolinearitas (VIF) ==="
print(vif_result)
## X1_Std X2_Std X3_Std
## 1.198497 1.135465 1.202865
#-------------------------------------------------------------------------------
#Estimasi Model Regresi Poisson (Model Awal)
# Estimasi Model Regresi Poisson
model_poisson <- glm(Jumlah_kejadian ~ X1_Std + X2_Std + X3_Std,
family = poisson(link = "log"),
data = data_karhutla)
# Tampilkan Ringkasan Model Poisson
print("=== Hasil Estimasi Regresi Poisson ===")
## [1] "=== Hasil Estimasi Regresi Poisson ==="
summary(model_poisson)
##
## Call:
## glm(formula = Jumlah_kejadian ~ X1_Std + X2_Std + X3_Std, family = poisson(link = "log"),
## data = data_karhutla)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.35789 0.08480 27.807 < 2e-16 ***
## X1_Std 0.15854 0.09399 1.687 0.091665 .
## X2_Std 0.26039 0.07339 3.548 0.000388 ***
## X3_Std 0.08885 0.07787 1.141 0.253878
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 98.103 on 13 degrees of freedom
## Residual deviance: 74.383 on 10 degrees of freedom
## AIC: 135.5
##
## Number of Fisher Scoring iterations: 5
# Hitung Incidence Rate Ratio (IRR)
print("Incidence Rate Ratios (IRR) Model Poisson:")
## [1] "Incidence Rate Ratios (IRR) Model Poisson:"
print(exp(coef(model_poisson)))
## (Intercept) X1_Std X2_Std X3_Std
## 10.568614 1.171797 1.297432 1.092919
#-------------------------------------------------------------------------------
#Uji Overdispersion
deviance_ratio <- model_poisson$deviance / model_poisson$df.residual
print(paste("Residual Deviance / DF (Rasio Dispersi):", round(deviance_ratio, 4)))
## [1] "Residual Deviance / DF (Rasio Dispersi): 7.4383"
#-------------------------------------------------------------------------------
#Regresi Binomial Negatif
# Estimasi Model Regresi Binomial Negatif (Koreksi Overdispersion)
model_nb <- glm.nb(Jumlah_kejadian ~ X1_Std + X2_Std + X3_Std, data = data_karhutla)
# summary Model Binomial Negatif
print("=== Hasil Estimasi Regresi Binomial Negatif ===")
## [1] "=== Hasil Estimasi Regresi Binomial Negatif ==="
summary(model_nb)
##
## Call:
## glm.nb(formula = Jumlah_kejadian ~ X1_Std + X2_Std + X3_Std,
## data = data_karhutla, init.theta = 1.832699577, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.32270 0.21686 10.711 <2e-16 ***
## X1_Std 0.19046 0.24767 0.769 0.4419
## X2_Std 0.48149 0.23000 2.093 0.0363 *
## X3_Std 0.08316 0.23787 0.350 0.7266
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8327) family taken to be 1)
##
## Null deviance: 19.819 on 13 degrees of freedom
## Residual deviance: 15.051 on 10 degrees of freedom
## AIC: 102.64
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.833
## Std. Err.: 0.849
##
## 2 x log-likelihood: -92.638
# Hitung Incidence Rate Ratio (IRR) Model Binomial Negatif
print("Incidence Rate Ratios (IRR) Model Binomial Negatif:")
## [1] "Incidence Rate Ratios (IRR) Model Binomial Negatif:"
print(exp(coef(model_nb)))
## (Intercept) X1_Std X2_Std X3_Std
## 10.203210 1.209811 1.618486 1.086721