---
title: "Guidelines for Composing Graphical Interfaces for Fast and Correct Perception of System State:  Data Analysis"
output: html_notebook
---


```{r echo=FALSE}
library(lme4)
library(dplyr)
library(ggplot2)
library(tidyr)
```
#Load data

Data was recorded with a dedicated software that showed the situations for the different designs. Data of all participants was afterwards concatenated into one file. 
Meta information about anonymized subjects was recorded in separated file. 
Load data and merge with subject meta data:

```{r}

# read dataset
data <- read.csv("DataKonectGuidelines_allSubjects.csv", header = T, sep = ",",stringsAsFactors=FALSE)

# read subject information
subjects <- read.csv("SubjectsKonectGuidelines.csv", header = T, sep = ",",stringsAsFactors=FALSE)

#join data with meta information about subjects. VPID is the identifier for subjects.
data <- merge(data, subjects, by="VPID")

#For simplification in the paper we use IDs "D1" - "D6" for the design concepts, instead of the HMI designer group ID, that were used in the workshops
data$BaseDesign <- mapvalues(data$BaseDesign, from=c('P1P2','P3P4','P5P6','P7P8','P9P10','P13P14') , to=c('D1','D2','D3','D4','D5','D6'))
#We use more readable situation desciptors 
data$DesignState <- mapvalues(data$DesignState, from=c('ok', 'lowSpeed', 'criticalTrim', 'relevantDuK', 'strongWind') , to=c('non-critical','critical ship speed','critical trim state','critical (low) depth under keel','critical wind condition'))
```

#Data selection.#

Dataset contains a trace of all events in the recording software. Reaction times and correctness of responses are recorded in datasample with "action = REACTION"

```{r}
data <- filter(data, action == 'REACTION')
```

Dataset is from a larger experiment adressing multiple issues. For analysing the effect of guidelines on reaction times and accuracy only data from subjects that conducted A-B-A-B blocks are of interest.

```{r}
data <- filter(data, experimentScenario == 'ABAB')
```
Block A and B were iterated twice. To eliminate strong training effects we only analyse data from the second iteration

```{r}
data <- filter(data, BlockRepetition == 2)

```
Cutoffs for reaction times: 
We consider reaction times over 10 seconds as failure. 
Reaction times include perception, decision and motor response. 
It is unlikely that this happens in less than 300 ms, which we consider a click-error. 
Reaction times exceeding these limits were excluded

```{r}
excluded <- filter(data, reaction_ms > 10000 | reaction_ms < 300)
# number of excluded samples
nrow(excluded)
# exclude according to cutoff times
dsta <- filter(data, reaction_ms <= 10000 & reaction_ms >= 300)

```


Relevant information for subsequent analysis are 
* VPID: Identifier for the subject from whom this data sample is, 
* critical: criticality of shown situation, 
* reaction_ms: reaction time in ms, 
* Guideline: whether the currently shown situation is from an initial design (Guideline==0) or from a design after applying the guidelines (Guideline == 1)
* BaseDesign: Identifier for the design concept. The same for the initial deisgn and the guideline design for the same design team.
* DesignState: currently shown situation. Referenced as S1 - S5 in the paper.
* block: Block A or B. See paper
* stimuli: Order of stimuli. In each block 108 stimuli were shown.
* correct: 1, if response of subject was correct


```{r}
data <- select (data, VPID, critical, reaction_ms, Guideline, BaseDesign, DesignState, block, stimuli, correct)
head(data)
```

# Analysis of reaction times
## Visualize training effect
Visualize training effect. Only for dataset with initial designs (guideline == 0) as an example.

```{r}
plotdata <- filter(data, Guideline == 0)
plot(plotdata$stimuli, plotdata$reaction_ms)
abline(lm(plotdata$reaction_ms ~ plotdata$stimuli))
summary(lm(plotdata$reaction_ms ~ plotdata$stimuli))
```
## Effect of guidelines
Fit mixed effects model using maximum likelihood approach. See paper for model argumentation.

```{r}
data.RT.model <- lmer(reaction_ms ~  Guideline + stimuli + (1|VPID) + (1|BaseDesign), data=data, REML=FALSE)
summary(data.RT.model)
```

Guidelines on average decreased reaction times (`r data.RT.model@beta[2]` ms)


Fit mixed effects model **without** guidelines as explanatory variable using maximum likelihood approach
```{r}
data.RT.nullGuideline <- lmer(reaction_ms ~ stimuli + (1|VPID) + (1|BaseDesign), data=data, REML=FALSE)
summary(data.RT.nullGuideline)

```

Use likelihood ratio to estimate significance of guideline usage
```{r}
data.RT.mGuideline <- anova(data.RT.nullGuideline, data.RT.model)
data.RT.mGuideline
```
Model fit improves with inclusion of guidelines as explanatory variable. Effect is significant ($\chi^2$(`r data.RT.mGuideline$'Chi Df'[2]`)=`r data.RT.mGuideline$Chisq[2]`, p=`r data.RT.mGuideline$'Pr(>Chisq)'`). 


### Random slopes for each design

Assumption is, that the effect of the guidelines strongly differs for each design.
Adding random slopes to the model:

``` {r}
#fit data
data.RT.modelRandomSlopes <- lmer(reaction_ms ~ Guideline + stimuli + (1|VPID) + (Guideline|BaseDesign), data=data, REML=FALSE)
# print summary
summary(data.RT.modelRandomSlopes)
# print parameters (including slopes for designs)
coef(data.RT.modelRandomSlopes)

# Likelihood ratio test with original model
data.RT.mSlopesGuideline <- anova(data.RT.model, data.RT.modelRandomSlopes)
data.RT.mSlopesGuideline
```
Adding random slopes for each design significantly improves model fit ($\chi^2$(`r data.RT.mSlopesGuideline$'Chi Df'[2]`)=`r data.RT.mSlopesGuideline$Chisq[2]`, p=`r data.RT.mSlopesGuideline$'Pr(>Chisq)'`). The effect of guidelines on reaction time improvement strongly differs between designs

Visualisation of random slopes:
```{r}
#Mean reaction time values for each design
dSlopes <- ddply(data, .(BaseDesign, Guideline), summarize,  reaction_ms=mean(reaction_ms))

# Plot slope for each design
plot(1,type='n',xlim=c(-0.5,1.5),ylim=c(900,1600),xlab='Guideline usage', ylab='reaction time', xaxt='n')
axis(1,at=seq(0,1), labels=c("initial design", "with guidelines"))
i=1
for (design in unique(dSlopes$BaseDesign)){
    i=i+1
    d <-dSlopes[which(dSlopes$BaseDesign==design), ]
    lines(d$Guideline, d$reaction_ms, col=i, type='o', lwd=2)
}
xrange <- range(dSlopes$Guideline)
yrange <- range(dSlopes$reaction_ms)

legend(1,1600,unique(dSlopes$BaseDesign), col=2:7, lty=1)


# For paper figures: Transform into wide format and use gnuplot
# Gnuplot figure is created by running 'design_slopes_reaction_times.gnuplot' with gnuplot

dSlopes$Guideline <- factor(dSlopes$Guideline)
dSlopesGnuplot <- spread(dSlopes, BaseDesign, reaction_ms)
dSlopesGnuplot$GuidelineLabel <- mapvalues(dSlopesGnuplot$Guideline, from=c('0','1') , to=c('initial','guideline'))

write.csv(dSlopesGnuplot, "design_slopes_reaction_times.csv")
```




## Effect of training
Fit mixed effects model **without stimulus order** as explanatory variable using maximum likelihood approach. And compare with full model.

```{r}
data.RT.nullStimuli <- lmer(reaction_ms ~ Guideline + (1|VPID) + (1|BaseDesign), data=data, REML=FALSE)
summary(data.RT.nullStimuli)
data.RT.mStimuli <- anova(data.RT.nullStimuli, data.RT.model)
data.RT.mStimuli
```

Model fit improves with inclusion of stimulus order as explanatory variable. Reaction time decreases on average by `r data.RT.model@beta[3]` ms per stimulus ($\hat{=}$`r (data.RT.model@beta[3] * max(data$stimuli))` ms reduction from first to last stimulus).
This training effect is significant ($\chi^2$(`r data.RT.mStimuli$'Chi Df'[2]`)=`r data.RT.mStimuli$Chisq[2]`, p=`r data.RT.mStimuli$'Pr(>Chisq)'`). 

# Analysis of accuracy of responses
## Effect of guidelines

Fit mixed effects model using maximum likelihood approach. See paper for model argumentation.

```{r}
data.correct.model <- lmer(correct ~ Guideline + (1|VPID) + (1|BaseDesign), data=data, REML=FALSE)
summary(data.correct.model)
```
Guidelines on average improve accuracy by `r data.correct.model@beta[2] * 100` \% (from `r data.correct.model@beta[1] * 100` \% to `r (data.correct.model@beta[1]+data.correct.model@beta[2])*100` \%).

Fit mixed effects model **without** guidelines as explanatory variable using maximum likelihood approach.
```{r}
data.correct.nullGuideline <- lmer(correct ~ (1|VPID) + (1|BaseDesign), data=data, REML=FALSE)
summary(data.correct.nullGuideline)
```

Use likelihood ratio to estimate significance of guideline usage
```{r}
data.correct.mGuideline <- anova(data.correct.nullGuideline, data.correct.model)
data.correct.mGuideline
```
Model fit improves with inclusion of guidelines as explanatory variable. Effect is significant ($\chi^2$(`r data.correct.mGuideline$'Chi Df'[2]`)=`r data.correct.mGuideline$Chisq[2]`, p=`r data.correct.mGuideline$'Pr(>Chisq)'`). 


### Random slopes for each design
Assumption is, that the effect of the guidelines strongly differs for each design.
Adding random slopes to the model:
``` {r}
# fit mixed effects model with different accuracy-slopes for each design
data.correct.modelRandomSlopes <- lmer(correct ~ Guideline + (1|VPID) + (Guideline|BaseDesign), data=data, REML=FALSE)
# print summary
summary(data.correct.modelRandomSlopes)
# print parameters (including slopes for designs)
coef(data.correct.modelRandomSlopes)


# Likelihood ratio test with original model
data.correct.mSlopesGuideline <- anova(data.correct.model, data.correct.modelRandomSlopes)
data.correct.mSlopesGuideline
```
Adding random slopes for each design significantly improves model fit ($\chi^2$(`r data.correct.mSlopesGuideline$'Chi Df'[2]`)=`r data.correct.mSlopesGuideline$Chisq[2]`, p=`r data.correct.mSlopesGuideline$'Pr(>Chisq)'`). The effect of guidelines on accuracy  differs between designs.

Visualisation of random slopes:

```{r}
# Mean accuracy values for each deisgn converted to percentage
dSlopes <- ddply(data, .(BaseDesign, Guideline), summarize,  accuracy=mean(correct)*100)

# Plot slope for each design
plot(1,type='n',xlim=c(-0.5,1.5),ylim=c(80,100),xlab='Guideline usage', ylab='accuracy', xaxt='n')
axis(1,at=seq(0,1), labels=c("initial design", "with guidelines"))
i=1
for (design in unique(dSlopes$BaseDesign)){
    i=i+1
    d <-dSlopes[which(dSlopes$BaseDesign==design), ]
    lines(d$Guideline, d$accuracy, col=i, type='o', lwd=2)
}
xrange <- range(dSlopes$Guideline)
yrange <- range(dSlopes$accuracy)

legend(1.1,100,unique(dSlopes$BaseDesign), col=2:7, lty=1)

# For paper figures: Transform into wide format and use gnuplot
# Gnuplot figure is created by running 'design_slopes_accuracy.gnuplot' with gnuplot
dSlopes$Guideline <- factor(dSlopes$Guideline)
dSlopesGnuplot <- spread(dSlopes, BaseDesign, accuracy)
dSlopesGnuplot$GuidelineLabel <- mapvalues(dSlopesGnuplot$Guideline, from=c('0','1') , to=c('initial','guideline'))

dSlopesGnuplot
write.csv(dSlopesGnuplot, paste(processing_dir, "design_slopes_accuracy.csv", sep=""))
```

## Effect of training

We didn't expect to observe a training effect for accuracy, because the vast majority of responses is already correct for the initial designs. Room for improvement is small and harder to measure with binary measurement. If we add stimulus order and compare model fit of our model with this extended one, we get:
```{r}
data.correct.modelStimuli <- lmer(correct ~ Guideline + stimuli + (1|VPID) + (1|BaseDesign), data=data, REML=FALSE)
summary(data.correct.modelStimuli)
data.correct.mStimuli <- anova(data.correct.model, data.correct.modelStimuli)
data.correct.mStimuli
```
There seems to be a small tendency for a training effect (accuracy improves by `r data.correct.modelStimuli@beta[3] *100` \% per stimulus $\hat{=}$ `r (data.correct.modelStimuli@beta[3] * max(data$stimuli) * 100)` \% improvement from first to last stimulus), but it is not significant ($\chi^2$(`r data.correct.mStimuli$'Chi Df'[2]`)=`r data.correct.mStimuli$Chisq[2]`, p=`r data.correct.mStimuli$'Pr(>Chisq)'`). 


# Exploration
## Slopes for each situation of each design
Visualize how close the changes of reaction times are for the different situations in each design. Group by designs to give an impression of closeness/distance.
```{r}
# Mean reaction timess for each situation of each design
dSituationSlopes <- ddply(data, .(BaseDesign, Guideline, DesignState), summarize,  reaction_ms=mean(reaction_ms))

# Select values from initial designs and guideline designs
initialData <- filter(dSituationSlopes, Guideline==0)
guidelineData <- filter(dSituationSlopes, Guideline==1)
# Join both values for each situation of each design and compare them
dSituationSlopes <- merge(initialData, guidelineData, by=c('BaseDesign','DesignState'))
dSituationSlopes$reaction_ms_Slope <- dSituationSlopes$reaction_ms.x - dSituationSlopes$reaction_ms.y

#Exclude irrelevant columns
dSituationSlopes <- select(dSituationSlopes, BaseDesign, DesignState, reaction_ms_Slope)
head(dSituationSlopes)

# Plot distribution of slopes for each design
dSituationSlopes$BaseDesign <- factor(dSituationSlopes$BaseDesign)
plot(dSituationSlopes$BaseDesign, dSituationSlopes$reaction_ms_Slope)

# Figure from paper is created with gnuplot from this data. For gnuplot transform data into wide format and assign numeric IDs to string descriptors of BaseDesign. 
dSituationSlopes$DesignState <- factor(dSituationSlopes$DesignState)
dSituationSlopes <- spread(dSituationSlopes, DesignState, reaction_ms_Slope)
dSituationSlopes$designIndex <- mapvalues(dSituationSlopes$BaseDesign, from=c('D2','D3','D4','D5','D6','D1') , to=c(1,2,3,4,5,6))

#Then write data to file. Run 'design_situation_slopes.gnuplot' manually with gnuplot to recompile figure.
write.csv(dSituationSlopes, "design_situation_slopes.csv", quote = FALSE)
```
