prep the data

# define aesthetic
palette = wesanderson::wes_palette("Zissou1", 15, "continuous")

# import data & check
food.dataW <- read.csv(file.path("~/Dropbox (University of Oregon)/Berkman Lab/CHIVES/Papers/ROCpredict/Analyses/AllData_wide_v4.csv", fsep=""))
food.dataL <- read.csv(file.path("~/Dropbox (University of Oregon)/Berkman Lab/CHIVES/Papers/ROCpredict/Analyses/AllData_long_v4.csv", fsep=""))

# center all IVs. month does not need to be centered since 0 is already baseline.
attach(food.dataL)
food.dataL$LCLNCvmpfc_c <- c(scale(LCLNC_vmpfc, center=TRUE, scale=FALSE))
food.dataL$LCLNCstria_c <- c(scale(LCLNC_striatum, center=TRUE, scale=FALSE))
food.dataL$RCLCdlpfc_c <- c(scale(RCLC_dlpfc, center=TRUE, scale=FALSE))
food.dataL$RCLCifg_c <- c(scale(RCLC_IFG, center=TRUE, scale=FALSE))
food.dataL$RCLCdacc_c <- c(scale(RCLC_dacc, center=TRUE, scale=FALSE))
food.dataL$RCLCvmpfc_c <- c(scale(RCLC_vmpfc, center=TRUE, scale=FALSE))
food.dataL$RCLCparhip_c <- c(scale(RCLC_Lparahip, center=TRUE, scale=FALSE))
food.dataL$RCLClifg_c <- c(scale(RCLC_Lifg, center=TRUE, scale=FALSE))
food.dataL$RCLClsupra_c <- c(scale(RCLC_Lsupra, center=TRUE, scale=FALSE))
food.dataL$RCLClcereb_c <- c(scale(RCLC_Lcereb, center=TRUE, scale=FALSE))
food.dataL$RCLCna_c <- c(scale(RCLC_na, center=TRUE, scale=FALSE))
food.dataL$RCLCrsupra_c <- c(scale(RCLC_Rsupra, center=TRUE, scale=FALSE))
food.dataL$RCLCrpmf_c <- c(scale(RCLC_RpostmedFront, center=TRUE, scale=FALSE))
food.dataL$BMI_c <- c(scale(BMI, center=TRUE, scale=FALSE))
food.dataL$age_c <- c(scale(age, center=TRUE, scale=FALSE))

# create regulation brain activity composites
food.dataL$regAnatROIAvg <- (food.dataL$RCLCdacc_c + food.dataL$RCLCdlpfc_c + food.dataL$RCLCifg_c) / 3
food.dataL$regPeakROIAvg <- (food.dataL$RCLCparhip_c + food.dataL$RCLClifg_c + food.dataL$RCLClsupra_c + 
                               food.dataL$RCLClcereb_c + food.dataL$RCLCna_c + food.dataL$RCLCrsupra_c + 
                               food.dataL$RCLCrpmf_c) / 7

generate tables

define functions

# define modeling function
fit_mod = function(data){
  mod = lmerTest::lmer(y_value ~ month * x_value + BMI_c + age_c + gender + condition + (1 | Participant), data = data)
  return(mod)
}

make_table = function(df, predictors = ".*", short = FALSE) {
  if (short == TRUE) {
    df %>%
      filter(effect == "fixed") %>%
      rename("p" = p.value,
             `neural predictor` = x_var) %>%
      mutate(term = gsub("\\(Intercept\\)", "intercept", term),
             term = gsub("x_value", "neural predictor", term),
             term = gsub("_c", "", term),
             term = gsub("month", "time", term),
             term = gsub(":", " x ", term),
             term = gsub("sd__", "", term),
             `neural predictor` = gsub("LCLNC", " LC > LNC ", `neural predictor`),
             `neural predictor` = gsub("RCLC", " RC > LC ", `neural predictor`),
             `neural predictor` = gsub("_c", "", `neural predictor`),
             `neural predictor` = gsub("pfc", "PFC", `neural predictor`),
             `neural predictor` = gsub("stria", "striatum", `neural predictor`),
             `neural predictor` = gsub("lifg", "l IFG", `neural predictor`),
             `neural predictor` = gsub("ifg", "IFG", `neural predictor`),
             `neural predictor` = gsub("dacc", "dACC", `neural predictor`),
             `neural predictor` = gsub("parhip", "parahippocampal gyrus", `neural predictor`),
             `neural predictor` = gsub("rsupra", "r supramarginal gyrus", `neural predictor`),
             `neural predictor` = gsub("lsupra", "l supramarginal gyrus", `neural predictor`),
             `neural predictor` = gsub("lcereb", "l cerebellum", `neural predictor`),
             `neural predictor` = gsub("na$", "midbrain", `neural predictor`),
             `neural predictor` = gsub("rpmf", "r posterior MFG", `neural predictor`),
             `neural predictor` = gsub("regAnatROIAvg", "average anatomical ROIs", `neural predictor`),
             `neural predictor` = gsub("regAnatROIAvg", "average peak ROIs", `neural predictor`),
             `b [95% CI]` = ifelse(effect == "fixed",
                                   sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high), 
                                   sprintf("%.02f (variance)", estimate)),
             p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))) %>%
      select(`neural predictor`, term, `b [95% CI]`, p) %>%
      mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
      mutate_if(is.character, funs(ifelse(. == "NA" | is.na(.), "--", .))) %>%
      filter(grepl(!!(predictors), term))
    
  } else {
    df %>%
      rename("SE" = std.error,
             "t" = statistic,
             "p" = p.value,
             `neural predictor` = x_var) %>%
      mutate(term = gsub("\\(Intercept\\)", "intercept", term),
             term = gsub("x_value", "neural predictor", term),
             term = gsub("_c", "", term),
             term = gsub("month", "time", term),
             term = gsub(":", " x ", term),
             term = gsub("sd__", "", term),
             term = ifelse(effect == "ran_pars", tolower(group), term),
             `neural predictor` = gsub("LCLNC", " LC > LNC ", `neural predictor`),
             `neural predictor` = gsub("RCLC", " RC > LC ", `neural predictor`),
             `neural predictor` = gsub("_c", "", `neural predictor`),
             `neural predictor` = gsub("pfc", "PFC", `neural predictor`),
             `neural predictor` = gsub("stria", "striatum", `neural predictor`),
             `neural predictor` = gsub("lifg", "l IFG", `neural predictor`),
             `neural predictor` = gsub("ifg", "IFG", `neural predictor`),
             `neural predictor` = gsub("dacc", "dACC", `neural predictor`),
             `neural predictor` = gsub("parhip", "parahippocampal gyrus", `neural predictor`),
             `neural predictor` = gsub("rsupra", "r supramarginal gyrus", `neural predictor`),
             `neural predictor` = gsub("lsupra", "l supramarginal gyrus", `neural predictor`),
             `neural predictor` = gsub("lcereb", "l cerebellum", `neural predictor`),
             `neural predictor` = gsub("na$", "midbrain", `neural predictor`),
             `neural predictor` = gsub("rpmf", "r posterior MFG", `neural predictor`),
             `neural predictor` = gsub("regAnatROIAvg", "average anatomical ROIs", `neural predictor`),
             `neural predictor` = gsub("regPeakROIAvg", "average peak ROIs", `neural predictor`),
             effect = gsub("ran_pars", "random", effect),
             `b [95% CI]` = ifelse(effect == "fixed",
                                   sprintf("%.02f [%.02f, %.02f]", estimate, conf.low, conf.high), 
                                   sprintf("%.02f (variance)", estimate)),
             SE = sprintf("%.02f", SE),
             t = sprintf("%.02f", t),
             df = sprintf("%.02f", df),
             p = ifelse(p < .001, "< .001", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))) %>%
      select(`neural predictor`, effect, term, `b [95% CI]`, SE, t, df, p) %>%
      mutate_if(is.numeric, funs(ifelse(is.na(.), "--", .))) %>%
      mutate_if(is.character, funs(ifelse(. == "NA" | is.na(.), "--", .))) %>%
        filter(grepl(!!(predictors), term))
  }
}

row_nums = c(1:10, 21:30, 41:50, 61:70, 81:90, 101:110, 121:130, 141:150)
print_table = function(tables, outcome, row_nums) {
  
  tables %>%
  filter(y_var == outcome) %>%
  select(table) %>%
  unnest() %>%
  group_by(term) %>%
  mutate(model = row_number()) %>%
  select(model, `neural predictor`, term, everything(), -y_var) %>%
  knitr::kable() %>%
  kableExtra::kable_styling() %>%
  kableExtra::row_spec(row_nums, background = "#f5f5f5")
  
}

run models

Run the following multilevel model for each neural predictor and dependent variable:

lmer(y_value ~ month * x_value + BMI_c + age_c + gender + condition + (1 | Participant))

  • y_value = dependent variable
  • month = time in months
  • x_value = grand-mean centered neural predictor (independent variable of interest)
  • BMI_c = grand-mean centered BMI
  • age_c = grand-mean centered age
  • condition = intervention condition (not relevant for the present analyses)
# tidy dataset
data_tables = food.dataL %>%
  gather(y_var, y_value, FCI_unhelCrv, FCI_helCrv, FCI_unhelLike, FCI_helLike, HEItot, Kcal, Fvavg, Empty) %>%
  gather(x_var, x_value, contains("_c"), contains("Avg"), -BMI_c, -age_c) %>%
  select(Participant, month, BMI_c, age_c, gender, condition, contains("y"), contains("x"))

# run models
models = data_tables %>%
  group_by(y_var, x_var) %>%
  nest() %>%
  mutate(test = map(data, fit_mod)) %>%
  mutate(draws = map(test, broom.mixed::tidy, conf.int = TRUE)) %>%
  select(-data, -test) %>%
  unnest()

# make tables
tables = models %>%
  group_by(y_var) %>%
  nest() %>%
  mutate(table = map(data, make_table, short = FALSE))