What variation are multiple fixed effects explaining? I wanted to share a thought experiment and get some reactions. I could not find this discussed in plain language. If I have a panel data set that has four variables a y, x, person1_id, and person2_id and I run the following regression, (y ~ x + person1_id + person2_id) what variation is the coefficient on x explaining? 
Imagine the data are the number of smiles (y) and person2's t-shirt color brightness (x). Is this regression explaining the variation in the number of smiles person1 gives person2 based on person2's t-shirt color brightness?
It is clear if we exclude the second set of dummies (person2_id) and have the following regression (y ~ x + person1_id) that the variation that this is explaining is how the same person (person1_id) smiles based on t-shirt color brightness. 
 A: Assumptions:


*

*Let assume, that observations are independent i.e. $Cov(\varepsilon_{i,j}, \varepsilon_{j,i}) = 0$ for $i \neq j$ and $y_{i,i}$ do not exist.
This is a strong assumption, but situation need somehow to be clarified. Are smiles mutually generated? How they are measured? Is it possible, that $\varepsilon_{i,j}$ and $\varepsilon_{j,i}$ are correlated? Maybe $y_{i,j} = y_{j,i}$? There are many possibilities, so let assume this one.

*Let assume then that precisely $y_{i,j}$ is the number of smiles person $i$ 'gives' to person $j$ based on the of brightness  of $j$'s shirt ($x_j$), the generosity of person $i$ of giving smiles away ($z_i$), and the $appeal$ of person $j$ ($v_j$). The observation is one-sided smilegiving fact.

*The true, full model is:
$y_{i,j} = \beta_0 + \beta_1 x_j + \beta_2v_j + \beta_3 z_i + \varepsilon_{i,j}$

*We can not observe $v_j$ and $z_i$, but we observe $x_j$. What are the consequences? The estimator may be inconsistent and biased. This situation happens only, when the variable $x_j$ is correlated with $v_j$, $z_i$ or both. 

*But hey, we have this smilegiving setting! We can think about it theoretically! How is it possible, that the generosity of person $i$ is correlated to brightness of person's $j$ shirt? Maybe it can somehow happen? How? We need another assumption - let assume that such thing do not happen, i.e. $x_j$ and $z_i$ are independent, and then uncorrelated.

*We can also think of second variable: is it possible that more applealing people wear brighter t-shirts? Maybe otherwise? Let assume, that there are two possibilities: independence and correlation.

Consequences of omitting of a variable:


*

*Instead of full model, we estimate model: 
$y_{i,j} = \beta_0 + \beta_1 x_j + \beta_2v_j + \varepsilon_{i,j}$
There are two consequences of such an action: First, nothing wrong happened to $\hat{\beta_1}$ - it is still consistent and unbiased. That is great thing in statistical (or causal) inference. Second, now variable $z_i$ is included in $\varepsilon_{i,j}$. The precision of prediction is lower, errors are bigger.

*Instead of full model we estimated model:
$y_{i,j} = \beta_0 + \beta_1 x_j + \varepsilon_{i,j}$
We assumed two possibilities: that $x_j$ and $v_j$ are independent, or that they are correlated. 


*

*If they are independent, then the situation is exactly the same as before: our estimator is unbiased and consistent, while prediction is less precise.

*What interests us, is other option: when the variables $x_j$ and $v_j$ are correlated. Now the estimator is biased and inconsistent. 


*Statistical (and dreamly causal) inference is dead in option 2, while it was perfectly fine in option 1.
The interpretation of $\hat{\beta_1}$ may be: When we observe change of $x_j$ by one unit, then we observe change in $y_{i,j}$ by $\hat{\beta_1}$, assuming other variables included in model do not change.
Some say ceteris paribus, but this seems wrong, because when error term is correlated with our variable $x_j$, then in fact 'other things' are not fixed - the error changes systematically with change of the $x_j$. Fixed remain only variables included in model.

*We can stil predict. The prediction is worse, than in any model including more variables, but it is still better, than no prediction.

*In second model, we suddenly encountered hidden problem - error clusterisation at the level of $j$. It will affect our estimators of $\beta_1$. I will not analyse this in sake of any simplicity.

Introduction of fixed effects:


*

*Now we get it: we observe people many times, both this who give smiles, and those, who get them. We can introduce fixed effects! But what are the consequences?

*At first, lets introduce dummy for $i$. This is not a very drastic case:
Instead of short model:
$y_{i,j} = \beta_0 + \beta_1 x_j + \varepsilon_{i,j}$
we estimate model: 
$y_{i,j} = \beta_0 + \gamma_i + \beta_1 x_j + \varepsilon_{i,j}$
Three things happen:


*

*R$^2$ increases to the exact of model:
$y_{i,j} = \beta_0 + \beta_1 x_j + \beta_3 z_i+ \varepsilon_{i,j}$. We explained variance of variable $z_i$. 

*We can not predict for an unknown person. If we want so make a prediction we must put there $i$. This is the cost of such model and such prediction.

*Estimator of $x_j$ is inconsistent and biased as it was before, because variable $v_j$ is still omitted.


*Let then introduce dummy for $j$. Instead of model:
$y_{i,j} = \beta_0 + \gamma_i + \beta_1 x_j + \varepsilon_{i,j}$
we estimate model:
$y_{i,j} = \beta_0 + \gamma_i + \theta_j + \beta_1 x_j + \varepsilon_{i,j}$
Four things happen:


*

*R$^2$ increases to the level of full model $y_{i,j} = \beta_0 + \beta_1 x_j + \beta_2v_j + \beta_3 z_i + \varepsilon_{i,j}$. We explained variance of full model (great!).

*As before, we can not predict for an unknown person, we must put $i$ and $j$ there.

*We get rid of clustered error problem. Error no longer contains $v_j$ (see 4.).

*Estimator of $\beta_1$ is not estimated due to full collinearity with dummies. AUCH.

Some playground
Below an R code, which allows to simulate, and believe this story by oneself:
n = 200 
d = data.frame(id = 1:(n*n-n))

z = rnorm(n)
v = rnorm(n)
x = v + rnorm(n)

k = 0
for(i in 1:n){
    z_i = z[i]
    for(j in 1:n){
        v_j = v[j]
        x_j = x[j]
        if(i!=j){
            k = k+1
            print(k)
            y_ij = x_j + v_j + z_i + rnorm(1)
            d[k, 'i'] = i
            d[k, 'j'] = j
            d[k, 'z_i'] = z_i
            d[k, 'v_j'] = v_j
            d[k, 'x_j'] = x_j
            d[k, 'y_ij'] = y_ij
        }
    }
}

d$i = as.factor(d$i)
d$j = as.factor(d$j)

m1 = lm(y_ij ~ x_j + v_j + z_i, data = d)
summary(m1)

m2 = lm(y_ij ~ x_j + v_j, data = d)
summary(m2)

m3 = lm(y_ij ~ x_j, data = d)
summary(m3)

m4 = lm(y_ij ~ i  + x_j, data = d)
summary(m4)

m5 = lm(y_ij ~ z_i  + x_j, data = d)
summary(m5)

m6 = lm(y_ij ~ j  + i + x_j, data = d)
summary(m6)

