# Analyses for The Attentional Cost of Receiving a Cell Phone Notification

# Author: Cary Stothart
# Date: 05/28/2015

library(lme4)
library(ggplot2)
library(scales)
library(extrafont)
library(plyr)
library(AICcmodavg)

# List of participants who correctly guessed manipulation/hypothesis.
part_know <- c(146)

# Load the trial- and participant-level datasets and create merged dataset.
trial.only.data<- read.delim("trial_level_data.txt", sep="\t", header=TRUE)
part.only.data <- read.delim("part_level_data.txt", sep="\t", header=TRUE)
colnames(part.only.data) <- tolower(colnames(part.only.data))
trial.data <- merge(trial.only.data, part.only.data, by="part_num")

# Remove unusable participants and trials.
trial.data <- trial.data[trial.data$bad_data==0, ] # Remove bad data (see data file
# for the reasons for removal).
trial.data <- trial.data[trial.data$block!=0, ] # Remove practice block trials.
trial.data <- trial.data[trial.data$ph_current_setting!="Off", ] # Remove participants
# who had their phone off.
trial.data <- trial.data[!trial.data$part_num %in% part_know, ] # Remove participants
# who had knowledge of the hypothesis.

# Create variables that are needed for the analyses.
trial.data$omit_trial <- 0
trial.data$omit_trial[trial.data$number==3] <- 1
trial.data$com_error <- 0
trial.data$com_error[trial.data$resp_acc==0 & trial.data$omit_trial==1] <- 1
trial.data$block_num <- factor(trial.data$block_num)
trial.data$part_num <- factor(trial.data$part_num)
trial.data$com_trial <- 0
trial.data$com_trial[trial.data$number!=3] <- 1
trial.data$omit_error <- 0
trial.data$omit_error[trial.data$resp_acc==0 & trial.data$com_trial==1] <- 1

# At this point, trial.data is the trial-level dataset with all exclusions applied  
# (see above). The dataset includes both omission and commission trials. Also contains 
# both trial-level and participant-level data.

# Fit the commission error data to a mixed effects logistic regression using group
# block number and their interaction as predictors.  Only uses omission trials.
omit.data <- trial.data[trial.data$omit_trial==1, ]
contrasts(omit.data$block_num) <- c(-1, 1)
alertsVsNoAlert <- c(1, -2, 1)
callVsText <- c(1, 0, -1)
contrasts(omit.data$group) <- cbind(alertsVsNoAlert, callVsText)
a.fit1 <- glmer(com_error~(1|part_num), family=binomial,
                data=omit.data)
a.fit2a <- glmer(com_error~block_num+(1|part_num), family=binomial,
                 data=omit.data)
a.fit2b <- glmer(com_error~group+(1|part_num), family=binomial,
                 data=omit.data)
a.fit3 <- glmer(com_error~block_num+group+(1|part_num), family=binomial,
                data=omit.data)
a.fit4 <- glmer(com_error~block_num*group+(1|part_num), family=binomial,
                data=omit.data)
anova(a.fit1, a.fit2a) # Block main effect.
anova(a.fit1, a.fit2b) # Group main effect.
anova(a.fit3, a.fit4) # Group by block interaction.
summary(a.fit4) # Group by block interaction contrasts.
#confint(a.fit4) # 95% confidence intervals for the coefficients (takes awhile to execute).
hist(ranef(a.fit4)$part_num[, 1]) # Check if intercepts are normally distributed.

# Theme for bar plots
font_import(pattern="[A/a]rial")
y
loadfonts(device="win")
font <- "Arial"
plot.theme <- theme(panel.background=element_rect(fill="#ffffff", color="#ffffff"),
                    axis.title.x=element_text(face="plain", size=14, color="#000000",
                                              vjust=0, family=font),                                              
                    axis.title.y=element_text(face="plain", size=14, color="#000000",
                                              vjust=1, family=font),
                    axis.text.x=element_text(face="plain", size=14, color="#000000",
                                             family=font),
                    axis.text.y=element_text(face="plain", size=14, color="#000000",
                                             family=font),
                    legend.title=element_text(face="plain", size=14, color="#000000",
                                              family=font),                     
                    legend.text=element_text(face="plain", size=14, color="#000000",
                                             family=font), 
                    axis.ticks=element_blank(),
                    panel.grid.minor=element_blank(),
                    panel.grid.major=element_blank())

# Create a bar graph displaying the probability of making a commission error for
# each cell.
plot1.data <- expand.grid(block_num=c("1", "2"), group=unique(omit.data$group))
plot1.data$block_num <- factor(plot1.data$block_num)
prob.and.se <- predictSE(a.fit4, newdata=plot1.data, type="response")
plot1.data$prob = prob.and.se$fit
plot1.data$prob_se = prob.and.se$se.fit
plot1.data$group2 <- c("Call", "Call", "Text message", "Text message",
                       "No notification", "No notification")
plot1.data$block2 <- rep(c("Block 1", "Block 2"), 3)
plot1 <- ggplot(plot1.data, aes(x=group2, y=prob)) +
  geom_bar(aes(fill=block2), position=position_dodge(0.9), stat = "identity") +
  geom_errorbar(aes(ymin=prob-prob_se, ymax=prob+prob_se, group=block2), 
                position=position_dodge(0.9), stat = "identity", width=0.5) +
  ylab("Probability of making a commission error") +
  xlab("") +
  scale_y_continuous(labels=percent, expand=c(0, 0)) +
  scale_x_discrete(expand=c(0, 0)) +  
  scale_fill_manual("", values=c("#12006f", "#8f9a9d")) +  
  plot.theme
plot1
ggsave(filename="Figure2_ComErrorProbPlot.tiff", dpi=600)

# Assess cell differences in mean number of commission errors made. Calculate
# Cohen's Dz for each block difference for each group.
omit.block.data <- ddply(omit.data, c("part_num", "block_num"), summarize,
                         part_num=part_num[1], group=group[1], block_num=block_num[1],
                         total_com_errors=sum(com_error))
cell.com.data <- data.frame(aggregate(total_com_errors~block_num:group, 
                                      data=omit.block.data, FUN=function(x) {
                                        c(mean=mean(x), sd=sd(x), n=length(x))})) 
call.test <- t.test(total_com_errors~block_num, data=omit.block.data, var.equal=TRUE,
                    paired=TRUE, subset=group=="Call")
sms.test <- t.test(total_com_errors~block_num, data=omit.block.data, var.equal=TRUE,
                   paired=TRUE, subset=group=="SMS")
noalert.test <- t.test(total_com_errors~block_num, data=omit.block.data, var.equal=TRUE,
                       paired=TRUE, subset=group=="No Alert")
call.test$statistic/sqrt(call.test$parameter["df"]+1) # Call
sms.test$statistic/sqrt(sms.test$parameter["df"]+1) # SMS
noalert.test$statistic/sqrt(noalert.test$parameter["df"]+1) # No Alert

# Compare groups on the number of fast responses made.
contrasts(trial.data$block_num) <- c(-1, 1)
alertsVsNoAlert <- c(1, -2, 1)
callVsText <- c(1, 0, -1)
contrasts(trial.data$group) <- cbind(alertsVsNoAlert, callVsText)
trial.data$fast_rt <- NA
trial.data$fast_rt[trial.data$resp_rt >= 0.1889059] <- 0
trial.data$fast_rt[trial.data$resp_rt < 0.1889059] <- 1
b.fit1 <- glmer(fast_rt~(1|part_num), family=binomial,
                data=trial.data)
b.fit2a <- glmer(fast_rt~block_num+(1|part_num), family=binomial,
                 data=trial.data)
b.fit2b <- glmer(fast_rt~group+(1|part_num), family=binomial,
                 data=trial.data)
b.fit3 <- glmer(fast_rt~block_num+group+(1|part_num), family=binomial,
                data=trial.data)
b.fit4 <- glmer(fast_rt~block_num*group+(1|part_num), family=binomial,
                data=trial.data)
anova(b.fit1, b.fit2a) # Block main effect.
anova(b.fit1, b.fit2b) # Group main effect.
anova(b.fit3, b.fit4) # Group by block interaction.
summary(b.fit4)
#confint(b.fit4) # 95% confidence intervals for the coefficients (takes awhile to execute).
hist(ranef(b.fit4)$part_num[, 1]) # Check if intercepts are normally distributed.

# Create a bar graph displaying the probability of making a fast response for
# each cell.
plot2.data <- expand.grid(block_num=c("1", "2"), group=unique(omit.data$group))
plot2.data$block_num <- factor(plot2.data$block_num)
prob.and.se2 <- predictSE(b.fit4, newdata=plot1.data, type="response")
plot2.data$prob = prob.and.se2$fit
plot2.data$prob_se = prob.and.se2$se.fit
plot2.data$group2 <- c("Call", "Call", "Text message", "Text message",
                       "No notification", "No notification")
plot2.data$block2 <- rep(c("Block 1", "Block 2"), 3)
plot2 <- ggplot(plot2.data, aes(x=group2, y=prob)) +
  geom_bar(aes(fill=block2), position=position_dodge(0.9), stat = "identity") +
  geom_errorbar(aes(ymin=prob-prob_se, ymax=prob+prob_se, group=block2), 
                position=position_dodge(0.9), stat="identity", width=0.5) +  
  ylab("Probability of making a fast response") +
  xlab("") +
  scale_y_continuous(labels=percent, expand=c(0, 0)) +
  scale_x_discrete(expand=c(0, 0)) +  
  scale_fill_manual("", values=c("#12006f", "#8f9a9d")) +  
  plot.theme
plot2
ggsave(filename="Figure3_FastRtProbPlot.tiff", dpi=600)

# Directly compare each notifcation group to the no alert control group 
# (commission errors).
no.call.omit.data <- omit.data[omit.data$group!="Call", ]
no.call.omit.data$group <- factor(no.call.omit.data$group)
no.sms.omit.data <- omit.data[omit.data$group!="SMS", ]
no.sms.omit.data$group <- factor(no.sms.omit.data$group)
callvscontrol.fit <- glmer(com_error~block_num*group+(1|part_num), 
                           family=binomial, data=no.sms.omit.data)
summary(callvscontrol.fit)
#confint(callvscontrol.fit) # 95% confidence intervals for the coefficients (takes awhile to execute).
hist(ranef(callvscontrol.fit)$part_num[, 1]) # Check if intercepts are normally distributed.
smsvscontrol.fit <- glmer(com_error~block_num*group+(1|part_num), 
                           family=binomial, data=no.call.omit.data)
summary(smsvscontrol.fit)
#confint(smsvscontrol.fit) # 95% confidence intervals for the coefficients (takes awhile to execute).
hist(ranef(smsvscontrol.fit)$part_num[, 1]) # Check if intercepts are normally distributed.

# Directly compare each notifcation group to the no alert control group 
# (fast responses).
no.call.trial.data <- trial.data[trial.data$group!="Call", ]
no.call.trial.data$group <- factor(no.call.trial.data$group)
no.sms.trial.data <- trial.data[trial.data$group!="SMS", ]
no.sms.trial.data$group <- factor(no.sms.trial.data$group)
callvscontrol.fit2 <- glmer(fast_rt~block_num*group+(1|part_num), family=binomial,
                      data=no.sms.trial.data)
summary(callvscontrol.fit2)
#confint(callvscontrol.fit2) # 95% confidence intervals for the coefficients (takes awhile to execute).
hist(ranef(callvscontrol.fit2)$part_num[, 1]) # Check if intercepts are normally distributed.
smsvscontrol.fit2 <- glmer(fast_rt~block_num*group+(1|part_num), family=binomial,
                            data=no.call.trial.data)
summary(smsvscontrol.fit2)
#confint(smsvscontrol.fit2) # 95% confidence intervals for the coefficients (takes awhile to execute).
hist(ranef(smsvscontrol.fit2)$part_num[, 1]) # Check if intercepts are normally distributed.