Based on the response to this question that I posed I have something like the following:
library(tidyverse)
library(dplyr)
library(broom)
library(tidyr)
library(purrr)
dataset <- tibble(
y1=rnorm(n=100),
y2=rnorm(n=100),
x1=rnorm(n=100),
x2=rnorm(n=100),)
outcomes <- dataset %>%
select(y1,y2) %>% colnames
covars <- dataset %>%
select(x1,x2) %>% colnames
paramlist <- expand_grid(outcomes, covars)
paramlist %>%
rowwise %>%
mutate(mod = list(lm(reformulate(outcomes, covars), data = dataset)),
res = list(broom::tidy(mod)),
predicted=list(predict(mod)),
data=list(cbind(dataset,predicted)))
# A tibble: 4 x 6
# Rowwise:
#> outcomes covars mod res predicted data
#> <chr> <chr> <list> <list> <list> <list>
#> 1 y1 x1 <lm> <tibble [2 x 5]> <dbl [100]> <df [100 x 5]>
#> 2 y1 x2 <lm> <tibble [2 x 5]> <dbl [100]> <df [100 x 5]>
#> 3 y2 x1 <lm> <tibble [2 x 5]> <dbl [100]> <df [100 x 5]>
#> 4 y2 x2 <lm> <tibble [2 x 5]> <dbl [100]> <df [100 x 5]>
What I would like to do now is - for each combination of outcomes and covars - I'd like to calculate the mean or sd of the predicted value in data conditional on some value of x1. For example, x1 might be a treatment variable, and I'd like the adjusted mean of the outcome for those with x1=0. The tricky part seems to be that the outcome and conditioning variable differ across rows.