##------------------------------------------------------
## load libraries
##------------------------------------------------------
if (require("stringr") == FALSE) install.packages("stringr"); library(stringr)     # string functions
if (require("reshape2") == FALSE) install.packages("reshape2"); library(reshape2)  # convert data from wide to long format
if (require("ez") == FALSE) install.packages("ez"); library(ez)                    # RM-ANOVA
if (require("car") == FALSE) install.packages("car"); library(car)                 # analysis-of-variance tables for statistical models
if (require("lme4") == FALSE) install.packages("lme4"); library(lme4)              # generalized linear mixed-effects model (GLMM)
if (require("ARTool") == FALSE) install.packages("ARTool"); library(ARTool)        # Aligned Rank Transform (ART) for nonparametric factorial ANOVAs
if (require("ggplot2") == FALSE) install.packages("ggplot2"); library(ggplot2)     # graphics and plots

source("Rallfun-v38.txt") # Rand Wilcox's R library for robust statistics
                          # available from https://dornsife.usc.edu/labs/rwilcox/software/
                          # imported for Brunner-Puri within groups rank-based ANOVA

##------------------------------------------------------
## Analyzes participants' preferences for wearable devices.
## Preferences are expressed as ratings on 5-point Likert scales.
##
## Note: multiple tests are conducted (RM-ANOVA, ART, and BPRM); RM-ANOVA is
## not appropriate for this data because the dependent variable is
## ordinal, but we run it as a first model; BPRM is an improvement
## over Friedman's nonparametric test.
##------------------------------------------------------
analyzePreferenceRatings <- function(filename) 
{
  # load data
  data = read.csv(filename)
  cat(filename, ": N=", nrow(data), " records read\r\n", sep="")
  data$Participant = as.factor(data$Participant)
  str(data)
  
  cat("\r\n=======================================================")
  cat("\r\nRunning descriptive statistics:")
  cat("\r\n=======================================================\r\n")
  for(idxWearable in seq(2,6)) 
  {
    test = shapiro.test(data[,idxWearable])
    cat(str_pad(colnames(data)[idxWearable], width=15), 
        " | M = ", sprintf("%.1f", mean(data[,idxWearable])), 
        " | Mdn = ", sprintf("%.1f", median(data[,idxWearable])),
        " | M20 = ", sprintf("%.1f", mean(data[,idxWearable], trim=.2)),
        " | Shapiro-Wilk W = ", sprintf("%.3f", test$statistic), " p = ", sprintf("%.3f", test$p.value), ifelse(test$p.value <.01, " ** ", ifelse(test$p.value < .05, " * ", "  ")), "\r\n", sep ="")
  } 
  
  cat("\r\n=======================================================")
  cat("\r\nRunning RM-ANOVA:")
  cat("\r\n=======================================================")
  data.long = melt(data, id.vars=c("Participant"), variable.name="Wearable", value.name="Rating")
  data.long$Participant = as.factor(data.long$Participant)
  anova = ezANOVA(data.long, type = 3, dv = .(Rating), wid = .(Participant), within =.(Wearable), detailed = TRUE)
  print(anova)
  cat("Partial eta squared:")
  print(anova$ANOVA$SSn/(anova$ANOVA$SSn+anova$ANOVA$SSd))
  
  n = length(colnames(data)) - 1
  alpha = 0.05/(n*(n-1)/2)
  cat("\r\n=======================================================")
  cat("\r\nRunning pairwise t-tests with Bonferroni corrected alpha = ", alpha, ":", sep = "")
  cat("\r\n=======================================================\r\n")
  for(idxWearable1 in seq(2,5)) 
  {
    for(idxWearable2 in seq(idxWearable1+1,6)) 
    {
      test = t.test(data[,idxWearable1], data[,idxWearable2], paired = TRUE)
      cat(str_pad(colnames(data)[idxWearable1], width=15), str_pad(colnames(data)[idxWearable2], width=15), 
          " | t(", test$parameter, ") = ", sprintf("%6.3f", test$statistic),
          " p = ", sprintf("%.3f", test$p.value), ifelse(test$p.value < alpha, " *", " "),
          "\r\n", sep = "")
    }
  }
  
  cat("\r\n=======================================================")
  cat("\r\nRunning the ART test:")
  cat("\r\n=======================================================\r\n")
  data.long = melt(data, id.vars = c("Participant"), variable.name="Wearable", value.name="Rating")
  data.long$Participant = as.factor(data.long$Participant)
  model = art(Rating ~ Wearable + Error(Participant), data = data.long)
  print(summary(model))
  print(anova(model))
  print(art.con(model, "Wearable"))
  
  cat("\r\n=======================================================")
  cat("\r\nRunning robust statistics (BPRM):")
  cat("\r\n=======================================================\r\n")
  data.long = melt(data, id.vars = c("Participant"), variable.name="Wearable", value.name="Rating")
  data.long$Participant = as.factor(data.long$Participant)
  groups = fac2list(data.long$Rating, data.long$Wearable, pr = FALSE)
  test = bprm(groups)
  cat(sprintf("F_{(%.3f,\\infty)}=%.3f, p=%.3f", test$nu1, test$test.stat, test$p.value))
}

analyzePreferenceRatings("perceived-ease-of-access.csv")
analyzePreferenceRatings("perceived-ease-of-wear.csv")
analyzePreferenceRatings("perceived-ease-of-use.csv")
analyzePreferenceRatings("perceived-desirability.csv")
analyzePreferenceRatings("perceived-fun.csv")
analyzePreferenceRatings("perceived-unwanted-attention.csv")


##------------------------------------------------------
## Analyzes participants' preferences for suitable input 
## modalities for wearable devices. 
##------------------------------------------------------
analyzePreferencesForInputModalities <- function(filename)
{
  # load data
  data = read.csv(filename)
  cat(filename, ": N=", nrow(data), " records read\r\n", sep="")
  data$Participant = as.factor(data$Participant)
  data$Wearable = as.factor(data$Wearable)
  data$Context = as.factor(data$Context)
  str(data)
  
  cat("\r\n=======================================================")
  cat("\r\nRunning the ART test:")
  cat("\r\n=======================================================\r\n")
  model = art(Rating ~ Wearable * Context + Error(Participant), data = data)
  print(summary(model))
  print(anova(model))
  print(art.con(model, "Wearable"))
  print(art.con(model, "Context"))
}

analyzePreferencesForInputModalities("input-modalities.csv")


##------------------------------------------------------
## Analyzes participants' willingness to use wearable
## devices in public contexts.
##------------------------------------------------------
analyzeSocialAcceptability <- function(filename)
{
  # load data
  data = read.csv(filename)
  cat(filename, ": N=", nrow(data), " records read\r\n", sep="")
  data$Participant = as.factor(data$Participant)
  data$Wearable = as.factor(data$Wearable)
  data$Context = as.factor(data$Context)
  str(data)
  
  cat("\r\n=======================================================")
  cat("\r\nRunning GLM:")
  cat("\r\n=======================================================\r\n")
  model = glmer(formula = Rating ~ Wearable * Context + (1|Participant), data = data, family = binomial, control=glmerControl(optimizer="bobyqa"))
  print(Anova(model, type ="III"))
}

analyzeSocialAcceptability("acceptability-location.csv")
analyzeSocialAcceptability("acceptability-audience.csv")