Return to Answer

fit_list <- setNames(lapply(paste0("Met", 1:790), function(i) {
    fit <- lmer(data_wide[[i]] ~ Diseasestatus + BB + WA + ACOG + 
                   Age + Adiposity + (1|participantID), 
                data = data_wide)
    summary(model)$coef["Diseasestatus", c("Estimate", "Pr(>|t|)")]
}), paste0("Met", 1:790))

dplyr::bind_rows(fit_list, .id = "Metabolite")

fit_list <- setNames(lapply(paste0("Met", 1:790), function(i) {
    fit <- lmer(data_wide[[i]] ~ Diseasestatus + BB + WA + ACOG + 
                   Age + Adiposity + (1|participantID), 
                data = data_wide)
    summary(fit)$coef["Diseasestatus", c("Estimate", "Pr(>|t|)")]
}), paste0("Met", 1:790))

dplyr::bind_rows(fit_list, .id = "Metabolite")

fit_list <- setNames(lapply(1:790, function(i) {
    fit <- lmer(data_wide[[i]] ~ Diseasestatus + BB + WA + ACOG + 
                   Age + Adiposity + (1|participantID), 
                data = data_wide)
    summary(model)$coef["Diseasestatus", c("Estimate", "Pr(>|t|)")]
}), names(data_wide)[1:790])

dplyr::bind_rows(fit_list, .id = "Metabolite")

fit_list <- setNames(lapply(paste0("Met", 1:790), function(i) {
    fit <- lmer(data_wide[[i]] ~ Diseasestatus + BB + WA + ACOG + 
                   Age + Adiposity + (1|participantID), 
                data = data_wide)
    summary(model)$coef["Diseasestatus", c("Estimate", "Pr(>|t|)")]
}), paste0("Met", 1:790))

dplyr::bind_rows(fit_list, .id = "Metabolite")

EDIT:

It looks like you want an estimate of the effect for each metabolite. It makes sense to run a separate model for each metabolite, though after doing so you should correct for multiple comparisons. You need to loop through each metabolite and fit a separate model to each one. You could do this with either the long or the wide dataset. We'll stick with the wide dataset since it's smaller.

fit_list <- setNames(lapply(1:790, function(i) {
    fit <- lmer(data_wide[[i]] ~ Diseasestatus + BB + WA + ACOG + 
                   Age + Adiposity + (1|participantID), 
                data = data_wide)
    summary(model)$coef["Diseasestatus", c("Estimate", "Pr(>|t|)")]
}), names(data_wide)[1:790])

dplyr::bind_rows(fit_list, .id = "Metabolite")

To run this model you must have the lmerTest package loaded. You didn't see any p-values previously because you didn't have it loaded before running the model. The p-values from these tests use the Satterthwaite degrees of freedom.