Response variable: Continuous, Explanatory variables: One factor (8 levels) and one multivariate. How to model the response variable? I am a biologist and keep little insights into statistics, apologies, if the question seems trivial, thus please provide your suggestions.
Let me explain the situation completely. From same plant, I have data for nitrogen content (y, continuous and normally distributed) as well as their genotype (x1, factor with 8 levels) and their microbiome profile (Multivariate with many features; dimension reduction can be performed). If I want to find which of this two variables (x1 and x2) better explains my nitrogen content, how do I model it in a statistically sound manner?
To summarise: y~x1+x2 and I want to find which variable, x1 or x2 is more important?
I am particularly concerned on how to incorporate x2 into the model as this variable has many features, like we see in omics datasets.
I would appreciate any lead, thanks!

Consider the linear model
$$y_{i} = \mu + \sum \limits_{j=1}^7 \alpha_{j} I_{x_1=j}(x_{1i}) + \sum \limits_{k=1}^M  \beta_k x_{2ki} + \epsilon_{i}$$
$i$ is the number of data points and j=1-7 (categories of $x_1$).
Thus $\alpha_j$ is the coefficient for the $j$th category. Note that R will restrict $\alpha_1=0$ (be careful about the interpretation).
So $\mu$ is the intercept which is essentially the effect of the first category of $x_1$.
The indicator, $I_{x_1=j}(x_{1i})=1$ is if $x_{1i}$ is in category $j$(the ith data point in $x_1$).
Also this assumes that $X_2=(x_{21}, \dots, x_{2M})$ has $M$ variables.
Finally: And $\epsilon_i \sim N(0, \sigma^2)$ are the errors.
Also the two sub-models:
$$M_1: y_{i} = \mu + \sum \limits_{j=1}^7 \alpha_{j} I_{x_1=j}(x_{1i})+\epsilon_{i}$$
$$M_2: y_{i} = \mu  + \sum \limits_{k=1}^M  \beta_k x_{2ki} + \epsilon_{i}$$
Models explaining y by x1 only or X2 only.
Objective: you basically want to determine which group of coefficients ($\alpha_i$s vs $\beta_k$s) are more "important".
Here's the R call:
full_model = lm(y~x1 + x2_1 + x2_2 + .... + x2_M)
m1 = lm(y~x1)
m2 = lm(y~x2_1 + x2_2 + .... + x2_M)

First test for significance: You can test the significance of each set of coefficients using anova
anova(m1, full_model)
anova(m2, full_model)

which tests if all $\alpha_i$ or $\beta_k$ are likely to be zero or not. In other words, are $x_1$ or  $X_2$ explain $y$ at all?
This does so by performing an F test of the full model vs dropping $x_1$ or $X_2$.
Compare AIC/Adjusted $R^2$ : If both are significant, you can proceed to compare $M_1$ and $M_2$.
You can compare the AIC & Adjusted $R^2$ of these two models.
In code:
AIC(M1); AIC(M2) 


All of this assumes that the residual analysis confirming the core assumptions of linear regression. See this link for details.
