This webpage was generated with R Markdown.

It is meant to give examples of data analysis and how to interpret it based on real study-data. This webpage is also very much linked to our Alt IHM 2017 paper: La Difference Significative entre Valeurs p et Intervalles de Confiance (available in French here https://hal.inria.fr/hal-01562281/document). An english translation of most of the article is available in the Appendix A of Lonni’s Besancon’s thesis. It is available here: https://tel.archives-ouvertes.fr/tel-01684210/document

If you want to credit this page and its explanations, you can use the following bibtex files:

Packages, environment and helper functions.

Packages and environment

R and RStudio have been used to generate this webpage and its data analysis. You will need to install at least the packages ggplot2 and PropCIs. They are available in the CRAN repository. You can simply type in the console:

install.packages("ggplot2")
install.packages("PropCIs")

Helper Functions

We first need to declare all the helper function that we are going to use in our example.

Confidence Intervals

Here are the helper function for the computation of Bootstrap Confidence Intervals. You can copy/paste into a file for BootstrapCI_macros.

################################################################
# Bootstrap CI Helper Functions.
# Computes bootstrap confidence intervals with the sensible defaults.
# 2014 Pierre Dragicevic
################################################################

# For more details
# See http://www.mayin.org/ajayshah/KB/R/documents/boot.html for boot
# and (Kirby and Gerlanc, 2013) http://web.williams.edu/Psychology/Faculty/Kirby/bootes-kirby-gerlanc-in-press.pdf for bootES and for referencing the bootstrap method in your paper

library(boot)
library(PropCIs)

# If set to TRUE, bootstrap intervals remain the same across executions.
deterministic <- TRUE

# The number of bootstrap resamples. Typically recommended is >= 1,000.
# The higher the number, the more precise the interval but also the slower
# the computation.
replicates <- 10000

# The method for computing intervals. The adjusted bootstrap
# percentile (BCa) method is recommended by (Kirby and Gerlanc, 2013)
# and should work in most cases. For other methods type help(boot.ci).

# FIXME: changing this does not work
intervalMethod <- "bca"

# Number of significant digits used for text output. Using many digits
# is not recommended as it gives a misleading impression of precision.
significantDigits <- 3

#####################################################################

# Statistics

samplemedian <- function(x, d) {
  return(median(x[d]))
}

samplemean <- function(x, d) {
  return(mean(x[d]))
}

# Data transformations

# No transformation, yields arithmetic means.
identityTransform <- c(
  function(x) {return (x)}, # the transformation
  function(x) {return (x)}, # the inverse transformation
  TRUE                      # TRUE if increasing, FALSE otherwise
)

# Log transformation, yields geometric means.
logTransform <- c(
  function(x) {return (log(x))}, # the transformation
  function(x) {return (exp(x))}, # the inverse transformation
  TRUE                         # TRUE if increasing, FALSE otherwise
)

# Inverse transformation, yields harmonic means.
inverseTransform <- c(
  function(x) {return (1/x)}, # the transformation
  function(x) {return (1/x)}, # the inverse transformation
  FALSE                     # TRUE if increasing, FALSE otherwise
)

# Returns the point estimate and confidence interval in an array of length 3
bootstrapCI <- function(statistic, datapoints) {
  # Compute the point estimate
  pointEstimate <- statistic(datapoints)
  # Make the rest of the code deterministic
  if (deterministic) set.seed(0)
  # Generate bootstrap replicates
  b <- boot(datapoints, statistic, R = replicates, parallel="multicore")
  # Compute interval
  ci <- boot.ci(b, type = intervalMethod)
  # Return the point estimate and CI bounds
  # You can print the ci object for more info and debug
  lowerBound <- ci$bca[4]
  upperBound <- ci$bca[5]
  return(c(pointEstimate, lowerBound, upperBound))
}

# Returns the mean and its confidence interval in an array of length 3
bootstrapMeanCI <- function(datapoints) {
  return(bootstrapCI(samplemean, datapoints))
}

# Returns the median and its confidence interval in an array of length 3
bootstrapMedianCI <- function(datapoints) {
  return(bootstrapCI(samplemedian, datapoints))
}

# Returns the mean and its "classical" confidence interval in an array of length 3.
# Uses the t-distribution confidence interval for samples that are approximately normally
# distributed and/or large. Not a bootstrap method but included here for convenience.
exactMeanCI <- function(datapoints, transformation = identityTransform) {
  datapoints <- transformation[[1]](datapoints)
  pointEstimate <- mean(datapoints)
  ttest <- t.test(datapoints)
  lowerBound <- ttest[4]$conf.int[1]
  upperBound <- ttest[4]$conf.int[2]
  if (transformation[[3]])
    return(transformation[[2]](c(pointEstimate, lowerBound, upperBound)))
  else
    return(transformation[[2]](c(pointEstimate, upperBound, lowerBound)))
}

# Returns a percentage and its confidence interval, for binary observations.
# Uses a confidence interval method for proportions. Not a bootstrap method but
# included here for convenience.
percentCI <- function(datapoints, value) {
  ncorrect <- length(datapoints[datapoints == value])
  total <- length(datapoints)
  percent <- 100 * ncorrect / total
  proportionCI <- midPci(x = ncorrect, n = total, conf.level = 0.95)
  percentCIlow <- 100 * proportionCI$conf.int[1]
  percentCIhigh <- 100 * proportionCI$conf.int[2]
  return (c(percent, percentCIlow, percentCIhigh))
}

# Returns the point estimate and confidence interval in a human-legible text format
bootstrapCI.text <- function(statistic, datapoints, unit) {
  result <- bootstrapCI(statistic, datapoints)
  point <- result[1]
  interval <- c(result[2], result[3])
  # Format results in human-legible format using APA style
  text <- paste(
    prettyNum(point, digits=significantDigits),
    unit,
    ", 95% CI [",
    prettyNum(interval[1], digits=significantDigits),
    unit,
    ", ",
    prettyNum(interval[2], digits=significantDigits),
    unit,
    "]",
    sep = "")
  return(text)
}

# Returns the mean and its confidence interval in a human-legible text format
bootstrapMeanCI.text <- function(datapoints, unit) {
  return(bootstrapCI.text(samplemean, datapoints, unit))
}

# Returns the median and its confidence interval in a human-legible text format
bootstrapMedianCI.text <- function(datapoints, unit) {
  return(bootstrapCI.text(samplemedian, datapoints, unit))
}

Plotting function helper

This function is designed to plot point estimates and confidence intervals for several factors. It is designed to be generic but has been thoroughly tested.

# 2015--2017 Wesley Willett, Petra Isenberg and Lonni Besancon

library(ggplot2)
library(reshape2)


barChart <- function(resultTable, techniques, nbTechs = -1, ymin, ymax, xAxisLabel = "I am the X axis", yAxisLabel = "I am the Y Label"){
  #tr <- t(resultTable)
  if(nbTechs <= 0){
    stop('Please give a positive number of Techniques, nbTechs');
  }
  
  tr <- as.data.frame(resultTable)
  nbTechs <- nbTechs - 1 ; # seq will generate nb+1
  
  #now need to calculate one number for the width of the interval
  tr$CI2 <- tr$upperBound_CI - tr$mean_time
  tr$CI1 <- tr$mean_time - tr$lowerBound_CI
  
  #add a technique column
  tr$technique <- factor(seq.int(0, nbTechs, 1));
  
  breaks <- c(as.character(tr$technique));
  print(tr)
  g <- ggplot(tr, aes(x=technique, y=mean_time)) + 
    geom_bar(stat="identity",fill = I("#CCCCCC")) +
    geom_errorbar(aes(ymin=mean_time-CI1, ymax=mean_time+CI2),
                  width=0,                    # Width of the error bars
                  size = 1.1
    ) +
    #labs(title="Overall time per technique") +
    labs(x = xAxisLabel, y = yAxisLabel) + 
    scale_y_continuous(limits = c(ymin,ymax)) +
    scale_x_discrete(name="",breaks,techniques)+
    coord_flip() +
    theme(panel.background = element_rect(fill = 'white', colour = 'white'),axis.title=element_text(size = rel(1.2), colour = "black"),axis.text=element_text(size = rel(1.2), colour = "black"),panel.grid.major = element_line(colour = "#DDDDDD"),panel.grid.major.y = element_blank(), panel.grid.minor.y = element_blank())+
    geom_point(size=4, colour="black")         # dots
  
  print(g)
}

The data

The data gathered are about 4 different techniques from the paper Pressure-Based Gain Factor Control for Mobile 3D Interaction using Locally-Coupled Devices which is available from: “https://hal.inria.fr/hal-01436172”. The idea was to compare different gain factor manipulations for a 3D docking task. The goal of a 3D docking task is to ask users to match an objet’s position and orientation ot a target position and orientation. We can then use the angular and Euclidean distances between the object and the target as a measure of accuracy.

We used 4 different techniques to manipulate the gain factor:

pressure-based manipulations, slider-based manipulations, rate-control manipulation, and speed-based manipulations.

The names of these techniques will be used in the following example code. Modify these names so that they fit your needs.

We will look at three different performance measures for these 4 techniques:

This data does not include time data (which is also very interesting), we will deal with this example later.

Data Analysis

We present here the data analysis for the Euclidean Distance, the Angular Distance, and the workload per technique. For all of these, the distribution of data is far from normal so we use the Bootstrap CIs helper functions presented above.

We first need to read the data with a line like:

################ BOOTSTRAP CI #################

data = read.table("./Info.csv", header=T, sep=",")
data = aggregate(data,by=list(data$pID), FUN=mean, na.rm=TRUE)

The data for the Euclidean Distance and the Angular Distance is summarized below:

##     Group.1      EucliPressure   AngularPressure     EucliRate      
##  Min.   : 1.00   Min.   : 3.13   Min.   :  7.234   Min.   :  7.987  
##  1st Qu.: 6.75   1st Qu.:11.75   1st Qu.: 25.299   1st Qu.: 29.167  
##  Median :12.50   Median :26.36   Median : 40.187   Median : 80.384  
##  Mean   :12.50   Mean   :31.73   Mean   : 45.614   Mean   :114.394  
##  3rd Qu.:18.25   3rd Qu.:49.43   3rd Qu.: 61.346   3rd Qu.:161.468  
##  Max.   :24.00   Max.   :86.47   Max.   :105.215   Max.   :427.770  
##   AngularRate       EucliSpeed      AngularSpeed      EucliSlider     
##  Min.   : 13.88   Min.   : 3.187   Min.   :  9.618   Min.   :  2.076  
##  1st Qu.: 42.92   1st Qu.:19.007   1st Qu.: 26.425   1st Qu.: 10.895  
##  Median : 73.53   Median :43.053   Median : 46.461   Median : 34.649  
##  Mean   : 66.90   Mean   :45.182   Mean   : 52.150   Mean   : 43.719  
##  3rd Qu.: 95.28   3rd Qu.:67.917   3rd Qu.: 74.324   3rd Qu.: 69.486  
##  Max.   :119.37   Max.   :99.690   Max.   :125.185   Max.   :119.997  
##  AngularSlider          pID       
##  Min.   :  5.613   Min.   : 1.00  
##  1st Qu.: 18.792   1st Qu.: 6.75  
##  Median : 35.933   Median :12.50  
##  Mean   : 47.440   Mean   :12.50  
##  3rd Qu.: 71.374   3rd Qu.:18.25  
##  Max.   :119.719   Max.   :24.00

Euclidean Distance

Code

Here is the code to compute the means and CIs for each of the 4 techniques.

data = aggregate(data,by=list(data$pID), FUN=mean, na.rm=TRUE)

resultPressure = bootstrapMeanCI(data$EucliPressure)
resultRate = bootstrapMeanCI(data$EucliRate)
resultSpeed = bootstrapMeanCI(data$EucliSpeed)
resultSlider = bootstrapMeanCI(data$EucliSlider)

analysisData = c()
analysisData$ratio = c("Slider","Speed","Rate","Pressure")

analysisData$pointEstimate = c(resultSlider[1],resultSpeed[1],resultRate[1],resultPressure[1])
analysisData$ci.max = c(resultSlider[3], resultSpeed[3], resultRate[3], resultPressure[3])
analysisData$ci.min = c(resultSlider[2], resultSpeed[2], resultRate[2], resultPressure[2])

datatoprint <- data.frame(factor(analysisData$ratio),analysisData$pointEstimate, analysisData$ci.max, analysisData$ci.min)
colnames(datatoprint) <- c("technique", "mean_time", "lowerBound_CI", "upperBound_CI ") #We use the name mean_time for the value of the mean even though it's not a time, it's just to parse the data for the plot
barChart(datatoprint,analysisData$ratio ,nbTechs = 4, ymin = 0, ymax = 170, "", "Euclidean Distance per technique. Error Bars, Bootstrap 95% CIs")
##   technique mean_time lowerBound_CI upperBound_CI         CI2        CI1
## 1         0  43.71918      60.49831       29.43647 -14.282712 -16.779135
## 2         1  45.18153      57.34625       33.60464 -11.576890 -12.164717
## 3         2 114.39352     166.24510       77.10322 -37.290293 -51.851581
## 4         3  31.72767      41.67720       23.08464  -8.643031  -9.949531

Interpretation

While these results can be interpreted using a binary approach based on the spacing of the confidence intervals, estimation calls for more nuance. Here are how the results were interpreted in the paper:

These results show strong evidence for a better performance of pressure-based, speed-based, and slider-based control compared to rate control. There is also weak evidence for a better performance of pressure-control over the slider-based and speed-based methods.

The second sentence is a bit risky given the rather large overlaps but we will see that it’s confirmed by an analysis of effect sizes.

Going further with a pair-wise difference to show effect sizes

While we managed to avoid binary analysis of the previous results, this is still not enough because we haven’t talked about effect sizes just yet.

We decide to compute the difference between pressure-based control (which seems to be our best technique) and all the other techniques. Since it’s a within-subject design, we compute the difference for every subject individually

data = aggregate(data,by=list(data$pID), FUN=mean, na.rm=TRUE)
data$EucliPR = (data$EucliRate-data$EucliPressure)
data$EucliPSp = (data$EucliSpeed-data$EucliPressure)
data$EucliPSl = (data$EucliSlider-data$EucliPressure)

resultPR = bootstrapMeanCI(data$EucliPR)
resultPSp = bootstrapMeanCI(data$EucliPSp)
resultPSl = bootstrapMeanCI(data$EucliPSl)

analysisData = c()
analysisData$ratio = c("Slider-Pressure","Speed-Pressure","Rate-Pressure")
analysisData$pointEstimate = c(resultPSl[1],resultPSp[1], resultPR[1])
analysisData$ci.max = c(resultPSl[3], resultPSp[3], resultPR[3])
analysisData$ci.min = c(resultPSl[2], resultPSp[2], resultPR[2])

datatoprint <- data.frame(factor(analysisData$ratio),analysisData$pointEstimate, analysisData$ci.max, analysisData$ci.min)
colnames(datatoprint) <- c("technique", "mean_time", "lowerBound_CI", "upperBound_CI ")
barChart(datatoprint,analysisData$ratio ,nbTechs = 3, ymin = 0 , ymax = 150, "", "Pair-wise differences. Error Bars, Bootstrap 95% CIs")
##   technique mean_time lowerBound_CI upperBound_CI         CI2        CI1
## 1         0  11.99151      25.40051       1.963988 -10.027519 -13.409003
## 2         1  13.45386      21.02514       6.187039  -7.266819  -7.571285
## 3         2  82.66585     132.93559      50.898956 -31.766890 -50.269748

We can now interpret these results.

The fact that none of the confidence intervals for these differences overlaps with 0 supports our finding that pressure-based control of the gain factor is more accurate than all other techniques. Analyzing the ratios between the techniques leads to errors about 4 times smaller on average (in terms of Euclidean distances) than rate-control, and about 2 times smaller than speed-based control.

The results are clearer than before (when we compared technique means) because it’s a within-subjects design.

Angular Distance

Code

Similarly we can compute the mean and CIs for each of the techniques for the angular distance with the following code.

resultPressure = bootstrapMeanCI(data$AngularPressure)
resultRate = bootstrapMeanCI(data$AngularRate)
resultSpeed = bootstrapMeanCI(data$AngularSpeed)
resultSlider = bootstrapMeanCI(data$AngularSlider)

analysisData = c()
analysisData$ratio = c("Slider","Speed","Rate","Pressure")

analysisData$pointEstimate = c(resultSlider[1],resultSpeed[1],resultRate[1],resultPressure[1])
analysisData$ci.max = c(resultSlider[3], resultSpeed[3], resultRate[3], resultPressure[3])
analysisData$ci.min = c(resultSlider[2], resultSpeed[2], resultRate[2], resultPressure[2])

datatoprint <- data.frame(factor(analysisData$ratio),analysisData$pointEstimate, analysisData$ci.max, analysisData$ci.min)
colnames(datatoprint) <- c("technique", "mean_time", "lowerBound_CI", "upperBound_CI ")
barChart(datatoprint,analysisData$ratio ,nbTechs = 4, ymin = 0, ymax = 85, "", "Angular Distances Per Technique. Error Bars, Bootstrap 95% CIs")
##   technique mean_time lowerBound_CI upperBound_CI         CI2       CI1
## 1         0  47.44002      61.61774       35.57512 -11.864900 -14.17772
## 2         1  52.14977      66.97297       39.94621 -12.203564 -14.82319
## 3         2  66.90453      79.29061       53.55713 -13.347401 -12.38608
## 4         3  45.61402      56.61330       36.03621  -9.577812 -10.99928

Additionally and following the previous example, we can compute the pair-wise comparisons between the techniques. The code can easily be derived from the previous example and the figure is: