You should have one row per measurement per individual. So, if for each of 100 individuals you took 790 measurements twice, you should have a dataset that is 100*790*2 rows long. Your columns should be the participant ID, the metabolite type, and your other variables. Because each participant is measured multiple times, participant-level characteristics like age will be duplicated across rows (i.e., you will have multiple rows with the same age because they all correspond to the same individual).

To get the data in this format, you'll need to turn your dataset from wide to long. There are a variety of tools in R to do that, but I think the `melt()` function in the `data.table` package is fairly easy to use. Your code will look something like the following, assuming the only variables in your dataset are the 790 metabolite columns, participant ID, and the variables you've included in the model:

    data_long <- melt(data, measure.vars = 1:790, variable.name = "met_type",
                      value.name = "expr_val")
         
This long dataset (which I've called `data_long`) should have a different structure from your current dataset. Rather than a column for each type of metabolite, there will be one column for the metabolite expression value (`expr_val`) and another for the type of metabolite (`met_type`). You should also have a column, which I've called `participantID` in the model, that identifies each participant (which you probably already have).

From here, how you specify the model will depend on your theories about the 790 metabolites. If each type of metabolite is important, you can add metabolite as a fixed, measurement-level covariate, which will be split into 789 dummy variables. If the specific type of the 790 metabolites is not important to your theory but needs to be accounted for, you can add a random effect for metabolite, which gives your model two random effects (one for participant and one for metabolite type). The model, with two random effects, should look like the following:

    fit <- lmerTest::lmer(expr_val ~ Diseasestatus + BB + WA + ACOG + Age + Adiposity + 
                           (1|participantID) + (1|met_type), data = data_long)

If you run `summary()` on the `fit` object, you will get the coefficient estimates and estimates of the variance components. If you want to get p-values to test whether the coefficient of `Diseasestatus` differs from zero, you need to use the `lmerTest` package.  The coefficient on `Diseasestatus` represents the difference between the average metabolite expression value for those with the disease and that for those without.


----------
**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.