Robustness-Methods-and-Results.Rmd

---
output: latex_fragment
editor_options:
  chunk_output_type: console
bibliography: references.bib
csl: apa.csl
---

```{r}
#| label: setup
#| include: false

# R Studio Clean-Up
cat("\014") # clear console
# rm(list=ls()) # clear workspace - use restart R instead [cmd/alt + shift + F10]
gc() # garbage collector

# Install and Load Packages
# !IMPORTANT!
# BEFORE FIRST RENDER:
# To install all relevant packages please run "renv::restore()" (or renv::init() and then initiate from lockfile) in the console before the first use to ensure that all packages are using the correct version.
# to store the packages in a contained library within the project folder: renv::settings$use.cache(FALSE) and add 'RENV_CONFIG_SANDBOX_ENABLED = FALSE' to an '.Renviron' file
lib <- c(
  "rmarkdown",
  "knitr",
  "berryFunctions",
  "remedy",
  "bookdown",
  "brms",
  "psych",
  "ggplot2",
  "ggthemes",
  "purrr",
  "haven",
  "RColorBrewer",
  "plotly",
  "gridExtra",
  "ggpattern",
  "lme4",
  "nlme",
  "jtools",
  "gtsummary",
  "sessioninfo",
  "tibble",
  "pander",
  "devtools",
  "mada",
  "data.table",
  "plyr",
  "dplyr",
  "tidyr",
  "Hmisc",
  "kableExtra",
  "papaja",
  "stringr",
  "stringi",
  "reshape2",
  "lubridate",
  "misty",
  "sjPlot",
  "sjmisc",
  "metafor"
)
invisible(lapply(lib, library, character.only = TRUE))
rm(lib)

# Load Custom Packages
source("./scripts/functions/fun.panel.R")
source("./scripts/functions/themes.R")
source("./scripts/functions/binaryCor.R")
source("./scripts/functions/MlCorMat.R")
source("./scripts/functions/MlTbl.R")
source("./scripts/functions/metaLmer.R")

# Markdown Options
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file()) # set working directory
knitr::opts_knit$get("root.dir") # check working directory
options(
  scipen = 999,
  digits = 4,
  width = 400
) # removes scientific quotation
# knitr::opts_chunk$set(echo = TRUE, cache = F, cache.path = rprojroot::find_rstudio_root_file('cache/')) # cache settings

# Global Chunk Options
knitr::opts_chunk$set(
  fig.width = 12,
  fig.height = 8,
  fig.path = "Figures/",
  echo = TRUE,
  warning = FALSE,
  message = FALSE
)
```

```{r}
#| label: import
#| include: false

load("data/wrangled.RData")
# remove accidental surplus to avoid conflicts
# rm(list = ls(pattern = "mdl|medicalSdtComp"))
```

As with the main analyses, full surveys are available in our OSF repository \citep{KreienkampMasked2022a} and the full data description is available in Online Supplementary Material A. Correlations and descriptive statistics of the included variables are available in \tblref{tab:descrFullWide} and \tblref{tab:descrOutWide}.

\subsection{Additional Materials}

In addition to the measurement of whether or not participants had an intergroup interaction and their situational core need fulfillment, we also included a number of variables that allowed us to assess the robustness of our results.  

\subsubsection{Specific Psychological Needs} In addition to the intergroup contact dummy and situational core need reported in the main text, we included a common measure of three self-determination theory needs \citep[see][]{Downie2008}. The measurement was identical in all three studies. The items were introduced either by “\textit{During the interaction:}” or “\textit{This morning [/afternoon]:}” and measured autonomy (“\textit{I was myself.}”), competence (“\textit{I felt competent.}”), and relatedness (without intergroup contact “\textit{I had a strong need to belong}”; with intergroup contact: “\textit{I shared information about myself.}” and “\textit{The other(s) shared information about themselves.}”). All items were rated on a continuous slider scale from very little (-50) to a great deal (+50).

\subsubsection{Interaction Intent} To assess whether an interaction was accidental (vs. planned), we asked participants with a single item to report the extend to which "The interaction with -X- was accidental". The respondents were asked to report this context variable for all interactions they reported on using a continuous slider ranging from "not at all" (0), through "very little" (33) and "somewhat" (66), to "a great deal" (100). In all studies the scale showed a right skew (\textit{mean} = `r format(round(mean(dt_goal_directedness$InteractionAccidental, na.rm = TRUE), 2), nsmall = 2)`, \textit{sd} = `r format(round(sd(dt_goal_directedness$InteractionAccidental, na.rm = TRUE), 2), nsmall = 2)`).

\subsubsection{Goal-directedness}
```{r}
#| label: goal_directedness_descr
#| include: false

kappa.full <- function(x, name, ...){
  xtab <- table(x)
  diagonal.counts <- diag(xtab)
  N <- sum(xtab)
  row.marginal.props <- rowSums(xtab)/N
  col.marginal.props <- colSums(xtab)/N
  # Compute kappa (k)
  Po <- sum(diagonal.counts)/N
  Pe <- sum(row.marginal.props*col.marginal.props)
  k <- (Po - Pe)/(1 - Pe)
  se <- sqrt((Po*(1-Po))/(N*(1-Pe)^2)) # according to PMID: 23092060
  lwr <- k-1.96*se
  upr <- k+1.96*se
  
  data.frame(name,Po,Pe,k,se,lwr,upr)
}

irr_practical <- dt_topic_coding %>%
  select(practical_need_01, 
         practical_need_02) %>%
  kappa.full(., "practical_need")

irr_practical_dum <- dt_topic_coding %>%
  transmute(
    practical_need_01_dum = ifelse(practical_need_01 == 2, 1, practical_need_01), 
    practical_need_02_dum = ifelse(practical_need_02 == 2, 1, practical_need_02)
  ) %>%
  kappa.full(., "practical_need_dum")
  

irr_psychological <- dt_topic_coding %>%
  select(psychological_need_01, 
         psychological_need_02) %>%
  kappa.full(., "psychological_need")

irr_psychological_dum <- dt_topic_coding %>%
  transmute(
    practical_need_01_dum = ifelse(psychological_need_01 == 2, 1, psychological_need_01), 
    practical_need_02_dum = ifelse(psychological_need_02 == 2, 1, psychological_need_02)
  ) %>%
  kappa.full(., "psychological_need_dum")

irr <- rbind(
  irr_practical,
  irr_practical_dum,
  irr_psychological, 
  irr_psychological_dum
)

dt_topic_coding %>%
  transmute(
    practical_need_dum = ifelse(practical_need == 2, 1, practical_need), 
    psychological_need_dum = ifelse(psychological_need == 2, 1, psychological_need)
  )

lab_mapping <- c("no need", "need")
goal_directedness_dist <- dt_goal_directedness %>%
  transmute(
    practical_need_dum = ifelse(practical_need == 2, 1, practical_need), 
    psychological_need_dum = ifelse(psychological_need == 2, 1, psychological_need)
  ) %>% 
  mutate(
    practical_need_label = lab_mapping[match(practical_need_dum, c(0, 1))],
    psychological_need_label = lab_mapping[match(psychological_need_dum, c(0, 1))]
  ) %>%
  group_by(
    practical_need_dum,
    practical_need_label,
    psychological_need_dum,
    psychological_need_label
  ) %>%
  summarise(count = n()) %>%
  ungroup() %>%
  na.omit %>%
  mutate(
    perc = count/sum(count)*100
  )

ggplot(
  goal_directedness_dist,
  aes (
    x = reorder(practical_need_label, practical_need_dum),
    y = reorder(psychological_need_label, psychological_need_dum),
    fill = count,
    label = paste0(format(round(perc, 2), nsmall = 2), "%")
  )) +
  geom_tile() +
  geom_text(col = "white") +
  labs(x = "Practical Motive",
       y = "Psychological Motive",
       title = "Goal Directedness Coding") +
  theme(legend.key.size = unit(0.75, "cm"))
```
To assess whether the need content (i.e., the motives) would impact the effect of the need fulfillment experiences, we manually coded the topics we extracted during the topic modeling on two dimensions of how much they reflect a practical and a psychological goal-directedness. We chose practical and psychological needs specifically as our dimensions to account for differences in the types of needs that participants commonly reported. With practical motives we refer to specific, tangible goals or tasks that participants aimed to accomplish during the interaction. These instrumental goals are usually observable, concrete, and often centered on external outcomes, such as acquiring resources, completing tasks, and addressing immediate challenges \citep[e.g.,][]{oduntan2019}. With psychological motives we refer to underlying motives or desires that are more abstract and relate to personal fulfillment and well-being. In contrast to practical needs, psychological needs delve into the subjective and internal aspects of human experiences. These needs pertain to emotions, social connections, and cognitive processes, reflecting individuals' quest for personal growth, well-being, and thriving in social relationships \citep[][]{dweck2017}. Note that with this approach any particular motive can include a practical and/or a psychological goal-directedness but can also be classified as not having any goal at all. The full coding protocol we developed with examples for each of the codes is available in our Online Supplemental Material D. After an initial training, each of the two coders independently coded the 47 topics on the two dimensions, using one of three options each (i.e., 0 = no goal, 1 = vague goal, 2 = clear goal). Inter-rater reliability assessments showed that for both the practical as well as the psychological needs, agreement was not optimal using the three answer options (agreement practical = `r irr$Po[irr$name == "practical_need"] %>% round(2) %>% format(nsmall = 2)`\%, agreement psychological = `r irr$Po[irr$name == "psychological_need"] %>% round(2) %>% format(nsmall = 2)`\%). However, most disagreements were, if a need was present, whether that need was vague (1) or concrete (2). We, thus, collapsed these two categories making the ratings binary (need absent vs. need present). With the simpler coding, inter-rater agreement for practical needs (`r irr$Po[irr$name == "practical_need_dum"] %>% round(2) %>% format(nsmall = 2)`\%) and the psychological need (`r irr$Po[irr$name == "psychological_need_dum"] %>% round(2) %>% format(nsmall = 2)`\%) were much more reliable. Using Cohen’s $\kappa$ as our measure of inter-rater reliability, we find that both the practical need coding (\textit{Cohen's} $\kappa$ = `r irr$k[irr$name == "practical_need_dum"] %>% round(2) %>% format(nsmall = 2)`, 95\%CI[`r irr$lwr[irr$name == "practical_need_dum"] %>% round(2) %>% format(nsmall = 2)`, `r ifelse(irr$upr[irr$name == "practical_need_dum"]>1, 1 %>% round(2) %>% format(nsmall = 2), irr$upr[irr$name == "practical_need_dum"] %>% round(2) %>% format(nsmall = 2))`]) as well as the psychological need coding (\textit{Cohen's} $\kappa$ = `r irr$k[irr$name == "psychological_need_dum"] %>% round(2) %>% format(nsmall = 2)`, 95\%CI[`r irr$lwr[irr$name == "psychological_need_dum"] %>% round(2) %>% format(nsmall = 2)`, `r ifelse(irr$upr[irr$name == "psychological_need_dum"]>1, 1 %>% round(2) %>% format(nsmall = 2), irr$upr[irr$name == "psychological_need_dum"] %>% round(2) %>% format(nsmall = 2))`]) were very good. We thus proceeded with this collapsed coding. After resolving coder disagreements and merging the codings back to the free-text responses, we found that a majority of responses showed both a practical as well as a psychological need (`r goal_directedness_dist$perc[goal_directedness_dist$practical_need_dum==1 & goal_directedness_dist$psychological_need_dum==1] %>% round(2) %>% format(nsmall=2)`\%) and only few responses had no goal at all (`r goal_directedness_dist$perc[goal_directedness_dist$practical_need_dum==0 & goal_directedness_dist$psychological_need_dum==0] %>% round(2) %>% format(nsmall=2)`\%) with the remaining `r goal_directedness_dist$perc[goal_directedness_dist$practical_need_dum+goal_directedness_dist$psychological_need_dum==1] %>% sum %>% round(2) %>% format(nsmall=2)`\% having either a practical or a psychological need only (see Online Supplemental Material A for more detailed tables and visualizations).

\subsubsection{Well-being} We measured experienced well-being using a visual analog scale adapted from \citet{davies2022}. Participants were asked to respond to the the question "How do you feel right now?" using a continuous visual slider ranging from "very sad"(-100) to "very happy" (100). The well-being ratings were generally normally distributed (\textit{mean} = `r format(round(mean(dt_goal_directedness$WellBeing, na.rm = TRUE), 2), nsmall = 2)`, \textit{sd} = `r format(round(sd(dt_goal_directedness$WellBeing, na.rm = TRUE), 2), nsmall = 2)`).

\subsection{Results}

To build further confidence in our results, we assessed a number of additional models that might offer alternative explanations. We will discuss the results in sequential order --- in every case first considering the a global test of the model across the three studies and only then assessing whether the global three-level regression model suppresses any important person-level variations within the studies.

\subsubsection{Contact specific}

```{r}
#| label: contact_overall
#| include: false

dt_exp_all <- dt_exp_all %>%
  mutate(
    OutgroupInteraction = (OutgroupInteraction-1.5)
  )

contact_overall <- select_best_lmer_model(
  data = dt_exp_all,
  prediction_form = "CoreNeed * OutgroupInteraction",
  structure = "study/PID",
  dependent_var = "Attitude",
  lmer_scale = TRUE
)

# contact_overall$anova
# contact_overall$choice
# contact_overall$lmer_model %>% summ(., confint = TRUE, digits = 3)
# contact_overall$lmer_model %>% summ(., confint = TRUE, digits = 3, scale = TRUE)

# interactions::sim_slopes(mdl_contact$lmerAttSlopeOverall, pred = CoreNeedC, modx = OutgroupInteraction, jnplot = TRUE)
# interactions::interact_plot(contact_overall$lmer_model, pred = CoreNeed, modx = OutgroupInteraction, interval = TRUE)
```

We begin our robustness analysis by testing whether the effect of core need fulfillment is specific to an actual outgroup contact, rather than need fulfillment in general. For this, we analyzed the generalized situational core need fulfillment (either during a contact or about the daytime in general) and tested whether the effect differed during experience sampling measurements with and without outgroup contacts. We start this test by assessing the effect across all three studies, using a three-level hierarchical model, where measurements are nested within participants, and participants are nested within studies. In this overall model, we found no main effect of core need fulfillment (random slopes model, grand-mean standardized to account for all levels of variance; `r get_latex_coef(contact_overall$lmer_coef, "CoreNeed")`) but a significant interaction effect of core need fulfillment and outgroup contact (`r get_latex_coef(contact_overall$lmer_coef, "CoreNeed:OutgroupInteraction")`; also see \tblref{tab:robustnessTblLong}). While the three-level hierarchical model can be sensitive to scaling issues, this already indicates that it is not key need fulfillment in general --- but only key need fulfillment during an outgroup contact that predicts more positive outgroup attitudes. 

```{r}
#| label: contact_s1
#| include: false

contact_s1 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 1),
  prediction_form = "CoreNeedC * OutgroupInteraction",
  structure = "PID",
  dependent_var = "Attitude"
)
# contact_s1$anova
# contact_s1$choice
# contact_s1$lmer_model %>% summ(., confint = TRUE, digits = 3)
# contact_s1$lmer_model_z %>% summ(., confint = TRUE, digits = 3)

# interactions::sim_slopes(mdl_contact$lmerAttSlopeOverall, pred = CoreNeedC, modx = OutgroupInteraction, jnplot = TRUE)
# interactions::interact_plot(contact_s1$lmer_model, pred = CoreNeedC, modx = OutgroupInteraction, interval = TRUE)
```

```{r}
#| label: contact_s2
#| include: false

contact_s2 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 2),
  prediction_form = "CoreNeedC * OutgroupInteraction",
  structure = "PID",
  dependent_var = "Attitude"
)
# contact_s2$anova
# contact_s2$choice
# contact_s2$lmer_model %>% summ(., confint = TRUE, digits = 3)
# contact_s2$lmer_model_z %>% summ(., confint = TRUE, digits = 3)

# interactions::sim_slopes(mdl_contact$lmerAttSlopeOverall, pred = CoreNeedC, modx = OutgroupInteraction, jnplot = TRUE)
# interactions::interact_plot(contact_s2$lmer_model, pred = CoreNeedC, modx = OutgroupInteraction, interval = TRUE)
```

```{r}
#| label: contact_s3
#| include: false

contact_s3 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 3),
  prediction_form = "CoreNeedC * OutgroupInteraction",
  structure = "PID",
  dependent_var = "Attitude"
)
# contact_s3$anova
# contact_s3$choice
# contact_s3$lmer_model %>% summ(., confint = TRUE, digits = 3)
# contact_s3$lmer_model_z %>% summ(., confint = TRUE, digits = 3)

# interactions::sim_slopes(mdl_contact$lmerAttSlopeOverall, pred = CoreNeedC, modx = OutgroupInteraction, jnplot = TRUE)
# interactions::interact_plot(contact_s3$lmer_model, pred = CoreNeedC, modx = OutgroupInteraction, interval = TRUE)
```

```{r}
#| label: contact_study_coef
#| include: false

contact_study_coef <- 
  rbind(
    contact_s1$lmer_coef %>% mutate(study = 1),
    contact_s2$lmer_coef %>% mutate(study = 2),
    contact_s3$lmer_coef %>% mutate(study = 3)
  )
```

To ensure that the results are not affected by scaling issues (e.g., study-level variances suppressing person-level variances) or a similar Simpson's paradox, we additionally assess the model within each of the three studies. Within each of the three studies, the effects are more pronounced, so that we also see a significant effect of core need fulfillment (all \textit{b} > `r min(abs(contact_study_coef[contact_study_coef$coef == "CoreNeedC", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r max(contact_study_coef[contact_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3)`) as well as outgroup contact itself (all |\textit{b}| > `r min(abs(contact_study_coef[contact_study_coef$coef == "OutgroupInteraction", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r max(contact_study_coef[contact_study_coef$coef == "OutgroupInteraction", "p"]) %>% round(3) %>% format(nsmall=3)`) but the interaction effect consistently remains the most reliable predictor of outgroup attitudes (all |\textit{b}| > `r min(abs(contact_study_coef[contact_study_coef$coef == "CoreNeedC:OutgroupInteraction", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r max(contact_study_coef[contact_study_coef$coef == "CoreNeedC:OutgroupInteraction", "p"]) %>% round(3) %>% format(nsmall=3)`, also see \tblref{tab:robustnessTblLong}). There is thus consistent evidence that need fulfillment relates to outgroup attitudes for outgroup contacts in particular but not need fulfillment in general. 


\subsubsection{Interaction intent}

```{r}
#| label: accidental_overall
#| include: false

accidental_overall <- select_best_lmer_model(
  data = dt_exp_all %>% filter(OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeed * InteractionAccidental",
  structure = "study/PID",
  dependent_var = "Attitude",
  lmer_scale = TRUE
)
# accidental_overall$anova
# accidental_overall$choice
# accidental_overall$lmer_model %>% summ(., confint = TRUE, digits = 3)
# accidental_overall$lmer_model %>% summ(., confint = TRUE, digits = 3, scale = TRUE)
```

Secondly, to assess whether the need fulfillment mechanism affected by whether the interaction was accidental or planned we ran an exploratory moderation analysis using the participants' ratings of how much they perceived the interaction as 'accidental'. It should be noted that we asked our participants to focus on most the important interaction (i.e., "\textit{The following questions will be about the interaction \underline{you consider most significant}.}"; emphasis as in original). We again start our analysis approach by assessing the model across all three studies, using a three-level hierarchical model. In this overall model, we retain the main effect of core need fulfillment (random slopes model, grand-mean standardized to account for all levels of variance; `r get_latex_coef(accidental_overall$lmer_coef, "CoreNeed")`) but neither contact intent nor the moderation effect affect the results (all |\textit{b}| < `r max(abs(accidental_overall$lmer_coef[accidental_overall$lmer_coef$coef != c("(Intercept)", "CoreNeed"), "b"])) %>% round(2) %>% format(nsmall=2)` and all \textit{p} > `r min(accidental_overall$lmer_coef[accidental_overall$lmer_coef$coef != c("(Intercept)", "CoreNeed"), "p"]) %>% round(3) %>% format(nsmall=3)`; see \tblref{tab:robustnessTblLong} for full results). 

```{r}
#| label: accidental_s1
#| include: false

accidental_s1 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 1, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC * InteractionAccidentalC",
  structure = "PID",
  dependent_var = "Attitude"
)
# accidental_s1$anova
# accidental_s1$choice
# accidental_s1$lmer_model %>% summ(., confint = TRUE, digits = 3)
# accidental_s1$lmer_model_z %>% summ(., confint = TRUE, digits = 3)

# interactions::sim_slopes(accidental_s1$lmer_model, pred = CoreNeedC, modx = InteractionAccidentalC, jnplot = TRUE)
# interactions::interact_plot(accidental_s1$lmer_model, pred = CoreNeedC, modx = InteractionAccidentalC, interval = TRUE)
```

```{r}
#| label: accidental_s2
#| include: false

accidental_s2 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC * InteractionAccidentalC",
  structure = "PID",
  dependent_var = "Attitude",
  return = "slope"
)
# accidental_s2$anova
# accidental_s2$choice
# accidental_s2$lmer_model %>% summ(., confint = TRUE, digits = 3)
# accidental_s2$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
# 
# interactions::sim_slopes(accidental_s2$lmer_model, pred = CoreNeedC, modx = InteractionAccidentalC, jnplot = TRUE)
# interactions::interact_plot(accidental_s2$lmer_model, pred = CoreNeedC, modx = InteractionAccidentalC, interval = TRUE)
```

```{r}
#| label: accidental_s3
#| include: false

accidental_s3 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 3, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC * InteractionAccidentalC",
  structure = "PID",
  dependent_var = "Attitude"
)
# accidental_s3$anova
# accidental_s3$choice
# accidental_s3$lmer_model %>% summ(., confint = TRUE, digits = 3)
# accidental_s3$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
# 
# interactions::sim_slopes(accidental_s3$lmer_model, pred = CoreNeedC, modx = InteractionAccidentalC, jnplot = TRUE)
# interactions::interact_plot(accidental_s3$lmer_model, pred = CoreNeedC, modx = InteractionAccidentalC, interval = TRUE)
```

```{r}
#| label: accidental_study_coef
#| include: false

accidental_study_coef <- 
  rbind(
    accidental_s1$lmer_coef %>% mutate(study = 1),
    accidental_s2$lmer_coef %>% mutate(study = 2),
    accidental_s3$lmer_coef %>% mutate(study = 3)
  )
```

We again sought to ensure that the results were not affected by scaling issues by additionally assessing the interaction intentionality model within each of the three studies. Within each of the three studies, the effect of core need fulfillment became even clearer (all |\textit{b}| > `r min(abs(accidental_study_coef[accidental_study_coef$coef == "CoreNeedC", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r ifelse(max(accidental_study_coef[accidental_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3)<.001, ".001", max(accidental_study_coef[accidental_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3))`). But in none of the studies neither outgroup contact intention nor the moderation effect explained a significant amount of variance in outgroup attitudes (all |\textit{b}| < `r max(abs(accidental_study_coef[accidental_study_coef$coef != c("(Intercept)", "CoreNeedC"), "b"])) %>% round(2) %>% format(nsmall=2)` and all \textit{p} > `r min(accidental_study_coef[accidental_study_coef$coef != c("(Intercept)", "CoreNeedC"), "p"]) %>% round(3) %>% format(nsmall=3)`; also see \tblref{tab:robustnessTblLong}). There is thus consistent evidence that need fulfillment is related to outgroup attitudes, even when taking the intentionality of the interaction into account --- at least in our three samples and with a focus on the most significant interactions.


\subsubsection{Well-being outcome}

```{r}
#| label: wellbeing_core_overall
#| include: false

wellbeing_core_overall <- select_best_lmer_model(
  data = dt_exp_all %>% filter(OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeed",
  structure = "study/PID",
  dependent_var = "WellBeing",
  lmer_scale = TRUE
)
# wellbeing_core_overall$anova
# wellbeing_core_overall$choice
# wellbeing_core_overall$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_core_overall$lmer_model %>% summ(., confint = TRUE, digits = 3, scale = TRUE)
```

```{r}
#| label: wellbeing_full_overall
#| include: false

wellbeing_full_overall <- select_best_lmer_model(
  data = dt_exp_all %>% filter(OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeed + Quality",
  structure = "study/PID",
  dependent_var = "WellBeing",
  lmer_scale = TRUE
)
# wellbeing_full_overall$anova
# wellbeing_full_overall$choice
# wellbeing_full_overall$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_full_overall$lmer_model %>% summ(., confint = TRUE, digits = 3, scale = TRUE)
```


```{r}
#| label: wellbeing_core_s1
#| include: false

wellbeing_core_s1 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 1, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC",
  structure = "PID",
  dependent_var = "WellBeing"
)
# wellbeing_core_s1$anova
# wellbeing_core_s1$choice
# wellbeing_core_s1$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_core_s1$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: wellbeing_full_s1
#| include: false

wellbeing_full_s1 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 1, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC + QualityC",
  structure = "PID",
  dependent_var = "WellBeing"
)
# wellbeing_full_s1$anova
# wellbeing_full_s1$choice
# wellbeing_full_s1$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_full_s1$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```


```{r}
#| label: wellbeing_core_s2
#| include: false

wellbeing_core_s2 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC",
  structure = "PID",
  dependent_var = "WellBeing"
)
# wellbeing_core_s2$anova
# wellbeing_core_s2$choice
# wellbeing_core_s2$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_core_s2$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: wellbeing_full_s2
#| include: false

wellbeing_full_s2 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC + QualityC",
  structure = "PID",
  dependent_var = "WellBeing"
)
# wellbeing_full_s2$anova
# wellbeing_full_s2$choice
# wellbeing_full_s2$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_full_s2$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```


```{r}
#| label: wellbeing_core_s3
#| include: false

wellbeing_core_s3 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 3, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC",
  structure = "PID",
  dependent_var = "WellBeing"
)
# wellbeing_core_s3$anova
# wellbeing_core_s3$choice
# wellbeing_core_s3$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_core_s3$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: wellbeing_full_s3
#| include: false

wellbeing_full_s3 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 3, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC + QualityC",
  structure = "PID",
  dependent_var = "WellBeing"
)
# wellbeing_full_s3$anova
# wellbeing_full_s3$choice
# wellbeing_full_s3$lmer_model %>% summ(., confint = TRUE, digits = 3)
# wellbeing_full_s3$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: wellbeing_core_study_coef
#| include: false

wellbeing_core_study_coef <- 
  rbind(
    wellbeing_core_s1$lmer_coef %>% mutate(study = 1),
    wellbeing_core_s2$lmer_coef %>% mutate(study = 2),
    wellbeing_core_s3$lmer_coef %>% mutate(study = 3)
  )
```

Thirdly, to build a stronger case for the relevance of need fulfillment to minority group members, we exploratorily assessed the effect of need fulfilling outgroup interactions on self-reported well-being. We, thus, re-ran our main analysis but substituted the outgroup attitudes outcome with situational well-being. As with the previous robustness analyses, we begin with a global three-level hierarchical model (across the three studies). We find that need fulfillment during outgroup contacts, indeed, has a similar effect on experienced well-being (random slopes model, grand-mean standardized to account for all levels of variance; `r get_latex_coef(wellbeing_core_overall$lmer_coef, "CoreNeed")`). We found the same result when we assessed each of the three studies individually. In each of the studies situational need fulfillment during the outgroup interaction was related with higher well-being ratings by the participants (random slopes model, centered within participants; all |\textit{b}| > `r min(abs(wellbeing_core_study_coef[wellbeing_core_study_coef$coef == "CoreNeedC", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r ifelse(max(wellbeing_core_study_coef[wellbeing_core_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3)<.001, ".001", max(wellbeing_core_study_coef[wellbeing_core_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3))`). We, thus, find consistent and meaningful evidence that need fulfilling outgroup interactions also relate to higher everyday well-being. 

\subsubsection{Need types}

```{r}
#| label: need_type_overall
#| include: false

need_type_overall <- select_best_lmer_model(
  data = dt_goal_directedness,
  prediction_form = "CoreNeed * practical_need_dum + CoreNeed * psychological_need_dum",
  structure = "study/PID",
  dependent_var = "Attitude",
  lmer_scale = TRUE
)
# need_type_overall$anova
# need_type_overall$choice
# need_type_overall$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_type_overall$lmer_model %>% summ(., confint = TRUE, digits = 3, scale = TRUE)
```

Fourthly, to assess the role of different types of motives reported by our participants, we added our coding of practical and psychological goal-directedness as additional predictors to our base model. We thus had core need fulfillment predicting outgroup attitudes while also accounting for whether the reported motives were capturing practical and/or psychological motives. We again ran a global model, across the three studies first. We found that core need fulfillment remain a core predictor of outgroup attitudes (random slopes model, grand-mean standardized to account for all levels of variance; `r get_latex_coef(need_type_overall$lmer_coef, "CoreNeed")`), even after accounting for different types of motives. None of the motive types nor the moderation effects reached statistical significance within the overall analysis (all |\textit{b}| < `r max(abs(need_type_overall$lmer_coef[need_type_overall$lmer_coef$coef != c("(Intercept)", "CoreNeed"), "b"])) %>% round(2) %>% format(nsmall=2)` and all \textit{p} > `r min(need_type_overall$lmer_coef[need_type_overall$lmer_coef$coef != c("(Intercept)", "CoreNeed"), "p"]) %>% round(3) %>% format(nsmall=3)`; see \tblref{tab:robustnessTblLong} for full results). 

```{r}
#| label: need_type_s1
#| include: false

need_type_s1 <- select_best_lmer_model(
  data = dt_goal_directedness %>% filter(study == 1),
  prediction_form = "CoreNeedC * practical_need_dum + CoreNeedC * psychological_need_dum",
  structure = "PID",
  dependent_var = "Attitude",
  return = "slope"
)
# need_type_s1$anova
# need_type_s1$choice
# need_type_s1$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_type_s1$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: need_type_s2
#| include: false

need_type_s2 <- select_best_lmer_model(
  data = dt_goal_directedness %>% filter(study == 2),
  prediction_form = "CoreNeedC * practical_need_dum + CoreNeedC * psychological_need_dum",
  structure = "PID",
  dependent_var = "Attitude"
)
# need_type_s2$anova
# need_type_s2$choice
# need_type_s2$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_type_s2$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: need_type_s3
#| include: false

need_type_s3 <- select_best_lmer_model(
  data = dt_goal_directedness %>% filter(study == 3),
  prediction_form = "CoreNeedC * practical_need_dum + CoreNeedC * psychological_need_dum",
  structure = "PID",
  dependent_var = "Attitude"
)
# need_type_s3$anova
# need_type_s3$choice
# need_type_s3$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_type_s3$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: need_type_study_coef
#| include: false

need_type_study_coef <- 
  rbind(
    need_type_s1$lmer_coef %>% mutate(study = 1),
    need_type_s2$lmer_coef %>% mutate(study = 2),
    need_type_s3$lmer_coef %>% mutate(study = 3)
  )
```

When looking at the individual studies, we again saw that the effect of core need fulfillment remained the on only clear effect (all |\textit{b}| > `r min(abs(need_type_study_coef[need_type_study_coef$coef == "CoreNeedC", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r ifelse(max(need_type_study_coef[need_type_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3)<.001, ".001", max(need_type_study_coef[need_type_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3))`). Additionally, in none of the studies neither motive type dummies nor the moderation effect explained a significant amount of variance in outgroup attitudes (all |\textit{b}| < `r max(abs(need_type_study_coef[need_type_study_coef$coef != c("(Intercept)", "CoreNeedC"), "b"])) %>% round(2) %>% format(nsmall=2)` and all \textit{p} > `r min(need_type_study_coef[need_type_study_coef$coef != c("(Intercept)", "CoreNeedC"), "p"]) %>% round(3) %>% format(nsmall=3)`; also see \tblref{tab:robustnessTblLong}). We, thus, find consistent evidence that need fulfillment is related to outgroup attitudes, even when taking the type of need into account --- at least in our three samples.

\subsubsection{Specific psychological needs}

```{r}
#| label: need_sdt_overall
#| include: false

need_sdt_overall <- select_best_lmer_model(
  data = dt_exp_all %>% filter(OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeed + Competence + Autonomy + Relatednesss",
  structure = "study/PID",
  dependent_var = "Attitude",
  lmer_scale = TRUE
)
# need_sdt_overall$anova
# need_sdt_overall$choice
# need_sdt_overall$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_sdt_overall$lmer_model %>% summ(., confint = TRUE, digits = 3, scale = TRUE)
```

In a final step, we checked whether during the interaction the core situational need remains a meaningful predictor even when taking other fundamental psychological needs into account. We again take a two-step approach, starting with cross-study global three-level test and then assessing the effects within the individual studies. Within the overall model we find that across the studies core need fulfillment remained a strong predictor of outgroup attitudes, even after controlling for the three self-determination theory need (random slopes model, grand-mean standardized to account for all levels of variance; `r get_latex_coef(need_sdt_overall$lmer_coef, "CoreNeed")`). Within this overall analysis, none of the self-determination theory needs independently predicted outgroup attitudes to a statistically significant extent (all \textit{p} > `r min(need_sdt_overall$lmer_coef[need_sdt_overall$lmer_coef$coef != c("(Intercept)", "CoreNeed"), "p"]) %>% round(3) %>% format(nsmall=3)`). However, some of the effect sizes were largely comparable to that of the core need fulfillment (all |\textit{b}| < `r max(abs(need_sdt_overall$lmer_coef[need_sdt_overall$lmer_coef$coef != c("(Intercept)", "CoreNeed"), "b"])) %>% round(2) %>% format(nsmall=2)`, particularly that of relatedness fulfillment; see \tblref{tab:robustnessTblLong} for full results). 

```{r}
#| label: need_sdt_s1
#| include: false

need_sdt_s1 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 1, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC + CompetenceC + AutonomyC + RelatednesssC",
  structure = "PID",
  dependent_var = "Attitude"
)
# need_sdt_s1$anova
# need_sdt_s1$choice
# need_sdt_s1$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_sdt_s1$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: need_sdt_s2
#| include: false

# need_sdt_s2 <- select_best_lmer_model(
#   data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes"),
#   prediction_form = "CoreNeedC + CompetenceC + AutonomyC + RelatednesssC",
#   structure = "PID",
#   dependent_var = "Attitude",
#   control = lmeControl(opt = "nlminb", msMaxIter = 500, tolerance = 1e-6)
# )
model0 <- lmer(
  Attitude ~ 1 + (1 | PID),
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes")
)
model1 <- lmer(
  Attitude ~ CoreNeedC + CompetenceC + AutonomyC + RelatednesssC + (1 | PID),
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes")
)
model2 <- lmer(
  Attitude ~ CoreNeedC + CompetenceC + AutonomyC + RelatednesssC + (1 + CoreNeedC + CompetenceC + AutonomyC + RelatednesssC | PID),
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes")
)
test <- anova(model0, model1, model2)
test$`Pr(>Chisq)`[3] < 0.05
rm(model0, model1, model2, test)

need_sdt_s2 <- list()
need_sdt_s2$lmer_model <- lmer(
  Attitude ~ CoreNeedC + CompetenceC + AutonomyC + RelatednesssC + (1 | PID),
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes")
)
need_sdt_s2$lmer_model_z <- lmer(
  AttitudeZ ~ CoreNeedZ + CompetenceZ + AutonomyZ + RelatednesssZ + (1 | PID),
  data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes")
)

need_sdt_s2$lmer_model %>% summ(., confint = TRUE, digits = 3)
need_sdt_s2$lmer_model_z %>% summ(., confint = TRUE, digits = 3)

need_sdt_s2$lmer_coef <- 
    summ(need_sdt_s2$lmer_model, confint = TRUE, scale = FALSE)$coeftable %>%
    as_tibble %>%
    transmute(
      coef = rownames(summ(need_sdt_s2$lmer_model)$coeftable),
      b = Est.,
      df = d.f.,
      t = `t val.`,
      p = p,
      lwr = `2.5%`,
      upr = `97.5%`,
      coef_latex = paste0(
        "\\textit{b} = ", format(round(b, 2), nsmall = 2), 
        ", \\textit{t}(", format(df, big.mark=","), ") = ", format(round(t, 2), nsmall = 2),
        ", \\textit{p} ", ifelse(p < .001, "< .001", paste0("= ", format(round(p, 3), nsmall = 3))),
        ", \\textit{95\\%CI}[",  format(round(lwr, 2), nsmall = 2), ", ",  format(round(upr, 2), nsmall = 2), "]"
      ),
      coef_html = paste0(
        "<i>b</i> = ", format(round(b, 2), nsmall = 2), 
        ", <i>t</i>(", format(df, big.mark=","), ") = ", format(round(t, 2), nsmall = 2),
        ", <i>p</i> ", ifelse(p < .001, "< .001", paste0("= ", format(round(p, 3), nsmall = 3))),
        ", <i>95%CI</i>[",  format(round(lwr, 2), nsmall = 2), ", ",  format(round(upr, 2), nsmall = 2), "]"
      )
    )
```

```{r}
#| label: need_sdt_s3
#| include: false

need_sdt_s3 <- select_best_lmer_model(
  data = dt_exp_all %>% filter(study == 3, OutgroupInteractionLab == "Yes"),
  prediction_form = "CoreNeedC + CompetenceC + AutonomyC + RelatednesssC",
  structure = "PID",
  dependent_var = "Attitude"
)
# need_sdt_s3$anova
# need_sdt_s3$choice
# need_sdt_s3$lmer_model %>% summ(., confint = TRUE, digits = 3)
# need_sdt_s3$lmer_model_z %>% summ(., confint = TRUE, digits = 3)
```

```{r}
#| label: need_sdt_study_coef
#| include: false

need_sdt_study_coef <- 
  rbind(
    need_sdt_s1$lmer_coef %>% mutate(study = 1),
    need_sdt_s2$lmer_coef %>% mutate(study = 2),
    need_sdt_s3$lmer_coef %>% mutate(study = 3)
  )

need_sdt_study_coef_other <- need_sdt_study_coef %>%
  filter(
    coef != "(Intercept)",
    coef != "CoreNeedC"
  ) 
```


When looking at the individual studies, we again saw that core need fulfillment remained a consistent predictor of outgroup attitudes, even after accounting for the self-determination theory need (all |\textit{b}| > `r min(abs(need_sdt_study_coef[need_sdt_study_coef$coef == "CoreNeedC", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r ifelse(max(need_sdt_study_coef[need_sdt_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3)<.001, ".001", max(need_sdt_study_coef[need_sdt_study_coef$coef == "CoreNeedC", "p"]) %>% round(3) %>% format(nsmall=3))`). However, across all three studies the fulfillment of relatedness motives also emerged as a consistent predictor of outgroup attitudes (all |\textit{b}| > `r min(abs(need_sdt_study_coef[need_sdt_study_coef$coef == "RelatednesssC", "b"])) %>% round(2) %>% format(nsmall=2)`, all \textit{p} < `r ifelse(max(need_sdt_study_coef[need_sdt_study_coef$coef == "RelatednesssC", "p"]) %>% round(3) %>% format(nsmall=3)<.001, ".001", max(need_sdt_study_coef[need_sdt_study_coef$coef == "RelatednesssC", "p"]) %>% round(3) %>% format(nsmall=3))`). Additionally, in the larger studies 2 and 3 competence fulfillment was also related to more positive outgroup attitudes (study 2: `r get_latex_coef(need_sdt_study_coef %>% filter(study == 2), "CompetenceC")`, study 3: `r get_latex_coef(need_sdt_study_coef %>% filter(study == 3), "CompetenceC")`). None of the autonomy fulfillment effects reached statistical significance nor did the competence fulfillment during study 1 (see \tblref{tab:robustnessTblLong} for the full results). In short, find that across our samples, relatedness fulfillment (and to a smaller extend competence fulfillment) are instrumental in understanding when an outgroup contact leads to more positive outgroup attitudes. Importantly, even when considering these effects situational core need fulfillment remains a strong and consistent predictor of outgroup attitudes. In some cases, we even find that core need fulfillment takes on some of the variance that would otherwise be explained by the self-determination theory needs (see Online Supplemental Material A).


<!-- TABLES -->

```{r}
#| label: descrFullWideAppB
#| include: false

vars_full <- c(
  "Core Need" = "CoreNeed",
  "Competence" = "Competence",
  "Autonomy" = "Autonomy",
  "Relatednesss" = "Relatednesss",
  #"Outgroup Interaction" = "OutgroupInteraction",
  "Accidental" = "InteractionAccidental",
  "Quality" = "Quality",
  "Attitudes NL" = "Attitude"
)

dt_exp_all <- dt_exp_all %>% 
  mutate(
    OutgroupInteraction = (OutgroupInteraction-1.5)
  )

s1_MlCor_full <- 
  MlCorMat(
    data = dt_exp_all %>% filter(study == 1),
    id = "PID",
    selection = setdiff(as.character(vars_full), c("InteractionAccidental", "Quality")),
    labels = setdiff(names(vars_full), c("Accidental", "Quality"))
  ) 
s2_MlCor_full <- 
  MlCorMat(
    data = dt_exp_all %>% filter(study == 2),
    id = "PID",
    selection = as.character(vars_full),
    labels = names(vars_full)
  ) 
s3_MlCor_full <- 
  MlCorMat(
    data = dt_exp_all %>% filter(study == 3),
    id = "PID",
    selection = as.character(vars_full),
    labels = names(vars_full)
  ) 
all_MlCor_full <- 
  MlCorMat(
    data = dt_exp_all,
    id = "PID",
    selection = as.character(vars_full),
    labels = names(vars_full)
  ) 


rbind(
  s1_MlCor_full %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable") %>%
    add_column(Accidental = "",
               .before = "Attitudes NL") %>%
    add_column(Quality = "",
               .before = "Attitudes NL"),
  s2_MlCor_full %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable"),
  s3_MlCor_full %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable"),
  all_MlCor_full %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable")
) %>%
  kbl(.,
      format = "latex",
      caption = "Full Sample: Correlation Table and Descriptive Statistics",
      booktabs = TRUE,
      align = c("l", rep("c", ncol(.) - 1)))  %>%
  kable_classic() %>%
  add_header_above(
    .,
    c(
      "",
      "Correlations" = ncol(s3_MlCor_full),
      "Descriptives" = nrow(s3_MlCor_full) - ncol(s3_MlCor_full)
    )
  ) %>%
  pack_rows("Study 1", 1, ncol(s1_MlCor_full)) %>%
  pack_rows("Study 2", ncol(s1_MlCor_full) + 1,
            sum(ncol(s1_MlCor_full), ncol(s2_MlCor_full))) %>%
  pack_rows("Study 3", sum(ncol(s1_MlCor_full), ncol(s2_MlCor_full)) + 1,  
            sum(ncol(s1_MlCor_full), ncol(s2_MlCor_full), ncol(s3_MlCor_full))) %>%
  pack_rows("Across Studies", sum(ncol(s1_MlCor_full), ncol(s2_MlCor_full), ncol(s3_MlCor_full)) + 1,  
            sum(ncol(s1_MlCor_full), ncol(s2_MlCor_full), ncol(s3_MlCor_full), ncol(all_MlCor_full))) %>% # fix table position
  footnote(
    general = c(
      "Upper triangle: Within-person correlations;",
      "Lower triangle: Between-person correlations;",
      "*** p < .001, ** p < .01,  * p < .05",
      "Interaction Intent and Interaction Quality were only measured for outgroup interactions in Study 1."
    )
  ) %>%
  kable_styling(latex_options = "scale_down") %>%
  gsub("\\begin{table}", "\\begin{table}\n\\begin{minipage}[t][\\textheight][t]{\\textwidth}", ., fixed = TRUE) %>% # fix table position
  gsub("\\end{table}", "\\end{minipage}\n\\end{table}", ., fixed = TRUE) %>% # fix table position
  save_kable("Tables/descrFullWideAppB.tex")

```

```{r}
#| label: descrOutWideAppB
#| include: false

vars_out <- c(
  "Core Need" = "CoreNeed",
  "Competence" = "Competence",
  "Autonomy" = "Autonomy",
  "Relatednesss" = "Relatednesss",
  # "Practical Need" = "practical_need_dum",
  # "Psychological Need" = "psychological_need_dum",
  "Accidental" = "InteractionAccidental",
  "Quality" = "Quality",
  "Attitudes NL" = "Attitude"
)

s1_MlCor_out <- 
  MlCorMat(
    data = dt_exp_all %>% filter(study == 1, OutgroupInteractionLab == "Yes", PID != 21), # PP 21 has inflated zero for Accidental
    id = "PID",
    selection = as.character(vars_out),
    labels = names(vars_out)
  ) 
s2_MlCor_out <- 
  MlCorMat(
    data = dt_exp_all %>% filter(study == 2, OutgroupInteractionLab == "Yes"),
    id = "PID",
    selection = as.character(vars_out),
    labels = names(vars_out)
  ) 
s3_MlCor_out <- 
  MlCorMat(
    data = dt_exp_all %>% filter(study == 3, OutgroupInteractionLab == "Yes"),
    id = "PID",
    selection = as.character(vars_out),
    labels = names(vars_out)
  ) 
all_MlCor_out <- 
  MlCorMat(
    data = dt_exp_all %>% filter(OutgroupInteractionLab == "Yes"),
    id = "PID",
    selection = as.character(vars_out),
    labels = names(vars_out)
  ) 


rbind(
  s1_MlCor_out %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable"),
  s2_MlCor_out %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable"),
  s3_MlCor_out %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable"),
  all_MlCor_out %>%
    t %>%
    as.data.frame %>%
    rownames_to_column(., "Variable")
) %>%
  kbl(.,
      format = "latex",
      caption = "Outgroup Interaction Sample: Correlation Table and Descriptive Statistics",
      booktabs = TRUE,
      align = c("l", rep("c", ncol(.) - 1)))  %>%
  kable_classic() %>%
  add_header_above(
    .,
    c(
      "",
      "Correlations" = ncol(s3_MlCor_out),
      "Descriptives" = nrow(s3_MlCor_out) - ncol(s3_MlCor_out)
    )
  ) %>%
  pack_rows("Study 1", 1, ncol(s1_MlCor_out)) %>%
  pack_rows("Study 2", ncol(s1_MlCor_out) + 1,
            sum(ncol(s1_MlCor_out), ncol(s2_MlCor_out))) %>%
  pack_rows("Study 3", sum(ncol(s1_MlCor_out), ncol(s2_MlCor_out)) + 1,  
            sum(ncol(s1_MlCor_out), ncol(s2_MlCor_out), ncol(s3_MlCor_out))) %>%
  pack_rows("Across Studies", sum(ncol(s1_MlCor_out), ncol(s2_MlCor_out), ncol(s3_MlCor_out)) + 1,  
            sum(ncol(s1_MlCor_out), ncol(s2_MlCor_out), ncol(s3_MlCor_out), ncol(all_MlCor_out))) %>% # fix table position
  footnote(
    general = c(
      "Upper triangle: Within-person correlations;",
      "Lower triangle: Between-person correlations;",
      "*** p < .001, ** p < .01,  * p < .05"
    )
  ) %>%
  kable_styling(latex_options = "scale_down") %>%
  gsub("\\begin{table}", "\\begin{table}\n\\begin{minipage}[t][\\textheight][t]{\\textwidth}", ., fixed = TRUE) %>% # fix table position
  gsub("\\end{table}", "\\end{minipage}\n\\end{table}", ., fixed = TRUE) %>% # fix table position
  save_kable("Tables/descrOutWideAppB.tex")

```


```{r}
#| label: robustnessTblLong
#| include: false

# collect models for table
mdlRobustTblElements <- list(
  # Contact specific
  contact_overall = c(contact_overall$lmer_model, contact_overall$lmer_model_z, TRUE),
  contact_s1 = c(contact_s1$lmer_model, contact_s1$lmer_model_z, FALSE),
  contact_s2 = c(contact_s2$lmer_model, contact_s2$lmer_model_z, FALSE),
  contact_s3 = c(contact_s3$lmer_model, contact_s3$lmer_model_z, FALSE),
  # Accidental
  accidental_overall = c(accidental_overall$lmer_model, accidental_overall$lmer_model_z, TRUE),
  accidental_s1 = c(accidental_s1$lmer_model, accidental_s1$lmer_model_z, FALSE),
  accidental_s2 = c(accidental_s2$lmer_model, accidental_s2$lmer_model_z, FALSE),
  accidental_s3 = c(accidental_s3$lmer_model, accidental_s3$lmer_model_z, FALSE),
  # well-being core
  wellbeing_core_overall = c(wellbeing_core_overall$lmer_model, wellbeing_core_overall$lmer_model_z, TRUE),
  wellbeing_core_s1 = c(wellbeing_core_s1$lmer_model, wellbeing_core_s1$lmer_model_z, FALSE),
  wellbeing_core_s2 = c(wellbeing_core_s2$lmer_model, wellbeing_core_s2$lmer_model_z, FALSE),
  wellbeing_core_s3 = c(wellbeing_core_s3$lmer_model, wellbeing_core_s3$lmer_model_z, FALSE),
  # # POTENTIALLT well-being full
  # c(wellbeing_full_overall$lmer_model, wellbeing_full_s1$lmer_model_z),
  # c(wellbeing_full_s1$lmer_model, wellbeing_full_s1$lmer_model_z),
  # c(wellbeing_full_s2$lmer_model, wellbeing_full_s2$lmer_model_z),
  # c(wellbeing_full_s3$lmer_model, wellbeing_full_s3$lmer_model_z),
  # Need type
  need_type_overall = c(need_type_overall$lmer_model, need_type_overall$lmer_model_z, TRUE),
  need_type_s1 = c(need_type_s1$lmer_model, need_type_s1$lmer_model_z, FALSE),
  need_type_s2 = c(need_type_s2$lmer_model, need_type_s2$lmer_model_z, FALSE),
  need_type_s3 = c(need_type_s3$lmer_model, need_type_s3$lmer_model_z, FALSE),
  # Need SDT
  need_sdt_overall = c(need_sdt_overall$lmer_model, need_sdt_overall$lmer_model_z, TRUE),
  need_sdt_s1 = c(need_sdt_s1$lmer_model, need_sdt_s1$lmer_model_z, FALSE),
  need_sdt_s2 = c(need_sdt_s2$lmer_model, need_sdt_s2$lmer_model_z, FALSE),
  need_sdt_s3 = c(need_sdt_s3$lmer_model, need_sdt_s3$lmer_model_z, FALSE)
  )



mdlRobustPrep <- list()
for (i in 1:length(mdlRobustTblElements)) {
  mdlRobustPrep[[names(mdlRobustTblElements)[[i]]]] <- lmerTblPrep(obj = mdlRobustTblElements[[i]][[1]], 
                                                                   objz = mdlRobustTblElements[[i]][[2]], 
                                                                   lmer_scale = mdlRobustTblElements[[i]][[3]],
                                                                   name = names(mdlRobustTblElements)[[i]])
}

# extract important elements
shortRobustPrep <- list()
for(i in 1:length(mdlRobustPrep)){
  tmp <- rbind(
    mdlRobustPrep[[i]]$coeficients[, c("coef", "B", "Beta")],
    data.frame(
      coef = "R2Marg_R2Cond",
      B = mdlRobustPrep[[i]]$fit$R2Marg_R2Cond,
      Beta = ""
    )
  )
  shortRobustPrep[[mdlRobustPrep[[i]]$fit$name]] <- tmp
}

analyses <- unique(sub("_(?!.*_).*", "", names(mdlRobustPrep), perl = TRUE))
studies <- unique(sub(".*_(.*)", "\\1", names(mdlRobustPrep)))

robustTblLong <- list()
for (analysis in analyses) {
  mylist <- shortRobustPrep[grepl(analysis, names(shortRobustPrep))]
  mylist <- lapply(names(mylist), function(x) transform(mylist[[x]], name = x))
  df <- do.call(rbind, mylist)
  df <- melt(df, id.vars = c("coef", "name")) %>%
    mutate(
      name = as.character(sub(".*_(.*)", "\\1", name)),
      coef = as.character(sub("CoreNeedC", "CoreNeed", coef)),
      variable = as.character(variable)
    )
  df <- dcast(df, coef ~ name + variable, value.var = "value")
  robustTblLong[[analysis]] <- df
}
robustTblLong <- do.call(rbind, robustTblLong) %>% 
  rownames_to_column(var = "group") %>%
  mutate(group = factor(sub("\\.\\d+$", "", group), levels = unique(sub("\\.\\d+$", "", group))),
         contains_colon = ifelse(grepl(":", coef), 1, 0),
         special = case_when(
           coef == "(Intercept)" ~ "first",
           coef == "CoreNeed" ~ "second",
           coef == "R2Marg_R2Cond" ~ "last",
           TRUE ~ "middle"
         ),
         special = factor(special, levels = c("first", "second", "middle", "last"))) %>% 
  arrange(group, special, contains_colon, coef) %>% 
  select(-contains_colon, -special)

replacement_mapping <- c("contact" = "Contact specific",
                         "accidental" = "Interaction Intent",
                         "wellbeing_core" = "Well-being Outcome",
                         "need_type" = "Need Type",
                         "need_sdt" = "Specific Psychological Needs")

robustTblLong <- robustTblLong %>%
  mutate(
    group = case_when(
      group %in% names(replacement_mapping) ~ replacement_mapping[group],
      TRUE ~ as.character(group)
    ),
    coef = mapvalues(
      robustTblLong$coef,
      from = c("R2Marg_R2Cond", "practical_need_dum", "psychological_need_dum", "CoreNeed:practical_need_dum", "CoreNeed:psychological_need_dum"),
      to = c("$R^2_{marginal}$ / $R^2_{conditional}$", "practical need", "psychological need", "CoreNeed:practical need", "CoreNeed:psychological need")
    )
  )

mdlTblRobustLong <- 
  robustTblLong %>%
  select(-group) %>%
  replace(is.na(.), "") %>%
  kbl(
    .,
    caption = "Robustness Analyses",
    #format = "latex",
    linesep = "",
    booktabs = T,
    align = c("l", rep("c", ncol(.)-1)),
    escape = FALSE,
    col.names = c("", rep(c("$B$", "$\\beta$"), (ncol(.)-1)/2))
  ) %>%
  kable_classic() %>%
  add_header_above(., c(" ", "Overall" = 2, "Study 1" = 2, "Study 2" = 2, "Study 3" = 2)) %>%
  pack_rows(index = table(factor(robustTblLong$group, levels = unique(robustTblLong$group)))) %>%
  kableExtra::footnote(
    .,
    general = c("**** $p$ < .0001, *** $p$ < .001, ** $p$ < .01, * $p$ < .05"), 
    escape = FALSE
  ) %>%
  kable_styling(latex_options = "scale_down") %>%
  #kableExtra::landscape() %>%
  gsub("\\begin{table}", "\\begin{table}\n\\begin{minipage}[t][\\textheight][t]{\\textwidth}", ., fixed = TRUE) %>% # fix table position
  gsub("\\end{table}", "\\end{minipage}\n\\end{table}", ., fixed = TRUE) %>% # fix table position
  save_kable("Tables/mdlRobustlong.tex")

```


<!-- ```{r} -->
<!-- #| label: robustnessTblWide -->
<!-- #| include: false -->

<!-- # Wide format -->
<!-- mdlRobustTbl <- data.frame( -->
<!--   coef = c( -->
<!--     sapply(mdlRobustPrep, function(x) x[1]$coeficients$coef) %>% unlist %>% unique, -->
<!--     names(mdlRobustPrep[[1]]$random), -->
<!--     names(mdlRobustPrep[[1]]$fit) -->
<!--   ) -->
<!-- ) -->
<!-- for (i in 1:length(mdlRobustPrep)){ -->
<!--   mdlRobustTbl <- join(mdlRobustTbl, mdlRobustPrep[[i]]$mdlTbl, by = "coef") -->
<!--   #mdlTbl <- merge(mdlRobustTbl, mdlRobustPrep[[i]]$mdlTbl, by = "coef", all = TRUE, sort = FALSE, no.dups = TRUE, suffixes=c("", "")) -->
<!-- } -->

<!-- mdlRobustTblHeader <- data.frame( -->
<!--   data = as.character(mdlRobustTbl[mdlRobustTbl$coef == "data",]), -->
<!--   formula = as.character(mdlRobustTbl[mdlRobustTbl$coef == "formula",]) -->
<!-- ) %>% -->
<!--   filter(data != "", -->
<!--          data != "data") %>% -->
<!--   mutate(study = str_extract(data, "Worker|Student|Medical")) %>% -->
<!--   mutate(study = gsub("Worker", "Study 1", study)) %>% -->
<!--   mutate(study = gsub("Student", "Study 2", study)) %>% -->
<!--   mutate(study = gsub("Medical", "Study 3", study)) %>% -->
<!--   mutate(dv = sub(" ~.*", "", formula)) %>% -->
<!--   mutate(dv = gsub("Raw", " ", dv)) %>% -->
<!--   mutate(pred = str_extract(formula, "CoreNeed|Allport")) %>% -->
<!--   mutate(header = paste(study, pred, sep = " - ")) -->

<!-- mdlRobustTblHeaderFreq <- rbind.data.frame(rle(mdlRobustTblHeader$header)) %>% -->
<!--   mutate(lengths = lengths*2) %>% -->
<!--   rename(names = values, span = lengths) %>% -->
<!--   rbind(c(1, ""), .) %>% -->
<!--   mutate(span = as.numeric(span)) %>% -->
<!--   select(names, span) -->

<!-- mdlRobustTblHeaderDv <- data.frame(names = mdlRobustTblHeader$dv, span = 2) %>%  -->
<!--   mutate(names = as.character(names)) %>% -->
<!--   rbind(c("", 1), .) %>% -->
<!--   mutate(span = as.numeric(span)) -->

<!-- mdlRobustTblHeaderStudy <- rbind.data.frame(rle(mdlRobustTblHeader$study)) %>% -->
<!--   mutate(lengths = lengths*2) %>% -->
<!--   rename(names = values, span = lengths) %>% -->
<!--   rbind(c(1, ""), .) %>% -->
<!--   mutate(span = as.numeric(span)) %>% -->
<!--   select(names, span) -->

<!-- mdlRobustTbl$coef <- mapvalues( -->
<!--   mdlRobustTbl$coef, -->
<!--   from = c( -->
<!--     "sigmaSq", -->
<!--     "tau00", -->
<!--     "tau11", -->
<!--     "rho01", -->
<!--     "rho00", -->
<!--     "R2_marginal", -->
<!--     "R2_conditional",  -->
<!--     "R2Marg_R2Cond" -->
<!--   ), -->
<!--   to = c( -->
<!--     "$\\sigma^2$", -->
<!--     "$\\tau_{00}$", -->
<!--     "$\\tau_{11}$", -->
<!--     "$\\rho_{01}$", -->
<!--     "$\\rho_{00}$", -->
<!--     "$R^2_{marginal}$", -->
<!--     "$R^2_{conditional}$", -->
<!--     "$R^2_{marginal}$ / $R^2_{conditional}$" -->
<!--   ) -->
<!-- ) -->

<!-- # Long Table -->
<!-- mdlTblRobustLong <- mdlRobustTbl %>% -->
<!--   replace(is.na(.), "") %>% -->
<!--   kbl( -->
<!--     ., -->
<!--     caption = "Robustness Analyses", -->
<!--     #format = "latex", -->
<!--     linesep = "", -->
<!--     booktabs = T, -->
<!--     align = c("l", rep("c", ncol(.) - 1)), -->
<!--     escape = FALSE, -->
<!--     digits = 3 -->
<!--   ) %>% -->
<!--   kable_classic() %>% -->
<!--   add_header_above(., mdlRobustTblHeaderStudy) %>% -->
<!--   add_header_above( -->
<!--     ., -->
<!--     c( -->
<!--       "", -->
<!--       "Contact [Full Sample]" = 6, -->
<!--       "Self-Determination Theory [Outgroup Interactions only]" = 6 -->
<!--     ) -->
<!--   ) %>% -->
<!--   pack_rows("Fixed", 1, 7) %>% -->
<!--   pack_rows("Random", 8, 12) %>% -->
<!--   pack_rows("Fit", 13, 25) %>% -->
<!--   footnote( -->
<!--     general = c("PID = Participant ID number (grouping variable)", -->
<!--                 "$\\sigma^2$: residual variance, $\\tau_{00}$: random intercept variance, $\\tau_{11}$: random slope variance(s), $\\rho_{01}$: random slope-intercept correlation(s) , $\\rho_{00}$: correlation(s) between random slopes", -->
<!--                 "**** $p$ < .0001, *** $p$ < .001, ** $p$ < .01, * $p$ < .05" -->
<!--     ), -->
<!--     escape = FALSE -->
<!--   ) %>% -->
<!--   gsub("_cwc", "", ., fixed = TRUE)  -->

<!-- # shortened Table -->
<!-- mdlTblRobustShort <- mdlRobustTbl %>% -->
<!--   replace(is.na(.), "") -->
<!-- names(mdlTblRobustShort) <- make.unique(names(mdlTblRobustShort)) -->
<!-- mdlTblRobustShort <- mdlTblRobustShort %>% -->
<!--   filter( -->
<!--     !coef %in% c( -->
<!--       "data", -->
<!--       "formula", -->
<!--       "groupId", -->
<!--       "$R^2_{marginal}$", -->
<!--       "$R^2_{conditional}$", -->
<!--       "logLik", -->
<!--       "deviance" -->
<!--     ) -->
<!--   ) %>% -->
<!--   mutate_all(linebreak) %>% -->
<!--   kbl( -->
<!--     ., -->
<!--     caption = "Robustness Analyses", -->
<!--     format = "latex", -->
<!--     linesep = "", -->
<!--     booktabs = TRUE, -->
<!--     col.names = c("", rep(c("$B$", "$\\beta$"), (ncol(.)-1)/2)), -->
<!--     align = c("l", rep("c", ncol(.) - 1)), -->
<!--     escape = FALSE, -->
<!--     digits = 3 -->
<!--   )  %>% -->
<!--   kable_classic() %>% -->
<!--   add_header_above(., mdlRobustTblHeaderStudy) %>% -->
<!--   add_header_above( -->
<!--     ., -->
<!--     c( -->
<!--       "", -->
<!--       "Contact [Full Sample]" = 6, -->
<!--       "Self-Determination Theory [Outgroup Interactions only]" = 6 -->
<!--     ) -->
<!--   ) %>% -->
<!--   pack_rows("Fixed", 1, 7) %>% -->
<!--   pack_rows("Random", 8, 12) %>% -->
<!--   pack_rows("Fit", 13, 18) %>% -->
<!--   footnote( -->
<!--     general = c("PID = Participant ID number (grouping variable)", -->
<!--                 "$\\sigma^2$: residual variance, $\\tau_{00}$: random intercept variance, $\\tau_{11}$: random slope variance(s), $\\rho_{01}$: random slope-intercept correlation(s) , $\\rho_{00}$: correlation(s) between random slopes", -->
<!--                 "**** $p$ < .0001, *** $p$ < .001, ** $p$ < .01, * $p$ < .05" -->
<!--     ),  -->
<!--     escape = FALSE -->
<!--   ) %>% -->
<!--   kable_styling(latex_options = "scale_down") %>% -->
<!--   kableExtra::landscape() %>% -->
<!--   gsub("_cwc", "", ., fixed = TRUE)  %>% -->
<!--   gsub("\\makecell[l]", "\\makecell[c]", ., fixed = TRUE) %>% -->
<!--   gsub("\\makecell[r]", "\\makecell[c]", ., fixed = TRUE) %>% -->
<!--   gsub("objVar$", "", ., fixed = TRUE) %>% -->
<!--   gsub("slope_intercept.", "slope.intercept", ., fixed = TRUE) %>% -->
<!--   gsub("\\begin{table}", "\\begin{table}\n\\begin{minipage}[t][\\textheight][t]{\\textwidth}", ., fixed = TRUE) %>% # fix table position -->
<!--   gsub("\\end{table}", "\\end{minipage}\n\\end{table}", ., fixed = TRUE) %>% # fix table position -->
<!--   save_kable("Tables/mdlRobustness.tex") -->
<!-- mdlTblRobustShort -->
<!-- ``` -->

<!-- FIGURES: -->

```{r}
#| label: metaRobust
#| include: false

analyses <- unique(sub("_(?!.*_).*", "", names(mdlRobustTblElements), perl = TRUE))
studies <- unique(sub(".*_(.*)", "\\1", names(mdlRobustTblElements)))

# i = 1
robustMetaTbl <- data.frame()
for (i in 1:length(mdlRobustTblElements)){
  Prep <- lmerTblPrep(mdlRobustTblElements[[i]][[1]], mdlRobustTblElements[[i]][[2]])
  lmerTblTmp <- Prep$beta %>%
  transmute(
    group = names(mdlRobustTblElements)[[i]],
    analysis = sub("_(?!.*_).*", "", names(mdlRobustTblElements)[[i]], perl = TRUE),
    coef = gsub(":", ".", coef),
    B = est,
    se,
    lwr,
    upr,
    tval,
    Rsq,
    n = Prep$fit$N,
    m = nrow(.), 
    #dv = gsub("Dutch", "", str_extract(lmerTblPrep(mdlRobustTblElements[[i]][[1]], mdlRobustTblElements[[i]][[2]])$fit$formula,".+?(?= ~)")),
    sample = sub(".*_(.*)", "\\1", names(mdlRobustTblElements)[[i]]),
    formula = gsub("\\s+", " ", Prep$fit$formula)
  )
  robustMetaTbl <- rbind(robustMetaTbl, lmerTblTmp)
}
robustMetaTblSorted <- robustMetaTbl %>% 
  mutate(group = factor(group, levels = unique(group)),
         contains_period = ifelse(grepl("\\.", coef), 1, 0),
         special = case_when(
           coef == "(Intercept)" ~ "first",
           coef == "CoreNeed" ~ "second",
           coef == "R2Marg_R2Cond" ~ "last",
           TRUE ~ "middle"
         ),
         special = factor(special, levels = c("first", "second", "middle", "last"))) %>% 
  arrange(group, special, contains_period, coef) %>% 
  select(-contains_period, -special) %>%
  filter(sample != "overall") %>% 
  mutate(formula = gsub(" \\+ \\(.*", "", formula))

replacement_mapping <- c("contact" = "Contact Specific",
                         "accidental" = "Interaction Intent",
                         "wellbeing_core" = "Well-being Outcome",
                         "need_type" = "Need Type",
                         "need_sdt" = "Specific Psychological Needs")

robustMeta <- list()
for (name in unique(robustMetaTblSorted$analysis)) {
  lab <- case_when(
    name %in% names(replacement_mapping) ~ replacement_mapping[name],
    TRUE ~ as.character(name)
  )
  robustMeta[[name]] <- metaLmerOut(
    lmerDataTbl = robustMetaTblSorted %>% filter(analysis == name),
    type = "REML", #"FE", 
    name = name,
    title = paste0("Robustness — ", lab), 
  )
}
```

```{r}
#| label: robustForestPltComb
#| include: false

process_element <- function(element) {
  # Apply dplyr filtering and mutating
  df_eff_adjusted <- element$eff %>%
    filter(coef != "(Intercept)") %>%
    mutate(metaId = paste(coef, analysis, sep = ".")) %>%
    as.data.frame
  
  return(df_eff_adjusted)
}
adjusted_dataframes <- lapply(robustMeta, process_element)

lookup_table <- data.frame(analysis = unique(names(adjusted_dataframes)), group_letter = LETTERS[seq_along(unique(names(adjusted_dataframes)))])

effRobust <- do.call(rbind, adjusted_dataframes) %>%
  left_join(lookup_table, by = "analysis") %>%
  mutate(
    analysis = paste0("(",group_letter, ") ",
      case_when(
        analysis %in% names(replacement_mapping) ~ replacement_mapping[analysis],
        TRUE ~ as.character(analysis)
      )
    ),
    contains_period = ifelse(grepl("\\.", coef), 1, 0),
    special = case_when(coef == "CoreNeedZ" ~ "first", 
                        TRUE ~ "last"),
    special = factor(special, levels = c("first", "last"))
  ) %>%
  arrange(analysis, special, contains_period, coef, sample) %>%
  select(-group_letter, -special, -contains_period)


metaSubForest(effects = effRobust, 
              title = "Meta Analysis: Forest Plot Robustness Analysis \n[Parametric]",
              filename = "RobustnessComb",
              addAbove = 5,
              width = 600,
              height = 1200)
```