---
title: "Human stress responses in office-like environments with wood furniture: Data processing and analysis"
author: | 
        Michael Burnard |
        InnoRenew CoE & University of Primorska |
date: "5 June 2019"
output:
  pdf_document:
    df_print: kable
    toc: true
    toc_depth: 3
documentclass: article
classoption: a4paper
geometry: margin=2cm 
---

\newpage   
# Introduction

This document contains the code and procedures for the analysis of salivary free cortisol levels used as an indicator of stress in "Human stress responses in office-like environments with wood furniture".

The analysis includes extracting cortisol levels and other information from the microtitre plates saliva samples were processed on and various statistical, summary, and visual analyses used to understand the data and determine the results of the experiment.

The data are not included because individuals and their health data may be identifiable. Sanitized, psuedoanonymous data may be available upon request in qualified circumstances.

## Acknowledgements

The authors gratefully acknowledge the European Commission for funding the InnoRenew CoE project (Grant Agreement #739574) under the Horizon2020 Widespread-Teaming program and the Slovenian infrastructure programme IO-0035.

# Setting the Environment

```{r setup, results="hide", message=FALSE}
knitr::opts_chunk$set(cache = TRUE)
library(drc)
library(tidyverse)
library(stringr)
library(ggforce)
library(scales)
library(zoo)
```



# Data

## Prepare storage data frames   

Here we:

  1. Create a data frames to store our important info   
    a. Results (readings, cortisol levels, participant id, test info)   
    b. Information about the model and its fit   
    c. Information about each processed tray (calibrators, controls, etc)   
  2. Create a list to store models. We want to come back to inspect them later.
```{r createDFs}
results <- tibble("Subject" = as.character(), "Test" = as.integer(),
                      "Interval" = as.numeric(), "Tray" = as.integer(),
                      "Reading.mu" = as.double(), "Cortisol.fitted.ng" = as.double(),
                      "Cortisol.fitted.nmol"=as.numeric())
control <- tibble("Name" = as.character(), "Reading.mu" = as.double(),
                      "Tray" = as.integer(), "Cortisol.fitted.ng" = as.double(),
                      "Cortisol.fitted.nmol"=as.numeric())
model.info <- tibble("Tray" = as.integer(), "b" = as.double(),
                         "c" = as.double(), "d" = as.double(),
                         "e" = as.double(), "Resid.var" = as.double())
tray.info <- tibble("Tray" = as.integer(), "Reading.name" = as.character(),
                        "Reading.value_1" = as.double(), 
                        "Reading.value_2" = as.double(), "Reading.mu" = as.double(), 
                        "Reading.absDiff"=as.double(),
                        "Fit.value.ng" = as.double(), "Fit.value.nmol"=as.numeric(),
                        "ExpectedRaw.value" = as.double(), 
                        "ExpectedFit.value" = as.numeric())
models.list <- vector("list", 25)
```

\newpage   

## Read, parse, and calculate cortisol concentrations   

Data files from the microtitre plate reader were stored in a single folder that was read, processed, and stored as follows. Files associating the readings with subject identifiers were also read and matched to the data for processing.

```{r warning=FALSE}
for(i in 1:25)
{
  #read data
  readings.path <- paste("data/cortisol/Readings/CortisolReadings_tray",
                         as.character(i), ".csv", sep="")
  layout.path <-  paste("data/cortisol/Layouts/CortisolLayout_tray",
                         as.character(i), ".csv", sep="")
 
  readings.all <- read.csv(readings.path, stringsAsFactors=FALSE)
  layout.all <- read.csv(layout.path, stringsAsFactors=FALSE)
  
  #prep data
  readings.all <- gather(readings.all, "Col", "Reading", 2:13)
  readings.all$Col <- as.numeric(str_sub(readings.all$Col, 2))
  layout.all <- gather(layout.all, "Col", "Name", 2:13)
  readings.all$Name <- layout.all$Name
  
  #Split data into useful groups
  readings.subject <- readings.all %>% filter(str_length(Name) > 5)
  readings.calib <- readings.all %>% filter(str_detect(Name, "C{1}[0123456]") == TRUE)
  readings.control_a <- readings.all %>% filter(str_detect(Name, "C{1}[LHM]") == TRUE)
  readings.control_b <- readings.all %>% filter(str_detect(Name, "C{1}[o]") == TRUE)

  #two temporary storage frames
  readings.calib_t <- readings.calib[,3:4]
  readings.control_t <- bind_rows(readings.control_a[,3:4], readings.control_b[,3:4])

  #finalise data prep
  readings.subject <- readings.subject %>% 
    separate(Name, c("Subject", "Test", "Interval", "Position"), convert=TRUE) %>% 
    group_by(Subject, Test, Interval) %>% summarise(Reading.mu = mean(Reading)) %>%
    mutate(Tray = i)
  readings.calib <- readings.calib %>% separate(Name, c("Name", "Position")) %>%
    group_by(Name) %>% summarise(Reading.mu = mean(Reading)) %>% mutate(Tray = i)
  readings.control_a <- readings.control_a %>% separate(Name, c("Name", "Position")) %>%
     group_by(Name) %>% summarise(Reading.mu = mean(Reading)) %>% mutate(Tray = i)
  readings.control_b <- readings.control_b %>% ungroup() %>% 
    mutate(Reading.mu = Reading, Tray = i)
  readings.control <- bind_rows(readings.control_a, readings.control_b[,c(4,5,6)])
  #the same except for summarise for the temp storage frames
  readings.calib_t <- readings.calib_t %>% separate(Name, c("Name", "Position"))
  readings.control_t <- readings.control_t %>% separate(Name, c("Name", "Position"))

  #add cortisol values for model fitting & storage
  readings.calib$Cortisol <- c(0.0, 0.5, 1, 5, 10, 20, 100)
  
  #fit a model and store it in our list of models to look at later.
  #there is a special case for tray 23 because of a robot error 
  #the robot missed a row (probably a pipette wasn't well attached)
  if(i == 23) {
    readings.calib <- readings.calib %>% filter(Name != "C1")
    m <- drm(Cortisol~Reading.mu, data=readings.calib, fct=LL.4())
    models.list[[i]] <- m
    names(models.list)[[i]] <- paste("Tray", i, sep="_")
    tmp.row <- data_frame("Name"="C1", "Reading.mu"=NA, "Tray"=i, "Cortisol"=0.5)
    readings.calib <- bind_rows(readings.calib, tmp.row) %>% arrange(Cortisol)
  } else {
    m <- drm(Cortisol~Reading.mu, data=readings.calib, fct=LL.4())
    models.list[[i]] <- m
    names(models.list)[[i]] <- paste("Tray", i, sep="_")
  }
  
  #store model info in a data frame to compare
  temp.row <- data.frame("Tray" = i, "b" = coef(m)[[1]], "c" = coef(m)[[2]],
                         "d" = coef(m)[[3]], "e" = coef(m)[[4]], 
                         "Resid.var" = summary(m)[[1]])
  model.info <- bind_rows(model.info, temp.row)
  
  #store calibration and control readings
  readings.calib_t <- spread(readings.calib_t, Position, Reading)
  temp.df <- as.data.frame(cbind(rep(i, 7), 
                         readings.calib_t$Name, 
                         readings.calib_t[2],
                         readings.calib_t[3], 
                         (readings.calib_t[2] + readings.calib_t[3]) / 2, 
                         abs(readings.calib_t[2] - readings.calib_t[3]),
                         rep(1,7), 
                         c(2.643, 2.229, 2.056, 1.295, 0.789, 0.451, 0.131), 
                         readings.calib$Cortisol))
  names(temp.df) <- c("Tray", "Reading.name", "Reading.value_1" ,
                      "Reading.value_2", "Reading.mu",
                      "Reading.absDiff", "Fit.value.ng", 
                      "ExpectedRaw.value", "ExpectedFit.value")
  temp.df$Fit.value.ng <- fitted(m, temp.df[,c(4,1)])
  temp.df$Fit.value.nmol <- temp.df$Fit.value.ng * 2.76 # add nmol/L value
  
  #another correction for tray 23
  if(i == 23) { temp.df <- temp.df[-2,]}
  tray.info <- bind_rows(tray.info, temp.df)
  tray.info$ExpectedFit.value.nmol <- tray.info$ExpectedFit.value * 2.76

  #use the fitted model to back calculate
  readings.subject$Cortisol.fitted.ng <- fitted(m, as.data.frame(readings.subject[,4]))
  readings.subject$Cortisol.fitted.nmol <- readings.subject$Cortisol.fitted.ng * 2.76
  results <- bind_rows(results, readings.subject[,c(1,2,3,5,7,4,6)])
  
  #and for controls
  readings.control$Cortisol.fitted.ng <- fitted(m, as.data.frame(readings.control[,2]))
  readings.control$Cortisol.fitted.nmol <- readings.control$Cortisol.fitted.ng * 2.76
  control <- bind_rows(control, readings.control)
  
}
```

Now, merge Test Record data (video, room, test sequence, etc.) with results. We also make an pseudoanonymous identifier that isn't associable with the subject identities. Because of anonymisation, this isn't run here. Instead we read a file with the saved results.
```{r, eval=FALSE}
#Read the file and select the columns of interest
tr <- read.csv("data/TestRecord.csv", stringsAsFactors=FALSE)
tr <- tr %>% select(Participant, TestType, Room, Video, SGL, WHO5_Total)
#rename for the merge
names(tr)[1] <- "Subject"
names(tr)[5] <- "Test"

#Merge them with the results
results <- left_join(results, tr, by=c("Subject", "Test"))

#Add a character string to further anonomise the results in output plots.
makeRandString <- function() {
  tmp = c(sample(LETTERS, 3, replace=TRUE), 
    sample(0:9, 2, replace=TRUE))
  return(paste0(tmp, collapse=""))
}


results <- results %>% group_by(Subject) %>% mutate(Anon=makeRandString())
# verify unique for each subject
length(unique(results$Anon))

```


```{r}

results <- read.csv("results_storage.csv", stringsAsFactors = FALSE, header = TRUE)
#to double check we have 14 total for each subject

verify <- data.frame(table(results$Subject, results$Tray))
verify <- spread(verify, Var2, Freq)
verify <- verify %>% mutate(RowSum = rowSums(verify[,2:26]))

#The pipette robot had a small malfuction (didn't seal a pipette tip) on one tray
#rendering the values from row unuseable.
results$Cortisol.fitted.ng <- ifelse(results$Tray == 23 & results$Interval == 2, NA,
                                     results$Cortisol.fitted.ng)
results$Cortisol.fitted.nmol <- ifelse(results$Tray == 23 & results$Interval == 2, NA,
                                       results$Cortisol.fitted.nmol)

#add time in minutes
results$Minutes <- ifelse(results$Interval == 1, 0, 
                    ifelse(results$Interval == 2, 15,
                    ifelse(results$Interval == 3, 25,
                    ifelse(results$Interval == 4, 35,
                    ifelse(results$Interval == 5, 45,
                    ifelse(results$Interval == 6, 60, 75))))))

#3 subject-test combos were duplicated on tray 25 for verification
# XX-2, XX-1, XX-2. Need to remove these from the general results 
#and keep them to check.
check <- results %>% filter ((Anon %in% c("XGD55","FQX13","KHO49") & Tray == 25))
results <- results %>% filter( !(Anon %in% c("XGD55","FQX13","KHO49") & Tray == 25) )



```


#Plate to Plate comparisons and results
```{r, warning=FALSE}
ggplot(data=tray.info, 
       aes(x=ExpectedFit.value * 2.76, y=Fit.value.ng * 2.76)) + theme_bw() +
  geom_point(alpha=.6, shape=1) +
  scale_x_continuous(limits=c(0,300), breaks=seq(0,300,25)) +
  scale_y_continuous(limits=c(0,300), breaks=seq(0,300,25)) +
  labs(x="Expected cortisol concentration (nmol/L)",
       y="Predicted cortisol concentration (nmol/L)")

tray.fit <- tray.info %>% select(Tray, Fit.value.nmol, ExpectedFit.value.nmol) %>% 
  group_by(Tray) %>%
  summarise(RMSE = sqrt(mean((Fit.value.nmol - ExpectedFit.value.nmol)^2)))

model.info$RMSE <- tray.fit$RMSE
#write.csv(file="Model_info.csv", model.info)
#write.csv(file="Calibrators.csv", tray.info %>% group_by(Reading.name)
#%>% summarise(mu = mean(Fit.value.nmol), med=median(Fit.value.nmol), 
#sd=sd(Fit.value.nmol), coef.v=sd/mu, min=min(Fit.value.nmol),max=max(Fit.value.nmol)) )

#these are articially high, because of what we are measuring
#the results are near perfect because we are measuring curves fit to a known value
#to other curves fit to the same known value.
r2 <- 1 - (sum((tray.info$ExpectedFit.value.nmol - tray.info$Fit.value.nmol)^2) / 
             sum((tray.info$Fit.value.nmol - mean(tray.info$Fit.value.nmol))^2))

rho <- cov(tray.info$Fit.value.nmol, tray.info$ExpectedFit.value.nmol) / 
  (sd(tray.info$Fit.value.nmol) * sd(tray.info$ExpectedFit.value.nmol))

tray.info %>% group_by(Reading.name) %>% summarise(min=min(Reading.mu))
```


## Duplicate testing   
Three sample groups were tested in duplicate (on different trays) to guage tray-to-tray variation.
```{r}
#plate to plate

#first, the duplicated trio

chck <- check %>% select(Subject, Anon, Test, Tray, Minutes, Reading.mu,
                         Cortisol.fitted.nmol)
r.chck <- results %>% select(Subject, Anon, Test, Tray, Minutes, Reading.mu, 
                             Cortisol.fitted.nmol) %>%
  filter( (Anon == "XGD55" & Test == 1) |
          (Anon == "FQX13" & Test == 2) |
          (Anon == "KHO49" & Test == 2))
r.chck$Type <- "Original"
chck$Type <- "Verify"
chck <- bind_rows(chck, r.chck)

ggplot(data=chck, aes(x=Minutes, y=Cortisol.fitted.nmol, colour=Type)) + theme_bw() +
  geom_point(shape=1) +
  facet_wrap(~Anon, nrow=3) +
  scale_x_continuous(limits=c(0,80), breaks=c(0,15,25,35,45,60,75)) +
  scale_colour_manual(values=c("#FF0000", "#0000FF")) +
  labs(x="Minutes", y="Mean Cortisol Level (nmol / L)")

#now values
chck.o <- chck %>% filter(Type == "Original") %>% 
  select(Subject, Cortisol.fitted.nmol)%>%  arrange(Subject) 
chck.v <- chck %>% filter(Type == "Verify") %>% 
  select(Subject, Cortisol.fitted.nmol) %>%  arrange(Subject) 
names(chck.o)[2] <- "Original"
names(chck.v)[2] <- "Verify"

chck.o$Verify <- chck.v$Verify

chck.o$RawDiff <- abs(chck.o$Verify - chck.o$Original)
chck.o$PerDiff <- abs((chck.o$Verify - chck.o$Original) / chck.o$Verify)

chck.o <- chck.o %>% group_by(Subject) %>% 
  summarise(raw.mu = mean(RawDiff), per.mu = mean(PerDiff))
```

## Model and fit parameters   
```{r}

#consistency of fit parameters
mod.con <- gather(model.info, "Parameter", "Value", 2:6)
mod.con$Parameter <- ifelse(mod.con$Parameter == "Resid.var",
                           "Residual Variance", mod.con$Parameter)

ggplot(data=mod.con[which(mod.con$Parameter != "Residual Variance"),], 
       aes(x=Parameter, y=Value)) + theme_bw() +
  geom_boxplot() +
  facet_wrap(~Parameter, scales="free") +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())

ggplot(data=mod.con[which(mod.con$Parameter == "Residual Variance"),], 
       aes(x=Parameter, y=Value)) + theme_bw() +
  geom_boxplot() +
  facet_wrap(~Parameter) +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())

#plot(models.list$Tray_23)

```

## Tray Controls   


```{r}
#Clear names for labels
control$Labels <- ifelse(control$Name == "CH", "External High",
                  ifelse(control$Name == "CL", "External Low",
                  ifelse(control$Name == "CM", "External Medium",
                  ifelse(control$Name == "Con-B", "Internal High",
                  ifelse(control$Name == "Con-A", "Internal Low", "Zed")))))
control$Labels <- as.factor(control$Labels)
control$Labels <- factor(control$Labels, levels=c("External Low", "External Medium",
                                                  "External High", "Internal Low",
                                                  "Internal High"))

ggplot(data=control %>% filter(Tray != 23), aes(x=Name, y=Cortisol.fitted.nmol)) + 
  theme_bw() +
  geom_boxplot() +
  facet_wrap(~Labels, scales="free") +
  theme(axis.text.x=element_blank(), axis.ticks.x=element_blank())
```

# Results inspection for comparisons   


Now we can graphically and statistically inspect the results.
```{r}

#change page=X to see different pages.
ggplot(data=results, aes(x=Interval, y=Cortisol.fitted.nmol, shape=TestType)) + 
  theme_bw() +
  geom_point() +
  facet_wrap_paginate(~Anon, ncol=2, nrow=5, page=1)

#generally all are fine except for DYG49. It's a very bizarre pattern for the
#treated test. Strange enough to warrant removal.
#The pattern doesn't follow expectations, and the values for tests are too
#different to provide much confidence. One situation or the other might be acceptable
#both are not.

r1 <- results %>% filter(Anon != "DYG49")

#first, check the responses to videos were not different


r1.vid <- r1 %>% group_by(Anon, Video) %>% filter(Interval %in% c(5,6,7)) %>%
  summarise(Test.mu=mean(Cortisol.fitted.nmol, na.rm=TRUE))
r1.vid <- r1.vid %>% group_by(Anon) %>% mutate(testDif=Test.mu - lag(Test.mu))
r1.vid$Outcome <- ifelse(is.na(r1.vid$testDif), "X", 
                 ifelse(r1.vid$testDif < 0, "Positive",
                  "Negative"))
r1.vid$Outcome <- ifelse(r1.vid$Outcome == "X", lead(r1.vid$Outcome), r1.vid$Outcome)
r1.vid <- r1.vid %>% select(Anon, Video, Test.mu) %>%
  ungroup() %>% spread(Video, Test.mu)

wilcox.test(r1.vid$A, r1.vid$B, paired=TRUE)
```

```{r}
r1.int <- r1 %>% select(Subject, Anon, Room, TestType, Interval,
                        Minutes, Cortisol.fitted.nmol) %>% 
  spread("TestType", Cortisol.fitted.nmol) %>% 
  mutate(Labels = ifelse(Room == "Oak", "Room A", "Room B"))

r1.int <- r1.int[complete.cases(r1.int),]

r1.int$Response <- r1.int$Control - r1.int$Treated
r1.int$Response.per <- r1.int$Response / r1.int$Control
r1.int$Outcome <- ifelse(r1.int$Response > 0, "Positive", "Negative")
r1.int.mu <- r1.int %>% group_by(Labels, Minutes) %>% 
  summarise(mu.raw = mean(Response),mu.per = mean(Response.per))
r1.int.mu$Outcome <- ifelse(r1.int.mu$mu.raw > 0, "Positive", "Negative")


ggplot(data=r1.int, aes(x=Minutes, y=Response, shape=Outcome)) + theme_bw() +
  geom_point(colour="#999999", alpha=0.5, size=5) +
  geom_point(data=r1.int.mu, aes(x=Minutes, y=mu.raw), shape=5, size=2.5) +
  geom_hline(yintercept=0, colour="#000000") +
  facet_wrap(~Labels, nrow=2) +
  scale_shape_manual(values=c(45,43)) +
  scale_x_continuous(limits=c(0,80), breaks=c(0,15,25,35,45,60,75)) +
  labs(y="Cortisol concentration difference in nmol/L (control - treated)")

#percent basis? Its difficult to interpret but offers another perspective
ggplot(data=r1.int, aes(x=Minutes, y=Response.per, group=Anon)) + theme_bw() +
  geom_line(colour="#BBBBBB", linetype=3) +
  geom_point(data=r1.int.mu, aes(x=Minutes, y=mu.per, group=1)) +
  geom_hline(yintercept=0, colour="#333333", size=0.2) +
  facet_wrap(~Labels, nrow=2) +
  scale_x_continuous(limits=c(0,80), breaks=c(0,15,25,35,45,60,75)) +
  scale_y_continuous(labels=percent_format())
```


## Mean overall stress response    

This is parametersied as the mean cortisol level in the control room minus the mean cortisol level in the treated room.
```{r}
#overall mean per subject

####
# Check mean of interval-interval differences (should be the same)
####

r1.mu <- r1 %>% group_by(Subject, Anon, TestType, Room) %>% 
  summarise(Test.mu=mean(Cortisol.fitted.nmol, na.rm=TRUE))

r1.mu <- r1.mu %>% group_by(Subject, Anon) %>% mutate(testDif=Test.mu - lag(Test.mu))
r1.mu$Outcome <- ifelse(is.na(r1.mu$testDif), "X", 
                 ifelse(r1.mu$testDif < 0, "Positive",
                  "Negative"))
r1.mu$Outcome <- ifelse(r1.mu$Outcome == "X", lead(r1.mu$Outcome), r1.mu$Outcome)

r1.mu$Label <- ifelse(r1.mu$Room == "A-Walnut", "Walnut", "Oak") 

ggplot(data=r1.mu, aes(x=TestType, y=Test.mu, colour=Outcome, group=Subject)) +
  theme_bw() +
  geom_line(alpha=.3) +
  geom_point(alpha=.7, shape=1) +
  facet_wrap(~Room) +
  scale_colour_manual(values=c("#FF0000", "#0000FF")) +
  labs(x="", y="Mean Cortisol Level (nmol / L)")

#create seperate data frames for each room
r1.mu.oak <- r1.mu %>% filter(Room == "Oak") %>% 
  select(Subject, Anon, Room, TestType, Test.mu) %>%
  ungroup() %>% spread(TestType, Test.mu)

r1.mu.wal <- r1.mu %>% filter(Room == "A-Walnut") %>% 
  select(Subject, Anon, Room, TestType, Test.mu) %>%
  ungroup() %>% spread(TestType, Test.mu)

#hist(r1.mu.wal$Control)
#hist(r1.mu.wal$Treated)
#hist(r1.mu.oak$Control)
#hist(r1.mu.oak$Treated)

#the data aren't normal enough, even with a log transform
#We use Wilcoxon rank sum, we're within-subjects so it's a paired test
#and with our parameterisation our alternative in this case is greater.
wilcox.test(r1.mu.oak$Control, r1.mu.oak$Treated, paired=TRUE,
            alternative="greater", conf.int=TRUE, conf.level=.95)
wilcox.test(r1.mu.wal$Control, r1.mu.wal$Treated, paired=TRUE, 
            alternative="greater", conf.int=TRUE, conf.level=.95)

#Summary info
r1.mu %>% filter(!is.na(testDif)) %>% group_by(Room, Outcome) %>% 
  mutate(realtd=testDif*-1) %>% summarise(mu=mean(realtd), sd=sd(realtd),
                                          min=min(realtd), max=max(realtd), n=n())

r1.mu.per <- bind_rows(r1.mu.oak, r1.mu.wal)
r1.mu.per$Response <- r1.mu.per$Control - r1.mu.per$Treated
r1.mu.per$Response.per <- (r1.mu.per$Control - r1.mu.per$Treated) / r1.mu.per$Control
r1.mu.per$Outcome <- ifelse(r1.mu.per$Response > 0, "Positive", "Negative")
r1.mu.per %>% group_by(Room) %>% summarise(mu=mean(Response.per),
                                           sd=sd(Response.per), min=min(Response.per), 
                                           max=max(Response.per), n=n())
```

## Response period stress level

```{r}
#mean per subject for periods intervals 4,5,6,7

r1.re <- r1 %>% group_by(Anon, TestType, Room) %>% filter(Interval %in% c(4,5,6,7)) %>%
  summarise(Test.mu=mean(Cortisol.fitted.nmol, na.rm=TRUE))

r1.re <- r1.re %>% group_by(Anon) %>% mutate(testDif=Test.mu - lag(Test.mu))
r1.re$Outcome <- ifelse(is.na(r1.re$testDif), "X", 
                 ifelse(r1.re$testDif < 0, "Positive",
                  "Negative"))
r1.re$Outcome <- ifelse(r1.re$Outcome == "X", lead(r1.re$Outcome), r1.re$Outcome)
r1.re$Label <- ifelse(r1.re$Room == "A-Walnut", "Walnut", "Oak") 
ggplot(data=r1.re, aes(x=TestType, y=Test.mu, colour=Outcome, group=Anon)) + 
  theme_bw() +
  geom_line(alpha=.3) +
  geom_point(alpha=.7, shape=1) +
  facet_wrap(~Label) +
  scale_colour_manual(values=c("#FF0000", "#0000FF")) +
  labs(x="", y="Mean Cortisol Level (nmol / L)")

#create seperate data frames for each room
r1.re.oak <- r1.re %>% filter(Room == "Oak") %>% 
  select(Anon, Room, TestType, Test.mu) %>%
  ungroup() %>% spread(TestType, Test.mu)

r1.re.wal <- r1.re %>% filter(Room == "A-Walnut") %>% 
  select(Anon, Room, TestType, Test.mu) %>%
  ungroup() %>% spread(TestType, Test.mu)

#hist(r1.re.wal$Control)
#hist(r1.re.wal$Treated)
#hist(r1.re.oak$Control)
#hist(r1.re.oak$Treated)
#the data aren't normal enough, even with a log transform
#We use Wilcoxon rank sum, we're within-subjects so it's a paired test
#and with our parameterisation our alternative in this case is greater.
wilcox.test(r1.re.oak$Control, r1.re.oak$Treated, paired=TRUE,
            alternative="greater", conf.int=TRUE, conf.level=.95)
wilcox.test(r1.re.wal$Control, r1.re.wal$Treated, paired=TRUE, 
            alternative="greater", conf.int=TRUE, conf.level=.95)

r1.re %>% filter(!is.na(testDif)) %>% group_by(Room, Outcome) %>% 
  mutate(realtd=testDif*-1) %>% 
  summarise(mu=mean(realtd), sd=sd(realtd), min=min(realtd), max=max(realtd), n=n())

r1.re.per <- bind_rows(r1.re.oak, r1.re.wal)
r1.re.per$Response <- r1.re.per$Control - r1.re.per$Treated
r1.re.per$Response.per <- (r1.re.per$Control - r1.re.per$Treated) / r1.re.per$Control
r1.re.per$Outcome <- ifelse(r1.re.per$Response > 0, "Positive", "Negative")
r1.re.per %>% filter(abs(Response.per) > 1) %>% group_by(Room, Outcome) %>% 
  summarise(mu=mean(Response.per), sd=sd(Response.per),
            min=min(Response.per), max=max(Response.per), n=n())
```
### Check for differences within-subjects in initial 3 samples   

If these are different this confuses the outcome slightly.
```{r}
#mean per subject for periods intervals 123

r1.re <- r1 %>% group_by(Anon, TestType, Room) %>% filter(Interval %in% c(1,2,3)) %>%
  summarise(Test.mu=mean(Cortisol.fitted.nmol, na.rm=TRUE))

r1.re <- r1.re %>% group_by(Anon) %>% mutate(testDif=Test.mu - lag(Test.mu))
r1.re$Outcome <- ifelse(is.na(r1.re$testDif), "X", 
                 ifelse(r1.re$testDif < 0, "Positive",
                  "Negative"))
r1.re$Outcome <- ifelse(r1.re$Outcome == "X", lead(r1.re$Outcome), r1.re$Outcome)
r1.re$Label <- ifelse(r1.re$Room == "A-Walnut", "Walnut", "Oak") 
ggplot(data=r1.re, aes(x=TestType, y=Test.mu, colour=Outcome, group=Anon)) +
  theme_bw() +
  geom_line(alpha=.3) +
  geom_point(alpha=.7, shape=1) +
  facet_wrap(~Label) +
  scale_colour_manual(values=c("#FF0000", "#0000FF")) +
  labs(x="", y="Mean Cortisol Level (nmol / L)")

#create seperate data frames for each room
r1.re.oak <- r1.re %>% filter(Room == "Oak") %>% 
  select(Anon, Room, TestType, Test.mu) %>%
  ungroup() %>% spread(TestType, Test.mu)

r1.re.wal <- r1.re %>% filter(Room == "A-Walnut") %>%
  select(Anon, Room, TestType, Test.mu) %>%
  ungroup() %>% spread(TestType, Test.mu)

#hist(r1.re.wal$Control)
#hist(r1.re.wal$Treated)
#hist(r1.re.oak$Control)
#hist(r1.re.oak$Treated)
#the data aren't normal enough, even with a log transform
#We use Wilcoxon rank sum, we're within-subjects so it's a paired test
#and with our parameterisation our alternative two-sided since this is a post-hoc test.
wilcox.test(r1.re.oak$Control, r1.re.oak$Treated, paired=TRUE,
            alternative="two.sided", conf.int=TRUE, conf.level=.95)
wilcox.test(r1.re.wal$Control, r1.re.wal$Treated, paired=TRUE, 
            alternative="two.sided", conf.int=TRUE, conf.level=.95)


r1.re %>% filter(!is.na(testDif)) %>% group_by(Room, Outcome) %>% 
  mutate(realtd=testDif*-1) %>% summarise(mu=mean(realtd), sd=sd(realtd), 
                                          min=min(realtd), max=max(realtd), n=n())

r1.re.per <- bind_rows(r1.re.oak, r1.re.wal)
r1.re.per$Response <- r1.re.per$Control - r1.re.per$Treated
r1.re.per$Response.per <- (r1.re.per$Control - r1.re.per$Treated) / r1.re.per$Control
r1.re.per$Outcome <- ifelse(r1.re.per$Response > 0, "Positive", "Negative")
r1.re.per %>% filter(abs(Response.per) > 1) %>% group_by(Room, Outcome) %>% 
  summarise(mu=mean(Response.per), sd=sd(Response.per), min=min(Response.per), 
            max=max(Response.per), n=n())
```
\newpage
## Response magnitude   

This is tricky. We find the minimum of either stage 4 or 5, then the maximum after 4 (could be the same as the minimum indicating no stress response detected).
```{r}
r1.min <-  r1 %>% group_by(Anon,TestType, Room) %>% filter(Interval %in% c(4,5)) %>%
  filter(Cortisol.fitted.nmol == min(Cortisol.fitted.nmol)) %>% arrange(Anon) %>% 
  mutate(Value="Min")
r1.max <- r1 %>% group_by(Anon,TestType, Room) %>% filter(Interval > 4) %>%
  filter(Cortisol.fitted.nmol == max(Cortisol.fitted.nmol)) %>% arrange(Anon) %>% 
  mutate(Value="Max")

r1.min.oak <- r1.min %>% filter(Room == "Oak") %>% 
  select(Anon, TestType, Room, Cortisol.fitted.nmol)
names(r1.min.oak)[4] <- "Min.cortisol.nmol"

r1.max.oak <- r1.max %>% filter(Room == "Oak") %>% 
  select(Anon, TestType, Room, Cortisol.fitted.nmol)

r1.mm.oak <- r1.min.oak
r1.mm.oak$Max.cortisol.nmol <- r1.max.oak$Cortisol.fitted.nmol
#r1.mm.oak$Interval.max <- r1.max.oak$Interval

r1.mm.oak$Response <- r1.mm.oak$Max.cortisol.nmol - r1.mm.oak$Min.cortisol.nmol
r1.mm.oak$Response.per <- (r1.mm.oak$Max.cortisol.nmol - r1.mm.oak$Min.cortisol.nmol) / 
  r1.mm.oak$Min.cortisol.nmol

r1.mm.oak.tst <- r1.mm.oak %>% ungroup() %>% select(Anon, TestType, Response.per, Room) %>% 
  spread(TestType, Response.per)

r1.min.wal <- r1.min %>% filter(Room == "A-Walnut") %>% 
  select(Anon, TestType, Room, Cortisol.fitted.nmol)
names(r1.min.wal)[4] <- "Min.cortisol.nmol"


r1.max.wal <- r1.max %>% filter(Room == "A-Walnut") %>% 
  select(Anon, TestType, Room, Cortisol.fitted.nmol)

r1.mm.wal <- r1.min.wal
r1.mm.wal$Max.cortisol.nmol <- r1.max.wal$Cortisol.fitted.nmol


r1.mm.wal$Response <- r1.mm.wal$Max.cortisol.nmol - r1.mm.wal$Min.cortisol.nmol
r1.mm.wal$Response.per <- (r1.mm.wal$Max.cortisol.nmol - r1.mm.wal$Min.cortisol.nmol) / 
  r1.mm.wal$Min.cortisol.nmol

r1.mm.wal.tst <- r1.mm.wal %>% ungroup() %>% select(Anon, TestType, Response.per, Room) %>% 
  spread(TestType, Response.per)

#visualisation
r1.mm <- bind_rows(r1.mm.oak.tst, r1.mm.wal.tst)
r1.mm$diff <- r1.mm$Control - r1.mm$Treated
r1.mm$Outcome <- ifelse(r1.mm$diff == 0, "Neutral",
                        ifelse(r1.mm$diff > 0, "Positive", "Negative"))
r1.mm <- r1.mm %>% gather("TestType", "Response", 3:4)
r1.mm$Labels <- ifelse(r1.mm$Room == "A-Walnut", "Walnut", "Oak")


ggplot(data=r1.mm, aes(x=TestType, y=Response, colour=Outcome, group=Anon)) + 
  theme_bw() +
  geom_line(alpha=.3) +
  geom_point(alpha=.7, shape=1) +
  facet_wrap(~Labels) +
  scale_y_continuous() +
  scale_colour_manual(values=c("#FF0000", "#000000", "#0000FF")) +
  labs(x="", y="Cortisol concentrtation (nmol / L )")


#the data aren't normal enough, even with a log transform
#We use Wilcoxon rank sum, we're within-subjects so it's a paired test
#and with our parameterisation our alternative in this case is greater.
wilcox.test(r1.mm.oak.tst$Control, r1.mm.oak.tst$Treated, paired=TRUE,
            alternative="greater", conf.int=TRUE, conf.level=.95)
wilcox.test(r1.mm.wal.tst$Control, r1.mm.wal.tst$Treated, paired=TRUE,
            alternative="greater", conf.int=TRUE, conf.level=.95)


#Percent base
r1.mm %>% filter(TestType == "Treated") %>% group_by(Room, Outcome) %>% 
  summarise(mu=mean(diff), sd=sd(diff), min=min(diff), 
            max=max(diff), n=n())

r1.mm %>% group_by(Room) %>% filter(TestType == "Treated") %>%
  summarise(mu=mean(diff), sd=sd(diff),
            min=min(diff), max=max(diff), n=n())

```
\newpage   

## Recovery   

To parameterise recovery, we take the max value of intervals 4, 5, 6, and 7, and the final value (interval 7). If those two are equal (i.e., the max value is at interval 7), that means recovery didn't take place in the test period. We compare as a percent of the max value returned, rather than raw differences.
```{r}
r1.rec <- r1 %>% group_by(Anon,TestType, Room) %>% filter(Interval > 3) %>%
  filter(Cortisol.fitted.nmol == max(Cortisol.fitted.nmol)) %>% 
  arrange(Anon, TestType) 

r1.rec7 <- r1 %>% group_by(Anon,TestType, Room) %>% filter(Interval == 7) %>% 
  arrange(Anon, TestType)

r1.rec <- r1.rec %>% select(Anon, Interval, TestType, Room, Cortisol.fitted.nmol) 
r1.rec$Final <- r1.rec7$Cortisol.fitted.nmol
names(r1.rec)[5] <- "Max"

r1.rec$RawDiff <- r1.rec$Max - r1.rec$Final
r1.rec$PerDiff <- (r1.rec$Max - r1.rec$Final) / r1.rec$Max

r1.rec.per.o <- r1.rec %>% filter(Room == "Oak") %>% 
  select(Anon, TestType, PerDiff, Room) %>% 
  spread(TestType, PerDiff)
r1.rec.per.w <- r1.rec %>% filter(Room == "A-Walnut") %>%
  select(Anon, TestType, PerDiff, Room) %>% 
  spread(TestType, PerDiff)

r1.rec.raw.o <- r1.rec %>% filter(Room == "Oak") %>% 
  select(Anon, TestType, RawDiff, Room) %>% 
  spread(TestType, RawDiff)
r1.rec.raw.w <- r1.rec %>% filter(Room == "A-Walnut") %>% 
  select(Anon, TestType, RawDiff, Room) %>% 
  spread(TestType, RawDiff)

#visualisation
r1.recv <- bind_rows(r1.rec.per.o, r1.rec.per.w)
r1.recv$diff <- r1.recv$Control - r1.recv$Treated
r1.recv$Outcome <- ifelse(r1.recv$diff == 0, "Neutral",
                        ifelse(r1.recv$diff > 0, "Positive", "Negative"))
r1.recv <- r1.recv %>% gather("TestType", "Response", 3:4)
r1.recv$Labels <- ifelse(r1.recv$Room == "A-Walnut", "Walnut", "Oak")

ggplot(data=r1.recv, aes(x=TestType, y=Response, colour=Outcome, group=Anon)) + 
  theme_bw() +
  geom_line(alpha=.3) +
  geom_point(alpha=.7, shape=1) +
  facet_wrap(~Labels) +
  scale_y_continuous() +
  scale_colour_manual(values=c("#FF0000", "#0000FF")) +
  labs(x="", y="Cortisol concentrtation (nmol / L )")

#parameterisation changes here:
#we're looking at recovery, and our hypothesis is that one recovers less
#in the control room. We keep greater, and reverse the order we specify control and
#treated
wilcox.test(r1.rec.per.o$Treated, r1.rec.per.o$Control, paired=TRUE,
            alternative="greater", conf.int = TRUE, conf.level = .95, correct=FALSE)
wilcox.test(r1.rec.per.w$Treated, r1.rec.per.w$Control, paired=TRUE,
            alternative="greater", conf.int = TRUE, conf.level = .95, correct=FALSE)

r1.rec.per <- bind_rows(r1.rec.per.w, r1.rec.per.o)
r1.rec.per$Resp.diff <- r1.rec.per$Treated - r1.rec.per$Control
r1.rec.per$Outcome <- ifelse(r1.rec.per$Resp.diff > 0, "Positive", "Negative")

r1.rec.per %>%  group_by(Room, Outcome) %>% 
  summarise(mu=mean(Resp.diff), sd=sd(Resp.diff), min=min(Resp.diff), 
            max=max(Resp.diff), n=n())
r1.rec.per %>%  group_by(Room) %>% 
  summarise(mu=mean(Resp.diff), sd=sd(Resp.diff), min=min(Resp.diff), 
            max=max(Resp.diff), n=n())

```


#WHO-5 Well-being index

```{r}
r.who.tst <- r1 %>% select(Anon, Test, TestType, Room, WHO5_Total) %>% 
  group_by(Anon, TestType, Room) %>% summarise(WHO5_Total = mean(WHO5_Total*4)) %>% 
  spread(TestType, WHO5_Total) %>% mutate(Diff = (Treated) - (Control))

r.who.tst$Outcome <- ifelse(r.who.tst$Diff == 0, "Neutral", 
                 ifelse(r.who.tst$Diff > 0, "Positive", "Negative"))

hist(r.who.tst$Diff)
t.test(r.who.tst$Control, r.who.tst$Treated, paired=TRUE, conf.level = 0.95)
r.who <- r1 %>% select(Anon, Test, TestType, Room, WHO5_Total) %>% 
  mutate(WHO5_100 = WHO5_Total *4)

r.who %>% group_by(Room, TestType) %>% 
  summarise(min=min(WHO5_100), max=max(WHO5_100), med=median(WHO5_100), 
            mean=mean(WHO5_100), sd=sd(WHO5_100))
r.who.tst %>% group_by(Room) %>% 
  summarise(min=min(Diff), max=max(Diff), med=median(Diff), 
            mean=mean(Diff), sd=sd(Diff))
```


# Heart rate

One example.
```{r}
px.c <- read.csv("data/garminxml/garmindats_ANON-Control.txt", header = FALSE, 
                 stringsAsFactors=FALSE)
px.w <- read.csv("data/garminxml/garmindats_ANON-Wood.txt", header = FALSE, 
                 stringsAsFactors=FALSE)

px.c$datetime <- as.POSIXct(px.c$V4)
px.w$datetime <- as.POSIXct(px.w$V4)
px.c <- px.c[complete.cases(px.c),]
px.w <- px.w[complete.cases(px.w),]
px.c <- px.c %>% select(V1, V2, V7, datetime) %>% 
  mutate(timedif = as.numeric(difftime(datetime, lag(datetime, default=0))))
px.c[1,5] <- 1
px.c <- px.c %>% mutate(elapsed = cumsum(timedif), elapsed.min = cumsum(timedif/60))

px.w <- px.w %>% select(V1, V2, V7, datetime) %>% 
  mutate(timedif = as.numeric(difftime(datetime, lag(datetime, default=0))))

px.w[1,5] <- 1
px.w <- px.w %>% mutate(elapsed = cumsum(timedif), elapsed.min = cumsum(timedif/60))

px <- bind_rows(px.c, px.w)
ggplot(data=px, aes(x=elapsed.min, y=V7, colour=V2)) + theme_bw() +
  geom_line() +
  scale_x_continuous(breaks=seq(0,75,5)) +
  labs(x="Time (min)",
       y=("Heart rate (BPM)"))

```

#Environment
```{r}
sessionInfo()
```

