Contact Experiment Data Analysis

First Step is to load the necessary package, If you dont have them just install them. For jjstatsplot you need to install it remotely. Just remove the dash and press enter. Then Press 3 (none package to be updated). ?group_by()

#remotes::install_github("sbalci/jjstatsplot") #Press 3 !!!!! i.e., installing/Updating none package! 
library(jmv)
library(datasets)
library(plyr)
library(readr)
library(dataframes2xls)
library(data.table)
library(plyr)
library(ggstatsplot)
library(jjstatsplot)
library(lme4)
library(lmerTest)
library(ggplot2)
library(rstatix)
library(coin)
library(ARTool)
library(ggpubr)
library(tidyverse)
library(dplyr)
library("afex")     
library("emmeans")  
library("multcomp") 

Import and merge all the csv files (x.csv where x = ID) of the folder. Note, that the ID is has already been added as column

Discard first X trials per interpenetration feedback condition and then create a summary table for each participant. You need to define nrTrialsPerBlockToRemove.

data$Part <- as.factor(data$Block < 4)
levels(data$Part)

data$Part <- factor(data$Part,levels = c("TRUE","FALSE"),
                  labels = c("Part 1","Part 2"))

nrTrialsPerBlockToRemove <- 1
#trialsToRemove <- seq(from = 1, to = nrTrialsPerBlockToRemove)

  data <-  data %>%
    group_by(ID, Part, InterpenetrationFeedback, FullyShaded) %>% # I have added here the fully shaded 
    slice(nrTrialsPerBlockToRemove+1:n())
  # to double check we are discarding the right rows
  #print(data[[i]]$Trial)

 # This df will be used to create the subsets for 1st part and 2nd part of the experiment. 
  
data$InterpenetrationFeedback  <- as.factor(data$InterpenetrationFeedback)
data$FullyShaded <- as.factor(data$FullyShaded)

Let’s visualize the data per interpenetration feedback and/or part of the experiment (part 1 & part 2).

Warning messages:
1: Unknown or uninitialised column: `Extreme1`. 
2: Unknown or uninitialised column: `Extreme2`. 

ParsiDFplots <- data
ParsiDFplots$ID[ParsiDFplots$ID == 9] <- NA
ParsiDFplots$ID[ParsiDFplots$ID == 17] <- NA
ParsiDFplots$ID[ParsiDFplots$ID == 20] <- NA
ParsiDFplots <- na.omit(ParsiDFplots)
ParsiDFplots <- aggregate(. ~ ID + Age + Gender + InterpenetrationFeedback + Part, ParsiDFplots, mean)
ParsiDFplots$AverageInterpenetration <-100 * ParsiDFplots$AverageInterpenetration
ParsiDFplots$MaxInterpenetration <- 100 * ParsiDFplots$MaxInterpenetration

p1 <- ggstatsplot::ggbetweenstats(
  data = ParsiDFplots,
  x = "InterpenetrationFeedback", #Indepedent Variable
  y = "MaxInterpenetration", # Depedent Variable
  grouping.var = "Part", # 2nd IV 
  type = "p", # parametric test i.e., p values
  pairwise.comparisons = FALSE, #compute pairwise comparisons
  pairwise.display = "significant", # show only the significant ones
  p.adjust.method = "bonferroni", # correction of p-value
  effsize.type = "unbiased", # Calculates the Hedge's g for t tests and the partial Omega for ANOVA
  results.subtitle = FALSE,
  xlab = "Type of Feedback", #label of X axis
  ylab = "Maximum Interpenetration", #label of y axis
  sample.size.label = FALSE,
  var.equal = TRUE, #Assuming Equal variances
  mean.plotting = FALSE,
  mean.ci = TRUE, #display the confidence interval of the mean
  paired = TRUE, #indicating that we have a within subject design
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
  ) 

p2 <- ggstatsplot::ggbetweenstats(
  data = ParsiDFplots,
  x = "InterpenetrationFeedback",
  y = "AverageInterpenetration",
  grouping.var = "Part",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Type of Feedback",
  ylab = "Average Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
)


# Replicating the above but this time we look on the effect of the type of feedback on the DVs in 1st and 2nd Part of the experiment individually
p3 <- ggstatsplot::grouped_ggbetweenstats(
  data = ParsiDFplots,
  x = "InterpenetrationFeedback",
  y = "MaxInterpenetration",
  grouping.var = "Part",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Type of Feedback",
  ylab = "Maximum Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
  ) 

p4 <- ggstatsplot::grouped_ggbetweenstats(
  data = ParsiDFplots,
  x = "InterpenetrationFeedback",
  y = "AverageInterpenetration",
  grouping.var = "Part",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Type of Feedback",
  ylab = "Average Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
) 

# Lets check the effect of shaded condition on DVs
p5 <- ggstatsplot:: ggbetweenstats(
  data = ParsiDFplots,
  x = "Part",
  y = "MaxInterpenetration",
  grouping.var = "InterpenetrationFeedback",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Order",
  ylab = "Maximum Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
  ) 

p6 <- ggstatsplot::ggbetweenstats(
  data = ParsiDFplots,
  x = "Part",
  y = "AverageInterpenetration",
  grouping.var = "InterpenetrationFeedback",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Order",
  ylab = "Average Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
) 

p7 <- ggstatsplot::grouped_ggbetweenstats(
  data = ParsiDFplots,
  x = "Part",
  y = "MaxInterpenetration",
  grouping.var = "InterpenetrationFeedback",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Order",
  ylab = "Maximum Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black") 

p8 <- ggstatsplot::grouped_ggbetweenstats(
  data = ParsiDFplots,
  x = "Part",
  y = "AverageInterpenetration",
  grouping.var = "InterpenetrationFeedback",
  type = "p",
  pairwise.comparisons = FALSE,
  pairwise.display = "significant",
  p.adjust.method = "bonferroni",
  effsize.type = "unbiased",
  results.subtitle = FALSE,
  xlab = "Order",
  ylab = "Average Interpenetration",
  sample.size.label = FALSE,
  var.equal = TRUE,
  mean.plotting = FALSE,
  mean.ci = TRUE,
  paired = TRUE,
  title.text = "Interpenetration Box-Violin Plots",
  caption.text = "Note: Interpenetration distance is displayed in cm.",
  title.color = "black",
  caption.color = "black"
  ) 
p1

p2

p3 

p4

p5 

p6

p7

p8

Let’s check the Two Way Repeated Measures ANOVA

aMax <- aov_ez("ID", "MaxInterpenetration", ParsiDF,

Warning messages:
1: Unknown or uninitialised column: `Extreme1`. 
2: Unknown or uninitialised column: `Extreme2`. 
3: Unknown or uninitialised column: `Extreme1`. 
4: Unknown or uninitialised column: `Extreme2`. 
5: Unknown or uninitialised column: `Extreme1`. 
6: Unknown or uninitialised column: `Extreme2`. 
7: Unknown or uninitialised column: `Extreme1`. 
8: Unknown or uninitialised column: `Extreme2`. 
9: Unknown or uninitialised column: `Extreme1`. 
10: Unknown or uninitialised column: `Extreme2`. 

             within = c("Part", "InterpenetrationFeedback"),
             anova_table = list(es = "pes"))

knitr::kable(nice(aMax$anova_table))

Effect	df	MSE	F	pes	p.value
Part	1, 20	0.06	31.62 ***	.613	<.001
InterpenetrationFeedback	1.91, 38.11	0.12	29.77 ***	.598	<.001
Part:InterpenetrationFeedback	2.12, 42.40	0.04	5.80 **	.225	.005

aAv <- aov_ez("ID", "AverageInterpenetration", ParsiDF,
             within = c("Part", "InterpenetrationFeedback"),
             anova_table = list(es = "pes"))

knitr::kable(nice(aAv$anova_table))

Effect	df	MSE	F	pes	p.value
Part	1, 20	0.07	34.72 ***	.634	<.001
InterpenetrationFeedback	1.98, 39.53	0.12	27.65 ***	.580	<.001
Part:InterpenetrationFeedback	2.17, 43.30	0.04	6.63 **	.249	.002

effectsize::omega_squared(aMax, partial = TRUE, ci = 0.95)

Parameter                     | Omega2 (partial) |        95% CI
----------------------------------------------------------------
Part                          |             0.58 | [ 0.26, 0.75]
InterpenetrationFeedback      |             0.57 | [ 0.40, 0.69]
Part:InterpenetrationFeedback |             0.18 | [ 0.00, 0.34]

effectsize::omega_squared(aAv, partial = TRUE, ci = 0.95)

Parameter                     | Omega2 (partial) |       95% CI
---------------------------------------------------------------
Part                          |             0.61 | [0.29, 0.77]
InterpenetrationFeedback      |             0.56 | [0.37, 0.67]
Part:InterpenetrationFeedback |             0.21 | [0.01, 0.37]

We can see that every type of feedback as well as the interrelationship with the part of the experiment have a large effect on DVs!!!!!!

Reference for interpreting Omega Squared Small effect: ω2 = 0.01; Medium effect: ω2 = 0.06; Large effect: ω2 = 0.14.

Let’s plot the main effects (Interpenetration Feedback OR Part of the experiment).

aMaxPlots <- aov_ez("ID", "MaxInterpenetration", ParsiDFplots,
             within = c("Part", "InterpenetrationFeedback"),
             anova_table = list(es = "pes"))

aAvPlots <- aov_ez("ID", "AverageInterpenetration", ParsiDFplots,
             within = c("Part", "InterpenetrationFeedback"),
             anova_table = list(es = "pes"))


afex_plot(aMaxPlots, x = "InterpenetrationFeedback", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5)) +
            ylim(0, 3)

NOTE: Results may be misleading due to involvement in interactions

afex_plot(aMaxPlots, x = "Part", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5))  +
            ylim(0, 3)

NOTE: Results may be misleading due to involvement in interactions

afex_plot(aAvPlots, x = "InterpenetrationFeedback", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5)) +
            ylim(0, 2.25)

NOTE: Results may be misleading due to involvement in interactions

afex_plot(aAvPlots, x = "Part", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5))  +
            ylim(0, 2.25)

NOTE: Results may be misleading due to involvement in interactions

?afex_plot()

Let’s plot the main interaction effects (Interpenetration Feedback AND Part of the experiment).

afex_plot(aMaxPlots, x = "InterpenetrationFeedback", trace = "Part", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5)) +
            ylim(0, 3)

afex_plot(aMaxPlots, x = "Part", trace = "InterpenetrationFeedback", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5))  +
            ylim(0, 3)

afex_plot(aAvPlots, x = "InterpenetrationFeedback",  trace = "Part", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol::geom_boxjitter, 
                data_arg = list(width = 0.5)) +
            ylim(0, 2.25)

afex_plot(aAvPlots, x = "Part", trace = "InterpenetrationFeedback", error = "within", 
                mapping = c("linetype", "shape", "fill"),
                data_geom = ggpol:: geom_boxjitter, 
                data_arg = list(width = 0.5))  +
            ylim(0, 2.25)

Post-hoc Tests

# Main Effects
lsmMaxFeedback <- lsmeans(aMax,~ InterpenetrationFeedback)

NOTE: Results may be misleading due to involvement in interactions

contrast(lsmMax, method="pairwise", interaction = TRUE)

 InterpenetrationFeedback_pairwise estimate     SE df t.ratio p.value
 Both - Electrotactile              -0.2308 0.0597 60 -3.865  0.0003 
 Both - NoFeedback                  -0.5275 0.0597 60 -8.835  <.0001 
 Both - Visual                      -0.0921 0.0597 60 -1.543  0.1281 
 Electrotactile - NoFeedback        -0.2967 0.0597 60 -4.969  <.0001 
 Electrotactile - Visual             0.1387 0.0597 60  2.322  0.0236 
 NoFeedback - Visual                 0.4354 0.0597 60  7.292  <.0001 

Results are averaged over the levels of: Part 

lsmMaxPart <- lsmeans(aMax,~ Part)

NOTE: Results may be misleading due to involvement in interactions

contrast(lsmMax, method="pairwise", interaction = TRUE)

 InterpenetrationFeedback_pairwise estimate     SE df t.ratio p.value
 Both - Electrotactile              -0.2308 0.0597 60 -3.865  0.0003 
 Both - NoFeedback                  -0.5275 0.0597 60 -8.835  <.0001 
 Both - Visual                      -0.0921 0.0597 60 -1.543  0.1281 
 Electrotactile - NoFeedback        -0.2967 0.0597 60 -4.969  <.0001 
 Electrotactile - Visual             0.1387 0.0597 60  2.322  0.0236 
 NoFeedback - Visual                 0.4354 0.0597 60  7.292  <.0001 

Results are averaged over the levels of: Part 

lsmMaxIntEff <- lsmeans(aMax,~ InterpenetrationFeedback:Part)
contrast(lsmMax, method="pairwise", interaction = TRUE)

 InterpenetrationFeedback_pairwise estimate     SE df t.ratio p.value
 Both - Electrotactile              -0.2308 0.0597 60 -3.865  0.0003 
 Both - NoFeedback                  -0.5275 0.0597 60 -8.835  <.0001 
 Both - Visual                      -0.0921 0.0597 60 -1.543  0.1281 
 Electrotactile - NoFeedback        -0.2967 0.0597 60 -4.969  <.0001 
 Electrotactile - Visual             0.1387 0.0597 60  2.322  0.0236 
 NoFeedback - Visual                 0.4354 0.0597 60  7.292  <.0001 

Results are averaged over the levels of: Part 

lsmAvFeedback <- lsmeans(aAv,~ InterpenetrationFeedback)

NOTE: Results may be misleading due to involvement in interactions

contrast(lsmAvIntEff, method="pairwise", interaction = TRUE)

 InterpenetrationFeedback_pairwise Part_pairwise   estimate     SE df t.ratio p.value
 Both - Electrotactile             Part.1 - Part.2  -0.0995 0.0702 60 -1.417  0.1617 
 Both - NoFeedback                 Part.1 - Part.2   0.1816 0.0702 60  2.585  0.0122 
 Both - Visual                     Part.1 - Part.2   0.1333 0.0702 60  1.898  0.0625 
 Electrotactile - NoFeedback       Part.1 - Part.2   0.2811 0.0702 60  4.002  0.0002 
 Electrotactile - Visual           Part.1 - Part.2   0.2329 0.0702 60  3.315  0.0016 
 NoFeedback - Visual               Part.1 - Part.2  -0.0483 0.0702 60 -0.687  0.4946 

lsmAvPart <- lsmeans(aAv,~ Part)

NOTE: Results may be misleading due to involvement in interactions

contrast(lsmAvIntEff, method="pairwise", interaction = TRUE)

 InterpenetrationFeedback_pairwise Part_pairwise   estimate     SE df t.ratio p.value
 Both - Electrotactile             Part.1 - Part.2  -0.0995 0.0702 60 -1.417  0.1617 
 Both - NoFeedback                 Part.1 - Part.2   0.1816 0.0702 60  2.585  0.0122 
 Both - Visual                     Part.1 - Part.2   0.1333 0.0702 60  1.898  0.0625 
 Electrotactile - NoFeedback       Part.1 - Part.2   0.2811 0.0702 60  4.002  0.0002 
 Electrotactile - Visual           Part.1 - Part.2   0.2329 0.0702 60  3.315  0.0016 
 NoFeedback - Visual               Part.1 - Part.2  -0.0483 0.0702 60 -0.687  0.4946 

lsmAvIntEff <- lsmeans(aAv,~ InterpenetrationFeedback:Part)
contrast(lsmAvIntEff, method="pairwise", interaction = TRUE)

 InterpenetrationFeedback_pairwise Part_pairwise   estimate     SE df t.ratio p.value
 Both - Electrotactile             Part.1 - Part.2  -0.0928 0.0738 60 -1.258  0.2132 
 Both - NoFeedback                 Part.1 - Part.2   0.1674 0.0738 60  2.268  0.0269 
 Both - Visual                     Part.1 - Part.2   0.1509 0.0738 60  2.046  0.0452 
 Electrotactile - NoFeedback       Part.1 - Part.2   0.2602 0.0738 60  3.527  0.0008 
 Electrotactile - Visual           Part.1 - Part.2   0.2437 0.0738 60  3.304  0.0016 
 NoFeedback - Visual               Part.1 - Part.2  -0.0164 0.0738 60 -0.223  0.8243 

?lsmeans()
# Only Interaction Effects
contrast(emmeans(aMax,~ Part:InterpenetrationFeedback), 
         method="pairwise", interaction=TRUE)

 Part_pairwise   InterpenetrationFeedback_pairwise estimate     SE df t.ratio p.value
 Part.1 - Part.2 Both - Electrotactile              -0.1201 0.0718 60 -1.672  0.0997 
 Part.1 - Part.2 Both - NoFeedback                   0.1262 0.0718 60  1.757  0.0841 
 Part.1 - Part.2 Both - Visual                       0.1153 0.0718 60  1.605  0.1137 
 Part.1 - Part.2 Electrotactile - NoFeedback         0.2463 0.0718 60  3.429  0.0011 
 Part.1 - Part.2 Electrotactile - Visual             0.2354 0.0718 60  3.277  0.0017 
 Part.1 - Part.2 NoFeedback - Visual                -0.0109 0.0718 60 -0.152  0.8800 

contrast(emmeans(aMax,~ InterpenetrationFeedback:Part), 
         method="pairwise", interaction=TRUE)

 InterpenetrationFeedback_pairwise Part_pairwise   estimate     SE df t.ratio p.value
 Both - Electrotactile             Part.1 - Part.2  -0.1201 0.0718 60 -1.672  0.0997 
 Both - NoFeedback                 Part.1 - Part.2   0.1262 0.0718 60  1.757  0.0841 
 Both - Visual                     Part.1 - Part.2   0.1153 0.0718 60  1.605  0.1137 
 Electrotactile - NoFeedback       Part.1 - Part.2   0.2463 0.0718 60  3.429  0.0011 
 Electrotactile - Visual           Part.1 - Part.2   0.2354 0.0718 60  3.277  0.0017 
 NoFeedback - Visual               Part.1 - Part.2  -0.0109 0.0718 60 -0.152  0.8800 

contrast(emmeans(aAv,~ Part:InterpenetrationFeedback), 
         method="pairwise", interaction=TRUE)

 Part_pairwise   InterpenetrationFeedback_pairwise estimate     SE df t.ratio p.value
 Part.1 - Part.2 Both - Electrotactile              -0.0928 0.0738 60 -1.258  0.2132 
 Part.1 - Part.2 Both - NoFeedback                   0.1674 0.0738 60  2.268  0.0269 
 Part.1 - Part.2 Both - Visual                       0.1509 0.0738 60  2.046  0.0452 
 Part.1 - Part.2 Electrotactile - NoFeedback         0.2602 0.0738 60  3.527  0.0008 
 Part.1 - Part.2 Electrotactile - Visual             0.2437 0.0738 60  3.304  0.0016 
 Part.1 - Part.2 NoFeedback - Visual                -0.0164 0.0738 60 -0.223  0.8243 

contrast(emmeans(aAv,~ InterpenetrationFeedback:Part), 
         method="pairwise", interaction=TRUE)

 InterpenetrationFeedback_pairwise Part_pairwise   estimate     SE df t.ratio p.value
 Both - Electrotactile             Part.1 - Part.2  -0.0928 0.0738 60 -1.258  0.2132 
 Both - NoFeedback                 Part.1 - Part.2   0.1674 0.0738 60  2.268  0.0269 
 Both - Visual                     Part.1 - Part.2   0.1509 0.0738 60  2.046  0.0452 
 Electrotactile - NoFeedback       Part.1 - Part.2   0.2602 0.0738 60  3.527  0.0008 
 Electrotactile - Visual           Part.1 - Part.2   0.2437 0.0738 60  3.304  0.0016 
 NoFeedback - Visual               Part.1 - Part.2  -0.0164 0.0738 60 -0.223  0.8243 

Both Part.1 - Electrotactile Part.1 -0.2960 0.0710 95.2 -4.171 0.0017 Both Part.1 - NoFeedback Part.1 -0.4434 0.0710 95.2 -6.248 <.0001 Both Part.1 - Visual Part.1 -0.0531 0.0710 95.2 -0.749 0.9952 Both Part.1 - Both Part.2 0.2904 0.0588 64.1 4.936 0.0002

Part.2 Both - Part.2 Electrotactile -0.1965 0.0710 95.2 -2.769 0.115