bone turnover X PA linear models
cecilia sena
2025-04-15
I. Assess distributions
plot(density(na.omit(dat$osteocalcin_pg_ml))) # lognormal
plot(density(na.omit(dat$ctx_1_ng_ml))) # lognormal
plot(density(na.omit(dat$ad_step_count_sum_0_24hr))) # normal
plot(density(na.omit(dat$dur_day_total_lig_min_pla))) # normal
plot(density(na.omit(dat$dur_day_total_mod_min_pla))) # normal
plot(density(na.omit(dat$dur_day_total_vig_min_pla))) # lognormal
plot(density(na.omit(dat$functional_status_n_y_0_1))) # binary
plot(density(na.omit(dat$pregnat_or_breastfeeding_0_1))) # binary
plot(density(na.omit(dat_full$left_tibia_sos))) # normal
plot(density(na.omit(dat_full$ad_step_count_sum_0_24hr))) # normal
II. Linear models
A. exposure: average daily steps
1. outcome: osteocalcin [bone formation]
i. Model summary
mod1<- lm(log(osteocalcin_pg_ml) ~ ad_step_count_sum_0_24hr * bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1, data = dat)
summary(mod1)##
## Call:
## lm(formula = log(osteocalcin_pg_ml) ~ ad_step_count_sum_0_24hr *
## bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 +
## functional_status_n_y_0_1, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.28788 -0.34024 0.00395 0.24898 1.68773
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 11.249879735520 0.568320079074 19.795
## ad_step_count_sum_0_24hr -0.000010848036 0.000046497241 -0.233
## bmi -0.030192558432 0.022260966634 -1.356
## age_years -0.001720612143 0.002132413593 -0.807
## sexmale -0.094848238768 0.072800563655 -1.303
## pregnat_or_breastfeeding_0_1 -0.147591227794 0.097328242150 -1.516
## functional_status_n_y_0_1 0.085049811041 0.115121617372 0.739
## ad_step_count_sum_0_24hr:bmi -0.000000001123 0.000001924238 -0.001
## Pr(>|t|)
## (Intercept) <0.0000000000000002 ***
## ad_step_count_sum_0_24hr 0.816
## bmi 0.176
## age_years 0.421
## sexmale 0.194
## pregnat_or_breastfeeding_0_1 0.131
## functional_status_n_y_0_1 0.461
## ad_step_count_sum_0_24hr:bmi 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4985 on 239 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1069, Adjusted R-squared: 0.08074
## F-statistic: 4.087 on 7 and 239 DF, p-value: 0.0002944ii. Raw data
ggplot(dat, aes(x = ad_step_count_sum_0_24hr, y = log(osteocalcin_pg_ml))) +
geom_point(alpha = 0.5) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(x="Average daily step count",
y="Log Serum Osteocalcin\n(pg/mL)") 
iii. Diagnostics
plot(mod1)
iv. Model prediction (marginal effects, averaging over covariates)
(compute predicted outcome for each data point, average those predicted outcomes for each exposure value)
plot_model(mod1, type='pred', terms=c("ad_step_count_sum_0_24hr")) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(title=NULL,
x="Average daily step count",
y="Predicted Serum Osteocalcin\n(pg/mL)") ## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the transformed
## scale.
v. Model prediction (conditional effects, average values for covariates) (sanity check)
dat_pred <- data.frame(
ad_step_count_sum_0_24hr = seq(min(dat$ad_step_count_sum_0_24hr, na.rm = TRUE), max(dat$ad_step_count_sum_0_24hr, na.rm = TRUE), length.out = 100),
bmi = mean(dat$bmi, na.rm = TRUE),
age_years = mean(dat$age_years, na.rm = TRUE),
sex = "female",
pregnat_or_breastfeeding_0_1 = mean(dat$pregnat_or_breastfeeding_0_1, na.rm = TRUE),
functional_status_n_y_0_1 = mean(dat$functional_status_n_y_0_1, na.rm = TRUE)
)
preds <- predict(mod1, newdata = dat_pred, interval = "confidence")
dat_pred <- cbind(dat_pred, preds)
ggplot(dat_pred, aes(x = ad_step_count_sum_0_24hr, y = fit)) +
geom_line() +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2)
They look identical (because none of the covariates have significant large effects) so the rest will use marginal effects
vi. predicted effect with actual data overlain for scale
plot_model(mod1, type='pred', terms=c("ad_step_count_sum_0_24hr"), show.data = TRUE) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(title=NULL,
x="Average daily step count",
y="Predicted Serum Osteocalcin\n(pg/mL)") ## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the transformed
## scale.## Data points may overlap. Use the `jitter` argument to add some amount of
## random variation to the location of data points and avoid overplotting.
2. outcome: ctx-1 [bone resorption]
i. Model summary
mod1<- lm(log(ctx_1_ng_ml) ~ ad_step_count_sum_0_24hr * bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1, data = dat)
summary(mod1)##
## Call:
## lm(formula = log(ctx_1_ng_ml) ~ ad_step_count_sum_0_24hr * bmi +
## bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1,
## data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2581 -0.3994 0.1270 0.4704 1.7138
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.100546674 0.904964962 -0.111 0.9116
## ad_step_count_sum_0_24hr -0.000027046 0.000075730 -0.357 0.7213
## bmi -0.057443173 0.035689003 -1.610 0.1089
## age_years -0.005788849 0.003367586 -1.719 0.0869 .
## sexmale -0.161866082 0.114965521 -1.408 0.1605
## pregnat_or_breastfeeding_0_1 -0.329694439 0.158905579 -2.075 0.0391 *
## functional_status_n_y_0_1 0.064988964 0.180783846 0.359 0.7196
## ad_step_count_sum_0_24hr:bmi -0.000000068 0.000003157 -0.022 0.9828
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7814 on 233 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.1504, Adjusted R-squared: 0.1248
## F-statistic: 5.89 on 7 and 233 DF, p-value: 0.000002584ii. Raw data
ggplot(dat, aes(x = ad_step_count_sum_0_24hr, y = log(ctx_1_ng_ml))) +
geom_point(alpha = 0.5) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(x="Average daily step count",
y="Log Serum CTX-1\n(ng/mL)") ## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
iii. Diagnostics
plot(mod1)
iv. Model prediction
plot_model(mod1, type='pred', terms=c("ad_step_count_sum_0_24hr")) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(title=NULL,
x="Average daily step count",
y="Predicted Serum CTX-1\n(ng/mL)") ## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the transformed
## scale.
3. outcome: tibial SOS [proxy BMD]
i. Model summary
mod1<- lm(log(left_tibia_sos) ~ ad_step_count_sum_0_24hr * bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1, data = dat_full)
summary(mod1)##
## Call:
## lm(formula = log(left_tibia_sos) ~ ad_step_count_sum_0_24hr *
## bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 +
## functional_status_n_y_0_1, data = dat_full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19348 -0.02382 0.00325 0.02522 0.09795
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 8.30021792587 0.03395128618 244.474
## ad_step_count_sum_0_24hr 0.00000171926 0.00000283203 0.607
## bmi -0.00139077969 0.00123089625 -1.130
## age_years -0.00054821777 0.00014984261 -3.659
## sexmale 0.01857990668 0.00450922218 4.120
## pregnat_or_breastfeeding_0_1 0.01095884335 0.00618456028 1.772
## functional_status_n_y_0_1 0.00104387633 0.00719715497 0.145
## ad_step_count_sum_0_24hr:bmi -0.00000008693 0.00000011577 -0.751
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## ad_step_count_sum_0_24hr 0.544154
## bmi 0.259216
## age_years 0.000288 ***
## sexmale 0.0000462 ***
## pregnat_or_breastfeeding_0_1 0.077181 .
## functional_status_n_y_0_1 0.884754
## ad_step_count_sum_0_24hr:bmi 0.453180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03991 on 390 degrees of freedom
## (329 observations deleted due to missingness)
## Multiple R-squared: 0.1847, Adjusted R-squared: 0.17
## F-statistic: 12.62 on 7 and 390 DF, p-value: 0.00000000000001304ii. Raw data
ggplot(dat_full, aes(x = ad_step_count_sum_0_24hr, y = left_tibia_sos)) +
geom_point(alpha = 0.5) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(x="Average daily step count",
y="Tibial SOS (m/s))") ## Warning: Removed 329 rows containing missing values or values outside the scale range
## (`geom_point()`).
iii. Diagnostics
plot(mod1)
iv. Model prediction
plot_model(mod1, type='pred', terms=c("ad_step_count_sum_0_24hr")) +
scale_x_continuous(labels = function(ad_step_count_sum_0_24hr) format_thousands(ad_step_count_sum_0_24hr)) +
labs(title=NULL,
x="Average daily step count",
y="Predicted Tibial SOS\n(m/s)") ## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the transformed
## scale.
B. exposure: average daily light- & moderate-intensity PA minutes
1. outcome: osteocalcin [bone formation]
i. Model summary
mod1<- lm(log(osteocalcin_pg_ml) ~ lig_mod_pooled_mins * bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1, data = dat)
summary(mod1)##
## Call:
## lm(formula = log(osteocalcin_pg_ml) ~ lig_mod_pooled_mins * bmi +
## bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1,
## data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.25262 -0.32571 0.02258 0.25676 1.67501
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 11.56659952 0.73396530 15.759
## lig_mod_pooled_mins -0.00110267 0.00168735 -0.653
## bmi -0.03971263 0.02919459 -1.360
## age_years -0.00121736 0.00212051 -0.574
## sexmale -0.11185087 0.07433708 -1.505
## pregnat_or_breastfeeding_0_1 -0.14283651 0.09700774 -1.472
## functional_status_n_y_0_1 0.07418778 0.11503324 0.645
## lig_mod_pooled_mins:bmi 0.00002411 0.00006815 0.354
## Pr(>|t|)
## (Intercept) <0.0000000000000002 ***
## lig_mod_pooled_mins 0.514
## bmi 0.175
## age_years 0.566
## sexmale 0.134
## pregnat_or_breastfeeding_0_1 0.142
## functional_status_n_y_0_1 0.520
## lig_mod_pooled_mins:bmi 0.724
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4975 on 239 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1106, Adjusted R-squared: 0.08454
## F-statistic: 4.245 on 7 and 239 DF, p-value: 0.000194ii. Raw data
ggplot(dat, aes(x = lig_mod_pooled_mins, y = log(osteocalcin_pg_ml))) +
geom_point(alpha = 0.5) +
labs(x="Average daily minutes of light- and moderate-intensity PA",
y="Log Serum Osteocalcin\n(pg/mL)") 
iii. Diagnostics
plot(mod1)
iv. Model prediction
plot_model(mod1, type='pred', terms=c("lig_mod_pooled_mins")) +
labs(title=NULL,
x="Average daily minutes of light- and moderate-intensity PA",
y="Predicted Serum Osteocalcin\n(pg/mL)") ## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the transformed
## scale.
2. outcome: ctx-1 [bone resorption]
i. Model summary
mod1<- lm(log(ctx_1_ng_ml) ~ lig_mod_pooled_mins * bmi + bmi + age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1, data = dat)
summary(mod1)##
## Call:
## lm(formula = log(ctx_1_ng_ml) ~ lig_mod_pooled_mins * bmi + bmi +
## age_years + sex + pregnat_or_breastfeeding_0_1 + functional_status_n_y_0_1,
## data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.2588 -0.3891 0.1351 0.4825 1.6942
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.29278075 1.18200272 0.248 0.8046
## lig_mod_pooled_mins -0.00184662 0.00276452 -0.668 0.5048
## bmi -0.07093807 0.04731028 -1.499 0.1351
## age_years -0.00463657 0.00335576 -1.382 0.1684
## sexmale -0.18582555 0.11739006 -1.583 0.1148
## pregnat_or_breastfeeding_0_1 -0.31312480 0.15862436 -1.974 0.0496 *
## functional_status_n_y_0_1 0.04298348 0.18108123 0.237 0.8126
## lig_mod_pooled_mins:bmi 0.00003477 0.00011242 0.309 0.7574
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7815 on 233 degrees of freedom
## (7 observations deleted due to missingness)
## Multiple R-squared: 0.1501, Adjusted R-squared: 0.1246
## F-statistic: 5.88 on 7 and 233 DF, p-value: 0.000002651ii. Raw data
ggplot(dat, aes(x = lig_mod_pooled_mins, y = log(ctx_1_ng_ml))) +
geom_point(alpha = 0.5) +
labs(x="Average daily minutes of light- and moderate-intensity PA",
y="Log Serum CTX-1\n(ng/mL)") ## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
iii. Diagnostics
plot(mod1)
iv. Model prediction
plot_model(mod1, type='pred', terms=c("lig_mod_pooled_mins")) +
labs(title=NULL,
x="Average daily minutes of light- and moderate-intensity PA",
y="Predicted Serum CTX-1\n(ng/mL)") ## Model has log-transformed response. Back-transforming predictions to
## original response scale. Standard errors are still on the transformed
## scale.