Search code examples
rggplot2nls

How can I ggplot a tibble with models and run F-test


I have a tibble where I have fitted a model to different subsets of data. I now want to plot each model over the data which contains all points before I do an F-test to see if there is any gain to the model from including the "Site_class" variable.

Data:

sitedata <- structure(list(Site_class = c("1", "1", "1", "1", "1", "1", "1", 
"1", "1", "1", "1", "1", "4", "4", "4", "4", "4", "4", "4", "4", 
"4", "4", "4", "4", "All", "All", "All", "All", "All", "All", 
"All", "All", "All", "All", "All", "All", "All", "All", "All", 
"All", "All", "All", "All", "All", "All", "All", "All", "All", 
"All", "All", "All", "All", "All", "All", "All", "All", "All", 
"All", "All", "All", "All", "All", "All", "All", "All", "All", 
"All", "All", "All", "All", "All", "All", "All", "All", "All", 
"All", "All", "All", "All", "All", "All", "All", "All"), VarA = c(2.4, 
6.5, 11.3, 16.1, 20.5, 24.3, 27.4, 30, 32, 33.6, 34.9, 36, 0.75, 
2.45, 4.75, 7.45, 10.3, 13.05, 15.55, 17.7, 19.4, 20.75, 21.9, 
22.8, 2.4, 6.5, 11.3, 16.1, 20.5, 24.3, 27.4, 30, 32, 33.6, 34.9, 
36, 1.85, 5.15, 9.1, 13.2, 17.1, 20.55, 23.45, 25.9, 27.8, 29.3, 
30.55, 31.6, 1.3, 3.8, 6.95, 10.35, 13.7, 16.8, 19.5, 21.8, 23.6, 
25.05, 26.25, 27.2, 0.75, 2.45, 4.75, 7.45, 10.3, 13.05, 15.55, 
17.7, 19.4, 20.75, 21.9, 22.8, 1.1, 2.6, 4.6, 6.9, 9.3, 11.6, 
13.6, 15.2, 16.5, 17.55, 18.4), VarB = c(12, 81, 220, 403, 605, 
806, 991, 1153, 1288, 1399, 1495, 1578, 1, 11, 45, 106, 189, 
283, 381, 473, 552, 619, 675, 723, 12, 81, 220, 403, 605, 806, 
991, 1153, 1288, 1399, 1495, 1578, 4, 51, 148, 286, 446, 609, 
760, 893, 1008, 1105, 1186, 1255, 2, 27, 93, 190, 307, 433, 557, 
670, 766, 848, 917, 975, 1, 11, 45, 106, 189, 283, 381, 473, 
552, 619, 675, 723, 2, 11, 42, 94, 161, 234, 304, 368, 423, 469, 
508)), row.names = c(NA, -83L), class = c("grouped_df", "tbl_df", 
"tbl", "data.frame"), groups = structure(list(Site_class = c("1", 
"4", "All"), .rows = list(1:12, 13:24, 25:83)), row.names = c(NA, 
-3L), class = c("tbl_df", "tbl", "data.frame"), .drop = FALSE))

This is how I've worked so far.

library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)

#Create model
modfit <- function(df){
  nls(VarB ~ a * (VarA^b), data=df, start=c(a=1,b=1))
}

#Nest the original data.frame
sitedata <- sitedata %>% group_by(Site_class) %>% nest()

#Fit the models
sitedata_model <- sitedata %>% mutate(
  Model= map(.x=data, .f= modfit)
)

#Attempt to plot the models:
ggplot(sitedata_model[[2]][[3]], aes(x=VarA,y=VarB)) + #All the data
  geom_point()+
  geom_function(fun=~sitedata_model[[3]][[1]]) + # I assume I will have to plot them separately?
  geom_function(fun=~sitedata_model[[3]][[2]]) +
  geom_function(fun=~sitedata_model[[3]][[3]]) 

It seems I have succeeded creating the models, but it then doesn't understand my call to plot them. I have also tried using predict, without success.

How can I:

  1. plot them on my graph &
  2. run an F-test between all 3 models to see if there is any difference?

Solution

  • geom_function expects a function. Not a formula, not a model, but a function.

    I haven't had much experience with this geom, & I'm not sure if the following is efficient, but it does seem to work for your use case:

    sitedata_model <- sitedata %>% 
      #Fit the models
      mutate(Model = purrr::map(.x=data, .f= modfit)) %>%
      
      # extract formula from each model, convert to one-sided form, &
      # replace coefficients with fitted values, & store in dataframe
      # as character string 
      rowwise() %>%
      mutate(func = formula(Model) %>% 
               as.character() %>% 
               magrittr::extract(3) %>%
               gsub("VarA", ".x", ., fixed = T) %>%
               gsub("a", Model$m$getPars()[1], .) %>%
               gsub("b", Model$m$getPars()[2], .) %>%
               paste("~", ., collapse = "")) %>%
      ungroup()
    
    # plot
    ggplot(data = sitedata_model$data[[3]]) +
      geom_point(aes(x = VarA, y = VarB)) + 
      
      # add formula in each row as a separate geom_function layer
      lapply(seq(1, nrow(sitedata_model)),
             function(i) geom_function(fun = rlang::as_function(formula(sitedata_model$func[i])),
                                       aes(colour = sitedata_model$Site_class[i]))) +
      
      # change legend name (can also change palette / labels / etc.)
      scale_colour_discrete(name = "Site class")
    

    plot

    As for running F-test, are you referring to ANOVA?

    > anova(sitedata_model$Model[[1]], sitedata_model$Model[[2]], sitedata_model$Model[[3]])
    Analysis of Variance Table
    
    Model 1: VarB ~ a * (VarA^b)
    Model 2: VarB ~ a * (VarA^b)
    Model 3: VarB ~ a * (VarA^b)
      Res.Df Res.Sum Sq  Df  Sum Sq F value Pr(>F)
    1     10      60.14                           
    2     10      90.36   0    0.00               
    3     57     741.38 -47 -651.02  1.5328 0.2386
    

    See here for an explanation of how to interpret ANOVA output.