# Experimental stimuli, results, and analysis code for:

# Four Experiments on the Perception of Bar Charts
# Justin Talbot, Vidya Setlur, and Anushka Anand
# Tableau Software
# InfoVis 2014

# To regenerate the estimates, confidence intervals,
# and figures from the paper, with R installed, run:
# Rscript analysis.R


rm(list=ls())

library(plyr)
library(ggplot2)
library(grid)
library(MASS)

# Make sure .Random.seed exists
invisible(runif(1))

# Options

trimamount <- 0.25


standard_theme <- function() {
    theme_set(theme_bw(base_size = 8))
    theme_update(
    panel.background= element_blank(),
    panel.grid.major.y=element_blank(),
    panel.grid.minor.y=element_blank(),
    panel.border=element_rect(color="#DDDDDD", fill=NA),
    plot.title=element_text( family="Times", size=11, vjust=1),
    axis.title.x= element_text(vjust=-0.25, family="Times", size=8),
    axis.text.x= element_text( family= "Times", size=7),
    axis.title.y= element_text( family="Times", size=8, angle=90, vjust=0.2),
    axis.text.y= element_text( family= "Times", hjust=1, size=7),
    strip.text.x= element_text( family= "Times", size=8),
    strip.text.y= element_text( family= "Times", size=8),
    legend.key = element_rect(color="#FFFFFF", fill="#FFFFFF"),
    legend.text = element_text( family= "Times", size=9),
    legend.title = element_text( color="#FFFFFF", family= "Times", size=9),
    strip.background= element_rect(color="#DDDDDD", fill="#F3F3F3"),
    plot.margin = unit(c(0.01,0.01,0.1,-0.1), "in")
    )
}

standard_xy_theme <- function() {
    theme_set(theme_bw(base_size = 8))
    theme_update(
    panel.background= element_blank(),
    panel.border=element_rect(color="#DDDDDD", fill=NA),
    plot.title=element_text( family="Times", size=11, vjust=1),
    axis.title.x= element_text(vjust=-0.25, family="Times", size=8),
    axis.text.x= element_text( family= "Times", size=7),
    axis.title.y= element_text( family="Times", size=8, angle=90, vjust=0.2),
    axis.text.y= element_text( family= "Times", hjust=1, size=7),
    strip.text.x= element_text( family= "Times", size=8),
    strip.text.y= element_text( family= "Times", size=8),
    legend.key = element_rect(color="#FFFFFF", fill="#FFFFFF"),
    legend.text = element_text( family= "Times", size=9),
    legend.title = element_text( color="#FFFFFF", family= "Times", size=9),
    strip.background= element_rect(color="#DDDDDD", fill="#F3F3F3"),
    plot.margin = unit(c(0.01,0.01,0.1,0.1), "in")
    )
}


# Some utilities

loadfiles <- function(dir) {
    filenames <- list.files(dir, pattern="*.csv", full.names=TRUE)
    df <- NULL
    for (f in filenames){
        data <- read.csv(f, header=TRUE)
        data$user <- basename(f)
        df<-rbind(df, data)
    }
    
    # compute error metrics
    df$err <- df$userratio-df$ratio
    df$abserr <- abs(df$err)
    
    df$ratio <- as.factor(df$ratio)

    # Add a UserID that is ordered by overall accuracy
    u <- by(df$abserr, df['user'], function(x) mean(x, trim=trimamount))
    df$user <- ordered(df$user, levels=names(u)[order(u)])
    df$userid <- factor(as.numeric(df$user))
    
    df
}

# Mean of trimmed means
# dims2 should be a subset of dims1
mean.of.midmeans <- function(x, dims1, dims2, ..., trim=trimamount) {
    cols <- as.list(substitute(list(...))[-1])
    trimmed <- lapply(cols, function(x) call('mean', x, trim=trim))
    overall <- lapply(names(cols), function(x) call('mean', as.symbol(x)))
    names(overall) <- names(cols)
    
    y <- do.call('ddply', c( list(x, dims1, summarize), trimmed) )
    y <- do.call('ddply', c( list(y, dims2, summarize), overall) )
    
    if(length(dims2)==0)
        y[-1]
    else
        y
}

# Bootstrap over users
# Randomly pick 50 users (with replacement)
# Compute f for each replication
# Puts the actual point estimate in the first entry
# For consistency, reset the random seed each time
bootstrap <- function(df, f, n=1000) {
    rs <- .Random.seed
    set.seed(1357)
    
    u <- unique(df$userid)

    s <- replicate(
            n-1,
            f(merge(df, data.frame(userid=sample(u, replace=TRUE)))),
            simplify=FALSE)
            
    .Random.seed <- rs
    
    do.call('rbind', c(list(f(df)), s))
}

# Get the point estimate and the 95% interval in the structure that ggplot wants
bootstrapped <- function(v) {
    x <- quantile(v, c(0.025, 0.975))
    data.frame(y=v[[1]], ymin=x[[1]], ymax=x[[2]])
}

labeller <- function(variable, value) {
    if(variable == 'height') {
        #ifelse(value==25, 'Short', ifelse(value==50, 'Medium', 'Tall'))
        ifelse(value==25, 'Height: 25', ifelse(value==50, 'Height: 50', 'Height: 75'))
    }
    else if(variable == 'distractors') {
        ifelse(value=='No', 'No distractors', ifelse(value=='Short', 'Short distractors', 'Tall distractors'))
    }
    else {
        as.character(value)
    }
}

pivot <- function(df, factors, by) {
    for(i in seq_along(by)) {
        df <- reshape(df, idvar=as.character(factors[!(factors %in% by[1:i])]), timevar=as.character(by[i]), direction="wide")
    }
    df[complete.cases(df),]
}

unpivot <- function(df, varying) {
    reshape(df, varying=varying, times=varying, direction="long", timevar="m", v.names="v")
}

compute <- function(df, factors, pivoted, at, ...) {
    unpivot(
        bootstrap(
            pivot(df, factors, pivoted),
                    function(x) mean.of.midmeans(x,
                factors[!(factors %in% c(pivoted, .(userid)))],
                at,
                ...
                )),
        names(substitute(list(...)))[-1])
}

printi <- function(label, data) {
    d <- getOption('digits')
    options(digits=3)
    interval <- bootstrapped(data$v)
    cat(label, ':    ', interval[[1]], '  [', interval[[2]], ', ', interval[[3]], ']\n', sep="")
    options(digits=d)
}


# Experiment 1
(function() {
    
    cat("Experiment 1 -----------------------------------\n\n")

    # load the files
    df <- loadfiles("data/Experiment1-3-25-14/")
    df$type <- as.factor(ifelse(df$distance==1, 'Adjacent', 'Distant'))
    df$height <- factor(ifelse(df$height==25, '125 pixels',
                        ifelse(df$height==50, '250 pixels', '375 pixels')), levels=c('375 pixels','250 pixels','125 pixels'))
    # the generation code assigned by factor id number, not text label, so fix it up
    df$distractors <- factor(ifelse(df$distractors=='short', 'Tall',
                                ifelse(df$distractors=='tall', 'No', 'Short')),
                                levels = c('Tall', 'Short', 'No'))
    
    # the factors in this experiment
    factors <- .(userid, ratio, height, distractors, type)


    d1 <- compute(df[df$distractors!='No',], factors, .(type), .(),
        'Total'=abserr.Distant-abserr.Adjacent)

    printi('Separation effect (only cases with distractors)', d1)
    cat("\n")


    d2 <- compute(df, factors, .(type, distractors), .(),
        'Separation' = (abserr.Distant.No-abserr.Adjacent.No),
        'Short Distractors' = (abserr.Distant.Short-abserr.Adjacent.Short) - (abserr.Distant.No-abserr.Adjacent.No),
        'Tall Distractors' = (abserr.Distant.Tall-abserr.Adjacent.Tall) - (abserr.Distant.No-abserr.Adjacent.No))

    printi('Contribution of separation only', d2[d2$m=='Separation',])
    printi('Contribution of short distractors only', d2[d2$m=='Short Distractors',])
    printi('Contribution of tall distractors only', d2[d2$m=='Tall Distractors',])
    cat("\n")
    
    
    # Pairwise comparison of contributions of distractors and separation
    
    d2p <- compute(df, factors, .(type, distractors), .(),
    'Short Distractors - Separation' = ((abserr.Distant.Short-abserr.Adjacent.Short) - (abserr.Distant.No-abserr.Adjacent.No)) - (abserr.Distant.No-abserr.Adjacent.No),
    'Tall Distractors - Separation' = ((abserr.Distant.Tall-abserr.Adjacent.Tall) - (abserr.Distant.No-abserr.Adjacent.No)) - (abserr.Distant.No-abserr.Adjacent.No),
    'Tall Distractors - Short Distractors' = ((abserr.Distant.Tall-abserr.Adjacent.Tall) - (abserr.Distant.No-abserr.Adjacent.No)) - ((abserr.Distant.Short-abserr.Adjacent.Short) - (abserr.Distant.No-abserr.Adjacent.No)))
    
    printi('Contribution of short distractors - separation', d2p[d2p$m=='Short Distractors - Separation',])
    printi('Contribution of tall distractors - separation', d2p[d2p$m=='Tall Distractors - Separation',])
    printi('Contribution of tall distractors - short distractors', d2p[d2p$m=='Tall Distractors - Short Distractors',])
    cat("\n")
    
    
    # Separation effect conditioned on reference bar height
    # Estimate only from distractor conditions as above
    
    d3 <- compute(df[df$distractors!='No',], factors, .(type), .(height),
        'Separation Effect'=abserr.Distant-abserr.Adjacent)
    d3 <- d3[-2]
    names(d3) <- c('m', 'v', 'id')
    
    printi('Separation effect conditioned on reference height=125 pixels', d3[d3$m=='125 pixels',])
    printi('Separation effect conditioned on reference height=250 pixels', d3[d3$m=='250 pixels',])
    printi('Separation effect conditioned on reference height=375 pixels', d3[d3$m=='375 pixels',])
    cat("\n")
    
    
    # Pairwise comparison of reference bar conditions
    d3p <- compute(df, factors, .(type, height), .(),
    '375-250' = (`abserr.Distant.375 pixels` - `abserr.Adjacent.375 pixels`) -
    (`abserr.Distant.250 pixels` - `abserr.Adjacent.250 pixels`),
    '375-125' = (`abserr.Distant.375 pixels` - `abserr.Adjacent.375 pixels`) -
    (`abserr.Distant.125 pixels` - `abserr.Adjacent.125 pixels`),
    '250-125' = (`abserr.Distant.250 pixels` - `abserr.Adjacent.250 pixels`) - (`abserr.Distant.125 pixels` - `abserr.Adjacent.125 pixels`))
    
    printi('375 - 250', d3p[d3p$m=='375-250',])
    printi('375 - 125', d3p[d3p$m=='375-125',])
    printi('250 - 125', d3p[d3p$m=='250-125',])
    cat("\n")
    
    
    # Post hoc pairwise analysis looking at the effect of distractors on both types
    d4p <- compute(df, factors, .(distractors), .(),
        'Short - None'=abserr.Short - abserr.No,
        'Tall - None'=abserr.Tall - abserr.No,
        'Tall - Short'=abserr.Tall - abserr.Short
        )
        
    printi('Short - None', d4p[d4p$m=='Short - None',])
    printi('Tall - None', d4p[d4p$m=='Tall - None',])
    printi('Tall - Short', d4p[d4p$m=='Tall - Short',])
    cat("\n")
    
    
    # Make the figures
    
    dx <- rbind(d1, d2)
    
    
    # Figure 3
    dx$m <- factor(dx$m, levels=c('Separation', 'Short Distractors', 'Tall Distractors', 'Total'))
    
    (p <- ggplot(dx, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Separation Effect (percentage points)")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment1_err_diff1.pdf",p,width=3.5,height=1.5)
	embedFonts('figures/experiment1_err_diff1.pdf')
    
    # Figure 4
    dx$m <- factor(dx$m, levels=c('Total', 'Separation', 'Short Distractors', 'Tall Distractors'))
    
    (p <- ggplot(dx, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Contribution to Separation Effect")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment1_err_diff2.pdf",p,width=3.5,height=1.5)
	embedFonts('figures/experiment1_err_diff2.pdf')
    
    # Figure 5
    (p <- ggplot(d3, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Separation Effect conditioned on Reference Bar Height")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment1_err_diff3.pdf",p,width=3.5,height=1.25)
    embedFonts('figures/experiment1_err_diff3.pdf')
    
    # Figure 6
    (p <- ggplot(d4p, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Distractor Effect (percentage points)")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment1_err_diff4.pdf",p,width=3.5,height=1.25)
    embedFonts('figures/experiment1_err_diff4.pdf')

	# All means (not included in paper)
    dall <- compute(df, factors, .(), .(type, height, distractors),
        'Abserr'=abserr)
    (p <- ggplot(dall, aes(distractors, v))
        + facet_grid(height~type, labeller=labeller)
        + stat_summary(fun.data='bootstrapped')
        + xlab("Distractors")
        + ylab("Absolute error (percentage points)")
        + coord_flip())

    standard_theme()
    ggsave("figures/experiment1_all_means.pdf",p,width=3.5,height=2.5)
    embedFonts('figures/experiment1_all_means.pdf')
    
})()


# Experiment 2
(function() {
    
    cat("Experiment 2 -----------------------------------\n\n")
    
    # load the files
    df <- loadfiles("data/Experiment2-3-26-14/")
    df$type <- as.factor(ifelse(df$distance1==1, 'Aligned', 'Unaligned'))
    df$height <- factor(df$height, levels=c(75,50,25))
    df$distractors <- factor(ifelse(df$distractors=='distractors', 'Yes', 'No'), levels=c('Yes', 'No'))
    
    # the factors in this experiment
    factors <- .(userid, ratio, height, distractors, type, dot)
    
    
    d1 <- compute(df[df$distractors!='No'&df$dot=='middle',], factors, .(type), .(),
    'Total'=abserr.Unaligned-abserr.Aligned)
    
    printi('Unalignment effect (only cases with distractors, dot in middle)', d1)
    cat("\n")
    
    
    d2 <- compute(df, factors, .(type, distractors), .(),
        'Unalignment' = (abserr.Unaligned.No - abserr.Aligned.No))
        
    printi('Contribution of unalignment only', d2)

    
    d3 <- compute(df, factors, .(type, distractors), .(),
    'Distractors' = (abserr.Unaligned.Yes - abserr.Aligned.Yes) - (abserr.Unaligned.No - abserr.Aligned.No))
    
    printi('Contribution of distractors only', d3)
    
    d4 <- compute(df, factors, .(type, dot), .(),
    'Dot In Middle' = (abserr.Unaligned.middle - abserr.Aligned.middle) - (abserr.Unaligned.bottom - abserr.Aligned.bottom))
    
    printi('Contribution of dot in middle only', d4)
    cat("\n")
    
    
    d5 <- compute(df[df$type=='Aligned',], factors, .(distractors), .(),
        'Distractor Effect on Aligned Comparisons'=abserr.Yes-abserr.No)
    printi('Distractor Effect on Aligned Comparisons', d5)
    cat("\n")
    
    
    d6 <- compute(df, factors, .(dot), .(),
        'Dot Effect' = abserr.middle - abserr.bottom)
    printi('Dot Effect', d6)
    
    # Make the figures
    
    dx <- rbind(d1, d2, d3, d4)
    
    # Figure 8
    dx$m <- factor(dx$m, levels=c('Unalignment', 'Dot In Middle', 'Distractors', 'Total'))
    
    (p <- ggplot(dx, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Unalignment Effect (percentage points)")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment2_err_diff1.pdf",p,width=3.5,height=1.5)
    embedFonts('figures/experiment2_err_diff1.pdf')
    
    # Figure 9
    dx$m <- factor(dx$m, levels=c('Total', 'Unalignment', 'Distractors', 'Dot In Middle'))
    
    (p <- ggplot(dx, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Contribution to Unalignment Effect")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment2_err_diff2.pdf",p,width=3.5,height=1.5)
    embedFonts('figures/experiment2_err_diff2.pdf')
    
    # Figure 10
    (p <- ggplot(d6, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Dot In Middle Effect")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip(ylim = c(-0.75, 0.1)))
    
    standard_theme()
    ggsave("figures/experiment2_err_diff3.pdf",p,width=3.5,height=0.7)
    embedFonts('figures/experiment2_err_diff3.pdf')
    
    # All means (not included in paper)
    dall <- compute(df, factors, .(), .(type, height, distractors, dot),
        'Abserr'=abserr)
    (p <- ggplot(dall, aes(distractors, v, color=dot))
    + facet_grid(height~type, labeller=labeller)
    + stat_summary(fun.data='bootstrapped')
    + xlab("Distractors")
    + ylab("Absolute error (percentage points)")
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment2_all_means.pdf",p,width=3.5,height=2.5)
    embedFonts('figures/experiment2_all_means.pdf')
})()


# Experiment 3
(function() {
    
    cat("Experiment 3 -----------------------------------\n\n")
        
    df <- loadfiles("data/Experiment3-3-26-14/")
    df$type <- factor(ifelse(df$distance1==1, 'Separated (2 bars)',
                         ifelse(df$distance1==2, 'Separated (1 bar)', 'Adjacent')),
                         levels=c('Adjacent', 'Separated (1 bar)', 'Separated (2 bars)'))
    df$height <- factor(df$height, levels=c(50,25,12.5))
    df$ratio <- as.factor(df$ratio)
    df$distractors <- factor(ifelse(df$distractors=='distractors', 'Yes', 'No'), levels=c('Yes', 'No'))
    
    # the factors in this experiment
    factors <- .(userid, ratio, height, distractors, type)
    

    d1 <- compute(df[df$distractors!='No'&df$type=='Adjacent',], factors, .(), .(),
    'Total'=abserr)
    
    printi('Overall absolute error for divided bar comparisons', d1)
    cat("\n")


    d2 <- compute(df, factors, .(distractors), .(), "Distractors"=abserr.Yes - abserr.No)
    printi('Contribution of distractors to divided bar comparisons', d2)
    
    d3 <- compute(df, factors, .(type), .(), "Separated (1 bar)"=`abserr.Separated (1 bar)` - abserr.Adjacent)
    printi('Contribution of separation (1 bar) to divided bar comparisons', d3)
    
    d4 <- compute(df, factors, .(type), .(), "Separated (2 bars)"=`abserr.Separated (2 bars)` - abserr.Adjacent)
    printi('Contribution of separation (2 bars) to divided bar comparisons', d4)


    # Check for an effect of reference bar height
    d5 <- compute(df, factors, .(), .(height), abserr=abserr)
    d5 <- d5[-2]
    names(d5) <- c('m', 'v', 'id')
    
    printi('Overall effect when height=50', d5[d5$m==50,])
    printi('Overall effect when height=25', d5[d5$m==25,])
    printi('Overall effect when height=12.5', d5[d5$m==12.5,])

    dx <- rbind(d2, d3, d4)
    
    # Figure 12
    (p <- ggplot(dx, aes(x=0.5, y=v))
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + facet_grid(m~., labeller=labeller)
    + xlab(NULL)
    + ylab("Marginal Effects in Divided Bar Comparison")
    + scale_x_continuous(labels=c(), breaks=c(), limits=c(0,1))
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment3_err_diff2.pdf",p,width=3.5,height=1.25)
    embedFonts('figures/experiment3_err_diff2.pdf')
    
    # Figure 13
    x <- bootstrap(df,
        function(x) mean.of.midmeans(x, .(ratio, height, distractors, type), .(type, ratio), err=err))
    
    standard_xy_theme()
    (p <- ggplot(x, aes(ratio, err))
    + facet_grid(~type, labeller=labeller)
    + geom_hline(yintercept=0, color="#AAAAAA")
    + stat_summary(fun.data='bootstrapped')
    + ylab("Bias (percentage points)")
    + xlab("True percent")
    + theme(axis.text.x = element_text(angle=90, vjust=0.5)))
    ggsave("figures/experiment3_bias.pdf",p,width=3.5,height=1.5)
    embedFonts('figures/experiment3_bias.pdf')
    
    # All means (not included in paper)
    dall <- compute(df, factors, .(), .(type, height, distractors),
    'Abserr'=abserr)
    (p <- ggplot(dall, aes(distractors, v))
    + facet_grid(height~type, labeller=labeller)
    + stat_summary(fun.data='bootstrapped')
    + xlab("Distractors")
    + ylab("Absolute error (percentage points)")
    + coord_flip())
    
    standard_theme()
    ggsave("figures/experiment3_all_means.pdf",p,width=3.5,height=2.5)
    embedFonts('figures/experiment3_all_means.pdf')
})()

# Experiment 4
(function() {
    
    cat("Experiment 4 -----------------------------------\n\n")
        
    df <- loadfiles("data/Experiment4-3-26-14/")
    
    # the factors in this experiment
    factors <- .(userid, ratio)
    
    d <- bootstrap(df,
            function(x) mean.of.midmeans(x, .(ratio), .(ratio),
                err=err, abserr=abserr, timetaken=timetaken))
                
    d$ratio <- as.numeric(d$ratio)

    zacks <- data.frame(
        ratio = c(1,8,15,22,29,36,43,50,57,64,71,78,85,92,99),
        abserr = c(1.9,2.6,3.95,4,3.9,4.4,4.25,2.25,4.25,4.15,3.9,3.5,3.7,3,1)
    )

    inverr <- function(x) (2^x)-0.125

    cleveland <- data.frame(
        ratio = c(17.8,26.1,38.3,46.4,56.2,68.2,82.5),
        err = c(0.5,-0.75,-5,-2,0.25,1,0.25),
        abserr = inverr(c(1.4,1.42,1.45,1.5,1.65,1.8,1.9))
        # this is the Type 2 lowess fit in Figure 13 of the Cleveland & McGill paper.
        )
        
    heer <- data.frame(
        ratio = c(18,22,26,32,39,47,55,56,69,85),
        abserr = inverr(c(0.98,1.01,1.05,1.15,1.35,1.6,1.78,1.8,2.1,2.4))
        # this is the Type 2 lowess fit in Figure 3 of the Heer and Bostock paper.
        )

    # Figure 15 (top)
    standard_xy_theme()
    (p <- ggplot(d, aes(ratio, abserr))
    + stat_summary(fun.data='bootstrapped', geom="smooth", fill='#C5C5C5')
    + geom_line(data=zacks, color="#A16058")
    + geom_line(data=heer, color="#17BECF")
    + geom_line(data=cleveland, color="#FF7F0E")
    + xlab("True percent")
    + ylab("Absolute error (percentage points)")
    + scale_x_continuous(breaks=seq(0,100,10), limits=c(0,100)))
    ggsave("figures/experiment4_abserr.pdf",p,width=3.5,height=2)
    embedFonts('figures/experiment4_abserr.pdf')
    
    # Figure 15 (middle)
    standard_xy_theme()
    (p <- ggplot(d, aes(ratio, timetaken))
    + stat_summary(fun.data='bootstrapped', geom="smooth", fill='#C5C5C5')
    + xlab("True percent")
    + ylab("Response time (s)")
    + scale_x_continuous(breaks=seq(0,100,10), limits=c(0,100)))
    ggsave("figures/experiment4_timetaken.pdf",p,width=3.5,height=2)
    embedFonts('figures/experiment4_timetaken.pdf')
    
    # Figure 15 (bottom)
    standard_xy_theme()
    (p <- ggplot(d, aes(ratio, err))
    + stat_summary(fun.data='bootstrapped', geom="smooth", fill='#C5C5C5')
    + xlab("True percent")
    + ylab("Error (percentage points)")
    + scale_x_continuous(breaks=seq(0,100,10), limits=c(0,100)))
    ggsave("figures/experiment4_err.pdf",p,width=3.5,height=2)
    embedFonts('figures/experiment4_err.pdf')
    
    # Figure 16
    df$userratiocnt <- 1:nrow(df)
    tt <- by(df, df[c('userid','userratio')], function(x) x$userratiocnt[order(x$err)])
    for(i in tt) df$userratiocnt[i] <- 1:length(i)
    x <- df[as.numeric(df$userid)<=10 | as.numeric(df$userid)>=41,]
    
    standard_theme()
    (p <- ggplot(x, aes(x=userratio,y=userratiocnt, color=sign(err)*(abs(err)>2.5)))
    + facet_grid(userid~.)
    + geom_point(size=0.75, shape=15)
    + scale_colour_gradient2(low="blue", mid="#B4B4B4", high="#FF7F0E", guide=FALSE)
    + scale_y_continuous(name="", breaks=NULL)
    + scale_x_continuous(name="Participant response", breaks=seq(0,100,10))
    + theme(panel.margin=unit(0.01,'in'))
    )
    ggsave("figures/experiment4_peruser.pdf",p,width=3.5,height=5.8)
    embedFonts('figures/experiment4_peruser.pdf')
    
    # Figure 17
    standard_xy_theme()
    x <- ddply(df, .(userid), summarize,
    err=mean(err, trim=trimamount), abserr=mean(abserr, trim=trimamount), timetaken=mean(timetaken, trim=trimamount))
    (p <- ggplot(x, aes(timetaken, abserr))
    + geom_point(size=1.25)
    + stat_smooth(method="rlm", fill='#C5C5C5')
    + xlab("Response time (s)")
    + ylab("Absolute error (percentage points)")
    )
    ggsave("figures/experiment4_timeerror.pdf",p,width=3,height=2)
    embedFonts('figures/experiment4_timeerror.pdf')
    
    # Figure 18
    standard_xy_theme()
    f <- sort(table(df$userratio))[77:101]
    f <- data.frame(ratio=factor(names(f), levels=names(f)), count=as.vector(f))
    (p <- ggplot(f, aes(count/50, ratio))
    + geom_point(size=1.25)
    + scale_x_continuous(breaks=0:5, labels=c('0', '1x', '2x', '3x', '4x', '5x'), limits=c(0,5.2))
    + xlab("Relative frequency")
    + ylab("Participant responses"))
    ggsave("figures/experiment4_responses.pdf",p,width=3.5,height=3.25)
    embedFonts('figures/experiment4_responses.pdf')
})()
