Search code examples
rggplot2model

Plot models in ggplot2

I want to replot a graph with ggplot instead of the plot function but have difficulty getting the same values.

The package growthcurver produces models I would like to plot with ggplot.

> data
   time         T1
1     0 0.01666667
2     1 0.06000000
3     2 0.12000000
4     3 0.34000000
5     4 0.29666667
6     5 0.34000000
7     6 1.23666667
8     7 2.21333333
9     8 2.64333333
10   24 3.89000000

model<- growthcurver:: SummarizeGrowth(data$time, data$T1)


plot(model.wt)

enter image description here

tibble(t = data$time, pred = model$model$m$fitted(),   N = data$T1) |> 
  ggplot(aes(t, pred)) +
  geom_point(aes(y = N)) +
  geom_line()

enter image description here

Solution

  • I had responded to your earlier question and with your current data I can see it doesn't quite work - especially if you have gaps. Borrowing some of the background code to the plotting function within the package you can create a new predict method for this class and use that:

    library(growthcurver)
library(tidyverse)

data <- tribble(
  ~time,  ~T1,
  0,      0.01666667,
  1,      0.06000000,
  2,      0.12000000,
  3,      0.34000000,
  4,      0.29666667,
  5,      0.34000000,
  6,      1.23666667,
  7,      2.21333333,
  8,      2.64333333,
  24,     3.89000000
  )

gc_mod <- SummarizeGrowth(data$time, data$T1)

    The NAtT function from the package produces the full set of fitted values when given a new set of independent variables:

    tibble(x = seq(min(data$time), max(data$time), length = 30),
       y = NAtT(gc_mod$vals$k, gc_mod$vals$n0, gc_mod$vals$r, x)) |> 
  ggplot(aes(x, y)) +
  geom_line()

    We can use this to make a new predict function tailored for the gcfit object type:

    
predict.gcfit <- function(gc_mod, new_data = NULL) {
  
  if(is.null(new_data)) {
    new_data <- gc_mod$data$t
  } 
  
  pred <- NAtT(gc_mod$vals$k, gc_mod$vals$n0, gc_mod$vals$r, new_data)
  
  pred
}

tibble(time = seq(0, 24),
       pred = predict(gc_mod, new_data = time)) |> 
  left_join(data, by = "time") |> 
  ggplot(aes(time)) +
  geom_point(aes(y = T1)) +
  geom_line(aes(y = pred), colour = "red")
#> Warning: Removed 15 rows containing missing values (`geom_point()`).