Do I need to use multivariate regression or several regression analyses? I have a data set of 45 participants with 96 variables each (although some measurements are missing). Some variables are simple such as age and disability while other measurements are scores on some given test (e.g. one test has 5 values as an outcome). I have data of 5 tests, given on 3 different points in time and as stated earlier sometimes multiple scores per test.
Since the dataset is so large (given the amount of features wrt the amount of participants), I decided to predict the scores on a test given all previous results (such as age, disability, and all the scores on the same previous test). So this basically boils down to that I want to predict 5 features given approximately 10 features in 45 participants using regression (I wish to view the exact coefficients, p-values and R squared measurements).
Should I do a regular regression on each of the features I wish to predict, or should I use multivariate regression on all of the features I wish to predict at once? What is the difference?
 A: Let $Y_i$ denote the vector of $i$th response, wehre $i = 1, \dots, r$. In your example $r$ is 5 since you have 5 test scores. Let $X$ be an $n \times p$ matrix of predictors. If you implement $r$ separate regressions (one for each $Y_i$),
$$Y_i = X\beta_i + \epsilon_i, $$
where $\epsilon_i \sim N_n(0, \sigma^2_iI_n)$. Using OLS, you get estimates for $\beta$. You can also do a multivariate regression,
$$\mathbf{Y} = X\mathbf{B} + \mathbf{E}, $$
where $\mathbf{Y}$ is the $n \times r$ matrix of responses, $\mathbf{B}$ is the $p \times r$ matrix of regression coefficients, and $\mathbf{E}$ is the error matrix such that the $i$th column, $\epsilon_i \overset{iid}{\sim} N_n(0, \sigma^2_iI_n)$. In this case, the OLS estimate for $\mathbf{B}$ is equivalent to the $r$ OLS estimates for $\beta_i$.
However, if you have reason to assume that conditioned on $X$, the 5 predictors are correlated (which seems like that would be a reasonable assumption in your case), then the rows of $\mathbf{E}$ are assumed to to be such that for $j = 1, 2, \dots, n, \epsilon_j \overset{iid}{\sim} N_r(0, \Sigma)$. Here $\Sigma$ now represents the correlation structure for the predictors as well.
It is important to note that even in this case, the estimate for $\mathbf{B}$ is the same as the OLS estimate, but the error structure of the estimates changes, and thus inference on the estimates changes. As a consequence, $p$-values change.
The MRCE R package allows for such model fits. This package also uses regularization methods for when $n$ is not large enough compared to $p$, so you might not be forced to reduce to a smaller number of predictors. You can also find more detailed theory here along with motivating examples. The authors state the following motivation

Applications of this general model arise in chemometrics, econometrics, psychometrics, and other quantitative disciplines where one predicts multiple responses with a single set of prediction variables. For example, predicting several measures of quality of paper with a set of variables relating to its production.

Similarly in your setup, you seem to have 5 sets of responses arising from the same predictors, with an inherent correlation structure between responses.
