[PADK] Regresi Logistik - Regresi Berganda

Angga Fathan Rofiqy

17 November, 2023

Kode di Hide dalam default, untuk menampilkan kode, klik Code .

#                      -=( Install & Load Package Function )=-
install_load <- function (package1, ...)  {   

   # convert arguments to vector
   packages <- c(package1, ...)

   # start loop to determine if each package is installed
   for(package in packages){

       # if package is installed locally, load
       if(package %in% rownames(installed.packages()))
          do.call('library', list(package))

       # if package is not installed locally, download, then load
       else {
          install.packages(package)
          do.call("library", list(package))
       }
   } 
}

#Path Function
path <- function(){
  gsub  ( "\\\\",  "/",  readClipboard ()  )
}
#Copy path, Panggil function di console
#Copy r path, paste ke var yang diinginkan

1 Pendahuluan

1.1 Kelompok 1

Nama NIM
Karimatu Ain G1401211001
Diva Nisfu Mustika G1401211002
Lutfi Syahreza Lubis G1401211003
Asfiah Adiba G1401211004
Angga Fathan Rofiqy G1401211006
Farhan Abdillah Harahap G1401211007
Nabil Naufal G1401211008
Muhammad Rizky Fajar G1401211009
Mutiara Andhini G1401211010
Muhammad Nafiz G1401211011
Syifa Khairunnisa G1401211012
Reynd Hamonangan Pasaribu G1401211013

1.2 Soal

1.3 Data Entry

install_load('dplyr', 'DT', 'tidyr', 'stringr')

# Mengubah Tabel Kontingensi ke dalam Data Frame R
data <- data.frame(
  City = rep(c("Beijing", "Shanghai", "Shenyang", "Nanjing", "Harbin", "Zehgzhou", "Taiyuan", "Nanchang"), each = 2),
  Smoking = rep(c("Yes", "No"), times = 8),
  LungCancer_Yes = c(126, 35, 908, 497, 913, 336, 235, 58, 402, 121, 182, 72, 60, 11, 104, 21),
  LungCancer_No = c(100, 61, 688, 807, 747, 598, 172, 121, 308, 215, 156, 98, 99, 43, 89, 36)
)
datatable(data)

1.4 Transform Data

data2 <- data %>%
  pivot_longer(cols = c("LungCancer_Yes", "LungCancer_No"), 
               names_to = "LungCancer", values_to = "Count") %>%
  mutate(LungCancer = ifelse(str_detect(LungCancer, "Yes"), "Yes", "No")) %>%
  uncount(weights = Count) %>%
  dplyr::select(City, Smoking, LungCancer) %>%
  mutate(
    City = as.factor(City),
    Smoking = as.factor(Smoking),
    LungCancer = as.factor(LungCancer)
  )
datatable(data2)

1.5 Data Checking

1.5.1 Tipe Data

str(data2)
## tibble [8,419 × 3] (S3: tbl_df/tbl/data.frame)
##  $ City      : Factor w/ 8 levels "Beijing","Harbin",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Smoking   : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ LungCancer: Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...

Agar dapat bisa dianalisis, data harus bertipe Factor. Terlihat bahwa tipe data di atas sudah sesuai.

1.5.2 Banyaknya data transform

data_restored <- data2 %>%
  group_by(City, Smoking, LungCancer, .drop = TRUE) %>%
  summarize(Count = n(), .groups = 'drop') %>%
  pivot_wider(names_from = LungCancer, values_from = Count) %>%
  ungroup()
datatable(data_restored)

Bisa di cek bahwa. Jumlah data yang di transform sudah benar.

2 Jawab Soal

2.1 No 1.

Hitunglah dugaan parameter regresi logistik berganda untuk Peubah Smoking

2.1.1 Rstudio

model <- glm(LungCancer ~ Smoking + City, data = data2, family = binomial)
summary(model)
## 
## Call:
## glm(formula = LungCancer ~ Smoking + City, family = binomial, 
##     data = data2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3001  -1.2637  -0.9444   1.0935   1.7530  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.548682   0.118022  -4.649 3.34e-06 ***
## SmokingYes    0.777062   0.046775  16.613  < 2e-16 ***
## CityHarbin    0.018187   0.129473   0.140    0.888    
## CityNanchang -0.054906   0.170996  -0.321    0.748    
## CityNanjing   0.005764   0.140911   0.041    0.967    
## CityShanghai  0.055618   0.119570   0.465    0.642    
## CityShenyang -0.027739   0.120071  -0.231    0.817    
## CityTaiyuan  -0.745683   0.185519  -4.019 5.83e-05 ***
## CityZehgzhou  0.028782   0.144755   0.199    0.842    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11663  on 8418  degrees of freedom
## Residual deviance: 11358  on 8410  degrees of freedom
## AIC: 11376
## 
## Number of Fisher Scoring iterations: 4

Bisa dilihat pada bagian bahwa nilai dugaan parameter regresi logistik berganda untuk Peubah Smoking adalah sekitar \(0.777062\).

2.1.2 SAS

Terlihat bahwa nilai yang diperoleh untuk dugaan parameter regresi logistik berganda pada Peubah Smoking adalah sebesar \(0.7704\).

2.1.3 Minitab

Terlihat bahwa nilai yang diperoleh untuk dugaan parameter regresi logistik berganda pada Peubah Smoking adalah sebesar \(0.7771\).

2.1.4 Perbandingan

no1 <- data.frame(R = summary(model)[["coefficients"]]["SmokingYes",1] %>%
                    round(.,4),
                  SAS = 0.7704, 
                  MiniTab = 0.7771)
datatable(no1)

Sehingga jawaban yang tepat adalah :

2.2 No 2.

Nilai odds, p-value, dan Statistik Uji Wald pada peubah Smoking

2.2.1 Rstudio

summary_model <- summary(model)
rasio.odds <- exp(coef(model)) %>% round(.,3)
p.value <- summary(model)$coefficients[, "Pr(>|z|)"] %>% 
  sprintf("%.2e", .)

datatable(data.frame(rasio.odds[1:2], p.value[1:2]))

Didapat bahwa nilai Rasio odds pada peubah Smoking diperoleh sebesar \(2.175\).

install_load('car')
anova <- Anova(model, type="II", test.statistic="Wald"); anova
## Analysis of Deviance Table (Type II tests)
## 
## Response: LungCancer
##         Df   Chisq Pr(>Chisq)    
## Smoking  1 275.988  < 2.2e-16 ***
## City     7  28.568  0.0001734 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Lalu nilai Statistik Uji Wald nya sebesar \(275.988\) dan p-value nya sangat kecil yakni kurang dari \(2.2\times10^{-16}\).

2.2.2 SAS

Didapatkan Nilai odds pada peubah Smoking sebesar \(2.161\).

Didapatkan Statistik Uji Wald pada peubah Smoking sebesar \(274.4587\). Dan p-value nya sangat kecil yakni kurang dari \(0.0001\).

2.2.3 Minitab

Didapatkan Nilai odds pada peubah Smoking sebesar \(2.1751\).

2.2.4 Perbandingan

no2 <- data.frame(R = c(rasio.odds[1:2][[2]], 
                        anova$Chisq[1] %>% round(.,4),
                        anova$`Pr(>Chisq)`[1] %>% sprintf("%.2e", .)),
                  SAS = c(2.161, 274.4587, "<0.0001"), 
                  MiniTab = c(2.1751, NA, NA))
rownames(no2) <- c("Rasio Odds", "Statistik Wald", "P-Value")
datatable(no2)

Sehingga jawaban yang tepat adalah

2.3 No 3.

Rasio Odds Lung Cancer

rasio.odds[1:2]
## (Intercept)  SmokingYes 
##       0.578       2.175

Bisa dilihat bahwa berdasarkan nilai rasio odds nya, kelompok yang memiliki kecenderungan lebih tinggi terkena kanker paru-paru adalah kelompok yang merokok. Yakni \(2.175\) kali lebih besar dari pada yang tidak merokok.

Pada dua sofware lainnya memiliki kesimpulan yang sama walaupun nilainya cukup berbeda. Yakni pada SAS sebesar \(2.161\) dan pada Minitab sebesar \(2.1751\).

Sehingga jawaban yang tepat adalah

2.4 No 4.

Kota yang memilikipengaruh yang signifikan

summary(model)
## 
## Call:
## glm(formula = LungCancer ~ Smoking + City, family = binomial, 
##     data = data2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3001  -1.2637  -0.9444   1.0935   1.7530  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.548682   0.118022  -4.649 3.34e-06 ***
## SmokingYes    0.777062   0.046775  16.613  < 2e-16 ***
## CityHarbin    0.018187   0.129473   0.140    0.888    
## CityNanchang -0.054906   0.170996  -0.321    0.748    
## CityNanjing   0.005764   0.140911   0.041    0.967    
## CityShanghai  0.055618   0.119570   0.465    0.642    
## CityShenyang -0.027739   0.120071  -0.231    0.817    
## CityTaiyuan  -0.745683   0.185519  -4.019 5.83e-05 ***
## CityZehgzhou  0.028782   0.144755   0.199    0.842    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11663  on 8418  degrees of freedom
## Residual deviance: 11358  on 8410  degrees of freedom
## AIC: 11376
## 
## Number of Fisher Scoring iterations: 4

Bisa dilihat pada bagian

Bahwa Kota yang memiliki pengaruh yang signifikan pada taraf nyata \(5\%\) adalah kota Taiyuan dan kota Beijing (Intercept, karena sebagai peubah referensi).

Untuk No 4, hanya bisa dilakukan di Rstudio saja.

Sehingga jawaban yang tepat adalah

2.5 No 5.

Lakukan Uji kebaikan suai Hosmer and Lemeshow terhadap model yang di tentukan. Tentukan khi-kuadrat, p-value, keputusan, dan kesimpulan.

2.5.1 Rstudio

Hipotesis

\(H_0\) : Model Cocok

\(H_1\) : Model Tidak cocok

install_load("ResourceSelection")
hoslem <- hoslem.test(model$y, fitted(model), g = 10); hoslem
## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  model$y, fitted(model)
## X-squared = 1.388, df = 8, p-value = 0.9944

Didapatkan Khi-kuadrat sebesar \(1.388\). Nilainya sangat kecil, sehingga dapat dinyatakan kurang dari \(0.001\). Didapat pula p-value yang besar yakni sebesar \(0.9944\). Karena p-value lebih besar dari taraf nyata \(5\%\), maka Tak Tolak \(H_0\). Artinya Model yang digunakan sudah cocok.

2.5.2 SAS

Didapatkan nilai khi-kuadrat sebesar \(1.6062\). Dan p-value sebesar \(0.9521\).

Hipotesis

\(H_0\) : Model Cocok

\(H_1\) : Model Tidak cocok

Karena p-value lebih besar dari taraf nyata \(5\%\), maka Tak Tolak \(H_0\). Artinya Model yang digunakan sudah cocok.

2.5.3 MiniTab

Didapatkan nilai khi-kuadrat sebesar \(0.70\). Dan p-value sebesar \(0.951\).

Hipotesis

\(H_0\) : Model Cocok

\(H_1\) : Model Tidak cocok

Karena p-value lebih besar dari taraf nyata \(5\%\), maka Tak Tolak \(H_0\). Artinya Model yang digunakan sudah cocok.

2.5.4 Perbandingan

no3 <- data.frame(R = c(hoslem[["statistic"]][[1]], 
                        hoslem[["parameter"]][["df"]], 
                        hoslem[["p.value"]][[1]]) %>% round(.,4),
                  SAS = c(1.6062, 6, 0.9521), 
                  MiniTab = c(0.70, 4, 0.951))
rownames(no3) <- c("khi-kuadrat", "Db", "P-Value")
datatable(no3)

Terdapat perbedaan cukup jauh pada nilai Khi-kuadrat nya. Namun itu juga bisa disebabkan oleh perbedaan Derajat bebasnya. Walaupun begitu p-value nya tidak berbeda jauh. Sehingga keputusan dan kesimpulannya tetap sama.

Karena ada dua software yang memiliki jawaban yang sama. Maka jawaban yang tepat adalah

2.6 No 6.

Jika ingin diuji apakah model dengan kovariat city dan smoking (L1) lebih baik daripada model dengan hanya intersep saja (L0), maka dilakukan LR test.

Hipotesis

\(H_0\) : City dan Smoking tidak berpengaruh terhadap model

\(H_1\) : City dan Smoking berpengaruh terhadap model

LR test

# Model hanya intersep (L0)
model.L0 <- glm(LungCancer ~ 1, data = data2, family = binomial)

# Model dengan City dan Smoking (L1)
model.L1 <- glm(LungCancer ~ City + Smoking, data = data2, family = binomial)

install_load("lmtest")
# Uji LR
lrtest(model.L0, model.L1)
## Likelihood ratio test
## 
## Model 1: LungCancer ~ 1
## Model 2: LungCancer ~ City + Smoking
##   #Df  LogLik Df Chisq Pr(>Chisq)    
## 1   1 -5831.7                        
## 2   9 -5678.8  8 305.7  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Terlihat bahwa p-value bernilai \(2.2\times10^{-16}\). Nilainya sangat kecil, sehingga dapat dinyatakan kurang dari \(0.001\). Sehingga Tolak \(H_0\), artinya City dan Smoking berpengaruh terhadap model.

Nilai L0, L1, L0-L1

# -2log L0
minus2logL0 <- -2 * logLik(model.L0)

# -2log L1
minus2logL1 <- -2 * logLik(model.L1)

# -2log (L0-L1)
minus2logL0.L1 <- minus2logL0 - minus2logL1

log <- data.frame(minus2logL0, minus2logL1, minus2logL0.L1) %>% round(.,2)
colnames(log) <- c("-2log L0", "-2log L1", "-2log (L0-L1)")
datatable(log)

Diperoleh nilai \(log{(\text{L}_0^{-2})}\) adalah sebesar \(11663.37\). Lalu \(log{(\text{L}_1^{-2})}\) bernilai \(11357.67\). Dan \(log{(\text{L}_0-\text{L}_1)^{-2}}\) sebesar \(305.7\).

Untuk No 6, hanya bisa dilakukan di Rstudio saja.

Sehingga jawaban yang tepat adalah

3 Pengerjaan di Software Lain

3.1 SAS

File.sas

3.1.1 Code

DATA tugas; 
INPUT city smoking lungcancer count;
DATALINES;
0 1 1 126
0 1 0 100
0 0 1 35
0 0 0 61
1 1 1 908
1 1 0 688
1 0 1 497
1 0 0 807
2 1 1 913
2 1 0 747
2 0 1 336
2 0 0 598
3 1 1 235
3 1 0 172
3 0 1 58
3 0 0 121
4 1 1 402
4 1 0 308
4 0 1 121
4 0 0 215
5 1 1 182
5 1 0 156
5 0 1 72
5 0 0 98
6 1 1 60
6 1 0 99
6 0 1 11
6 0 0 43
7 1 1 104
7 1 0 89
7 0 1 21
7 0 0 36
;
RUN;

PROC LOGISTIC DATA = tugas descending; weight count;
MODEL lungcancer = city smoking/rsquare lackfit;
RUN;

3.1.2 Output

3.2 Minitab

File.mpx

3.2.1 Import data

Dimana :

  • Cancer = LungCancer

  • 1 = Yes. 0 = No

  • Trial = LungCancer_Yes + LungCancer_No

  • Event = LungCancer_Yes

3.2.2 Steps

Pastikan selain City peubah lainnya adalah numerik. Bisa dibuah di excel atau di minitab nya langsung. Cara ubah di minitab

  1. Pada Menu, pilih Stat ➡️ Regression ➡️Binary Logistic Regression ➡️ Fit Binary Logistic Model

  2. Memasukkan peubah-peubah seperti pada gambar

  3. Pada bagian Result pastikan tampilannya seperti ini (Bagian Display of results dipilih yang Expanded tables)

3.2.3 Output