In bsaul/elktoeChemistry: Data and code for Elktoe Chemistry project

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  echo    = FALSE,
  message = FALSE,
  warning = FALSE
)

library(elktoeChemistry)
library(gt)
library(ggplot2)
# library(geex)
library(geepack)
library(lme4)
library(quantreg)
# library(MASS)
library(dplyr)
analysis_data  <- readRDS(params$inputs$analysis)
# analysis_data <- readRDS("data/analysis_data.rds")

site_colors <- 
  c("LiTN 1" = "#99d8c9",
    "LiTN 2" = "#41ae76",
    "LiTN 3" = "#005824",
    "Tuck 1" = "#fdae6b",
    "Tuck 2" = "#f16913",
    "Tuck 3" = "#8c2d04",
    "Baseline" = "#525252")

dt <-
  analysis_data(
    contrast = "all",
    test_data_FUN = function(data) data,
    elements  = "all", 
    signals   = c("mad_10"),
    group_by_valve = TRUE,
    transect_opts  = list(
      .layers    = c("ipx", "ncr", "psm", "pio", "opx"),
      .min_n_obs = 5L,
      .outer_buffer = 15,
      .inner_buffer = 15
    )
  )

layer2 <- function(layer, annuli){
   case_when(
      layer == "ncr" & annuli == "A" ~ "ncr_new",
      layer == "ncr" & annuli != "A" ~ "ncr_old",
      # layer %in% c("ipx", "opx")     ~ "epx",
      TRUE ~ layer
  )
}

Summary of elements by species/site/layer

specimen_summaries <- 
  summarize_specimens(dt, modify_layer = layer2) %>%
  filter(site != "Baseline") 

adt <- 
  specimen_summaries %>%
  filter(species == "Arav") %>%
  left_join(
    dt$data[[1]] %>%
      filter(site != "Baseline") %>%
      distinct(id, river, dead, moribund, baseline_weight, baseline_volume, 
              final_status, prop_weight_lost),
    by = "id"
  ) %>%
  filter(!is.na(final_status)) %>%
  mutate(
    element =  case_when(
      stringr::str_detect(element, "^(Zn|Cu|Mg)") ~
        paste0(stringr::str_replace(element, "_ppm_m", "["), "]"),
      TRUE ~ stringr::str_remove(element, "_ppm_m.*")
    )
  )

yy <- adt %>% 
  filter(!is.na(baseline_weight), !is.na(baseline_volume)) 

xx <- 
  yy %>%
  filter(!(element %in% c("Cd", "U"))) %>%
  group_by(element, layer) %>%
  mutate(
    zz = residuals(lm(log(l1) ~ baseline_weight + baseline_volume + site)),
    v = log10(l1),
    z = (v - mean(v))/sd(v)
  )

ggplot(
  data = xx,
  aes(x = final_status, y = z)
) + 
  geom_jitter(
    size = 0.2,
    width = 0.25,
    color = "grey50"
  ) + 
  # geom_point(size = 0.2) +
  geom_point(
    data = xx %>%
      group_by(element, layer, final_status) %>%
      summarize(z = median(z)),
    aes(color = final_status)
  ) +
  geom_line(
    data = xx %>%
      group_by(element, layer, final_status) %>%
      summarize(z = median(z)) %>%
      ungroup(),
    aes(x = final_status, y = z, group = layer)
  ) +
  # stat_smooth(method = "lm", se = TRUE) +
  facet_grid(layer ~ element, scales = "free")

Models

mdt <- 
  adt %>%
  select(river, drawer, site, layer, dead, final_status, l1, prop_weight_lost,
         baseline_weight, baseline_volume) %>%
  arrange(id) %>%
  ungroup() %>%
  left_join(
    select(mussels_wide, id, wd = buoyant_weight_g_d),
    by = "id"
  ) %>%
  group_by(id) %>%
  group_nest() %>%
  mutate(id = 1:n()) %>%
  tidyr::unnest(cols = data) %>%
  group_by(element, layer) %>%
  mutate(
    # id     = as.integer(as.factor(id)),
    site   = factor(site),
    drawer = factor(drawer),
    layer  = factor(layer, levels = c("ipx", "ncr_new", "ncr_old", "psm", "pio", "opx"),
                    ordered = TRUE),
    l1s    = (l1 - mean(l1))/sd(l1)
  ) %>%
  ungroup() %>%
  arrange(id)

mdat <- 
  mdt %>%
  group_by(element) %>%
  group_nest()

summarise_gee <- function(m){
  summary(m) %>%
  purrr::pluck("coefficients") %>% 
  { x <- .
    x %>% 
      mutate(parameter = rownames(x)) %>%
      select(
        estimate = Estimate,
        std_error = Std.err,
        p_value   = `Pr(>|W|)`,
        everything()
      )
  }
}

xx <- mdat$data[[14]]

xx <- xx %>% 
  group_by(site, layer) %>%
  mutate(
    wd = if_else(dead == 1, -baseline_weight, wd),
    l1_c_s = l1/exp(mean(log(l1))),
    l1c  = l1 - mean(l1)
    # l1_c_s = l1 - mean(l1)
    # l1_c_s = l1_c_s - min(log(l1))
  )

ggplot(
  xx,
  aes(x = dead, y = log(l1))) +
  geom_point() +
  stat_smooth(method = "lm") + 
  facet_grid(layer ~ site, scales = "free")


library(qif)
summary(geeglm(dead ~ site*layer + layer:log(l1), data = xx, id = id, family = binomial))
summary(qif(dead ~ site + layer:log(l1), data = xx))
summary(qif(wd ~ site + layer + layer:log(l1), data = xx))
summary(qif(wd ~ site + layer + layer:log(l1), data = xx, corstr = "AR-1", 
            invfun = "ginv"))

xx %>%
  ggplot(
    aes(x = log(l1))
  ) + geom_histogram() +
  facet_grid(site~ layer)

models <-
  mdat %>%
  mutate(
    m1 = purrr::map(
      .x = data,
      .f = ~ {
        geeglm(
          dead ~ log(l1):layer + site + site:layer,
          data = .x,
          family = binomial(link = "logit"),
          id = id)
      }),
    m2 = purrr::map(
      .x = data,
      .f = ~ {
        geeglm(
          dead ~ log(l1):layer + site + site:layer,
          data = .x %>% 
            filter(site != "Tuck 1") %>% 
            mutate(site = factor(site)),
          family = binomial(link = "logit"),
          id = id)
      }),
  ) %>%
  select(element, m1, m2) 

res <- 
  models %>%
  tidyr::pivot_longer(
    names_to = "model",
    cols = c(m1, m2)
  ) %>%
  mutate(
    m_summary = purrr::map(value, summarise_gee)
  ) %>%
  select(-value) %>%
  tidyr::unnest(cols = m_summary) 

pdt <-
  res %>% 
  filter(grepl("l1", parameter)) %>%
  mutate(parameter = stringr::str_replace(parameter, "log\\(l1\\):layer", "")) %>%
  filter(!(element %in% c("Sr"))) %>%
  group_by(model, element, parameter) %>%
  mutate(
    layer = factor(parameter, levels = c("ipx", "ncr_new", "ncr_old", "psm", "pio", "opx"),
                    ordered = TRUE),
    conf_lo   = estimate - 2*std_error,
    conf_hi   = estimate + 2*std_error,
    cross0    = p_value < 0.05
  ) 

Estimated coefficients from a logistic GEE model predicting death

ggplot(
  data = filter(pdt, model == "m1", std_error < 0.5),
  aes(x = element, y = estimate, color = cross0)
) + 
  geom_hline(yintercept = 0) +
  geom_point(size = 0.5) + 
  geom_segment(
    aes(xend = element, y = conf_lo, yend = conf_hi)) + 
  scale_color_manual(
    values = c("black", "red"),
    guide = FALSE
  ) +
  facet_grid(layer ~ ., scale = "free") +
  theme(
    axis.title.x = element_blank(),
    axis.text.x  = element_text(size = 8, angle = 45, hjust = 1) 
  )

summary_by_site <-
  mdt %>%
  mutate(
    wd = if_else(dead == 1, -baseline_weight, wd)
  ) %>%
  group_by(river, element, layer, site) %>%
  summarise(
    l1m  = exp(mean(log(l1))),
    dead = mean(dead),
    bw   = sum(baseline_weight),
    lw   = sum(wd),
    pw   = bw - lw,
    rl   = (bw - pw)/bw
  ) 

plot_summary <- function(layer){
  dt <- 
  summary_by_site %>%
    filter(layer == !! layer)

  ggplot(
    data = dt,
    aes(x = l1m, y = rl, color = river)
  ) +
  ggtitle(layer) +
  geom_point() +
  stat_smooth(method ="lm", se = FALSE, color = "blue") +
  stat_smooth(
    data = dt %>% filter(site != "Tuck 1"),
    method ="lm", se = FALSE,
    color = "black"
  ) +
  scale_x_continuous("Geometric mean of means of specimens") +
  scale_y_continuous("Relative loss of site biomass") +
  facet_wrap(. ~ element, ncol = 5, scales = "free") +
  theme(
    axis.text.y = element_text(size = 8),
    axis.text.x = element_text(size = 8)
  )
}

Geometric mean of means within site by loss of biomass within site

(Beware of ecological fallacy with this view!)

• Black = best fit excluding Tuck 1
• Blue= best fit

plot_summary("ipx")

plot_summary("ncr_new")

plot_summary("ncr_old")

plot_summary("psm")

plot_summary("pio")

plot_summary("opx")

mdt2 <- 
  mdt %>%
  mutate(
    wd = if_else(dead == 1, -baseline_weight, wd)
  ) 

summarize_rlm <- function(m){
  summary(m) %>%
  purrr::pluck("coefficients") %>%
  { x <- .
    x %>% 
      as.data.frame() %>%
      mutate(parameter = rownames(x)) %>%
      select(
        estimate = Value,
        std_error = `Std. Error`,
        t   = `t value`,
        everything()
      )
  }
}


summarize_lmer <- function(m){
  summary(m) %>%
  purrr::pluck("coefficients") %>%
  { x <- .
    x %>% 
      as.data.frame() %>%
      mutate(parameter = rownames(x)) %>%
      select(
        estimate = Estimate,
        std_error = `Std. Error`,
        t   = `t value`,
        everything()
      )
  }
}

summarize_rq <- function(m){
  m %>%
  purrr::pluck("coefficients") %>%
  {
    x <- .
    as.data.frame(x) %>%
      mutate(parameter = rownames(x),
             tau = m$tau)
  }
}

mods <- 
  mdt2 %>%
  group_by(element, layer) %>%
  group_nest() %>%
  mutate(
    mlm = purrr::map(
      .x = data,
      .f = ~ {
        MASS::rlm(wd ~ site + log(l1),
                  data = .x, 
                  maxit = 100)
      }
    ),
    rlm = purrr::map(
      .x = data,
      .f = ~ {
        lmer(
          wd ~ site + log(l1) + (1|site),
          data = .x)
      }
    ),
    rqm = purrr::map(
      .x = data,
      .f = ~ {
        rq(
          wd ~ log(l1) + site, tau = seq(0.2, 0.8, by = 0.1),
          data = .x
        )
      }
    ),
    ms = purrr::map(mlm, summarize_rlm),
    re = purrr::map(rlm, summarize_lmer),
    rqs = purrr::map(
      .x = rqm, 
      .f = ~ { summary(.x) %>%  purrr::map_dfr(summarize_rq) }
    )
  )

rq_res <- 
  mods %>% 
  select(element, layer, rqs) %>%
  tidyr::unnest(cols = rqs) %>%
  filter(parameter == "log(l1)") %>%
  select(
    estimate = coefficients,
    conf_lo  = `lower bd`,
    conf_hi = `upper bd`,
    everything()
  )

plot_rq <- function(layer){
  rq_res %>%
  filter(layer == !! layer) %>%
  ggplot(
    aes(x = tau, y = estimate)
  ) + 
  ggtitle(layer) +
  geom_hline(yintercept = 0) +
  geom_point(size = 0.5) +
  geom_ribbon(
    aes(ymin = conf_lo, ymax = conf_hi),
    alpha = 0.3
  ) + 
  geom_line() +
  scale_y_continuous("Quantreg Estimated Coefficient of mean concentration") +
  scale_x_continuous("Quantile") +
  facet_wrap(
    element ~ ., ncol = 5,  scales = "free"
  )
}

Quantile regression estimating effect of log(l1) on weight lost

plot_rq("ipx")

plot_rq("ncr_new")

plot_rq("ncr_old")

plot_rq("psm")

plot_rq("pio")

plot_rq("opx")

rq_by_site <- 
 mdt2 %>%
 group_by(element, layer) %>%
  group_nest() %>%
  mutate(
   rqm = purrr::map(
      .x = data,
      .f = ~ {
        rq(
          wd ~ -1 + site + log(l1):site, tau = c(0.25, 0.5, 0.75),
          data = .x
        )
      }
    ),
    rqs = purrr::map(
      .x = rqm, 
      .f = ~ { summary(.x) %>%  purrr::map_dfr(summarize_rq) }
    )
  )

rq_res <- 
  rq_by_site %>% 
  select(element, layer, rqs) %>%
  tidyr::unnest(cols = rqs) %>%
  filter(grepl("log\\(l1\\)", parameter)) %>%
  mutate(
    site     = substr(parameter, 13, 18),
  ) %>%
  select(
    estimate = coefficients,
    conf_lo  = `lower bd`,
    conf_hi = `upper bd`,
    everything()
  )

rq_res %>%
  filter(layer == "psm") %>%
  ggplot(
    aes(x = tau, y = estimate)
  ) + geom_point() +
  facet_wrap(element ~ ., ncol = 5, scale = "free")

plot_l1_by_wd <- function(layer){
  ggplot(
    data = mdt2 %>% filter(layer == !! layer),
    aes(x = log(l1), y = wd, color = site)
  ) + 
  ggtitle(layer) +
  geom_point(size = 0.5) +
  scale_color_manual(
    values = site_colors
  ) +
  scale_y_continuous(
    "Pre to post difference in weight"
  ) + 
  scale_x_continuous(
    "Log(mean in layer)"
  ) +
  geom_quantile(
    method = "rq", color = "blue",
    quantiles = c(0.5) ) + 
  facet_wrap(~ element, ncol = 5, scales = "free")
}

Scatterplots of log(l1) by weight loss

plot_l1_by_wd("ipx")

plot_l1_by_wd("ncr_new")

plot_l1_by_wd("ncr_old")

plot_l1_by_wd("psm")

plot_l1_by_wd("pio")

plot_l1_by_wd("opx")

ee <- function(data, md, mw){
  # browser()
  psi1 <- grab_psiFUN(md, data = data)
  psi2 <- grab_psiFUN(mw, data = data)
  p1 <- length(coef(md))
  p2 <- length(coef(mw))
  d  <- all(data$dead == 1)

  function(theta){
    # print(data$id)
    # browser()
    c(psi1(theta[1:p1]), 
      { if(!d) psi2(theta[(p1+1):(p1 + p2)]) else rep(0, p2) } )
      # (!d * 1) * ( )

  }
}

library(mediation)

summarise_mediation <- function(x){
  els <- c("d0", "d0.ci", "d0.p", "d1", "d1.ci", "d1.p",
           "d.avg", "d.avg.ci", "d.avg.p",
           "z.avg", "z.avg.ci", "z.avg.p",
           "n0", "n1", "tau.coef", "tau.ci", "tau.p")
  tibble(!!! unlist(x[els]), use.names = TRUE)
}

med_dat <- 
  mdat %>%
  tidyr::unnest(cols = data) %>%
  mutate(
    wd1   = if_else(dead == 1, -baseline_weight, wd),
    ll1   = log(l1),
    treat = (river == "Little Tennessee"),
    # treat = !(site == "Tuck 1"),
    comparison = "rivers"
  ) %>%
  {
    dt <- .
    purrr::map_dfr(
      .x = c("Tuck 2", "Tuck 3", "LiTN 1", "LiTN 2", "LiTN 3"),
      .f = ~ dt %>% filter(site %in% c("Tuck 1", .x)) %>%
        mutate(
          treat      = !(site == "Tuck 1"),
          comparison = paste("Tuck1 v ", .x)
        )
    ) %>%
    bind_rows(dt)
  } %>%
  group_by(comparison, element, layer) %>%
  group_nest()

medf <- function(dt, sims = 250){ 
  model.m <- lm(ll1 ~ treat, data = dt)
  # model.y <- lm(wd1 ~ ll1*treat, data = dt)
  # model.y <- glm(dead ~ ll1*treat, data = dt, family = binomial)
  model.y <- glm(dead ~ ll1 + treat, data = dt)
  med <- mediate(model.m, model.y, treat = "treat", mediator = "ll1", 
                 covariates = "site", 
                 robustSE = TRUE, sims = sims, dropobs = FALSE) 
  summary(med)
}





med_res <- 
  med_dat %>%
  # filter(comparison == "Tuck1 v  LiTN 3") %>%
  mutate(
    med = purrr::map(data, ~ medf(.x, sims = 1000)),
  )

xx <- 
med_res %>%
  mutate(
    res = purrr::map(med, summarise_mediation)
  ) %>%
  dplyr::select(comparison, element, layer, res) %>%
  tidyr::unnest(cols = res)

View(xx)
zz <- 
xx %>%
  filter(d.avg.p < .2 | d0.p < .2 | d1.p < .2) %>% 
  distinct(comparison, element, layer) %>%
  left_join(
    med_dat, by = c("comparison", "element", "layer")
  ) %>%
  mutate(
    med = purrr::map(data, ~ medf(.x, sims = 1000)),
  )

yy <- 
  zz %>%
  mutate(
    res = purrr::map(med, summarise_mediation)
  ) %>%
  dplyr::select(comparison, element, layer, res) %>%
  tidyr::unnest(cols = res)

View(yy)

med_dat2 <-  
  mdat %>%
  mutate(
    # treat = !(site == "Tuck 1"),
    comparison = "rivers"
  ) %>%
  mutate(
    data = purrr::map(
      .x = data,
      .f = ~.x %>%
          dplyr::select(-l1s) %>%
          filter(layer %in% c("ipx", "ncr_new", "ncr_old")) %>%
          mutate(
            treat = (river == "Little Tennessee"),
            l1 = log(l1),
            wd1   = if_else(dead == 1, -baseline_weight, wd)
          ) %>%
          tidyr::pivot_wider(
            names_from = layer,
            values_from = l1
          )
    )
  ) %>%
  tidyr::unnest(cols = data) %>%
  {
    dt <- .
    purrr::map_dfr(
      .x = c("Tuck 2", "Tuck 3", "LiTN 1", "LiTN 2", "LiTN 3"),
      .f = ~ dt %>% filter(site %in% c("Tuck 1", .x)) %>%
        mutate(
          treat      = !(site == "Tuck 1"),
          comparison = paste("Tuck1 v ", .x)
        )
    ) %>%
    bind_rows(dt)
  } %>%
  group_by(comparison, element) %>%
  group_nest()

medf2 <- function(dt, sims = 250){ 
  model.m <- lm(ncr_new ~ treat + ncr_old + ipx, data = dt)
  # model.y <- lm(wd1 ~ ll1*treat, data = dt)
  # model.y <- glm(dead ~ ll1*treat, data = dt, family = binomial)
  model.y <- glm(dead ~ ncr_new + treat, data = dt)
  med <- mediate(model.m, model.y, treat = "treat", mediator = "ncr_new", 
                 covariates = "site", 
                 robustSE = TRUE, sims = sims, dropobs = TRUE) 
  summary(med)
}


med_res <- 
  med_dat2 %>%
  # filter(comparison == "Tuck1 v  LiTN 3") %>%
  mutate(
    med = purrr::map(data, ~ medf2(.x, sims = 1000)),
  )

xx <- 
med_res %>%
  mutate(
    res = purrr::map(med, summarise_mediation)
  ) %>%
  dplyr::select(comparison, element, res) %>%
  tidyr::unnest(cols = res)

View(xx)

ggplot(
  xx %>% filter(comparison == "rivers"),
  aes(x = d.avg, y = -log10(d.avg.p))) +
  geom_point()

bsaul/elktoeChemistry documentation built on Nov. 17, 2022, 8:10 a.m.