How to solve polynomial equation in tidyverse framework?

Question

I want to apply groupwise polynomial model and solve it to get the x value when the y is highest individually for each group. I could able to do it for one group at a time like

library(tidyverse)

#Create soem data
df = data.frame(Group = c("A", "A","A","A","A","A","A","A","A","A","A","A",
                            "B", "B","B","B","B","B","B","B","B","B","B","B"),
       price = c(448, 560, 575, 500, 612, 600, 610, 590, 589, 532, 577, 560,
                 454, 568, 584, 506, 617, 608, 617, 599, 598, 537, 583, 567),
       score = c(83, 86, 89, 85, 89, 90, 91, 91, 91, 85, 88, 93,
                 89, 94, 98, 91, 94, 98, 98, 100, 100, 90, 94, 100))

#Subset group A
data <- df %>% 
  subset(Group %in% "A")

#Plot it
with(data, plot(x = score, y = price,
                pch=19, xlab = "Score", ylab = "Price"))

#Develop the polynomial model
model <- lm(price ~ score + I(score^2), data = data)

xx <- seq(min(data$score, na.rm = T), max(data$score, na.rm = T), 0.01)

lines(xx,  model$coefficients[[3]]*xx^2 + model$coefficients[[2]]*xx + 
        model$coefficients[[1]], col='blue', lwd = 3)

#Optimum score calculation
xoptimum <- -model$coefficients[[2]]/2/model$coefficients[[3]]
abline(v=xoptimum, col="red")

Now, how can I fit and solve the polynomial model for both groups using tidyverse framework? I have many groups in my original dataset.

Answer 1

You can use nest and map to work on groups of data:

library(tidyverse)

df %>%
  nest(data = -Group) %>%
  mutate(model = map(data, ~ lm(price ~ score + I(score^2), data = .x))) %>%
  mutate(x_opt = unlist(map(model, ~-.x$coef[[2]]/2/.x$coef[[3]]))) %>%
  mutate(y_opt = unlist(map2(x_opt, model, ~predict(.y, list(score = .x))))) %>%
  select(Group, x_opt, y_opt)
#> # A tibble: 2 x 3
#> Group x_opt y_opt
#> <chr> <dbl> <dbl>
#> 1 A   89.9  597.
#> 2 B   97.0  607.

To plot it with ggplot, we can instead do

df %>%
  nest(data = -Group) %>%
  mutate(model = map(data, ~ lm(price ~ score + I(score^2), data = .x))) %>%
  mutate(x_opt = unlist(map(model, ~-.x$coef[[2]]/2/.x$coef[[3]]))) %>%
  mutate(y_opt = unlist(map2(x_opt, model, ~predict(.y, list(score = .x))))) %>%
  unnest(data) %>%
  ggplot(aes(score, price)) +
  geom_point() +
  geom_smooth(method = lm, formula = y ~ x + I(x^2), se = FALSE) +
  geom_vline(aes(xintercept = x_opt), color = "red") +
  geom_hline(aes(yintercept = y_opt), linetype = 2) +
  theme_bw() +
  facet_grid(.~Group)