Loading and preparing data for analysis

data = read.table("XP4ComparativeStudy_logs.csv", header=TRUE, sep=",")

data = data %>% arrange(participantId, tech, block, targetSize, density)

# filtering first trial of each set, because we don't know where in the scene, the cursor come from.
data = data %>% filter(trialId > 1)

Outliers stats

Filtering data that are above mean + 3 * std_dev for each technique:

outliersStats = data %>%
    filter(success == "True") %>% # ignoring error trials for the value of outliers threshold
    group_by(tech) %>% # the outlier threshold is computed for each technique
    summarize(outlierThreshold = mean(finalTimer) + 3 * sd(finalTimer))

# contains only outliers records
dataOutliers = data %>% left_join(outliersStats, by = c("tech")) %>%
    filter(success == "True" & finalTimer > outlierThreshold)

# number of outliers trials for each technique
kable(dataOutliers %>% count(tech))
tech n
1_RayCasting 16
2_RayCursor 13
3_SemiAuto 24
4_RoEtAl 22
# filter data to remove outliers (don't filter if success == false)
data = data %>% left_join(outliersStats, by = c("tech")) %>%
    filter(success == "False" | finalTimer <= outlierThreshold)
# aggregating by set of trials, and computing error rate:
dataErrorMeans = data %>%
    group_by(participantId, tech, block, targetSize, density) %>%
    summarize(errorRate = computeErrorRate(success))

# same aggregation as above, but we compute the mean time:
dataSuccessAggr = data %>%
    filter(success == "True") %>%
    group_by(participantId, tech, block, targetSize, density) %>%
    summarize(meanFinalTimer = mean(finalTimer))

# aggregated data, with mean time and error rate in the same data table:
dataAggr = left_join(dataSuccessAggr, dataErrorMeans, by = c("participantId", "tech", "block", "targetSize", "density"))

# useful if we exclude block 1 in some analysis:
dataAggrBlock23 = dataAggr %>% filter(block > 1)

techniques = unique(data$tech)
blocks = unique(data$block)
targetSizes = unique(data$targetSize)
densities = unique(data$density)

1 Selection time

1.1 Effect of blocks on selection time

1.1.1 Check the normality of data

data.long = melt(dataAggr, id = c("meanFinalTimer", "participantId", "tech", "block", "targetSize", "density"))

data.long$tech = factor(data.long$tech)
data.long$block = factor(data.long$block)
data.long$targetSize = factor(data.long$targetSize)
data.long$density = factor(data.long$density)

m <- aov(meanFinalTimer ~ tech*block, data=data.long)
pander(normalCheck(m))

Shapiro-Wilk normality test

data: res W = 0.90592, p-value < 2.2e-16

Shapiro-Wilk normality test: res
Test statistic P value
0.9059 2.238e-18 * * *

1.1.2 Determine boxcox transform lambda to make residuals normal

boxcox(meanFinalTimer ~ tech*block*targetSize*density, data=data.long, plotit=T)

A lambda of around -0.3 is effective.

1.1.3 Transform data

datatr = data.long %>%
    mutate(meanFinalTimer = meanFinalTimer^(-0.3))
m <- aov(meanFinalTimer ~ tech*block, data=datatr)
pander(normalCheck(m))

Shapiro-Wilk normality test

data: res W = 0.99682, p-value = 0.3192

Shapiro-Wilk normality test: res
Test statistic P value
0.9968 0.3192

1.1.4 Repeated measures ANOVA on transformed data

anova = ezANOVA(datatr, dv=.(meanFinalTimer), wid=.(participantId), within=.(tech,block), detailed=TRUE)

kable(anova_apa(anova, sph_corr ="gg", es = "ges", print=FALSE))
effect text
(Intercept) F(1, 11) = 3292.29, p < .001, getasq > .99
tech F(1.45, 15.93) = 120.06, p < .001, getasq = .71
block F(2, 22) = 21.26, p < .001, getasq = .05
tech:block F(6, 66) = 0.86, p = .529, getasq < .01

1.1.5 Repeated measures ANOVA on non-transformed data

anova = ezANOVA(data.long, dv=.(meanFinalTimer), wid=.(participantId), within=.(tech,block), detailed=TRUE)

kable(anova_apa(anova, sph_corr ="gg", es = "ges", print=FALSE))
effect text
(Intercept) F(1, 11) = 213.93, p < .001, getasq = .92
tech F(1.61, 17.72) = 66.73, p < .001, getasq = .66
block F(2, 22) = 11.31, p < .001, getasq = .03
tech:block F(2.17, 23.91) = 1.57, p = .229, getasq = .01

1.1.6 Selection time per block

kable(data.long %>% group_by(block) %>% summarise(mean = mean(meanFinalTimer)))
block mean
1 2.965367
2 2.723390
3 2.592280

1.1.7 Post-hoc analysis with Bonferroni correction

attach(datatr)
pw <- pairwise.t.test(meanFinalTimer, interaction(block), p.adj = "bonferroni")
detach(datatr)
kable(pw$p.value)
1 2
2 0.1605040 NA
3 0.0168319 1

1.2 Analysis with all blocks

1.2.1 Repeated measures ANOVA on transformed data

anova = ezANOVA(datatr, dv=.(meanFinalTimer), wid=.(participantId), within=.(tech,targetSize,density), detailed=TRUE)

kable(anova_apa(anova, sph_corr ="gg", es = "ges", print=FALSE))
effect text
(Intercept) F(1, 11) = 3292.29, p < .001, getasq > .99
tech F(1.45, 15.93) = 120.06, p < .001, getasq = .70
targetSize F(1, 11) = 11.04, p = .007, getasq = .02
density F(1, 11) = 5.81, p = .035, getasq < .01
tech:targetSize F(3, 33) = 11.80, p < .001, getasq = .03
tech:density F(3, 33) = 3.39, p = .030, getasq = .01
targetSize:density F(1, 11) = 18.78, p = .001, getasq < .01
tech:targetSize:density F(3, 33) = 0.97, p = .417, getasq < .01

1.2.2 Repeated measures ANOVA Selection time on non-transformed data

dataSuccess.long = melt(dataAggr, id = c("meanFinalTimer", "participantId", "tech", "block", "targetSize", "density"))

dataSuccess.long$tech = factor(dataSuccess.long$tech)
dataSuccess.long$block = factor(dataSuccess.long$block)
dataSuccess.long$targetSize = factor(dataSuccess.long$targetSize)
dataSuccess.long$density = factor(dataSuccess.long$density)

anova = ezANOVA(dataSuccess.long, dv=.(meanFinalTimer), wid=.(participantId), within=.(tech,targetSize,density), detailed=TRUE)

kable(anova_apa(anova, sph_corr ="gg", es = "ges", print=FALSE))
effect text
(Intercept) F(1, 11) = 213.93, p < .001, getasq = .92
tech F(1.61, 17.72) = 66.73, p < .001, getasq = .66
targetSize F(1, 11) = 22.80, p < .001, getasq = .02
density F(1, 11) = 0.14, p = .717, getasq < .01
tech:targetSize F(3, 33) = 27.46, p < .001, getasq = .04
tech:density F(1.32, 14.50) = 3.09, p = .092, getasq = .01
targetSize:density F(1, 11) = 6.26, p = .029, getasq < .01
tech:targetSize:density F(1.39, 15.32) = 0.27, p = .689, getasq < .01

1.2.3 Selection time per technique

tech meanFinalTimer
1_RayCasting 1.770374
2_RayCursor 2.721099
3_SemiAuto 1.876999
4_RoEtAl 4.672911

1.2.3.1 Post-hoc analysis with Bonferroni correction

attach(datatr)
pw <- pairwise.t.test(meanFinalTimer, interaction(tech), p.adj = "bonferroni")
detach(datatr)
kable(pw$p.value)
1_RayCasting 2_RayCursor 3_SemiAuto
2_RayCursor 0 NA NA
3_SemiAuto 1 0 NA
4_RoEtAl 0 0 0

1.2.4 Selection time per technique * size

tech targetSize meanFinalTimer
1_RayCasting Large 1.663111
2_RayCursor Large 2.707410
3_SemiAuto Large 1.911865
4_RoEtAl Large 4.270223
1_RayCasting Small 1.877638
2_RayCursor Small 2.734788
3_SemiAuto Small 1.842134
4_RoEtAl Small 5.075599

1.2.4.1 Post-hoc analysis with Bonferroni correction

attach(datatr)
pw <- pairwise.t.test(meanFinalTimer, interaction(tech,targetSize), p.adj = "bonferroni")
detach(datatr)
kable(pw$p.value)
1_RayCasting.Large 2_RayCursor.Large 3_SemiAuto.Large 4_RoEtAl.Large 1_RayCasting.Small 2_RayCursor.Small 3_SemiAuto.Small
2_RayCursor.Large 0.0000000 NA NA NA NA NA NA
3_SemiAuto.Large 0.3864316 0 NA NA NA NA NA
4_RoEtAl.Large 0.0000000 0 0 NA NA NA NA
1_RayCasting.Small 0.0798407 0 1 0.0000000 NA NA NA
2_RayCursor.Small 0.0000000 1 0 0.0000000 0 NA NA
3_SemiAuto.Small 1.0000000 0 1 0.0000000 1 0 NA
4_RoEtAl.Small 0.0000000 0 0 0.0579343 0 0 0

1.2.5 Selection time per technique * density

tech density meanFinalTimer
1_RayCasting High 1.828241
2_RayCursor High 2.824155
3_SemiAuto High 1.894124
4_RoEtAl High 4.535939
1_RayCasting Low 1.712508
2_RayCursor Low 2.618043
3_SemiAuto Low 1.859875
4_RoEtAl Low 4.809883

1.2.5.1 Post-hoc analysis with Bonferroni correction

attach(datatr)
pw <- pairwise.t.test(meanFinalTimer, interaction(tech,density), p.adj = "bonferroni")
detach(datatr)
kable(pw$p.value)
1_RayCasting.High 2_RayCursor.High 3_SemiAuto.High 4_RoEtAl.High 1_RayCasting.Low 2_RayCursor.Low 3_SemiAuto.Low
2_RayCursor.High 0 NA NA NA NA NA NA
3_SemiAuto.High 1 0 NA NA NA NA NA
4_RoEtAl.High 0 0 0 NA NA NA NA
1_RayCasting.Low 1 0 1 0 NA NA NA
2_RayCursor.Low 0 1 0 0 0 NA NA
3_SemiAuto.Low 1 0 1 0 1 0 NA
4_RoEtAl.Low 0 0 0 1 0 0 0

2 Error rate

2.0.1 Check the normality of data

data.long = melt(dataAggr, id = c("errorRate", "participantId", "tech", "block", "targetSize", "density"))

data.long$tech = factor(data.long$tech)
data.long$block = factor(data.long$block)
data.long$targetSize = factor(data.long$targetSize)
data.long$density = factor(data.long$density)

m <- aov(errorRate ~ tech*block, data=data.long)
pander(normalCheck(m))

Shapiro-Wilk normality test

data: res W = 0.84744, p-value < 2.2e-16

Shapiro-Wilk normality test: res
Test statistic P value
0.8474 3.963e-23 * * *

Error rate does not follow a normal distribution.

2.1 Effect of blocks on error rate

2.1.1 Running ANOVA on ART

m = art(errorRate ~ block + (1|participantId), data=data.long)
kable(anova(m)) 
Term F Df Df.res Pr(>F)
block block 2.886354 2 562 0.0566066

2.1.2 Repeated measures ANOVA on non-transformed data

anova = ezANOVA(data.long, dv=.(errorRate), wid=.(participantId), within=.(tech,block), detailed=TRUE)

kable(anova_apa(anova, sph_corr ="gg", es = "ges", print=FALSE))
effect text
(Intercept) F(1, 11) = 47.37, p < .001, getasq = .54
tech F(3, 33) = 14.04, p < .001, getasq = .32
block F(2, 22) = 1.60, p = .224, getasq = .02
tech:block F(6, 66) = 2.31, p = .044, getasq = .05

No learning effect.

2.2 Analysis with all blocks

2.2.1 Running ANOVA on ART

m = art(errorRate ~ tech*targetSize*density + (1|participantId), data=data.long)
kable(anova(m)) 
Term F Df Df.res Pr(>F)
tech tech 25.5157414 3 549 0.0000000
targetSize targetSize 11.9173815 1 549 0.0005990
density density 0.0132082 1 549 0.9085451
tech:targetSize tech:targetSize 2.6162936 3 549 0.0503224
tech:density tech:density 2.7323168 3 549 0.0431323
targetSize:density targetSize:density 3.6180970 1 549 0.0576773
tech:targetSize:density tech:targetSize:density 0.3485347 3 549 0.7902296

2.2.2 Repeated measures ANOVA on non-transformed data

data.long = melt(dataAggr, id = c("errorRate", "participantId", "tech", "block", "targetSize", "density"))

data.long$tech = factor(data.long$tech)
data.long$block = factor(data.long$block)
data.long$targetSize = factor(data.long$targetSize)
data.long$density = factor(data.long$density)

anova = ezANOVA(data.long, dv=.(errorRate), wid=.(participantId), within=.(tech,targetSize,density), detailed=TRUE)

kable(anova_apa(anova, sph_corr ="gg", es = "ges", print=FALSE))
effect text
(Intercept) F(1, 11) = 47.37, p < .001, getasq = .48
tech F(3, 33) = 14.04, p < .001, getasq = .27
targetSize F(1, 11) = 1.93, p = .192, getasq = .01
density F(1, 11) = 1.28, p = .283, getasq < .01
tech:targetSize F(3, 33) = 1.49, p = .236, getasq = .01
tech:density F(3, 33) = 1.34, p = .279, getasq = .02
targetSize:density F(1, 11) = 0.01, p = .906, getasq < .01
tech:targetSize:density F(1.65, 18.16) = 0.20, p = .776, getasq < .01

2.2.3 Error rate per technique

tech errorRate
1_RayCasting 5.737434
2_RayCursor 3.405258
3_SemiAuto 2.434689
4_RoEtAl 11.887401

2.2.3.1 Post-hoc analysis with Bonferroni correction

pw = lsmeans(artlm(m, "tech"), pairwise ~ tech)
## Loading required namespace: lmerTest
## NOTE: Results may be misleading due to involvement in interactions
kable(summary(pw$contrasts))
## Warning in ptukey(sqrt(2) * abst, fam.size, zapsmall(df), lower.tail =
## FALSE): NaNs produced
contrast estimate SE df t.ratio p.value
1_RayCasting - 2_RayCursor 60.67361 16.75024 34.3125000 3.622252 0.0049363
1_RayCasting - 3_SemiAuto 37.14583 16.75024 34.3125000 2.217629 0.1387060
1_RayCasting - 4_RoEtAl -76.12500 16.75024 0.2835744 -4.544710 NaN
2_RayCursor - 3_SemiAuto -23.52778 16.75024 34.3125000 -1.404623 0.5051509
2_RayCursor - 4_RoEtAl -136.79861 16.75024 0.2835744 -8.166962 NaN
3_SemiAuto - 4_RoEtAl -113.27083 16.75024 0.2835744 -6.762339 NaN
attach(data.long)
pw <- pairwise.t.test(errorRate, interaction(tech), p.adj = "bonferroni")
detach(data.long)
kable(pw$p.value)
1_RayCasting 2_RayCursor 3_SemiAuto
2_RayCursor 0.1804213 NA NA
3_SemiAuto 0.0130348 1 NA
4_RoEtAl 0.0000001 0 0

2.2.4 Error rate per technique * size

tech targetSize errorRate
1_RayCasting Large 5.464616
2_RayCursor Large 3.462302
3_SemiAuto Large 2.041997
4_RoEtAl Large 9.995040
1_RayCasting Small 6.010251
2_RayCursor Small 3.348214
3_SemiAuto Small 2.827381
4_RoEtAl Small 13.779762

2.2.5 Error rate per size

kable(data.long %>% group_by(targetSize) %>% summarise(mean = mean(errorRate)))
targetSize mean
Large 5.240989
Small 6.491402

2.2.6 Error rate per technique * density

tech density errorRate
1_RayCasting High 7.076720
2_RayCursor High 4.379960
3_SemiAuto High 3.224206
4_RoEtAl High 11.210318
1_RayCasting Low 4.398148
2_RayCursor Low 2.430556
3_SemiAuto Low 1.645172
4_RoEtAl Low 12.564484

3 Preferences

participant 1_RayCasting 2_RayCursor 3_SemiAuto 4_RoEtAl
0 3 2 1 4
1 2 3 1 4
2 1 2 3 4
3 2 3 1 4
4 3 2 1 4
5 3 2 1 4
6 3 1 2 4
7 2 3 1 4
8 2 3 1 4
9 3 2 1 4
10 1 2 3 4
11 2 3 1 4

Median

technique score
1_RayCasting 2
2_RayCursor 2
3_SemiAuto 1
4_RoEtAl 4

3.0.1 Friedman

## 
##  Friedman rank sum test
## 
## data:  data.matrix(data_tr)
## Friedman chi-squared = 25.3, df = 3, p-value = 1.336e-05

3.0.2 Wilcoxon post-hoc analysis with Bonferroni correction

1_RayCasting 2_RayCursor 3_SemiAuto
2_RayCursor 1.0000000 NA NA
3_SemiAuto 0.0927879 0.0405306 NA
4_RoEtAl 0.0000499 0.0000463 3.22e-05