2
$\begingroup$

I would like to use the glm() function in R to run a multivariate logistic regression. I have also run bi-variate statistics for each variable but want a test that controls for all variables at once (what I understand the glm function in R does - let me know if I should use a different function). My problem is that I don't know what to report in an academic paper and can't find anything online. There are clear guidelines for reporting other tests (example: t (34.17) = 2.22, p = 0.033) but I don't know what to report for the glm. Here is the output of my glm:

Call:
glm(formula = EPI ~ premed + kg + sim + bup + premed.min + ami.type + breed, family = "binomial", data = df)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7316  -0.2980  -0.2351  -0.1573   2.9108  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)        -1.364678   1.239446  -1.101   0.2709  
premedT             1.103479   0.499234   2.210   0.0271 *
kg                 -0.032467   0.014633  -2.219   0.0265 *
simS                0.193276   0.431850   0.448   0.6545  
bupN                0.037000   1.048603   0.035   0.9719  
premed.min          0.002211   0.002542   0.870   0.3845  
ami.typeCRI       -12.604613 781.507743  -0.016   0.9871  
ami.typeDRIP       -0.226335   0.466006  -0.486   0.6272  
breedYorkshireMix  -0.852239   0.489624  -1.741   0.0818 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 272.44  on 938  degrees of freedom
Residual deviance: 255.21  on 930  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 273.21

Number of Fisher Scoring iterations: 15

It doesn't help that I cant figure out how I would report a bi-variate logistic regression either (I've seen people report odds ratios, χ2 values etc). I have found out how to get confidence intervals and odds ratios in R if I need those:

> exp(cbind(OR = coef(LM), confint(LM)))
Waiting for profiling to be done...
                            OR         2.5 %       97.5 %
(Intercept)       2.554629e-01  2.222514e-02 2.873313e+00
premedT           3.014635e+00  1.198440e+00 8.735251e+00
kg                9.680544e-01  9.402870e-01 9.957784e-01
simS              1.213217e+00  5.144926e-01 2.832683e+00
bupN              1.037693e+00  5.650816e-02 5.366901e+00
premed.min        1.002213e+00  9.966528e-01 1.006788e+00
ami.typeCRI       3.356495e-06 9.742338e-102 1.000098e+12
ami.typeDRIP      7.974511e-01  3.090303e-01 1.957706e+00
breedYorkshireMix 4.264589e-01  1.724568e-01 1.213685e+00
There were 36 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: glm.fit: fitted probabilities numerically 0 or 1 occurred

I want to say something like "The results of the multivariate logistic regression indicated that subjects given pre-medication "T" were more likely to have outcome "EPI" than subjects given pre-medication "X" (blah, blah, p = 0.027) and subjects that had lower "kg" were also more likely to have outcome "EPI" than subjects that had higher "kg" (blah, blah, p = 0.0265)."

EDIT: EPI is a binary variable that is assigned 0 or 1

$\endgroup$
2
$\begingroup$

Is your EPI variable a binary variable taking the values 0 or 1? Then the glm() function the way you used it here will fit a binary logistic regression model relating this binary variable to the predictors of interest.

Before you report the results from this model, note that R posts a concerning warning message that fitted probabilities numerically 0 or 1 have occurred. If you examine the standard errors for your estimated regression coefficients, you will note that the standard error for ami.typeCRI is huge compared to the other standard errors. Same for the confidence interval limits! This could be a sign that you have very few observations in the category CRI of your ami.type predictor variable. You may have no choice but to combine this category with one of the other categories: DRIP or the reference category used in your model. After combining, you can refit the model with the new version of ami.type to see if R will stop posting the pesky error message.

Once you get a glm model with no more warnings, you can check the model diagnostics (e.g., using the DHARMa package). Presuming the model diagnostics look fine, you can report the odds ratios produced by the model and their corresponding 95% confidence intervals:

round(exp(cbind(OR = coef(LM), confint(LM))),2)

Then you can say things like:

The results of the multiple binary logistic regression indicated that, all else being equal, subjects given pre-medication "T" had higher odds of having the outcome "EPI" than subjects given pre-medication "X" (OR = 1.92 ; 95% CI: 1.15 to 2.45; p = 0.027). Subjects who had a larger body weight (kg) had lower odds of having the outcome "EPI" than subjects with a lower body weight (OR = 0.59, 95% CI: 0.26 to 0.78, p = 0.0265).

I made up the numbers just to give you an idea of how you would report your findings in terms of odds ratios, thereby using language involving odds.

The medical literature is full of papers which report this type of model.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Thanks for your very clear answer. Yes EPI is a binary variable that is assigned 0 or 1. I'll look into the DARMAa package. I think I can justify collapsing the ami.type predictor, so I'll do that. Thanks again. $\endgroup$ – B.Kenobi Jul 24 at 14:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.