How to perform a what-if study using observational data? 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.
 A: This is linked to causal inference. The potential outcomes/counterfactual model will be useful for conceptualization. One easy way to do this is to use the g-formula. The g-formula is essentially the same as do-calculus. For an introductory paper, I would recommend this article by Snowden et al. 2011
The key part of the application of the g-formula is the identifiability assumptions. Common ones are exchangeability, causal consistency, and positivity. This paper describes the conditions. Exchangeability is the assumption that there is no common cause of $x_j$ and $y$. This is the d-separation rule of do-calculus. Causal consistency implies that you have sufficiently specified $x_j$ so that there is no difference in $y$ by residual differences in $x_j$. You can think about this as assuming taking an aspirin in the morning versus night has no difference on mortality. Positivity is the assumption that all individuals have a non-zero probability of having all values of $x_j$. Westreich & Cole provide further details.
Having said that, I would be careful with the usage of machine learning for causal inference. In general, confidence intervals for the g-formula are calculated by bootstrapping. There has been a lot of work showing that confidence intervals based on bootstrapping procedures is invalid for machine learning algorithms, like random forests. For an example, see Naimi & Edwards.
Lastly, since you mentioned that $x_j$ is a continuous variable, you will need to consider how $x_j$ is changed. You can imagine that everyone's $x_j$ is shifted by some quantity $\delta$. Alternatively, you can set a threshold, $\gamma$, where those below (above) are shifted to $\gamma$. The remainder of individuals retain their $x_j$ value
