0

I've data that look like this.

Probability <- c(1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0)
Score2 <- c(2,3,4,2,3,4,2,3,4,4,3,2,3,2,3,2,3,2,3,4,2,3,4,2,3,3,4,3,2,3,2,3,2,3)
Data2 <- data.frame(Probability, Score2)

Data
   Probability Score2
1            1      2
2            1      3
3            1      4
4            1      2
5            1      3
6            1      4
7            1      2
8            1      3
9            1      4
10           0      4
11           0      3
12           0      2
13           0      3
14           0      2
15           0      3
16           0      2
17           0      3
18           1      2
19           1      3
20           1      4
21           1      2
22           1      3
23           1      4
24           1      2
25           1      3
26           1      3
27           0      4
28           0      3
29           0      2
30           0      3
31           0      2
32           0      3
33           0      2
34           0      3

I need to plot not just a correlation plot, but also add odds ratios comparing the odds of Y at 1 standard deviation above and below the mean of X (with the corresponding p-value).

The following gets me everything but the odds ratio of +/- 1 SD and its p-value.

ggplot(Data2, aes(Score2, Probability))+ 
  geom_smooth(method='lm', alpha = .3, color = "black")+
  stat_cor(method = "pearson", label.x.npc = "left", label.y.npc= "top", label.sep = "
")+
  scale_colour_grey(start = .6, end = .2)+
  scale_fill_grey(start = 0.6, end = 1)+
  theme_classic()+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1), limits = c(0, 1))

image of correlation plot without odds ratio

Question

How can I add the odds ratio (comparing Probability at 1 SD above and below the mean of X)?

2 Answers 2

1

You could estimate the CI using contrast package for the linear probability model also.

library(ggpubr)

mean_x = mean(Score2)
sd_x = sd(Score2)

Linear Probability Model

lm_prob <- lm(Probability ~ Score2,data=Data2)

pred_probs = predict(lm_prob,newdata = data.frame(Score2 = c(mean_x - sd_x, mean_x+sd_x)))
or_pred = (pred_probs[2]/(1-pred_probs[2]))/(pred_probs[1]/(1-pred_probs[1]))

Logistic Regression

glm_prob <- glm(Probability ~ Score2,data=Data2,family=binomial())

glm_pred_probs = predict(glm_prob,newdata = data.frame(Score2 = c(mean_x - sd_x, mean_x+sd_x)),type = "response",se.fit = TRUE)
glm_or_pred = (glm_pred_probs$fit[2]/(1-glm_pred_probs$fit[2]))/(glm_pred_probs$fit[1]/(1-glm_pred_probs$fit[1]))

The contrasts statements can be updated with mean_x + sd_x and mean_x - sd_x. It doesn't change the result.

library(contrast)
glm_contrast <- 
  contrast(glm_prob, 
           list(Score2 = sd_x),
           list(Score2 = -sd_x)
  )
print(glm_contrast, X = TRUE)

or_ci = paste0(round(exp(glm_contrast$Contrast),2), 
               ", 95% CI:",
               round(exp(glm_contrast$Lower),2),
               ", ",
               round(exp(glm_contrast$Upper),2),
               ", p = ",
               round(glm_contrast$Pvalue,3)
               )

Plot

ggplot(Data2, aes(Score2, Probability))+ 
  geom_smooth(method='lm', alpha = .3, color = "black")+
  stat_cor(method = "pearson", label.x.npc = "left", label.y.npc= "top", label.sep = "
")+
  annotate("text", x=3.0, y=0.9, label= paste0("OR = ",or_ci)) +
  scale_colour_grey(start = .6, end = .2)+
  scale_fill_grey(start = 0.6, end = 1)+
  theme_classic()+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1), limits = c(0, 1))
0

@jvargh7 got me most of the way to this solution (minus the standard deviation bracket) here

library(ggpubr)

### @jVargh7's code ###    
library(ggplot2)  
library(ggpubr)
  
mean_x = mean(Score2)
sd_x = sd(Score2)

lm_prob <- lm(Probability ~ Score2,data=Data2)

pred_probs = predict(lm_prob,newdata = data.frame(Score2 = c(mean_x - sd_x, mean_x+sd_x)))
or_pred = (pred_probs[2]/(1-pred_probs[2]))/(pred_probs[1]/(1-pred_probs[1]))

glm_prob <- glm(Probability ~ Score2,data=Data2,family=binomial())

glm_pred_probs = predict(glm_prob,newdata = data.frame(Score2 = c(mean_x - sd_x, mean_x+sd_x)),type = "response",se.fit = TRUE)
glm_or_pred = (glm_pred_probs$fit[2]/(1-glm_pred_probs$fit[2]))/(glm_pred_probs$fit[1]/(1-glm_pred_probs$fit[1]))

install.packages("contrast")
library(contrast)
glm_contrast <- contrast(glm_prob, list(Score2 = sd_x), list(Score2 = -sd_x))
print(glm_contrast, X = TRUE)

or_ci = paste0(round(exp(glm_contrast$Contrast),2), 
               ", 95% CI: ",
               round(exp(glm_contrast$Lower),2),
               ", ",
               round(exp(glm_contrast$Upper),2),
               ", p = ",
               round(glm_contrast$Pvalue,3)
)   

ggplot(Data2, aes(Score2, Probability))+ 
  geom_smooth(method='lm', alpha = .3, color = "black")+
  stat_cor(method = "pearson", label.x = 2.75, label.y= .6, label.sep = ",")+
  scale_colour_grey(start = .6, end = .2)+
  scale_fill_grey(start = 0.6, end = 1)+
  theme_classic()+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1), limits = c(0, 1))+ 
### new code ###
  geom_bracket(xmin = mean(Score2)-sd(Score2), 
               xmax = mean(Score2)+sd(Score2), 
               y.position = .9, 
               label = paste0("OR = ", or_ci), 
               tip.length = c(0.42, 0.29),
               vjust = -1)

Example output below

image of linear regression plot with odds ratio of Y at 1 S.D. above and below the mean of X

Again, thanks to @jvargh7 for getting me most of the way to this solution here.

1
  • Updated the solution for contrast estimation.
    – jsv
    Commented Mar 4, 2021 at 21:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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