data_import.qmd

---
title: "Data import & descriptives"
title-block-banner: "#870052"
title-block-banner-color: "#FFFFFF"
author: 
  name: Magnus Johansson
  affiliation: RISE Research Institutes of Sweden
  affiliation-url: https://ri.se/shic
  orcid: 0000-0003-1669-592X
date: 2023-06-12
date-format: YYYY-MM-DD
format: 
  html:
    toc: true
    embed-resources: true
    standalone: true
    mainfont: 'Arial'
    code-overflow: wrap
    code-fold: show
    code-link: true
    number-sections: true
    fig-dpi: 96
    layout-align: left
    linestretch: 1.6
    theme: materia
    
  pdf:
    papersize: a4
    documentclass: article #article, report or book
    #classoption: [twocolumn, portrait]
  docx:
    toc: false
    number-sections: false
    highlight-style: github
    
always_allow_html: true
execute:
  echo: true
  warning: false
  message: false
  cache: false
css: style.css
editor_options: 
  chunk_output_type: console
---

## Overview of this file

- importing data
- looking at data
- descriptives
  - tables
  - figures
  - basic statistical analyses
- missing data
- checking for outliers
- wrangling data (wide to long, etc)


## Data cleaning

![Image credit: <https://kbroman.org/Talk_DataCleaning/data_cleaning.pdf>](images/RulesDataCleaning.jpg){fig-align="center"}

## Setting up for data analysis {background-color="#99D9DD"}

Let's load packages/libraries.

```{r}
#| code-fold: show
# these are mostly for data management/wrangling and visualization
library(tidyverse) # for most things
library(foreign) # for reading SPSS files
library(readxl) # read MS Excel files
library(showtext) # get fonts
library(glue) # simplifies mixing text and code in figures and tables
library(arrow) # support for efficient file formats
library(grateful) # create table+references for packages used in a project
library(styler) # only a one-time installation (it is an Rstudio plugin)
library(car) # for car::recode only
library(skimr) # data skimming
library(lubridate) # for handling dates in data
library(janitor) # for many things in data cleaning

# these are mostly for data analysis and visualization
library(gtsummary)
library(scales)
library(visdat)
library(psych)
library(lme4)
library(nlme)
library(broom.mixed)
library(patchwork)
library(easystats)
library(mice)
library(modelsummary)
library(ggrain)
library(ggdist)
library(kableExtra)
library(formattable)
library(ggrepel)
library(GGally)

```

Define a ggplot theme `theme_ki()`, a standard table function, `kbl_ki()`, and a color palette based on KI's design guide, `ki_color_palette`.

```{r}
#| code-fold: show

source("ki.R") # this reads an external file and loads whatever is in it
```

### Adaptions {.smaller}

Some functions exist in multiple packages, which can be a source of headaches and confusion. Loading `library(conflicted)` will provide errors every time you use a function that is available in multiple loaded packages, which can be helpful to avoid problems (but also annoying if you already have things under control).

Below we define preferred functions that are frequently used. If desired, we can still use specific functions by using their package prefix, for instance `dplyr::recode()`.

```{r}
#| code-fold: show

#library(conflicted)
select <- dplyr::select
count <- dplyr::count
recode <- car::recode
rename <- dplyr::rename
filter <- dplyr::filter
clean_names <- janitor::clean_names
```

## Importing data

The open dataset we will use for our experiments was retrieved from <https://doi.org/10.26180/13240304> and is available in the `data` subfolder of the R project folder we are currently working in. The description of the dataset on Figshare is:

>De-identified dataset from a randomised controlled trial of Mindfulness-integrated cognitive behaviour therapy (MiCBT) versus a treatment-as-usual waitlist control. All participants completed the measures one week before the start of the MiCBT group intervention (T0), after week 4 (T1), at week 8 (T2, post-intervention), and then again after a 6-month follow up period (T3). A full description of the project methodology including the measures used in the trial is provided in the protocol paper (see References).

And from the study protocol:

> The intent of this study is to examine the effectiveness of MiCBT to create changes in clinical measures of depression, anxiety and stress. It is hypothesized that these changes will occur during the program in stages 1,2 and 3 and be enhanced in stage 4 because of the additional practice time. Compassion and ethics are taught in Stage 4 for relapse prevention which is not the focus of the current study.

Looking at [the abbreviations section of the study protocol](https://bmcpsychiatry.biomedcentral.com/articles/10.1186/s12888-019-2411-1#abbreviations), we can hopefully get some variable name explanations:

- EQ: Experiences Questionnaire
- FS: Flourishing scale
- K10: Kessler Psychological Distress Scale
- MAIS: Multidimensional Assessment of Interoceptive Awareness
- MB-EAT: Mindfulness-based Eating Program
- MBI: Mindfulness-based Intervention
- MBRE: Mindfulness-based Relationship Enhancement
- MBSR: Mindfulness-based Stress Reducyion
- MiCBT: Mindfulness integrated Cognitive Behavior Therapy
- MSES: Mindfulness-based Self-efficacy Scale
- NAS: Non-attachment Scale
- SWLS: Satisfaction with Life Scale

Let's read the datafile and have a first look.

```{r}
df <- read_excel("data/MiCBT RCT data_Bridges repository.xlsx")
```
Look in the Environment quadrant (upper right). How many observations and variables do we have in the `df` object?

Press the circle to the left of `df` to get a quick look at the data. We can see the word "missing" noted in several fields. Anything else you notice about the variables?

Let's re-import the data and tell `read_excel()` to code missing correctly.

```{r}
df <- read_excel("data/MiCBT RCT data_Bridges repository.xlsx",
                 na = "missing")
```

Have another look at the data now and see what happened. You can go back and run the previous chunk to see the difference more clearly.

Also, have a look at the naming scheme and see what pattern you find?

The K10 questionnaire is used for pre-intervention measurement and screening, as well as follow-up measurement. Let's look at the variables containing "K10".

```{r}
df %>% 
  select(contains("K10"))
```

`K10_di_t0` is a categorical variable created from `TK10_t0`, and it does not repeat for other time points. As such, it is mislabeled and we want to fix this. While we are at it, we can rename some other variables too.

The syntax for `dplyr::rename()` is `newname = oldname`.

```{r}
df <- df %>% 
  rename(id = BridgesID,
         Group = GROUP,
         K10preCat = K10_di_t0)
```

## Demographics

The dataset does not include any demographics. Just for fun, we'll add randomly assigned age and gender variables. `mutate()` helps us create or modify variables.

```{r}
set.seed(1234)

df <- df %>% 
  mutate(age = rnorm(nrow(df), 
                     mean = 44, 
                     sd = 8),
         age = as.integer(age),
         age = recode(age,"0:18=19"),
         sex = sample(1:2, nrow(df), replace=TRUE))
```

We made a gender/sex variable that is numeric, since this is often the case. You may want to turn it into a factor instead, with labels.

```{r}

df <- df %>% 
  mutate(sex = factor(sex, 
                 levels = c(1,2),
                 labels = c("Female","Male")))
```


`summary()` is a general function that can be used with many types of objects, including model outputs.
```{r}
summary(df$age)
```


### Tables

::: panel-tabset
#### All participants

```{r}
df %>% 
  select(age) %>% 
  tbl_summary()
```

#### Control/intervention group
```{r}
df %>% 
  select(age,Group) %>% 
  tbl_summary(by = Group,
              statistic = list(all_continuous() ~ "{mean} ({sd}, {min}-{max})"))
```
#### Difference?
We can see that there isn't a difference, and we expect this since random sampling was used. But if you wanted to test the difference, this is one way.

```{r}
age.ttest <- t.test(age ~ Group, data = df)
age.ttest
```

```{r}
tidy(age.ttest) %>% 
  mutate_if(is.double, round, digits = 2) %>% 
  kbl_ki()
```

:::



### Figures

Base R examples.
```{r}
hist(df$age)

hist(df$age, 
     col = "lightblue", 
     main = "Histogram of participant age",
     xlab = "Age",
     breaks = 24)
```

With ggplot we have a lot more flexibility. Note that as soon as `ggplot()` has been called, the line ends with `+` when we add plot configurations.
```{r}
df %>% 
  ggplot(aes(x = age)) +
  geom_histogram(fill = "lightblue",
                 color = "black") +
  labs(title = "Histogram of participant age",
       x = "Age",
       y = "Count") +
  theme_ki()
```

Let's look separately at the Control group.

```{r}
df %>% 
  filter(Group == "Control") %>% 
  ggplot(aes(x = age)) +
  geom_histogram(fill = "darkgreen",
                 color = "white") +
  labs(title = "Histogram of participant age",
       x = "Age",
       y = "Count",
       subtitle = "Control group only") +
  theme_ki()
```

::: {.callout-tip}
#### Practice 

- Make a new plot for intervention group age.
:::


### tidy filter/select

::: {.callout-note}
Two key things to learn:

- `filter()` works on rows, based on their column content
- `select()` works on columns, based on their names
:::

#### Examples

This filters for inclusion of only one municipality:
``` r
df %>% 
  filter(municipality == "Vallentuna")
```

This filters out (excludes) multiple municipalities:
``` r
df %>% 
  filter(!municipality %in% c("Vallentuna","Vaxholm"))
```

#### Fill/color based on Group variable.

We can have both in the same plot.

- dynamic color/fill (based on a data variable) needs to be defined within the `aes()` function (aesthetics).

```{r}
df %>% 
  ggplot(aes(x = age, fill = Group)) +
  geom_histogram(color = "black") +
  labs(title = "Histogram of participant age",
       x = "Age",
       y = "Count") +
  theme_ki() +
  scale_y_continuous(breaks = c(0,4,8,12)) +
  scale_color_manual(values = ki_color_palette, 
                     aesthetics = c("fill","color"))

```

Or we can use `facet_wrap()` to make parallel plots.
```{r}
#| fig-height: 4
df %>% 
  ggplot(aes(x = age, fill = Group)) +
  geom_histogram(color = "white",
                 binwidth = 3) +
  labs(title = "Histogram of participant age",
       x = "Age",
       y = "Count") +
  theme_ki() +
  scale_y_continuous(breaks = c(0,4,8,12)) +
  scale_color_manual(values = ki_color_palette, 
                     aesthetics = c("fill","color"),
                     guide = "none") +
  facet_wrap(~Group,
             ncol = 2)
```

::: {.callout-tip}
#### Practice 

- Plot age grouped by sex 
  - bonus task: try facet_grid() to plot age by sex and Group
:::

## Variable names

Generally we should have systematic naming of variables, avoiding things like spaces (" "). There is an amazing function called `janitor::clean_names()` which defaults to using snake_case. It also offers options for things like camelCase and others. This functions is primarily useful when you get a dataset that someone else collected and you need to bring order to variables names.

```{r}
df %>% 
  select(ends_with("t0")) %>% 
  names()
```

Why is there a "T" at the beginning of most variables? What happens if we remove it?

```{r}
df %>% 
  select(ends_with("t0")) %>% 
  rename_all(~ str_replace(.x, "^T","")) %>% 
  names()
```

A lot more readable. Please note that we did not "save" our changes in the previous code chunk. Let's rename all variables in the dataframe that start with a capital "T".

```{r}
df <- df %>% 
  rename_all(~ str_replace(.x, "^T",""))
```

## Skim

We can use `skim()` to get a more detailed overview of the variables.

```{r}
df %>% 
  select(starts_with(c("ANX","DEP","STRESS"))) %>% 
  skim()
```

And grouped:

```{r}
df %>% 
  select(Group,starts_with(c("ANX","DEP","STRESS"))) %>% 
  group_by(Group) %>% 
  skim() %>% 
  arrange(Group) #try other sorting variables, and also try to reverse sorting
```


## Correlation & visualization

### Correlation matrix


```{r}
df %>% 
  select(ANX_t0,DEP_t0,STRESS_t0) %>% 
  ggpairs()
```

### Correlation grouped

```{r}
df %>% 
  select(sex,ANX_t0,DEP_t0,STRESS_t0) %>% 
  ggpairs(aes(color = sex, alpha = 0.85)) +
  scale_color_manual(values = ki_color_palette,
                     aesthetics = c("color", "fill")) +
  theme_ki()

```

### easystats

::: {.callout-note}
From here, we will practice to "break down" the pieces in the tidy/ggplot code by selecting pieces of it and running it, adding one row/function/layour at a time.
:::

```{r}
df %>% 
  cor_test("ANX_t0", "DEP_t0") %>% 
  plot()
```

With some styling/theming.
```{r}
df %>% 
  cor_test("ANX_t0", "DEP_t0") %>% 
  plot() +
  theme_ki() +
  geom_point(data = df, 
             aes(ANX_t0, DEP_t0), 
             size = 2.4,
             color = "#870052") +
  geom_smooth(data = df, 
              aes(ANX_t0, DEP_t0),
              method = "lm",
              fill = "#FF876F",
              color = "#4F0433",
              alpha = 0.4) +
  labs(y = "Depression at pre (t0)",
       x = "Anxiety at pre (t0)",
       title = "Correlation between DASS-D and DASS-A at time 0.")
```

**Exercise:** create separate correlation plots for gender.
Bonus points if you can get both in the same plot!

## Missing data
```{r}
df %>% 
  select(starts_with(c("ANX","DEP","STRESS"))) %>% 
  vis_dat(palette = "qual")
```

### Imputation etc?

Please see <https://vincentarelbundock.github.io/marginaleffects/articles/multiple_imputation.html> for some options.

And notes on missing data and LMM's: <https://rpsychologist.com/lmm-slope-missingness>

## Outliers

```{r}
df %>% 
  select(id, starts_with(c("ANX","DEP","STRESS"))) %>% 
  na.omit() %>% 
  check_outliers(ID = "id")
```

Should we try other methods? See `?check_outliers`.


## Wide to long format

Almost everything in R likes long format. Let's look at the variables ending with "t".

```{r}
df %>% 
  select(ends_with("t0")) %>% 
  names()

df %>% 
  select(ends_with("t1")) %>% 
  names()
```

For simplicity, we'll initially focus on the main outcome variables, DEP, ANX, and STRESS. Let's pivot the data to long format.

```{r}
df.long <- df %>% 
  select(id,Group,starts_with(c("ANX","DEP","STRESS"))) %>% 
  pivot_longer(starts_with(c("ANX","DEP","STRESS")),
               names_to = c("measure", "time"),
               names_sep = "_t")
```

And look at the df.

```{r}
glimpse(df.long)
```

## Done!

We'll write `df.long` to file for use in analysis (soon).

```{r}
write_csv(df.long,"data/dataLong.csv")
```

::: {.callout-tip}
If you have a large dataset, I highly recommend using `library(arrow)` and the function `write_parquet()`, since it is incredibly fast and produces a small file. As an example, a 450mb SPSS datafile became 8mb when saved in .parquet format.
:::

## Bonus - dealing with questionnaire data

::: {.callout-note}
This section will mostly be skimmed and the various plots and solutions are included for you to use as a reference if/when you would like to implement something similar in your future analyses. We can also go back to this section if we have time at the end.
:::

We'll use a dataset that actually includes raw response data from SurveyMonkey.

```{r}
df2 <- read.spss("data/2023-04-26 Prevent OSA-enkat.sav", to.data.frame = T) %>% 
  select(starts_with("q0010"))

itemlabels <- read_excel("data/Itemlabels.xlsx") %>% 
  filter(str_detect(itemnr, pattern = "ls")) %>% 
  select(!Dimension)

names(df2) <- itemlabels$itemnr
  
```

```{r}
glimpse(df2)
```

```{r}
skim(df2)
```

### Ordering response categories

Here is a figure with categories as they are.

```{r}
#| fig-height: 8
df2 %>% 
  pivot_longer(everything(),
               values_to = "category",
               names_to = "itemnr") %>% 
  group_by(itemnr) %>% 
  count(category) %>% 
  left_join(.,itemlabels, # this adds the item description to the dataset
            by = "itemnr") %>% 
  ggplot(aes(x = category, y = n, fill = item)) +
  geom_col() +
  facet_wrap(~item, # makes a separate facet/plot for each item
             ncol = 1) +
  theme_ki() +
  scale_fill_manual(values = ki_color_palette,
                     guide = "none")

```

Reversing the response categories with `ggplot(aes(x = fct_rev(category)`

```{r}
#| fig-height: 8

df2 %>% 
  na.omit() %>% 
  pivot_longer(everything(),
               values_to = "category",
               names_to = "itemnr") %>% 
  group_by(itemnr) %>% 
  count(category) %>% 
  left_join(.,itemlabels, # this adds the item description to the dataset
            by = "itemnr") %>% 
  ### reverse response categories
  ggplot(aes(x = fct_rev(category), y = n, fill = item)) +
  geom_col() +
  facet_wrap(~item, # makes a separate facet/plot for each item
             ncol = 1) +
  theme_ki() +
  scale_fill_manual(values = ki_color_palette,
                     guide = "none")
```

And, just for reference, manually ordering the response categories. Here, we also added the category names to each plot facet by adding `scales = "free"` to the `facet_wrap()` call. Note that this frees the y axis to vary for each facet too, which can be less desirable. This can be easily solved by adding `scale_y_continuous(limits = c(0,150))` to have all facets range from 0 to 150.

```{r}
#| fig-height: 8

df2 %>% 
  na.omit() %>% 
  pivot_longer(everything(),
               values_to = "category",
               names_to = "itemnr") %>% 
  group_by(itemnr) %>% 
  count(category) %>% 
  left_join(.,itemlabels, # this adds the item description to the dataset
            by = "itemnr") %>% 
  mutate(category = factor(category, levels = c("Aldrig","Ibland","Ganska ofta",
                                                "Sällan","Mycket ofta","Alltid"))) %>%  ### order response categories
  ggplot(aes(x = category, y = n, fill = item)) +
  geom_col() +
  facet_wrap(~item, # makes a separate facet/plot for each item
             ncol = 1,
             scales = "free") +
  theme_ki() +
  scale_fill_manual(values = ki_color_palette,
                     guide = "none")
```

### Stacked bar plot

```{r}
df2 %>%
  na.omit() %>%
  pivot_longer(everything()) %>%
  dplyr::count(name, value) %>%
  mutate(Item = factor(name, levels = rev(names(df2))),
         value = factor(value)) %>%
  ggplot(aes(x = n, y = Item, fill = value)) +
  geom_col() +
  scale_fill_viridis_d("Category",
                       direction = -1) +
  labs(title = "Item responses",
       x = "Number of responses") +
  theme_ki()
```

### Recoding response categories to integers

Using `car::recode()`

```{r}
df2.recoded <- df2 %>% 
  mutate(across(ls1:ls6, ~ recode(.x,"'Aldrig'=0;
                                   'Sällan'=1;
                                   'Ibland'=2;
                                   'Ganska ofta'=3;
                                   'Mycket ofta'=4;
                                   'Alltid'=5",
                                  as.factor = FALSE)))
  
```


### Tile plot

```{r}
df2.recoded %>% 
  na.omit() %>% 
  pivot_longer(everything()) %>%
  dplyr::count(name, value) %>%
  mutate(name = factor(name, levels = rev(names(df2.recoded)))) %>%
  ggplot(aes(x = value, y = name, fill = n)) +
  geom_tile() +
  scale_fill_viridis_c(expression(italic(n)), limits = c(0, NA)) +
  scale_x_continuous("Response category", expand = c(0, 0), breaks = 0:max(df2.recoded, na.rm = T)) + # change breaks to fit number of response categories
  labs(y = "Items") +
  theme(axis.text.x = element_text(size = 8)) +
  geom_text(aes(label = n), colour = "orange") +
  theme_ki()
```


## Exercise?

Plot/table of age+gender for df2's spss data. Start with importing the whole dataset and figure out which variables are age and gender.

```
df3 <- read.spss("data/2023-04-26 Prevent OSA-enkat.sav", to.data.frame = T)
```