In bsaul/tafp: THN Foundation Taste Adaptation Study Pilot
knitr::opts_chunk$set(echo = FALSE)
library(ggplot2)
library(geepack)
# library(lme4)
library(geex)
source("4a_create_analysis_datasets.r")
PRELIMARY RESULTS
Methods
Outcomes
- primary: changes in suprathreshold intensity
- secondary: changes in detection and recognition thresholds
Sample Size and Power
The sample size for this pilot study was based on a model-based simulation study of changes in suprathreshold intensity from day 0 to day 5 of fasting. The target of the simulation was estimating 90\% confidence intervals instead of power. This approach for justifying pilot study sample size is recommended when little is known about the parameter(s) of interest over a hypothesis testing approach [@hertzog2008; @moore2011; @lee2014].
Statistical Analyses
Demographics, baseline characteristics, and fasting lengths were described by appropriate statistical and graphical summaries.
For analysis of suprathreshold intensities, we excluded the zero concentrations as only r supra_analysis %>% filter(conc ==0, response > 0 ) %>% nrow() of the r supra_analysis %>% filter(conc ==0) %>% nrow() tests had a value other than zero at this concentration. The mean and standard deviation of responses were computed for each assay, concentration, and time point. To estimate the marginal effect of days since start of fast on suprathrehold intensity, a Generalized Estimating Equation (GEE) [@liang1986] model of $log(1 + \text{response})$ was fit with terms for days since start of fast, $log_{10}(\text{concentration})$, and their interaction. Models were fit separately within each taste. Model estimates were transformed back to the original scale are presented with point-wise 90\% confidence intervals.
Recognition and detection thresholds were summarized by mean and standard deviations for each assay and time point. TODO: modelling
Statistical analyses were done in r R.version.string [@rcore] using the geepack [@geepack], geex [@geex], and ggplot2 [@ggplot2] packages.
Results
TODO: describe subject characteristics
counts <-
supra_analysis %>%
ungroup() %>%
distinct(record_id, time) %>%
group_by(time) %>%
count()
All subjects completed the assays at the start of fast, day 5, and end of fast + 3 days. Only r counts$n[3] subjects completed the assays at the end of fast timepoint.
fast_length <-
supra_analysis %>%
ungroup() %>%
select(record_id, time, days) %>%
filter(time == 4) %>%
group_by(record_id) %>%
summarise(fast_length = max(days) - 3)
clinical_summary <-
clinical %>%
filter(time == 1) %>%
select(
record_id, height, weight, bmi,
systolic_blood_pressure, diastolic_blood_pressure
) %>%
left_join(
fast_length, by = "record_id"
) %>%
left_join(
demographics, by = "record_id"
) %>%
mutate(
record_id = as.character(record_id)
) %>%
mutate(
male = (gender == 1) * 1
)
fvars <-
clinical_summary %>%
select_if(is.factor) %>%
summarise_all(
.funs = list(tab = ~ list(table(.)))
)
clinical_summary %>%
select_if(is.numeric) %>%
tidyr::gather(key = "key", value = "value") %>%
group_by(key) %>%
summarise(
count = sum(value, na.rm = TRUE),
mean = mean(value, na.rm = TRUE),
median = median(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
iqr = IQR(value, na.rm = TRUE)
) %>%
knitr::kable(
digits = 2
)
clinical_summary %>%
filter(fast_length >= 10) %>%
select_if(is.numeric) %>%
tidyr::gather(key = "key", value = "value") %>%
group_by(key) %>%
summarise(
count = sum(value, na.rm = TRUE),
mean = mean(value, na.rm = TRUE),
median = median(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
iqr = IQR(value, na.rm = TRUE)
) %>%
knitr::kable(
digits = 2,
)
clinical %>%
select(
record_id, time, height, weight, bmi,
sbp = systolic_blood_pressure, dbp = diastolic_blood_pressure
) %>%
tidyr::gather(key = "key", value = "value", -record_id, -time) %>%
group_by(time, key) %>%
summarise(
x = sprintf("%.1f (%.2f)", mean(value, na.rm = TRUE), sd(value, na.rm = TRUE))
) %>%
tidyr::spread(
key = time,
value = x
) %>%
knitr::kable(
align = c("l", rep("c", 4)),
col.names = c("Measure",
paste(c("Start of Fast",
"Day 5", "End of fast", "End of fast + 3 days"),
paste0("(n=", counts$n, ")"))),
caption = "Mean (Standard deviation) of clinical measures over study"
)
counts2 <-
supra_analysis %>%
ungroup() %>%
distinct(record_id, time) %>%
group_by(record_id) %>%
filter(sum(time == 3) > 0) %>%
ungroup() %>%
group_by(time) %>%
count()
clinical %>%
select(
record_id, time, height, weight, bmi,
sbp = systolic_blood_pressure, dbp = diastolic_blood_pressure
) %>%
group_by(record_id) %>%
filter(sum(time == 3) > 0) %>%
filter(!is.na(height[time == 3])) %>%
tidyr::gather(key = "key", value = "value", -record_id, -time) %>%
group_by(time, key) %>%
summarise(
x = sprintf("%.1f (%.2f)", mean(value, na.rm = TRUE), sd(value, na.rm = TRUE))
) %>%
tidyr::spread(
key = time,
value = x
) %>%
knitr::kable(
align = c("l", rep("c", 4)),
col.names = c("Measure",
paste(c("Start of Fast",
"Day 5", "End of fast", "End of fast + 3 days"),
paste0("(n=", counts2$n, ")"))),
caption = "Mean (Standard deviation) of clinical measures over study among those with end of fast record"
)
Changes in Supratheshold Intensity
supra_analysis %>%
ungroup %>%
# filter(conc > 0) %>%
group_by(assay_taste, conc, time) %>%
summarise(
x = sprintf("%.1f (%.2f)", mean(response), sd(response))
) %>%
tidyr::spread(
key = time,
value = x
) %>%
ungroup() %>%
mutate(
assay_taste = factor(
assay_taste,
levels = c("nacl", "sucr"),
labels = c("NaCL", "Sucrose")),
) %>%
knitr::kable(
align = c("l", "l", rep("c", 4)),
col.names = c("Assay", "Concentration",
paste(c("Start of Fast",
"Day 5", "End of fast", "End of fast + 3 days"),
paste0("(n=", counts$n, ")"))),
caption = "Mean (Standard deviation) of Suprathreshold repsonses at timepoints"
)
supra_model_dt <-
supra_analysis %>%
ungroup %>%
filter(conc > 0) %>%
mutate(
y = log1p(response),
) %>%
select(y, days, assay_taste, time, c = conc, id = record_id) %>%
arrange(id) %>%
group_nest(assay_taste) %>%
mutate(
data_pre_to_5 = purrr::map(data, ~ .x %>% filter(time %in% c(1, 2)) )
)
geex_correct <- function(data, model){
f <- geex::grab_psiFUN(object = model, data = data)
function(theta){
f(theta)
}
}
gee_fit1 <- function(data){
geepack::geeglm(
formula = y ~ days * log10(c),
id = id,
data = data,
family = gaussian(link = 'identity'))
}
mest <- function(fitter, data){
m <- fitter(data)
m_estimate(
data = data,
units = 'id',
estFUN = geex_correct,
outer_args = list(model = m),
corrections = list(fay = correction(
FUN = fay_bias_correction,
b = 0.75)),
compute_roots = FALSE,
roots = coef(m))
}
get_est <- function(fit, ndays = 17, conc = c(50, 100, 200)){
entries <-
expand.grid(
days = 0:ndays,
conc = conc
)
# Contrast matrix
# TODO: This is hardcoded based on the model specification!!
make_L <- function(d, c){
c(1, d, log10(c), d * log10(c))
}
L <- t(apply(entries, 1, function(x) { make_L(x[1], x[2]) }))
cbind(
entries,
est = L %*% coef(fit),
se = sqrt(diag(L %*% vcov(fit) %*% t(L)))
)
}
get_est_diff <- function(fit, ndays = 10, conc = c(50, 100, 200)){
# Contrast matrix
# TODO: This is hardcoded based on the model specification!!
make_L <- function(d, c){
c(1, d, log10(c), d * log10(c))
}
# Use delta method to obtain standard errors:
# h(theta) = expm1(X_10 %*% theta) - expm1(X_0 %*% theta)
# h'(theta) =
# (
# exp(X_10 %*% theta) - exp(X_0 %*% theta)
# 10 * exp(X_10 %*% theta)
# conc * (exp(X_10 %*% theta) - exp(X_0 %*% theta))
# 10 * conc * exp(X_10 %*% theta)
# )
theta <- coef(fit)
sigma <- vcov(fit)
est <- numeric(length(conc))
se <- numeric(length(conc))
# Obtaining for each concentration independently; though this could be done
# all at once
for(i in seq_along(conc)){
Ln <- make_L(ndays, conc[i])
L0 <- make_L(0, conc[i])
J <- (Ln %*% exp(Ln %*% theta)) - (L0 %*% exp(L0 %*% theta))
V <- t(J) %*% vcov(fit) %*% J
est[i] <- expm1(Ln %*% theta) - expm1(L0 %*% theta)
se[i] <- sqrt(V)
}
tibble::tibble(
conc,
est = est,
se = se
)
}
supra_model_dt %>%
mutate(
fit_alldays = purrr::map(
.x = data,
.f = ~ mest(gee_fit1, .x)
),
est_alldays = purrr::map2(
.x = fit_alldays,
.y = assay_taste,
.f = ~ {
if(.y == "sucr") {
get_est(.x, conc = c(100, 200, 400))
} else {
get_est(.x)
} }
),
est_diff_day10 = purrr::map2(
.x = fit_alldays,
.y = assay_taste,
.f = ~ {
if(.y == "sucr") {
get_est_diff(.x, ndays = 10, conc = c(100, 200, 400))
} else {
get_est_diff(.x, ndays = 10)
} })
) ->
supra_fits
supra_fits %>%
select(
assay_taste, est_alldays
) %>%
tidyr::unnest(cols = c(est_alldays)) %>%
mutate(
conf_lo = est - qnorm(.9) * se,
conf_hi = est + qnorm(.9) * se,
est = expm1(est),
conf_lo = expm1(conf_lo),
conf_hi = expm1(conf_hi)
) ->
supra_model_res
supra_changes <-
supra_fits %>%
select(
assay_taste, est_diff_day10
) %>%
tidyr::unnest(cols = c(est_diff_day10)) %>%
mutate(
conf_lo = est - qnorm(.9) * se,
conf_hi = est + qnorm(.9) * se,
est = est
) %>%
mutate(label = sprintf("%.1f (%.1f, %.1f)", est, conf_lo, conf_hi))
supra_nacl_changes <- filter(supra_changes, assay_taste == "nacl") %>% pull(label)
supra_sucr_changes <- filter(supra_changes, assay_taste == "sucr") %>% pull(label)
Figure \@ref(fig:supra_days) displays changes in suprathreshold intensity for salt and sweet taste during the study period for all subjects along with the estimated GEE model. Responses are highly variable, especially between subjects (Table \@ref(table1)), but show a general trend of increasing response to salt and descreasing response to sweet. For example, the GEE model estimates that after 10 days of fasting the average salt taste response increases by r supra_nacl_changes %>% glue::glue_collapse(sep = ", ", last = ", and ") for the 50, 100, and 200 mM concentrations, respectively. Over the same time span, the average sweet taste response is estimated to change by r supra_sucr_changes %>% glue::glue_collapse(sep = ", ", last = ", and ") for each concentration respectively.
supra_plot <- supra_analysis %>%
filter(conc != 0) %>%
ungroup() %>%
mutate(
assay_taste = factor(
assay_taste,
levels = c("nacl", "sucr"),
labels = c("NaCL", "Sucrose")),
conc_label = factor(
conc,
levels = sort(unique(conc)),
labels = sprintf("%s mM", sort(unique(conc))),
ordered = TRUE)
)
conc_labs <- distinct(supra_plot, assay_taste, conc, conc_rank, conc_label)
supra_model_plot <-
supra_model_res %>%
mutate(
assay_taste = factor(
assay_taste,
levels = c("nacl", "sucr"),
labels = c("NaCL", "Sucrose"))
) %>%
left_join(
conc_labs,
by = c("assay_taste", "conc")
)
supra_changes_plot <-
supra_changes %>%
mutate(
assay_taste = factor(
assay_taste,
levels = c("nacl", "sucr"),
labels = c("NaCL", "Sucrose"))
) %>%
left_join(
conc_labs,
by = c("assay_taste", "conc")
) %>%
left_join(
y = supra_model_plot %>% filter(days == 10) %>%
select(assay_taste, conc, ypos = est),
by = c("assay_taste", "conc")
) %>%
left_join(
y = supra_model_plot %>% filter(days == 0) %>%
select(assay_taste, conc, y0pos = est),
by = c("assay_taste", "conc")
)
ggplot(
data = supra_plot,
mapping = aes(x = days, y = response)
) +
geom_hline(
yintercept = 0, color = "grey50", size = 0.5
) +
geom_line(
aes(group = record_id),
color = "grey80",
size = .3) +
geom_point(
color = "grey80",
size = 0.3
) +
geom_line(
data = supra_model_plot,
aes(x = days, y = est, color = assay_taste),
inherit.aes = FALSE
) +
geom_ribbon(
data = supra_model_plot,
aes(x = days, ymin = conf_lo, ymax = conf_hi, fill = assay_taste),
alpha = 0.25,
inherit.aes = FALSE
) +
geom_segment(
data = supra_changes_plot,
aes(x = 0, xend = 10, y = y0pos, yend = y0pos),
linetype = "dotted",
size = .2, color = "grey20"
) +
geom_segment(
data = supra_changes_plot,
aes(x = 10, xend = 10, y = ypos, yend = y0pos),
linetype = "dotted",
size = .2, color = "grey20"
) +
geom_text(
data = supra_changes_plot,
aes(x = 10, y = ypos, label = label),
hjust = 0, vjust = 0, size = 2, color = "grey10"
) +
geom_text(
data = conc_labs,
aes(x = 2, y = 5, label = conc_label),
hjust = 0, size = 2.5, color = "grey40"
) +
facet_grid(assay_taste ~ conc_rank) +
scale_y_continuous(
expand = c(0, 0)
) +
scale_color_manual(
guide = FALSE,
values = c("#7b4244", "#60e9dc")
) +
scale_fill_manual(
guide = FALSE,
values = c("#7b4244", "#60e9dc")
) +
labs(
title = "Changes in suprathreshold intensity",
x = "Days since start of fast",
y = "Response",
caption = "Each line is a single subject's trajectory.\nColored lines show GEE model fit with pointwise 90% confidence intervals.\nValue is estimated average change from Day 0 to 10 with 90% confidence interval."
) +
theme_classic() +
theme(
strip.background = element_blank(),
# strip.background.y = element_blank(),
strip.text.x = element_blank(),
strip.text.y = element_text(angle = 0),
axis.line.x = element_blank(),
axis.line.y = element_line(color = "grey50", size = 0.5),
axis.ticks.y = element_line(color = "grey50", size = 0.1),
axis.ticks.x = element_line(color = "grey50", size = 0.1),
axis.text = element_text(size = 6),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
plot.caption = element_text(size = 6)
)
Changes in Detection/Recognition Thresholds
dtrt_dt <- dtrt_thresholds %>%
ungroup() %>%
select(-dtp, -rtp, -taste_order) %>%
tidyr::gather(key = "type", value = "conc",
-record_id, -time, -days, -assay_taste) %>%
mutate(
assay_taste = factor(
assay_taste,
levels = c("nacl", "sucr"),
labels = c("NaCL", "Sucrose")),
type = factor(
type,
levels = c("dt", "rt"),
labels = c("Detection", "Recognition")
)
) %>%
group_by(assay_taste, type, record_id) %>%
arrange(time, .by_group = TRUE) %>%
mutate(
conc2 = log1p(conc) - lag(log1p(conc), default = NA)
)
dtrt_dt %>%
group_by(assay_taste, time, type) %>%
summarise(
x = sprintf("%.1f (%.2f)", mean(conc), sd(conc))
) %>%
tidyr::spread(
key = time,
value = x
) %>%
knitr::kable(
align = c("l", "l", rep("c", 4)),
col.names = c("Assay", "Threshold",
"Start of Fast",
"Day 5", "End of fast", "End of fast + 3 days"),
caption = "Mean (Standard deviation) of thresholds at study timepoints"
)
A simple paired t-test between day 0 and end of fast:
dtrt_dt %>%
ungroup() %>%
filter(time %in% c(1, 3)) %>%
select(-days, -conc2) %>%
group_by(assay_taste, type) %>%
tidyr::pivot_wider(
names_from = time,
names_prefix = "time",
values_from = conc
) %>%
tidyr::nest() %>%
mutate(
p_value = purrr::map_dbl(
.x = data,
.f = ~ t.test(x = .x[["time1"]], y = .x[["time3"]], paired = TRUE)[["p.value"]])
) %>%
select(
assay_taste, type, p_value
)
dtrt_model_dt <-
dtrt_dt %>%
ungroup %>%
mutate(
y = log1p(conc),
) %>%
select(type, y, days, assay_taste, time, c = conc, id = record_id) %>%
arrange(id) %>%
group_nest(type, assay_taste) %>%
mutate(
data_pre_to_5 = purrr::map(data, ~ .x %>% filter(time %in% c(1, 2)) )
)
geex_correct <- function(data, model){
f <- geex::grab_psiFUN(object = model, data = data)
L <- matrix(
c(rep(1, 20), 0:19), ncol = 2
)
function(theta){
# browser()
m <- f(theta[1:2])
c(m, # GEE model
expm1(L %*% theta[1:2]) - theta[-c(1:2, 23)], # estimate for days 0:19
theta[23] - (theta[13] - theta[3]) # difference between day 10 and day 0
)
}
}
gee_fit1 <- function(data){
geepack::geeglm(
formula = y ~ days,
id = id,
data = data,
family = gaussian(link = 'identity'))
}
mest <- function(fitter, data){
m <- fitter(data)
m_estimate(
data = data,
units = 'id',
estFUN = geex_correct,
outer_args = list(model = m),
corrections = list(fay = correction(
FUN = fay_bias_correction,
b = 0.75)),
compute_roots = TRUE,
root_control = setup_root_control(start = c(coef(m), rep(3, 20), 0))
# roots = c(coef(m), rep(3, 20))
)
}
# get_est <- function(fit, ndays = 17){
#
# entries <-
# expand.grid(
# days = 0:ndays
# )
#
# # Contrast matrix
# # TODO: This is hardcoded based on the model specification!!
# make_L <- function(d, c){
# c(1, d)
# }
#
# L <- t(apply(entries, 1, function(x) { make_L(x[1], x[2]) }))
#
# cbind(
# entries,
# est = L %*% coef(fit),
# se = sqrt(diag(L %*% vcov(fit) %*% t(L)))
# )
# }
get_est <- function(fit){
tibble::tibble(
days = c(0:19, -10),
est = fit@estimates[-c(1,2)],
se = sqrt(diag(vcov(fit)))[-c(1,2)]
)
}
dtrt_model_dt %>%
mutate(
fit_alldays = purrr::map(
.x = data,
.f = ~ mest(gee_fit1, .x)
),
est_alldays = purrr::map(
.x = fit_alldays,
.f = ~ get_est(.x)
)
) ->
dtrt_fits
dtrt_fits %>%
select(
type, assay_taste, est_alldays
) %>%
tidyr::unnest(cols = c(est_alldays)) %>%
mutate(
conf_lo = est - qnorm(.9) * se,
conf_hi = est + qnorm(.9) * se,
est = est,
conf_lo = conf_lo,
conf_hi = conf_hi
) ->
dtrt_model_res
dtrt_changes <-
dtrt_model_res %>%
filter(days == -10) %>%
group_by(type, assay_taste) %>%
mutate(label = sprintf("%.1f (%.1f, %.1f)", est, conf_lo, conf_hi))
# dtrt_changes <-
# dtrt_model_res %>%
# filter(days %in% c(0, 10)) %>%
# group_by(type, assay_taste) %>%
# summarise(
# diff = expm1(est[days == 10]) - expm1(est[days == 0])
# )
# dtrt_changes_orig_scale <- dtrt_model_res %>%
# filter(days == 10) %>%
# mutate(
# est = expm1(est)
# )
dt_nacl_changes <- filter(
dtrt_changes,
assay_taste == "NaCL",
type == "Detection") %>% pull(label)
rt_nacl_changes <- filter(
dtrt_changes,
assay_taste == "NaCL",
type == "Recognition") %>% pull(label)
dt_sucr_changes <- filter(
dtrt_changes,
assay_taste == "Sucrose",
type == "Detection") %>% pull(label)
rt_sucr_changes <- filter(
dtrt_changes,
assay_taste == "Sucrose",
type == "Recognition") %>% pull(label)
dtrt_changes %>%
select(type, assay_taste, label) %>%
knitr::kable(
caption = "Estimated change in threshold from day 0 to day 10 (90% confidence interval)."
)
dtrt_plot <- dtrt_model_dt %>%
select(type, assay_taste, data) %>%
tidyr::unnest(cols = data) %>%
filter(days >= 0)
dtrt_model_plot <- dtrt_model_res %>%
filter(days >= 0 )
dtrt_changes_plot <-
dtrt_changes %>%
left_join(
y = dtrt_model_plot %>% filter(days == 10) %>%
select(type, assay_taste, ypos = est),
by = c("type", "assay_taste")
) %>%
left_join(
y = dtrt_model_plot %>% filter(days == 0) %>%
select(type, assay_taste, y0pos = est),
by = c("type", "assay_taste")
) %>%
mutate(days = 10)
Figure \@ref(fig:dtrt_days) displays changes in suprathreshold intensity for salt and sweet taste during the study period for all subjects along with the estimated GEE model. Responses are highly variable, especially between subjects (Table \@ref(table1)), but show a general trend of decreasing recognition response for salt and sweet response.
ggplot(
data = dtrt_plot,
mapping = aes(x = days, y = expm1(y))
) +
geom_hline(
yintercept = 0, color = "grey50", size = 0.5
) +
geom_line(
aes(group = id),
color = "grey80",
size = .3) +
geom_point(
color = "grey80",
size = 0.3
) +
geom_line(
data = dtrt_model_plot,
aes(x = days, y = est, color = assay_taste),
inherit.aes = FALSE
) +
geom_ribbon(
data = dtrt_model_plot,
aes(x = days, ymin = conf_lo, ymax = conf_hi, fill = assay_taste),
alpha = 0.25,
inherit.aes = FALSE
) +
geom_segment(
data = dtrt_changes_plot,
aes(x = 0, xend = 10, y = y0pos, yend = y0pos),
linetype = "dotted",
size = .2, color = "grey20"
) +
geom_segment(
data = dtrt_changes_plot,
aes(x = 10, xend = 10, y = ypos, yend = y0pos),
linetype = "dotted",
size = .2, color = "grey20"
) +
geom_text(
data = dtrt_changes_plot,
aes(x = 10, y = ypos, label = label),
hjust = 0, vjust = 0, size = 2, color = "grey10"
) +
facet_grid(assay_taste ~ type, scales = "free") +
scale_y_continuous(
expand = c(0, 0)
) +
scale_color_manual(
guide = FALSE,
values = c("#7b4244", "#60e9dc")
) +
scale_fill_manual(
guide = FALSE,
values = c("#7b4244", "#60e9dc")
) +
labs(
title = "Changes in detection/recognition thresholds",
x = "Days since start of fast",
y = "Threshold (mM)",
caption = "Each line is a single subject's trajectory.\nColored lines show GEE model fit with pointwise 90% confidence intervals.\nValue is estimated average change from Day 0 to 10 with 90% confidence interval."
) +
theme_classic() +
theme(
strip.background = element_blank(),
# strip.background.y = element_blank(),
# strip.text.x = element_blank(),
strip.text.y = element_text(angle = 0),
axis.line.x = element_blank(),
axis.line.y = element_line(color = "grey50", size = 0.5),
axis.ticks.y = element_line(color = "grey50", size = 0.1),
axis.ticks.x = element_line(color = "grey50", size = 0.1),
axis.text = element_text(size = 6),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
plot.caption = element_text(size = 6)
)
dtrt_dt %>%
filter(time %in% c(1, 2)) %>%
ggplot(
aes(x = factor(days), y = conc, group = record_id)
) +
geom_point(size = 0.2) +
geom_line(size = 0.2) +
facet_grid(
assay_taste ~ type, scales = "free_y"
)
baseline_vars <-
clinical_summary %>%
select(record_id, bmi, age, male)
hold_thresholds <-
dtrt_thresholds %>%
filter(time %in% c(1, 4)) %>%
ungroup() %>%
select(-taste_order, -days, -dtp, -rtp) %>%
tidyr::pivot_longer(cols = c(dt, rt)) %>%
group_by(record_id, assay_taste, name) %>%
summarise(
value = value[time == 4] - value[time == 1]
) %>%
mutate(
time = 41
) %>%
bind_rows(
select(ungroup(dtrt_thresholds), -taste_order, -days, -dtp, -rtp) %>%
tidyr::pivot_longer(cols = c(dt, rt))
)
hold_for_cors <-
hold_thresholds %>%
filter(time %in% c(1, 41)) %>%
left_join(
baseline_vars, by = "record_id"
)
hold_for_cors %>%
tidyr::pivot_longer(
cols = c(age, male, bmi),
names_to = "baseline_var",
values_to = "baseline_value"
) %>%
ungroup() %>%
group_nest(
time, assay_taste, name, baseline_var
) %>%
mutate(
cort = purrr::map(
.x = data,
.f = ~ cor.test(.x$value, .x$baseline_value)),
cor_est = purrr::map_dbl(cort, "estimate"),
cor_pval = purrr::map_dbl(cort, "p.value"),
value = sprintf("%.3f (%.2f)", cor_est, cor_pval)
) %>%
select(-data, -cort, -cor_est, -cor_pval) %>%
tidyr::pivot_wider(
id_cols = c(time, assay_taste, name),
names_from = baseline_var
) %>%
mutate(
time = if_else(time == 1, "Baseline",
"Change from baseline to re-feed")
) %>%
knitr::kable(
caption = "Correlations between baseline covariates and dectection and recognition thresholds."
)
Changes in food-liking
flq_completed <-
foodliking %>%
filter(complete.cases(.)) %>%
group_by(record_id) %>%
filter(n() == 2) %>%
ungroup() %>%
mutate(record_id = as.character(record_id))
There were r length(unique(flq_completed$record_id)) subjects who completed the food-liking questionnaire both before and after the exposure time.
TODO: describe the scoring
Figure XX shows the changes in food-liking.
flq_plot <-
flq_completed %>%
tidyr::pivot_longer(
cols = ends_with("flq")
) %>%
mutate(
name = stringr::str_remove(name, "_flq"),
name = factor(name, levels = c("salt", "sweet", "fatsalt", "fatsweet"),
labels = c("Salt", "Sweet", "Fat & Salt", "Fat & Sweet"))
)
hold <-
flq_plot %>%
tidyr::pivot_wider(
id_cols = c("record_id", "name"),
names_from = "time"
) %>%
ungroup() %>%
group_nest(name) %>%
mutate(
ttest = purrr::map(
.x = data,
.f = ~ t.test(.x$`5`, .x$`0`, paired = TRUE)
),
tpval = purrr::map_dbl(ttest, "p.value"),
ttext = sprintf("Paired t-test p-value: %.3f", tpval)
)
hold %>%
mutate(
estimate = purrr::map(ttest, ~ .x$estimate)
) %>%
select(name, estimate) %>%
knitr::kable(caption = "Paired Changes in FLQ")
flq_plot %>%
ggplot(
data = .,
aes(x = time, y = value, color = name)
) +
geom_hline(
yintercept = 0,
color = "grey50"
) +
geom_point(
alpha = 0.5
) +
geom_line(
aes(group = record_id),
size = 0.1,
alpha = 0.5
) +
geom_text(
data = hold,
aes(label = ttext, x = 2, y = 1),
color = "black",
size = 2
) +
stat_smooth(
# aes(color = assay_taste),
method = "lm", se = FALSE
) +
scale_y_continuous(
"FLQ Score"
) +
scale_x_continuous(
"Time",
breaks = c(0, 5),
labels = c("Pre", "Post")
) +
scale_color_manual(
"",
values = c("#8c75d2", "#11677e", "#d547ca", "#3e3c8d"),
guide = FALSE
) +
facet_wrap(
~ name,
ncol = 2
) +
theme_classic() +
theme(
strip.background = element_blank(),
# strip.background.y = element_blank(),
# strip.text.x = element_blank(),
strip.text.y = element_text(angle = 0),
axis.line.x = element_blank(),
axis.line.y = element_line(color = "grey50", size = 0.5),
axis.ticks.y = element_line(color = "grey50", size = 0.1),
axis.ticks.x = element_line(color = "grey50", size = 0.1),
axis.text = element_text(size = 6),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
plot.caption = element_text(size = 6)
)
Food-liking and Adaptation
dtrt_hold <-
dtrt_thresholds %>%
ungroup() %>%
filter(time %in% c(1, 4), record_id %in% flq_completed$record_id) %>%
select(record_id, time, assay_taste, dt, rt) %>%
tidyr::pivot_longer(
cols = c(dt, rt),
names_to = "threshold"
) %>%
group_by(record_id, assay_taste, threshold) %>%
summarise(
dtrt_del = value[time ==4] - value[time == 1]
) %>%
ungroup() %>%
mutate(
threshold = factor(threshold, levels = c("dt", "rt"),
labels = c("Detection", "Recognition")),
assay_taste = factor(assay_taste, levels = c("nacl", "sucr"),
labels = c("NaCL", "Sucrose"))
)
flq_hold <-
flq_completed %>%
ungroup() %>%
filter(time %in% c(0, 5)) %>%
tidyr::pivot_longer(
cols = ends_with("_flq"),
names_to = "flq_name"
) %>%
group_by(record_id, flq_name) %>%
summarise(
flq_del = value[time == 5] - value[time == 0]
) %>%
mutate(
flq_name = stringr::str_remove(flq_name, "_flq"),
flq_name = factor(flq_name, levels = c("salt", "sweet", "fatsalt", "fatsweet"),
labels = c("Salt", "Sweet", "Fat & Salt", "Fat & Sweet"))
)
plot_hold <-
full_join(
dtrt_hold,
flq_hold,
by = "record_id"
)
ggplot(
data = plot_hold,
aes(x = dtrt_del, y = flq_del)
) +
stat_smooth(
method = "lm",
size = 0.5,
se = TRUE
) +
geom_point(
aes(color = flq_name),
alpha = 0.5
) +
scale_color_manual(
"",
values = c("#8c75d2", "#11677e", "#d547ca", "#3e3c8d"),
guide = FALSE
) +
scale_y_continuous(
"Change in food liking score"
) +
scale_x_continuous(
"Change in threshold"
) +
facet_grid(
flq_name ~ threshold + assay_taste,
scales = "free"
) +
theme_classic() +
theme(
panel.background = element_rect(fill = "grey90"),
panel.grid.major = element_line(color = "white"),
strip.background = element_blank(),
# strip.background.y = element_blank(),
# strip.text.x = element_blank(),
strip.text.y = element_text(angle = 0),
# axis.line.x = element_blank(),
axis.line.x = element_line(color = "grey50", size = 0.5),
axis.line.y = element_line(color = "grey50", size = 0.5),
axis.ticks.y = element_line(color = "grey50", size = 0.1),
axis.ticks.x = element_line(color = "grey50", size = 0.1),
axis.text = element_text(size = 6),
axis.title.y = element_text(size = 8),
axis.title.x = element_text(size = 8),
plot.caption = element_text(size = 6)
)
References
bsaul/tafp documentation built on Jan. 28, 2022, 10:16 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = FALSE) library(ggplot2) library(geepack) # library(lme4) library(geex) source("4a_create_analysis_datasets.r")
PRELIMARY RESULTS
Methods
Outcomes
- primary: changes in suprathreshold intensity
- secondary: changes in detection and recognition thresholds
Sample Size and Power
The sample size for this pilot study was based on a model-based simulation study of changes in suprathreshold intensity from day 0 to day 5 of fasting. The target of the simulation was estimating 90\% confidence intervals instead of power. This approach for justifying pilot study sample size is recommended when little is known about the parameter(s) of interest over a hypothesis testing approach [@hertzog2008; @moore2011; @lee2014].
Statistical Analyses
Demographics, baseline characteristics, and fasting lengths were described by appropriate statistical and graphical summaries.
For analysis of suprathreshold intensities, we excluded the zero concentrations as only r supra_analysis %>% filter(conc ==0, response > 0 ) %>% nrow() of the r supra_analysis %>% filter(conc ==0) %>% nrow() tests had a value other than zero at this concentration. The mean and standard deviation of responses were computed for each assay, concentration, and time point. To estimate the marginal effect of days since start of fast on suprathrehold intensity, a Generalized Estimating Equation (GEE) [@liang1986] model of $log(1 + \text{response})$ was fit with terms for days since start of fast, $log_{10}(\text{concentration})$, and their interaction. Models were fit separately within each taste. Model estimates were transformed back to the original scale are presented with point-wise 90\% confidence intervals.
Recognition and detection thresholds were summarized by mean and standard deviations for each assay and time point. TODO: modelling
Statistical analyses were done in r R.version.string [@rcore] using the geepack [@geepack], geex [@geex], and ggplot2 [@ggplot2] packages.
Results
TODO: describe subject characteristics
counts <- supra_analysis %>% ungroup() %>% distinct(record_id, time) %>% group_by(time) %>% count()
All subjects completed the assays at the start of fast, day 5, and end of fast + 3 days. Only r counts$n[3] subjects completed the assays at the end of fast timepoint.
fast_length <- supra_analysis %>% ungroup() %>% select(record_id, time, days) %>% filter(time == 4) %>% group_by(record_id) %>% summarise(fast_length = max(days) - 3) clinical_summary <- clinical %>% filter(time == 1) %>% select( record_id, height, weight, bmi, systolic_blood_pressure, diastolic_blood_pressure ) %>% left_join( fast_length, by = "record_id" ) %>% left_join( demographics, by = "record_id" ) %>% mutate( record_id = as.character(record_id) ) %>% mutate( male = (gender == 1) * 1 ) fvars <- clinical_summary %>% select_if(is.factor) %>% summarise_all( .funs = list(tab = ~ list(table(.))) ) clinical_summary %>% select_if(is.numeric) %>% tidyr::gather(key = "key", value = "value") %>% group_by(key) %>% summarise( count = sum(value, na.rm = TRUE), mean = mean(value, na.rm = TRUE), median = median(value, na.rm = TRUE), sd = sd(value, na.rm = TRUE), iqr = IQR(value, na.rm = TRUE) ) %>% knitr::kable( digits = 2 )
clinical_summary %>% filter(fast_length >= 10) %>% select_if(is.numeric) %>% tidyr::gather(key = "key", value = "value") %>% group_by(key) %>% summarise( count = sum(value, na.rm = TRUE), mean = mean(value, na.rm = TRUE), median = median(value, na.rm = TRUE), sd = sd(value, na.rm = TRUE), iqr = IQR(value, na.rm = TRUE) ) %>% knitr::kable( digits = 2, )
clinical %>% select( record_id, time, height, weight, bmi, sbp = systolic_blood_pressure, dbp = diastolic_blood_pressure ) %>% tidyr::gather(key = "key", value = "value", -record_id, -time) %>% group_by(time, key) %>% summarise( x = sprintf("%.1f (%.2f)", mean(value, na.rm = TRUE), sd(value, na.rm = TRUE)) ) %>% tidyr::spread( key = time, value = x ) %>% knitr::kable( align = c("l", rep("c", 4)), col.names = c("Measure", paste(c("Start of Fast", "Day 5", "End of fast", "End of fast + 3 days"), paste0("(n=", counts$n, ")"))), caption = "Mean (Standard deviation) of clinical measures over study" )
counts2 <- supra_analysis %>% ungroup() %>% distinct(record_id, time) %>% group_by(record_id) %>% filter(sum(time == 3) > 0) %>% ungroup() %>% group_by(time) %>% count() clinical %>% select( record_id, time, height, weight, bmi, sbp = systolic_blood_pressure, dbp = diastolic_blood_pressure ) %>% group_by(record_id) %>% filter(sum(time == 3) > 0) %>% filter(!is.na(height[time == 3])) %>% tidyr::gather(key = "key", value = "value", -record_id, -time) %>% group_by(time, key) %>% summarise( x = sprintf("%.1f (%.2f)", mean(value, na.rm = TRUE), sd(value, na.rm = TRUE)) ) %>% tidyr::spread( key = time, value = x ) %>% knitr::kable( align = c("l", rep("c", 4)), col.names = c("Measure", paste(c("Start of Fast", "Day 5", "End of fast", "End of fast + 3 days"), paste0("(n=", counts2$n, ")"))), caption = "Mean (Standard deviation) of clinical measures over study among those with end of fast record" )
Changes in Supratheshold Intensity
supra_analysis %>% ungroup %>% # filter(conc > 0) %>% group_by(assay_taste, conc, time) %>% summarise( x = sprintf("%.1f (%.2f)", mean(response), sd(response)) ) %>% tidyr::spread( key = time, value = x ) %>% ungroup() %>% mutate( assay_taste = factor( assay_taste, levels = c("nacl", "sucr"), labels = c("NaCL", "Sucrose")), ) %>% knitr::kable( align = c("l", "l", rep("c", 4)), col.names = c("Assay", "Concentration", paste(c("Start of Fast", "Day 5", "End of fast", "End of fast + 3 days"), paste0("(n=", counts$n, ")"))), caption = "Mean (Standard deviation) of Suprathreshold repsonses at timepoints" )
supra_model_dt <- supra_analysis %>% ungroup %>% filter(conc > 0) %>% mutate( y = log1p(response), ) %>% select(y, days, assay_taste, time, c = conc, id = record_id) %>% arrange(id) %>% group_nest(assay_taste) %>% mutate( data_pre_to_5 = purrr::map(data, ~ .x %>% filter(time %in% c(1, 2)) ) )
geex_correct <- function(data, model){ f <- geex::grab_psiFUN(object = model, data = data) function(theta){ f(theta) } } gee_fit1 <- function(data){ geepack::geeglm( formula = y ~ days * log10(c), id = id, data = data, family = gaussian(link = 'identity')) } mest <- function(fitter, data){ m <- fitter(data) m_estimate( data = data, units = 'id', estFUN = geex_correct, outer_args = list(model = m), corrections = list(fay = correction( FUN = fay_bias_correction, b = 0.75)), compute_roots = FALSE, roots = coef(m)) } get_est <- function(fit, ndays = 17, conc = c(50, 100, 200)){ entries <- expand.grid( days = 0:ndays, conc = conc ) # Contrast matrix # TODO: This is hardcoded based on the model specification!! make_L <- function(d, c){ c(1, d, log10(c), d * log10(c)) } L <- t(apply(entries, 1, function(x) { make_L(x[1], x[2]) })) cbind( entries, est = L %*% coef(fit), se = sqrt(diag(L %*% vcov(fit) %*% t(L))) ) } get_est_diff <- function(fit, ndays = 10, conc = c(50, 100, 200)){ # Contrast matrix # TODO: This is hardcoded based on the model specification!! make_L <- function(d, c){ c(1, d, log10(c), d * log10(c)) } # Use delta method to obtain standard errors: # h(theta) = expm1(X_10 %*% theta) - expm1(X_0 %*% theta) # h'(theta) = # ( # exp(X_10 %*% theta) - exp(X_0 %*% theta) # 10 * exp(X_10 %*% theta) # conc * (exp(X_10 %*% theta) - exp(X_0 %*% theta)) # 10 * conc * exp(X_10 %*% theta) # ) theta <- coef(fit) sigma <- vcov(fit) est <- numeric(length(conc)) se <- numeric(length(conc)) # Obtaining for each concentration independently; though this could be done # all at once for(i in seq_along(conc)){ Ln <- make_L(ndays, conc[i]) L0 <- make_L(0, conc[i]) J <- (Ln %*% exp(Ln %*% theta)) - (L0 %*% exp(L0 %*% theta)) V <- t(J) %*% vcov(fit) %*% J est[i] <- expm1(Ln %*% theta) - expm1(L0 %*% theta) se[i] <- sqrt(V) } tibble::tibble( conc, est = est, se = se ) } supra_model_dt %>% mutate( fit_alldays = purrr::map( .x = data, .f = ~ mest(gee_fit1, .x) ), est_alldays = purrr::map2( .x = fit_alldays, .y = assay_taste, .f = ~ { if(.y == "sucr") { get_est(.x, conc = c(100, 200, 400)) } else { get_est(.x) } } ), est_diff_day10 = purrr::map2( .x = fit_alldays, .y = assay_taste, .f = ~ { if(.y == "sucr") { get_est_diff(.x, ndays = 10, conc = c(100, 200, 400)) } else { get_est_diff(.x, ndays = 10) } }) ) -> supra_fits supra_fits %>% select( assay_taste, est_alldays ) %>% tidyr::unnest(cols = c(est_alldays)) %>% mutate( conf_lo = est - qnorm(.9) * se, conf_hi = est + qnorm(.9) * se, est = expm1(est), conf_lo = expm1(conf_lo), conf_hi = expm1(conf_hi) ) -> supra_model_res supra_changes <- supra_fits %>% select( assay_taste, est_diff_day10 ) %>% tidyr::unnest(cols = c(est_diff_day10)) %>% mutate( conf_lo = est - qnorm(.9) * se, conf_hi = est + qnorm(.9) * se, est = est ) %>% mutate(label = sprintf("%.1f (%.1f, %.1f)", est, conf_lo, conf_hi)) supra_nacl_changes <- filter(supra_changes, assay_taste == "nacl") %>% pull(label) supra_sucr_changes <- filter(supra_changes, assay_taste == "sucr") %>% pull(label)
Figure \@ref(fig:supra_days) displays changes in suprathreshold intensity for salt and sweet taste during the study period for all subjects along with the estimated GEE model. Responses are highly variable, especially between subjects (Table \@ref(table1)), but show a general trend of increasing response to salt and descreasing response to sweet. For example, the GEE model estimates that after 10 days of fasting the average salt taste response increases by r supra_nacl_changes %>% glue::glue_collapse(sep = ", ", last = ", and ") for the 50, 100, and 200 mM concentrations, respectively. Over the same time span, the average sweet taste response is estimated to change by r supra_sucr_changes %>% glue::glue_collapse(sep = ", ", last = ", and ") for each concentration respectively.
supra_plot <- supra_analysis %>% filter(conc != 0) %>% ungroup() %>% mutate( assay_taste = factor( assay_taste, levels = c("nacl", "sucr"), labels = c("NaCL", "Sucrose")), conc_label = factor( conc, levels = sort(unique(conc)), labels = sprintf("%s mM", sort(unique(conc))), ordered = TRUE) ) conc_labs <- distinct(supra_plot, assay_taste, conc, conc_rank, conc_label) supra_model_plot <- supra_model_res %>% mutate( assay_taste = factor( assay_taste, levels = c("nacl", "sucr"), labels = c("NaCL", "Sucrose")) ) %>% left_join( conc_labs, by = c("assay_taste", "conc") ) supra_changes_plot <- supra_changes %>% mutate( assay_taste = factor( assay_taste, levels = c("nacl", "sucr"), labels = c("NaCL", "Sucrose")) ) %>% left_join( conc_labs, by = c("assay_taste", "conc") ) %>% left_join( y = supra_model_plot %>% filter(days == 10) %>% select(assay_taste, conc, ypos = est), by = c("assay_taste", "conc") ) %>% left_join( y = supra_model_plot %>% filter(days == 0) %>% select(assay_taste, conc, y0pos = est), by = c("assay_taste", "conc") )
ggplot( data = supra_plot, mapping = aes(x = days, y = response) ) + geom_hline( yintercept = 0, color = "grey50", size = 0.5 ) + geom_line( aes(group = record_id), color = "grey80", size = .3) + geom_point( color = "grey80", size = 0.3 ) + geom_line( data = supra_model_plot, aes(x = days, y = est, color = assay_taste), inherit.aes = FALSE ) + geom_ribbon( data = supra_model_plot, aes(x = days, ymin = conf_lo, ymax = conf_hi, fill = assay_taste), alpha = 0.25, inherit.aes = FALSE ) + geom_segment( data = supra_changes_plot, aes(x = 0, xend = 10, y = y0pos, yend = y0pos), linetype = "dotted", size = .2, color = "grey20" ) + geom_segment( data = supra_changes_plot, aes(x = 10, xend = 10, y = ypos, yend = y0pos), linetype = "dotted", size = .2, color = "grey20" ) + geom_text( data = supra_changes_plot, aes(x = 10, y = ypos, label = label), hjust = 0, vjust = 0, size = 2, color = "grey10" ) + geom_text( data = conc_labs, aes(x = 2, y = 5, label = conc_label), hjust = 0, size = 2.5, color = "grey40" ) + facet_grid(assay_taste ~ conc_rank) + scale_y_continuous( expand = c(0, 0) ) + scale_color_manual( guide = FALSE, values = c("#7b4244", "#60e9dc") ) + scale_fill_manual( guide = FALSE, values = c("#7b4244", "#60e9dc") ) + labs( title = "Changes in suprathreshold intensity", x = "Days since start of fast", y = "Response", caption = "Each line is a single subject's trajectory.\nColored lines show GEE model fit with pointwise 90% confidence intervals.\nValue is estimated average change from Day 0 to 10 with 90% confidence interval." ) + theme_classic() + theme( strip.background = element_blank(), # strip.background.y = element_blank(), strip.text.x = element_blank(), strip.text.y = element_text(angle = 0), axis.line.x = element_blank(), axis.line.y = element_line(color = "grey50", size = 0.5), axis.ticks.y = element_line(color = "grey50", size = 0.1), axis.ticks.x = element_line(color = "grey50", size = 0.1), axis.text = element_text(size = 6), axis.title.y = element_text(size = 8), axis.title.x = element_text(size = 8), plot.caption = element_text(size = 6) )
Changes in Detection/Recognition Thresholds
dtrt_dt <- dtrt_thresholds %>% ungroup() %>% select(-dtp, -rtp, -taste_order) %>% tidyr::gather(key = "type", value = "conc", -record_id, -time, -days, -assay_taste) %>% mutate( assay_taste = factor( assay_taste, levels = c("nacl", "sucr"), labels = c("NaCL", "Sucrose")), type = factor( type, levels = c("dt", "rt"), labels = c("Detection", "Recognition") ) ) %>% group_by(assay_taste, type, record_id) %>% arrange(time, .by_group = TRUE) %>% mutate( conc2 = log1p(conc) - lag(log1p(conc), default = NA) )
dtrt_dt %>% group_by(assay_taste, time, type) %>% summarise( x = sprintf("%.1f (%.2f)", mean(conc), sd(conc)) ) %>% tidyr::spread( key = time, value = x ) %>% knitr::kable( align = c("l", "l", rep("c", 4)), col.names = c("Assay", "Threshold", "Start of Fast", "Day 5", "End of fast", "End of fast + 3 days"), caption = "Mean (Standard deviation) of thresholds at study timepoints" )
A simple paired t-test between day 0 and end of fast:
dtrt_dt %>% ungroup() %>% filter(time %in% c(1, 3)) %>% select(-days, -conc2) %>% group_by(assay_taste, type) %>% tidyr::pivot_wider( names_from = time, names_prefix = "time", values_from = conc ) %>% tidyr::nest() %>% mutate( p_value = purrr::map_dbl( .x = data, .f = ~ t.test(x = .x[["time1"]], y = .x[["time3"]], paired = TRUE)[["p.value"]]) ) %>% select( assay_taste, type, p_value )
dtrt_model_dt <- dtrt_dt %>% ungroup %>% mutate( y = log1p(conc), ) %>% select(type, y, days, assay_taste, time, c = conc, id = record_id) %>% arrange(id) %>% group_nest(type, assay_taste) %>% mutate( data_pre_to_5 = purrr::map(data, ~ .x %>% filter(time %in% c(1, 2)) ) )
geex_correct <- function(data, model){ f <- geex::grab_psiFUN(object = model, data = data) L <- matrix( c(rep(1, 20), 0:19), ncol = 2 ) function(theta){ # browser() m <- f(theta[1:2]) c(m, # GEE model expm1(L %*% theta[1:2]) - theta[-c(1:2, 23)], # estimate for days 0:19 theta[23] - (theta[13] - theta[3]) # difference between day 10 and day 0 ) } } gee_fit1 <- function(data){ geepack::geeglm( formula = y ~ days, id = id, data = data, family = gaussian(link = 'identity')) } mest <- function(fitter, data){ m <- fitter(data) m_estimate( data = data, units = 'id', estFUN = geex_correct, outer_args = list(model = m), corrections = list(fay = correction( FUN = fay_bias_correction, b = 0.75)), compute_roots = TRUE, root_control = setup_root_control(start = c(coef(m), rep(3, 20), 0)) # roots = c(coef(m), rep(3, 20)) ) } # get_est <- function(fit, ndays = 17){ # # entries <- # expand.grid( # days = 0:ndays # ) # # # Contrast matrix # # TODO: This is hardcoded based on the model specification!! # make_L <- function(d, c){ # c(1, d) # } # # L <- t(apply(entries, 1, function(x) { make_L(x[1], x[2]) })) # # cbind( # entries, # est = L %*% coef(fit), # se = sqrt(diag(L %*% vcov(fit) %*% t(L))) # ) # } get_est <- function(fit){ tibble::tibble( days = c(0:19, -10), est = fit@estimates[-c(1,2)], se = sqrt(diag(vcov(fit)))[-c(1,2)] ) } dtrt_model_dt %>% mutate( fit_alldays = purrr::map( .x = data, .f = ~ mest(gee_fit1, .x) ), est_alldays = purrr::map( .x = fit_alldays, .f = ~ get_est(.x) ) ) -> dtrt_fits dtrt_fits %>% select( type, assay_taste, est_alldays ) %>% tidyr::unnest(cols = c(est_alldays)) %>% mutate( conf_lo = est - qnorm(.9) * se, conf_hi = est + qnorm(.9) * se, est = est, conf_lo = conf_lo, conf_hi = conf_hi ) -> dtrt_model_res dtrt_changes <- dtrt_model_res %>% filter(days == -10) %>% group_by(type, assay_taste) %>% mutate(label = sprintf("%.1f (%.1f, %.1f)", est, conf_lo, conf_hi)) # dtrt_changes <- # dtrt_model_res %>% # filter(days %in% c(0, 10)) %>% # group_by(type, assay_taste) %>% # summarise( # diff = expm1(est[days == 10]) - expm1(est[days == 0]) # ) # dtrt_changes_orig_scale <- dtrt_model_res %>% # filter(days == 10) %>% # mutate( # est = expm1(est) # ) dt_nacl_changes <- filter( dtrt_changes, assay_taste == "NaCL", type == "Detection") %>% pull(label) rt_nacl_changes <- filter( dtrt_changes, assay_taste == "NaCL", type == "Recognition") %>% pull(label) dt_sucr_changes <- filter( dtrt_changes, assay_taste == "Sucrose", type == "Detection") %>% pull(label) rt_sucr_changes <- filter( dtrt_changes, assay_taste == "Sucrose", type == "Recognition") %>% pull(label)
dtrt_changes %>% select(type, assay_taste, label) %>% knitr::kable( caption = "Estimated change in threshold from day 0 to day 10 (90% confidence interval)." )
dtrt_plot <- dtrt_model_dt %>% select(type, assay_taste, data) %>% tidyr::unnest(cols = data) %>% filter(days >= 0) dtrt_model_plot <- dtrt_model_res %>% filter(days >= 0 ) dtrt_changes_plot <- dtrt_changes %>% left_join( y = dtrt_model_plot %>% filter(days == 10) %>% select(type, assay_taste, ypos = est), by = c("type", "assay_taste") ) %>% left_join( y = dtrt_model_plot %>% filter(days == 0) %>% select(type, assay_taste, y0pos = est), by = c("type", "assay_taste") ) %>% mutate(days = 10)
Figure \@ref(fig:dtrt_days) displays changes in suprathreshold intensity for salt and sweet taste during the study period for all subjects along with the estimated GEE model. Responses are highly variable, especially between subjects (Table \@ref(table1)), but show a general trend of decreasing recognition response for salt and sweet response.
ggplot( data = dtrt_plot, mapping = aes(x = days, y = expm1(y)) ) + geom_hline( yintercept = 0, color = "grey50", size = 0.5 ) + geom_line( aes(group = id), color = "grey80", size = .3) + geom_point( color = "grey80", size = 0.3 ) + geom_line( data = dtrt_model_plot, aes(x = days, y = est, color = assay_taste), inherit.aes = FALSE ) + geom_ribbon( data = dtrt_model_plot, aes(x = days, ymin = conf_lo, ymax = conf_hi, fill = assay_taste), alpha = 0.25, inherit.aes = FALSE ) + geom_segment( data = dtrt_changes_plot, aes(x = 0, xend = 10, y = y0pos, yend = y0pos), linetype = "dotted", size = .2, color = "grey20" ) + geom_segment( data = dtrt_changes_plot, aes(x = 10, xend = 10, y = ypos, yend = y0pos), linetype = "dotted", size = .2, color = "grey20" ) + geom_text( data = dtrt_changes_plot, aes(x = 10, y = ypos, label = label), hjust = 0, vjust = 0, size = 2, color = "grey10" ) + facet_grid(assay_taste ~ type, scales = "free") + scale_y_continuous( expand = c(0, 0) ) + scale_color_manual( guide = FALSE, values = c("#7b4244", "#60e9dc") ) + scale_fill_manual( guide = FALSE, values = c("#7b4244", "#60e9dc") ) + labs( title = "Changes in detection/recognition thresholds", x = "Days since start of fast", y = "Threshold (mM)", caption = "Each line is a single subject's trajectory.\nColored lines show GEE model fit with pointwise 90% confidence intervals.\nValue is estimated average change from Day 0 to 10 with 90% confidence interval." ) + theme_classic() + theme( strip.background = element_blank(), # strip.background.y = element_blank(), # strip.text.x = element_blank(), strip.text.y = element_text(angle = 0), axis.line.x = element_blank(), axis.line.y = element_line(color = "grey50", size = 0.5), axis.ticks.y = element_line(color = "grey50", size = 0.1), axis.ticks.x = element_line(color = "grey50", size = 0.1), axis.text = element_text(size = 6), axis.title.y = element_text(size = 8), axis.title.x = element_text(size = 8), plot.caption = element_text(size = 6) )
dtrt_dt %>% filter(time %in% c(1, 2)) %>% ggplot( aes(x = factor(days), y = conc, group = record_id) ) + geom_point(size = 0.2) + geom_line(size = 0.2) + facet_grid( assay_taste ~ type, scales = "free_y" )
baseline_vars <- clinical_summary %>% select(record_id, bmi, age, male) hold_thresholds <- dtrt_thresholds %>% filter(time %in% c(1, 4)) %>% ungroup() %>% select(-taste_order, -days, -dtp, -rtp) %>% tidyr::pivot_longer(cols = c(dt, rt)) %>% group_by(record_id, assay_taste, name) %>% summarise( value = value[time == 4] - value[time == 1] ) %>% mutate( time = 41 ) %>% bind_rows( select(ungroup(dtrt_thresholds), -taste_order, -days, -dtp, -rtp) %>% tidyr::pivot_longer(cols = c(dt, rt)) ) hold_for_cors <- hold_thresholds %>% filter(time %in% c(1, 41)) %>% left_join( baseline_vars, by = "record_id" ) hold_for_cors %>% tidyr::pivot_longer( cols = c(age, male, bmi), names_to = "baseline_var", values_to = "baseline_value" ) %>% ungroup() %>% group_nest( time, assay_taste, name, baseline_var ) %>% mutate( cort = purrr::map( .x = data, .f = ~ cor.test(.x$value, .x$baseline_value)), cor_est = purrr::map_dbl(cort, "estimate"), cor_pval = purrr::map_dbl(cort, "p.value"), value = sprintf("%.3f (%.2f)", cor_est, cor_pval) ) %>% select(-data, -cort, -cor_est, -cor_pval) %>% tidyr::pivot_wider( id_cols = c(time, assay_taste, name), names_from = baseline_var ) %>% mutate( time = if_else(time == 1, "Baseline", "Change from baseline to re-feed") ) %>% knitr::kable( caption = "Correlations between baseline covariates and dectection and recognition thresholds." )
Changes in food-liking
flq_completed <- foodliking %>% filter(complete.cases(.)) %>% group_by(record_id) %>% filter(n() == 2) %>% ungroup() %>% mutate(record_id = as.character(record_id))
There were r length(unique(flq_completed$record_id)) subjects who completed the food-liking questionnaire both before and after the exposure time.
TODO: describe the scoring
Figure XX shows the changes in food-liking.
flq_plot <- flq_completed %>% tidyr::pivot_longer( cols = ends_with("flq") ) %>% mutate( name = stringr::str_remove(name, "_flq"), name = factor(name, levels = c("salt", "sweet", "fatsalt", "fatsweet"), labels = c("Salt", "Sweet", "Fat & Salt", "Fat & Sweet")) ) hold <- flq_plot %>% tidyr::pivot_wider( id_cols = c("record_id", "name"), names_from = "time" ) %>% ungroup() %>% group_nest(name) %>% mutate( ttest = purrr::map( .x = data, .f = ~ t.test(.x$`5`, .x$`0`, paired = TRUE) ), tpval = purrr::map_dbl(ttest, "p.value"), ttext = sprintf("Paired t-test p-value: %.3f", tpval) )
hold %>% mutate( estimate = purrr::map(ttest, ~ .x$estimate) ) %>% select(name, estimate) %>% knitr::kable(caption = "Paired Changes in FLQ")
flq_plot %>% ggplot( data = ., aes(x = time, y = value, color = name) ) + geom_hline( yintercept = 0, color = "grey50" ) + geom_point( alpha = 0.5 ) + geom_line( aes(group = record_id), size = 0.1, alpha = 0.5 ) + geom_text( data = hold, aes(label = ttext, x = 2, y = 1), color = "black", size = 2 ) + stat_smooth( # aes(color = assay_taste), method = "lm", se = FALSE ) + scale_y_continuous( "FLQ Score" ) + scale_x_continuous( "Time", breaks = c(0, 5), labels = c("Pre", "Post") ) + scale_color_manual( "", values = c("#8c75d2", "#11677e", "#d547ca", "#3e3c8d"), guide = FALSE ) + facet_wrap( ~ name, ncol = 2 ) + theme_classic() + theme( strip.background = element_blank(), # strip.background.y = element_blank(), # strip.text.x = element_blank(), strip.text.y = element_text(angle = 0), axis.line.x = element_blank(), axis.line.y = element_line(color = "grey50", size = 0.5), axis.ticks.y = element_line(color = "grey50", size = 0.1), axis.ticks.x = element_line(color = "grey50", size = 0.1), axis.text = element_text(size = 6), axis.title.y = element_text(size = 8), axis.title.x = element_text(size = 8), plot.caption = element_text(size = 6) )
Food-liking and Adaptation
dtrt_hold <- dtrt_thresholds %>% ungroup() %>% filter(time %in% c(1, 4), record_id %in% flq_completed$record_id) %>% select(record_id, time, assay_taste, dt, rt) %>% tidyr::pivot_longer( cols = c(dt, rt), names_to = "threshold" ) %>% group_by(record_id, assay_taste, threshold) %>% summarise( dtrt_del = value[time ==4] - value[time == 1] ) %>% ungroup() %>% mutate( threshold = factor(threshold, levels = c("dt", "rt"), labels = c("Detection", "Recognition")), assay_taste = factor(assay_taste, levels = c("nacl", "sucr"), labels = c("NaCL", "Sucrose")) ) flq_hold <- flq_completed %>% ungroup() %>% filter(time %in% c(0, 5)) %>% tidyr::pivot_longer( cols = ends_with("_flq"), names_to = "flq_name" ) %>% group_by(record_id, flq_name) %>% summarise( flq_del = value[time == 5] - value[time == 0] ) %>% mutate( flq_name = stringr::str_remove(flq_name, "_flq"), flq_name = factor(flq_name, levels = c("salt", "sweet", "fatsalt", "fatsweet"), labels = c("Salt", "Sweet", "Fat & Salt", "Fat & Sweet")) ) plot_hold <- full_join( dtrt_hold, flq_hold, by = "record_id" )
ggplot( data = plot_hold, aes(x = dtrt_del, y = flq_del) ) + stat_smooth( method = "lm", size = 0.5, se = TRUE ) + geom_point( aes(color = flq_name), alpha = 0.5 ) + scale_color_manual( "", values = c("#8c75d2", "#11677e", "#d547ca", "#3e3c8d"), guide = FALSE ) + scale_y_continuous( "Change in food liking score" ) + scale_x_continuous( "Change in threshold" ) + facet_grid( flq_name ~ threshold + assay_taste, scales = "free" ) + theme_classic() + theme( panel.background = element_rect(fill = "grey90"), panel.grid.major = element_line(color = "white"), strip.background = element_blank(), # strip.background.y = element_blank(), # strip.text.x = element_blank(), strip.text.y = element_text(angle = 0), # axis.line.x = element_blank(), axis.line.x = element_line(color = "grey50", size = 0.5), axis.line.y = element_line(color = "grey50", size = 0.5), axis.ticks.y = element_line(color = "grey50", size = 0.1), axis.ticks.x = element_line(color = "grey50", size = 0.1), axis.text = element_text(size = 6), axis.title.y = element_text(size = 8), axis.title.x = element_text(size = 8), plot.caption = element_text(size = 6) )
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.