Pertemuan 9 ADK
aldi8402
2024-10-09
MODEL REGRESI LOGISTIK ORDINAL
Tujuan
- Mempelajari program R untuk membuat model regresi logistik ordinal.
- Mengaplikasikan program R untuk membuat model regresi logistik
ordinal.
Pendahuluan
Regresi logistik ordinal digunakan untuk memodelkan hubungan antara peubah respon yang berskala ordinal dengan peubah peubah penjelasnya. Jika diasumsikan terdapat peubah respon Y berskala ordinal dengan J kategori dan \(x'\)=(\(x_1\),\(x_2\),…,\(x_p\) ) adalah vektor peubah penjelas, maka peluang dari peubah respon kategori ke-j pada peubah penjelas X tertentu dapat dinyatakan dengan \(P[Y=j| x]\)\(=\)\(π_j\)(x) dan peluang kumulatifnya adalah (Hosmer & Lemeshow 2000):
Metode yang digunakan untuk mengestimasi parameter model regresi logistik ordinal pada penulisan ini adalah metode maksimum likelihood (maximum likelihood methods). Estimasi parameter model regresi logistik ordinal menggunakan algoritma Newton Raphson. Kemudian setelah diperoleh estimasi parameter, dilakukan uji taraf nyata parameter menggunakan Uji rasio likelihood dan uji Wald.
Program R
Contoh kasus:
Untuk analisis data di bawah ini, kita akan
mengembangkan contoh 3 tentang pendaftaran ke sekolah pascasarjana.
Disimulasikan beberapa data untuk contoh ini, data dapat diperoleh dari
situs web:
Panggil pagkage
library(foreign)
library(ggplot2)
library(MASS)
#install.packages("Hmisc")
library(Hmisc)##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(reshape2)
dat <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")
head(dat)## apply pared public gpa
## 1 very likely 0 0 3.26
## 2 somewhat likely 1 0 3.21
## 3 unlikely 1 1 3.94
## 4 somewhat likely 0 0 2.81
## 5 somewhat likely 0 0 2.53
## 6 unlikely 0 1 2.59
Set data hipotetis ini memiliki variabel tiga tingkat yang disebut apply, dengan tingkat “unlikely”, “somewhat likely”, dan “very likely”, masing-masing diberi kode 1, 2, dan 3, yang akan digunakan sebagai variabel hasil. Juga terdapat tiga variabel yang akan digunakan sebagai prediktor: pared, yang merupakan variabel 0/1 yang menunjukkan apakah setidaknya salah satu orang tua memiliki gelar sarjana; public, yang merupakan variabel 0/1 di mana 1 menunjukkan bahwa lembaga pendidikan tinggi tersebut adalah negeri dan 0 swasta, dan GPA, yang merupakan indeks prestasi kumulatif siswa. Selanjutnya dalah melihat statistic deskriptif dari variabel-variabel ini.
## one at a time, table apply, pared, and public
lapply(dat[, c("apply", "pared", "public")], table)## $apply
##
## unlikely somewhat likely very likely
## 220 140 40
##
## $pared
##
## 0 1
## 337 63
##
## $public
##
## 0 1
## 343 57## three way cross tabs (xtabs) and flatten the table
ftable(xtabs(~ public + apply + pared, data = dat))## pared 0 1
## public apply
## 0 unlikely 175 14
## somewhat likely 98 26
## very likely 20 10
## 1 unlikely 25 6
## somewhat likely 12 4
## very likely 7 3summary(dat$gpa)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.900 2.720 2.990 2.999 3.270 4.000sd(dat$gpa)## [1] 0.3979409
Kita juga dapat memeriksa distribusi GPA di setiap tingkat apply yang dipecah berdasarkan public dan pared. Hal ini akan menciptakan grid 2 x 2 dengan boxplot GPA untuk setiap level pelamar, untuk nilai tertentu dari public dan pared. Untuk melihat data dengan lebih baik, ditambahkan titik data mentah di atas boxplot, dengan sedikit noise (sering disebut “jitter”) dan transparansi 50% sehingga tidak memenuhi boxplot. Terakhir, selain sel, kita memplot semua hubungan marjinal. Margin membuat plot akhir menjadi kotak 3 x 3. Di pojok kanan bawah, adalah hubungan keseluruhan antara apply dan GPA yang terlihat sedikit positif. Untuk melakukan hal ini, kita menggunakan package ggplot2.
ggplot(dat, aes(x = apply, y = gpa)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
facet_grid(pared ~ public, margins = TRUE) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
Untuk memahami bagaimana menginterpretasikan koefisien, pertama-tama mari kita membuat beberapa notasi dan meninjau konsep-konsep yang terlibat dalam regresi logistik ordinal. Misalkan \(Y\) menjadi hasil ordinal dengan \(J\) kategori. Kemudian \(P(Y≤j)\)adalah probabilitas kumulatif dari \(Y\) kurang dari atau sama dengan kategori tertentu \(j=1,…,j-1\). Peluang kurang dari atau sama dengan kategori tertentu dapat didefinisikan sebagai \(\frac{P(Y≤j)}{P(Y>j)}\) untuk \(j=1,…,j-1\) karena \(P(Y>j)\) dan membagi dengan nol tidak didefinisikan. Log odds juga dikenal sebagai logit, sehingga $log \[\frac{P(Y≤j)}{P(Y>j)}\] =logit( P(Y≤j))$ Dalam polr R, model regresi logistik ordinal diparameterkan sebagai \(logit(P(Y≤j))=β_j0-η_1 x_1-…-η_p x_p\)
## fit ordered logit model and store results 'm'
m <- polr(apply ~ pared + public + gpa, data = dat, Hess=TRUE)
## view a summary of the model
summary(m)## Call:
## polr(formula = apply ~ pared + public + gpa, data = dat, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## pared 1.04769 0.2658 3.9418
## public -0.05879 0.2979 -0.1974
## gpa 0.61594 0.2606 2.3632
##
## Intercepts:
## Value Std. Error t value
## unlikely|somewhat likely 2.2039 0.7795 2.8272
## somewhat likely|very likely 4.2994 0.8043 5.3453
##
## Residual Deviance: 717.0249
## AIC: 727.0249Model estimasi dapat ditulis sebagai:
# store table
(ctable <- coef(summary(m)))## Value Std. Error t value
## pared 1.04769010 0.2657894 3.9418050
## public -0.05878572 0.2978614 -0.1973593
## gpa 0.61594057 0.2606340 2.3632399
## unlikely|somewhat likely 2.20391473 0.7795455 2.8271792
## somewhat likely|very likely 4.29936315 0.8043267 5.3452947# calculate and store p values
p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
# combined table
(ctable <- cbind(ctable, "p value" = p))## Value Std. Error t value p value
## pared 1.04769010 0.2657894 3.9418050 8.087072e-05
## public -0.05878572 0.2978614 -0.1973593 8.435464e-01
## gpa 0.61594057 0.2606340 2.3632399 1.811594e-02
## unlikely|somewhat likely 2.20391473 0.7795455 2.8271792 4.696004e-03
## somewhat likely|very likely 4.29936315 0.8043267 5.3452947 9.027008e-08
Kita juga bisa mendapatkan interval kepercayaan untuk estimasi parameter. Ini dapat diperoleh dengan membuat profil fungsi likelihood atau dengan menggunakan kesalahan standar dan mengasumsikan distribusi normal. Perhatikan bahwa profil CI tidak simetris (meskipun biasanya mendekati simetris). Jika 95% CI tidak melewati 0, maka estimasi parameter signifikan secara statistik.
(ci <- confint(m)) # default method gives profiled CIs## Waiting for profiling to be done...## 2.5 % 97.5 %
## pared 0.5281768 1.5721750
## public -0.6522060 0.5191384
## gpa 0.1076202 1.1309148confint.default(m) # CIs assuming normality## 2.5 % 97.5 %
## pared 0.5267524 1.5686278
## public -0.6425833 0.5250119
## gpa 0.1051074 1.1267737
CI untuk pared dan GPA tidak
menyertakan 0; public menyertakan 0. Estimasi dalam
output diberikan dalam satuan logit terurut, atau peluang log odds. Jadi
untuk pared, dapat dinyatakan bahwa untuk peningkatan
satu unit dalam pared (yaitu, dari 0 ke 1), kami
mengharapkan peningkatan 1,05 dalam nilai yang diharapkan dari
apply pada skala odds log, mengingat semua variabel
lain dalam model dianggap konstan. Untuk GPA,
dinyatakan bahwa untuk peningkatan satu unit GPA, kita
mengharapkan peningkatan 0,62 dalam nilai yang diharapkan dari
apply pada skala odds log, mengingat semua variabel
lain dalam model dianggap konstan.
Koefisien dari model ini bisa
jadi agak sulit untuk ditafsirkan karena diskalakan dalam bentuk log.
Cara lain untuk menginterpretasikan model regresi logistik adalah dengan
mengubah koefisien menjadi rasio odds. Untuk mendapatkan OR dan interval
kepercayaan, kita hanya perlu mengeksponensikan estimasi dan interval
kepercayaan.
## odds ratios
exp(coef(m))## pared public gpa
## 2.8510579 0.9429088 1.8513972## OR and CI
exp(cbind(OR = coef(m), ci))## OR 2.5 % 97.5 %
## pared 2.8510579 1.6958376 4.817114
## public 0.9429088 0.5208954 1.680579
## gpa 1.8513972 1.1136247 3.098490Koefisien ini disebut rasio odds proporsional dan kita akan
menginterpretasikannya sama seperti kita menginterpretasikan rasio odds
dari regresi logistik biner.
Perintah terakhir meminta R untuk
mengembalikan konten ke objek s, yang merupakan sebuah tabel.
sf <- function(y) { c('Y>=1' = qlogis(mean(y >= 1)), 'Y>=2' = qlogis(mean(y >= 2)), 'Y>=3' = qlogis(mean(y >= 3))) }
(s <- with(dat, summary(as.numeric(apply) ~ pared + public + gpa, fun=sf)))## as.numeric(apply) N= 400
##
## +-------+-----------+---+----+-----------+---------+
## | | | N|Y>=1| Y>=2| Y>=3|
## +-------+-----------+---+----+-----------+---------+
## | pared| No|337| Inf|-0.37833644|-2.440735|
## | | Yes| 63| Inf| 0.76546784|-1.347074|
## +-------+-----------+---+----+-----------+---------+
## | public| No|343| Inf|-0.20479441|-2.345006|
## | | Yes| 57| Inf|-0.17589067|-1.547563|
## +-------+-----------+---+----+-----------+---------+
## | gpa|[1.90,2.73)|102| Inf|-0.39730180|-2.772589|
## | |[2.73,3.00)| 99| Inf|-0.26415158|-2.302585|
## | |[3.00,3.28)|100| Inf|-0.20067070|-2.090741|
## | |[3.28,4.00]| 99| Inf| 0.06062462|-1.803594|
## +-------+-----------+---+----+-----------+---------+
## |Overall| |400| Inf|-0.20067070|-2.197225|
## +-------+-----------+---+----+-----------+---------+glm(I(as.numeric(apply) >= 2) ~ pared, family="binomial", data = dat)##
## Call: glm(formula = I(as.numeric(apply) >= 2) ~ pared, family = "binomial",
## data = dat)
##
## Coefficients:
## (Intercept) pared
## -0.3783 1.1438
##
## Degrees of Freedom: 399 Total (i.e. Null); 398 Residual
## Null Deviance: 550.5
## Residual Deviance: 534.1 AIC: 538.1glm(I(as.numeric(apply) >= 3) ~ pared, family="binomial", data = dat)##
## Call: glm(formula = I(as.numeric(apply) >= 3) ~ pared, family = "binomial",
## data = dat)
##
## Coefficients:
## (Intercept) pared
## -2.441 1.094
##
## Degrees of Freedom: 399 Total (i.e. Null); 398 Residual
## Null Deviance: 260.1
## Residual Deviance: 252.2 AIC: 256.2s[, 4] <- s[, 4] - s[, 3]
s[, 3] <- s[, 3] - s[, 3]
s # print## as.numeric(apply) N= 400
##
## +-------+-----------+---+----+----+---------+
## | | | N|Y>=1|Y>=2| Y>=3|
## +-------+-----------+---+----+----+---------+
## | pared| No|337| Inf| 0|-2.062399|
## | | Yes| 63| Inf| 0|-2.112541|
## +-------+-----------+---+----+----+---------+
## | public| No|343| Inf| 0|-2.140211|
## | | Yes| 57| Inf| 0|-1.371672|
## +-------+-----------+---+----+----+---------+
## | gpa|[1.90,2.73)|102| Inf| 0|-2.375287|
## | |[2.73,3.00)| 99| Inf| 0|-2.038434|
## | |[3.00,3.28)|100| Inf| 0|-1.890070|
## | |[3.28,4.00]| 99| Inf| 0|-1.864219|
## +-------+-----------+---+----+----+---------+
## |Overall| |400| Inf| 0|-1.996554|
## +-------+-----------+---+----+----+---------+plot(s, which=1:3, pch=1:3, xlab='logit', main=' ', xlim=range(s[,3:4]))
newdat <- data.frame(
pared = rep(0:1, 200),
public = rep(0:1, each = 200),
gpa = rep(seq(from = 1.9, to = 4, length.out = 100), 4))
newdat <- cbind(newdat, predict(m, newdat, type = "probs"))
##show first few rows
head(newdat)## pared public gpa unlikely somewhat likely very likely
## 1 0 0 1.900000 0.7376186 0.2204577 0.04192370
## 2 1 0 1.921212 0.4932185 0.3945673 0.11221424
## 3 0 0 1.942424 0.7325300 0.2244841 0.04298593
## 4 1 0 1.963636 0.4866885 0.3984676 0.11484395
## 5 0 0 1.984848 0.7273792 0.2285470 0.04407383
## 6 1 0 2.006061 0.4801630 0.4023098 0.11752712
Sekarang kita dapat membentuk ulang data yang panjang dengan paket reshape2 dan memplot semua probabilitas yang diprediksi untuk kondisi yang berbeda. Kita memplot probabilitas yang diprediksi, dihubungkan dengan sebuah garis, diwarnai berdasarkan tingkat hasil, berlaku, dan di-faceted berdasarkan tingkat pared dan publik. Kita juga menggunakan fungsi label khusus, untuk menambahkan label yang lebih jelas yang menunjukkan apa yang diwakili oleh setiap kolom dan baris plot.
lnewdat <- melt(newdat, id.vars = c("pared", "public", "gpa"),
variable.name = "Level", value.name="Probability")
## view first few rows
head(lnewdat)## pared public gpa Level Probability
## 1 0 0 1.900000 unlikely 0.7376186
## 2 1 0 1.921212 unlikely 0.4932185
## 3 0 0 1.942424 unlikely 0.7325300
## 4 1 0 1.963636 unlikely 0.4866885
## 5 0 0 1.984848 unlikely 0.7273792
## 6 1 0 2.006061 unlikely 0.4801630ggplot(lnewdat, aes(x = gpa, y = Probability, colour = Level)) +
geom_line() + facet_grid(pared ~ public, labeller="label_both")