---
layout:  post
title: "Bowel movements, losing weight while sleeping, and other questions 💩"
comments:  true
published:  true
author: "Zach Burchill"
date: 2020-01-20 10:00:00
permalink: /weight_analysis_pt_3/
categories: ["raspberry pi",R,health,weight,poop,sleep,hypotheses,data,"data science","data analysis"]
output:
  html_document:
    mathjax:  default
    fig_caption:  true
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE, warning=FALSE, message=FALSE}
knitr::opts_chunk$set(echo=FALSE, autocaption=TRUE,
                      warning = FALSE, message = FALSE,
                      dpi = 200)
options(kableExtra.latex.load_packages = FALSE)

library(tidyverse)
library(lubridate)
# library(testthat)
library(mgcv)
library(gratia)
 
```

```{r functions, echo=FALSE} 

is_empty_val <- function(v, ev=EMPTY_VALS) {
  map_lgl(v, ~.x %in% ev)
}
 
time_between <- function(first,second,units="seconds") {
  interval(first, second) %>%
    as.duration() %>%
    as.numeric(units=units)
}

are_close <- function(first, second, t) {
  as.duration(interval(first, second)) < t
}

# Currently unused
find_post_vals <- function(time_val, post_df, time_col, nt) {
  tcol <- enquo(time_col)
  
  vals <- post_df %>% 
    filter(!!tcol > time_val) %>%
    filter(are_close(time_val, !!tcol, nt))
  if (nrow(vals) > 2)
    warning("Multiple values in range, taking closest")
  if (nrow(vals) == 0) tibble()
  else vals %>% arrange(!!tcol) %>% {.[1,]}
}

# Right now, this doesn't take into account that calibrations can change over time
make_calibration_fn <- function(gt, measured) {
  ratio <- gt/measured
  fn <- approxfun(measured, ratio, rule=2)
  return(fn)
}

# Turn a cleaned df into a df for the GAMs
add_time_vars <- function(df) {
  df %>%
    mutate(ToD = hms::as_hms(time) %>% 
             as.numeric(),
           ToD_c = ToD %>% scale() %>% as.vector(),
           Hour = hour(time),
           DoW = wday(time, week_start = 1),
           DoW_fac = wday(time, label=T, week_start = 1) %>%
             `class<-`("factor"),
           ToW = DoW*3600*24+ToD,
           Date = date(time) %>% as.numeric(),
           Date_c = scale(Date) %>% as.vector(),
           Date_fac = as.factor(date(time)),
           Time_n = as.numeric(time),
           Time_c = scale(Time_n) %>% as.vector(),
           Time_n = Time_n - min(Time_n),
           WoY = week(time),
           FromStart = interval(min(time), time),
           nWeeks = trunc(time_length(FromStart,"week")))
} 

# Return the important predictions of a gam model in a nice format
get_gam_predictions <- function(model, df, time_term=Time_n) {
  tt <- enquo(time_term)
  # print(as_label(tt))
  smooth_names <- smooths(model)
  smooth_cols <- syms(smooth_names)
  time_smooth <- model %>%
    get_smooths_by_id(which_smooths(model, as_label(tt))) %>%
    pluck(1) %>% pluck("label") %>% sym()
  
  pred_df <- model %>% 
    predict(newdata = df, type = "terms") %>% 
    {mutate(as.data.frame(.), Intercept=attr(., "constant")[[1]]) } %>%
    as.data.frame() %>%
    transmute(all_pred = rowSums(.),
              smooth_pred = rowSums(select(., !!!smooth_cols, Intercept)),
              time_smooth_pred = !!time_smooth + Intercept)
  
  no_char_df <- df %>% drop_touchy_columns()
  
  derivs <- derivatives(model, newdata=no_char_df,
                        interval="simultaneous") %>%
    filter(smooth == as_label(time_smooth)) %>%
    transmute(sig_deriv = !map2_lgl(upper, lower, ~between(0,.y,.x)),
              pos_deriv = derivative > 0)
  
  pred_df %>%
    bind_cols(derivs)
}

# Get the regions of a smooth with a significant derivative
get_contiguous_regions <- function(df, n = 3, val = TRUE,
                                   time_col = Time_n,
                                   return_sum_df=FALSE) {
  tc <- enquo(time_col)
  df_sum <- df %>% 
    mutate(switch = ifelse(sig_deriv == lead(sig_deriv), 0, 1),
           ba = cumsum(switch)) %>%
    group_by(ba) %>%
    summarise(max = max(!!tc),
              min = min(!!tc) + 1e-7,
              is_pos = pos_deriv[2]==TRUE,
              numgood = sum(ifelse(sig_deriv==val, 1, 0))) %>%
    filter(numgood >= n) %>%
    mutate(group_id = seq_along(max))
  
  if (return_sum_df)
    return(df_sum)
  
  df %>%
    mutate(group_id = map(
      !!tc, 
      ~which(.x >= df_sum$min & .x <= df_sum$max)) %>%
        modify_if(is_empty, ~NA_integer_) %>%
        map_int(~.)) %>%
    filter(!is.na(group_id))
}

# For dealing with bugs in `gratia`
is_fac_or_can_inc <- function(x, by=1e-7) {
  is.factor(x) || can_be_incremented(x)
}
can_be_incremented <- function(x, by=1e-7) {
  tryCatch({x+by;TRUE}, error=function(e) FALSE)
}
drop_touchy_columns <- function(df) {
  good_cols <- vapply(df, is_fac_or_can_inc, logical(1L))
  df[,good_cols]
}

get_terms_data <- function(m, ...) {
  predict(m, ..., type="terms") %>% 
  {
    mutate(as.data.frame(.),
           Intercept=attr(., "constant")[[1]])
  }
}

get_all_pred_data <- function(m, og_data, newdata, 
                            time_term=Time_n,
                            min_contig=3) {
  tt <- enquo(time_term)
  smooth_names <- smooths(m)
  smooth_cols <- syms(smooth_names)
  time_smooth <- m %>%
    get_smooths_by_id(which_smooths(m, as_label(tt))) %>%
    pluck(1) %>% pluck("label") %>% sym()
  
  m_preds <-  get_terms_data(m) %>%
    select(tidyselect::matches("(s|ti|te|t2)\\(.+\\)"), `Intercept`) %>%
    transmute(Pred = rowSums(.)) %>%
    mutate(Resids = unlist(residuals(m))) %>%
    bind_cols(select(og_data, time)) %>%
    mutate(fuller_pred = Pred+Resids)
  
  m_even_time_se <- m %>%
    predict(newdata=newdata, type="terms", se.fit = TRUE) %>%
    pluck("se.fit") %>%
    as.data.frame() %>% 
    select(!!time_smooth) %>%
    transmute(se2 = !!time_smooth*2)
  
  m_even_preds <- m %>%
    get_gam_predictions(df=newdata, time_term=!!tt) %>%
    bind_cols(newdata) %>%
    bind_cols(m_even_time_se) %>%
    mutate(upper=time_smooth_pred+se2,
           lower=time_smooth_pred-se2)
  
  sig_regions_w_time <- get_contiguous_regions(m_even_preds,
                                               n = min_contig) %>%
    group_by(group_id) %>%
    summarise(max=max(time),
              min=min(time),
              is_pos=pos_deriv[2] == TRUE)
  
  list(actual_preds = m_preds,
       spline_data = m_even_preds,
       sig_regions = sig_regions_w_time)
}


prep_simple_df <- function(id_df, col,
                           pre_val, post_val,
                           og_df=df.clean) {
  col = enquo(col)
  id_df %>%
    left_join(og_df) %>%
    group_by(PID) %>%
    summarise(when = min(time),
              how_long = time_between(min(time), max(time), "minutes"),
              diff = median[!!col == pre_val] -
                median[!!col == post_val],
              weight = mean(median)) %>%
    mutate(when = scale(when) %>% as.vector()) %>%
    select(-PID)
}

```

```{r gathering functions}
# Breaks up an array of booleans into lists, 
#   chunking them where `bools==TRUE`.
break_at_true <- function(bools, 
                          break_before=FALSE,
                          exclude_singletons=TRUE) {
  
  if (exclude_singletons) comp <- `<`
  else                    comp <- `<=`
  
  if (break_before) adder <- 0
  else              adder <- 1
  
  nvi = which(bools)
  len = length(bools)
  
  if (length(nvi) == 0)
    return(list(1:len))
  
  nums <- map2(c(1, nvi + adder), c(nvi-1, len), ~ list(.x, .y))
  
  nums %>%
    keep(~ comp(.[[1]], .[[2]])) %>%
    map(~ .[[1]] : .[[2]])
}

# Instead of trying to match all possible before-and-after measurements in a data frame,
#   I can use my real-world knowledge of how these measurements work to narrow
#   the possible matches for each value.  This function breaks up a data frame
#   into smaller data frames, each containing only the possible matches for those rows.
#----------------------------
# @param is_special_val If there can be no intervening measurements between before-and-after
#   pairs, then set this parameter to an bool that evaluates TRUE if a row is the same type 
#   of measurement. E.g., for getting all the 'poo' measurements I use `df.clean$poo != "",`
#   since the tags that aren't marked for being about #2s have a value of `""` for the `poo`
#   column.
# @param duration The maximum duration (ie a `lubridate` `duration` object) apart two
#   measurements can be and still be linked
# @param time_col_name A string of the name of the time column in the dataframe, if 
#   `duration` is used
break_df <- function(df, is_special_val = NA,
                     duration = NA, time_col_name = NA) {
  
  if (!rlang::is_na(is_special_val)) {
    dfl <- break_at_true(!is_special_val,
                         break_before = FALSE,
                         exclude_singletons = TRUE) %>%
      map(~df[.x,])
  } else {
    dfl <- list(df)
  }
  
  if (!rlang::is_na(duration)) {
    dfl <- dfl %>%
      map(function(tiny_df) {
        within_bools <- are_close(lag(tiny_df[[time_col_name]]),
                                  tiny_df[[time_col_name]],
                                  duration)
        break_at_true(!within_bools,
                      break_before = TRUE,
                      exclude_singletons = TRUE) %>%
          map(~tiny_df[.x,])
      })
  } else {
    dfl <- dfl %>% map(list)
  }
  
  dfl %>%
    unlist(recursive=FALSE) %>%
    keep(~!is_empty(.))
}

# Given a score calculator function that calculates a score 
#   for two rows and an optional constraint function, it 
#   calculates the scores for each combination of row pairs, 
#   removing the ones that violate the constraint. 
# It returns a list of data frames with the two ids of the 
#   rows and their score.
calculate_scores <- function(test_df, calculator,
                             constraint=NA,
                             temp_index_col = Index,
                             ID_col = X) {
  ID_col = enexpr(ID_col) %>% as.character()
  # Adds a new index column
  index_col_name = enexpr(temp_index_col)
  test_df <- mutate(test_df, !!index_col_name := 1:nrow(test_df))
  
  # just making things shorter
  constr <- as_mapper(constraint)
  calc <- as_mapper(calculator)
  
  # First separate df into a list of its rows
  slices <- map(1:nrow(test_df), ~slice(test_df,.x))
  pairs <- combn(slices, 2, simplify = F)

  # Applies the constraints
  if (!rlang::is_na(constr))
    pairs <- pairs %>% keep(~constr(.x[[1]], .x[[2]]))
  
  # Calculate scores
  map(pairs, ~data.frame(FirstID  = .x[[1]][[ID_col]],
                         SecondID = .x[[2]][[ID_col]],
                         score = calc(.x[[1]], .x[[2]])))
}

# This function takes the output of `calculate_scores` and attempts
#   to return the paired rows that maximize the total score such that
#   each row is only ever present in one pair, removing those rows that
#   are not present in any pair.
# It then returns a data frame of IDs with pair ids.
# Note that the output is not at ALL guaranteed to be globally optimal,
#   nor is the speed optimal. I just wanted to get something quick
get_best_scores <- function(pairs, hi_scores_bad=FALSE,
                            ID_col = X) {
  ID_col = enexpr(ID_col)
  # Order the pairs by their score
  pairs <- pairs[order(map_dbl(pairs, ~.$score),
                       decreasing = !hi_scores_bad)]
  
  good_pairs <- list()
  # Go through the lists, picking the best pairs first
  #   and removing the ones that were already used
  while (length(pairs) > 0) {
    best <- pairs[[1]]
    newIDs <- c(best$FirstID, best$SecondID)
    good_pairs <- append(good_pairs, list(best))
    pairs <- pairs %>%
      keep(~!(.$FirstID %in% newIDs | .$SecondID %in% newIDs))
  }

  good_pairs %>%
    imap_dfr(~mutate(.x, PairID=.y)) %>%
    tidyr::gather("IDthing", !!ID_col,
                  FirstID:SecondID) %>%
    select(!!ID_col, PairID)
}

# test_df %>%
#   calculate_scores(~abs(time_between(.x$time, .y$time)),
#               ~.x$poo!=.y$poo & .x$clothing==.y$clothing) %>%
#   get_best_scores()

```

```{r constants}
# SIMPLE CONSTANTS
DATA_PATH <- "~/Desktop/"
DATA_FILE <- paste0(DATA_PATH, "for_r.csv")
EMPTY_VALS <- c("", NA) 
KG2LB <- 2.20462

```

```{r plotting constants and functions}
base_fig_height = 6
base_fig_width = 8

knitr::opts_chunk$set(fig.width = base_fig_width, 
                      fig.height = base_fig_height)

kg_lb_scale <- function(s, in_kg=TRUE) {
  if (in_kg) {
    trans <- ~KG2LB * .
    left  <- "(kg)"
    right <- "(lb)"
  } else {
    trans <- ~./KG2LB 
    right <- "(kg)"
    left  <- "(lb)"
  }
  scale_y_continuous(
    paste(s, left),
    sec.axis = dup_axis(trans=trans, name=paste(s, right)))
}

base_text_size = 20
shrinking_factor = 0.8
geom_text_size = base_text_size * 5/15 * shrinking_factor
theme_set(theme_bw(base_size = base_text_size))

stat_cols <- quos(sd, median, mean)
neg_stat_cols <- quos(-sd, -median, -mean)

med_y_scale  <- kg_lb_scale("median weight")
mean_y_scale <- kg_lb_scale("mean weight")

```

Although I started [this whole personal weight analysis thing]({{ site.baseurl }}{% post_url 2019-09-07-bluetooth_scale_intro %}) because I thought I was going to start dieting, I quickly realized that frequent fine-grained weight measurements could serve other, more _experimental_ purposes.

Finally, the scatological blog post you've all been waiting for! 💩💩💩

### Part 4: Bowel movements, losing weight while sleeping, and other questions

<!--more-->

_For background on the project this post stems from, check out [Part 1]({{ site.baseurl }}{% post_url 2019-09-07-bluetooth_scale_intro %}), [Part 2]({{ site.baseurl }}{% post_url 2020-01-02-weight_analysis_pt_1 %}), and [Part 3]({{ site.baseurl }}{% post_url 2020-01-04-weight_analysis_pt_2 %})  of this series._

```{r load and process data}
raw_df <- read.csv(DATA_FILE) %>% 
  as_tibble() %>%
  mutate(time = as.POSIXct(time, origin="1970-01-01")) %>%
  # HACK! should be changed soon (GG1)
  mutate(time = with_tz(time, tzone="America/New_York")) %>%
  select(-weight_values, everything(), weight_values) %>%
  filter(matched=="True") %>% 
  # This is for backwards compatability:
  { 
    if ("username" %in% names(.))
      filter(., grepl("zach", tolower(username)) | username=="")
    else
      .
  } %>%
  { 
    if (!("is_metric" %in% names(.)))
      mutate(., is_metric="")
    else
      .
  } %>% 
  filter(error=="") %>%
  mutate(clothing = factor(clothing,
                           levels = c("naked", "some_clothes","full_clothes")))

# Some of the older data needs to be cleaned up a bit.
# My measurement code now avoids this necessity
mezzies_need_cleaning <- raw_df %>%
  filter(git_hash=="") %>%
  tidyr::separate_rows(weight_values) %>%
  mutate(w=as.numeric(weight_values)) %>%
  group_by(X, time) %>%
  filter(any(w < 75)) %>%
  summarise() %>% pluck("X")

df.clean <- raw_df %>%
  filter(X %in% mezzies_need_cleaning) %>%
  tidyr::separate_rows(weight_values) %>%
  mutate(w=as.numeric(weight_values)) %>%
  filter(w > 60) %>%
  group_by(X, time) %>%
  filter(abs(scale(w)) < 3) %>%
  group_by_at(vars(-w, -weight_values, -sd:-mean)) %>%
  summarise(mean   =   mean(w),
            sd     =     sd(w),
            median = median(w),
            weight_values = paste(w, collapse=", ")) %>%
  bind_rows(filter(raw_df, !(X %in% mezzies_need_cleaning))) %>%
  ungroup()

calib_df <- df.clean %>%
  filter(!is_empty_val(is_calib)) %>%
  select(X,median,mean,sd,
         calib_weight, is_metric,
         ID, time, weight_values) %>%
  mutate(calib_weight = ifelse(is_metric==TRUE,
                               calib_weight, 
                               calib_weight/KG2LB))

calibrate <- make_calibration_fn(calib_df$calib_weight,
                                 calib_df$median)

df.clean <- df.clean %>%
  mutate(calib_ratio = calibrate(median)) %>%
  mutate_at(vars(mean,median,sd), ~ . * calib_ratio) %>%
  filter(is_empty_val(is_calib) | is_calib==FALSE) %>%
  filter(!is_empty_val(clothing))


```

```{r group tags}

# An example of how I used these functions to look at sleep/wake data:
sleepytime_ids <- df.clean %>%
  arrange(time) %>%
  break_df(is_special_val = df.clean$sleepwake != "", 
 # I decided instances of 'good' sleep/wake pairs shouldn't have measurements 
 #  between them, so using the `is_special_val` will chunk the data frame on 
 #  cases where there aren't sleep/wake measurements, narrowing down the possible
 #  'pair' space immensely
           duration = duration(24, "hours"),
 # Sleep/wake measurements need to be less than 24 hours apart. Including this as
 #  the `duration` parameter also helps reduce the possible pair space by breaking
 #' up the data frame
           time_col_name = "time") %>%
 # Now that the data frame has been broken up into a list of more manageable
 #  smaller dfs, I calculate the pairing scores for the rows in each
  map(function(z) 
    calculate_scores(z, 
                     # Penalize for being longer than 8 hours or less than six
                     ~max(time_between(.x$time, .y$time, "minutes")-8*60, 0) +
                       max(6*60-time_between(.x$time, .y$time, "minutes"), 0),
                     # They have to be different sleep values, starting with 'bedtime'
                     ~.x$sleepwake != .y$sleepwake & 
                       .x$sleepwake == "bedtime" &
                       # They have to have the same clothing on
                       .x$clothing == .y$clothing &
                       # There can't be intervening sleep values between a pair
                       .x$Index == .y$Index - 1 &
                       # It can't be longer than 16 hours or shorter than 3
                       time_between(.x$time, .y$time, "hours") < 16 &
                       time_between(.x$time, .y$time, "hours") > 3
    )) %>%
 # I throw away empty data frames
  keep(~length(.x) > 0) %>%
 # I get the best-scoring pair out of each sublist
  map(~get_best_scores(.x, hi_scores_bad = TRUE)) %>%
  imap_dfr(~mutate(.x, PID = paste0(.y,"-",PairID,"-sleep")) %>% select(-PairID))

poo_ids <- df.clean %>%
  ungroup() %>%
  break_df(df.clean$poo != "", 
           duration = duration(120, "minutes"),
           "time") %>%
  map(function(z) 
    calculate_scores(z, 
                     ~abs(time_between(.x$time, .y$time)),
                     ~.x$poo != .y$poo & .x$clothing == .y$clothing &
                       .x$poo == "pre_poop")) %>%
  keep(~length(.x) > 0) %>%
  map(~get_best_scores(.x, hi_scores_bad = TRUE)) %>%
  imap_dfr(~mutate(.x, PID = paste0(.y,"-",PairID,"-poo")) %>% select(-PairID))

pee_ids <- df.clean %>%
  ungroup() %>%
  break_df(df.clean$pee != "", 
           duration = duration(120, "minutes"),
           "time") %>%
  map(function(z) 
    calculate_scores(z, 
                     ~abs(time_between(.x$time, .y$time)),
                     ~.x$pee != .y$pee & .x$clothing == .y$clothing &
                       .x$pee == "pre_pee")) %>%
  keep(~length(.x) > 0) %>%
  map(~get_best_scores(.x, hi_scores_bad = TRUE)) %>%
  imap_dfr(~mutate(.x, PID = paste0(.y,"-",PairID,"-pee")) %>% select(-PairID))



shower_ids <- df.clean %>%
  arrange(time) %>%
  break_df(df.clean$shower != "", 
           duration = duration(1, "hours"),
           time_col_name = "time") %>%
  map(function(z) 
    calculate_scores(z, 
                     # Penalize for being longer than 8 hours or less than six
                     ~max(time_between(.x$time, .y$time, "minutes")-45, 0) +
                       max(7-time_between(.x$time, .y$time, "minutes"), 0),
                       # They have to be different sleep values, starting with pre-shower
                       ~.x$shower != .y$shower & 
                         .x$shower == "pre_shower" &
                         # They have to have the same clothing
                         .x$clothing == .y$clothing &
                         # There can't be intervening sleep values
                         .x$Index == .y$Index - 1
      )) %>%
    keep(~length(.x) > 0) %>%
    map(~get_best_scores(.x, hi_scores_bad = TRUE)) %>%
    imap_dfr(~mutate(.x, PID = paste0(.y,"-",PairID,"-shower")) %>% select(-PairID))

all_ids <- map(list(poo_ids, pee_ids, sleepytime_ids, shower_ids), ~.$X) %>% unlist()

collapsed_df <- map_dfr(list(poo_ids, pee_ids, sleepytime_ids, shower_ids),
                  ~left_join(.x, df.clean)) %>%
                  { bind_rows(., anti_join(df.clean, ., by="X"))} %>%
  mutate(PID = ifelse(is.na(PID), paste(seq_along(PID), "ZA"), PID)) %>%
  group_by(PID, clothing) %>% #filter(n() > 2)
  summarise_at(vars(time, !!!stat_cols, X), mean) %>%
  ungroup()

```


# Note to all possible future employers:

**Stop reading now.** 

By all that is good and true in this world, close out of this tab and just forget about this post.  

I promise you that there is nothing _offensive_ in this post, nothing even politically incorrect.  But I cannot promise you that you will not be grossed out by this post.  You have every right to be. You might even come to the conclusion that my undertaking here crosses the borders of "juvenile" or that I have the mind of an eleven year-old.

And if you did, you would be right.

But in my defense, I believe the taboo nature of the topics I will be touching upon[^1] is increasingly a relic of the past. For, in the immortal words of that hallowed book: *"Everyone poops."*

![Everyone Poops, by Tarō Gomi.](https://upload.wikimedia.org/wikipedia/en/a/a4/Everyone_Poops.jpg){:width="400px"}

<p class = "figcaption">The book I remember reading in my grandparents' house: <a href="https://en.wikipedia.org/wiki/Everyone_Poops"><i>Everyone Poops</i>, by Tarō Gomi</a>. (Taken from Wikipedia)</p>

## Why am I like this?

About a day after planning my weight loss measuring journey, I realized that the scale system I was coming up with could also measure very _short_ term weight losses/fluctuations.  Within minutes, I realized that I could use a before-and-after system to weigh my own bodily functions.

Yes, _bowel_ movements.

I mean, it makes sense! The heaviest short-term weight fluctuations are going to be from eating food and getting _rid_ of that food. It's only natural for my mind to go there!

While, yes, I _did_ consider the utility of developing a more... "detailed" scatological corpus for wider medical use, after speaking with a few friends involved in medical fields, I decided that the need was just not there.

But understanding more about one's excretory habits has vital importance for understanding one's own personal health, and _haven't you ever wondered how much your poop weighs?_

Anyway, I came up with roughly three areas of interest about my weight:

## My questions:

1. How much weight, if any, do I lose in my _sleep_?  How do my sleeping habits affect my overall weight?

2. How much does my average bowel movement weigh? How do the different _types_ of excreta compare? Do these habits have any bearing on my overall weight?

3. How do other habits impact fluctuations in my weight?

See? Only *one* of my questions has to do with poop!

### Technical background

I have developed a local [Flask](https://en.wikipedia.org/wiki/Flask_(web_framework)) website that interfaces with the [Bluetooth monitoring system]({{ site.baseurl }}{% post_url 2019-09-07-bluetooth_scale_intro %}) that's running on a RaspPi in my house, letting me "tag" different weight recordings with meta-data via browser.[^2]  I'm pretty dang satisfied with what I've built[^3], but the one major kink in my pipeline is the fact that the Bluetooth/scale part can only run in Python 2, while everything else runs in Python 3.[^4]

For better or for worse, the system I've created essentially work with two stacks: a stack of weight data and a stack for the meta-data, and as weights and tags are put onto these stacks, it attempts to link weights to tags. So far this has worked great.

But to measure short-term weight changes, I needed to record my weight both before and after a change, and supply meta-data that would let me know what kind of change I was measuring. Finding the difference in the before-and-after measurements would let me know exactly how much weight I lost. While my system for linking meta-data to the individual weights works great, my system for linking the two before-and-after measurements is a little more of a pain. Mostly because of how forgetful I am. 

I forget to tag some weight measurements, I forget to weigh myself again to complete a "before-and-after" measurement, etc. When I came up with the design that I did, I figured that these problems would arise, but that I would have the chops to deal with them in "post" and clean the data up afterwards. Leaving all the hard work for future me. But let's put that on hold for now.

## Question 1: Do I lose weight in my sleep? {#question_1}

One of the questions I thought would be interesting to ask with this system is whether my weight changes from when I go to sleep to when I wake up, and if so, by how much.  

To measure this change, I measured my weight immediately before and after I slept, and tagged each of these weights with that information. According to [science blogger Derek Muller](https://www.youtube.com/watch?v=lL2e0rWvjKI) and [NPR](https://www.npr.org/sections/krulwich/2013/06/19/193556929/every-night-you-lose-more-than-a-pound-while-youre-asleep-for-the-oddest-reason), we lose about a pound of weight in our sleep per night. Is that what _my_ data says?

### The baby solution

The simplest solution would be to filter all the weight sleep/wake measurements and do a _t_-test to see if the two category means are different. 

This is the easy way---I don't have to link 'before' and 'after' measurements---I can count on the fact that they're mostly probably paired up, and if there's a difference, the noise might not be _too_ much of a problem.  But before we do that, let's look at the data:

```{r, collapse=TRUE, echo=FALSE}
temp_df <- df.clean %>%
  filter(clothing=="naked", sleepwake!="") %>%
  group_by(sleepwake)

# temp_df %>% 
#   summarise(lbs = mean(median)*KG2LB)

temp_df %>%
  ggplot(aes(x=sleepwake, y=median, color=sleepwake)) +
  geom_boxplot() +
  geom_point(position = position_jitter(0.1), alpha=0.2, size=3) +
  scale_color_discrete(guide=FALSE) +
  med_y_scale + xlab("Time") +
  ggtitle("Before and after sleeping","naked measurements only") 

bedtimes <- temp_df[temp_df$sleepwake=="bedtime",]$median * KG2LB
wakeups <-  temp_df[temp_df$sleepwake=="wakeup",]$median * KG2LB
```

Although the data is certainly noisy, it looks like I'm _relatively_ close to that mark (~`r temp_df %>% summarise(lbs = mean(median)*KG2LB) %>% summarise(diff=abs(diff(lbs))) %>% pluck("diff") %>% round(3)` lbs).  But is this weight loss significant?

```{r}
t.test(bedtimes, wakeups, paired=FALSE)
```

_Womp womp!_  Looks like it isn't!  But can we do better?  What if I'm able to link the 'before' and 'after measurements?

### Fixing it in 'post'

I wanted to create an entirely separate blog post about how I approached this problem, but I think I've gone on enough about this project as-is. I'll just say that my approach to how best to link these values in post was a lesson in how "worth it" it is to approach problems in a principled, thoughtful way: by considering the problem as a whole, I was able to whip up some very general code that solved all my problems in one go.

I really wanted ~~to brag about it~~ use it as a introduction to the design choices that programmers have to make, but that will have to be sometime later.  If you're interested in the code, you can look at the R (commented) code in the [source file of this post](https://raw.githubusercontent.com/burchill/burchill.github.io/master/_source/2020-01-20-weight_analysis_pt_3.Rmd).[^5] 

Long story short, I was able to break the problem up into easier, more manageable pieces, and use what I know about how these measurements were taken to link them together.

### Big boy solution

```{r}
ya <- sleepytime_ids %>%
  left_join(df.clean) %>%
  arrange(PID)
paired_bedtimes <- ya[ya$sleepwake=="bedtime",]$median * KG2LB
paired_wakeups <- ya[ya$sleepwake=="wakeup",]$median * KG2LB
```

Ok, so what happens if we use my code to match up the corresponding sleep/wake measurements? We can essentially turn it into a paired _t_-test, for one. Each 'before' measurement (i.e., when I went to bed) has a corresponding 'after' measurement (when I woke up).

```{r}
ttt <- t.test(paired_bedtimes, paired_wakeups, paired=TRUE)
ttt
```

Not only does using these paired measurements make the difference significant, but we see that the difference is almost exactly what was previously reported: `r round(ttt$estimate, 3) %>% as.vector()` lbs!  Props to Derek!


## Question 2: Potty talk

Alright, this is more of a _series_ of questions than a single one. I'm _generally_ curious about what the data looks like, but there are some hard questions as well, like comparing the weights of the type of bowel movements, etc. But before we get to the specific questions, let's just see what the data looks like, starting with the data for the 'poop' measurements. 

**Note on my methodology:** I measured the entire weight of a bowel movement as a single measurement, so a "#2" included both the weight of the solid _and_ liquids released. I'm interested in the topic, but not enough to do one after the other, measuring myself in between!


```{r get simple dfs}

poo_df <- prep_simple_df(poo_ids, poo, "pre_poop","post_poop")
pee_df <- prep_simple_df(pee_ids, pee, "pre_pee","post_pee")
sleep_df <- prep_simple_df(sleepytime_ids, sleepwake, "bedtime", "wakeup")
shower_df <- prep_simple_df(shower_ids, shower, "pre_shower", "post_shower")

```

```{r huh, echo=TRUE, fig.cap="A quick visual summary of the 'poo' data. <code>when</code>, <code>how_long</code>, <code>diff</code>, and <code>weight</code> correspond to the date, the time between measurements, the difference in kg, and the mean kg, respectively. The numbers in the upper-right panels are the Pearson correlations."} 
psych::pairs.panels(poo_df)
```

Nothing too weird, none of the correlations are really that striking.

### Pizzles vs. deuces

So let's see how the number ones compare to the number twos. Before running this experiment, I expected that the average doodoo would outweigh the average piddle.  Although I'm pretty sure more liquid weight makes its way through our bodies each day, I generally micturate more than I defecate. And since  the way I measured a number two _included_ the accompanying number one, I thought this hypothesis would be a slam dunk.

```{r, fig.cap="Boxplot of weight of my recorded BMs by type. The thick dark lines represent the medians and the red lines represent the means."}

combo_df <- bind_rows("number one" = pee_df,
                      "number two" = poo_df,
                      .id = "Type") %>%
  filter(diff >= 0) %>%
  group_by(Type) 

stat_combo_df <- combo_df %>%
  summarise(median=median(diff),
            mean = mean(diff)) %>%
  tidyr::gather(stat,value,median:mean)

combo_df %>%
  ggplot(aes(y = diff, x = Type)) + 
  # geom_violin() +
  geom_boxplot( outlier.alpha=0) +
  geom_point(position = position_jitter(0.1)) +
  ggtitle("How heavy are my toilet trips?") +
  stat_summary(fun.data = function(x) data.frame(y = mean(x), ymax = mean(x), ymin = mean(x)), 
               geom="errorbar", color="red",
               width=0.75) +
  kg_lb_scale("weight change") +
  xlab("bowel movement type") +
  ggrepel::geom_label_repel(aes(label = paste(round(value, 3), "kg"), 
        y = value,
        color = stat),
    data = stat_combo_df,
    min.segment.length=0,
    force = 500,
    box.padding = 4,
    max.iter = 10000,
    seed=2,
    size = geom_text_size) +
  scale_color_manual(guide=FALSE,
                     values=c("red","black"))
```

Surprisingly, my recorded tinkles outweigh my recorded dookies in both mean _and_ median!  This difference isn't statistically significant, but it is numerically intriguing.

Thinking about it, I have a few potential explanations.

First, it might be true. Water (essentially what urine is) by volume is pretty heavy. Try picking up a bucket of water, it's heavier than you might think. Maybe a full bladder is heavier than a full lower intestine.  What makes this theory pretty weak is that, as we have all seen, feces most often sink to the bottom of the toilet bowl---they must be denser than water/urine.

More likely, it could be a problem with how I _collected_ the data: I have the feeling that I was more likely to remember to weigh myself when I had to pee really badly, but had an easier time remembering to weigh myself before any poop. With a poop on its way, I would know that I might have to 'settle' in for a longer toilet visit, so I would take the time to measure beforehand.

To solve this last issue, I could be *incredibly* diligent about weighing myself every time I use the restroom. I'd have to stay home for a while, but it could be possible to pull off. If _every_ bowel movement is measured, this bias should decrease.  However, this line of reasoning led me to a different thought.

### Maybe the wrong question?

The frequency with which we poop is not (fully) determined by how much we eat: many studies say that what is healthy ranges from [three times a day to three times a week](https://www.vice.com/en_us/article/j54xv7/how-often-should-you-poop).[^6] In my layman's opinion, this frequency is much more stable than the frequency of urination, which I feel is more determined by things such as how busy we are, etc. And given that the weight per widdle is a function of consumption by frequency, does talking about "individual" urinations make that much sense?   

Sure, the frequency of these bowel movements could be interesting in its own right, but I would have to be much more diligent about recording all of them. This just comes back to the second issue I brought up. I approached this question with the lazy understanding that I could just record whichever bowel movements I happened to remember about, and that these ~~stool~~ samples would be representative of an "average" bowel movement. But in reality, thinking about an "average" BM isn't that useful, especially for urination. Overall "output" rates and frequency seem much more applicable.


## Question 3: Lots of questions

Time for a lightning round!

### Question 3.1: Does a shower make me lighter?

I'm pretty sure scientists say that you lose dead skin cells, hair, grime, etc. in the process of showering.  But is that change noticeable? I recently started weighing myself before and after a shower to see.

Below are the changes in weight, before-and-after, in pounds:

```{r}

shower_df %>%
  arrange(diff) %>% 
  pluck("diff") %>% {.*KG2LB} %>% round(2)

```

**Answer: Inconclusive.**  While the mean is slightly positive, there seems to be enough noise in the data that I can't tell as of now!


### Question 3.2: Longer bathroom breaks = more poo/pee? 

Do I stay in the bathroom longer for bigger BMs?
 
```{r q32, fig.height = base_fig_height * 0.8 }

diff_dfs <- bind_rows("poo" = poo_df,
                      "pee" = pee_df %>% filter(diff > 0),
                      "sleep" = sleep_df,
                      "shower" = shower_df,
                      .id = "df_type") %>%
  tidyr::nest(-df_type) %>%
  mutate(p.value = data %>% 
           map_dbl(~cor.test(.$how_long, .$diff)$p.value),
         cor = data %>% 
           map_dbl(~cor.test(.$how_long, .$diff)$estimate),
         R2 = cor^2,
         lm = data %>% 
           map(~lm(diff ~ how_long, data=.x) %>%
                 broom::tidy() %>% select(term, estimate)),
         slope = map_dbl(lm, ~.[.$term=="how_long",]$estimate),
         int = map_dbl(lm, ~.[.$term=="(Intercept)",]$estimate)) %>%
    tidyr::unnest(data)


diff_dfs %>%
  filter(df_type %in% c("poo","pee")) %>%
  ggplot(aes(x = how_long, y = diff)) +
  geom_point() +
  # geom_abline(aes(intercept=int, slope=slope),
  #             color="red") +
  geom_smooth(method="lm") +
  facet_wrap(~df_type, scales="free") +
  zplyr::geom_abs_label(ypos=0.75, xpos=0.7,
                        aes(label = paste0("p=", format.pval(p.value, digits = 2))),
                        size=geom_text_size) +
  zplyr::geom_abs_label(ypos=0.85, xpos=0.7,
                        aes(label = paste0("R=", round(cor,2))),
                        size=geom_text_size) +
  kg_lb_scale("BM weight") +
  xlab("minutes spent on toilet")

```

**Answer: No.** In fact, it looks like I stay longer for <em>smaller</em> poops and pees. Must be messing around on my phone!

### Question 3.3: Longer sleep = more weight lost?

If Derek and NPR are right, we lose most of the weight lost in sleep through respiration. If this hypothesis is correct, we would expect that longer sleep = more respiration = more weight lost. Is this borne out in the data?

```{r, out.width="40%", fig.height = base_fig_height * 0.8, fig.width = base_fig_width/2 }

diff_dfs %>%
  filter(df_type %in% c("sleep")) %>%
  ggplot(aes(x = how_long/60, y = diff)) +
  geom_point() +
  geom_smooth(method="lm") +
  facet_wrap(~df_type, scales="free") +
  zplyr::geom_abs_label(ypos=0.75, xpos=0.7,
                        aes(label = paste0("p=", format.pval(p.value, digits = 2))),
                        size=geom_text_size) +
  zplyr::geom_abs_label(ypos=0.85, xpos=0.7,
                        aes(label = paste0("R=", round(cor,2))),
                        size=geom_text_size) +
  kg_lb_scale("weight lost in sleep") +
  xlab("hours slept")
  
```

**Answer: Probably!** The effect is close to significance, which I'll take, and the trend is in the right direction!    
          
### Question 3.4: How well can GAMs predict before-and-after effects?

In my [previous post]({{ site.baseurl }}{% post_url 2020-01-04-weight_analysis_pt_2 %}) about using GAMs on this data set, I pondered how well a GAM could capture the weight changes from the before-and-after measurements when you throw everything into a model together.

I decided to test this, with a kind of jury-rigged model comparison. Assuming that the methods I've used earlier in this post capture the "ground truth" for the 'average' number 1, number 2, and weight lost in sleep,[^7] I can compare how well different GAMs accurately capture these average changes.

I compared 25 values of `k` for the time smooth, for GAMs made with both the `gamm()` and `gam()` functions. I also compared models with time as a linear predictor ("k=0"), models with the time smooth fit to the "collapsed" data residualized out ("k=c"), and a simple linear regression model.


```{r grouped_gam_stuff}

collapsed_gam_df <- collapsed_df %>% 
  add_time_vars() %>%
  mutate(Y = median) %>%
  mutate(clothing = clothing %>% 
           as.character() %>% 
           factor(levels = 
                    c("naked", "some_clothes", "full_clothes"))) %>%
  arrange(time) %>%
  filter(time < max(time)) # just removes a weird point

collapsed_dummy_df <- data.frame(
  time = seq(from = min(collapsed_gam_df$time), 
             to = max(collapsed_gam_df$time), 
             length = nrow(collapsed_gam_df))) %>%
  add_time_vars() %>%
  mutate(clothing = collapsed_gam_df$clothing[[1]])

m <- gamm(Y ~ 1 + s(Time_n, k=15) + 
            s(DoW, k=4, bs="cc") + 
            s(ToD, k=12, bs="cc") + 
            clothing,
          data = collapsed_gam_df)
# print(plot(mgcViz::getViz(m$gam), allTerms = T), pages = 1)

```

```{r regular_gam_data}

big_gam_df <- df.clean %>%
  add_time_vars() %>%
  mutate(Y = median) %>%
  mutate_at(vars(clothing, poo, pee, sleepwake, shower),
            ~as.character(.)) %>%
  mutate(clothing = factor(clothing, 
                           levels = 
                             c("naked","some_clothes","full_clothes")),
         poo = ifelse(poo=="", "normal", poo) %>%
           factor(levels=c("normal", "pre_poop", "post_poop")),
         pee = ifelse(pee=="", "normal", pee) %>%
           factor(levels=c("normal","pre_pee","post_pee")),
         sleepwake = ifelse(sleepwake=="","normal", sleepwake) %>%
           factor(levels=c("normal","bedtime","wakeup")),
         shower = ifelse(shower=="","normal", shower) %>%
           factor(levels=c("normal","pre_shower","post_shower"))) %>%
  arrange(time) %>%
  filter(time < max(time)) # just removes a weird point

pred_d <- m$gam %>%
  get_all_pred_data(collapsed_gam_df, big_gam_df, Time_n)

big_gam_df$NewY <- big_gam_df$Y - pred_d$spline_data$time_smooth_pred

```

```{r total_trash_code, cache=TRUE}
# This is all garbage. Please forgive me

get_coefs <- function(m) {
  gam_res_data <- data.frame(
    estimate = summary(m)$p.coeff,
    terms = names(summary(m)$p.coeff)
  ) %>%
    as_tibble() %>%
    mutate(se = summary(m)$se[1:nrow(.)])
  gam_res_data %>%
  {
    bind_rows(
    "sleep" = get_coeff_df(., "sleepwake","bedtime","wakeup"),
    "poo"  = get_coeff_df(., "poo","pre_poop","post_poop"),
    "pee"  = get_coeff_df(., "pee","pre_pee","post_pee"),
    "shower"  = get_coeff_df(., "shower","pre_shower","post_shower"),
    .id="coef"
  ) %>%
    mutate(edf = summary(m)$edf[[1]])
  }
}

get_coeff_df <- function(coef_df,name,pre,post) {
  pre_name = paste0(name,pre)
  post_name = paste0(name,post)
  coef_df %>%
    summarise(est = estimate[terms==pre_name] - 
                estimate[terms==post_name],
              se = se[terms==pre_name] + 
                se[terms==post_name]) %>%
    mutate(upper = est+se*1.96,
           lower = est-se*1.96)
  
}

# for the linear model
get_coeff_df2 <- function(coef_df,name,pre,post) {
  pre_name = paste0(name,pre)
  post_name = paste0(name,post)
  coef_df %>%
    summarise(est = estimate[term==pre_name] - 
                estimate[term==post_name],
              se = std.error[term==pre_name] + 
                std.error[term==post_name]) %>%
    mutate(upper = est+se*1.96,
           lower = est-se*1.96)
  
}

different_ks <- seq(3,50,by=2)#  c(3,5,8,10,12,13,14,15,20,30,50)

lowest_gamm <- gamm(Y ~ 1 + s(Time_n, k=3) + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)
noest_gamm <- gamm(Y ~ 1 + Time_n + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)

collapsed_gamm <- gamm(NewY ~ 1 + 
                        s(DoW, k=4, bs="cc") + 
                        
                        clothing + sleepwake + poo + pee + shower,
                      select=TRUE,
                      data = big_gam_df)
collapsed_gam <- gam(NewY ~ 1 + 
                        s(DoW, k=4, bs="cc") + 
                        
                        clothing + sleepwake + poo + pee + shower,
                      select=TRUE,
                      data = big_gam_df)
lm_model <- lm(Y ~ 1 + Time_n + clothing + sleepwake + poo + pee + shower, data=big_gam_df)

lm_dat <- broom::tidy(lm_model) %>%
{
  bind_rows(
    "sleep" = get_coeff_df2(., "sleepwake", "bedtime", "wakeup"),
    "pee" = get_coeff_df2(., "pee", "pre_pee", "post_pee"),
    "poo" = get_coeff_df2(., "poo", "pre_poop", "post_poop"),
    "shower" = get_coeff_df2(., "shower", "pre_shower", "post_shower"),
    .id="coef"
    )
}

lowest_gam <- gam(Y ~ 1 + s(Time_n, k=3) + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)
noest_gam <- gam(Y ~ 1 + Time_n + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)

reg_gamm_mods <- map(different_ks, function(x) {
  gamm(Y ~ 1 + s(Time_n, k=x, fx=T) + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)
})

reg_gam_mods <- map(different_ks, function(x) {
  gam(Y ~ 1 + s(Time_n, k=x, fx=T) + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)
})


free_reg_gamm_mods <- map(different_ks, function(x) {
  gamm(Y ~ 1 + s(Time_n, k=x, fx=F) + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)
})

free_reg_gam_mods <- map(different_ks, function(x) {
  gam(Y ~ 1 + s(Time_n, k=x, fx=F) + 
            s(DoW, k=4, bs="cc") + 
            
            clothing + sleepwake + poo + pee + shower,
            select=TRUE,
          data = big_gam_df)
})

all_gamm_data <- reg_gamm_mods %>%
  map2_dfr(different_ks,
           ~get_coefs(.x$gam) %>% 
             mutate(k=as.character(.y))) %>%
  bind_rows(get_coefs(lowest_gamm$gam) %>%
              mutate(k="1"),
            get_coefs(noest_gamm$gam) %>%
              mutate(k="0"),
            get_coefs(collapsed_gamm$gam) %>%
              mutate(k="c"),
            lm_dat %>% mutate(k="lm")
            ) 
all_gam_data <- reg_gam_mods %>%
  map2_dfr(different_ks,
           ~get_coefs(.x) %>% 
             mutate(k=as.character(.y))) %>%
  bind_rows(get_coefs(lowest_gam) %>%
              mutate(k="1"),
            get_coefs(noest_gam) %>%
              mutate(k="0"),
            get_coefs(collapsed_gam) %>%
              mutate(k="c")
            ) 

all_free_gamm_data <- free_reg_gamm_mods %>%
  map2_dfr(different_ks,
           ~get_coefs(.x$gam) %>% 
             mutate(k=as.character(.y))) %>%
  bind_rows(get_coefs(lowest_gamm$gam) %>%
              mutate(k="1"),
            get_coefs(noest_gamm$gam) %>%
              mutate(k="0"),
            get_coefs(collapsed_gamm$gam) %>%
              mutate(k="c"),
            lm_dat %>% mutate(k="lm")
            ) 
all_free_gam_data <- free_reg_gam_mods %>%
  map2_dfr(different_ks,
           ~get_coefs(.x) %>% 
             mutate(k=as.character(.y))) %>%
  bind_rows(get_coefs(lowest_gam) %>%
              mutate(k="1"),
            get_coefs(noest_gam) %>%
              mutate(k="0"),
            get_coefs(collapsed_gam) %>%
              mutate(k="c")
            ) 

actual_vals <- bind_rows(
  "poo" = poo_df %>% 
    summarise(mean=mean(diff), sd=sd(diff)),
  "pee" = pee_df %>% 
    summarise(mean=mean(diff), sd=sd(diff)),
  "sleep" = sleep_df %>% 
    summarise(mean=mean(diff), sd=sd(diff)),
  "shower" = shower_df %>% 
    summarise(mean=mean(diff), sd=sd(diff)),
  .id = "coef"
)

all_dfs <- bind_rows("gamm()"=all_gamm_data,
          "gam()" =all_gam_data,
          .id="type") %>%
  left_join(actual_vals) %>%
  # arrange(k) %>%
  mutate(in_95 = pmap_lgl(list(mean,lower,upper),
                          ~between(..1,..2,..3)),
         in_1se = pmap_lgl(list(mean,est,se),
                          ~between(..1,..2-..3,..2+..3))) %>%
  mutate(k = factor(k, 
    levels = c("lm","0", "c","1", as.character(different_ks)))) %>%
  mutate(edf = ifelse(k %in% c("0","c","lm"),
                      NA, edf)) %>%
  mutate(g = ifelse(k=="lm",1,0)) %>%
  mutate(type=factor(type,levels=c("gamm()","gam()"))) %>%
  mutate(fake_k = case_when(!(k %in% c("0","c","lm", "1")) ~ as.numeric(as.character(k)),
                            k == "lm" ~ -20,
                            k == "0" ~ -10,
                            k == "c" ~ -5,
                            TRUE ~ 0))

all_free_dfs <- bind_rows("gamm()"=all_free_gamm_data,
          "gam()" =all_free_gam_data,
          .id="type") %>%
  left_join(actual_vals) %>%
  mutate(in_95 = pmap_lgl(list(mean,lower,upper),
                          ~between(..1,..2,..3)),
         in_1se = pmap_lgl(list(mean,est,se),
                          ~between(..1,..2-..3,..2+..3))) %>%
  mutate(k = factor(k, 
    levels = c("lm","0", "c","1", as.character(different_ks)))) %>%
  mutate(edf = ifelse(k %in% c("0","c","lm"),
                      NA, edf)) %>%
  mutate(g = ifelse(k=="lm",1,0)) %>%
  mutate(type=factor(type,levels=c("gamm()","gam()"))) 
```


```{r, fig.cap="'lm' indicates results from a linear model, '0' is from models with time as a linear predictor, '1' is from models with the fewest possible basis functions, and 'c' is from models with the effect of time from a different model residualized out. For the time smooths of these GAMS, <code>fx</code> is set to <code>TRUE</code>."}

all_dfs %>%
  filter(coef != "shower") %>%
  ggplot(aes(x=fake_k, y=abs(1-est/mean), color=coef, group=coef)) +
  geom_point() +
  geom_rect(ymin = -Inf, ymax=Inf,
            xmin=-Inf,   xmax=-15, 
            fill="grey90",
            data = data.frame(type = "gamm()", g = 0),
            inherit.aes = FALSE) +
  geom_point() +
  geom_line() +
  facet_wrap(~type,scales="free_x") +
  scale_size_continuous(guide=FALSE) +
  # theme(axis.text.x = element_text(angle = 45,hjust=1)) +
  scale_color_discrete("effect") +
  scale_y_continuous("percent error\n(1-predicted/gt)",
                     labels = scales::percent_format(accuracy = 1)) +
  ggtitle("How well do GAMs capture the effects?") +
  scale_x_continuous(
    breaks = c( -20,  -10, -5,  0,  10,20,30,40,50),
    labels = c("lm", "0","c","1","10","20","30","40","50")) +
  xlab("'k' of the time smooth")
    
```



First, I should say that "ground truth" estimates for each effect are all within the 95% CI for each GAM. So in some sense, the answer seems to be "yes."  However, as we see, we can get estimates that are ~40% off in some models, so it's worth being cautious.

Perhaps unsurprisingly, as the number of basis functions increases, the accuracy of the estimations seem to get better and better, up to a point.[^8]

However, in the figure above, I fixed the `k` term for each smooth, instead of letting `mgcv()` automatically determine the best number of basis functions with `k`-1 as the upper bound. 

When I repeat what I did above, but let `mgcv` decide the optimal number of basis functions (represented below via the estimated degrees of freedom), we immediately see that there is *much* less 'thrashing' around of the estimated effect sizes.
          
```{r, fig.cap="How well the GAMs capture the effects versus the number of basis functions of the time smooth (represented via the estimated degrees of freedom for the smooth) as <code>k</code> increases. Note how the <code>gamm()</code> function tops out around 25."}

all_free_dfs %>%
  filter(coef != "shower") %>%
  filter(k %in% as.character(3:50)) %>%
  mutate(`function`=type) %>%
  ggplot(aes(x=edf, y=abs(1-est/mean), color=coef,
             linetype=`function`)) +
  geom_point() +
  geom_line() +
  scale_size_continuous(guide=FALSE) +
  scale_color_discrete("effect") +
  scale_y_continuous("percent error",
                     labels = scales::percent_format(accuracy = 1)) +
  ggtitle("How well do GAMs capture the effects?",
          "smooth parameters estimated automatically") +
  xlab("e.d.f. of time smooth")
```

**Lesson learned: don't fix `k` unless you really know what you're doing!** 

Interestingly though, it seems that `gamm()` is much more conservative with the amount of wiggliness it bestows the time spline---even when `k` is set to `50`, the edf barely goes above 25, while the edf of `gam()` go much higher.  If anybody has any clue _why_, let me know!


<hr />
<br />

## Source Code:

> [`2020-01-20-weight_analysis_pt_3.Rmd`](https://raw.githubusercontent.com/burchill/burchill.github.io/master/_source/2020-01-20-weight_analysis_pt_3.Rmd)

If you want to see how I did what I did in this post, check out the source code of the R Markdown file!

### Footnotes:

[^1]: Eh, bad phrasing.

[^2]: If you want access to the git repo for this project, drop me a line on Twitter and I'll add you as a collaborator. Otherwise, I'm keeping it private for the time being. 😉

[^3]: I want to brag about it so much. I've set up a subdomain to redirect me to the landing page of my local site: `scale.zachburchill.ml` (only available when you're connected to my home WiFi network). I also have a webhook system that will send push notifications directly to my phone. I feel so 133t.

[^4]: See the issues in [my previous post]({{ site.baseurl }}{% post_url 2019-09-07-bluetooth_scale_intro %}). If not for these issues, I could get a _really_ sexy system going easy-peasy with some real-time dashboards and whatnot.

[^5]: The important functions are defined in the `gathering functions` chunk and there's a commented example in the `group tags` chunk.

[^6]: I'll go on record that I think both of those frequencies are _bizarre_, by the way.

[^7]: I'm making a lot of assumptions here.

[^8]: I originally had a 'time-of-day' spline in the GAMs, but it ended up eating up the variance of the sleep effect (since I have a somewhat regular bedtime), and led to very unstable estimates.