Regresión logística en R (Odds Ratio)

Question

Estoy tratando de realizar un análisis de regresión logística en R . He asistido a cursos sobre este material utilizando STATA. Me resulta muy difícil reproducir la funcionalidad en R . ¿Está maduro en esta zona? Parece que hay poca documentación u orientación disponible. La producción de resultados de odds ratio parece requerir la instalación de epicalc y/o epitools y/u otros, ninguno de los cuales puedo hacer funcionar, están desactualizados o carecen de documentación. He utilizado glm para hacer la regresión logística. Cualquier sugerencia será bienvenida.

Será mejor que haga una pregunta real. ¿Cómo puedo ejecutar una regresión logística y producir raciones de probabilidades en R ?

Esto es lo que he hecho para un análisis univariante:

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

Y para las multivariantes:

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

Luego he mirado x , y , summary(x) y summary(y) .

Es x$coefficients de algún valor?

Answer 1

Si quiere interpretar los efectos estimados como odds ratios relativas, simplemente haga exp(coef(x)) (te da $e^\beta$ el cambio multiplicativo de la razón de momios para $y=1$ si la covariable asociada a $\beta$ aumenta en 1). Para los intervalos de probabilidad del perfil de esta cantidad, se puede hacer

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

EDITAR: @caracal fue más rápido...

Answer 2

Tienes razón en que la salida de R suele contener sólo la información esencial, y hay que calcular más por separado.

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() le da los parámetros de regresión estimados $b_{j}$ . Es más fácil interpretar $exp(b_{j})$ sin embargo (a excepción de la intercepción).

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

Para obtener la razón de probabilidades, necesitamos la tabla cruzada de clasificación de la VD dicotómica original y la clasificación predicha según algún umbral de probabilidad que hay que elegir primero. También se puede ver la función ClassLog() en el paquete QuantPsyc (como chl mencionó en una pregunta relacionada ).

# 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

Para la proporción de probabilidades, puede utilizar el paquete vcd o hacer el cálculo manualmente.

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

Answer 3

La página de estadísticas de la UCLA tiene un buen recorrido de realizar una regresión logística en R. Incluye una breve sección sobre el cálculo de odds ratios.

Answer 4

El paquete epiDisplay lo hace muy fácilmente.

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

Answer 5

R ha madurado con respecto a los cálculos de odds ratio más de dos décadas. Es mejor pensar en esto en términos generales. Por ejemplo, ¿qué pasa si x1 y x2 son continuos y tienen efectos no lineales e interactúan entre sí? Aquí hay un código de ejemplo en el que se calcula el efecto del rango intercuartil de x1, ajustado a x2=1,5.

require(rms)
dd <- datadist(mydata); options(datadist='dd')
# restricted cubic spline in x1, quadratic in x2, interactions
f <- lrm(y ~ rcs(x1, 4) * pol(x2, 2), data=mydata)
summary(f, x2=1.5)  # gives IQR odds ratio for x1
summary(f, x1=c(1, 3), x2=1.5)  # OR for x1=1 vs. x1=3, still x2=1.5

Puedes ver que tratar con coeficientes individuales no es la solución general.