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

subcohort <- fread("data/manifest.csv", colClasses = classes)[has_close_microbiome == T]
subcohort[, "public_client_id" := paste0("0", public_client_id)]

cohort[, "additional_assays" := FALSE]
cohort[public_client_id %chin% subcohort$public_client_id, "additional_assays" := TRUE]
library(arivale.data.interface)
chem <- get_snapshot("chemistries", clean=T)
chem <- chem[cohort, on=c("public_client_id", "days_in_program"), roll="nearest"]
tt <- cohort[, t.test(bmi[since_baseline > 0], bmi[since_baseline == 0], paired=TRUE)]

ggplot(cohort, aes(x=since_baseline, y=bmi, group=public_client_id, 
                   fill=additional_assays, color=subset)) +
    geom_line(alpha=0.5) +
    geom_point(shape=21, stroke=0, size=2.5) +
    scale_fill_manual(values=c("black", "mediumseagreen")) +
    labs(x="time in program [days]", y="BMI [kg/m²]") +
    guides(color="none", fill="none")

ggsave("figures/bmi_trajectories.png", width=4, height=3, dpi=300)

cohort[, t.test(bmi[since_baseline == 0], bmi[since_baseline > 0], paired=TRUE)]

    Paired t-test

data:  bmi[since_baseline == 0] and bmi[since_baseline > 0]
t = 5.2091, df = 104, p-value = 9.647e-07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.655356 1.461042
sample estimates:
mean of the differences 
               1.058199 
measures <- c("ADIPONECTIN__SERUM", "HDL_CHOL_DIRECT", "LDL_CHOL_CALCULATION", 
              "GLUCOSE", "HOMA_IR", "INSULIN", "CRP_HIGH_SENSITIVITY", 
              "GLYCOHEMOGLOBIN_A1C")
legends <- c("serum adiponectin [mg/L]", "serum HDL [mg/L]", "serum LDL [mg/L]", 
             "serum glucose [mg/L]", "HOMA IR [index]", "serum insulin [mIU/L]", 
             "serum CRP [mg/L]", "glycated hemoglobin [%]")
select <- chem[, c("public_client_id", "since_baseline", "additional_assays", measures), with=F]
names(select) <- c("public_client_id", "since_baseline", "additional_assays", legends)
select <- melt(select, id.vars=c("public_client_id", "since_baseline", "additional_assays"), value.name="value", variable.name="measure")

stats <- function(value, since_baseline) {
    s <- t.test(value[since_baseline > 0], value[since_baseline == 0], paired=TRUE)
    list(t=s$statistic, p=s$p.value, delta_mu=s$estimate)
}

res <- select[, stats(value, since_baseline), by="measure"]
res[, "label" := sprintf("Δ=%.2g, p=%.2g", delta_mu, p)]

pl <- ggplot(select, aes(x=since_baseline, y=value, group=public_client_id, color=additional_assays)) +
    geom_line(alpha=0.5, color="seagreen4") +
    geom_point() +
    labs(x="time in program [days]", y="") +
    facet_wrap(~ measure, scales="free_y") +
    guides(color="none") +
    scale_color_manual(values=c("black", "tomato"))
pl <- pl + geom_text(data = res, 
                     mapping = aes(x = Inf, y = Inf, label = label, group=1),
                     color = "black",
                     hjust = 1,
                     vjust = 2.5)
ggsave("figures/other_trajectories.png", width=8, height=8, dpi=300)
pl

LS0tCnRpdGxlOiAiQ29ob3J0IGNoYXJhY3RlcmlzdGljcyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkobWJ0b29scykKdGhlbWVfc2V0KHRoZW1lX21pbmltYWwoKSkKYGBgCgpgYGB7cn0KY29ob3J0IDwtIHJiaW5kKAogICAgZnJlYWQoIm5vX3dlaWdodF9sb3NzLmNzdiIsIGNvbENsYXNzZXM9YyhwdWJsaWNfY2xpZW50X2lkPSJjaGFyYWN0ZXIiKSksCiAgICBmcmVhZCgic3VjY2Vzc2Z1bF93ZWlnaHRfbG9zcy5jc3YiLCBjb2xDbGFzc2VzPWMocHVibGljX2NsaWVudF9pZD0iY2hhcmFjdGVyIikpCikKY29ob3J0WywgInN1YnNldCIgOj0gImNvbnRyb2xzIl0KY29ob3J0W3dlaWdodF9jaGFuZ2VfcmVsYXRpdmUgPCAwLCAic3Vic2V0IiA6PSAid2VpZ2h0IGxvc3MiXQpjb2hvcnRbLCBzdWJzZXQgOj0gZmFjdG9yKHN1YnNldCldCgpzdWJjb2hvcnQgPC0gZnJlYWQoImRhdGEvbWFuaWZlc3QuY3N2IiwgY29sQ2xhc3NlcyA9IGNsYXNzZXMpW2hhc19jbG9zZV9taWNyb2Jpb21lID09IFRdCnN1YmNvaG9ydFssICJwdWJsaWNfY2xpZW50X2lkIiA6PSBwYXN0ZTAoIjAiLCBwdWJsaWNfY2xpZW50X2lkKV0KCmNvaG9ydFssICJhZGRpdGlvbmFsX2Fzc2F5cyIgOj0gRkFMU0VdCmNvaG9ydFtwdWJsaWNfY2xpZW50X2lkICVjaGluJSBzdWJjb2hvcnQkcHVibGljX2NsaWVudF9pZCwgImFkZGl0aW9uYWxfYXNzYXlzIiA6PSBUUlVFXQpgYGAKCmBgYHtyfQpsaWJyYXJ5KGFyaXZhbGUuZGF0YS5pbnRlcmZhY2UpCmNoZW0gPC0gZ2V0X3NuYXBzaG90KCJjaGVtaXN0cmllcyIsIGNsZWFuPVQpCmNoZW0gPC0gY2hlbVtjb2hvcnQsIG9uPWMoInB1YmxpY19jbGllbnRfaWQiLCAiZGF5c19pbl9wcm9ncmFtIiksIHJvbGw9Im5lYXJlc3QiXQpgYGAKCmBgYHtyLCBmaWcud2lkdGg9NCwgZmlnLmhlaWdodD0zfQp0dCA8LSBjb2hvcnRbLCB0LnRlc3QoYm1pW3NpbmNlX2Jhc2VsaW5lID4gMF0sIGJtaVtzaW5jZV9iYXNlbGluZSA9PSAwXSwgcGFpcmVkPVRSVUUpXQoKZ2dwbG90KGNvaG9ydCwgYWVzKHg9c2luY2VfYmFzZWxpbmUsIHk9Ym1pLCBncm91cD1wdWJsaWNfY2xpZW50X2lkLCAKICAgICAgICAgICAgICAgICAgIGZpbGw9YWRkaXRpb25hbF9hc3NheXMsIGNvbG9yPXN1YnNldCkpICsKICAgIGdlb21fbGluZShhbHBoYT0wLjUpICsKICAgIGdlb21fcG9pbnQoc2hhcGU9MjEsIHN0cm9rZT0wLCBzaXplPTIuNSkgKwogICAgc2NhbGVfZmlsbF9tYW51YWwodmFsdWVzPWMoImJsYWNrIiwgIm1lZGl1bXNlYWdyZWVuIikpICsKICAgIGxhYnMoeD0idGltZSBpbiBwcm9ncmFtIFtkYXlzXSIsIHk9IkJNSSBba2cvbcKyXSIpICsKICAgIGd1aWRlcyhjb2xvcj0ibm9uZSIsIGZpbGw9Im5vbmUiKQpnZ3NhdmUoImZpZ3VyZXMvYm1pX3RyYWplY3Rvcmllcy5wbmciLCB3aWR0aD00LCBoZWlnaHQ9MywgZHBpPTMwMCkKCmNvaG9ydFssIHQudGVzdChibWlbc2luY2VfYmFzZWxpbmUgPT0gMF0sIGJtaVtzaW5jZV9iYXNlbGluZSA+IDBdLCBwYWlyZWQ9VFJVRSldCmBgYAoKYGBge3IsIGZpZy53aWR0aD04LCBmaWcuaGVpZ2h0PTh9Cm1lYXN1cmVzIDwtIGMoIkFESVBPTkVDVElOX19TRVJVTSIsICJIRExfQ0hPTF9ESVJFQ1QiLCAiTERMX0NIT0xfQ0FMQ1VMQVRJT04iLCAKICAgICAgICAgICAgICAiR0xVQ09TRSIsICJIT01BX0lSIiwgIklOU1VMSU4iLCAiQ1JQX0hJR0hfU0VOU0lUSVZJVFkiLCAKICAgICAgICAgICAgICAiR0xZQ09IRU1PR0xPQklOX0ExQyIpCmxlZ2VuZHMgPC0gYygic2VydW0gYWRpcG9uZWN0aW4gW21nL0xdIiwgInNlcnVtIEhETCBbbWcvTF0iLCAic2VydW0gTERMIFttZy9MXSIsIAogICAgICAgICAgICAgInNlcnVtIGdsdWNvc2UgW21nL0xdIiwgIkhPTUEgSVIgW2luZGV4XSIsICJzZXJ1bSBpbnN1bGluIFttSVUvTF0iLCAKICAgICAgICAgICAgICJzZXJ1bSBDUlAgW21nL0xdIiwgImdseWNhdGVkIGhlbW9nbG9iaW4gWyVdIikKc2VsZWN0IDwtIGNoZW1bLCBjKCJwdWJsaWNfY2xpZW50X2lkIiwgInNpbmNlX2Jhc2VsaW5lIiwgImFkZGl0aW9uYWxfYXNzYXlzIiwgbWVhc3VyZXMpLCB3aXRoPUZdCm5hbWVzKHNlbGVjdCkgPC0gYygicHVibGljX2NsaWVudF9pZCIsICJzaW5jZV9iYXNlbGluZSIsICJhZGRpdGlvbmFsX2Fzc2F5cyIsIGxlZ2VuZHMpCnNlbGVjdCA8LSBtZWx0KHNlbGVjdCwgaWQudmFycz1jKCJwdWJsaWNfY2xpZW50X2lkIiwgInNpbmNlX2Jhc2VsaW5lIiwgImFkZGl0aW9uYWxfYXNzYXlzIiksIHZhbHVlLm5hbWU9InZhbHVlIiwgdmFyaWFibGUubmFtZT0ibWVhc3VyZSIpCgpzdGF0cyA8LSBmdW5jdGlvbih2YWx1ZSwgc2luY2VfYmFzZWxpbmUpIHsKICAgIHMgPC0gdC50ZXN0KHZhbHVlW3NpbmNlX2Jhc2VsaW5lID4gMF0sIHZhbHVlW3NpbmNlX2Jhc2VsaW5lID09IDBdLCBwYWlyZWQ9VFJVRSkKICAgIGxpc3QodD1zJHN0YXRpc3RpYywgcD1zJHAudmFsdWUsIGRlbHRhX211PXMkZXN0aW1hdGUpCn0KCnJlcyA8LSBzZWxlY3RbLCBzdGF0cyh2YWx1ZSwgc2luY2VfYmFzZWxpbmUpLCBieT0ibWVhc3VyZSJdCnJlc1ssICJsYWJlbCIgOj0gc3ByaW50ZigizpQ9JS4yZywgcD0lLjJnIiwgZGVsdGFfbXUsIHApXQoKcGwgPC0gZ2dwbG90KHNlbGVjdCwgYWVzKHg9c2luY2VfYmFzZWxpbmUsIHk9dmFsdWUsIGdyb3VwPXB1YmxpY19jbGllbnRfaWQsIGNvbG9yPWFkZGl0aW9uYWxfYXNzYXlzKSkgKwogICAgZ2VvbV9saW5lKGFscGhhPTAuNSwgY29sb3I9InNlYWdyZWVuNCIpICsKICAgIGdlb21fcG9pbnQoKSArCiAgICBsYWJzKHg9InRpbWUgaW4gcHJvZ3JhbSBbZGF5c10iLCB5PSIiKSArCiAgICBmYWNldF93cmFwKH4gbWVhc3VyZSwgc2NhbGVzPSJmcmVlX3kiKSArCiAgICBndWlkZXMoY29sb3I9Im5vbmUiKSArCiAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzPWMoImJsYWNrIiwgInRvbWF0byIpKQpwbCA8LSBwbCArIGdlb21fdGV4dChkYXRhID0gcmVzLCAKICAgICAgICAgICAgICAgICAgICAgbWFwcGluZyA9IGFlcyh4ID0gSW5mLCB5ID0gSW5mLCBsYWJlbCA9IGxhYmVsLCBncm91cD0xKSwKICAgICAgICAgICAgICAgICAgICAgY29sb3IgPSAiYmxhY2siLAogICAgICAgICAgICAgICAgICAgICBoanVzdCA9IDEsCiAgICAgICAgICAgICAgICAgICAgIHZqdXN0ID0gMi41KQpnZ3NhdmUoImZpZ3VyZXMvb3RoZXJfdHJhamVjdG9yaWVzLnBuZyIsIHdpZHRoPTgsIGhlaWdodD04LCBkcGk9MzAwKQpwbApgYGAK