TLMM_JointActionTask
Giovanna Cuomo
2024-08-29
CODE BOOK & DESCRIPTION
This script contains the analyses conducted on data related to participants who took part in our study “TMS B&W” for the Joint Action Task. In three different experiments, participants were presented with two virtual partners during the joint action task: an out-group (Black) or an in-group (White) avatar (the control condition). The instruction was to synchronize with the partner to gasp together a bottle-shaped object in front of them, performing the Same or Opposite movement relative to the avatar. A recent study has already proven how a partner’s racial membership can modulate interactive behaviors (Sacheli et al., 2015). To investigate the role of mPFC (vs Vertex, vs vPM) in coordinating with an out-group partner after a virtual inhibition of the area, before the joint action task, a continuous theta burst stimulation was administered (cTBS; Huang et al., 2005).
Hypotheses 1) Partner’s ethnicity modulates people’s ability to coordinate with them for shared goals. 2) mPFC inhibition reduces predictive simulation during joint action with an ethnic out-group member, deteriorating the interactive performance. 3) mPFC inhibition increases explicit and implicit racial bias.
Preprocessing The raw data were cleaned as follows: 1) trials in which participants failed to grasp the bottle-shaped object correctly were excluded from the analyses (i.e., when the touch-sensitive plates were not correctly pressed); 2) trials when participants did not follow the instructions, and grasped the object on the wrong half the trial was discarded; 3) in a similar way, we did not include in the analyses those trials in which participants started their movement before the audio instruction was given; 4) catch trials (i.e., when the virtual partner switched movement direction) from the analyses. 5) outliers trial values 2.5 standard deviations above or below the mean for each subject within each condition
Independent Variables: Site (mPFC vs vPM vs Vertex) Virtual partner Ethnicity (Black vs White) Movement (Same vs Opposite) Experiment (1, 2 or 3)
Dependent Variables: Grasping Asynchrony (GA), i.e., the absolute value of the difference between the time of the thumb-index contact on the bottle-shaped object for the virtual partner and the participant. Movement Time (MT), i.e., the difference between the time the participant released the starting button and the time of their thumb-index contact with the object. Reaction Time (RT), i.e., the time passed between the end of the audio instruction and the moment the participant released the stating button.
MLMMs In constructing the MLMMs that follows, we performed a “stepwise backward” selection process for the random effects. Starting from a structure that included all the relevant effects and interactions, we systematically removed terms based on singularity to optimize model fit and ensure the convergence of the model. Post-hoc comparisons were conducted utilizing the ‘Estimated Marginal Means’ R package (Lenth, 2020), employing the ‘emmeans’ and ‘emtrends’ functions to further investigate significant interactions. Bonferroni correction for multiple comparisons was applied.
First, we load packages we need for analyses #PACKAGES
library (tidyverse) #for data handling## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 0.3.5
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.4.1
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(readxl) #to read excel files
library(writexl) #to export data as excel file ## Warning: il pacchetto 'writexl' è stato creato con R versione 4.2.3library(lme4) ## Caricamento del pacchetto richiesto: Matrix
##
## Caricamento pacchetto: 'Matrix'
##
## I seguenti oggetti sono mascherati da 'package:tidyr':
##
## expand, pack, unpacklibrary(emmeans) ## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'library(ggplot2)
library(rstatix)## Warning: il pacchetto 'rstatix' è stato creato con R versione 4.2.3##
## Caricamento pacchetto: 'rstatix'
##
## Il seguente oggetto è mascherato da 'package:stats':
##
## filterlibrary(stats)
library(dplyr)
library(yarrr)## Warning: il pacchetto 'yarrr' è stato creato con R versione 4.2.3## Caricamento del pacchetto richiesto: jpeg
## Caricamento del pacchetto richiesto: BayesFactor## Warning: il pacchetto 'BayesFactor' è stato creato con R versione 4.2.3## Caricamento del pacchetto richiesto: coda
## ************
## Welcome to BayesFactor 0.9.12-4.4. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
##
## Type BFManual() to open the manual.
## ************
## Caricamento del pacchetto richiesto: circlize## Warning: il pacchetto 'circlize' è stato creato con R versione 4.2.3## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================
##
## yarrr v0.1.5. Citation info at citation('yarrr'). Package guide at yarrr.guide()
## Email me at Nathaniel.D.Phillips.is@gmail.com
##
## Caricamento pacchetto: 'yarrr'
##
## Il seguente oggetto è mascherato da 'package:ggplot2':
##
## diamondslibrary(xlsx)## Warning: il pacchetto 'xlsx' è stato creato con R versione 4.2.3options(scipen=999)##STUDY 1 & 3
dati_tms = read.xlsx("dati_tms.xlsx", sheetIndex = 1)Specifying the factors
dati_tms$sito = as.factor(dati_tms$sito)
dati_tms$bw = as.factor(dati_tms$bw)
dati_tms$Movement = as.factor(dati_tms$Movement)
dati_tms$Grasp = as.factor(dati_tms$Grasp)
dati_tms$Subject = as.factor(dati_tms$Subject)#Grasping Asynchrony
ASY_model = lmer(Asynchrony ~ Exp * bw * sito * Movement +
(1+ sito*bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dati_tms)
#summary(ASY_model)
isSingular(ASY_model)## [1] FALSEAnova(ASY_model)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Asynchrony
## Chisq Df Pr(>Chisq)
## Exp 1.7769 1 0.18253
## bw 0.1392 1 0.70911
## sito 0.5381 1 0.46320
## Movement 24.8172 1 0.0000006303 ***
## Exp:bw 0.6884 1 0.40670
## Exp:sito 4.2156 1 0.04005 *
## bw:sito 2.9614 1 0.08528 .
## Exp:Movement 3.3220 1 0.06836 .
## bw:Movement 0.0067 1 0.93494
## sito:Movement 0.0475 1 0.82747
## Exp:bw:sito 0.0797 1 0.77767
## Exp:bw:Movement 0.1429 1 0.70546
## Exp:sito:Movement 0.4671 1 0.49434
## bw:sito:Movement 0.0122 1 0.91189
## Exp:bw:sito:Movement 0.0382 1 0.84510
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residuals check
res = residuals(ASY_model)
plot(density(res))
plot(ASY_model)
qqnorm(residuals(ASY_model))
qqline(residuals(ASY_model))
significant interactions post-hocs
#emmeans/posthoc
emmeans_interaction <- emmeans(ASY_model, ~ Exp * sito)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 8356' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 8356)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 8356' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 8356)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsposthoc_interaction <- pairs(emmeans_interaction, adjust = "bonferroni")
summary(posthoc_interaction)## contrast estimate SE df z.ratio p.value
## Exp1 MPFC - Exp3 MPFC 16.37 11.17 Inf 1.467 0.8549
## Exp1 MPFC - Exp1 Vertex 11.49 6.56 Inf 1.752 0.4783
## Exp1 MPFC - Exp3 Vertex 10.48 12.00 Inf 0.873 1.0000
## Exp3 MPFC - Exp1 Vertex -4.89 12.38 Inf -0.395 1.0000
## Exp3 MPFC - Exp3 Vertex -5.89 5.36 Inf -1.099 1.0000
## Exp1 Vertex - Exp3 Vertex -1.01 13.14 Inf -0.077 1.0000
##
## Results are averaged over the levels of: bw, Movement
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 testsplot dataset
pirateplot(formula = Asynchrony ~ bw + sito + Exp, #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
main effect of Movement - plot
pirateplot(formula = Asynchrony ~ Movement, #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
#Since data were not normal: Alternative - GA log
dati_tms$Asynchrony_log = log(dati_tms$Asynchrony + 225)run LMM log
ASY_model3 = lmer(Asynchrony_log ~ Exp * bw * sito * Movement +
(1+bw*sito|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dati_tms)
summary(ASY_model3)## Linear mixed model fit by REML ['lmerMod']
## Formula: Asynchrony_log ~ Exp * bw * sito * Movement + (1 + bw * sito |
## Subject)
## Data: dati_tms
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 2416.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8437 -0.7972 -0.0586 0.6999 4.3998
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 0.005007 0.07076
## bwWhite 0.001245 0.03528 -0.09
## sitoVertex 0.001161 0.03407 -0.14 0.32
## bwWhite:sitoVertex 0.001264 0.03555 0.27 -0.06 0.50
## Residual 0.075788 0.27530
## Number of obs: 8356, groups: Subject, 40
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.928590 0.022395 264.724
## Exp3 -0.047240 0.028821 -1.639
## bwWhite -0.014499 0.021201 -0.684
## sitoVertex -0.036363 0.020992 -1.732
## MovementImitative -0.033922 0.019131 -1.773
## Exp3:bwWhite 0.007847 0.027265 0.288
## Exp3:sitoVertex 0.032818 0.027159 1.208
## bwWhite:sitoVertex 0.020897 0.028489 0.734
## Exp3:MovementImitative 0.005940 0.024654 0.241
## bwWhite:MovementImitative -0.006611 0.027069 -0.244
## sitoVertex:MovementImitative -0.003908 0.026875 -0.145
## Exp3:bwWhite:sitoVertex 0.012533 0.036799 0.341
## Exp3:bwWhite:MovementImitative 0.014454 0.034856 0.415
## Exp3:sitoVertex:MovementImitative 0.013385 0.034877 0.384
## bwWhite:sitoVertex:MovementImitative 0.004195 0.038104 0.110
## Exp3:bwWhite:sitoVertex:MovementImitative -0.014617 0.049205 -0.297##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itisSingular(ASY_model3)## [1] FALSEAnova(ASY_model3)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Asynchrony_log
## Chisq Df Pr(>Chisq)
## Exp 1.7123 1 0.19069
## bw 0.0634 1 0.80123
## sito 0.4471 1 0.50371
## Movement 22.2574 1 0.000002384 ***
## Exp:bw 0.9910 1 0.31950
## Exp:sito 4.9048 1 0.02678 *
## bw:sito 3.8262 1 0.05046 .
## Exp:Movement 1.7383 1 0.18736
## bw:Movement 0.0003 1 0.98646
## sito:Movement 0.0215 1 0.88341
## Exp:bw:sito 0.0362 1 0.84919
## Exp:bw:Movement 0.0837 1 0.77229
## Exp:sito:Movement 0.0603 1 0.80603
## bw:sito:Movement 0.0359 1 0.84964
## Exp:bw:sito:Movement 0.0882 1 0.76642
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1significant interactions post-hocs
#emmeans/posthoc
emmeans_interaction <- emmeans(ASY_model3, ~ Exp * sito)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 8356' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 8356)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 8356' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 8356)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsposthoc_interaction <- pairs(emmeans_interaction, adjust = "bonferroni")
summary(posthoc_interaction)## contrast estimate SE df z.ratio p.value
## Exp1 MPFC - Exp3 MPFC 0.03673 0.0246 Inf 1.491 0.8153
## Exp1 MPFC - Exp1 Vertex 0.02682 0.0149 Inf 1.803 0.4287
## Exp1 MPFC - Exp3 Vertex 0.02143 0.0266 Inf 0.805 1.0000
## Exp3 MPFC - Exp1 Vertex -0.00991 0.0275 Inf -0.360 1.0000
## Exp3 MPFC - Exp3 Vertex -0.01530 0.0122 Inf -1.258 1.0000
## Exp1 Vertex - Exp3 Vertex -0.00539 0.0293 Inf -0.184 1.0000
##
## Results are averaged over the levels of: bw, Movement
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 testsSignificant effects - plot
pirateplot(formula = Asynchrony_log ~ Movement, #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
mean and ds
library(dplyr)
summary_stats_asy <- dati_tms %>%
group_by(Movement) %>%
summarise(
mean_Asy = mean(Asynchrony_log, na.rm = TRUE), # Calcola la media
sd_Asy = sd(Asynchrony_log, na.rm = TRUE) # Calcola la deviazione standard
)
summary_stats_asy1 <- dati_tms %>%
group_by(Movement) %>%
summarise(
mean_Asy = mean(Asynchrony, na.rm = TRUE), # Calcola la media
sd_Asy = sd(Asynchrony, na.rm = TRUE) # Calcola la deviazione standard
)residuals check
res = residuals(ASY_model3)
plot(density(res))
plot(ASY_model3)
qqnorm(residuals(ASY_model3))
qqline(residuals(ASY_model3))
#Movement Time
MT_model = lmer(MovTime ~ Exp * bw * sito * Movement +
(1+sito*bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dati_tms)
summary(MT_model)## Linear mixed model fit by REML ['lmerMod']
## Formula: MovTime ~ Exp * bw * sito * Movement + (1 + sito * bw | Subject)
## Data: dati_tms
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 112493.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0744 -0.6967 -0.0999 0.6334 5.6873
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 14791 121.62
## sitoVertex 6987 83.59 -0.44
## bwWhite 2820 53.10 -0.13 0.59
## sitoVertex:bwWhite 7247 85.13 0.41 -0.66 -0.71
## Residual 36526 191.12
## Number of obs: 8415, groups: Subject, 40
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1281.4345 31.8410 40.245
## Exp3 -69.6006 41.0854 -1.694
## bwWhite -19.7585 18.7861 -1.052
## sitoVertex -4.2750 24.7532 -0.173
## MovementImitative -23.7644 13.2235 -1.797
## Exp3:bwWhite 28.5220 24.2157 1.178
## Exp3:sitoVertex 0.3497 32.0491 0.011
## bwWhite:sitoVertex 37.8366 28.3325 1.335
## Exp3:MovementImitative 8.2832 17.0579 0.486
## bwWhite:MovementImitative 6.3111 18.6848 0.338
## sitoVertex:MovementImitative 14.3453 18.5935 0.772
## Exp3:bwWhite:sitoVertex -44.5937 36.6565 -1.217
## Exp3:bwWhite:MovementImitative 1.9502 24.0864 0.081
## Exp3:sitoVertex:MovementImitative -6.1128 24.1637 -0.253
## bwWhite:sitoVertex:MovementImitative -8.6684 26.3108 -0.329
## Exp3:bwWhite:sitoVertex:MovementImitative -0.7008 34.0238 -0.021##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itisSingular(MT_model)## [1] FALSEAnova(MT_model)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MovTime
## Chisq Df Pr(>Chisq)
## Exp 3.3224 1 0.068341 .
## bw 0.3259 1 0.568096
## sito 0.2354 1 0.627531
## Movement 8.2647 1 0.004042 **
## Exp:bw 0.3626 1 0.547069
## Exp:sito 1.8662 1 0.171910
## bw:sito 0.1763 1 0.674577
## Exp:Movement 0.5018 1 0.478726
## bw:Movement 0.1241 1 0.724625
## sito:Movement 0.5317 1 0.465878
## Exp:bw:sito 1.9240 1 0.165418
## Exp:bw:Movement 0.0088 1 0.925112
## Exp:sito:Movement 0.1445 1 0.703851
## bw:sito:Movement 0.2968 1 0.585925
## Exp:bw:sito:Movement 0.0004 1 0.983567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residuals check
res = residuals(MT_model)
plot(density(res))
plot(MT_model)
qqnorm(residuals(MT_model))
qqline(residuals(MT_model))
plot dataset
pirateplot(formula = MovTime ~ bw + sito + Movement, #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Movement Time", xlab = "",#titles for x and y axes
)
plot main effect
pirateplot(formula = MovTime ~ Movement, #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Movement Time", xlab = "",#titles for x and y axes
)
mean and ds
library(dplyr)
summary_stats_asy <- dati_tms %>%
group_by(Movement) %>%
summarise(
mean_MovTime = mean(MovTime, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(MovTime, na.rm = TRUE) # Calcola la deviazione standard
)#Start - RT
RT_model = lmer(Start ~ Exp * bw * sito * Movement +
(1+sito*bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dati_tms)
summary(RT_model)## Linear mixed model fit by REML ['lmerMod']
## Formula: Start ~ Exp * bw * sito * Movement + (1 + sito * bw | Subject)
## Data: dati_tms
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 103746
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6279 -0.6037 -0.0732 0.5382 5.1626
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 10971 104.74
## sitoVertex 6497 80.60 -0.35
## bwWhite 2913 53.97 -0.25 0.48
## sitoVertex:bwWhite 5080 71.27 0.03 -0.40 -0.57
## Residual 14541 120.58
## Number of obs: 8324, groups: Subject, 40
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 395.3584 26.8647 14.717
## Exp3 -5.4322 34.6725 -0.157
## bwWhite -24.4025 15.9225 -1.533
## sitoVertex -6.5585 21.8341 -0.300
## MovementImitative -12.5495 8.3804 -1.497
## Exp3:bwWhite 0.2307 20.5363 0.011
## Exp3:sitoVertex -3.3633 28.2844 -0.119
## bwWhite:sitoVertex -17.8724 21.4117 -0.835
## Exp3:MovementImitative 16.9921 10.8140 1.571
## bwWhite:MovementImitative 14.8005 11.8656 1.247
## sitoVertex:MovementImitative -3.7189 11.7818 -0.316
## Exp3:bwWhite:sitoVertex 38.7627 27.7364 1.398
## Exp3:bwWhite:MovementImitative -18.7756 15.2894 -1.228
## Exp3:sitoVertex:MovementImitative 8.8443 15.3182 0.577
## bwWhite:sitoVertex:MovementImitative 2.7282 16.6904 0.163
## Exp3:bwWhite:sitoVertex:MovementImitative -6.8240 21.5863 -0.316##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itisSingular(RT_model)## [1] FALSEAnova(RT_model)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Start
## Chisq Df Pr(>Chisq)
## Exp 0.1131 1 0.736630
## bw 7.3784 1 0.006601 **
## sito 0.2118 1 0.645338
## Movement 0.0075 1 0.930968
## Exp:bw 0.1641 1 0.685370
## Exp:sito 0.5092 1 0.475499
## bw:sito 0.1305 1 0.717898
## Exp:Movement 3.5998 1 0.057786 .
## bw:Movement 0.3036 1 0.581622
## sito:Movement 0.0266 1 0.870349
## Exp:bw:sito 1.9170 1 0.166191
## Exp:bw:Movement 4.2304 1 0.039707 *
## Exp:sito:Movement 0.2511 1 0.616320
## bw:sito:Movement 0.0163 1 0.898407
## Exp:bw:sito:Movement 0.0999 1 0.751908
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1check residuals
res = residuals(RT_model)
plot(density(res))
plot(RT_model)
qqnorm(residuals(RT_model))
qqline(residuals(RT_model))
posthocs
#emmeans/posthoc
emmeans_interaction1 <- emmeans(RT_model, ~ Exp * Movement)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 8324' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 8324)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 8324' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 8324)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsemmeans_interaction2 <- emmeans(RT_model, ~ Movement * Exp * bw)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 8324' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 8324)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 8324' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 8324)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsposthoc_interaction1 <- pairs(emmeans_interaction1, adjust = "bonferroni")
posthoc_interaction2 <- pairs(emmeans_interaction2, adjust = "bonferroni")
summary(posthoc_interaction1)## contrast estimate SE df z.ratio p.value
## Exp1 Complementary - Exp3 Complementary -2.69 31.48 Inf -0.086 1.0000
## Exp1 Complementary - Exp1 Imitative 6.33 4.17 Inf 1.516 0.7767
## Exp1 Complementary - Exp3 Imitative -6.69 31.47 Inf -0.212 1.0000
## Exp3 Complementary - Exp1 Imitative 9.02 31.47 Inf 0.287 1.0000
## Exp3 Complementary - Exp3 Imitative -3.99 3.42 Inf -1.167 1.0000
## Exp1 Imitative - Exp3 Imitative -13.01 31.47 Inf -0.413 1.0000
##
## Results are averaged over the levels of: bw, sito
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 testssummary(posthoc_interaction2)## contrast estimate SE df z.ratio
## Complementary Exp1 Black - Imitative Exp1 Black 14.409 5.89 Inf 2.446
## Complementary Exp1 Black - Complementary Exp3 Black 7.114 32.21 Inf 0.221
## Complementary Exp1 Black - Imitative Exp3 Black 0.109 32.21 Inf 0.003
## Complementary Exp1 Black - Complementary Exp1 White 33.339 12.68 Inf 2.629
## Complementary Exp1 Black - Imitative Exp1 White 31.583 12.67 Inf 2.492
## Complementary Exp1 Black - Complementary Exp3 White 20.840 32.45 Inf 0.642
## Complementary Exp1 Black - Imitative Exp3 White 19.858 32.45 Inf 0.612
## Imitative Exp1 Black - Complementary Exp3 Black -7.295 32.20 Inf -0.227
## Imitative Exp1 Black - Imitative Exp3 Black -14.300 32.20 Inf -0.444
## Imitative Exp1 Black - Complementary Exp1 White 18.930 12.66 Inf 1.496
## Imitative Exp1 Black - Imitative Exp1 White 17.174 12.65 Inf 1.358
## Imitative Exp1 Black - Complementary Exp3 White 6.431 32.44 Inf 0.198
## Imitative Exp1 Black - Imitative Exp3 White 5.449 32.44 Inf 0.168
## Complementary Exp3 Black - Imitative Exp3 Black -7.005 4.89 Inf -1.431
## Complementary Exp3 Black - Complementary Exp1 White 26.225 32.60 Inf 0.804
## Complementary Exp3 Black - Imitative Exp1 White 24.469 32.60 Inf 0.751
## Complementary Exp3 Black - Complementary Exp3 White 13.727 10.42 Inf 1.318
## Complementary Exp3 Black - Imitative Exp3 White 12.744 10.41 Inf 1.224
## Imitative Exp3 Black - Complementary Exp1 White 33.230 32.60 Inf 1.019
## Imitative Exp3 Black - Imitative Exp1 White 31.475 32.60 Inf 0.965
## Imitative Exp3 Black - Complementary Exp3 White 20.732 10.41 Inf 1.991
## Imitative Exp3 Black - Imitative Exp3 White 19.750 10.41 Inf 1.898
## Complementary Exp1 White - Imitative Exp1 White -1.756 5.91 Inf -0.297
## Complementary Exp1 White - Complementary Exp3 White -12.498 32.84 Inf -0.381
## Complementary Exp1 White - Imitative Exp3 White -13.480 32.84 Inf -0.411
## Imitative Exp1 White - Complementary Exp3 White -10.743 32.84 Inf -0.327
## Imitative Exp1 White - Imitative Exp3 White -11.725 32.83 Inf -0.357
## Complementary Exp3 White - Imitative Exp3 White -0.982 4.78 Inf -0.205
## p.value
## 0.4045
## 1.0000
## 1.0000
## 0.2396
## 0.3553
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
##
## Results are averaged over the levels of: sito
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 28 testsmeans, ds
summary_statsRT <- dati_tms %>%
group_by(bw) %>%
summarise(
mean_start = mean(Start, na.rm = TRUE), # Calcola la media
sd_start = sd(Start, na.rm = TRUE) # Calcola la deviazione standard
)dataset plot
pirateplot(formula = Start ~ bw + sito + Movement, #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Times", xlab = "",#titles for x and y axes
)
main effect plot
pirateplot(formula = Start ~ bw , #dependent variable ~ independent
data = dati_tms, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Times", xlab = "",#titles for x and y axes
)
##STUDIO 2
EXP2 = read.xlsx("dati_tms_EXP2.xlsx", sheetIndex = 1)Factor Specification
EXP2$sito = as.factor(EXP2$sito)
EXP2$bw = as.factor(EXP2$bw)
EXP2$Grasp= as.factor(EXP2$Grasp)
EXP2$Movement = as.factor(EXP2$Movement)
EXP2$Subject = as.factor (EXP2$Subject)#Grasping Asynchrony
ASY_EXP2 = lmer (Asynchrony ~ bw * sito * Movement +
(1|Subject), control = lmerControl(optimizer = "bobyqa"),
data = EXP2)
isSingular(ASY_EXP2)## [1] FALSEsummary(ASY_EXP2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Asynchrony ~ bw * sito * Movement + (1 | Subject)
## Data: EXP2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 39010.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8506 -0.7425 -0.1798 0.5811 7.1245
##
## Random effects:
## Groups Name Variance Std.Dev.
## Subject (Intercept) 821.8 28.67
## Residual 13518.8 116.27
## Number of obs: 3160, groups: Subject, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 171.829 9.426 18.230
## bwWhite -17.747 8.295 -2.139
## sitovpm 1.320 8.238 0.160
## MovementImitative -28.946 8.270 -3.500
## bwWhite:sitovpm 12.457 11.721 1.063
## bwWhite:MovementImitative 15.445 11.710 1.319
## sitovpm:MovementImitative 10.838 11.673 0.928
## bwWhite:sitovpm:MovementImitative -16.648 16.549 -1.006
##
## Correlation of Fixed Effects:
## (Intr) bwWhit sitvpm MvmntI bwWht: bwW:MI stv:MI
## bwWhite -0.436
## sitovpm -0.439 0.498
## MovmntImttv -0.437 0.496 0.500
## bwWht:stvpm 0.308 -0.708 -0.703 -0.351
## bwWht:MvmnI 0.308 -0.708 -0.353 -0.706 0.501
## stvpm:MvmnI 0.310 -0.352 -0.706 -0.708 0.496 0.500
## bwWht:st:MI -0.218 0.501 0.498 0.500 -0.708 -0.708 -0.705Anova(ASY_EXP2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Asynchrony
## Chisq Df Pr(>Chisq)
## bw 3.6852 1 0.05490 .
## sito 4.4940 1 0.03401 *
## Movement 23.3007 1 0.000001385 ***
## bw:sito 0.2462 1 0.61975
## bw:Movement 0.7383 1 0.39020
## sito:Movement 0.0954 1 0.75740
## bw:sito:Movement 1.0120 1 0.31441
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1check residuals
res = residuals(ASY_EXP2)
plot(density(res))
plot(ASY_EXP2)
qqnorm(residuals(ASY_EXP2))
qqline(residuals(ASY_EXP2))
data plot
pirateplot(formula = Asynchrony ~ bw + sito + Movement, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "Experiment 2",#titles for x and y axes
)
Main effect Site
pirateplot(formula = Asynchrony ~ sito, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "Experiment 2",#titles for x and y axes
)
Main effect Movement
pirateplot(formula = Asynchrony ~ Movement, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "Experiment 2",#titles for x and y axes
)
mean and ds
statsGA_EXP2 <- EXP2 %>%
group_by(Movement) %>%
summarise(
mean_asy = mean(Asynchrony, na.rm = TRUE), # Calcola la media
sd_asy = sd(Asynchrony, na.rm = TRUE) # Calcola la deviazione standard
)
statsGA_EXP2_sito <- EXP2 %>%
group_by(sito) %>%
summarise(
mean_sito = mean(Asynchrony, na.rm = TRUE), # Calcola la media
sd_sito = sd(Asynchrony, na.rm = TRUE) # Calcola la deviazione standard
)#Since data were not normal: Alternative - GA log
EXP2$Asynchrony_log = log(EXP2$Asynchrony + 225)run LMM log
ASY_model_studio2 = lmer(Asynchrony_log ~ bw * Movement * sito +
(1|Subject), control = lmerControl(optimizer = "bobyqa"),
data = EXP2)
#summary(ASY_model_studio2)
isSingular(ASY_model_studio2)## [1] FALSEAnova(ASY_model_studio2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Asynchrony_log
## Chisq Df Pr(>Chisq)
## bw 4.4719 1 0.03446 *
## Movement 22.5408 1 0.000002057 ***
## sito 2.5316 1 0.11159
## bw:Movement 0.1595 1 0.68962
## bw:sito 0.2306 1 0.63108
## Movement:sito 0.1169 1 0.73240
## bw:Movement:sito 1.6254 1 0.20234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residuals check
res = residuals(ASY_model_studio2)
plot(density(res))
plot(ASY_model_studio2)
qqnorm(residuals(ASY_model_studio2))
qqline(residuals(ASY_model_studio2))
data plots
pirateplot(formula = Asynchrony_log ~ bw, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
pirateplot(formula = Asynchrony_log ~ Movement, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
Means, ds
log_statsASY2 <- EXP2 %>%
group_by(bw) %>%
summarise(
mean_bw = mean(Asynchrony_log, na.rm = TRUE), # Calcola la media
sd_bw = sd(Asynchrony_log, na.rm = TRUE) # Calcola la deviazione standard
)
log_statsASY_m2 <- EXP2 %>%
group_by(Movement) %>%
summarise(
mean_m = mean(Asynchrony_log, na.rm = TRUE), # Calcola la media
sd_m = sd(Asynchrony_log, na.rm = TRUE) # Calcola la deviazione standard
)#Movement Time
MT_EXP2 = lmer(MovTime ~ bw * sito * Movement +
(1+sito*bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = EXP2)
summary(MT_EXP2)## Linear mixed model fit by REML ['lmerMod']
## Formula: MovTime ~ bw * sito * Movement + (1 + sito * bw | Subject)
## Data: EXP2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 41817.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6970 -0.6989 -0.0768 0.6399 4.5265
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 29141 170.71
## sitovpm 13913 117.95 -0.68
## bwWhite 1661 40.75 -0.67 0.56
## sitovpm:bwWhite 6872 82.90 0.67 -0.94 -0.44
## Residual 32035 178.98
## Number of obs: 3160, groups: Subject, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1304.871 44.992 29.002
## bwWhite 9.786 16.567 0.591
## sitovpm 5.373 33.007 0.163
## MovementImitative 7.162 12.766 0.561
## bwWhite:sitovpm -28.003 28.021 -0.999
## bwWhite:MovementImitative -10.326 18.042 -0.572
## sitovpm:MovementImitative -12.379 17.986 -0.688
## bwWhite:sitovpm:MovementImitative 23.157 25.484 0.909
##
## Correlation of Fixed Effects:
## (Intr) bwWhit sitvpm MvmntI bwWht: bwW:MI stv:MI
## bwWhite -0.523
## sitovpm -0.666 0.479
## MovmntImttv -0.142 0.386 0.194
## bwWht:stvpm 0.568 -0.565 -0.837 -0.228
## bwWht:MvmnI 0.100 -0.548 -0.137 -0.707 0.324
## stvpm:MvmnI 0.101 -0.274 -0.273 -0.710 0.321 0.502
## bwWht:st:MI -0.071 0.388 0.192 0.501 -0.458 -0.708 -0.706isSingular(MT_EXP2)## [1] FALSEAnova(MT_EXP2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MovTime
## Chisq Df Pr(>Chisq)
## bw 0.0001 1 0.9903
## sito 1.5579 1 0.2120
## Movement 0.0604 1 0.8059
## bw:sito 0.4304 1 0.5118
## bw:Movement 0.0101 1 0.9200
## sito:Movement 0.0044 1 0.9472
## bw:sito:Movement 0.8257 1 0.3635Posthocs
#emmeans/posthoc
emmeans_interaction <- emmeans(MT_EXP2, ~ bw * sito)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3160' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3160)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3160' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3160)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsemmeans_interaction1 <- emmeans(MT_EXP2, ~ bw * Movement)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3160' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3160)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3160' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3160)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsemmeans_interaction2 <- emmeans(MT_EXP2, ~ Movement * sito)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3160' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3160)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3160' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3160)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsposthoc_interaction <- pairs(emmeans_interaction, adjust = "bonferroni")
posthoc_interaction1 <- pairs(emmeans_interaction1, adjust = "bonferroni")
posthoc_interaction2 <- pairs(emmeans_interaction2, adjust = "bonferroni")
summary(posthoc_interaction)## contrast estimate SE df z.ratio p.value
## Black v - White v -4.623 13.9 Inf -0.334 1.0000
## Black v - Black vpm 0.816 31.8 Inf 0.026 1.0000
## Black v - White vpm 12.618 22.7 Inf 0.555 1.0000
## White v - Black vpm 5.439 27.5 Inf 0.197 1.0000
## White v - White vpm 17.241 15.6 Inf 1.107 1.0000
## Black vpm - White vpm 11.802 21.3 Inf 0.555 1.0000
##
## Results are averaged over the levels of: Movement
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 testssummary(posthoc_interaction1)## contrast estimate SE df z.ratio p.value
## Black Complementary - White Complementary 4.216 14.44 Inf 0.292 1.0000
## Black Complementary - Black Imitative -0.973 8.99 Inf -0.108 1.0000
## Black Complementary - White Imitative 1.991 14.41 Inf 0.138 1.0000
## White Complementary - Black Imitative -5.188 14.44 Inf -0.359 1.0000
## White Complementary - White Imitative -2.225 9.03 Inf -0.246 1.0000
## Black Imitative - White Imitative 2.964 14.40 Inf 0.206 1.0000
##
## Results are averaged over the levels of: sito
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 testssummary(posthoc_interaction2)## contrast estimate SE df z.ratio p.value
## Complementary v - Imitative v -2.00 9.02 Inf -0.222 1.0000
## Complementary v - Complementary vpm 8.63 22.62 Inf 0.382 1.0000
## Complementary v - Imitative vpm 7.43 22.60 Inf 0.329 1.0000
## Imitative v - Complementary vpm 10.63 22.60 Inf 0.470 1.0000
## Imitative v - Imitative vpm 9.43 22.59 Inf 0.417 1.0000
## Complementary vpm - Imitative vpm -1.20 9.00 Inf -0.133 1.0000
##
## Results are averaged over the levels of: bw
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 testsControl residuals distributions
res = residuals(MT_EXP2)
plot(density(res))
plot(MT_EXP2)
qqnorm(residuals(MT_EXP2))
qqline(residuals(MT_EXP2))
data plots
pirateplot(formula = MovTime ~ bw + sito + Movement, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Movement Time", xlab = "Experiment 2",#titles for x and y axes
)
#Start - RT
RT_EXP2 = lmer(Start ~ bw * sito * Movement +
(1+sito*bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = EXP2)
summary(RT_EXP2)## Linear mixed model fit by REML ['lmerMod']
## Formula: Start ~ bw * sito * Movement + (1 + sito * bw | Subject)
## Data: EXP2
## Control: lmerControl(optimizer = "bobyqa")
##
## REML criterion at convergence: 37649.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8586 -0.6389 -0.0617 0.5418 4.9565
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 12071.8 109.87
## sitovpm 9977.2 99.89 -0.69
## bwWhite 763.7 27.64 -0.23 0.55
## sitovpm:bwWhite 6082.0 77.99 0.29 -0.71 -0.83
## Residual 10081.8 100.41
## Number of obs: 3114, groups: Subject, 15
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 274.408 28.821 9.521
## bwWhite -7.326 10.163 -0.721
## sitovpm 21.584 26.771 0.806
## MovementImitative -18.340 7.197 -2.548
## bwWhite:sitovpm 16.362 22.583 0.725
## bwWhite:MovementImitative 10.445 10.200 1.024
## sitovpm:MovementImitative 5.852 10.153 0.576
## bwWhite:sitovpm:MovementImitative -27.739 14.406 -1.926
##
## Correlation of Fixed Effects:
## (Intr) bwWhit sitvpm MvmntI bwWht: bwW:MI stv:MI
## bwWhite -0.244
## sitovpm -0.687 0.467
## MovmntImttv -0.125 0.354 0.134
## bwWht:stvpm 0.293 -0.749 -0.695 -0.159
## bwWht:MvmnI 0.088 -0.505 -0.095 -0.706 0.227
## stvpm:MvmnI 0.088 -0.251 -0.190 -0.709 0.225 0.500
## bwWht:st:MI -0.062 0.358 0.134 0.500 -0.321 -0.708 -0.705isSingular(RT_EXP2)## [1] FALSEAnova(RT_EXP2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Start
## Chisq Df Pr(>Chisq)
## bw 0.0558 1 0.81327
## sito 2.0054 1 0.15674
## Movement 22.5716 1 0.000002025 ***
## bw:sito 0.0125 1 0.91104
## bw:Movement 0.2308 1 0.63092
## sito:Movement 1.2110 1 0.27114
## bw:sito:Movement 3.7078 1 0.05416 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Control residual distribution
res = residuals(RT_EXP2)
plot(density(res))
plot(RT_EXP2)
qqnorm(residuals(RT_EXP2))
qqline(residuals(RT_EXP2))
Posthocs
#emmeans/posthoc
emmeans_interaction <- emmeans(RT_EXP2, ~ bw * sito * Movement)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 3114' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 3114)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 3114' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 3114)' or larger];
## but be warned that this may result in large computation time and memory use.posthoc_interaction <- pairs(emmeans_interaction, adjust = "bonferroni")
summary(posthoc_interaction)## contrast estimate SE df z.ratio
## Black v Complementary - White v Complementary 7.326 10.16 Inf 0.721
## Black v Complementary - Black vpm Complementary -21.584 26.77 Inf -0.806
## Black v Complementary - White vpm Complementary -30.620 19.99 Inf -1.532
## Black v Complementary - Black v Imitative 18.340 7.20 Inf 2.548
## Black v Complementary - White v Imitative 15.220 10.13 Inf 1.503
## Black v Complementary - Black vpm Imitative -9.096 26.77 Inf -0.340
## Black v Complementary - White vpm Imitative -0.839 19.97 Inf -0.042
## White v Complementary - Black vpm Complementary -28.909 23.79 Inf -1.215
## White v Complementary - White vpm Complementary -37.946 19.66 Inf -1.931
## White v Complementary - Black v Imitative 11.014 10.16 Inf 1.084
## White v Complementary - White v Imitative 7.895 7.23 Inf 1.092
## White v Complementary - Black vpm Imitative -16.422 23.79 Inf -0.690
## White v Complementary - White vpm Imitative -8.164 19.64 Inf -0.416
## Black vpm Complementary - White vpm Complementary -9.037 16.41 Inf -0.551
## Black vpm Complementary - Black v Imitative 39.924 26.77 Inf 1.491
## Black vpm Complementary - White v Imitative 36.804 23.77 Inf 1.548
## Black vpm Complementary - Black vpm Imitative 12.488 7.16 Inf 1.744
## Black vpm Complementary - White vpm Imitative 20.745 16.39 Inf 1.266
## White vpm Complementary - Black v Imitative 48.960 19.99 Inf 2.449
## White vpm Complementary - White v Imitative 45.841 19.64 Inf 2.335
## White vpm Complementary - Black vpm Imitative 21.524 16.41 Inf 1.312
## White vpm Complementary - White vpm Imitative 29.782 7.22 Inf 4.122
## Black v Imitative - White v Imitative -3.120 10.13 Inf -0.308
## Black v Imitative - Black vpm Imitative -27.436 26.77 Inf -1.025
## Black v Imitative - White vpm Imitative -19.179 19.97 Inf -0.960
## White v Imitative - Black vpm Imitative -24.316 23.77 Inf -1.023
## White v Imitative - White vpm Imitative -16.059 19.62 Inf -0.819
## Black vpm Imitative - White vpm Imitative 8.257 16.39 Inf 0.504
## p.value
## 1.0000
## 1.0000
## 1.0000
## 0.3031
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 0.4011
## 0.5479
## 1.0000
## 0.0011
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 28 testsMean and ds
summary_statsRT_EXP2 <- EXP2 %>%
group_by(Movement) %>%
summarise(
mean_MovTime = mean(Start, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(Start, na.rm = TRUE) # Calcola la deviazione standard
)Dataplots
pirateplot(formula = Start ~ bw + sito + Movement, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Time", xlab = "Experiment 2",#titles for x and y axes
)
Main Effect Plot
pirateplot(formula = Start ~ Movement, #dependent variable ~ independent
data = EXP2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Time", xlab = "Experiment 2",#titles for x and y axes
)
##ALL STUDIES
dataset_2 = read.xlsx("dataset_2.xlsx", sheetIndex = 1)factors declaration
dataset_2$sito = as.factor(dataset_2$sito)
dataset_2$bw = as.factor(dataset_2$bw)
dataset_2$Movement = as.factor(dataset_2$Movement)
dataset_2$Grasp = as.factor(dataset_2$Grasp)
dataset_2$Subject = as.factor(dataset_2$Subject)#Grasping Asynchrony
ASY_model_3studi = lmer(Asynchrony ~ Exp * bw * Movement + sito +
(1+bw+sito|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dataset_2)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?#summary(ASY_model_3studi)
isSingular(ASY_model_3studi)## [1] FALSEAnova(ASY_model_3studi)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Asynchrony
## Chisq Df Pr(>Chisq)
## Exp 2.2764 2 0.32039
## bw 0.0514 1 0.82066
## Movement 46.1702 1 0.00000000001084 ***
## sito 0.5868 2 0.74573
## Exp:bw 2.5650 2 0.27734
## Exp:Movement 5.8548 2 0.05354 .
## bw:Movement 0.2921 1 0.58890
## Exp:bw:Movement 0.6230 2 0.73234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residuals check
res = residuals(ASY_model_3studi)
plot(density(res))
plot(ASY_model_3studi)
qqnorm(residuals(ASY_model_3studi))
qqline(residuals(ASY_model_3studi))
means, ds
summary_statsASY_ALL <- dataset_2 %>%
group_by(Movement) %>%
summarise(
mean_Mov = mean(Asynchrony, na.rm = TRUE), # Calcola la media
sd_Mov = sd(Asynchrony, na.rm = TRUE) # Calcola la deviazione standard
)data plots
pirateplot(formula = Asynchrony ~ Movement, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
#Since data were not normal: Alternative - GA log
dataset_2$Asynchrony_log = log(dataset_2$Asynchrony + 225)run LMM log
ASY_model_3studi2 = lmer(Asynchrony_log ~ Exp * bw * Movement + sito +
(1+bw+sito|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dataset_2)
#summary(ASY_model_3studi2)
isSingular(ASY_model_3studi2)## [1] FALSEAnova(ASY_model_3studi2)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Asynchrony_log
## Chisq Df Pr(>Chisq)
## Exp 2.7353 2 0.2547
## bw 0.1667 1 0.6830
## Movement 42.4888 1 0.00000000007109 ***
## sito 0.3596 2 0.8354
## Exp:bw 3.0046 2 0.2226
## Exp:Movement 4.4030 2 0.1106
## bw:Movement 0.0390 1 0.8435
## Exp:bw:Movement 0.2080 2 0.9012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1residuals check
res = residuals(ASY_model_3studi2)
plot(density(res))
plot(ASY_model_3studi2)
qqnorm(residuals(ASY_model_3studi2))
qqline(residuals(ASY_model_3studi2))
data plots
pirateplot(formula = Asynchrony_log ~ Movement, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Grasping Asynchrony", xlab = "",#titles for x and y axes
)
means, ds
log_statsASY_ALL <- dataset_2 %>%
group_by(Movement) %>%
summarise(
mean_Mov = mean(Asynchrony_log, na.rm = TRUE), # Calcola la media
sd_Mov = sd(Asynchrony_log, na.rm = TRUE) # Calcola la deviazione standard
)#Movement Time
MT_model_3studi = lmer(MovTime ~ Exp * bw * Movement +
(1+bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dataset_2)
#summary(MT_model_3studi)
isSingular(MT_model_3studi)## [1] FALSEAnova(MT_model_3studi)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: MovTime
## Chisq Df Pr(>Chisq)
## Exp 8.1235 2 0.01722 *
## bw 0.1436 1 0.70470
## Movement 5.3885 1 0.02027 *
## Exp:bw 0.5771 2 0.74936
## Exp:Movement 3.1673 2 0.20522
## bw:Movement 0.1746 1 0.67608
## Exp:bw:Movement 0.0105 2 0.99476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1check residuals
res = residuals(MT_model_3studi)
plot(density(res))
plot(MT_model_3studi)
qqnorm(residuals(MT_model_3studi))
qqline(residuals(MT_model_3studi))
means, ds
summary_statsMT_Exp <- dataset_2 %>%
group_by(Exp) %>%
summarise(
mean_MovTime = mean(MovTime, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(MovTime, na.rm = TRUE) # Calcola la deviazione standard
)
summary_statsMT <- dataset_2 %>%
group_by(Movement) %>%
summarise(
mean_MovTime = mean(MovTime, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(MovTime, na.rm = TRUE) # Calcola la deviazione standard
)posthocs
#emmeans/posthoc
emmeans_Experiment <- emmeans(MT_model_3studi, ~ Exp)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 11575' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 11575)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 11575' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 11575)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsposthoc_Experiment <- contrast(emmeans_Experiment, adjust = "bonferroni")
summary(posthoc_Experiment)## contrast estimate SE df z.ratio p.value
## Exp1 effect 9.6 24.6 Inf 0.389 1.0000
## Exp2 effect 44.6 25.1 Inf 1.778 0.2263
## Exp3 effect -54.2 22.4 Inf -2.419 0.0467
##
## Results are averaged over the levels of: bw, Movement
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 3 testsmain effects plots
pirateplot(formula = MovTime ~ Movement, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Movement Time", xlab = "",#titles for x and y axes
)
pirateplot(formula = MovTime ~ Exp, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Movement Time", xlab = "",#titles for x and y axes
)
#Start - RT
RT_model_3studi = lmer(Start ~ Exp * bw * Movement +
(1+bw|Subject), control = lmerControl(optimizer = "bobyqa"),
data = dataset_2)
#summary(RT_model_3studi)
isSingular(RT_model_3studi)## [1] FALSEAnova(RT_model_3studi)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Start
## Chisq Df Pr(>Chisq)
## Exp 12.9470 2 0.001544 **
## bw 7.2134 1 0.007236 **
## Movement 3.8781 1 0.048919 *
## Exp:bw 2.5901 2 0.273891
## Exp:Movement 13.0285 2 0.001482 **
## bw:Movement 0.0809 1 0.776099
## Exp:bw:Movement 4.2210 2 0.121176
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1check residuals
res = residuals(RT_model_3studi)
plot(density(res))
plot(RT_model_3studi)
qqnorm(residuals(RT_model_3studi))
qqline(residuals(RT_model_3studi))
posthocs
#emmeans/posthoc
RT_emmeans_interaction <- emmeans(RT_model_3studi, ~ Exp * Movement)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 11438' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 11438)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 11438' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 11438)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsRT_emmeans_Experiment <- emmeans(MT_model_3studi, ~ Exp)## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'pbkrtest.limit = 11575' (or larger)
## [or, globally, 'set emm_options(pbkrtest.limit = 11575)' or larger];
## but be warned that this may result in large computation time and memory use.## Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
## To enable adjustments, add the argument 'lmerTest.limit = 11575' (or larger)
## [or, globally, 'set emm_options(lmerTest.limit = 11575)' or larger];
## but be warned that this may result in large computation time and memory use.## NOTE: Results may be misleading due to involvement in interactionsposthoc_interaction <- pairs(RT_emmeans_interaction, adjust = "bonferroni")
posthoc_Experiment <- contrast(RT_emmeans_Experiment, adjust = "bonferroni")
summary(posthoc_interaction)## contrast estimate SE df z.ratio p.value
## Exp1 Complementary - Exp2 Complementary 90.34 33.90 Inf 2.664 0.1157
## Exp1 Complementary - Exp3 Complementary -2.67 30.45 Inf -0.088 1.0000
## Exp1 Complementary - Exp1 Imitative 6.39 4.24 Inf 1.506 1.0000
## Exp1 Complementary - Exp2 Imitative 106.42 33.90 Inf 3.139 0.0254
## Exp1 Complementary - Exp3 Imitative -6.71 30.45 Inf -0.220 1.0000
## Exp2 Complementary - Exp3 Complementary -93.01 31.05 Inf -2.995 0.0411
## Exp2 Complementary - Exp1 Imitative -83.94 33.90 Inf -2.476 0.1992
## Exp2 Complementary - Exp2 Imitative 16.08 4.40 Inf 3.655 0.0039
## Exp2 Complementary - Exp3 Imitative -97.05 31.05 Inf -3.126 0.0266
## Exp3 Complementary - Exp1 Imitative 9.06 30.44 Inf 0.298 1.0000
## Exp3 Complementary - Exp2 Imitative 109.09 31.05 Inf 3.514 0.0066
## Exp3 Complementary - Exp3 Imitative -4.04 3.48 Inf -1.161 1.0000
## Exp1 Imitative - Exp2 Imitative 100.03 33.90 Inf 2.951 0.0476
## Exp1 Imitative - Exp3 Imitative -13.11 30.44 Inf -0.431 1.0000
## Exp2 Imitative - Exp3 Imitative -113.13 31.05 Inf -3.644 0.0040
##
## Results are averaged over the levels of: bw
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 15 testssummary(posthoc_Experiment)## contrast estimate SE df z.ratio p.value
## Exp1 effect 9.6 24.6 Inf 0.389 1.0000
## Exp2 effect 44.6 25.1 Inf 1.778 0.2263
## Exp3 effect -54.2 22.4 Inf -2.419 0.0467
##
## Results are averaged over the levels of: bw, Movement
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 3 testsmeans and ds
statsRT_Exp <- dataset_2 %>%
group_by(Exp) %>%
summarise(
mean_MovTime = mean(Start, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(Start, na.rm = TRUE) # Calcola la deviazione standard
)
statsRT_Mov <- dataset_2 %>%
group_by(Movement) %>%
summarise(
mean_MovTime = mean(Start, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(Start, na.rm = TRUE) # Calcola la deviazione standard
)
statsRT_bw <- dataset_2 %>%
group_by(bw) %>%
summarise(
mean_MovTime = mean(Start, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(Start, na.rm = TRUE) # Calcola la deviazione standard
)
statsRT_interaction <- dataset_2 %>%
group_by(Exp, Movement) %>%
summarise(
mean_MovTime = mean(Start, na.rm = TRUE), # Calcola la media
sd_MovTime = sd(Start, na.rm = TRUE) # Calcola la deviazione standard
)## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.significant results plots
pirateplot(formula = Start ~ Movement, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Times", xlab = "",#titles for x and y axes
)
pirateplot(formula = Start ~ Exp, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Times", xlab = "",#titles for x and y axes
)
pirateplot(formula = Start ~ bw, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Times", xlab = "",#titles for x and y axes
)
pirateplot(formula = Start ~ Exp + Movement, #dependent variable ~ independent
data = dataset_2, #data
theme = 3, #set theme
gl.col = "white",
ylab = "Reaction Times", xlab = "",#titles for x and y axes
)