Partial least squares regression (PLSR) is a statistical technique that allows you to predict multiple response variables from multiple predictor variables. It works by essentially running separate principal components analyses (PCAs) on the matrix of the response variables and the matrix of the predictor variables, and then it regresses the equivalent axes from each. For example, the first PCA axis from the predictor variables is regressed with the first PCA axis from the response variables, the second PCA axis from the predictor variables is regressed with the second PCA axis from the response variables, and so on. It isn't exactly PCA that is being performed on the matrices of the predictor and the response variables, though, because the goal is not to capture the maximum variability within predictor and response matrices. Instead, the goal is to capture the maximum variability between the predictor and the response matrices. Thus, what I have been calling PCA axes may not actually be true PCA axes - they may be rotated to allow the subsequent regressions between the predictor and the response variable axes to be better.
Since PLSR is based on a similar technique as PCA, it assumes that there are linear relationships between variables. In practice, these relationships are not always linear.
Non-metric multidimensional scaling (NMDS) is a more robust technique than PCA as it does not assume relationships between variables are linear, but NMDS axis scores may be somewhat meaningless, especially when there are high stress values.
Is it reasonable to use a technique like NMDS to generate axis scores to use in PLSR regression? In other words, instead of using a technique like PCA to generate axis scores for points to be used in subsequent regression analysis (where we regress the first axis from the predictor variable ordination with the first axis of the response variable ordination), could we use something like NMDS and avoid making the assumption that relationships between data are linear?