Home > other >  Adjusted odds ratios using the or_glm() function?
Adjusted odds ratios using the or_glm() function?

Time:06-06

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:

  1. How can I also extract a p-value for each odds ratio?
  2. 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)

  • Related