How high must logistic covariates' predictive accuracy be for a reversal effect to show up?

Question

I am modeling an outcome for hospital patients, 'RA' (whether readmitted; 0=No, 1=Yes). My predictor of interest is 'HHS' (whether referred to Home Health Services such as from a visiting nurse; 0=No, 1=Yes). Those referred readmit at a 15.2% rate; others, 9.2%, but the former are needier, sicker patients. Conventional thinking is that if we controlled for severity of illness this difference would not only be washed out but would reverse itself. In other words, holding constant the severity of illness, having HHS should mean a lower RA rate.

With HHS as the sole predictor, its coefficient (B) in a logistic regression = 0.6 (N ~ 25k). B is reduced to 0.2 with a group of covariates controlled, each accounting for some aspect of severity of illness, but B doesn't fall below zero.

HHS alone explains only about 1% of the variance in RA; with the other predictors, this becomes 4%.* Perhaps this is the problem--that these covariates are not explaining enough variance to "succeed" in reversing the sign of the coefficient of interest. If this is true, is there a way to estimate how high their explained variance needs to be for such a reversal to show up?

EDIT: Alecos Papadopoulos has come up with an impressive solution that answers this question, soon to be published in The American Statistician. See https://www.tandfonline.com/doi/full/10.1080/00031305.2019.1704873

---

*Using either of 2 pseudo-R-squared formulas; Cox & Snell's or Menard's [-2LL0 - (-2LL1)] / [-2LL0.]

Answer 1

(This answer uses results from W.H. Greene (2003), Econometric Analysis, 5th ed. ch.21)

I will answer the following modified version, which I believe accomplishes the goals of the OP's question : "If we only estimate a logit model with one binary regressor of interest and some (dummy or continuous) control variables, can we tell whether dropping the control variables will result in a change of sign for the (coefficient of) the regressor of interest?"

Notation: Let $RA\equiv Y$ be the dependent variable, $HHS \equiv X$ the binary regressor of interest and $\mathbf Z$ a matrix of control variables. The size of the sample is $n$. Denote $n_0$ the number of zero-realizations of $X$ and $n_1$ the number of non-zero realizations, $n_0+n_1=n$. Denote $\Lambda()$ the cdf of the logistic distribution.
Let the model including the control variables (the "unrestricted" model) be

$$M_U : \begin{align} &P(Y=1\mid X,\mathbf Z)=\Lambda(X, \mathbf Z,b,\mathbf c)\\ &P(Y=0\mid X,\mathbf Z)=1-\Lambda(X, \mathbf Z,b,\mathbf c) \end{align}$$

where $b$ is the coefficient on the regressor of interest.
Let the model including only the regressor of interest (the "restricted" model) be

$$M_R : \begin{align} &P(Y=1\mid X)=\Lambda(X, \beta)\\ &P(Y=0\mid X)=1-\Lambda(X,\beta) \end{align}$$

STEP 1

Consider the unrestricted model. The first-derivative of the log-likelihood w.r.t to $b$ and the condition for a maximum is

$$\frac {\partial \ln L_U}{\partial b}= \sum_{i=1}^n\left[(y_i-\Lambda_i(x_i, \mathbf z_i,b,\mathbf c)\right]x_i=0 \Rightarrow b^*: \sum_{i=1}^ny_ix_i=\sum_{i=1}^n\Lambda_i(x_i, \mathbf z_i,b^*,\mathbf c^*)x_i \;[1]$$

The analogous relations for the restricted model is $$\frac {\partial \ln L_R}{\partial \beta}= \sum_{i=1}^n\left[(y_i-\Lambda_i(x_i,\beta)\right]x_i=0 \Rightarrow \beta^*: \sum_{i=1}^ny_ix_i=\sum_{i=1}^n\Lambda_i(x_i, \beta^*)x_i \qquad[2]$$

We have

$$\Lambda_i(X,\beta^*) = \frac {1}{1+e^{-x_i\beta^*}}$$

and since $X$ is a zero/one binary variable relation $[2]$ can be written

$$\beta^*: \sum_{i=1}^ny_ix_i=\frac {n_1}{1+e^{-\beta^*}} \qquad[2a]$$

Combining $[1]$ and $[2a]$ and using again the fact that $X$ is binary we obtain the following equality relation between the estimated coefficients of the two models:

$$\frac {n_1}{1+e^{-\beta^*}} = \sum_{i=1}^n\Lambda_i(x_i, \mathbf z_i,b^*,\mathbf c^*)x_i $$

$$\Rightarrow \frac {1}{1+e^{-\beta^*}} = \frac {1}{n_1}\sum_{x_i=1}\Lambda_i(x_i=1, \mathbf z_i,b^*,\mathbf c^*) \qquad [3]$$

$$\Rightarrow \hat P_R(Y=1\mid X=1) = \hat {\bar P_U}(Y=1\mid X=1,\mathbf Z) \qquad [3a]$$

or in words, that the estimated probability from the restricted model will equal the restricted average estimated probability from the model that includes the control variables.

STEP 2
For a sole binary regressor in a logistic regression, its marginal effect $m_R(X)$ is

$$ \hat m_R(X)= \hat P_R(Y=1\mid X=1) - \hat P_R(Y=1\mid X=0)$$

$$ \Rightarrow \hat m_R(X) = \frac {1}{1+e^{-\beta^*}} - \frac 12$$

and using $[3]$

$$ \hat m_R(X) = \frac {1}{n_1}\sum_{x_i=1}\Lambda_i(x_i=1, \mathbf z_i,b^*,\mathbf c^*) - \frac 12 \qquad [4]$$

For the unrestricted model that includes the control variables we have

$$ \hat m_U(X)= \hat P_U(Y=1\mid X=1, \bar {\mathbf z}) - \hat P_U(Y=1\mid X=0, \bar {\mathbf z})$$

$$\Rightarrow \hat m_U(X) = \frac {1}{1+e^{-b^*-\bar {\mathbf z}'\mathbf c^*}} - \frac {1}{1+e^{-\bar {\mathbf z}'\mathbf c^*}} \qquad [5]$$

where $\bar {\mathbf z}$ contains the sample means of the control variables.

It is easy to see that the marginal effect of $X$ has the same sign as its estimated coefficient. Since we have expressed the marginal effect of $X$ from both models in terms of the estimated coefficients from the unrestricted model, we can estimated only the latter, and then calculate the above two expressions ($[4]$ and $[5]$) which will tell us whether we will observe a sign reversal for the coefficient of $X$ or not, without the need to estimate the restricted model.

Answer 2

This is for OLS regression. Consider a geometric representation of three variables -- two predictors, $X_1$ and $X_2$, and a dependent variable, $Y$. Each variable is represented by a vector from the origin. The length of the vector equals the standard deviation of the corresponding variable. The cosine of the angle between any two vectors equals the correlation of the corresponding two variables. I will take all the standard deviations to be 1.

The picture shows the plane determined by the $X_1$ and $X_2$ when they correlate positively with one another. $Y$ is a vector coming out of the screen; the dashed line is its projection into the predictor space and is the regression estimate of $Y$, $\hat{Y}$. The length of the dashed line equals the multiple correlation, $R$, of $Y$ with $X_1$ and $X_2$.

If the projection is in any of the colored sectors then both predictors correlate positively with $Y$. The signs of the regression coefficients $\beta_1$ and $\beta_2$ are immediately apparent visually, because $\hat{Y}$ is the vector sum of $\beta_1 X_1$ and $\beta_2 X_2$. If the projection is in the yellow sector then both $\beta_1$ and $\beta_2$ are positive, but if the projection is in either the red or the blue sector then we have what appears to be suppression; that is, the sign of one of the regression weights is opposite to the sign of the corresponding simple correlation with $Y$. In the picture, $\beta_1$ is positive and $\beta_2$ is negative.

Since the length of the projection can vary between 0 and 1 no matter where it is in the predictor space, there is no minimum $R^2$ for suppression.

Answer 3

There is no obvious relationship between $R^2$ and reversal of the sign of a regression coefficient. Assume you have data for which the true model is for example $$ y_i = 0+5x_i -z_z + \epsilon_i $$ with $\epsilon_i \sim N(0, sd_\text{error}^2)$. I show the zero to make explicit that the intercept of the true model is zero, this is just a simplification.

When x and y are highly correlated and centered about zero then the coefficient of z when regressing over just z will be positive instead of negative. Note that the true model coefficients do not change with $sd_\text{error}$ but you can make $R^2$ vary between zero and one by changing the magnitude of the residual error. Look for example at the following R-code:

require(MASS)
sd.error <- 1
x.and.z <- mvrnorm(1000, c(0,0) , matrix(c(1, 0.9,0.9,1),nrow=2)) # set correlation to 0.9
x <- x.and.z[, 1]
z <- x.and.z[, 2]
y <- 5*x - z + rnorm(1000, 0, sd.error) # true model
modell1 <- lm(y~x+z)
modell2 <- lm(y~z)
print(summary(modell1)) # coefficient of z should be negative
print(summary(modell2)) # coefficient of z should be positive   

and play a bit with sd.error. Look for example at $sd_\text{error}=50$.

Note that with a very large sd.error the coefficient estimation will become more unstable and the reversal might not show up every time. But that's a limitation of the sample size.

A short summary would be that the variance of the error does not affect the expectations and thus reversal. Therefore neither does $R^2$.