I'm using or_glm() to calculate odds ratios, using this reproducible example:
library(oddsratio)
or_glm(data = data_glm,
model = glm(admit ~ gre gpa rank,
data = data_glm,
family = "binomial"),
incr = list(gre = 1, gpa = 1, rank = 1))
I have two questions:
- How can I also extract a p-value for each odds ratio?
- How can I get an odds ratio for "gre" adjusted for for "gpa" and "rank"?
CodePudding user response:
I would try as follows:
library(oddsratio)
library(mfx)
model = glm(admit ~ gre gpa rank,
data = data_glm,
family = "binomial")
logitor(admit ~ gre gpa rank,data=data_glm)
Call:
logitor(formula = admit ~ gre gpa rank, data = data_glm)
Odds Ratio:
OddsRatio Std. Err. z P>|z|
gre 1.0022670 0.0010965 2.0699 0.0384651 *
gpa 2.2345448 0.7414651 2.4231 0.0153879 *
rank2 0.5089310 0.1610714 -2.1342 0.0328288 *
rank3 0.2617923 0.0903986 -3.8812 0.0001039 ***
rank4 0.2119375 0.0885542 -3.7131 0.0002047 ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef(model))
(Intercept) gre gpa rank2
0.0185001 1.0022670 2.2345448 0.5089310
rank3 rank4
0.2617923 0.2119375
exp(cbind(OR=coef(model), confint(model)))
Waiting for profiling to be done...
OR 2.5 % 97.5 %
(Intercept) 0.0185001 0.001889165 0.1665354
gre 1.0022670 1.000137602 1.0044457
gpa 2.2345448 1.173858216 4.3238349
rank2 0.5089310 0.272289674 0.9448343
rank3 0.2617923 0.131641717 0.5115181
rank4 0.2119375 0.090715546 0.4706961
CodePudding user response:
The coefficients of a logistic regression (i.e. a binomial glm) are simply the natural log of the odds ratios, so you get the odds ratio by doing exp(coefficient)
. Similarly, the 95% confidence intervals are the exponents of the coefficients /- 1.96 * the standard error. The p values for the odds ratio are already given in the coefficient table, and these don't need modified at all.
It is therefore easy to make the table yourself. Let's start by getting the coefficient table from the model:
library(oddsratio)
mod <- glm(admit ~ gre gpa rank,
data = data_glm,
family = "binomial")
tab <- summary(mod)$coef
tab
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) -3.989979073 1.139950936 -3.500132 0.0004650273
#> gre 0.002264426 0.001093998 2.069864 0.0384651284
#> gpa 0.804037549 0.331819298 2.423119 0.0153878974
#> rank2 -0.675442928 0.316489661 -2.134171 0.0328288188
#> rank3 -1.340203916 0.345306418 -3.881202 0.0001039415
#> rank4 -1.551463677 0.417831633 -3.713131 0.0002047107
We can make the odds ratio table like this:
data.frame(variable = rownames(tab),
oddsratio = round(exp(tab[,1]), 3),
ci_low = round(exp(tab[,1] - 1.96 * tab[,2]), 3),
ci_high = round(exp(tab[,1] 1.96 * tab[,2]), 3),
pval = scales::pvalue(tab[,4]),
row.names = NULL)[-1,]
#> variable oddsratio ci_low ci_high pval
#> 2 gre 1.002 1.000 1.004 0.038
#> 3 gpa 2.235 1.166 4.282 0.015
#> 4 rank2 0.509 0.274 0.946 0.033
#> 5 rank3 0.262 0.133 0.515 <0.001
#> 6 rank4 0.212 0.093 0.481 <0.001
Compare this to the output of your or_glm
function:
#> predictor oddsratio ci_low (2.5) ci_high (97.5) increment
#> 1 gre 1.002 1.000 1.004 1
#> 2 gpa 2.235 1.174 4.324 1
#> 3 rank2 0.509 0.272 0.945 Indicator variable
#> 4 rank3 0.262 0.132 0.512 Indicator variable
#> 5 rank4 0.212 0.091 0.471 Indicator variable
Created on 2022-06-05 by the reprex package (v2.0.1)