Logistische Regression in R (Odds Ratio)

Ich versuche eine logistische Regressionsanalyse durchzuführen R. Ich habe Kurse besucht, die dieses Material mit STATA behandeln. Ich finde es sehr schwierig, die Funktionalität in zu replizieren R. Ist es in diesem Bereich ausgereift? Es scheint wenig Dokumentation oder Anleitung zu geben. Die Erstellung der Odds Ratio-Ausgabe erfordert anscheinend die Installation epicalcund / oder epitoolsund / oder andere, von denen keine zur Arbeit kommen kann, veraltet ist oder keine Dokumentation vorliegt. Ich habe glmdie logistische Regression durchgeführt. Anregungen wären willkommen.

Ich mache das besser zu einer echten Frage. Wie führe ich eine logistische Regression durch und erzeuge Quotenrationen in R?

Folgendes habe ich für eine univariate Analyse getan:

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

Und für multivariate:

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

Ich habe dann sieht x, y, summary(x)und summary(y).

Ist x$coefficientsirgendetwas wert?

Wenn Sie die geschätzten Effekte als relative Quotenverhältnisse interpretieren möchten, tun exp(coef(x))Sie einfach Folgendes (ergibt , die multiplikative Änderung des Quotenverhältnisses für wenn die mit verbundene Kovariate um 1 zunimmt). Für Profilwahrscheinlichkeitsintervalle für diese Menge können Sie Folgendes tun y = 1 $e^\beta$$y=1$$\beta$

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

EDIT: @caracal war schneller ...

Sie haben Recht, dass die Ausgabe von R normalerweise nur wichtige Informationen enthält und weitere Informationen separat berechnet werden müssen.

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()gibt Ihnen die geschätzten Regressionsparameter . Es ist jedoch einfacher, zu interpretieren (mit Ausnahme des Intercept).$b_{j}$$exp(b_{j})$

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

Um die Odds Ratio zu erhalten, benötigen wir die Klassifikationskreuztabelle des ursprünglichen dichotomen DV und die vorhergesagte Klassifikation gemäß einer Wahrscheinlichkeitsschwelle, die zuerst ausgewählt werden muss. Sie können die Funktion auch ClassLog()im Paket sehen QuantPsyc(wie in einer verwandten Frage erwähnt ).

# 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

Für das Odds Ratio können Sie entweder das Paket verwenden vcdoder die Berechnung manuell durchführen.

> 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 ***

Die UCLA-Statistikseite enthält eine ausführliche Beschreibung der Durchführung der logistischen Regression in R. Sie enthält einen kurzen Abschnitt zur Berechnung der Quotenverhältnisse.

Das epiDisplay-Paket erledigt dies sehr einfach.

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