1
$\begingroup$

This is the desired graph

enter image description here

Here is the data: https://github.com/UnlimitedR/share/blob/main/mydata.csv

If so, could anyone please tell me how should I implement the R codes for it? Including drawing the CI area.

This is currently my code

df <- read.csv("mydata.csv")
print(sum(df$PRS))

library(ggplot2)
df$incidence <- 100 * df$PRS / sum(df$PRS,na.rm=T)
g <- ggplot(df, aes(x = CIT, y = incidence)) +
  geom_line()
print(g)

But the result is

enter image description here

btw, what does the incidence of PRS means? I think PRS is a binary variable?

I dont know if I have understood the meaning of incidence of PRS correctly?


Thanks for the help from Dave, I rewrote the code with ggplot

dat <- rio::import("https://raw.githubusercontent.com/UnlimitedR/share/main/mydata.csv")
library(ggplot2)
library(ggeffects)

mod <- glm(PRS ~ CIT, data=dat, family=binomial)
g <- ggpredict(mod, terms="CIT [all]") 

ggplot(g, aes(x, y = predicted, color="yellow")) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.5) +
  ggtitle("") +
  labs(x = "CIT", y = "Incidence of PRS, %") +
  scale_y_continuous(labels = ~sprintf("%.0f", .x*100)) +
  #scale_y_continuous(labels = percent_format())+
  theme_bw()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.2,0.9))+
  scale_color_discrete(name="",
                       labels="95% CI")

However, the color seems to have a little bit problems.

enter image description here

$\endgroup$
1
  • 3
    $\begingroup$ It looks like you're plotting the raw data, whereas the authors of the graphic you want to emulate are plotting a model fit. Do you know what kind of model they might be using? $\endgroup$ Dec 30, 2022 at 16:44

1 Answer 1

2
$\begingroup$

Maybe something like this:

dat <- rio::import("https://raw.githubusercontent.com/UnlimitedR/share/main/mydata.csv")
library(ggplot2)
library(ggeffects)

mod <- glm(PRS ~ CIT, data=dat, family=binomial)
g <- ggpredict(mod, terms="CIT [all]") 
plot(g) + 
  ggtitle("") + 
  labs(x="CIT", y="Incidence of PRS")

Created on 2022-12-30 by the reprex package (v2.0.1)

The code above assumes that the underlying model predicting PRS with CIT is a logistic regression.

To print without the % symbol, using the code in the comments:

plot(g) + 
  ggtitle("") + 
  labs(x="CIT", y="Incidence of PRS, %") + 
  scale_y_continuous(labels = ~sprintf("%.0f", .x*100)) + 
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.border = element_rect(fill = NA))
#> Scale for y is already present.
#> Adding another scale for y, which will replace the existing scale.


Edit: Modifying code in question

ggplot(g, aes(x, y = predicted)) +
  geom_line(aes(colour="Predictions\nw/95% CI")) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill="Predictions\nw/95% CI"), alpha = 0.5, colour="transparent") +
  ggtitle("") +
  labs(x = "CIT", y = "Incidence of PRS, %") +
  scale_y_continuous(labels = ~sprintf("%.0f", .x*100)) +
  #scale_y_continuous(labels = percent_format())+
  theme_bw()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = c(0.15,0.9)) +
  scale_fill_manual(values="gray50") + 
  scale_colour_manual(values="black") + 
  labs(colour="", fill = "")

$\endgroup$
6
  • $\begingroup$ Thank you! Im sorry Im new to R. I still have some minor problems about ploting. I checked it in RStudio and found that the g object is like a dataframe with three columns CIT, Predicted and 95% CI. Does that mean the plot function can automatically detect the predicted value and CI area and draw it? (like take the second column to draw the line and third column to draw the band?) $\endgroup$
    – user900476
    Dec 30, 2022 at 17:32
  • 1
    $\begingroup$ The first class for g is ggeffects, which means when you use plot() on the object g, it will first try to find plot.ggeffects() which is defined in the ggeffects package. That does know to use the predicted values and the confidence intervals to make the line and confidence envelope. To see the code for the function, you could type: ggeffects:::plot.ggeffects into R and press enter (without the parentheses after the function). $\endgroup$ Dec 30, 2022 at 17:38
  • $\begingroup$ Thank you! one last minor request is that, I want to make the graph as close as to the example. So I modified the plot statement a little bit ``` plot(g) + ggtitle("") + labs(x="CIT", y="Incidence of PRS, %") + scale_y_continuous(labels = percent_format())+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ theme(panel.border = element_rect(fill = NA)) ``` But it still shows the percentage symbol, how should I modify it to not show it? $\endgroup$
    – user900476
    Dec 30, 2022 at 17:43
  • $\begingroup$ @user900476 I updated the answer to remove the percent sign. $\endgroup$ Dec 30, 2022 at 17:48
  • $\begingroup$ Hello, I rewrote the code with ggplot and tried to add a legend like in the example, would you mind taking a look at it? $\endgroup$
    – user900476
    Dec 30, 2022 at 19:23

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.