Yvonne Jansen - R Scripts Analysis of experiment data

This is an R HTML document in which we document the analysis scripts for our study on the size of physical bars and spheres. An author's version of the corresponding publication can be accessed here.

Helper Functions

# load all required libraries
pacman::p_load(stats, nls2, nlstools, boot, ggplot2, grid, gridExtra, plyr, robustbase)




##########################################################################################
# Here are our helper functions 

# multmerge, a function that takes a directory and extension as strings and loads all files in the indicated directory with that extension  
multmerge = function(mypath, ext){
  filenames=list.files(path=mypath, pattern=paste("*\\.", ext, "$", sep=""), full.names=TRUE)
  datalist = lapply(filenames, function(x){read.csv(file=x,header=T)})
  Reduce(function(x,y) {merge(x,y, all=T)}, datalist)
}

# adjusting the ggplot theme
mytheme <-     theme_bw() +
  theme( # remove the vertical grid lines
    plot.background = element_blank(),
    legend.position = "none",
    axis.ticks.margin=unit(c(1.5,-1.0),'mm'),
    axis.title.y = element_text(size=12),
    axis.title.x = element_text(size=12, vjust=-0.7),
    axis.text.y = element_text(size=10, colour="#888888"),
    axis.text.x = element_text(size=10, colour="#888888"),
    panel.border = element_rect(fill=NULL, colour="#888888"), # element_blank(),
    panel.grid.minor.x = element_line(size=0.01, color="#000000"),
    panel.grid.major.x = element_line(size=0.03, color="#000000"),
    panel.margin = unit(2, "lines"),
    axis.line = element_line(size=0.1, color="#000000"),
    axis.ticks = element_line(size=0.1, color="#000000"),
    axis.line.x = element_blank(),
    panel.grid.minor.x = element_line(size=0.01, color="#000000"),
    panel.grid.major.x = element_line(size=0.03, color="#000000"),
    strip.background = element_rect(colour="#EEEEEE", fill="#EEEEEE"),
    strip.text.y = element_text(size=12),
    strip.text.x = element_text(size=12)
  )
myaxes <- coord_fixed(ratio = 1, xlim = c(0,100), ylim = c(0,100)) # expand_limits(y=c(0,100), x=c(0,100))


# an inc function for convenience since R doesn't have one
inc <- function(x)
{
 eval.parent(substitute(x <- x + 1))
}

#end of helper functions
#########################################################################################

Loading, Checking, and Preparing Data

# load data

data <- multmerge("logs/" , "csv")

# data sanity check
sanity <- c(logical())

numParticipants  <- length(unique(data$participant))
okParticipants  <- numParticipants == 10
sanity  <- append(sanity, okParticipants)

numConditions <- length(unique(data$condition))
okConditions  <- numConditions == 4
sanity  <- append(sanity, okConditions)

numStimuli  <-  length(unique(data$code))
okStimuli  <-  numStimuli == 56
sanity  <- append(sanity, okStimuli)

numMethods  <-  length(unique(data$method))
okMethods  <- numMethods == 2
sanity  <- append(sanity, okMethods)

numTasks  <- numStimuli * numMethods
numTasksPerParticipant  <- daply(data, .(participant), nrow)
okTasks  <- all(numTasksPerParticipant == numTasks)
sanity  <- append(sanity, okTasks)

cat(ifelse(all(sanity == TRUE), "Successfully loaded data.\n", paste("There was a problem loading data: \n",
  "Number of participants:", numParticipants, " -> ", ifelse(okParticipants, "OK", "WRONG") , "\n",
  "Number of conditions:", numConditions, " -> ", ifelse(okConditions, "OK", "WRONG") , "\n",
  "Number of stimuli:", numConditions, " -> ", ifelse(okStimuli, "OK", "WRONG") , "\n",
  "Number of methods:", numMethods, " -> ", ifelse(okMethods, "OK", "WRONG") , "\n",
  "Number of experimental tasks:", numTasks, " -> ", ifelse(okTasks, "OK", "WRONG") , "\n")))

## Successfully loaded data.

# Compute the estimate for the constant sum method to make the values between the two estimation methods comparable

data  <- ddply(data, .(), transform, estimate = ifelse(data$method == "ratioPercentage", estimateRatio.Percentage, (100 - data$estimateConstant.Sum)/data$estimateConstant.Sum[29]*100))

# create a column indicating whether it's a bar or a sphere stimuli

data <- ddply(data, .(), transform, shape = ifelse(grepl("bar", condition), "bar", "sphere"))
data  <- ddply(data, .(), transform, method = ifelse(grepl("ratio", method), "ratio estimation", "constant sum"))

# compute the true value from the actual physical size of the stimuli.  

data  <- ddply(data, .(), transform, trueP = size1 / size2 * 100)


# compute the accuracy, i.e., the difference between the estimate and true ratio  

data  <- ddply(data, .(), transform, accuracy = abs(estimate - trueP) )

data  <- ddply(data, .(), transform, rawError = estimate - trueP )

Rationale for computing the "true" ratio from 1D measures

The physical sizes of both our physical marks can be described through a single number: length for bars and diameter for spheres. We use these to compute the "true ratio" between our stimuli. By using linear measures to describe the ratios between our stimulis, we are able to analyze the "perception dimensionality". To illustrate our rationale, consider the charts below showing randomly drawn numbers plotting on the x-axis and linear, quadratic, and cubic transformations respectively on the y-axis. Each transformation results in a distinctively different shape of the point cloud. 

testdata <- data.frame(relation=character(), trueP=numeric(), estimate=numeric())

as.numeric(Sys.time())-> t; set.seed((t - floor(t)) * 1e8 -> seed);
tsample <- sample(1:100, 112, replace=T)
lest  <- tsample
testdata <- rbind(testdata, data.frame(relation=rep("linear", times=112), trueP=tsample, estimate=lest))
qest <- (lest/100)^2 * 100
testdata <- rbind(testdata, data.frame(relation=rep("quadratic", times=112), trueP=tsample, estimate=qest))
cest <- (lest/100)^3 * 100
testdata <- rbind(testdata, data.frame(relation=rep("cubic", times=112), trueP=tsample, estimate=cest))

qplot(data=testdata, trueP, estimate, facets=~relation, position=position_jitter(w=0, h=3), size=I(1.5), xlab="true ratio") + mytheme + myaxes
Random data to illustrate different curves for linear, quadratic, and cubic relationships between the true ratio and an estimate.

Analysis

Now the data is ready so that we can take a first look. We are interested in comparing the estimated values against the true values for both shapes and both estimation methods:

qplot(data=data, trueP, estimate, facets = method ~ shape, alpha = I(0.3), xlab="true ratio", ylab="estimated ratio") + mytheme + myaxes
Basic scatterplot matrix similar to Figure 8 in the paper

Here we can see that estimates made using the constant sum method (top row) have a higher spread than the ratio estimates. Furthermore, the estimates for bars lie more on a straight line than the estimates for spheres. While we expected that bars, who grow linearly in only one direction no matter whether it's a 2D or 3D bar, are more accurately perceived than spheres, we did not expect the high variation with the constant sum method.

Our next analysis step is now to take a closer look at the observed variances. What we want to know is:

  1. How accurate are people overall?
  2. Is the variance due to variance between or within people? That is, if a person is asked to make estimates for two similar ratios, how similar are the two estimates? If two people are asked to make an estimate for one set of stimuli, how similar are their estimates?
  3. Are estimates systematically biased?
We now investigate these three questions in turn.

Accuracy

We want to compare how accurately people estimated shapes in our experiment against data from previous work. Our article compares our findings against those by Spence [1990] who only studied two-dimensional shapes with and without extraneous dimensions, but we might include more here in the future. 

# compute accuracies


# we want bootstrapped confidence intervals, so we need a function to compute them conveniently
# (these convenience functions are borrowed from Pierre Dragicevic)
samplemean <- function(x, d) {
  return(mean(x[d]))
}

# Returns the point estimate and confidence interval in an array of length 3
bootstrapCI <- function(datapoints) {
  datapoints <- datapoints[!is.na(datapoints)]
  # Compute the point estimate
        pointEstimate <- samplemean(datapoints)
        # Make the rest of the code deterministic
        # if (deterministic) set.seed(0)
        # Generate bootstrap replicates
        b <- boot(datapoints, samplemean, R = 10000, parallel="multicore")
        # Compute interval
        ci <- boot.ci(b, type = "bca")
        # 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))
}

# a data frame to hold our confidence intervals
accuracies  <- data.frame(pub=character(), shape=character(), measure=character(), estimate=numeric(), ci.low=numeric(), ci.high=numeric(), order=integer())

ord <- 10
cis <- bootstrapCI(subset(data, shape == "bar" & method == "ratio estimation")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="bar (RE)", measure="height", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))

cis <- bootstrapCI(subset(data, shape == "bar" & method == "constant sum")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="bar (CS)", measure="height", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))

cis <- bootstrapCI(subset(data, shape == "sphere" & method == "ratio estimation")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (RE)", measure="diameter", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))

cis <- bootstrapCI(subset(data, shape == "sphere" & method == "constant sum")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (CS)", measure="diameter", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))


ord <- 0
# manually entered data from previous work
prevD  <- data.frame(pub=character(), shape=character(), measure=character(), estimate=numeric(), ci.low=numeric(), ci.high=numeric(), order=integer())

prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="line (H)", measure="length", estimate=3.1519, ci.low=2.84, ci.high=3.4479, order=inc(ord)))
prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="line (V)", measure="length", estimate=3.75, ci.low=3.026, ci.high=4.4479, order=inc(ord)))
prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="bar", measure="height", estimate=2.75, ci.low=2.5, ci.high=3.05, order=inc(ord)))

prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="box", measure="height", estimate=3.2434, ci.low=2.8339, ci.high=3.6527, order=inc(ord)))
prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="cylinder", measure="height", estimate=3.344, ci.low=3.0338, ci.high=3.656, order=inc(ord)))

# merge the two data frames into one for plotting
cidf  <-rbind(prevD, accuracies)

# Cleveland dot plot

# R doesn't have a built-in function to create a dot plot similar to Cleveland & McGill [1984]. We therefore create a fake data set containing all the locations of the dots we want to see. We then plot those first and draw our confidence intervals on top.

# the fake data for the dots
dotData <- data.frame(x = rep(0:20, times=nrow(cidf)), y = rep(paste(cidf$shape, cidf$measure), each=21), order = rep(cidf$order, each = 21))

# we start by plotting the dots
qplot(data=dotData, reorder(y, order), x, alpha=I(0.2)) +
  expand_limits(x=c(1,5), y=c(0,20)) +
  mytheme +
  coord_flip() +
  # and now we plot the dataframe containing our confidence intervals on top
  geom_pointrange(data=cidf, aes(y=estimate, x=reorder(paste(shape, measure), order), ymin=ci.low, ymax=ci.high), width=1, size=0.7) +
  labs(x = "", y = "Average accuracy (absolute discrepancy in percent)") +
   annotate("text", x= 3.2, y = 16, label = cidf$pub[1]) +
   annotate("segment", x=0.5, xend=5.5, y = 12, yend = 12) +
   annotate("segment", x=5.5, xend=5.5, y = 11.5, yend = 12) +
   annotate("segment", x=0.5, xend=0.5, y = 11.5, yend = 12)
Dot plot showing accuracies similar to Figure 11 in the paper
# compute accuracies


# we want bootstrapped confidence intervals, so we need a function to compute them conveniently
# (these convenience functions are borrowed from Pierre Dragicevic)
samplemean <- function(x, d) {
  return(mean(x[d]))
}

# Returns the point estimate and confidence interval in an array of length 3
bootstrapCI <- function(datapoints) {
  datapoints <- datapoints[!is.na(datapoints)]
        # Compute the point estimate
        pointEstimate <- samplemean(datapoints)
        # Make the rest of the code deterministic
        # if (deterministic) set.seed(0)
        # Generate bootstrap replicates
        b <- boot(datapoints, samplemean, R = 10000, parallel="multicore")
        # Compute interval
        ci <- boot.ci(b, type = "bca")
        # 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))
}

# a data frame to hold our confidence intervals
accuracies  <- data.frame(pub=character(), shape=character(), measure=character(), estimate=numeric(), ci.low=numeric(), ci.high=numeric(), order=integer())

ord <- 10
cis <- bootstrapCI(subset(data, shape == "bar" & method == "ratio estimation")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="bar (RE)", measure="height", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))

cis <- bootstrapCI(subset(data, shape == "bar" & method == "constant sum")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="bar (CS)", measure="height", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))

cis <- bootstrapCI(subset(data, shape == "sphere" & method == "ratio estimation")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (RE)", measure="diameter", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))

cis <- bootstrapCI(subset(data, shape == "sphere" & method == "constant sum")$accuracy)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (CS)", measure="diameter", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order=inc(ord)))


ord <- 0
# manually entered data from previous work
prevD  <- data.frame(pub=character(), shape=character(), measure=character(), estimate=numeric(), ci.low=numeric(), ci.high=numeric(), order=integer())

prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="line (H)", measure="length", estimate=3.1519, ci.low=2.84, ci.high=3.4479, order=inc(ord)))
prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="line (V)", measure="length", estimate=3.75, ci.low=3.026, ci.high=4.4479, order=inc(ord)))
prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="bar", measure="height", estimate=2.75, ci.low=2.5, ci.high=3.05, order=inc(ord)))

prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="box", measure="height", estimate=3.2434, ci.low=2.8339, ci.high=3.6527, order=inc(ord)))
prevD  <- rbind(prevD, data.frame(pub="Spence [1990]", shape="cylinder", measure="height", estimate=3.344, ci.low=3.0338, ci.high=3.656, order=inc(ord)))

# merge the two data frames into one for plotting
cidf  <-rbind(prevD, accuracies)

# Cleveland dot plot

# R doesn't have a built-in function to create a dot plot similar to Cleveland & McGill [1984]. We therefore create a fake data set containing all the locations of the dots we want to see. We then plot those first and draw our confidence intervals on top.

# the fake data for the dots
dotData <- data.frame(x = rep(0:20, times=nrow(cidf)), y = rep(paste(cidf$shape, cidf$measure), each=21), order = rep(cidf$order, each = 21))

# we start by plotting the dots
qplot(data=dotData, reorder(y, order), x, alpha=I(0.2)) +
  expand_limits(x=c(1,5), y=c(0,20)) +
  mytheme +
  coord_flip() +
  # and now we plot the dataframe containing our confidence intervals on top
  geom_pointrange(data=cidf, aes(y=estimate, x=reorder(paste(shape, measure), order), ymin=ci.low, ymax=ci.high), width=1, size=0.7) +
  labs(x = "", y = "Average accuracy (absolute discrepancy in percent)") +
   annotate("text", x= 3.2, y = 16, label = cidf$pub[1]) +
   annotate("segment", x=0.5, xend=5.5, y = 12, yend = 12) +
   annotate("segment", x=5.5, xend=5.5, y = 11.5, yend = 12) +
   annotate("segment", x=0.5, xend=0.5, y = 11.5, yend = 12)
Dot plot showing accuracies similar to Figure 11 in the paper

As we illustrated earlier, the "true" ratio between two spheres can be computed in different ways. Since the computation of accuracy depends on the definition of what is considered the "true" ratio, we now extend the above figure to include the estimation accuracy for spheres under the assumption that the spheres' surfaces or the spheres' volumes represent their "true" size. 

# compute the true value from the actual physical size of the stimuli.  

sphereSur <- function(d){ # d is the sphere's diameter
  return(4 * pi * (d/2)^2)
}

sphereVol <-function(d){ # d is the sphere's diameter
    return(4 / 3 * pi * (d/2)^3)
  }

data  <- ddply(data, .(), transform, truePSur = ifelse(shape == "sphere", sphereSur(size1) / sphereSur(size2) * 100, trueP))

data  <- ddply(data, .(), transform, truePVol = ifelse(shape == "sphere", sphereVol(size1) / sphereVol(size2) * 100, trueP))



# compute the accuracy, i.e., the difference between the estimate and true ratio  

data  <- ddply(data, .(), transform, accuracySur = abs(estimate - truePSur) )
data  <- ddply(data, .(), transform, accuracyVol = abs(estimate - truePVol) )

data  <- ddply(data, .(), transform, rawErrorSur = estimate - truePSur)
data  <- ddply(data, .(), transform, rawErrorVol = estimate - truePVol)


# now we compute bootstrapped confidence intervals for the accuracies as we did above

ord <- max(accuracies$order)

cis <- bootstrapCI(subset(data, shape == "sphere" & method == "ratio estimation")$accuracySur)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (RE)", measure="surface", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order = inc(ord)))




cis <- bootstrapCI(subset(data, shape == "sphere" & method == "constant sum")$accuracySur)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (CS)", measure="surface", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order = inc(ord)))



cis <- bootstrapCI(subset(data, shape == "sphere" & method == "constant sum")$accuracyVol)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (CS)", measure="volume", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order = inc(ord)))

cis <- bootstrapCI(subset(data, shape == "sphere" & method == "ratio estimation")$accuracyVol)

accuracies  <- rbind(accuracies, data.frame(pub="Jansen [2016]", shape="sphere (RE)", measure="volume", estimate=cis[1], ci.low=cis[2], ci.high=cis[3], order = inc(ord)))


cidf <- rbind(prevD, accuracies)

dotData <- data.frame(x = rep(0:20, times=nrow(cidf)), y = rep(paste(cidf$shape, cidf$measure), each=21), order = rep(cidf$order, each = 21))

qplot(data=dotData, reorder(y, order), x, alpha=I(0.2)) + expand_limits(x=c(1,5), y=c(0,20)) +
  mytheme +
  coord_flip() +
  geom_pointrange(data=cidf, aes(y=estimate, x=reorder(paste(shape, measure), order), ymin=ci.low, ymax=ci.high), width=1, size=0.7) +
  labs(x = "", y = "Average accuracy (absolute discrepancy in percent)") +
   annotate("text", x= 3.2, y = 16, label = cidf$pub[1]) +
   annotate("segment", x=0.5, xend=5.5, y = 12, yend = 12) +
   annotate("segment", x=5.5, xend=5.5, y = 11.5, yend = 12) +
   annotate("segment", x=0.5, xend=0.5, y = 11.5, yend = 12)
Dot plot showing accuracies to compare different 'true' ratios between spheres similar to Figure 11 in the paper

As we can see in the figure above, people's accuracy is a lot better when we calculate the ratio between two spheres using their surface areas. In comparison to the accuracy based on volume ratios, basing it on surface area cuts the average error in half, at least when using ratio estimation.

Variance Between People - Individual Differences

Next we investigate individual differences. To get a first overview, we plot data for each participant separately. 

#sData <- ddply(data, .(), transform, estimate = estimate / 100.0, trueP = trueP / 100.0)

coeffs <- ddply(data, .(method, shape, participant), function(dat){
  df <- data.frame(x = dat$trueP, y = dat$estimate)
  #print(head(df))
  regSN <- nls(y ~ 100 * (x/100)^a, data = df, start = list(a=1.0), control = list(maxiter = 500))
  coefSN <- coef(regSN)
  a <- coefSN[[1]]

  # compute the residual standard error as a goodness-of-fit metric
  errorfn <- function(reg) sqrt(deviance(reg)/df.residual(reg))

  # compute bootstrapped confidence intervals for the regression
  bootsJ <- nlsJack(regSN)
  # boots <- nlsBoot(regSN, niter=500)
  # get ci
  a_lower <- mean(bootsJ$jackCI$Low)
  a_upper <- mean (bootsJ$jackCI$Up)

  data.frame(participant=dat$participant[1], method=dat$method[1], shape=dat$shape[1], coef = a, coef.low = a_lower, coef.high = a_upper, regError = errorfn(regSN))
})

# add the coefficients as a column in our main data file
data <- merge(data, coeffs)

# compute lci and hci values for each ratio from the coeffs
data  <- ddply(data, .(), transform, ci.low = (trueP/100)^coef.low, ci.high = (trueP/100)^coef.high)
regFun <- function(x, a) {100*(x/100)^a}
plots <- dlply(data, .(method, shape, participant), function(dat){

p <- ggplot(data=dat, aes(x=trueP, y=estimate)) +
  stat_function(fun=regFun, args=list(a = dat$coef[1]), col="red") +
  #geom_smooth(method="nls2", formula = y ~ (x/100)^a, #dat$estimate ~ 100*((dat$trueP/100)^a), 
   #           data = dat, start = list(a=1), control = list(maxiter = 500), aes(ymin = (pred$fit)^coef.low, ymax = (pred$fit)^coef.high)) +
  geom_point()  +
  labs(title = dat$participant, x = "", y = "") +
  mytheme + myaxes + theme(plot.margin = unit(c(0,0,0,0), "lines"))
return(p)
})
label = textGrob("Bars (constant sum)", just = c("centre", "top"), gp = gpar(fontsize = 14))
label2 = textGrob("Spheres (constant sum)", just = c("centre", "top"), gp = gpar(fontsize = 14))
label3 = textGrob("Bars (ratio estimation)", just = c("centre", "top"), gp = gpar(fontsize = 14))
label4 = textGrob("Spheres (ratio estimation)", just = c("centre", "top"), gp = gpar(fontsize = 14))

blank<-rectGrob(gp=gpar(col="white")) # make a white spacer grob

grid.arrange(do.call("arrangeGrob", c(plots[21:30], ncol = 5, top="Bars (ratio estimation)", left="estimated ratio (in percent)", bottom="ratio between heights (in percent)")),
            blank, do.call("arrangeGrob", c(plots[31:40], ncol = 5, top="Spheres (ratio estimation)", left="estimated ratio (in percent)", bottom="ratio between heights (in percent)")),
            blank, do.call("arrangeGrob", c(plots[1:10], ncol = 5, top="Bars (constant sum)", left="estimated ratio (in percent)", bottom="ratio between heights (in percent)")),
            blank, do.call("arrangeGrob", c(plots[11:20], ncol = 5, top="Spheres (constant sum)", left="estimated ratio (in percent)", bottom="ratio between heights (in percent)")), heights=c(4, 0.7, 4, 0.7, 4, 0.7, 4),
            ncol=1)
plot of chunk compute-coefficients

The individual curces indicate that, for estimates made using the ratio estimation method, the nonlinear regression assuming a exponential function, fits rather well. For estimates made using the constant sum method however, some curves are a bad fit for our obtained data. Since we were unaware of this phenomenon at the time of the experiment, and even expected the method to outperform ratio estimation, we were unable to collect further information on this issue from our participants. Since for most people the individual judgments still seem to lie on an approximative curve, we assume that the mental models of participants for how to interpret the task varied immensely. Further studies will be necessary to investigate this question further. Next, to be able to better compare the individually fit curves, we now overplot the regression curves similarly to figure 10 in our article. 

par(mfrow=c(2,2), pty="s")
# p <- plot()
minA <- 100
maxA <- 0
curve((x), xlim=c(0,100), ylim=c(0,100), col="red", ylab="predicted estimated ratio (in percent)", xlab="actual ratio (in percent)", main="Spheres - ratio estimation")
d_ply(subset(data, shape=="sphere" & method== "ratio estimation"), .(participant), function(d){
  curve(regFun(x,a=d$coef[1]), add = TRUE, col=rgb(0,0,0,0.5))
  #p  <- p + geom_smooth(data=d, aes(x = d$trueP, y = d$estimate), method="loess", alpha=I(0.1))
})
  text(x=c(20,80), y=c(70, 30), labels = c(prettyNum(min(subset(data, shape=="sphere" & method== "ratio estimation")$coef), digits=2), prettyNum(max(subset(data, shape=="sphere" & method== "ratio estimation")$coef), digits=2)))

curve((x), xlim=c(0,100), ylim=c(0,100), col="red", ylab="predicted estimated ratio (in percent)", xlab="actual ratio (in percent)", , main="Spheres - constant sum")
d_ply(subset(data, shape=="sphere" & method== "constant sum"), .(participant), function(d){
  curve(regFun(x,a=d$coef[1]), add = TRUE, col=rgb(0,0,0,0.5))
  #p  <- p + geom_smooth(data=d, aes(x = d$trueP, y = d$estimate), method="loess", alpha=I(0.1))
})
  text(x=c(20,80), y=c(70, 30), labels = c(prettyNum(min(subset(data, shape=="sphere" & method== "constant sum")$coef), digits=2), prettyNum(max(subset(data, shape=="sphere" & method== "constant sum")$coef), digits=2)))


curve((x), xlim=c(0,100), ylim=c(0,100), col="red", ylab="predicted estimated ratio (in percent)", xlab="actual ratio (in percent)", main="Bars - ratio estimation")
d_ply(subset(data, shape=="bar" & method== "ratio estimation"), .(participant), function(d){
  curve(regFun(x,a=d$coef[1]), add = TRUE, col=rgb(0,0,0,0.5))
  #p  <- p + geom_smooth(data=d, aes(x = d$trueP, y = d$estimate), method="loess", alpha=I(0.1))
})
  text(x=c(20,80), y=c(70, 30), labels = c(prettyNum(min(subset(data, shape=="bar" & method== "ratio estimation")$coef), digits=2), prettyNum(max(subset(data, shape=="bar" & method== "ratio estimation")$coef), digits=2)))


curve((x), xlim=c(0,100), ylim=c(0,100), col="red", ylab="predicted estimated ratio (in percent)", xlab="actual ratio (in percent)", main="Bars - constant sum")
d_ply(subset(data, shape=="bar" & method== "constant sum"), .(participant), function(d){
  curve(regFun(x,a=d$coef[1]), add = TRUE, col=rgb(0,0,0,0.5))
  #p  <- p + geom_smooth(data=d, aes(x = d$trueP, y = d$estimate), method="loess", alpha=I(0.1))
})
  text(x=c(20,80), y=c(70, 30), labels = c(prettyNum(min(subset(data, shape=="bar" & method== "constant sum")$coef), digits=2), prettyNum(max(subset(data, shape=="bar" & method== "constant sum")$coef), digits=2)))
Overplotted individual regression curves similar to Figure 10 in the paper
par(mfrow=c(1,1), pty="s")

The overplotted curves further confirm that variations between participants are much higher when comparing the constant sum estimations. The fact that this can be observed for both bars and spheres confirms that these findings are likely due to an "implementation" problem of the method. While we initially assumed that a visual response format - the constant sum method - should be more intuitive than one that requires a conversion to numbers, our results show otherwise. However, previous studies, such as Spence (1990) report high accuracies using this method. This discrepancy might be due to Spence having participants go through a more extensive training phase and requiring them to repeat the training phase until they performed sufficiently well. We enforced no such protocol and had participants continue with the experiment when they reported understanding the task after two training stimuli. Since we could only speculate on possible reasons for our data concerning the constant sum method, we focus from here on our analysis on estimates made using the ratio estimation method.

Residuals

Next, we investigate the question whether people make consistent judgments across the range of different ratios. If some ratios, for example, small ones, are judged quite accuractely while others show large errors, for example large ones, then the tested shape would be rarely suitable to encode data. Similarly, overestimations at end end of the range and underestimations on the other end of the scale would render a shape equally of little use. We therefore now analyze residuals.

volResidualPlots <- dlply(subset(data, method=="ratio estimation"), .(shape), function(d){
   p  <-  ggplot(data=d, aes(x=truePVol, y = rawErrorVol)) + geom_point(alpha=I(0.2))  +
  labs(title = " ", x = " ", y = "residuals") +
  mytheme + coord_fixed(ratio = 1, xlim = c(0,100), ylim = c(-50,50))
  return(p)


})

volResidualHist <- dlply(subset(data, method=="ratio estimation"), .(method, shape), function(d){
  p  <- qplot(data=d, rawErrorVol, geom="histogram", binwidth = 1)  + coord_equal(ratio = 6) + coord_flip()  + mytheme  + scale_y_continuous(limits = c(0, 30), expand=c(0,0)) + scale_x_continuous(limits = c(-50, 50)) + geom_vline(xintercept=0, col="red") +
    theme(
    panel.margin = element_blank(),
    axis.ticks = element_blank(),
    panel.border = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x  = element_blank(),
    axis.text.y = element_blank(),
    panel.margin.x = NULL,
    panel.margin.y = NULL,
    panel.border = element_blank()) + labs(main=" ", x=" ")
  return(p)
})

idx <- order(c(seq_along(volResidualPlots), seq_along(volResidualHist)))
mergedplots <- (c(volResidualPlots, volResidualHist))[idx]



grid.arrange(do.call("arrangeGrob", c(mergedplots[1:2], list(ncol = 2, widths=c(4,1.5)))), top = "Bars (ratio estimation)", bottom = "ratios between heights")
plot of chunk residual_bars
grid.arrange(do.call("arrangeGrob", c(mergedplots[3:4], list(ncol = 2, widths=c(4,1.5)))), top = "Spheres (ratio estimation)", bottom = "ratios between volumes")
plot of chunk raw_error_spheres

The residuals for estimates of ratios between heights of bars are centered around the origin (marked in red). For spheres, we find consistent overestimations if we would encode data in their volume. As our prior analyses of accuracy indicated that error rates are lower when encoding data in other measures, we now investigate residuals for diameters, surfaces, and a globally fitted regression. 

# diameter
diaResidualPlot <- ggplot(data=subset(data, method=="ratio estimation" & shape == "sphere"), aes(x=trueP, y = rawError)) + geom_point(alpha=I(0.2))  +
  labs(title = " ", x = " ", y = "residuals") +
  mytheme + coord_fixed(ratio = 1, xlim = c(0,100), ylim = c(-50,50))


diaResidualHist <- qplot(data=subset(data, method=="ratio estimation" & shape == "sphere"), rawError, geom="histogram", binwidth = 1)  + coord_equal(ratio = 6) + coord_flip()  + mytheme  + scale_y_continuous(limits = c(0, 30), expand=c(0,0)) + scale_x_continuous(limits = c(-50, 50)) + geom_vline(xintercept=0, col="red") +
    theme(
    panel.margin = element_blank(),
    axis.ticks = element_blank(),
    panel.border = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x  = element_blank(),
    axis.text.y = element_blank(),
    panel.margin.x = NULL,
    panel.margin.y = NULL,
    panel.border = element_blank()) + labs(main=" ", x=" ")


mergedplots <- list(diaResidualPlot, diaResidualHist)


grid.arrange(do.call("arrangeGrob", c(mergedplots, list(ncol = 2, widths=c(4,1.5)))), top = "Spheres (ratio estimation)", bottom = "ratios between diameters")
plot of chunk residuals_diameter

While ratios between diameters lead to lower error rates than ratios between volumes, the distribution of residuals is still unsatisfying. 

# surfaces
surResidualPlot <- ggplot(data=subset(data, method=="ratio estimation" & shape == "sphere"), aes(x=truePSur, y = rawErrorSur)) + geom_point(alpha=I(0.2))  +
  labs(title = " ", x = " ", y = "residuals") +
  mytheme + coord_fixed(ratio = 1, xlim = c(0,100), ylim = c(-50,50))


surResidualHist <- qplot(data=subset(data, method=="ratio estimation" & shape == "sphere"), rawErrorSur, geom="histogram", binwidth = 1)  + coord_equal(ratio = 6) + coord_flip()  + mytheme  + scale_y_continuous(limits = c(0, 30), expand=c(0,0)) + scale_x_continuous(limits = c(-50, 50)) + geom_vline(xintercept=0, col="red") +
    theme(
    panel.margin = element_blank(),
    axis.ticks = element_blank(),
    panel.border = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x  = element_blank(),
    axis.text.y = element_blank(),
    panel.margin.x = NULL,
    panel.margin.y = NULL,
    panel.border = element_blank()) + labs(main=" ", x=" ")


mergedplots <- list(surResidualPlot, surResidualHist)


grid.arrange(do.call("arrangeGrob", c(mergedplots, list(ncol = 2, widths=c(4,1.5)))), top = "Spheres (ratio estimation)", bottom = "ratios between surface areas")
plot of chunk residuals_surfaces

As we saw earlier, ratios based on surface areas cut error rates in half. However, they still lead overall to overestimations as residuals are still not centered around the origin.

As a last step we now plot residuals for our fitted residual. 

dat <- subset(data, method == "ratio estimation" & shape == "sphere")
  df <- data.frame(x = dat$trueP, y = dat$estimate)
  regSN <- nls(y ~ 100 * (x/100)^a, data = df, start = list(a=1.0), control = list(maxiter = 500))

  regs <- nlsResiduals(regSN)


names(regs) <- gsub(" ", ".", names(regs))

  resi  <-  as.data.frame(regs$resi1)
  names(resi)[1] <- "Fitted"
  #print(head(resi))
  fitResidualPlot  <-  ggplot(data=resi, aes_string(x="Fitted", y = "Residuals")) + geom_point(alpha=I(0.3))  +
  labs(title = "", x = "fitted ratios", y = "residuals") +
  mytheme + coord_fixed(ratio = 1, xlim = c(0,100), ylim = c(-50,50))


fitResidualHist <- qplot(data=resi, Residuals, geom="histogram", binwidth = 1)  + coord_equal(ratio = 6) + coord_flip()  + mytheme  + scale_y_continuous(limits = c(0, 30), expand=c(0,0)) + scale_x_continuous(limits = c(-50, 50)) + geom_vline(xintercept=0, col="red") +
    theme(
    panel.margin = element_blank(),
    axis.ticks = element_blank(),
    panel.border = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.x  = element_blank(),
    axis.text.y = element_blank(),
    panel.margin.x = NULL,
    panel.margin.y = NULL,
    panel.border = element_blank()) + labs(main=" ", x=" ")


mergedplots <- list(fitResidualPlot, fitResidualHist)

grid.arrange(do.call("arrangeGrob", c(mergedplots, list(ncol = 2, widths=c(4,1.5)))), top = "Spheres (ratio estimation)", bottom = "ratios according to fitted regression")
plot of chunk residuals_fitted

This last plot indicates that a perceptual correction for spheres might be possible. Error rates are still twice as high as with bar charts but the profound misperception of the encoded data observed when assuming ratios between volumes, is corrected.

However, further studies will be necessary to determine how including manual exploration affects people's estimates. Previous work on haptic size perception of physical shapes [Kahrimanovic 2010] indicates that people's estimates are there also close to a surface area based scale. It thus might be that encoding data in surface area would minimize errors when both visual and manual exploration of shapes is possible.

Kahrimanovic, M., Tiest, W. M. B., & Kappers, A. M. L. (2010). Haptic perception of volume and surface area of 3-D objects. Attention, Perception & Psychophysics, 72(2), 517–527.

Kahrimanovic, M., Bergmann Tiest, W. M., & Kappers, A. M. L. (2010). Seeing and feeling volumes: The influence of shape on volume perception. Acta Psychologica, 134(3), 385–390.