Search code examples
rdplyrlinear-regressionsegmentpiecewise

How to use segmented package when working with data frames with dplyr package to perform piecewise linear regression?


My data frame is separated by groups. I want to perform piecewise linear regression on each group and for that I intend to use the segmented package.

First I created the linear models for each group using the dplyr package. The next step is to segment these models, however this is where I'm stuck. Any tips or other way to do this? The ultimate goal is to use these segments to make a graph.

library(dplyr)
library(segmented)

Group <- c("A", "B")
x <- 0:10
y <- c(0, 0.4, 0.6, 0.8, 0.9, 0.9, 0.95, 0.97, 0.98, 0.99, 1,
       0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1)

df <- expand.grid(x = x,
                  Group = Group)

df$y <- y

Segment <- df %>%
  group_by(Group) %>%
  do(my.lm = lm(data = .,
                formula = y ~ x))

Unsuccessful attempts:

Segment <- df %>%
  group_by(Group) %>%
  do(my.lm = lm(data = .,
                formula = y ~ x),
     my.seg = segmented(my.lm,
                        seg.Z = x))


Segment <- df %>%
  group_by(Group) %>%
  do(my.lm = lm(data = .,
                formula = y ~ x)) %>%
  do(my.seg = segmented(my.lm,
                        seg.Z = x))


Segment <- df %>%
  group_by(Group) %>%
  do(my.lm = lm(data = .,
                formula = y ~ x)) %>%
  mutate(my.seg = segmented(my.lm,
                        seg.Z = x))

Solution

  • One option is to wrap with tryCatch and return a NA for possible errors

    library(dplyr)
    out <- df %>% 
        nest_by(Group) %>%
        mutate(my.lm = list(lm(y ~ x, data = data)),
            my.seg = list(tryCatch(segmented(my.lm, seg.Z = ~ x),
             error = function(e) list(NA))))
    

    -output

    > out
    # A tibble: 2 x 4
    # Rowwise:  Group
      Group               data my.lm  my.seg    
      <fct> <list<tibble[,2]>> <list> <list>    
    1 A               [11 × 2] <lm>   <segmentd>
    2 B               [11 × 2] <lm>   <list [1]>
    > out$my.seg
    [[1]]
    Call: segmented.lm(obj = my.lm, seg.Z = ~x)
    
    Meaningful coefficients of the linear terms:
    (Intercept)            x         U1.x  
        0.03333      0.30000     -0.27488  
    
    Estimated Break-Point(s):
    psi1.x  
     2.691  
    
    [[2]]
    [[2]][[1]]
    [1] NA