Régression logistique en R (rapport de cotes)

J'essaie d'entreprendre une analyse de régression logistique en format R. J'ai suivi des cours sur ce matériel avec STATA. Je trouve très difficile de reproduire la fonctionnalité dans R. Est-il mature dans ce domaine? Il semble y avoir peu de documentation ou de conseils disponibles. La production du rapport de cotes semble nécessiter une installation epicalcet / ou epitools/ et d’autres, sans que je puisse me mettre au travail, sont périmées ou manquent de documentation. J'ai l'habitude glmde faire la régression logistique. Toute suggestion serait la bienvenue.

Je ferais mieux d'en faire une vraie question. Comment exécuter une régression logistique et générer des cotes de cotes dans R?

Voici ce que j'ai fait pour une analyse univariée:

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

Et pour les variables multiples:

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

Je l' ai ensuite regardé x, y, summary(x)et summary(y).

Est x$coefficientsde n'importe quelle valeur?

si vous voulez interpréter les effets estimés sous forme de rapports de cotes relatifs, faites simplement exp(coef(x))(vous donne , le changement multiplicatif du rapport de cotes pour si la covariable associée à augmente de 1). Pour les intervalles de probabilité de profil pour cette quantité, vous pouvez faire y = 1 $e^\beta$$y=1$$\beta$

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

EDIT: @caracal était plus rapide ...

Vous avez raison de dire que la sortie de R ne contient généralement que des informations essentielles et que d’autres doivent être calculées séparément.

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()vous donne les paramètres de régression estimés . Il est cependant plus facile d’interpréter (à l’exception de l’interception). e x p ( b j$b_{j}$$exp(b_{j})$

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

Pour obtenir le rapport de cotes, nous avons besoin du tableau croisé de classification du DV d'origine dichotomique et de la classification prédite en fonction d'un seuil de probabilité qu'il convient de choisir en premier. Vous pouvez également voir la fonction ClassLog()dans le package QuantPsyc(comme chl mentionné dans une question connexe ).

# 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

Pour le rapport de cotes, vous pouvez utiliser un package vcdou effectuer le calcul manuellement.

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

La page de statistiques de UCLA propose une régression logistique en R. Elle comprend une brève section sur le calcul des rapports de cotes.

Le paquet epiDisplay le fait très facilement.

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