TPDT Analysis

Tidy data
Descriptive plots
Assumption tests
Two way mixed ANOVA results
- Post hoc tests
Interaction plot
Shortcuts

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

library(readxl)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggplot2)
library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(rstatix)

## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter

library(ggpubr)

df <- read_excel("TPDT.xlsx")

## New names:
## • `` -> `...16`

df

Tidy data

df_long <- df %>%
  pivot_longer(
    cols = 7:11,              # these are the trial columns
    names_to = "Trial",
    values_to = "Measurement"
  ) %>%
  mutate(
    Trial = as.integer(Trial)   # Convert trial number to numeric
  ) %>%
  select(ID, Location, Pregnant, Trial, Measurement)  # Keep relevant columns

df_long %>%
  group_by(Location, Pregnant) %>%
  summarise(
    n = sum(!is.na(Measurement)),
    mean = mean(Measurement, na.rm = TRUE),
    sd = sd(Measurement, na.rm = TRUE),
    median = median(Measurement, na.rm = TRUE),
    min = min(Measurement, na.rm = TRUE),
    max = max(Measurement, na.rm = TRUE)
  )

## `summarise()` has grouped output by 'Location'. You can override using the
## `.groups` argument.

df_long %>% #group by participant
  group_by(ID, Location) %>%
  summarise(
    trials_completed = sum(!is.na(Measurement)),
    mean = mean(Measurement, na.rm = TRUE),
    sd = sd(Measurement, na.rm = TRUE)
  )

## `summarise()` has grouped output by 'ID'. You can override using the `.groups`
## argument.

df_long$Pregnant <- factor(df_long$Pregnant,
                          levels = c(0, 1),
                          labels = c("Not Pregnant", "Pregnant"))

df_long %>%
  distinct(ID, Pregnant) %>%  # remove duplicate rows per participant
  count(Pregnant)             # count participants per group

threshold_summary <- df_long %>%
  filter(!is.na(Measurement), !Trial >=4) %>%                # Exclude missing trials and trials above 3
  group_by(ID, Pregnant, Location) %>%
  summarise(
    n_trials = n(),
    mean_threshold = mean(Measurement),
    median_threshold = median(Measurement),
    .groups = "drop"
  )

# View the result
print(threshold_summary)

## # A tibble: 92 × 6
##       ID Pregnant Location n_trials mean_threshold median_threshold
##    <dbl> <fct>    <chr>       <int>          <dbl>            <dbl>
##  1     1 Pregnant Abdomen         3           46.3             50  
##  2     1 Pregnant Forearm         3           38.3             40  
##  3     3 Pregnant Abdomen         3           32.7             32  
##  4     3 Pregnant Forearm         3           15.3             15  
##  5     4 Pregnant Abdomen         3           33.7             34  
##  6     4 Pregnant Forearm         3           34               35  
##  7     6 Pregnant Abdomen         3           31.7             30  
##  8     6 Pregnant Forearm         2           31.5             31.5
##  9     8 Pregnant Abdomen         3           36.3             34  
## 10     8 Pregnant Forearm         2           57               57  
## # ℹ 82 more rows

threshold_forearm <- threshold_summary %>%
  filter(Location == "Forearm") 
  
threshold_abdomen <- threshold_summary %>%
  filter(Location == "Abdomen")

Descriptive plots

# Scatter plot: Mean threshold by participant and location
ggplot(threshold_forearm) +
  aes(x = factor(Pregnant), y = mean_threshold, color = Pregnant) +
  geom_jitter() +
    labs(
    title = "Mean Two-Point Discrimination Threshold - Forearm",
    x = "Pregnancy status",
    y = "Mean Threshold (mm)"
  ) +
  theme(legend.position = "none")

# Scatter plot: Mean threshold by participant and location
ggplot(threshold_abdomen) +
  aes(x = factor(Pregnant), y = mean_threshold, color = Pregnant) +
  geom_jitter() +
    labs(
    title = "Mean Two-Point Discrimination Threshold - Abdomen",
    x = "Pregnancy status",
    y = "Mean Threshold (mm)"
  ) +
  theme(legend.position = "none")

Assumption tests

# Outliers
threshold_summary %>%
  group_by(Location, Pregnant) %>%
  identify_outliers(mean_threshold) #no extreme values

# Normality
threshold_summary %>%
  group_by(Location, Pregnant) %>%
  shapiro_test(mean_threshold) #can assume normality if P > 0.05 for each cell

ggqqplot(threshold_summary, "mean_threshold", ggtheme = theme_bw()) +
  facet_grid(Location ~ Pregnant)

Note on Outliers - Values above Q3 + 1.5xIQR or below Q1 - 1.5xIQR are considered as outliers. Values above Q3 + 3xIQR or below Q1 - 3xIQR are considered as extreme points (or extreme outliers).

Two way mixed ANOVA results

anova_results <- threshold_summary %>%
  anova_test(dv = mean_threshold, wid = ID, within = Location, between = Pregnant) # Two-way mixed ANOVA
get_anova_table(anova_results)

Post hoc tests

## Post-hoc tests
# Simple main effect of group variable
# Effect of group at each location
one.way <- threshold_summary %>%
  group_by(Location) %>%
  anova_test(dv = mean_threshold, wid = ID, between = Pregnant) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")
one.way

# Pairwise comparison between group levels
between_pwc <- threshold_summary %>%
  group_by(Location) %>%
  pairwise_t_test(mean_threshold ~ Pregnant, paired = FALSE, p.adjust.method = "bonferroni")
between_pwc

# Simple main effect of location variable
# Effect of location for each group
one.way2 <- threshold_summary %>%
  group_by(Pregnant) %>%
  anova_test(dv = mean_threshold, wid = ID, between = Location) %>%
  get_anova_table() %>%
  adjust_pvalue(method = "bonferroni")

## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning:
## ! The 'wid' column contains duplicate ids across between-subjects variables. Automatic unique id will be created
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.

one.way2

# Pairwise comparison between locations at each group level
within_pwc <- threshold_summary %>%
  group_by(Pregnant) %>%
  pairwise_t_test(mean_threshold ~ Location, paired = TRUE, p.adjust.method = "bonferroni")
within_pwc

Interaction plot

threshold_summary %>%
  group_by(Pregnant, Location) %>%
  summarise(mean = mean(mean_threshold), se = sd(mean_threshold)/sqrt(n())) %>%
  ggplot(aes(x = Location, y = mean, color = Pregnant, group = Pregnant)) +
  geom_line() +
  geom_point() +
  geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) +
      labs(
    title = "Two Point Discrimination Test - Inteaction Plot",
    x = "Location",
    y = "Mean Threshold (mm)"
  ) +
  theme_minimal()

## `summarise()` has grouped output by 'Pregnant'. You can override using the
## `.groups` argument.

Shortcuts

Execute Chunk Cmd+Shift+Enter.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I or Ctrl+Alt+I with keyboard.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.