Suppose we are trying to model by linear regression $Y$ on $X$. However, we only observe $(Y,X)$ for $R = 1$. We know that $Pr(R|Y,X) = Pr(R|X)$, meaning that the missingness has nothing to do with the outcome. Why can we ignore the missing data in this case? Why will we have unbiased estimators?
I have seen many resources on econometrics claiming this without giving a proof.
In the context of regression analysis, all of our inferences are conditional on the explanatory variables. Consequently, if the missingness is random when conditioning on these explanatory variables (which is the case for MAR) this is sufficient to ensure that estimators retain the properties they would usually have under MCAR. The basic intuition for this is that if the missingness distribution only depends on the explanatory variables then it is essentially "fixed" when we condition on these variables.
The theory of missingness is usually shown in a more generalised setting than linear regression, and its main theorems are usually established directly with respect to the likelihood function in a highly general form of analysis. A good expository paper on the concept (with theorems and proofs) is Seaman et al (2013) (version without paywall here). Below I will adapt some material from this paper to show how MAR works in the context of regression analysis. This will show that ---under the assumption of MAR--- the likelihood function for the regression parameters is proportionate to the version that "ignores" the missingness mechanism.
Setting up the problem: Consider a problem involving a set of values $(\mathbf{x}_1, \mathbf{y}_1),...,(\mathbf{x}_n, \mathbf{y}_n)$ where the conditional distribution for $\mathbf{y}|\mathbf{x}$ uses the density $f_\theta$ which depends on a parameter $\theta$. Suppose we observe this process subject to MAR with missingness indicators $\mathbf{m}_1,...,\mathbf{m}_n$ that have a conditional distribution that uses the density $g_\phi$ that depends on a parameter $\phi$. Since the sampling and missingness determine the observed outcomes, there is a mapping $H$ from any vectors of explanatory and response variables and any vector of corresponding missingness indicators to the observed vectors of explanatory and response variables.
To facilitate our analysis, let $\mathbf{x}_*$ and $\mathbf{y}_*$ denote the actual vectors of explanatory and response variables and let $\mathbf{m}_*$ denote the missingness vector that leads to the observed values. Then let $\mathbf{x}_\text{obs}$ and $\mathbf{y}_\text{obs}$ denote the observed vectors of explanatory and response variables. These latter vectors are obtained by:
$$(\mathbf{x}_\text{obs}, \mathbf{y}_\text{obs}) = H(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*).$$
Without loss of generality, we will suppose that there are $k = n-\sum_i m_{*i}$ observed values in the observe vectors and $n-k = \sum_i m_{*i}$ missing values that do not appear in these vectors.
Likelihood analysis: If we were to ignore the missingness mechanism in the sampling method, and just treat this as if we sampled the observed outcomes at random without any missing data, our likelihood function for inferences about $\theta$ would be:
$$L_{\mathbf{x}_\text{obs}, \mathbf{y}_\text{obs}}(\theta) = \prod_{i=1}^k f_\theta(\mathbf{y}_{\text{obs}, i}|\mathbf{x}_{\text{obs}, i}).$$
However, with the missingness mechanism the actual likelihood function is:
$$\tilde{L}_{\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*}(\theta, \phi) = \int \limits_{\mathscr{H}(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*)} \bigg( \prod_{i=1}^n f_\theta(\mathbf{y}_i|\mathbf{x}_{*i})^{1-m_{*i}} \bigg) \cdot g_\phi(\mathbf{m}_*|\mathbf{x}_*, \mathbf{y}) \ d \mathbf{y},$$
where we define the space $\mathscr{H}(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*) \equiv \{ \mathbf{y} | H(\mathbf{x}_*, \mathbf{y}, \mathbf{m}_*) = H(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*) \}$.
Now, obviously these are very different likelihood functions. However, if missingness-at-random (MAR) holds then the distribution of the missingness vector depends on $(\mathbf{x},\mathbf{y})$ only through $\mathbf{x}$. Under this assumption we have $g_\phi(\mathbf{m}|\mathbf{x}, \mathbf{y}) = g_\phi(\mathbf{m}|\mathbf{x})$ for all $\mathbf{x},\mathbf{y},\mathbf{m}$, which then gives:
$$\begin{align} \tilde{L}_{\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*}(\theta, \phi) &= \int \limits_{\mathscr{H}(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*)} \bigg( \prod_{i=1}^n f_\theta(\mathbf{y}_i|\mathbf{x}_{*i})^{1-m_{*i}} \bigg) \cdot g_\phi(\mathbf{m}_*|\mathbf{x}_*, \mathbf{y}) \ d \mathbf{y} \\[6pt] &= g_\phi(\mathbf{m}_*|\mathbf{x}_*) \int \limits_{\mathscr{H}(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*)}\bigg( \prod_{i=1}^n f_\theta(\mathbf{y}_i|\mathbf{x}_{*i})^{1-m_{*i}} \bigg) \ d \mathbf{y} \\[6pt] &= g_\phi(\mathbf{m}_*|\mathbf{x}_*) \int \limits_{\mathscr{H}(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*)}\bigg( \prod_{i=1}^k f_\theta(\mathbf{y}_{\text{obs}, i}|\mathbf{x}_{\text{obs}, i}) \bigg) \ d \mathbf{y} \\[6pt] &= g_\phi(\mathbf{m}_*|\mathbf{x}) \cdot \bigg( \prod_{i=1}^k f_\theta(\mathbf{y}_{\text{obs}, i}|\mathbf{x}_{\text{obs}, i}) \bigg) \cdot |\mathscr{H}(\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*)| \\[6pt] &\propto g_\phi(\mathbf{m}_*|\mathbf{x}) \cdot \bigg( \prod_{i=1}^k f_\theta(\mathbf{y}_{\text{obs}, i}|\mathbf{x}_{\text{obs}, i}) \bigg) \\[12pt] &= g_\phi(\mathbf{m}_*|\mathbf{x}) \cdot L_{\mathbf{x}_\text{obs}, \mathbf{y}_\text{obs}}(\theta). \\[6pt] \end{align}$$
If we are only concerned with making inferences about $\theta$ (i.e., if $\phi$ s a nuisance parameter) then we can make use of the fact that:
$$\tilde{L}_{\mathbf{x}_*, \mathbf{y}_*, \mathbf{m}_*}(\theta, \phi) \overset{\theta}{\propto} L_{\mathbf{x}_\text{obs}, \mathbf{y}_\text{obs}}(\theta).$$
As you can see, this means that the true likelihood/profile likelihood function for the paramater $\theta$ (the one that does not ignore the missingness mechanism) is proportionate to the sipler likelihood function obtained by ignoring the missingness mechanism. If you look through the mathematics above you will see that this occurs because the equation $g_\phi(\mathbf{m}|\mathbf{x}, \mathbf{y}) = g_\phi(\mathbf{m}|\mathbf{x})$ allows us to "separate" this part from the part that depends on $\theta$ ---i.e., the MAR assumption makes the true likelihood function "seperable" with respect to the parameters $\theta$ and $\phi$.
I really like the following figure from Daniel, Kenward, Cousens, & De Stavola (2012) that gives an intuitive view of why Missing At Random (MAR) given X is not a problem, but MAR given Y is a problem.
(NOTE that this figure uses the symbol A instead of X for the predictor.)
For the sake of illustration this figure shows an exaggerated situation where the missingness is either 0% or 100%, depending on the value of an observed variable. But the same principles hold in more realistic cases where missingness depends in a more fuzzy way on the observed variables.
The top panel shows MAR given X. In this situation, the efficiency/precision of the slope estimate is certainly lower, but the average estimate is the same as it would be if there were no missing data.
The bottom panel shows MAR given Y. Here the average slope estimate is not the same as it would be if there were no missing data -- that is, it's biased. That's because the missing large Y values lead the OLS slope estimates to be too shallow.
This figure can also be used to illustrate why even certain kinds of Missing Not At Random (MNAR) are not a problem either, although I won't get into that right now.
If you want to understand this better, you may be interested in this old blog post I wrote on the topic:
Maybe a couple of examples will help convince you.
(1) Suppose I have a random sample of size $m = 50$ from the distribution $\mathsf{Binom}(n = 10, p = .8)$ with $\mu = 8.$
set.seed(2022) # for reproducibility
x = rbinom(50, 10, 0.8); x
[1] 7 8 9 8 9 8 10 10 9 7
[11] 10 9 9 8 8 9 7 8 8 7
[21] 7 8 8 9 9 6 8 4 8 7
[31] 8 8 9 8 9 10 7 8 7 7
[41] 8 8 8 8 7 10 6 8 9 7
The sample mean is $\bar X = 8.04,$ reasonably close to $\mu = 8.$
mean(x)
[1] 8.04
Now, I sample indexes of four values at random to be 'missing'. The respective observations are $7, 9, 8, 7.$
sample(1:50, 4)
20 5 32 37
The mean of the remaining (non-missing) 46 observations is $8.0652,$ which is not much different from the mean $8.04$ of all 50 observations.
(sum(x) - 7 - 9 - 8 - 7)/46
[1] 8.065217
This is really too small a sample for a guaranteed good example, especially with a relatively large proportion (4 of 50) missing, but it worked pretty well.
(2) Suppose I have 1000 observations from $\mathsf{Norm}(100, 15)$ and randomly get rid of 50 observations.
This time I will randomize the order of the observations and leave off the last 50, because it would be impractical to show 1000 observations.
I'll also show that a few observations missing at random doesn't change the standard deviation by much.
set.seed(1234)
y = sample(rnorm(1000, 100, 15)) # scranble
mean(y); mean(y[1:950])
[1] 99.60104 # mean of all 1000
[1] 99.71998 # mean with 50 missing @ rand
sd(y); sd(y[1:950])
[1] 14.96007 # SD of all 1000
[1] 14.86055 # SD with 50 missing @ rand
Boxplots of the entire sample (bottom) and of the sample with 50 missing observations look much the same.
boxplot(y, y[1:950], horizontal=T, col="skyblue2")
Note: By contrast, if I delete observations above 124, instead of deleting at random, then it makes a larger difference in means, SDs, and the appearance of the boxplot. The at random part of "missing at random" is crucial.
mean(y); mean(y[y<124])
[1] 99.60104
[1] 97.6231
sd(y); sd(y[y<124])
[1] 14.96007
[1] 13.17083
