A team fit a Random Forest model to a dataset $S=\{\mathbf{x}_i,y_i\}_{i=1}^N$, where $\mathbf{x}$ is a vector of continuous and categorical variables, and $y$ is a binary response. The model has a low CV-error, and the precision-recall curve looks good.
Now, for new samples they would like to perform what-if studies, i.e. for a fixed variable $x_j$ they would like to see how $p(y=1|x_j)$ changes as a function of $x_j$, all other things being equal. The goal of this activity is the following: for new samples $\mathbf{x}$, they want to modify the value of $x_j$ to $x'_j$ , so that $x'_jp(y=1|(x_1,\dots,x'_j,\dots,x_n))$ is maximized. However, by doing so, it seems to me that they modify the distribution $p(\mathbf{x},y)$, thus the results are not reliable anymore. I think this is related to causal inference. Is there a way to modify this process so that it actually works? Which methods should I study to help this team?
EDIT: An important point I forgot to mention is that $x_j$ is one of the continuous variables, unfortunately.