Supplemental-Material-B-Power-Simulation_between.Rmd

---
title: "Supplemental Information C: Power Simulations"
subtitle: "Supplemental Material for 'Need Fulfillment During Intergroup Contact: Three Experience Sampling Studies'"
author:
- ██████ ██████████^1^
- ██████████ ████████^1^
- █████ █████████^1^
- █████ ██ █████^1^
- ███ ███████^1^
- ^1^██████████ ██ █████████, Department of Psychology
- "Author Information:"
- "Correspondence concerning this article should be addressed to ██████ ██████████, ██████████ ██ ██████████, █████████ ██ ██████████, █████ ██████████ ███, ████ ██ █████████ ████████████████. E-mail: █████████████████"
- 'The main manuscript is available at <a href="https://www.doi.org/ToBePublished" target="_blank">doi.org/ToBePublished</a>'
- 'The data repository for this manuscript is available at <a href="https://osf.io/pr9zs/?view_only=1ea47bb646694632a764dead807ef970" target="_blank">osf.io/pr9zs/?view_only=1ea47bb646694632a764dead807ef970</a>'
- 'The GitHub repository for this manuscript is available at <a href="https://github.com/maskedForPeerReview" target="_blank">github.com/maskedForPeerReview</a>'
date: "Last updated: `r format(Sys.time(), '%d %B, %Y')`"
output: 
  bookdown::html_document2:
    fig_caption: yes
    md_extensions: +footnotes
    code_folding: hide
    mathjax: default
    theme: yeti
    toc: yes
    toc_float: yes
    number_sections: false
    css: style.css
    includes:
      in_header: "_includes/head-custom-rmd.html" 
editor_options:
  chunk_output_type: console
bibliography: references.bib
csl: apa.csl
header-includes:
   - \usepackage{amsmath, nccmath}
---

<style type="text/css">
.main-container {
  max-width: 1300px;
  margin-left: auto;
  margin-right: auto;
}
.table {
  margin-left:auto; 
  margin-right:auto;
}
</style>


```{r 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",
  "remedy",
  "bookdown",
  "brms",
  "psych",
  "ggplot2",
  "ggthemes",
  # "haven",
  "RColorBrewer",
  # "plotly",
  "gridExtra",
  "ggpattern",
  "ggridges",
  "binom",
  "iterators",
  "pbkrtest",
  "plotrix",
  "RLRsim",
  "stats",
  "methods",
  "utils",
  "graphics",
  "grDevices",
  "car",
  "testthat",
  "simr",
  "data.table",
  "lme4",
  "lmerTest",
  # "nlme",
  # "jtools",
  # "gtsummary",
  "sessioninfo",
  # "tibble",
  "pander",
  # "devtools",
  # "mada",
  #"tidyr",
  #"tidyverse",
  "plyr",
  "dplyr",
  # "Hmisc",
  "kableExtra",
  # "papaja",
  "stringr"#,
  # "stringi",
  # "reshape2",
  # "lubridate",
  # "purrr",
  # "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
knitr::knit_hooks$set(
  error = function(x, options) {
    paste('\n\n<div class="alert alert-danger">',
      gsub("##", "\n", gsub("^##\ Error", "**Error**", x)),
      "</div>",
      sep = "\n"
    )
  },
  warning = function(x, options) {
    paste('\n\n<div class="alert alert-warning">',
      gsub("##", "\n", gsub("^##\ Warning:", "**Warning**", x)),
      "</div>",
      sep = "\n"
    )
  },
  message = function(x, options) {
    paste('\n\n<div class="alert alert-info">',
      gsub("##", "\n", x),
      "</div>",
      sep = "\n"
    )
  }
)
htmltools::tagList(rmarkdown::html_dependency_font_awesome())

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

# set ggplot theme
theme_set(theme_Publication())
```

<br/>

<i class="fas fa-exclamation-circle"></i> Note. Boxplots display the interquartile range (IQR, center box), and the whiskers extend 1.5*IQR from the lower and upper hinge. The white point indicates the mean and the white center line indicates the median.   

<br/>

Generally, power considerations of mixed effects model such as repeated measures ANCOVAs and multi-level regressions are difficult to estimate because of the complex covariance structures. We ran simulation-studies based on our initial results following the first study. We found that sufficient sample sizes heavily depended on the type and size of the effect. When we simulated models with the marginalization sub-scale (which was the best distributed with mostly small to medium effect sizes), we found that for between-subjects effects, a sample size of 80 participants would have been sufficient (power above .8) with two measurements (see Figure 9 for an illustration). However, for the same between-subject effects we did not reach a power above .8 with any increase of repeated follow-up measurements (based on model extrapolations with a fixed _n_ = 20; tested up to 12 months). The within time effect, on the other hand, was estimated well above .8 by the seventh measurement (i.e., after six months, see Figure 9). The cross-level interaction effects (within-between) were largely a mixture of their individual effects and changes in measurement numbers and sample size affected the models accordingly. We, therefore, urge researchers to rely on pilot data and to run extensive simulation studies to maximize effort and feasibility (i.e., sufficient power).

# Data Import  

```{r}
#| label: import data

load("data/S1_Workers/processed/df.btw.RData")
```

# Extract Effects for Simulation

## Prepare Data

```{r}
#| label: data prep

# prepare data for regression (i.e, put into long format and change variable types)
marginalization = melt(df.btw[!is.na(df.btw$assimilation.post),c("marginalization.pre","marginalization.post",
                                                                 "CtContactNL_c","AvKeyNeedInt_c",
                                                                 "ExternalReference")],
                       id=c("ExternalReference","CtContactNL_c","AvKeyNeedInt_c"))
marginalization$time = ifelse(grepl(".pre", marginalization$variable), "pre", "post")
marginalization$variable = gsub(".pre$|.post$", "", marginalization$variable)
marginalization$time = factor(marginalization$time)
marginalization$time = factor(marginalization$time,levels(marginalization$time)[c(2,1)])
marginalization$id = rep(seq(1:20), 2)
```


# ML Regression

Run multilevel regression to extract the effect sizes and covariance matrices

```{r}
#| label: run regression 

# run regresion 
lm.marg <- lmer(value~
                  AvKeyNeedInt_c*time+
                  CtContactNL_c*time+
                  (1|id), 
                data=marginalization)
summary(lm.marg)
fixef(lm.marg)
```


# Post-Hoc Power

Check post-hoc power for comparison.

```{r}
#| label: post hoc simulation

# Post-hoc power analysis for single effects (example: key need)
post_hoc <- powerSim(lm.marg, fixed("AvKeyNeed"), nsim=1000)
post_hoc
```

# Power Simulation Participant Numbers

Simulate different participant numbers for all fixed effects.

## Power of contact effect for different numbers of participants

```{r}
#| label: participant simulation contact

# extrapolate data:
sim.N <- extend(lm.marg, along="id", n=500)

# Run and plot simulations
plt.pwr.contact = powerCurve(fit = sim.N, 
                             test = simr::fixed("CtContactNL_c", "t"), 
                             along="id", 
                             progress = FALSE,
                             breaks = seq(10,100,10), 
                             nsim=1000)
plot(plt.pwr.contact) +
  title("Outgroup Contact Count [participant numbers]")
```

## Power of time effect for different numbers of participants

```{r}
#| label: participant simulation time

plt.pwr.time = powerCurve(sim.N, simr::fixed("time"), 
                          along="id", 
                          progress = FALSE, 
                          breaks = seq(10,100,10), 
                          nsim=1000)
plot(plt.pwr.time) +
  title("Time [participant numbers]")
```

## Power of need fulfillment effect for different numbers of participants

```{r}
#| label: participant simulation avg need fulfillment

plt.pwr.Need = powerCurve(sim.N, simr::fixed("AvKeyNeedInt_c", "t"), 
                          along="id", 
                          progress = FALSE, 
                          breaks = seq(10,100,10), 
                          nsim=1000)
plot(plt.pwr.Need) +
  title("Avergage Need Fulfillment [participant numbers]")
```

## Power of time by contact effect for different numbers of participants

```{r}
#| label: participant simulation contact by time

plt.pwr.timeCont = powerCurve(sim.N, simr::fixed("time:CtContactNL_c"), 
                             along="id", 
                             progress = FALSE, 
                             breaks = seq(10,500,50), 
                             nsim=1000)
plot(plt.pwr.timeCont) +
  title("Time X Outgroup Contact Count [participant numbers]")
```

## Power of time by need fulfillment effect for different numbers of participants

```{r}
#| label: participant simulation need fulfillment by time

plt.pwr.timeNeed = powerCurve(sim.N, simr::fixed("AvKeyNeedInt_c:time"), 
                              along="id", 
                              progress = FALSE, 
                              breaks = seq(10,50,5), 
                              nsim=1000)
plot(plt.pwr.timeNeed) +
  title("Time X Avergage Need Fulfillment [participant numbers]")
```


# Power Simulation Measurement Numbers

Simulate different number of within person measurements.

## Power of time effect for different numbers of measurements

```{r}
#| label: measurement simulation time

# extrapolate data:
sim.within <- extend(lm.marg, within="id", n=12)

# Run and plot simulations
plt.pwr.t.time = powerCurve(sim.within, simr::fixed("time"), 
                          within="id", 
                          progress = FALSE, 
                          breaks = seq(2,12,1), 
                          nsim=1000)
plot(plt.pwr.t.time) +
  title("Time [measurement numbers]")
```

## Power of contact effect for different numbers of measurements

```{r}
#| label: measurement simulation contact count

plt.pwr.t.contact = powerCurve(sim.within, simr::fixed("CtContactNL_c", "t"), 
                               within="id", 
                               progress = FALSE, 
                               breaks = seq(2,12,1), 
                               nsim=1000)
plot(plt.pwr.t.contact) +
  title("Outgroup Contact Count [measurement numbers]")
```

## Power of need fulfillment effect for different numbers of measurements

```{r}
#| label: measurement simulation contact need fulfillment

plt.pwr.t.Need = powerCurve(sim.within, simr::fixed("AvKeyNeedInt_c", "t"), 
                            within="id", 
                            progress = FALSE, 
                            breaks = seq(2,12,1), 
                            nsim=1000)
plot(plt.pwr.t.Need) +
  title("Average Need Fulfillment [measurement numbers]")
```

## Power of time by contact effect for different numbers of measurements

```{r}
#| label: measurement simulation contact time by contact count

plt.pwr.t.timeCont = powerCurve(sim.within, simr::fixed("time:CtContactNL_c"), 
                          within="id", 
                          progress = FALSE, 
                          breaks = seq(2,12,1), 
                          nsim=1000)
plot(plt.pwr.t.timeCont) +
  title("Time X Outgroup Contact Count [measurement numbers]")
```

## Power of time by need fulfillment effect for different numbers of measurements

```{r}
#| label: measurement simulation need fulfillment by time

plt.pwr.t.timeNeed = powerCurve(sim.within, simr::fixed("AvKeyNeedInt_c:time"), 
                               within="id", 
                               progress = FALSE, 
                               breaks = seq(2,12,1), 
                               nsim=1000)
plot(plt.pwr.t.timeNeed) +
  title("Time X Average Need Fulfillment [measurement numbers]")
```

# Export data and plots

```{r}
#| label: export

# save simulations, because this takes for freakin ever.      
save(plt.pwr.contact,plt.pwr.Need,plt.pwr.t.contact,plt.pwr.t.Need,
     plt.pwr.t.time,plt.pwr.t.timeCont,plt.pwr.t.timeNeed,plt.pwr.time,
     plt.pwr.timeCont,plt.pwr.timeNeed, 
     file = "data/S1_Workers/processed/PowerSimPlots.RData")

#load("PowerSimPlots.RData")
```

# Software Information  

The full session information with all relevant system information and all loaded and installed packages is available in the collapsible section below.  

<details>
  <summary>System Info</summary>
  
  \renewcommand{\arraystretch}{0.8} <!-- decrease line spacing for the table -->
```{r Reproducibility-SessionInfo-R-environment, echo=FALSE, message=FALSE, warning=FALSE, fig.align="center", out.width='100%', results='asis'}
df_session_platform <- devtools::session_info()$platform %>%
  unlist(.) %>%
  as.data.frame(.) %>%
  rownames_to_column(.)

colnames(df_session_platform) <- c("Setting", "Value")

kbl(
  df_session_platform,
  booktabs = T,
  align = "l",
  caption = "R environment session info for reproducibility of results" # complete caption for main document
) %>%
  kable_classic(
    full_width = F,
    lightable_options = "hover",
    html_font = "Cambria"
  )
```
  \renewcommand{\arraystretch}{1} <!-- reset row height/line spacing -->
 </details>
 <br>
 <details>
  <summary>Package Info</summary>

\renewcommand{\arraystretch}{0.6} <!-- decrease line spacing for the table -->
```{r Reproducibility-SessionInfo-R-packages, echo=FALSE, message=FALSE, warning=FALSE, fig.align="center", out.width='100%', results='asis'}
df_session_packages <- devtools::session_info()$packages %>%
  as.data.frame(.) %>%
  filter(attached == TRUE) %>%
  dplyr::select(loadedversion, date, source) %>%
  rownames_to_column()

colnames(df_session_packages) <- c("Package", "Loaded version", "Date", "Source")

kbl(
  df_session_packages,
  booktabs = T,
  align = "l",
  caption = "Package info for reproducibility of results" # complete caption for main document
) %>%
  kable_classic(
    full_width = F,
    lightable_options = "hover",
    html_font = "Cambria"
  )
```
\renewcommand{\arraystretch}{1} <!-- reset row height/line spacing -->
</details>
<br>
<details>
  <summary>Full Session Info (including loaded but unattached packages --- for troubleshooting only)</summary>
    `r pander(sessionInfo(), compact = FALSE)`
</details>


--------------------------------------------------------------------

</br>