Diet covariates

library(arivale.data.interface)

list_snapshot_contents()

qs <- get_snapshot("assessments", clean=T)
diet_variables <- names(qs)[grepl("diet_freq", names(qs))]
lifestyle_variables <- c(
    "assessment_lifestyle_cruciferous_vegetables_enum",
    "assessment_lifestyle_fruits_enum",                      
    "assessment_lifestyle_vegetables_enum",                  
    "assessment_lifestyle_sugary_drinks_enum",               
    "assessment_lifestyle_water_enum",                       
    "assessment_lifestyle_alcohol_drinks_a_day_enum",        
    "assessment_lifestyle_grains_enum"
)
variables <- c(diet_variables, lifestyle_variables)
qs <- qs[, c("public_client_id", variables), with=FALSE]
qs[, "complete" := apply(qs, 1, function(x) sum(!is.na(x))/length(x))]
qs <- qs[order(-complete)][!duplicated(public_client_id)]

scales <- apply(qs[, variables, with=F], 2, function(x) unique(x, na.rm=T) %>% 
          sort() %>% paste(collapse = ", ")) %>% data.table()
scales <- data.table(quantity=variables, scale=scales)
fwrite(scales, "data/diet_scales.csv")

All of those should be on a ordinal scale so we will conver them to the corresponding rank values.

library(stringr)
library(magrittr)

for (v in variables) {
    qs[[v]] <- str_match(qs[[v]], "^\\((\\d)\\)")[, 2] %>% as.numeric()
}

qs

Now let’s merge in our selected cohort.

cohort <- rbind(
    fread("no_weight_loss.csv", colClasses=c(public_client_id="character")),
    fread("successful_weight_loss.csv", colClasses=c(public_client_id="character"))
)[since_baseline == 0]
cohort[, "subset" := "controls"]
cohort[weight_change_relative < 0, "subset" := "weight loss"]
cohort[, subset := factor(subset)]

data <- qs[cohort, on="public_client_id"]

Let’s visualize what we have.

library(pheatmap)
library(viridisLite)

df <- as.data.frame(data[, variables, with=F])
rownames(df) <- data$public_client_id
anns <- data.frame(BMI=data$bmi, subset=data$subset, row.names=data$public_client_id)
pheatmap(t(df), color=viridis(256), annotation_col = anns, show_colnames=FALSE)

# Save image
pheatmap(t(df), color=viridis(256), annotation_col = anns, show_colnames=FALSE,
         filename="figures/diets.png", width=16, height=6, dpi=300)

Finally we can run the models for BMI.

models <- lapply(variables, function(v){
    model <- reformulate(c("age", "sex", "bmi"), v) %>% glm(data=data)
    attr(model, "feature") <- v
    model
})

bmi_tests <- lapply(models, function(m) {
    coefs <- summary(m)$coefficients
    data.table(
        feature=attr(m, "feature"), 
        coef=coefs[4, 1], 
        se=coefs[4, 2], 
        t=coefs[4, 3], 
        p=coefs[4, 4]
    )
}) %>% rbindlist()
bmi_tests[, "padj" := p.adjust(p, method="fdr")]
bmi_tests[, "variable" := "bmi"]
bmi_tests[order(p)]

And for weight loss.

models <- lapply(variables, function(v){
    model <- reformulate(c("age", "sex", "bmi", "subset"), v) %>% 
             glm(data=data)
    attr(model, "feature") <- v
    model
})

wl_tests <- lapply(models, function(m) {
    coefs <- summary(m)$coefficients
    data.table(
        feature=attr(m, "feature"), 
        coef=coefs[5, 1], 
        se=coefs[5, 2], 
        t=coefs[5, 3], 
        p=coefs[5, 4]
    )
}) %>% rbindlist()
wl_tests[, "padj" := p.adjust(p, method="fdr")]
wl_tests[, "variable" := "subset"]
wl_tests[order(p)]

tests <- rbind(bmi_tests, wl_tests)

library(ggplot2)
theme_minimal() %>% theme_set()

wide <- dcast(tests, feature ~ variable, value.var=c("t", "padj"), fill=0)
wide[, "sig" := "none"]
wide[padj_bmi < 0.05, "sig" := "BMI"]
wide[padj_subset < 0.05, "sig" := "weight loss"]
wide[padj_bmi < 0.05 & padj_subset < 0.05, "sig" := "both"]
wide[, "sig" := factor(sig, levels=c("none", "BMI", "weight loss", "both"))]
ggplot(wide, aes(x=t_bmi, y=t_subset, color=sig)) +
    geom_vline(xintercept = 0, color="gray30", lty="dashed") +
    geom_hline(yintercept = 0, color="gray30", lty="dashed") +
    geom_point() + stat_smooth(method="glm", aes(group = 1)) +
    labs(x = "t statistic BMI", y = "t statistic weight loss") +
    scale_color_manual(values=c(none = "black", BMI = "royalblue", `weight loss`="orange", both = "red"))

ggsave("figures/diet_t.png", dpi=300, width=5, height=3)
cor.test(~ t_bmi + t_subset, data=wide)

    Pearson's product-moment correlation

data:  t_bmi and t_subset
t = -1.1438, df = 37, p-value = 0.26
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4727561  0.1388057
sample estimates:
      cor 
-0.184805 

Overall explained variance

library(vegan)

m <- as(data[, variables, with=F], "matrix")
m[is.na(m)] <- min(as.numeric(m), na.rm=T)
perm <- adonis(m ~ age + sex + bmi + subset, data=data, method="euclidean")
perm

Call:
adonis(formula = m ~ age + sex + bmi + subset, data = data, method = "euclidean") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
age         1     324.9  324.89 1.96175 0.01883  0.124
sex         1      94.0   94.01 0.56765 0.00545  0.538
bmi         1     159.0  158.96 0.95982 0.00921  0.320
subset      1     118.8  118.76 0.71709 0.00688  0.447
Residuals 100   16561.3  165.61         0.95963       
Total     104   17257.9                 1.00000