Search code examples
rtidyversenls

How to return NA when nls model is unable to fit data?


I am applying groupwise nlsLM model and solving it to get the x value when the y is highest individually for each group. But some of the studies are giving following error

Error in mutate(): ℹ In argument: model = map(data, ~nlsLM(y ~ SSlinp(x, a, b, xs), data = .x)). Caused by error in map(): ℹ In index: 2. Caused by error in nlsModel(): ! singular gradient matrix at initial parameter estimates

How to avoid those studies and run for those studies where nlsLM model can be run?

Here is the code I am using with sample data

library(tidyverse)
library(gslnls)
library(broom)
library(minpack.lm)
library(nlraa)

df %>%
  nest(data = -SN) %>%
  mutate(model = map(data, ~ nlsLM(y ~ SSlinp(x, a, b, xs), data = .x))) %>% 
  mutate(x_opt = map(model, tidy)) %>%
  unnest(x_opt) %>% 
  filter(term %in% "xs") %>% 
  mutate(y_opt = unlist(map2(estimate, model, ~predict(.y, list(x = .x))))) %>%
  select(SN, estimate, y_opt)

Data

df = structure(list(SN = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L), x = c(0L, 60L, 90L, 120L, 150L, 0L, 60L, 
90L, 120L, 150L, 0L, 60L, 90L, 120L, 150L), y = c(3569.5, 6698.1, 
6595.2, 6834.4, 6966.7, 3649.7, 4895.1, 5162, 4443.5, 5038.9, 
4097.7, 5664.9, 5518.1, 7608.1, 8081.6)), class = "data.frame", row.names = c(NA, 
-15L))

Solution

  • This is a clunky, partial solution. You'd have to add a similar my_predict function that will return the right null/NA results if the model is bad.

    Alternatively, you could break the pipeline after the model-fitting stage and get rid of the models that didn't work (e.g. filter(!inherits(model, "try-error")) or some such), then run tidy/predict without hassles on the remaining elements ...

    empty <- tibble(
        term = c("a", "b", "xs"),
        estimate = NA_real_,
        std.error = NA_real_,
        p.value = NA_real_
    )
    my_tidy <- function(m) {
        if (inherits(m, "try-error")) empty else broom::tidy(m)
    }
    
    df %>%
      nest(data = -SN) %>%
        mutate(model = map(data, ~ try(nlsLM(y ~ SSlinp(x, a, b, xs), data = .x)))) %>% 
        mutate(x_opt = map(model, my_tidy)) %>%
        unnest(x_opt) %>%
        filter(term %in% "xs")