#load packages
library(googlesheets)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggthemes)
library(RColorBrewer)
library(scales)

## Multiplot function: ####
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

# Import ####

log <- read.csv("emma_log.csv", stringsAsFactors = F)
log$date <- as.Date(log$date, "%m/%d/%Y")
log[is.na(log)] <- 0
tot <- sum(log$hour)

# Breast Preferences ####

# Most popular breast by time
round((sum(log$right.time, na.rm = T))/60,2)
round((sum(log$left.time, na.rm = T))/60,2)

dat2 = data.frame(count=c(107.37, 147.73), category=c("Right", "Left"))

dat2$fraction = dat2$count / sum(dat2$count)
dat2 = dat2[order(dat2$fraction), ]
dat2$ymax = cumsum(dat2$fraction)
dat2$ymin = c(0, head(dat2$ymax, n=-1))

d1 <- ggplot(dat2, aes(fill=category, ymax=ymax, ymin=ymin, xmax=4, xmin=3)) +
  geom_rect(colour="grey30") +
  coord_polar(theta="y") +
  xlim(c(0, 4)) +
  labs(fill="") +
  theme_bw() +
  theme(panel.grid=element_blank()) +
  theme(axis.text=element_blank(),
        legend.position="bottom", 
        legend.direction="horizontal") +
  theme(axis.ticks=element_blank()) +
  labs(title="Total time spent on each side",
       subtitle = "Baby spent 147.7 hours on the left side, \n58% of her total time feeding.")
ggsave("total_time_side.jpg", width = 5)

# Popular breast by count
sum(log$right, na.rm = T)
sum(log$left, na.rm = T)

dat = data.frame(count=c(390, 478), category=c("Right", "Left"))

dat$fraction = dat$count / sum(dat$count)
dat = dat[order(dat$fraction), ]
dat$ymax = cumsum(dat$fraction)
dat$ymin = c(0, head(dat$ymax, n=-1))

d2 <- ggplot(dat, aes(fill=category, ymax=ymax, ymin=ymin, xmax=4, xmin=3)) +
  geom_rect(colour="grey30") +
  coord_polar(theta="y") +
  xlim(c(0, 4)) +
  labs(fill="") +
  theme_few() +
  theme(panel.grid=element_blank()) +
  theme(axis.text=element_blank(),
        legend.position="bottom", 
        legend.direction="horizontal") +
  theme(axis.ticks=element_blank()) +
  labs(title="Total Feeds on each side",
       subtitle = "478 (55%) of baby's feeds were on the left side. \nLeft side average: 18.5; right side average: 16.5.")
ggsave("total_feeds_side.jpg", width = 5)

# Average time per breast
round(mean(log$right.time[log$right.time!=0], na.rm = T),2)
round(mean(log$left.time[log$left.time!=0], na.rm = T),2)

median(log$right.time[log$right.time!=0], na.rm = T)
median(log$left.time[log$left.time!=0], na.rm = T)

# Feeding Trends ####

# Total number of hours feeding
log$feed.time <- log$right.time + log$left.time
fd.t <- round((sum(log$feed.time))/60,2)
fd.t / tot

dat3 = data.frame(count=c(255.1, 1185), category=c("Feeding", "Not Feeding"))

dat3$fraction = dat3$count / sum(dat3$count)
dat3 = dat3[order(dat3$fraction), ]
dat3$ymax = cumsum(dat3$fraction)
dat3$ymin = c(0, head(dat3$ymax, n=-1))

d3 <- ggplot(dat3, aes(fill=category, ymax=ymax, ymin=ymin, xmax=4, xmin=3)) +
  geom_rect(colour="grey30") +
  coord_polar(theta="y") +
  xlim(c(0, 4)) +
  labs(fill="") +
  theme_few() +
  theme(panel.grid=element_blank()) +
  theme(axis.text=element_blank(),
        legend.position="bottom", 
        legend.direction="horizontal") +
  theme(axis.ticks=element_blank()) +
  labs(title="Proportion of time spent feeding",
       subtitle = "Baby spent 255 hours feeding, 18% of total time \nduring the period tracked.")
ggsave("total_time_feeding.jpg", width = 5)

mplot <- multiplot(d1,d2,d3,cols = 3)


# Total feed time per day

daily <- aggregate(feed.time ~ date, data = log, sum)
daily$hours <- round(daily$feed.time/60,2)

ggplot(filter(daily, hours != 0), aes(date, hours)) +
  geom_bar(stat = "identity", fill = "pink1") +
  labs(x = "", y = "Feeding hours per day") +
  ggtitle("Total feeding time per day", 
          subtitle="After some initial issues, baby got better at feeding. Good job baby! \nMore recently, I have it on good authority that feeds are getting more efficient.") +
  theme_few()
ggsave("daily_feed_time.jpg")

# Most popular feeding times: total hours feeding
hourly <- aggregate(feed.time ~ time.s, data = log, sum)
hourly <- hourly %>% arrange(-feed.time)
hourly$hours <- hourly$feed.time/60
head(hourly,3)
tail(hourly,5)

ggplot(hourly, aes(time.s, hours)) +
  geom_bar(stat = "identity", fill = "pink1") +
  labs(x = "Hour ", y = "Total hours across all days") +
  scale_x_continuous(breaks = c(1:24),labels = c(1:24)) +
  ggtitle(
    "Total time feeding during each hour of the day", 
    subtitle="Most popular feed times are 7pm, 4pm, and 11pm; least popular feed times are 1am, 6am, and 2am. \nBaby seems to like to feed most in the afternoon and evenings. Good job baby!") +
  theme_few()
ggsave("hourly_feed_time_total.jpg")

# Most popular feeding times: average length of feed

hourly.avg <- log %>% filter(feed.time != 0)
hourly.avg <- round(aggregate(feed.time ~ time.s, data = hourly.avg, mean),2)

ggplot(hourly.avg, aes(time.s, feed.time)) +
  geom_bar(stat = "identity", fill = "pink1") +
  labs(x = "Hour", y = "Average feed length in minutes") +
  scale_x_continuous(breaks = c(1:24),labels = c(1:24)) +
  ggtitle(
    "Average feed length in minutes", 
    subtitle="Baby's favorite times of day are less apparent here: in general, average feed length is fairly consistent \nacross hours of the day. 4pm, 9pm, and 10pm show the longest average feeds.") +
  theme_few()
ggsave("hourly_feed_time_average.jpg")

# Count of feeds per day
feed.count <- log %>% filter(feed.time !=0) %>%
  count(date, hour) %>% select(-hour)

# Length of feeds over time

r.feeds <- log %>%
  filter(right.time != 0) %>%
  select(date, right.time)

l.feeds <- log %>%
  filter(left.time != 0) %>%
  select(date, left.time)

#feeds <- merge(r.feeds, l.feeds, by = "date")

ggplot() +
  geom_point(data = r.feeds, aes(date, right.time), color = "blue") +
  geom_point(data = l.feeds, aes(date, left.time), color = "red")

# Diaper Trends ####
# Total number of poops
sum(log$poo, na.rm = T)
sum(log$pee, na.rm = T)

# Total number of diaper changes (poop or pee)
log <- mutate(log, diaper = poo + pee)
log$diaper[log$diaper==2] <- 1
sum(log$diaper)

# Total number of feeds
log <- mutate(log, feed = right + left)
log$feed[log$feed==2] <- 1
sum(log$feed)

# Total events
log <- mutate(log, event = diaper + feed)
log$event[log$event==2] <- 1
sum(log$event)
round(sum(log$event) / tot,2)

dat4 = data.frame(count=c(572, 647), category=c("Event Hours", "No Event Hours"))

dat4$fraction = dat4$count / sum(dat4$count)
dat4 = dat4[order(dat4$fraction), ]
dat4$ymax = cumsum(dat4$fraction)
dat4$ymin = c(0, head(dat4$ymax, n=-1))

ggplot(dat4, aes(fill=category, ymax=ymax, ymin=ymin, xmax=4, xmin=3)) +
  geom_rect(colour="grey30") +
  coord_polar(theta="y") +
  xlim(c(0, 4)) +
  labs(fill="") +
  theme_few() +
  theme(panel.grid=element_blank()) +
  theme(axis.text=element_blank(),
        legend.position="bottom", 
        legend.direction="horizontal") +
  theme(axis.ticks=element_blank()) +
  labs(title="Proportion of time doing something",
       subtitle = "Event hours defined as containing a diaper change or a feed. \nWe have had 572 event hours, approximately 45% of our total hours.")
ggsave("event_totals.jpg", width = 5)

# Average diaper changes per day (poop or pee)

d.diaper <- aggregate(diaper ~ date, data = log, sum)
d.poo <- aggregate(poo ~ date, data = log, sum)
d.pee <- aggregate(pee ~ date, data = log, sum)

ggplot() +
  geom_point(data = d.diaper, aes(date, diaper), color = "black", size = 4) +
  geom_point(data = d.poo, aes(date, poo), color = "tan4", size = 3) +
  geom_point(data = d.pee, aes(date, pee), color = "yellow2", size = 2) +
  scale_y_continuous(breaks = c(0,2,4,6,8,10,12)) +
  labs(x = "", y = "Total number per day") +
  ggtitle("Diaper events per day", 
          subtitle="Events are color coded; black is total diaper changes, you can figure out the rest") +
  theme_bw()
ggsave("daily_diapers.jpg")

# Number of uninterrupted hours of sleep per night

night <- log %>% filter(time.s %in% c(11,12,1,2,3,4,5,6)) %>%
  filter(feed.time == 0) %>%
  count(date, hour)

ggplot(filter(night, n <= 7), aes(date, n)) +
  geom_bar(stat = "identity", fill = "navyblue") +
  scale_y_continuous(breaks = c(1:8)) +
  labs(x = "", y = "Hours of sleep per night") +
  ggtitle("Hours of sleep per night", 
          subtitle="Night defined as 11pm to 7am. Estimates based on the number of hours each night with no feed or diaper event. \nWe had some pretty rough nights in the middle of April. Mean hours of sleep per night is approximately 5.2. \nNote: not continuous/uninterrupted sleep!") +
  theme_few()
ggsave("nightly_sleep.jpg")

# Average uninterrupted hours per night
round((mean(night$n, na.rm = T)),2)
round((median(night$n, na.rm = T)),2)

night2 <- log %>% filter(time.s %in% c(11,12,1,2,3,4,5,6))
y <- rle(night2$feed.time)
uninterrupted.sleep <- y$lengths[y$values==0]
uninterrupted.sleep
round(mean(uninterrupted.sleep),2)
median(uninterrupted.sleep)

# Additional analysis ####

daily.feed <- aggregate(feed ~ date, data = log, sum)
daily.night.feed <- log %>% filter(time.s %in% c(11,12,1,2,3,4,5,6))
daily.night.feed <- aggregate(feed ~ date, data = daily.night.feed, sum)

ggplot() +
  geom_line(data = filter(daily.feed, feed > 0), aes(date,feed), color = "pink1") +
  geom_line(data = filter(daily.night.feed, feed > 0), aes(date,feed), color = "navyblue") +
  scale_y_continuous(breaks = c(0:16)) +
  labs(x = "", y = "Count of feeds") +
  ggtitle("Feeding trends, day and night", 
          subtitle="The pink line represents a count of all feeds per day; the blue line represents a count of nightly feeds.") +
  theme_bw()
ggsave("feedingtrends.jpg")

weekly.feed <- aggregate(feed.time ~ X, data = log, sum)
weekly.feed$hours <- round(weekly.feed$feed.time/60,2)
weekly.night.feed <- log %>% filter(time.s %in% c(11,12,1,2,3,4,5,6))
weekly.night.feed <- aggregate(feed.time ~ X, data = weekly.night.feed, sum)
weekly.night.feed$hours <- round(weekly.night.feed$feed.time/60,2)

ggplot() +
  geom_line(data = weekly.feed, aes(X,hours), color = "pink1") +
  geom_line(data = weekly.night.feed, aes(X,hours), color = "navyblue") +
  scale_x_continuous(breaks = c(0:10)) +
  scale_y_continuous(breaks = c(5,10,15,20,25,30,35,40,45)) +
  labs(x = "Week", y = "Hours feeding") +
  ggtitle("Weekly feeding trends", 
          subtitle="The pink line shows the total number of hours spent feeding per week; the blue line shows the nightly hours \nspent feeding per week. A significant downward trend is evident (even bearing in mind that week 9 is incomplete.") +
  theme_bw()
ggsave("feedingtrends.jpg")