# Load the necessary libraries
library(ggplot2)
library(tidyverse)
library(ggsci)
library(tidyr)
library(gridExtra)

# read in the data set
s <- read.csv(file.choose())

# make sure the right variables are formatted as dates
s$Date <- as.Date(s$Date, format = "%m/%d/%y")
s$Start.Time <- as.POSIXct(s$Start.Time, format = "%Y-%m-%d %H:%M%p")
s$End.Time <- as.POSIXct(s$End.Time, format = "%Y-%m-%d %H:%M%p")
s$Weekday <- weekdays(s$Date)

# let R know the activity variables are factors
s$Exercise <- as.factor(s$Exercise)
s$Write.cry <- as.factor(s$Write.cry)
s$Paint <- as.factor(s$Paint)
s$Meditate <- as.factor(s$Meditate)

# Plot activity assignment by day of the experiment
s <- s[order(s$Date),]
s$dayNumber <- c(1:nrow(s))
t <- pivot_longer(s, cols = c(12:15), values_to = "Assignment")
t$name <- factor(t$name, levels = (c("Meditate","Write.cry","Paint","Exercise")))

# First half of experiment
g1 <- ggplot(t[t$dayNumber <= 60,], aes(x = dayNumber, y = name, fill = Assignment)) +
  geom_tile(color = "black") +
  ylab("Activity") +
  xlab("Day Number") +
  ggtitle("Days 1 - 60") +
  scale_fill_manual(values=c("#b8b8b8", "#619CFF"))

# Second half of experiment
g2 <- ggplot(t[t$dayNumber > 60,], aes(x = dayNumber, y = name, fill = Assignment)) +
  geom_tile(color = "black") +
  ylab("Activity") +
  xlab("Day Number") +
  ggtitle("Days 61 - 113") +
  scale_fill_manual(values=c("#b8b8b8", "#619CFF"))

# Plot first and second half of experiment together in one plot
grid.arrange(g1, g2, nrow = 2, ncol = 1)

# Histogram of minutes awake
ggplot(s, aes(x=Minutes.Awake)) + 
  geom_histogram(fill = "#E77577", binwidth = 5) +
  xlab("Minutes Awake") +
  ylab("Count")

# Make weekday variable to plot according to weekday
s$Weekday <- factor(s$Weekday, levels = c("Saturday","Friday","Thursday",
"Wednesday","Tuesday","Monday","Sunday"))

# Plot minutes awake by weekday
ggplot(s, aes(x=Minutes.Awake, y = Weekday, color = Weekday)) + 
  xlab("Minutes Awake") +
  ylab("Weekay") +
  theme(legend.position = "none") +
  geom_jitter(position = position_jitter(seed = 2019, width = 0.05), size = 2, alpha = 0.4) +
  stat_summary(fun = mean, geom = "point", size = 5) +
  guides(fill = guide_legend(reverse = TRUE))

# make awakeasleep variable which is ratio of time awake to time asleep
s$awakeasleep <- s$Minutes.Awake/s$Minutes.Asleep

# Plot awakeasleep (awake:asleep ratio) by weekday
ggplot(s, aes(x=awakeasleep, y = Weekday, color = Weekday)) + 
  xlab("Minutes Awake / Minutes Asleep") +
  #scale_color_uchicago() +
  ylab("Weekay") +
  theme(legend.position = "none") +
  geom_jitter(position = position_jitter(seed = 2019, width = 0.05), size = 2, alpha = 0.4) +
  stat_summary(fun = mean, geom = "point", size = 5)

# Histogram of number of steps by assignment to exercise
ggplot(s, aes(x=Steps, fill=Exercise)) +
  geom_histogram( color="#e9ecef", position = 'identity', binwidth = 1000) +
  scale_fill_manual(values=c("#B8B8B8", "#619CFF")) +
  facet_grid(Exercise~.) +
  ylab("Count")

# Univariate analyses of exercise
mEx <- lm(Minutes.Awake~Exercise, data = s)
summary(mEx)

# Plot minutes awake by assignment to exercise
ggplot(s, aes(x = Exercise, y = Minutes.Awake, color = Exercise)) +
  geom_violin(trim = TRUE) + 
  coord_flip() +
  geom_jitter(shape=16, position=position_jitter(0.2), size = .5) +
  stat_summary(fun = mean, geom = "point", size = 4) +
  theme(legend.position = "none") +
  scale_color_manual(values=c("#B8B8B8", "#619CFF")) +
  xlab("Exercise Assignment") +
  ylab("Minutes Awake")

# Univariate analyses of paint
summary(lm(Minutes.Awake~Paint, data = s))

# Plot minutes awake by assignment to paint
ggplot(s, aes(x = Paint, y = Minutes.Awake, color = Paint)) +
  geom_violin(trim = TRUE) + 
  coord_flip() +
  geom_jitter(shape=16, position=position_jitter(0.2), size = .5) +
  stat_summary(fun = mean, geom = "point", size = 4) +
  theme(legend.position = "none") +
  scale_color_manual(values=c("#B8B8B8", "#619CFF")) +
  xlab("Painting Assignment") +
  ylab("Minutes Awake")

# Univariate analyses of write/cry
mWr <- lm(Minutes.Awake~Write.cry, data = s)
summary(mWr)
confint(mWr)

# Plot minutes awake by assignment to write/cry
ggplot(s, aes(x = Write.cry, y = Minutes.Awake, color = Write.cry)) +
  geom_violin(trim = TRUE, size = .5) + 
  coord_flip() +
  geom_jitter(shape=16, position=position_jitter(0.2), size = .5) +
  stat_summary(fun = mean, geom = "point", size = 4) +
  theme(legend.position = "none") +
  scale_color_manual(values=c("#B8B8B8", "#619CFF")) +
  xlab("Cry/Write.cry Assignment") +
  ylab("Minutes Awake")

# Univariate analyses of meditation
summary(lm(Minutes.Awake~Meditate, data = s)) 

# Plot minutes awake by assignment to meditate
ggplot(s, aes(x = Meditate, y = Minutes.Awake, color = Meditate)) +
  geom_violin(trim = TRUE) + 
  coord_flip() +
  geom_jitter(shape=16, position=position_jitter(0.2), size = .5) +
  stat_summary(fun = mean, geom = "point", size = 4) +
  theme(legend.position = "none") +
  scale_color_manual(values=c("#B8B8B8", "#619CFF")) +
  xlab("Meditation Assignment") +
  ylab("Minutes Awake")

# Checking univariate relationships between outcome and covariates to include in model
summary(lm(s$Minutes.Awake~s$Minutes.Asleep))
summary(lm(s$Minutes.Awake~s$Weekday))

# calculate lag correlations
# need the lag to be actual number of days between data points, and some days were skipped
s$diffjan1 <- difftime(s$Date, as.Date("01-01-2020", format = "%m-%d-%Y"), unit = "days")
daylag <- s[,which(colnames(s) %in% c("Minutes.Awake","diffjan1"))]
# figure out which days were skipped and add them to the data frame
a <- c(0:137)
extras <- as.data.frame(matrix(data = NA, nrow = length(a[!a %in% daylag$diffjan1]), ncol = 2))
extras$V2 <- a[!a %in% daylag$diffjan1]
colnames(extras) <- colnames(daylag)
daylag <- rbind(daylag, extras)
daylag <- daylag[order(daylag$diffjan1),]

# store the p-values, correlation coefficient, and lag for each lag value
ps <- rep(NA, 6)
cors <- rep(NA, 6)
lag <- c(1:6)
# store these 3 variables in a data frame
cordata <- as.data.frame(cbind(ps,cors,lag))
cortestlag <- function(lag) {
  x1 <- daylag$Minutes.Awake[-c(1:lag)]
  x2 <- daylag$Minutes.Awake[-c((nrow(daylag)-(lag-1)):nrow(daylag))]
  ct <- cor.test(x1,x2, use = "complete.obs")
  cordata$ps[lag] <- ct$p.value
  cordata$cors[lag] <- ct$estimate
  return(cordata)
}

# calculate these 3 values for each lag
for(i in c(1:6)) {
  cordata <- cortestlag(i)
}

which(cordata$ps < .05) # no significant correlations
cordata <- rbind(c(1,1,0),cordata)

# Plot the correlation data
with(cordata, plot(lag, cors, type = "h", lwd = 1, xlab = "Lag",
                   ylab = "Correlation"))
# Add line at x-axis
abline(h = 0)

# The final model!
fm <- lm(Minutes.Awake ~ Write.cry + Weekday + WFH + Minutes.Asleep, data = s)
# Plot the residuals vs. the fitted values to check model assumptions
plot(fm$residuals~fm$fitted.values, xlab = "Fitted Values", ylab = "Residuals")
# summarize the model and find the confidence intervals
summary(fm)
confint(fm)