Get predicted value out of a betaRegression model when using `nest()`

Question

I'm running some regression analyses using a nested file. The code that I have is similar to this:

library(tidyverse)
library(purrr)
library(betareg)

# 1. Create dataframe ----
dt <- data.frame(Marker = as.factor(paste0('m', rep(seq(1,10), 10))),
                 Year = rep(1990:1999, each = 10),
                 Ahat = rnorm(100, 0.5, 0.1)) %>% 
  mutate(Group = case_when(
    Marker %in% c("m1", "m2", "m3") ~ "A", 
    Marker %in% c("m4", "m5", "m6") ~ "B", 
    Marker %in% c("m7", "m8")       ~ "C",
    TRUE ~ "D"))

# 2 Nesting ----
nested_dt <- dt %>% 
  group_by(Group) %>% 
  nest()  

# 3 Beta Regression Function
marker_model <- function(dt) {
  betareg(Ahat ~ Year, data = dt)
}
# 4 Run Reg Model
models <- nested_dt %>% mutate(mod = map(data, marker_model))

It works fine, but I'd like now to get predicted values from a different time frame (e.g., 2000 to 2010). When using a normal file (not nested), I could easily do this using:

fit <- betareg(formula = Ahat ~ Year, dt)
overtime <- data.frame(Year = seq(2000, 2010)) 
predict(fit, type = "response", newdata = overtime)

So, does anyone know how to use this predict function on nested files (models in this case)?

Thank you!

Answer 1

You can get predicted values using map, as you create the mod column with it. But extracting the predicted in a human readable manner is not easy.

models %>%
    mutate(predicted_values = map(mod, predict, type = 'response', newdata = overtime)) %>%
    select(Group, predicted_values) %>%
    unnest(predicted_values) %>%
    mutate(rn = row_number(Group)) %>%
    pivot_wider(names_from = Group, values_from = predicted_values)
# # A tibble: 11 x 5
#       rn     A     B     C     D
#    <int> <dbl> <dbl> <dbl> <dbl>
#  1     1 0.509 0.487 0.506 0.484
#  2     2 0.510 0.486 0.512 0.487
#  3     3 0.512 0.485 0.519 0.489
#  4     4 0.514 0.483 0.525 0.492
#  5     5 0.515 0.482 0.531 0.495
#  6     6 0.517 0.481 0.537 0.498
#  7     7 0.519 0.479 0.544 0.501
#  8     8 0.520 0.478 0.550 0.504
#  9     9 0.522 0.477 0.556 0.507
# 10    10 0.524 0.476 0.562 0.510
# 11    11 0.525 0.474 0.568 0.513

Or just:

x <- models %>%
    mutate(pred = map(mod, predict, type = 'response', newdata = overtime))
as.data.frame(setNames(x$pred, x$Group))
#            A         B         C         D
# 1  0.5087058 0.4872116 0.5060286 0.4836964
# 2  0.5103726 0.4859097 0.5123000 0.4865974
# 3  0.5120391 0.4846079 0.5185676 0.4894992
# 4  0.5137054 0.4833064 0.5248293 0.4924018
# 5  0.5153713 0.4820050 0.5310832 0.4953049
# 6  0.5170369 0.4807040 0.5373274 0.4982083
# 7  0.5187021 0.4794031 0.5435598 0.5011118
# 8  0.5203670 0.4781026 0.5497787 0.5040153
# 9  0.5220313 0.4768024 0.5559821 0.5069185
# 10 0.5236952 0.4755024 0.5621680 0.5098212
# 11 0.5253585 0.4742029 0.5683347 0.5127233