age X bone turnover

A. outcome: osteocalcin [bone formation]

1. GAM with smooth age, sex, age*sex interaction

plot gam-predicted lines by average daily minutes of light- and moderate-intensity activity tertiles facetted by sex

# create tertiles for minutes
dat <- dat %>%
  filter(!is.na(lig_mod_pooled_mins)) %>%
  mutate(tertile = ntile(lig_mod_pooled_mins, 3))

#mod1 <- gam(log(osteocalcin_pg_ml) ~ s(age_years) + sex + age_years*sex, data = dat)

age_seq <- seq(min(dat$age_years), max(dat$age_years), length.out = 100)
dat$sex  <- as.factor(dat$sex)

pred_grid <- dat %>%
  group_by(tertile) %>%
  group_split() %>%
  map_df(~ {
    data_tertile <- .
    model <- gam(log(osteocalcin_pg_ml) ~ s(age_years) + sex + age_years*sex, data = data_tertile)
    pred_grid <- expand.grid(age_years = age_seq, sex = unique(data_tertile$sex))
    pred <- predict(model, newdata = pred_grid, se.fit=TRUE)
    pred_grid$fit <- exp(pred$fit)
    upper <- pred$fit + 1.96 * pred$se.fit
    lower <- pred$fit - 1.96 * pred$se.fit
    pred_grid$lwr  <- exp(lower)
    pred_grid$upr  <- exp(upper)
    pred_grid$tertile <- unique(data_tertile$tertile)
    pred_grid
  })

facet_labs <- c("female" = "Women", "male" ="Men")
tert_labs <- c("1" = "1st", "2" = "2nd", "3" = "3rd")

ggplot(pred_grid, aes(x = age_years, y = fit, color = factor(tertile))) +
  geom_point(data = dat, aes(x = age_years, y = osteocalcin_pg_ml, color = factor(tertile)), size = 1, alpha = 0.5) +
  geom_line(aes(x = age_years, y = fit, color = factor(tertile)), size = 1) +
  geom_ribbon(aes(x = age_years, ymin = lwr, ymax = upr, fill=factor(tertile)), color=NA) + 
  facet_wrap(~ sex, labeller = labeller(sex= facet_labs)) +
  labs(title = NULL,
       x = "Age (years)",
       y = "Osteocalcin (pg/mL)",
       color = "Activity Tertile", fill = "Activity Tertile")+
  scale_color_manual(name = "Activity (Minutes) Tertile",
                     values = c("1" = "#CC3D12", "2" = "#283845", "3" = "#C5A359"),
                     labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_fill_manual(name = "Activity (Minutes) Tertile",
                    values = c("1" = "#CC3D1233", "2" = "#28384533", "3" = "#C5A35933"),
                    labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_y_continuous(labels = format_thousands) +
  coord_cartesian(ylim = c(0, 200000)) +
  theme(legend.position="top", legend.title = element_text(size = 12, face = "bold", color="#000000"), 
    legend.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot gam-predicted lines by average daily steps tertiles facetted by sex

# create tertiles for steps
dat <- dat %>%
  filter(!is.na(ad_step_count_sum_0_24hr)) %>%
  mutate(tertile = ntile(ad_step_count_sum_0_24hr, 3))

#mod1 <- gam(log(osteocalcin_pg_ml) ~ s(age_years) + sex + age_years*sex, data = dat)

age_seq <- seq(min(dat$age_years), max(dat$age_years), length.out = 100)
dat$sex  <- as.factor(dat$sex)

pred_grid <- dat %>%
  group_by(tertile) %>%
  group_split() %>%
  map_df(~ {
    data_tertile <- .
    model <- gam(log(osteocalcin_pg_ml) ~ s(age_years) + sex + age_years*sex, data = data_tertile)
    pred_grid <- expand.grid(age_years = age_seq, sex = unique(data_tertile$sex))
    pred <- predict(model, newdata = pred_grid, se.fit=TRUE)
    pred_grid$fit <- exp(pred$fit)
    upper <- pred$fit + 1.96 * pred$se.fit
    lower <- pred$fit - 1.96 * pred$se.fit
    pred_grid$lwr  <- exp(lower)
    pred_grid$upr  <- exp(upper)
    pred_grid$tertile <- unique(data_tertile$tertile)
    pred_grid
  })

facet_labs <- c("female" = "Women", "male" ="Men")
tert_labs <- c("1" = "1st", "2" = "2nd", "3" = "3rd")

ggplot(pred_grid, aes(x = age_years, y = fit, color = factor(tertile))) +
  geom_point(data = dat, aes(x = age_years, y = osteocalcin_pg_ml, color = factor(tertile)), size = 1, alpha = 0.5) +
  geom_line(aes(x = age_years, y = fit, color = factor(tertile)), size = 1) +
  geom_ribbon(aes(x = age_years, ymin = lwr, ymax = upr, fill=factor(tertile)), color=NA) + 
  facet_wrap(~ sex, labeller = labeller(sex= facet_labs)) +
  labs(title = NULL,
       x = "Age (years)",
       y = "Osteocalcin (pg/mL)",
       color = "Steps Tertile", fill = "Steps Tertile")+
 scale_color_manual(name = "Activity (Steps) Tertile",
                     values = c("1" = "#CC3D12", "2" = "#283845", "3" = "#C5A359"),
                     labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_fill_manual(name = "Activity (Steps) Tertile",
                    values = c("1" = "#CC3D1233", "2" = "#28384533", "3" = "#C5A35933"),
                    labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  coord_cartesian(ylim = c(0, 200000)) +
  scale_y_continuous(labels = format_thousands) +
  theme(legend.position="top", legend.title = element_text(size = 12, face = "bold", color="#000000"), 
    legend.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"))

B. outcome: ctx-1 [bone resorption]

1. GAM with smooth age, sex, age*sex interaction

plot gam-predicted lines by average daily minutes of light- and moderate-intensity activity tertiles facetted by sex

# create tertiles for minutes
dat <- dat %>%
  filter(!is.na(lig_mod_pooled_mins)) %>%
  mutate(tertile = ntile(lig_mod_pooled_mins, 3))

#mod1 <- gam(log(ctx_1_ng_ml) ~ s(age_years) + sex + age_years*sex, data = dat)

age_seq <- seq(min(dat$age_years), max(dat$age_years), length.out = 100)
dat$sex  <- as.factor(dat$sex)

pred_grid <- dat %>%
  group_by(tertile) %>%
  group_split() %>%
  map_df(~ {
    data_tertile <- .
    model <- gam(log(ctx_1_ng_ml) ~ s(age_years) + sex + age_years*sex, data = data_tertile)
    pred_grid <- expand.grid(age_years = age_seq, sex = unique(data_tertile$sex))
    pred <- predict(model, newdata = pred_grid, se.fit=TRUE)
    pred_grid$fit <- exp(pred$fit)
    upper <- pred$fit + 1.96 * pred$se.fit
    lower <- pred$fit - 1.96 * pred$se.fit
    pred_grid$lwr  <- exp(lower)
    pred_grid$upr  <- exp(upper)
    pred_grid$tertile <- unique(data_tertile$tertile)
    pred_grid
  })

facet_labs <- c("female" = "Women", "male" ="Men")
tert_labs <- c("1" = "1st", "2" = "2nd", "3" = "3rd")

ggplot(pred_grid, aes(x = age_years, y = fit, color = factor(tertile))) +
  geom_point(data = dat, aes(x = age_years, y = ctx_1_ng_ml, color = factor(tertile)), size = 1, alpha = 0.5) +
  geom_line(aes(x = age_years, y = fit, color = factor(tertile)), size = 1) +
  geom_ribbon(aes(x = age_years, ymin = lwr, ymax = upr, fill=factor(tertile)), color=NA) + 
  facet_wrap(~ sex, labeller = labeller(sex= facet_labs)) +
  labs(title = NULL,
       x = "Age (years)",
       y = "CTX-1 (ng/mL)",
       color = "Activity Tertile", fill = "Activity Tertile")+
   scale_color_manual(name = "Activity (Minutes) Tertile",
                     values = c("1" = "#CC3D12", "2" = "#283845", "3" = "#C5A359"),
                     labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_fill_manual(name = "Activity (Minutes) Tertile",
                    values = c("1" = "#CC3D1233", "2" = "#28384533", "3" = "#C5A35933"),
                    labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  coord_cartesian(ylim = c(0, 1)) +
  theme(legend.position="top", legend.title = element_text(size = 12, face = "bold", color="#000000"), 
    legend.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"))

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).

plot gam-predicted lines by average daily steps tertiles facetted by sex

# create tertiles for steps
dat <- dat %>%
  filter(!is.na(ad_step_count_sum_0_24hr)) %>%
  mutate(tertile = ntile(ad_step_count_sum_0_24hr, 3))

#mod1 <- gam(log(ctx_1_ng_ml) ~ s(age_years) + sex + age_years*sex, data = dat)

age_seq <- seq(min(dat$age_years), max(dat$age_years), length.out = 100)
dat$sex  <- as.factor(dat$sex)

pred_grid <- dat %>%
  group_by(tertile) %>%
  group_split() %>%
  map_df(~ {
    data_tertile <- .
    model <- gam(log(ctx_1_ng_ml) ~ s(age_years) + sex + age_years*sex, data = data_tertile)
    pred_grid <- expand.grid(age_years = age_seq, sex = unique(data_tertile$sex))
    pred <- predict(model, newdata = pred_grid, se.fit=TRUE)
    pred_grid$fit <- exp(pred$fit)
    upper <- pred$fit + 1.96 * pred$se.fit
    lower <- pred$fit - 1.96 * pred$se.fit
    pred_grid$lwr  <- exp(lower)
    pred_grid$upr  <- exp(upper)
    pred_grid$tertile <- unique(data_tertile$tertile)
    pred_grid
  })

facet_labs <- c("female" = "Women", "male" ="Men")
tert_labs <- c("1" = "1st", "2" = "2nd", "3" = "3rd")

ggplot(pred_grid, aes(x = age_years, y = fit, color = factor(tertile))) +
  geom_point(data = dat, aes(x = age_years, y = ctx_1_ng_ml, color = factor(tertile)), size = 1, alpha = 0.5) +
  geom_line(aes(x = age_years, y = fit, color = factor(tertile)), size = 1) +
  geom_ribbon(aes(x = age_years, ymin = lwr, ymax = upr, fill=factor(tertile)), color=NA) + 
  facet_wrap(~ sex, labeller = labeller(sex= facet_labs)) +
  labs(title = NULL,
       x = "Age (years)",
       y = "CTX-1 (ng/mL)",
       color = "Steps Tertile", fill = "Steps Tertile")+
   scale_color_manual(name = "Activity (Steps) Tertile",
                     values = c("1" = "#CC3D12", "2" = "#283845", "3" = "#C5A359"),
                     labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_fill_manual(name = "Activity (Steps) Tertile",
                    values = c("1" = "#CC3D1233", "2" = "#28384533", "3" = "#C5A35933"),
                    labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  coord_cartesian(ylim = c(0, 1)) +
  theme(legend.position="top", legend.title = element_text(size = 12, face = "bold", color="#000000"), 
    legend.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"))

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).

C. outcome: left tibial SOS [proxy BMD]

1. GAM with smooth age terms, sex, and age*sex interaction

plot gam-predicted lines by average daily light- and moderate-intensity activity minute tertiles facetted by sex

# create tertiles for minutes
dat_full <- dat_full %>%
  filter(!is.na(lig_mod_pooled_mins)) %>%
  mutate(tertile = ntile(lig_mod_pooled_mins, 3))

#mod1 <- gam(log(left_tibia_sos) ~ s(age_years) + sex + age_years*sex, data = dat_full)

age_seq <- seq(min(dat_full$age_years), max(dat_full$age_years), length.out = 100)
dat_full$sex  <- as.factor(dat_full$sex)

pred_grid <- dat_full %>%
  group_by(tertile) %>%
  group_split() %>%
  map_df(~ {
    data_tertile <- .
    model <- gam(log(left_tibia_sos) ~ s(age_years) + sex + age_years*sex, data = data_tertile)
    pred_grid <- expand.grid(age_years = age_seq, sex = unique(data_tertile$sex))
    pred <- predict(model, newdata = pred_grid, se.fit=TRUE)
    pred_grid$fit <- exp(pred$fit)
    upper <- pred$fit + 1.96 * pred$se.fit
    lower <- pred$fit - 1.96 * pred$se.fit
    pred_grid$lwr  <- exp(lower)
    pred_grid$upr  <- exp(upper)
    pred_grid$tertile <- unique(data_tertile$tertile)
    pred_grid
  })

facet_labs <- c("female" = "Women", "0" ="Men")
tert_labs <- c("1" = "1st", "2" = "2nd", "3" = "3rd")

ggplot(pred_grid, aes(x = age_years, y = fit, color = factor(tertile))) +
  geom_point(data = dat_full, aes(x = age_years, y = left_tibia_sos, color = factor(tertile)), size = 1, alpha = 0.5) +
  geom_line(aes(x = age_years, y = fit, color = factor(tertile)), size = 1) +
  geom_ribbon(aes(x = age_years, ymin = lwr, ymax = upr, fill=factor(tertile)), color=NA) + 
  facet_wrap(~ sex, labeller = labeller(sex= facet_labs)) +
  labs(title = NULL,
       x = "Age (years)",
       y = "Tibial SOS (m/s)",
       color = "Activity Tertile", fill = "Activity Tertile")+
   scale_color_manual(name = "Activity (Minutes) Tertile",
                     values = c("1" = "#CC3D12", "2" = "#283845", "3" = "#C5A359"),
                     labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_fill_manual(name = "Activity (Minutes) Tertile",
                    values = c("1" = "#CC3D1233", "2" = "#28384533", "3" = "#C5A35933"),
                    labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  #coord_cartesian(ylim = c(0, 1)) +
  theme(legend.position="top", legend.title = element_text(size = 12, face = "bold", color="#000000"), 
    legend.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"))

## Warning: Removed 282 rows containing missing values or values outside the scale range
## (`geom_point()`).

plot gam-predicted lines by average daily steps tertiles facetted by sex

# create tertiles for steps
dat_full <- dat_full %>%
  filter(!is.na(ad_step_count_sum_0_24hr)) %>%
  mutate(tertile = ntile(ad_step_count_sum_0_24hr, 3))

#mod1 <- gam(log(left_tibia_sos) ~ s(age_years) + sex + age_years*sex, data = dat_full)

age_seq <- seq(min(dat_full$age_years), max(dat_full$age_years), length.out = 100)
dat_full$sex  <- as.factor(dat_full$sex)

pred_grid <- dat_full %>%
  group_by(tertile) %>%
  group_split() %>%
  map_df(~ {
    data_tertile <- .
    model <- gam(log(left_tibia_sos) ~ s(age_years) + sex + age_years*sex, data = data_tertile)
    pred_grid <- expand.grid(age_years = age_seq, sex = unique(data_tertile$sex))
    pred <- predict(model, newdata = pred_grid, se.fit=TRUE)
    pred_grid$fit <- exp(pred$fit)
    upper <- pred$fit + 1.96 * pred$se.fit
    lower <- pred$fit - 1.96 * pred$se.fit
    pred_grid$lwr  <- exp(lower)
    pred_grid$upr  <- exp(upper)
    pred_grid$tertile <- unique(data_tertile$tertile)
    pred_grid
  })

facet_labs <- c("0" = "Women", "1" ="Men")
tert_labs <- c("1" = "1st", "2" = "2nd", "3" = "3rd")

ggplot(pred_grid, aes(x = age_years, y = fit, color = factor(tertile))) +
  geom_point(data = dat_full, aes(x = age_years, y = left_tibia_sos, color = factor(tertile)), size = 1, alpha = 0.5) +
  geom_line(aes(x = age_years, y = fit, color = factor(tertile)), size = 1) +
  geom_ribbon(aes(x = age_years, ymin = lwr, ymax = upr, fill=factor(tertile)), color=NA) + 
  facet_wrap(~ sex, labeller = labeller(sex= facet_labs)) +
  labs(title = NULL,
       x = "Age (years)",
       y = "Tibial SOS (m/s)",
       color = "Steps Tertile", fill = "Steps Tertile")+
   scale_color_manual(name = "Activity (Steps) Tertile",
                     values = c("1" = "#CC3D12", "2" = "#283845", "3" = "#C5A359"),
                     labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  scale_fill_manual(name = "Activity (Steps) Tertile",
                    values = c("1" = "#CC3D1233", "2" = "#28384533", "3" = "#C5A35933"),
                    labels = c("1" = "Low", "2" = "Middle", "3" = "High")) +
  #coord_cartesian(ylim = c(0, 1)) +
  theme(legend.position="top", legend.title = element_text(size = 12, face = "bold", color="#000000"), 
    legend.text = element_text(size = 12), legend.key.size = unit(0.5, "cm"))

## Warning: Removed 277 rows containing missing values or values outside the scale range
## (`geom_point()`).