Regresi Logistik dalam R (Odds Ratio)

41

Saya mencoba melakukan analisis regresi logistik di R. Saya telah mengikuti kursus yang membahas materi ini menggunakan STATA. Saya merasa sangat sulit untuk meniru fungsi di R. Apakah sudah matang di area ini? Tampaknya ada sedikit dokumentasi atau panduan yang tersedia. Memproduksi output odds ratio tampaknya membutuhkan instalasi epicalcdan / atau epitoolsdan / atau yang lainnya, yang tidak ada yang bisa saya gunakan, sudah ketinggalan zaman atau kurang dokumentasi. Saya sudah terbiasa glmmelakukan regresi logistik. Setiap saran akan diterima.

Saya lebih baik membuat ini pertanyaan nyata. Bagaimana cara menjalankan regresi logistik dan menghasilkan rasio odds R?

Inilah yang saya lakukan untuk analisis univariat:

x = glm(Outcome ~ Age, family=binomial(link="logit"))

Dan untuk multivarian:

y = glm(Outcome ~ Age + B + C, family=binomial(link="logit"))

Saya telah lalu memandang x, y, summary(x)dan summary(y).

Apakah x$coefficientsada nilainya?

r logistic odds-ratio SabreWolfy
sumber

Jawaban:

36

jika Anda ingin menginterpretasikan efek yang diperkirakan sebagai rasio odds relatif, lakukan saja exp(coef(x))(memberi Anda , perubahan multiplikasi dalam rasio odds untuk jika kovariat yang terkait dengan meningkat sebesar 1). Untuk interval kemungkinan profil untuk jumlah ini, Anda dapat melakukannya y = 1 βeβy=1β

require(MASS)
exp(cbind(coef(x), confint(x)))

Sunting: @caracal lebih cepat ...

orang fabian
sumber
1
+1 untuk saran @ fabian. Cara lebih rendah untuk melakukan hal ini yang biasanya menghasilkan interval yang sama adalah untuk menghitung interval pada skala logit dan kemudian mengubah dengan skala kemungkinan: cbind( exp(coef(x)), exp(summary(x)$coefficients[,1] - 1.96*summary(x)$coefficients[,2]), exp(summary(x)$coefficients[,1] + 1.96*summary(x)$coefficients[,2]) ). Ada juga metode delta: ats.ucla.edu/stat/r/faq/deltamethod.htm
terkunci
42

Anda benar bahwa output R biasanya hanya berisi informasi penting, dan lebih banyak yang perlu dihitung secara terpisah.

N  <- 100               # generate some data
X1 <- rnorm(N, 175, 7)
X2 <- rnorm(N,  30, 8)
X3 <- abs(rnorm(N, 60, 30))
Y  <- 0.5*X1 - 0.3*X2 - 0.4*X3 + 10 + rnorm(N, 0, 12)

# dichotomize Y and do logistic regression
Yfac   <- cut(Y, breaks=c(-Inf, median(Y), Inf), labels=c("lo", "hi"))
glmFit <- glm(Yfac ~ X1 + X2 + X3, family=binomial(link="logit"))

coefficients()memberi Anda estimasi parameter regresi . Lebih mudah untuk menafsirkan (kecuali untuk intersep).bjexhal(bj)

> exp(coefficients(glmFit))
 (Intercept)           X1           X2           X3 
5.811655e-06 1.098665e+00 9.511785e-01 9.528930e-01

Untuk mendapatkan odds ratio, kita perlu tabel silang klasifikasi DV dikotomus asli dan klasifikasi yang diprediksi berdasarkan beberapa ambang batas probabilitas yang perlu dipilih terlebih dahulu. Anda juga dapat melihat fungsi ClassLog()dalam paket QuantPsyc(seperti chl disebutkan dalam pertanyaan terkait ).

# predicted probabilities or: predict(glmFit, type="response")
> Yhat    <- fitted(glmFit)
> thresh  <- 0.5  # threshold for dichotomizing according to predicted probability
> YhatFac <- cut(Yhat, breaks=c(-Inf, thresh, Inf), labels=c("lo", "hi"))
> cTab    <- table(Yfac, YhatFac)    # contingency table
> addmargins(cTab)                   # marginal sums
     YhatFac
Yfac   lo  hi Sum
  lo   41   9  50
  hi   14  36  50
  Sum  55  45 100

> sum(diag(cTab)) / sum(cTab)        # percentage correct for training data
[1] 0.77

Untuk rasio odds, Anda dapat menggunakan paket vcdatau melakukan perhitungan secara manual.

> library(vcd)                       # for oddsratio()
> (OR <- oddsratio(cTab, log=FALSE)) # odds ratio
[1] 11.71429

> (cTab[1, 1] / cTab[1, 2]) / (cTab[2, 1] / cTab[2, 2])
[1] 11.71429

> summary(glmFit)  # test for regression parameters ...

# test for the full model against the 0-model
> glm0 <- glm(Yfac ~ 1, family=binomial(link="logit"))
> anova(glm0, glmFit, test="Chisq")
Analysis of Deviance Table
Model 1: Yfac ~ 1
Model 2: Yfac ~ X1 + X2 + X3
  Resid. Df Resid. Dev Df Deviance P(>|Chi|)    
1        99     138.63                          
2        96     110.58  3   28.045 3.554e-06 ***
caracal
sumber
2
Terima kasih - saya harus memeriksa jawaban Anda dengan cermat. Di STATA kita bisa menjalankan logitdan logisticdan mendapatkan odds rasio dan interval kepercayaan dengan mudah. Saya agak frustrasi karena ini tampaknya sangat rumit dan tidak standar R. Bisakah saya menggunakan exp(cbind(coef(x), confint(x)))jawaban orang fabian di bawah ini untuk mendapatkan OD dan CI? Saya tidak jelas apa jawaban Anda?
SabreWolfy
3
cTabexhal(bj)exhal(bj)Xj
4
sebenarnya @SabreWolfy Saya merasa frustasi karena orang dapat mengklik satu tombol di stata / sas / spss dll, dan mendapatkan odds rasio (masukkan statistik kecocokan, tipe III SS, apa pun yang Anda suka di sini) tanpa memiliki petunjuk tentang apa artinya / bagaimana cara menghitungnya / apakah itu bermakna dalam situasi tertentu / dan (mungkin lebih penting) tanpa memiliki pengetahuan tentang bahasa itu sendiri.
rawr
5

Paket epiDisplay melakukan ini dengan sangat mudah.

library(epiDisplay)
data(Wells, package="carData")
glm1 <- glm(switch~arsenic+distance+education+association, 
            family=binomial, data=Wells)
logistic.display(glm1)
Logistic regression predicting switch : yes vs no 

                       crude OR(95%CI)         adj. OR(95%CI)         P(Wald's test) P(LR-test)
arsenic (cont. var.)   1.461 (1.355,1.576)     1.595 (1.47,1.731)     < 0.001        < 0.001   

distance (cont. var.)  0.9938 (0.9919,0.9957)  0.9911 (0.989,0.9931)  < 0.001        < 0.001   

education (cont. var.) 1.04 (1.021,1.059)      1.043 (1.024,1.063)    < 0.001        < 0.001   

association: yes vs no 0.863 (0.746,0.999)     0.883 (0.759,1.027)    0.1063         0.1064    

Log-likelihood = -1953.91299
No. of observations = 3020
AIC value = 3917.82598
Edward
sumber
Apakah ada cara untuk menggabungkan tampilan logistik dengan pembungkus lateks seperti outregatau xtable?
Kesalahan Nama Terkemuka