Can I make my for loop in a for loop more efficient?

The slider package specializes in helping with these kinds of rolling operations. When used with a data frame, the slide() functions apply a function for each window of rows in the data. You could use that to power the model fitting process.

First read the example data:

library(tidyverse)

df <- read.table(header = TRUE, text = "
Time   A   B   C
 330 102 179 303
 340 103 194 308
 350 101 198 348
 360 114 199 347
 370 120 214 371
 380 131 224 420
 390 128 226 430
 400 128 246 481
 410 138 260 541
 420 146 277 583
 430 155 290 640
 440 154 315 653")

We can then replace the outer loop over outcomes with map_dfr() which combines the results in a data frame, and the inner loop with slide_dfr() which does the same, but for a window of rows:

window_width <- 10
outcomes <- c("A", "B", "C")

# Estimate coefficients for each outcome
coefs <- map_dfr(set_names(outcomes), function(outcome) {
  # Fit model for each complete window in `df`
  slider::slide_dfr(df, function(data) {
    # Build model formula for outcome and fit
    f <- reformulate("Time", outcome)
    m <- lm(formula = f, data = data)
    set_names(coef(m), c("Inter", "Slope"))
  }, .after = window_width - 1, .complete = TRUE, .names_to = "Window")
}, .id = "Outcome")

coefs
#>   Outcome Window       Inter     Slope
#> 1       A      1   -67.30909 0.5024242
#> 2       A      2   -89.20000 0.5600000
#> 3       A      3   -87.54545 0.5545455
#> 4       B      1  -158.98182 1.0151515
#> 5       B      2  -191.86667 1.1030303
#> 6       B      3  -265.72727 1.2927273
#> 7       C      1  -749.07273 3.0993939
#> 8       C      2  -939.80000 3.6018182
#> 9       C      3 -1019.60000 3.8000000

Note that we kept the coefficients from all windows above, not just the top 5. It’s best to save the summarization of results to be done outside the loops, if you have enough memory for it. That way, if you decide that you wanted a different summary after all, you won’t have to redo the expensive part of fitting the models.

With the coefficients in a neat data frame, it’s straightforward to summarize:

n_top_slopes <- 5

mean_top_coefs <- coefs %>% 
  group_by(Outcome) %>% 
  slice_max(Slope, n = n_top_slopes) %>%
  summarise(across(c(Inter, Slope), mean))
#> `summarise()` ungrouping output (override with `.groups` argument)

mean_top_coefs
#> # A tibble: 3 x 3
#>   Outcome  Inter Slope
#>   <chr>    <dbl> <dbl>
#> 1 A        -81.4 0.539
#> 2 B       -206.  1.14 
#> 3 C       -903.  3.50

A quick way to plot these linear fits with the data would be with geom_abline():

df_long <- df %>% 
  pivot_longer(
    cols = all_of(outcomes),
    names_to = "Outcome",
    values_to = "Value"
  )

scatterplot <- function() {
  ggplot(df_long, aes(Time, Value)) +
    facet_wrap(~ Outcome) + geom_point()
}

scatterplot() +
  geom_abline(data = mean_top_coefs, aes(intercept = Inter, slope = Slope))

A slightly more involved but more general approach would be to consturct a new data frame of predictions first, and plot with geom_line():

mean_top_fits <- df_long %>% 
  inner_join(mean_top_coefs, by = "Outcome") %>% 
  mutate(Value = Inter + Slope * Time)

scatterplot() + geom_line(data = mean_top_fits)

Using this method, we could e.g. plot the fits for each individual window:

time_range <- function(window) {
  range(df$Time[seq(window, window + window_width - 1)])
}

# Calculate predicted lines for each window
window_fits <- coefs %>% 
  mutate(Time = map(Window, time_range)) %>% 
  unnest_longer(Time) %>% 
  mutate(Value = Inter + Slope * Time)

scatterplot() + geom_line(data = window_fits, aes(group = Window))

I'm not sure the code you posted works exactly as you want: you are using the same index for your nested for loops (i in both cases), so I'd expect the internal loop may rewrite the external one [edit: actually not: see comments]. Also, when I try to execute this code, the definition of dat fails: when i is equal to the number of columns, then i+1 is outside of the data frame. And finally the code doesn't correspond exactly to your description, e.g. taking 10 elements or 15, etc.

I tried to write a pure dplyr version based on your description, I'm not sure I understood everything though. It should be faster than the version with the for loop, but it's hard to be sure since we're currently not doing the same thing.

Let's first process a single column. To slice the column into groups of 10 observations, I define a function:

library(tidyverse)

make_groups <- function(n, k){
  # return a factor of length n defining groups of size k, except the last one
  nb_groups <- floor(n/k)
  overhang <- n - k*nb_groups
  c(rep(1:nb_groups, each = k), rep(nb_groups+1, overhang))
}

And some parameters

# how many rows to group together
group_size <- 10
# how many top coefs to use for mean
mean_top <- 2

Now, I'll have several steps: first, add a column to group the rows 10 by 10. Then run the lm(). Since the result of lm() is a complex object, I store it as a list-column. Then, I can simply do a mutate() to extract the coefficients of interest, and a summarize() to group the observations and only keep the mean of the top ones.

df %>%
  mutate(group = make_groups(n(), k=group_size)) %>%
  group_by(group) %>%
  summarize(mod = list(lm(A~Time))) %>%
  mutate(intercept = map_dbl(mod, ~.x$coefficients[[1]]),
         slope = map_dbl(mod, ~ .x$coefficients[[2]])) %>%
  summarize(mean_intercept = mean(sort(intercept, decreasing = TRUE)[1:mean_top]),
            mean_slope = mean(sort(slope,decreasing = TRUE)[1:mean_top]))
# A tibble: 1 x 2
#   mean_intercept mean_slope
#            <dbl>      <dbl>
# 1           65.3      0.201

This seems to work, so now how do we run it for every column? One possibility would be to use across(), but it gets complicated and unreadable, so I'll go with a simpler loop approach. Rather than a for loop, I use the equivalent map_df(), which does the looping and formatting by itself:

map_df(colnames(df)[2:4],
       function(cur_col_name){
         df %>%
           mutate(group = make_groups(n(), k=group_size)) %>%
           group_by(group) %>%
           summarize(mod = list(lm(formula(paste0(cur_col_name,"~Time")))),
                     .groups = "drop") %>%
           mutate(intercept = map_dbl(mod, ~.x$coefficients[[1]]),
                  slope = map_dbl(mod, ~ .x$coefficients[[2]])) %>%
           summarize(mean_intercept = mean(sort(intercept, decreasing = TRUE)[1:mean_top]),
                     mean_slope = mean(sort(slope,decreasing = TRUE)[1:mean_top]),
                     .groups = "drop")
       }
)
# A tibble: 3 x 2
#   mean_intercept mean_slope
#            <dbl>      <dbl>
# 1           65.3      0.201
# 2         -472.       1.76 
# 3         -334.       2.20 

So, finally we can compare performance (but since the two functions are not doing the same thing, that's not very meaningful). After resetting everything (ctrl+shift+F10 on RStudio) I run this:

library(tidyverse)
ex_df <- read.table(text="Time    A   B   C
330 102 179 303
340 103 194 308
350 101 198 348
360 114 199 347
370 120 214 371
380 131 224 420
390 128 226 430
400 128 246 481
410 138 260 541
420 146 277 583
430 155 290 640
440 154 315 653",header=TRUE)




with_for_loop <- function(df){
  regressors<-setdiff(names(df),"Time")
  
  coef2 <- data.frame()
  coef3 <- data.frame()
  
  #dat <- data.frame(df[[1]],df[[n]])
  
  for (i in 1:length(regressors)) {
    dat <- data.frame(df[[1]], df[[(i + 1)]])
    for (i in 1:length(dat[[1]])) {
      temp_data <- dat[i:(i + 15),]
      linm <- lm(temp_data[[2]] ~ temp_data[[1]], data = temp_data)
      Inter <- summary(linm)$coefficients[1]
      Slope <- summary(linm)$coefficients[2]
      coef2 <-
        data.table::rbindlist(list(coef2, data.frame(Inter, Slope)), use.names = T)
    }
    mean <- coef2 %>%
      arrange(desc(Slope)) %>%
      slice(2:7)
    meanx <-  sapply(mean, FUN = mean)
    meanx <-
      data.frame(lapply(meanx, type.convert), stringsAsFactors = FALSE)
    coef3 <-  data.table::rbindlist(list(coef3, meanx), use.names = T)
    coef2 <- data.frame()
    # z <- function(x)
    #   (meanx[[2]] * x + meanx[[1]])
    # p <- qplot(dat[[1]],dat[[2]], data=dat, xlab="X-axis", ylab="Y-axis")+ylim(0,(max(dat[[2]], na.rm = TRUE)+100))
    # b <- stat_function(fun=z)
    # v <- print(p + b)
  }
  coef3
}


make_groups <- function(n, k){
  # return a factor of length n defining groups of size k, except the last one
  nb_groups <- floor(n/k)
  overhang <- n - k*nb_groups
  c(rep(1:nb_groups, each = k), rep(nb_groups+1, overhang))
}

# how many rows to group together
group_size <- 10
# how many top coefs to use for mean
mean_top <- 2

with_map <- function(df){
  map_df(colnames(df)[2:4],
         function(cur_col_name){
           df %>%
             mutate(group = make_groups(n(), k=group_size)) %>%
             group_by(group) %>%
             summarize(mod = list(lm(formula(paste0(cur_col_name,"~Time")))),
                       .groups = "drop") %>%
             mutate(intercept = map_dbl(mod, ~.x$coefficients[[1]]),
                    slope = map_dbl(mod, ~ .x$coefficients[[2]])) %>%
             summarize(mean_intercept = mean(sort(intercept, decreasing = TRUE)[1:mean_top]),
                       mean_slope = mean(sort(slope,decreasing = TRUE)[1:mean_top]),
                       .groups = "drop")
         }
  )
}

microbenchmark::microbenchmark(with_for_loop(ex_df),
                               with_map(ex_df))
# Unit: milliseconds
#                  expr     min       lq     mean  median       uq      max neval
#  with_for_loop(ex_df) 65.2575 73.36925 81.36406 77.4504 83.82025 223.8993   100
#       with_map(ex_df) 39.2748 42.54920 47.53608 45.7630 49.27740  92.4213   100

So the dplyr approach is 2x as fast, that's not very impressive. But I expect it might get better on a bigger dataset (because you have assignments in your for loop that I suspect might choke on big datasets).