Data import & descriptives

1 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)

2 Data cleaning

Image credit: https://kbroman.org/Talk_DataCleaning/data_cleaning.pdf

3 Setting up for data analysis

Let’s load packages/libraries.

Code

# 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.

Code

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

3.1 Adaptions

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().

Code

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

4 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, 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.

Code

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.

Code

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”.

Code

df %>% 
  select(contains("K10"))

# A tibble: 106 × 7
   K10_Score_GP GPK10_coded TK10_t0 K10_di_t0    TK10_t1 TK10_t2 TK10_t3
          <dbl>       <dbl>   <dbl> <chr>          <dbl>   <dbl>   <dbl>
 1           33           1      30 30+               25      30      21
 2           24           0      30 30+               27      31      32
 3           36           1      38 30+               28      18      23
 4           29           0      27 less than 30      21      24      27
 5           27           0      22 less than 30      32      NA      23
 6           29           0      16 less than 30      15      16      17
 7           23           0      23 less than 30      20      15      17
 8           27           0      30 30+               42      27      19
 9           29           0      25 less than 30      17      16      15
10           30           1      29 less than 30      21      17      21
# ℹ 96 more rows

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.

Code

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

5 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.

Code

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.

Code

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.

Code

summary(df$age)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  25.00   37.00   40.50   42.24   47.00   64.00

5.1 Tables

Code

df %>% 
  select(age) %>% 
  tbl_summary()

Characteristic	N = 106¹
age	41 (37, 47)
¹ Median (IQR)

Code

df %>% 
  select(age,Group) %>% 
  tbl_summary(by = Group,
              statistic = list(all_continuous() ~ "{mean} ({sd}, {min}-{max})"))

Characteristic	Control, N = 55¹	MiCBT, N = 51¹
age	43 (8, 26-63)	42 (8, 25-64)
¹ Mean (SD, Minimum-Maximum)

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.

Code

age.ttest <- t.test(age ~ Group, data = df)
age.ttest


    Welch Two Sample t-test

data:  age by Group
t = 0.64071, df = 101.07, p-value = 0.5232
alternative hypothesis: true difference in means between group Control and group MiCBT is not equal to 0
95 percent confidence interval:
 -2.061734  4.028936
sample estimates:
mean in group Control   mean in group MiCBT 
             42.70909              41.72549

Code

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

estimate	estimate1	estimate2	statistic	p.value	parameter	conf.low	conf.high	method	alternative
0.98	42.71	41.73	0.64	0.52	101.07	-2.06	4.03	Welch Two Sample t-test	two.sided

5.2 Figures

Base R examples.

Code

hist(df$age)

Code

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.

Code

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.

Code

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()

Practice

Make a new plot for intervention group age.

5.3 tidy filter/select

Note

Two key things to learn:

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

5.3.1 Examples

This filters for inclusion of only one municipality:

df %>% 
  filter(municipality == "Vallentuna")

This filters out (excludes) multiple municipalities:

df %>% 
  filter(!municipality %in% c("Vallentuna","Vaxholm"))

5.3.2 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).

Code

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.

Code

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)

Practice

Plot age grouped by gender (bonus: facet_wrap by Group)

6 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.

Code

df %>% 
  select(ends_with("t0")) %>% 
  names()

 [1] "TK10_t0"           "Temotreg_t0"       "Tequanimity_t0"   
 [4] "Tsocialskills_t0"  "Tdistresstol_t0"   "Ttakerespons_t0"  
 [7] "Tinterpersonal_t0" "Tnoticing_t0"      "Tnotdisract_t0"   
[10] "Tnotworry_t0"      "Tattenreg_t0"      "Temotaware_t0"    
[13] "Tbodylisten_t0"    "Ttrusting_t0"      "Tselfreg_t0"      
[16] "TNAS_t0"           "TNAS_MEAN_t0"      "EQ_t0"            
[19] "TFs_t0"            "TSWLS_t0"          "TDASS_t0"         
[22] "TANX_t0"           "TSTRESS_t0"        "TDEP_t0"          
[25] "TMSES_t0"          "TMAIA_t0"          "MSES_Et0"         
[28] "MSES_It0"

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

Code

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

 [1] "K10_t0"           "emotreg_t0"       "equanimity_t0"    "socialskills_t0" 
 [5] "distresstol_t0"   "takerespons_t0"   "interpersonal_t0" "noticing_t0"     
 [9] "notdisract_t0"    "notworry_t0"      "attenreg_t0"      "emotaware_t0"    
[13] "bodylisten_t0"    "trusting_t0"      "selfreg_t0"       "NAS_t0"          
[17] "NAS_MEAN_t0"      "EQ_t0"            "Fs_t0"            "SWLS_t0"         
[21] "DASS_t0"          "ANX_t0"           "STRESS_t0"        "DEP_t0"          
[25] "MSES_t0"          "MAIA_t0"          "MSES_Et0"         "MSES_It0"

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”.

Code

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

7 Skim

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

Code

df %>% 
  select(starts_with(c("ANX","DEP","STRESS"))) %>% 
  skim()

Data summary
Name	Piped data
Number of rows	106
Number of columns	12
_______________________
Column type frequency:
numeric	12
________________________
Group variables	None

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
ANX_t0	1	0.99	11.20	8.02	0	6.0	8	16	36	▇▅▃▂▁
ANX_t1	13	0.88	10.06	8.34	0	4.0	8	14	36	▇▆▂▂▁
ANX_t2	22	0.79	8.33	7.51	0	2.0	6	14	32	▇▂▃▁▁
ANX_t3	20	0.81	9.07	8.87	0	2.0	8	12	38	▇▆▁▁▁
DEP_t0	1	0.99	16.46	10.49	0	8.0	16	22	42	▇▇▆▂▂
DEP_t1	13	0.88	14.45	11.26	0	4.0	12	22	42	▇▅▃▂▂
DEP_t2	22	0.79	12.57	10.21	0	4.0	10	18	42	▇▅▂▂▁
DEP_t3	20	0.81	12.91	9.48	0	6.0	12	18	42	▇▇▃▂▁
STRESS_t0	1	0.99	20.70	8.46	2	14.0	22	26	42	▂▇▇▃▁
STRESS_t1	13	0.88	18.22	9.16	0	12.0	16	22	42	▃▇▆▂▂
STRESS_t2	22	0.79	16.19	8.93	0	10.0	14	22	42	▃▇▆▂▁
STRESS_t3	20	0.81	16.79	8.80	0	10.5	14	22	42	▂▇▅▂▁

And grouped:

Code

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

Data summary
Name	Piped data
Number of rows	106
Number of columns	13
_______________________
Column type frequency:
numeric	12
________________________
Group variables	Group

Variable type: numeric

skim_variable	Group	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
ANX_t0	Control	0	1.00	11.38	8.35	0	4.0	8	17.0	36	▇▆▃▂▁
ANX_t1	Control	2	0.96	11.66	9.62	0	4.0	10	16.0	36	▇▇▂▃▁
ANX_t2	Control	10	0.82	10.22	8.30	0	4.0	8	16.0	32	▇▃▅▂▁
ANX_t3	Control	9	0.84	9.74	9.05	0	4.0	8	13.5	38	▇▇▂▁▁
DEP_t0	Control	0	1.00	16.76	10.07	0	11.0	16	22.0	40	▆▇▇▂▂
DEP_t1	Control	2	0.96	16.38	12.04	0	6.0	14	24.0	42	▇▅▅▂▃
DEP_t2	Control	10	0.82	15.16	10.14	2	8.0	14	20.0	42	▇▇▃▂▂
DEP_t3	Control	9	0.84	14.04	9.94	0	6.5	12	19.5	42	▆▇▅▂▁
STRESS_t0	Control	0	1.00	21.13	8.10	6	16.0	22	26.0	42	▅▇▇▃▁
STRESS_t1	Control	2	0.96	20.79	10.31	0	14.0	18	26.0	42	▂▇▆▅▃
STRESS_t2	Control	10	0.82	19.69	8.90	4	14.0	20	24.0	42	▃▇▇▃▁
STRESS_t3	Control	9	0.84	19.39	8.64	4	14.0	18	23.5	42	▂▇▆▂▁
ANX_t0	MiCBT	1	0.98	11.00	7.72	0	6.0	9	16.0	34	▇▃▅▁▁
ANX_t1	MiCBT	11	0.78	7.95	5.72	0	4.0	6	10.0	26	▆▇▁▂▁
ANX_t2	MiCBT	12	0.76	6.15	5.86	0	2.0	6	9.0	20	▇▅▂▁▂
ANX_t3	MiCBT	11	0.78	8.30	8.72	0	2.0	6	10.5	38	▇▃▁▁▁
DEP_t0	MiCBT	1	0.98	16.12	11.02	0	8.0	15	20.0	42	▇▇▅▂▂
DEP_t1	MiCBT	11	0.78	11.90	9.70	0	4.0	11	18.0	34	▇▃▅▂▂
DEP_t2	MiCBT	12	0.76	9.59	9.57	0	2.0	6	13.0	34	▇▂▁▂▁
DEP_t3	MiCBT	11	0.78	11.60	8.88	0	4.0	10	16.0	34	▇▇▅▂▂
STRESS_t0	MiCBT	1	0.98	20.24	8.90	2	14.0	20	26.0	40	▃▆▇▆▁
STRESS_t1	MiCBT	11	0.78	14.80	5.92	6	10.0	15	20.0	28	▇▅▅▆▂
STRESS_t2	MiCBT	12	0.76	12.15	7.16	0	8.0	12	15.0	32	▃▇▅▂▁
STRESS_t3	MiCBT	11	0.78	13.80	8.10	0	9.5	12	16.5	36	▃▇▂▂▁

8 Correlation & visualization

8.1 Correlation matrix

Code

df %>% 
  select(ANX_t0,DEP_t0,STRESS_t0) %>% 
  ggpairs()

8.2 Correlation grouped

Code

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()

8.3 easystats

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.

Code

df %>% 
  cor_test("ANX_t0", "DEP_t0") %>% 
  plot()

With some styling/theming.

Code

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!

9 Missing data

Code

df %>% 
  select(starts_with(c("ANX","DEP","STRESS"))) %>% 
  vis_dat(palette = "qual")

9.1 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

10 Outliers

Code

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

OK: No outliers detected.
- Based on the following method and threshold: mahalanobis (34.528).
- For variables: id, ANX_t0, ANX_t1, ANX_t2, ANX_t3, DEP_t0, DEP_t1, DEP_t2, DEP_t3, STRESS_t0, STRESS_t1, STRESS_t2, STRESS_t3

Should we try other methods? See ?check_outliers.

11 Wide to long format

Almost everything in R likes long format. Let’s look at the variables ending with “t”.

Code

df %>% 
  select(ends_with("t0")) %>% 
  names()

 [1] "K10_t0"           "emotreg_t0"       "equanimity_t0"    "socialskills_t0" 
 [5] "distresstol_t0"   "takerespons_t0"   "interpersonal_t0" "noticing_t0"     
 [9] "notdisract_t0"    "notworry_t0"      "attenreg_t0"      "emotaware_t0"    
[13] "bodylisten_t0"    "trusting_t0"      "selfreg_t0"       "NAS_t0"          
[17] "NAS_MEAN_t0"      "EQ_t0"            "Fs_t0"            "SWLS_t0"         
[21] "DASS_t0"          "ANX_t0"           "STRESS_t0"        "DEP_t0"          
[25] "MSES_t0"          "MAIA_t0"          "MSES_Et0"         "MSES_It0"

Code

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

 [1] "emotreg_t1"       "equanimity_t1"    "socialskills_t1"  "distresstol_t1"  
 [5] "takerespons_t1"   "interpersonal_t1" "MSES_t1"          "noticing_t1"     
 [9] "notdisract_t1"    "notworry_t1"      "attenreg_t1"      "emotaware_t1"    
[13] "bodylisten_t1"    "trusting_t1"      "selfreg_t1"       "NAS_t1"          
[17] "NAS_MEAN_t1"      "EQ_t1"            "Fs_t1"            "SWLS_t1"         
[21] "K10_t1"           "DASS_t1"          "ANX_t1"           "STRESS_t1"       
[25] "DEP_t1"           "MAIA_t1"          "MSES_Et1"         "MSES_It1"

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

Code

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.

Code

glimpse(df.long)

Rows: 1,272
Columns: 5
$ id      <dbl> 830, 830, 830, 830, 830, 830, 830, 830, 830, 830, 830, 830, 57…
$ Group   <chr> "Control", "Control", "Control", "Control", "Control", "Contro…
$ measure <chr> "ANX", "ANX", "ANX", "ANX", "DEP", "DEP", "DEP", "DEP", "STRES…
$ time    <chr> "0", "1", "2", "3", "0", "1", "2", "3", "0", "1", "2", "3", "0…
$ value   <dbl> 16, 22, 22, 8, 22, 24, 22, 14, 22, 30, 28, 16, 8, 6, 0, 8, 14,…

12 Done!

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

Code

write_csv(df.long,"data/dataLong.csv")

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.

13 Bonus - dealing with questionnaire data

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.

Code

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

Code

glimpse(df2)

Rows: 581
Columns: 6
$ ls1 <fct> Ganska ofta, Sällan, Ganska ofta, Alltid, NA, NA, NA, Ibland, Gans…
$ ls2 <fct> Mycket ofta, Ganska ofta, Mycket ofta, Mycket ofta, NA, NA, NA, Sä…
$ ls3 <fct> Ganska ofta, Ganska ofta, Ganska ofta, Alltid, NA, NA, NA, Aldrig,…
$ ls4 <fct> Ganska ofta, Ganska ofta, Alltid, Alltid, NA, NA, NA, Ibland, Gans…
$ ls5 <fct> Ganska ofta, Ganska ofta, Alltid, Mycket ofta, NA, NA, NA, Ibland,…
$ ls6 <fct> Ganska ofta, Sällan, Mycket ofta, Alltid, NA, NA, NA, Ibland, Ibla…

Code

skim(df2)

Data summary
Name	df2
Number of rows	581
Number of columns	6
_______________________
Column type frequency:
factor	6
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
ls1	96	0.83	FALSE	6	Gan: 139, Myc: 106, Ibl: 102, Säl: 69
ls2	97	0.83	FALSE	6	Myc: 144, Gan: 124, Ibl: 78, Säl: 60
ls3	97	0.83	FALSE	6	Gan: 105, Ibl: 104, Säl: 94, Myc: 86
ls4	97	0.83	FALSE	6	Myc: 147, All: 125, Gan: 96, Ibl: 63
ls5	98	0.83	FALSE	6	Myc: 118, Gan: 110, All: 86, Ibl: 85
ls6	96	0.83	FALSE	6	All: 132, Myc: 123, Gan: 95, Ibl: 66

13.1 Ordering response categories

Here is a figure with categories as they are.

Code

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)

Code

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.

Code

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")

13.2 Stacked bar plot

Code

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()

13.3 Recoding response categories to integers

Using car::recode()

Code

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)))

13.4 Tile plot

Code

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()

14 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)