How can I overlay linear and non-linear curve?

[gilbertlzrus]

Hi, fellow statisticians!

I am very new to R.

I am currently running a dose-response meta-analysis, following the tutorial provided by Shim et al

In short, the graph function uses the following command:

 with(predict(varname, ref, order=TRUE/FALSE, exp=TRUE/FALSE), 
     {plot() lines(lower 95% CI) lines(upper 95% CI)})

Based on the tutorial, this function was used for both linear and non-linear curves. However, I would like to combine both linear and non-linear curve, along with their 95% CIs. My question is:

1. What command should I use to combine these plots by using the with() command. Is there any manual/reference I can look up to?
2. I would also like to color the 95% CI area. I noticed that I can do that with Stata's rarea attribute in twoway function, is there any equivalent command in R?

Thank you very much

[MyKo101]

The best way to do something like this, (IMO) would be using ggplot2, rather than with() and plot(). The main gist of the difference (in a psuedo-code kinda way) is that with() would look like this:

with(predict(make_model(y~x,Data)),plot(...))

But ggplot() and tidyverse would look like this:

Data %>%
  make_model(y ~ x)
  predict() %>%
  ggplot() +
  geom_*(...)

The ggplot version (to me at least) looks a bit tidier. There are less surrounding brackets and the code flows easier. This is big thing amongst R users (base vs tidy), but I say use whichever makes your code easier for others to read.

So, you want to plot the main estimate from your data and a confidence band, I'd recommend using the geom_ribbon() function for the confidence interval and then when you get better at ggplot stuff, exploring the different themes available

I'll use the dosresmeta() function from Shim et al for making the model. I'll also say that this paper is definitely aimed at teaching the dose-response methods much more than teaching R because the way it teaches R is (quite frankly) appalling. If you're wanting to learn R from this, I'd suggest finding a much better resource.

library(tidyverse)
read.csv(“drma_bin.csv”, header=TRUE) %>%
  dosresmeta(logrr ~ dose, id=id, type=type, se=se, cases=cases, n=n, data=.) %>%
  predict() %>%
  ggplot() + #This starts a new canvas with your data as the main source
  geom_line(aes(dose,pred)) + #draw a line for your main estimate
 geom_ribbon(aes(ymin = ci.lb, ymax= ci.ub))# this draws your 95% CI

In base R (without the tidyverse stuff), this would best be written as:

data_bin <- read.csv(“drma_bin.csv”, header=TRUE)
model <- dosresmeta(logrr ~ dose, id=id, type=type, se=se, cases=cases, n=n,data = data_bin)
pred <- predict(model)
plot(pred$dose, pred$pred, type="l")
lines(pred$dose,pred$ci.lb)
lines(pred$dose,pred$ci.ub)

or using with() as suggested by Shim et al:

data_bin <- read.csv(“drma_bin.csv”, header=TRUE)
model <- dosresmeta(logrr ~ dose, id=id, type=type, se=se, cases=cases, n=n,data = data_bin)

with(
  predict(model),
  {
    plot(dose, pred);
    lines(dose, ci.lb);
    lines(dos, ci.ub)
  })

This is of course assuming you get the data out of the supp. material they provided which is just more effort than it should be. Learn R elsewhere, this paper is not useful for learning it.

Edit:
Also, if you ever need to know something about a function, almost all functions in R have their own help page. You can access it using ? so run ?with and it will pop up with a help file. This is much easier if using RStudio since it organises the windows in panes for you.

[gilbertlzrus]

Dear @MyKo101,

Thank you very much for the help! I managed to overlay and shade the 95% CIs using polygon(). Once again, thank you very much for the help!

Best regards,

[MyKo101]

I'm not 100% sure what you mean, however my interpretation is that you have groups within your data and you want separate lines and CIs for each group. Is this correct?

If so, you can use the group= argument within the aes() function to split your data by a variable, then use both col= and fill= for the same variable to colour them distinctly, col= assigns a border (for lines, this is the whole line, but for ribbons, it's just the outline) and fill= assigns the variable to fill your ribbon with.

Here's an example of them in action:

library(tidyverse)

df <- tibble() %>%
  expand(time = 1:10,grp=1:2) %>%
  mutate(Estimate = grp*time + rnorm(n()), #Time-trend dependent on the group
         Error = abs(rnorm(n()))+0.2,
         Lower_bound = Estimate - 1.96*Error, #Normal CI
         Upper_bound = Estimate + 1.96*Error,
         grp = c("A","B")[grp]) #Easier to understand output with discrete groups
         
ggplot(df,
       aes(x = time,      #aesthetics can be supplied to the
           y = Estimate,  #ggplot() function as the plain canvas
           group = grp,   #of our plot, and then specific ones
           col = grp,     #can be given laterwhen needed
           fill = grp)) +
  geom_line() +
  geom_ribbon(aes(ymin=Lower_bound, #Add the extra aesthetics here
                  ymax=Upper_bound),
              alpha=0.2) + #alpha makes it see-through
  theme_bw() #Add a nice theme to the plot