Is repeated propensity score matching over many 0-1-features a valid procedure? I would like to do a simple linear model where the outcome $y$ is real-valued, but my data matrix $X$ consists of several hundred features that all are $0$-$1$-valued. The number of observations $n$ is approximately of the order $10^5$, while the number of features $p$ is between $200$ and $300$.
The $0$-$1$-distributions are very different between features (column sums of $X$ range between 1 and several thousand) and some of them are correlated (positively or negatively) to different degrees. I did some pruning by only considering features with at least 30 (100, ...) $1$'s or cases.
My aim is to assign to every feature a number of how much it changes the value of the outcome $y$ (in units of $y$). As such, I'm not nailed down on regression coefficients $\beta$, but it seemed the most natural way to do. I tried (penalized) regression, but the results are so-so and sometimes uninterpretable because of the imbalances of the features.
Another idea is inspired by propensity score matching, where I think about every feature as a "treatment":
I would go through features one-by-one (let's fix feature $x_j$) and do a propensity score matching for the "treatment" $x_j$ on all other features $x_i$, $i\neq j$ (wihtout considering $y$). Then I have a matched sample for the $0$- and $1$-cases of $x_j$ and could analyze how much different the $y(1)$-values are compared to the $y(0)$-cases (sorry for abuse of notation). Afterwards I would go on to the next feature and repeat.
Is this a valid procedure? It could be a problem that I repeatedly use features for the conditional matching?
 A: You could do this, but the results you would get would not be any more meaningfully interpretable than regression coefficients. The interpretation of a treatment effect in a linear regression model that controls for covariates is the same as the interpretation of a treatment effect resulting from a propensity score analysis. This method would not buy you anything and is not what propensity score analysis was designed to do. In addition, treatment effect estimates from propensity score analysis have more uncertainty than those from regression, so if your regression results are unstable, your propensity score results would be even more so.
It might make sense to use a flexible machine learning technique that allows you to assess variable importance in a non-parametric way. I'm no expert on this so I can't make a more concrete recommendation than this. It's not immediately clear to me why penalized regression (e.g., elastic net) would not work well in this scenario. One thing I can think of is that because all the features are on the same scale, you should use the unstandardized coefficients in the loss function (rather than standardized coefficients, which is what's usually used in penalized regression to ensure the scaling of the variables doesn't impact their regression weights). But propensity score analysis is definitely not the right way to go about solving this problem.
