motion production
2023-03-06
knitr::opts_chunk$set(echo = TRUE)
library(afex)
library(emmeans)
library(sjPlot)
library(tidyverse)
library(ggstatsplot)data <- read.csv("data_clean2.csv", header = TRUE)
data <- subset(data, Participant_ID != 4 & Participant_ID != 12) Quick look
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.1.0
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(ggstatsplot)## You can cite this package as:
## Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
## Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167data %>% group_by(Participant_ID,trialType,indicatorType) %>% count() %>% pivot_wider(names_from = indicatorType, values_from = n)## # A tibble: 56 × 5
## # Groups: Participant_ID, trialType [56]
## Participant_ID trialType `0.2` `0.4` `0.6`
## <int> <chr> <int> <int> <int>
## 1 1 left 15 15 15
## 2 1 right 15 15 15
## 3 2 left 15 15 15
## 4 2 right 15 15 15
## 5 3 left 15 15 15
## 6 3 right 15 15 15
## 7 5 left 15 15 15
## 8 5 right 15 15 15
## 9 6 left 15 15 15
## 10 6 right 15 15 15
## # … with 46 more rows# create a plot
ggwithinstats(
data = data %>% mutate(sideXindicator=interaction(trialType,indicatorType)) %>% na.omit(),
x = sideXindicator,
y = amount_accuracy,
type = "np",
p.adjust.method="none",
pairwise.comparisons = FALSE,
)
# Participant 11 and 22 might be suspect (>1.5 IQR)
ggwithinstats(
data = data %>% na.omit() %>% mutate(sideXindicator=interaction(trialType,indicatorType)) %>% group_by(Participant_ID,sideXindicator) %>% summarise(maa=mean(amount_accuracy)),
x = sideXindicator,
y = maa,
type = "np",
p.adjust.method="none",
pairwise.comparisons = FALSE,
outlier.tagging = TRUE,
outlier.label=Participant_ID
)## `summarise()` has grouped output by 'Participant_ID'. You can override using
## the `.groups` argument.
grouped_ggwithinstats(
data = data %>% na.omit() %>% group_by(Participant_ID,trialType,indicatorType) %>% summarise(maa=mean(amount_accuracy)),
x = trialType,
y = maa,
grouping.var = indicatorType,
type = "np",
)## `summarise()` has grouped output by 'Participant_ID', 'trialType'. You can
## override using the `.groups` argument.
Continous, ordered factor or factor
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(sjPlot)## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!# continous
a0 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a1 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 + trialType | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a2 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 + indicatorType | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
tab_model(a0,a1,a2)| amount_accuracy | amount_accuracy | amount_accuracy | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 0.33 | 0.22 – 0.43 | <0.001 | 0.33 | 0.21 – 0.44 | <0.001 | 0.33 | 0.26 – 0.39 | <0.001 |
| trialType [right] | 0.00 | -0.05 – 0.06 | 0.909 | 0.00 | -0.06 – 0.07 | 0.919 | 0.00 | -0.05 – 0.05 | 0.906 |
| indicatorType | 0.64 | 0.55 – 0.73 | <0.001 | 0.64 | 0.55 – 0.72 | <0.001 | 0.64 | 0.35 – 0.92 | <0.001 |
|
trialType [right] × indicatorType |
-0.03 | -0.16 – 0.10 | 0.622 | -0.03 | -0.16 – 0.09 | 0.620 | -0.03 | -0.15 – 0.08 | 0.583 |
| Random Effects | |||||||||
| σ2 | 0.07 | 0.07 | 0.06 | ||||||
| τ00 | 0.07 Participant_ID | 0.08 Participant_ID | 0.02 Participant_ID | ||||||
| τ11 | 0.01 Participant_ID.trialTyperight | 0.53 Participant_ID.indicatorType | |||||||
| ρ01 | -0.45 Participant_ID | -0.38 Participant_ID | |||||||
| ICC | 0.51 | 0.52 | 0.61 | ||||||
| N | 28 Participant_ID | 28 Participant_ID | 28 Participant_ID | ||||||
| Observations | 2516 | 2516 | 2516 | ||||||
| Marginal R2 / Conditional R2 | 0.066 / 0.544 | 0.066 / 0.554 | 0.066 / 0.637 | ||||||
plot_model(a2, type = "pred", terms = c("indicatorType","trialType"))
plot_model(a2, type="re",sort.est = "indicatorType")## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.5.3
## Current Matrix version is 1.5.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package## Warning in checkDepPackageVersion(dep_pkg = "TMB"): Package version inconsistency detected.
## glmmTMB was built with TMB version 1.9.1
## Current TMB version is 1.9.2
## Please re-install glmmTMB from source or restore original 'TMB' package (see '?reinstalling' for more information)
data$indicatorType <- factor(data$indicatorType, levels = c("0.2", "0.4", "0.6"), ordered = TRUE)
a0 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a1 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 + trialType | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a2 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 + indicatorType | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
tab_model(a0,a1,a2)| amount_accuracy | amount_accuracy | amount_accuracy | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 0.58 | 0.48 – 0.68 | <0.001 | 0.58 | 0.47 – 0.69 | <0.001 | 0.58 | 0.48 – 0.68 | <0.001 |
| trialType [right] | -0.01 | -0.03 – 0.01 | 0.366 | -0.01 | -0.05 – 0.03 | 0.609 | -0.01 | -0.03 – 0.01 | 0.264 |
| indicatorType [linear] | 0.18 | 0.15 – 0.21 | <0.001 | 0.18 | 0.15 – 0.21 | <0.001 | 0.18 | 0.10 – 0.26 | <0.001 |
| indicatorType [quadratic] | 0.00 | -0.02 – 0.03 | 0.804 | 0.00 | -0.02 – 0.03 | 0.804 | 0.00 | -0.05 – 0.06 | 0.900 |
|
trialType [right] × indicatorType [linear] |
-0.01 | -0.05 – 0.03 | 0.623 | -0.01 | -0.04 – 0.03 | 0.620 | -0.01 | -0.04 – 0.02 | 0.556 |
|
trialType [right] × indicatorType [quadratic] |
-0.00 | -0.04 – 0.03 | 0.944 | -0.00 | -0.04 – 0.03 | 0.946 | -0.00 | -0.03 – 0.03 | 0.935 |
| Random Effects | |||||||||
| σ2 | 0.07 | 0.07 | 0.05 | ||||||
| τ00 | 0.07 Participant_ID | 0.08 Participant_ID | 0.07 Participant_ID | ||||||
| τ11 | 0.01 Participant_ID.trialTyperight | 0.04 Participant_ID.indicatorType.L | |||||||
| 0.02 Participant_ID.indicatorType.Q | |||||||||
| ρ01 | -0.45 Participant_ID | 0.88 | |||||||
| -0.08 | |||||||||
| ICC | 0.51 | 0.52 | 0.65 | ||||||
| N | 28 Participant_ID | 28 Participant_ID | 28 Participant_ID | ||||||
| Observations | 2516 | 2516 | 2516 | ||||||
| Marginal R2 / Conditional R2 | 0.066 / 0.544 | 0.066 / 0.554 | 0.066 / 0.675 | ||||||
data$indicatorType <- factor(data$indicatorType, ordered = FALSE)
a0 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a1 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 + trialType | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a2 <- lmer(amount_accuracy ~ 1 + trialType * indicatorType + (1 + indicatorType | Participant_ID), data, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
tab_model(a0,a1,a2)| amount_accuracy | amount_accuracy | amount_accuracy | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 0.45 | 0.35 – 0.56 | <0.001 | 0.45 | 0.34 – 0.56 | <0.001 | 0.45 | 0.39 – 0.52 | <0.001 |
| trialType [right] | -0.00 | -0.04 – 0.03 | 0.839 | -0.00 | -0.05 – 0.04 | 0.878 | -0.00 | -0.03 – 0.03 | 0.793 |
| indicatorType [0.4] | 0.12 | 0.09 – 0.16 | <0.001 | 0.12 | 0.09 – 0.16 | <0.001 | 0.12 | 0.03 – 0.21 | 0.006 |
| indicatorType [0.6] | 0.25 | 0.22 – 0.29 | <0.001 | 0.25 | 0.22 – 0.29 | <0.001 | 0.25 | 0.14 – 0.37 | <0.001 |
|
trialType [right] × indicatorType [0.4] |
-0.00 | -0.06 – 0.05 | 0.853 | -0.00 | -0.06 – 0.05 | 0.850 | -0.00 | -0.05 – 0.04 | 0.823 |
|
trialType [right] × indicatorType [0.6] |
-0.01 | -0.06 – 0.04 | 0.623 | -0.01 | -0.06 – 0.04 | 0.620 | -0.01 | -0.06 – 0.03 | 0.556 |
| Random Effects | |||||||||
| σ2 | 0.07 | 0.07 | 0.05 | ||||||
| τ00 | 0.07 Participant_ID | 0.08 Participant_ID | 0.03 Participant_ID | ||||||
| τ11 | 0.01 Participant_ID.trialTyperight | 0.05 Participant_ID.indicatorType0.4 | |||||||
| 0.09 Participant_ID.indicatorType0.6 | |||||||||
| ρ01 | -0.45 Participant_ID | 0.19 | |||||||
| 0.55 | |||||||||
| ICC | 0.51 | 0.52 | 0.65 | ||||||
| N | 28 Participant_ID | 28 Participant_ID | 28 Participant_ID | ||||||
| Observations | 2516 | 2516 | 2516 | ||||||
| Marginal R2 / Conditional R2 | 0.066 / 0.544 | 0.066 / 0.554 | 0.066 / 0.675 | ||||||
VE
data.aggr <- data %>% na.omit() %>% group_by(Participant_ID,indicatorType,trialType) %>% mutate(VE=sd(amount_accuracy))
a0 <- lmer(VE ~ 1 + trialType * indicatorType + (1 | Participant_ID), data.aggr, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a1 <- lmer(VE ~ 1 + trialType * indicatorType + (1 + trialType | Participant_ID), data.aggr, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
a2 <- lmer(VE ~ 1 + trialType * indicatorType + (1 + indicatorType | Participant_ID), data.aggr, REML = FALSE, control = lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 30000000)))
tab_model(a0,a1,a2)| VE | VE | VE | |||||||
|---|---|---|---|---|---|---|---|---|---|
| Predictors | Estimates | CI | p | Estimates | CI | p | Estimates | CI | p |
| (Intercept) | 0.16 | 0.13 – 0.18 | <0.001 | 0.16 | 0.13 – 0.19 | <0.001 | 0.16 | 0.14 – 0.18 | <0.001 |
| trialType [right] | 0.01 | -0.00 – 0.02 | 0.108 | 0.01 | -0.01 – 0.03 | 0.503 | 0.01 | 0.00 – 0.01 | 0.030 |
| indicatorType [0.4] | 0.04 | 0.03 – 0.05 | <0.001 | 0.04 | 0.03 – 0.05 | <0.001 | 0.04 | 0.01 – 0.07 | 0.004 |
| indicatorType [0.6] | 0.07 | 0.06 – 0.08 | <0.001 | 0.07 | 0.06 – 0.08 | <0.001 | 0.07 | 0.04 – 0.10 | <0.001 |
|
trialType [right] × indicatorType [0.4] |
-0.02 | -0.03 – -0.01 | <0.001 | -0.02 | -0.03 – -0.01 | <0.001 | -0.02 | -0.03 – -0.01 | <0.001 |
|
trialType [right] × indicatorType [0.6] |
-0.00 | -0.01 – 0.01 | 0.874 | -0.00 | -0.01 – 0.01 | 0.875 | -0.00 | -0.01 – 0.01 | 0.832 |
| Random Effects | |||||||||
| σ2 | 0.00 | 0.00 | 0.00 | ||||||
| τ00 | 0.00 Participant_ID | 0.01 Participant_ID | 0.00 Participant_ID | ||||||
| τ11 | 0.00 Participant_ID.trialTyperight | 0.01 Participant_ID.indicatorType0.4 | |||||||
| 0.01 Participant_ID.indicatorType0.6 | |||||||||
| ρ01 | -0.47 Participant_ID | 0.12 | |||||||
| 0.03 | |||||||||
| ICC | 0.55 | 0.63 | 0.77 | ||||||
| N | 28 Participant_ID | 28 Participant_ID | 28 Participant_ID | ||||||
| Observations | 2516 | 2516 | 2516 | ||||||
| Marginal R2 / Conditional R2 | 0.096 / 0.596 | 0.096 / 0.666 | 0.096 / 0.796 | ||||||
plot_model(a2, type = "pred", terms = c("indicatorType","trialType"))
fitting a linear regression to ..
reg_results <- data %>% na.omit() %>%
group_by(Participant_ID, trialType) %>%
summarize(intercept = coef(lm(amount_accuracy ~ indicatorType))[1],
slope = coef(lm(amount_accuracy ~ indicatorType))[2],
mean_dependent_variable = mean(amount_accuracy))## `summarise()` has grouped output by 'Participant_ID'. You can override using
## the `.groups` argument.ggwithinstats(
data = reg_results,
x = trialType,
y = slope,
type = "p"
)## Adding missing grouping variables: `Participant_ID`
ggplot(data = data,
aes(x=indicatorType,y=amount_accuracy, group = trialType, color = as.factor(Participant_ID)))+
geom_smooth(aes(linetype=trialType),method = "lm") + facet_wrap(~Participant_ID)## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).