I used to think that "random effects model" in econometrics corresponds to a "mixed model with random intercept" outside of econometrics, but now I am not sure. Does it?

Econometrics uses terms like "fixed effects" and "random effects" somewhat differently from the literature on mixed models, and this causes a notorious confusion. Let us consider a simple situation where $y$ linearly depends on $x$ but with a different intercept in different groups of measurements:

$$y_{it} = \beta x_{it} + u_i + \epsilon_{it}.$$

Here each unit/group $i$ is observed at different timepoints $t$. Econometricians call it "panel data".

  • In mixed models terminology, we can treat $u_i$ as a fixed effect or as a random effect (in this case, it's random intercept). Treating it as fixed means fitting $\hat \beta$ and $\hat u_i$ to minimize squared error (i.e. running OLS regression with dummy group variables). Treating it as random means that we additionally assume that $u_i\sim\mathcal N(u_0,\sigma^2_u)$ and use maximum likelihood to fit $u_0$ and $\sigma^2_u$ instead of fitting each $u_i$ on its own. This leads to the "partial pooling" effect, where the estimates $\hat u_i$ get shrunk toward their mean $\hat u_0$.

    R formula when treating group as fixed:    y ~ x + group
    R formula when treating group as random:   y ~ x + (1|group)
    
  • In econometrics terminology, we can treat this whole model as a fixed effects model or as a random effects model. The first option is equivalent to the fixed effect above (but econometrics has its own way of estimating $\beta$ in this case, called "within" estimator). I used to think that the second option is equivalent to the random effect above; e.g. @JiebiaoWang in his highly upvoted answer to What is a difference between random effects-, fixed effects- and marginal model? says that

    In econometrics, the random-effects model may only refer to random intercept model as in biostatistics

Okay --- let us test if this understanding is correct. Here is some random data generated by @ChristophHanck in his answer to What is the difference between fixed effect, random effect and mixed effect models? (I put the data here on pastebin for those who do not use R):

enter image description here

@Christoph does two fits using econometrics approaches:

fe <- plm(stackY~stackX, data = paneldata, model = "within")
re <- plm(stackY~stackX, data = paneldata, model = "random")

The first one yields the estimate of beta equal to -1.0451, the second one 0.77031 (yes, positive!). I tried to reproduce it with lm and lmer:

l1 = lm(stackY ~ stackX + as.factor(unit), data = paneldata)
l2 = lmer(stackY ~ stackX + (1|as.factor(unit)), data = paneldata)

The first one yields -1.045 in perfect agreement with the within estimator above. Cool. But the second yields -1.026, which is miles away from the random effects estimator. Heh? What is going on? In fact, what is plm even doing, when called with model = "random"?

Whatever it is doing, can one somehow understand it via the mixed models perspective?

And what is the intuition behind whatever it is doing? I read in a couple of econometrics places that random effects estimator is a weighted average between the fixed effects estimator and the "between" estimator which is more or less regression slope if we do not include group identity in the model at all (this estimate is strongly positive in this case, around 4.) E.g. @Andy writes here:

The random effects estimator then uses a matrix weighted average of the within and between variation of your data. [...] This makes random effects more efficient[.]

Why? Why would we want this weighted average? And in particular, why would we want it instead of running a mixed model?

up vote 13 down vote accepted
+100

Summary: the "random-effects model" in econometrics and a "random intercept mixed model" are indeed the same models, but they are estimated in different ways. The econometrics way is to use FGLS, and the mixed model way is to use ML. There are different algorithms of doing FGLS, and some of them (on this dataset) produce results that are very close to ML.


1. Differences between estimation methods in plm

I will answer with my testing on plm(..., model = "random") and lmer(), using the data generated by @ChristophHanck.

According to the plm package manual, there are four options for random.method: the method of estimation for the variance components in the random effects model. @amoeba used the default one swar (Swamy and Arora, 1972).

For random effects models, four estimators of the transformation parameter are available by setting random.method to one of "swar" (Swamy and Arora (1972)) (default), "amemiya" (Amemiya (1971)), "walhus" (Wallace and Hussain (1969)), or "nerlove" (Nerlove (1971)).

I tested all the four options using the same data, getting an error for amemiya, and three totally different coefficient estimates for the variable stackX. The ones from using random.method='nerlove' and 'amemiya' are nearly equivalent to that from lmer(), -1.029 and -1.025 vs -1.026. They are also not very different from that obtained in the "fixed-effects" model, -1.045.

# "amemiya" only works using the most recent version:
# install.packages("plm", repos="http://R-Forge.R-project.org")

re0 <- plm(stackY~stackX, data = paneldata, model = "random") #random.method='swar'
re1 <- plm(stackY~stackX, data = paneldata, model = "random",  random.method='amemiya')
re2 <- plm(stackY~stackX, data = paneldata, model = "random",  random.method='walhus')
re3 <- plm(stackY~stackX, data = paneldata, model = "random",  random.method='nerlove')
l2  <- lmer(stackY~stackX+(1|as.factor(unit)), data = paneldata)

coef(re0)     #    (Intercept)   stackX    18.3458553   0.7703073 
coef(re1)     #    (Intercept)   stackX    30.217721   -1.025186 
coef(re2)     #    (Intercept)   stackX    -1.15584     3.71973 
coef(re3)     #    (Intercept)   stackX    30.243678   -1.029111 
fixef(l2)     #    (Intercept)   stackX    30.226295   -1.026482 

Unfortunately I do not have time right now, but interested readers can find the four references, to check their estimation procedures. It would be very helpful to figure out why they make such a difference. I expect that for some cases, the plm estimation procedure using the lm() on transformed data should be equivalent to the maximum likelihood procedure utilized in lmer().

2. Comparison between GLS and ML

The authors of plm package did compare the two in Section 7 of their paper: Yves Croissant and Giovanni Millo, 2008, Panel Data Econometrics in R: The plm package.

Econometrics deal mostly with non-experimental data. Great emphasis is put on specification procedures and misspecification testing. Model specifications tend therefore to be very simple, while great attention is put on the issues of endogeneity of the regressors, dependence structures in the errors and robustness of the estimators under deviations from normality. The preferred approach is often semi- or non-parametric, and heteroskedasticity-consistent techniques are becoming standard practice both in estimation and testing.

For all these reasons, [...] panel model estimation in econometrics is mostly accomplished in the generalized least squares framework based on Aitken’s Theorem [...]. On the contrary, longitudinal data models in nlme and lme4 are estimated by (restricted or unrestricted) maximum likelihood. [...]

The econometric GLS approach has closed-form analytical solutions computable by standard linear algebra and, although the latter can sometimes get computationally heavy on the machine, the expressions for the estimators are usually rather simple. ML estimation of longitudinal models, on the contrary, is based on numerical optimization of nonlinear functions without closed-form solutions and is thus dependent on approximations and convergence criteria.


3. Update on mixed models

I appreciate that @ChristophHanck provided a thorough introduction about the four random.method used in plm and explained why their estimates are so different. As requested by @amoeba, I will add some thoughts on the mixed models (likelihood-based) and its connection with GLS.

The likelihood-based method usually assumes a distribution for both the random effect and the error term. A normal distribution assumption is commonly used, but there are also some studies assuming a non-normal distribution. I will follow @ChristophHanck's notations for a random intercept model, and allow unbalanced data, i.e., let $T=n_i$.

The model is \begin{equation} y_{it}= \boldsymbol x_{it}^{'}\boldsymbol\beta + \eta_i + \epsilon_{it}\qquad i=1,\ldots,m,\quad t=1,\ldots,n_i \end{equation} with $\eta_i \sim N(0,\sigma^2_\eta), \epsilon_{it} \sim N(0,\sigma^2_\epsilon)$.

For each $i$, $$\boldsymbol y_i \sim N(\boldsymbol X_{i}\boldsymbol\beta, \boldsymbol\Sigma_i), \qquad\boldsymbol\Sigma_i = \sigma^2_\eta \boldsymbol 1_{n_i} \boldsymbol 1_{n_i}^{'} + \sigma^2_\epsilon \boldsymbol I_{n_i}.$$ So the log-likelihood function is $$const -\frac{1}{2} \sum_i\mathrm{log}|\boldsymbol\Sigma_i| - \frac{1}{2} \sum_i(\boldsymbol y_i - \boldsymbol X_{i}\boldsymbol\beta)^{'}\boldsymbol\Sigma_i^{-1}(\boldsymbol y_i - \boldsymbol X_{i}\boldsymbol\beta).$$

When all the variances are known, as shown in Laird and Ware (1982), the MLE is $$\hat{\boldsymbol\beta} = \left(\sum_i\boldsymbol X_i^{'} \boldsymbol\Sigma_i^{-1} \boldsymbol X_i \right)^{-1} \left(\sum_i \boldsymbol X_i^{'} \boldsymbol\Sigma_i^{-1} \boldsymbol y_i \right),$$ which is equivalent to the GLS $\hat\beta_{RE}$ derived by @ChristophHanck. So the key difference is in the estimation for the variances. Given that there is no closed-form solution, there are several approaches:

  • directly maximization of the log-likelihood function using optimization algorithms;
  • Expectation-Maximization (EM) algorithm: closed-form solutions exist, but the estimator for $\boldsymbol \beta$ involves empirical Bayesian estimates of the random intercept;
  • a combination of the above two, Expectation/Conditional Maximization Either (ECME) algorithm (Schafer, 1998; R package lmm). With a different parameterization, closed-form solutions for $\boldsymbol \beta$ (as above) and $\sigma^2_\epsilon$ exist. The solution for $\sigma^2_\epsilon$ can be written as $$\sigma^2_\epsilon = \frac{1}{\sum_i n_i}\sum_i(\boldsymbol y_i - \boldsymbol X_{i} \hat{\boldsymbol\beta})^{'}(\hat\xi \boldsymbol 1_{n_i} \boldsymbol 1_{n_i}^{'} + \boldsymbol I_{n_i})^{-1}(\boldsymbol y_i - \boldsymbol X_{i} \hat{\boldsymbol\beta}),$$ where $\xi$ is defined as $\sigma^2_\eta/\sigma^2_\epsilon$ and can be estimated in an EM framework.

In summary, MLE has distribution assumptions, and it is estimated in an iterative algorithm. The key difference between MLE and GLS is in the estimation for the variances.

Croissant and Millo (2008) pointed out that

While under normality, homoskedasticity and no serial correlation of the errors OLS are also the maximum likelihood estimator, in all the other cases there are important differences.

In my opinion, for the distribution assumption, just as the difference between parametric and non-parametric approaches, MLE would be more efficient when the assumption holds, while GLS would be more robust.

  • I would suspect that the issue with the error message is somehow related to me generating the variables as vectors? Maybe plm prefers the data to be stored differently? – Christoph Hanck Oct 4 '16 at 7:56
  • 1
    nerlove works well here, but is not implemented for unbalanced panels, as I found out by deleting 1 observation from the last panel and trying to run all methods. – amoeba Oct 4 '16 at 9:55
  • 2
    @ChristophHanck @amoeba The plm error for random.method="amemiya" occurs to me that they probably should use X[, -1, drop=FALSE] instead of X[, -1] to retain the matrix format of X[, -1] when there is only one covariate in the model. Anyway, I tried to overcome this by adding a standard normal variable to the formula. amemiya reproduces the result with an estimate of -1.02 and it works for unbalanced data as well. – Randel Oct 7 '16 at 22:38
  • 3
    @jiebiao-wang @ChristophHanck @amoeba the current development version of plm runs fine with random.method="amemiya": var std.dev share idiosyncratic 0.6360 0.7975 0.002 individual 313.6510 17.7102 0.998 theta: 0.9841 – Helix123 Mar 5 '17 at 22:24
  • 1
    Hello @JiebiaoWang. I figured that after your Update, your answer does answer my question satisfactorily. I took the liberty to make some edits and to insert an update about amemiya and a quote on the ML vs GLS. I am marking it as accepted and am going to award it a bounty. Cheers. – amoeba Jun 9 '17 at 9:58

This answer doesn't comment on mixed models, but I can explain what the random-effects estimator does and why it screws up on that graph.

Summary: the random-effects estimator assumes $E[u_i \mid x ] = 0$, which is not true in this example.


What is the random effects estimator doing?

Assume we have the model:

$$ y_{it} = \beta x_{it} + u_i + \epsilon_{it}$$

We have two dimensions of variation: groups $i$ and time $t$. To estimate $\beta$ we could:

  1. Only use time-series variation within a group. This is what the fixed-effect estimator does (and this is why it's also often called the within estimator.)
  2. If $u_i$ is random, we could use only cross-sectional variation between the time-series means of groups. This is known as the between estimator.

    Specifically, for each group $i$, take the average over time of the above panel data model to get:

    $$ \bar{y}_{i} = \beta \bar{x}_{i} + v_i \quad \quad \text{ where } v_i = u_i + \bar{\epsilon}_i$$

    If we run this regression, we get the between estimator. Observe that it is a consistent estimator if the effects $u_i$ are random white noise, uncorrelated with $x$! If this is the case, then completely tossing the between group variation (as we do with the fixed effects estimator) is inefficient.

The random-effects estimator of econometrics combines the (1) within estimator (i.e. fixed effects estimator) and (2) the between estimator in a way to maximize efficiency. It is an application of generalized least squares and the basic idea is inverse variance weighting. To maximize efficiency, the random-effects estimator calculates $\hat \beta$ as a weighted average of the within estimator and the between estimator.

What's going on in that graph...

The effects $u_i$ are not random! $E[u_i \mid x ] \neq 0$. The group effects $u_i$ are not orthogonal to $x$ (in a statistical sense), rather, the group effects have a clear positive relationship with $x$.

The between estimator assumes $E[u_i \mid x ] = 0$. The between estimator says, "sure I can impose $E[u_i \mid x ] = 0$, by making $\hat \beta$ positive!"

Then in turn, the random-effects estimator is off because it's a weighted average of the within estimator and the between estimator.

  • +1, thanks Matthew. Not sure why somebody downvoted your answer. I am looking for an answer establishing the connection to mixed models so I won't accept yours, but I still find it helpful for this discussion. If you can expand a little bit on how GLS and inverse-variance weighting are applied and computed here, it would be very useful. – amoeba Oct 3 '16 at 18:47

In this answer, I would like to elaborate a little on Matthew's +1 answer regarding the GLS perspective on what the econometrics literature calls the random effects estimator.

GLS perspective

Consider the linear model \begin{equation} y_{it}=\alpha + X_{it}\beta+u_{it}\qquad i=1,\ldots,m,\quad t=1,\ldots,T \end{equation} If it held that $E(u_{it}\vert X_{it})=0$ we could simply estimate the model by pooled OLS, which amounts to ignoring the panel data structure and simply lump all $n=mT$ observations together.

We model the $u_{it}$ using the error-component model \begin{equation} u_{it}=\eta_i+\epsilon_{it} \end{equation}

In matrix notation, the model can be written as \begin{equation} y=\alpha \iota_{mT}+X\beta+D\eta+\epsilon \end{equation} where $y$ and $\epsilon$ are $n$-vectors with typical elements $y_{it}$ and $\epsilon_{it}$, and $D$ is an $n \times m$ (one column per unit) matrix of dummy variables. $D$ is such that if a row corresponds to an observation belonging to unit $i$, then $D$ has a one in column $i$ and 0 else, $i=1,\ldots,m$.

We furthermore assume $$ E(\epsilon\epsilon^\prime)=\sigma_\epsilon^2I $$

The individual-specific effects $\eta$ must be independent of the $\epsilon_{it}$. The random-effects estimator, unlike the fixed effects (again, econometrics terminology) one, however additionally requires the stronger assumption that \begin{equation} E(\eta_i\vert X)=0 \end{equation} Under this assumption, pooled OLS would be unbiased, but we can derive a GLS estimator. Assume that the $\eta_i$ are IID with mean zero and variance $\sigma^2_\eta$.

This assumption accounts for the term random effects. Assuming, moreover, that the two error components are independent, it is easy to see that \begin{align*} \operatorname{Var}(u_{it})&=\sigma^2_\eta+\sigma^2_\epsilon\\ \operatorname{Cov}(u_{it},u_{is})&=\sigma^2_\eta\\ \operatorname{Cov}(u_{it},u_{js})&=0\qquad\text{for all } i\neq j \end{align*}

We then get the following $n\times n$ variance-covariance matrix $\Omega$: $$ \Omega= \begin{pmatrix} \Sigma&O&\cdots&O\\ O&\Sigma&\cdots&O\\ \vdots&\vdots&&\vdots\\ O&O&\cdots&\Sigma \end{pmatrix} $$ Here, $$ \Sigma=\sigma^2_\eta \iota\iota^\prime+\sigma^2_\epsilon I_T $$ with $\iota$ a $T$-vector of ones. We may hence write $$ \Omega=\sigma^2_\eta (I_m\otimes\iota\iota^\prime)+\sigma^2_\epsilon (I_m\otimes I_T) $$ For the GLS estimator $$ \hat\beta_{RE}=(X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y $$ we require $\Omega^{-1}$. To this end, let $J_T=\iota\iota^\prime$, $\bar J_T=J_T/T$ and $E_T=I_T-\bar J_T$. Then, write $$ \Omega=T\sigma^2_\eta (I_m\otimes\bar J_T)+\sigma^2_\epsilon (I_m\otimes E_T)+\sigma^2_\epsilon (I_m\otimes \bar J_T) $$ or, collecting terms with the same matrices, $$ \Omega=(T\sigma^2_\eta+\sigma^2_\epsilon) (I_m\otimes\bar J_T)+\sigma^2_\epsilon (I_m\otimes E_T) $$ Idempotency of $P=I_m\otimes\bar J_T$ and $Q=I_m\otimes E_T$ then allows us to show that $$\Omega^{-1}=\frac{1}{\sigma^2_1}P+\frac{1}{\sigma^2_\epsilon}Q= -\frac{\sigma^2_\eta}{\sigma^2_1\sigma^2_\epsilon}(I_m\otimes\iota\iota^\prime) + \frac{1}{\sigma^2_\epsilon}(I_m\otimes I_T),$$ where $\sigma^2_1=T\sigma^2_\eta+\sigma^2_\epsilon$.

Gauss-Markov logic then explains why the random effects estimator may be useful, as it is a more efficient estimator than pooled OLS or fixed effects under the given assumptions (provided, which is a very big if in many panel data applications, that the $\eta_i$ are indeed uncorrelated with the regressors). In short, GLS is more efficient because the error covariance matrix is not homoskedastic in this model.

One can show that GLS estimate can be obtained by running OLS on the partially demeaned data: $$(y_{it}-\theta \bar y_{i\cdot}) = (X_{it} - \theta \bar X_{i\cdot})\beta + (u_{it} - \theta u_{i\cdot}),$$ where $\theta = 1-\sigma_\eta/\sigma_1$. For $\theta=1$ one gets the fixed effect ("within") estimator. For $\theta\to -\infty$ one gets the "between" estimator. The GLS estimator is a weighted average between the two. (For $\theta=0$ one gets the pooled OLS estimator.)

Feasible GLS

To make an FGLS approach practical, we require estimators of $\sigma^2_1$ and $\sigma^2_\epsilon$. Baltagi, Econometric Analysis of Panel Data, p. 16 (quoting from the 3rd edition), discusses the following options on how to proceed.

Assume first we observe $u_{it}$. Then,

$$\hat\sigma^2_1=T\frac{1}{m}\sum_{i=1}^m\bar{u}_{i\cdot}^2$$ and $$\hat\sigma^2_\epsilon=\frac{1}{m(T-1)}\sum_{i=1}^m\sum_{t=1}^T\left(u_{it}-\frac{1}{m}\sum_{i=1}^m\bar{u}_{i\cdot}\right)^2$$ would be good estimators of their parameters, with $\bar{u}_{i\cdot}$ the time-average corresponding to the obseravations of unit $i$.

The Wallace and Hussein (1969) approach consists of replacing $u$ with residuals of a pooled OLS regression (which, after all, still is unbiased and consistent under the present assumptions).

The Amemiya (1971) approach suggests using FE (or LSDV) residuals instead. As a computational matter, we impose the restriction that $\sum_i\eta_i=0$ to circumvent the dummy variable trap so as to be able to get $\hat\alpha=\bar y_{\cdot\cdot}-\bar X_{\cdot\cdot}'\hat\beta_{FE}$ with $\cdot\cdot$ denoting grand averages over $i$ and $t$ for the LSDV residuals $\hat u=y-\hat\alpha-X\hat\beta_{FE}$.

The default Swamy and Arora (1972) approach estimates $$ \hat\sigma^2_\epsilon=[y'Q(I-X(X'QX)^{-1}X'Q)y]/[m(T-1)-K] $$ and $$ \hat\sigma^2_1=[y'P(I-Z(Z'PX)^{-1}Z'P)y]/[m-K-1] $$ Here, $Z=(\iota_{mT}\quad X)$.

The Nerlove (1971) approach estimates $\sigma_\eta^2$ from $\sum_{i=1}^m(\hat\eta_i-\bar{\hat\eta})^2/(m-1)$ where the $\hat\eta_i$ are dummies from a fixed effects regression and $\hat\sigma^2_\epsilon$ is estimated from the within residual sums of squares from this regression, with $mT$ in the denominator.

I am also very surprised that these make such a big difference as shown by Randel's calculations!

EDIT:

Regarding the differences, the estimates of the error components may be retrived in the plm package, and indeed return vastly different results, explaining the difference in the point estimates for $\beta$ (as per @Randel's answer, amemiya throws an error that I did not attempt to fix):

> ercomp(stackY~stackX, data = paneldata, method = "walhus")
                  var std.dev share
idiosyncratic 21.0726  4.5905 0.981
individual     0.4071  0.6380 0.019
theta:  0.06933  
> ercomp(stackY~stackX, data = paneldata, method = "swar")
                 var std.dev share
idiosyncratic 0.6437  0.8023 0.229
individual    2.1732  1.4742 0.771
theta:  0.811  
> ercomp(stackY~stackX, data = paneldata, method = "nerlove")
                   var  std.dev share
idiosyncratic   0.5565   0.7460 0.002
individual    342.2514  18.5000 0.998
theta:  0.9857  

I suspect that the estimators of the error components are also not consistent in my example in the sister thread where I aim to demonstrate differences between FE and RE using data where the individual effects and $X$ are correlated. (In fact, they cannot be, because they ultimately drive away the RE estimate from the FE estimate as per the fact that RE is a weighted average of FE and between estimation with weights determined by the error component estimates. So, if RE is not consistent, that must ultimately be due to these estimates.)

If you replace the "offending" feature of that example,

alpha = runif(n,seq(0,step*n,by=step),seq(step,step*n+step,by=step))

by simply, say,

alpha = runif(n)

so random effects that are uncorrelated with $X$, you get RE point estimates for $\beta$ very close to the true value $\beta=-1$ for all variantes of estimating the error components.


References

Amemiya, T., 1971, The estimation of the variances in a variance-components model, International Economic Review 12, 1–13.

Baltagi, B. H., Econometric Analysis of Panel Data, Wiley.

Nerlove, M., 1971a, Further evidence on the estimation of dynamic economic relations from a time-series of cross-sections, Econometrica 39, 359–382.

Swamy, P.A.V.B. and S.S. Arora, 1972, The exact finite sample properties of the estimators of coefficients in the error components regression models, Econometrica 40, 261–275.

Wallace, T.D. and A. Hussain, 1969, The use of error components models in combining cross-section and time-series data, Econometrica 37, 55–72.

  • 4
    +1. Thanks Christoph, this is helpful and I am happy to finally see some mathematical details in this thread. It would be great to look up how the four methods implemented in plm and listed by Randel work and update your answer with some comments about it. If not the detailed descriptions, then at least some brief notes about what is going on. Do you think you would be able to look it up? I am happy to offer a bounty for that :-) My naive approach would be to estimate both sigmas from the fixed effects solution. Does it correspond to one of the "named" methods? – amoeba Oct 4 '16 at 8:57
  • @amoeba, I included some comments on how to estimate the variances in the error component model. Your suggestion then appears to be closely related to that of Amemiya. – Christoph Hanck Oct 4 '16 at 12:23
  • Very nice, thanks. Isn't Nerlove also using regression with dummies? In fact, I don't quite understand what is the difference between Amemiya and Nerlove. My "naive" suggestion was to fit dummy regression, use the residual variance as an estimate of $\sigma_\epsilon$ and use the variance of dummy coefficients as an estimate of $\sigma_\eta$. It looks like that's what Nerlove is doing. I am not sure I understand what Amemiya is doing and how it's different. (And I agree that the huge over-aching question remains as to why these methods yield such a difference in this case.) – amoeba Oct 4 '16 at 13:39
  • Yes, both use regression with dummies. As far as I understand, one difference between Amemiya and Nerlove is the denominator for the degrees of freedom correction. Another is that I am not sure that the variance of the estimated dummy coefficients is the same as the variance of the residuals. Another crucial one is that Nerlove directly aims to estimate $\sigma_\eta^2$, whereas you would have to back out the estimate via $(\hat\sigma_1^2-\hat\sigma_\epsilon^2)/T$ for the three others, and one known disadvantage of these is that there is no guarantee that these are nonnegative. – Christoph Hanck Oct 4 '16 at 13:58
  • 1
    Thanks. I made an edit to supply a more explicit formula for $\Omega^{-1}$, you might want to double-check (but I think it's correct). I started a bounty that I am going to award to your answer. However, I am still looking for an answer that will draw a connection to mixed models, contrast GLS to MLE, and explain why and when one should prefer which approach (none of the current answers does it, so there is currently no answer I would tick as "accepted"). It's interesting that MLE (as implemented by lmer) yields variance estimates that are very close to Nerlove ones. – amoeba Oct 5 '16 at 21:23

I am not really familiar enough with R to comment on your code, but the simple random intercept mixed model should be identical to the RE MLE estimator, and very close to the RE GLS estimator, except when total $N = \sum_i T_i$ is small and the data are unbalanced. Hopefully, this will be useful in diagnosing the problem. Of course, this is all assuming that the RE estimator is appropriate.

Here is some Stata showing the equivalence (requires esttab and eststo from SSC):

set more off
estimates clear
webuse nlswork, clear
eststo, title(mixed): mixed ln_w grade age c.age#c.age ttl_exp tenure c.tenure#c.tenure || id: // Mixed estimator
eststo, title(MLE): xtreg ln_w grade age c.age#c.age ttl_exp tenure c.tenure#c.tenure, i(id) mle // MLE RE estimator 
eststo, title(GLS): xtreg ln_w grade age c.age#c.age ttl_exp tenure c.tenure#c.tenure, i(id) re // GLS RE estimato
esttab *, b(a5) se(a5) mtitle 

Here's the output of the last line:

. esttab *, b(a5) se(a5) mtitle 

------------------------------------------------------------
                      (1)             (2)             (3)   
                    mixed             MLE             GLS   
------------------------------------------------------------
main                                                        
grade            0.070790***     0.070790***     0.070760***
              (0.0017957)     (0.0017957)     (0.0018336)   

age              0.031844***     0.031844***     0.031906***
              (0.0027201)     (0.0027202)     (0.0027146)   

c.age#c.age   -0.00065130***  -0.00065130***  -0.00065295***
             (0.000044965)    (0.000044971)    (0.000044880)   

ttl_exp          0.035228***     0.035228***     0.035334***
              (0.0011382)     (0.0011392)     (0.0011446)   

tenure           0.037134***     0.037134***     0.037019***
              (0.0015715)     (0.0015723)     (0.0015681)   

c.tenure#c~e   -0.0018382***   -0.0018382***   -0.0018387***
             (0.00010128)    (0.00010128)    (0.00010108)   

_cons             0.14721***      0.14721***      0.14691** 
               (0.044725)      (0.044725)      (0.044928)   
------------------------------------------------------------
lns1_1_1                                                    
_cons            -1.31847***                                
               (0.013546)                                   
------------------------------------------------------------
lnsig_e                                                     
_cons            -1.23024***                                
              (0.0046256)                                   
------------------------------------------------------------
sigma_u                                                     
_cons                             0.26754***                
                              (0.0036240)                   
------------------------------------------------------------
sigma_e                                                     
_cons                             0.29222***                
                              (0.0013517)                   
------------------------------------------------------------
N                   28099           28099           28099   
------------------------------------------------------------
Standard errors in parentheses
* p<0.05, ** p<0.01, *** p<0.001

In your data, the assumptions for using the RE estimator are not satisfied since the group effect is clearly correlated with x, so you get very different estimates. The GLS RE estimator actually a generalized method of moments (GMM) estimator that is a matrix-weighted average of the between and within estimators. The within estimator is going to be OK here, but the between is going to be profoundly screwed, showing large positive effects of X. So GLS will be mostly the between estimator. The MLE RE is an MLE that maximizes the likelihood of the random-effects model. They are no longer expected to produce the same answer. Here the mixed estimator is giving something very close to FE "Within" estimator:

. esttab *, b(a5) se(a5) mtitle 

----------------------------------------------------------------------------
                      (1)             (2)             (3)             (4)   
                    mixed             GLS             MLE          Within   
----------------------------------------------------------------------------
main                                                                        
x                -1.02502***      0.77031**       3.37983***     -1.04507***
               (0.092425)       (0.26346)       (0.20635)      (0.093136)   

_cons             30.2166***      18.3459***      0.49507         30.3492***
                (5.12978)       (2.31566)             (.)       (0.62124)   
----------------------------------------------------------------------------
lns1_1_1                                                                    
_cons             2.87024***                                                
                (0.20498)                                                   
----------------------------------------------------------------------------
lnsig_e                                                                     
_cons            -0.22598**                                                 
               (0.077195)                                                   
----------------------------------------------------------------------------
sigma_u                                                                     
_cons                                             2.40363                   
                                                (1.28929)                   
----------------------------------------------------------------------------
sigma_e                                                                     
_cons                                             4.23472***                
                                                (0.37819)                   
----------------------------------------------------------------------------
N                      96              96              96              96   
----------------------------------------------------------------------------
Standard errors in parentheses
* p<0.05, ** p<0.01, *** p<0.001

Here is the Stata code for the above table:

clear
set more off
estimates clear

input int(obs id t) double(y x)
1      1           1  2.669271  0.5866982
2      1           2  1.475540  1.3500454
3      1           3  4.430008  0.6830919
4      1           4  2.162789  0.5845966
5      1           5  2.678108  1.0038879
6      1           6  3.456636  0.5863289
7      1           7  1.769204  2.3375403
8      1           8  3.413790  0.9640034
9      2           1  4.017493  1.5084121
10     2           2  4.218733  2.8982499
11     2           3  4.509530  3.2141335
12     2           4  6.106228  2.0317799
13     2           5  5.161379  2.1231733
14     2           6  2.724643  4.3369017
15     2           7  4.500306  1.9141065
16     2           8  4.119322  2.8667938
17     3           1  9.987779  2.3961969
18     3           2  7.768579  3.5509275
19     3           3  9.379788  3.3284869
20     3           4 10.035937  2.2997389
21     3           5 11.752360  2.8143474
22     3           6  9.500264  2.1825704
23     3           7  8.921687  5.0126462
24     3           8  8.269932  3.4046339
25     4           1 12.101253  3.2928033
26     4           2 11.482337  3.1645218
27     4           3 10.648010  4.8073987
28     4           4  9.687320  5.3394193
29     4           5 12.796925  3.1197431
30     4           6  9.971434  4.6512983
31     4           7 10.239717  4.7709378
32     4           8 12.245207  2.7952426
33     5           1 18.473320  5.8421967
34     5           2 19.097212  4.9425391
35     5           3 19.460495  4.9166172
36     5           4 18.642305  4.9856035
37     5           5 17.723912  5.0594425
38     5           6 16.783248  4.8615618
39     5           7 16.100984  6.2069167
40     5           8 18.851351  3.8856152
41     6           1 19.683171  7.5568816
42     6           2 21.104231  6.7441900
43     6           3 22.115529  6.4486514
44     6           4 22.061362  5.3727434
45     6           5 22.457905  5.8665798
46     6           6 21.424413  6.0578997
47     6           7 23.475946  4.4024323
48     6           8 24.884950  4.1596914
49     7           1 25.809011  7.6756255
50     7           2 25.432828  7.7910756
51     7           3 26.790387  7.3858301
52     7           4 24.640850  8.2090606
53     7           5 26.050086  7.3779219
54     7           6 25.297148  6.8098617
55     7           7 26.551229  7.6694272
56     7           8 26.669760  6.4425772
57     8           1 26.409669  8.3040894
58     8           2 26.570003  8.4686087
59     8           3 29.018818  7.2476785
60     8           4 30.342613  4.5207729
61     8           5 26.819959  8.7935557
62     8           6 27.147711  8.3141224
63     8           7 26.168568  9.0148308
64     8           8 27.653552  8.2081808
65     9           1 34.120485  7.8415520
66     9           2 31.286463  9.7234259
67     9           3 35.763403  6.9202442
68     9           4 31.974599  9.0078286
69     9           5 32.273719  9.4954288
70     9           6 29.666208 10.2525763
71     9           7 30.949857  9.4751679
72     9           8 33.485967  8.1824810
73    10           1 36.183128 10.7891587
74    10           2 37.706116  9.7119548
75    10           3 38.582725  8.6388290
76    10           4 35.876781 10.8259279
77    10           5 37.111179  9.9805046
78    10           6 40.313149  7.7487456
79    10           7 38.606329 10.2891107
80    10           8 37.041938 10.3568765
81    11           1 42.617586 12.1619185
82    11           2 41.787495 11.1420338
83    11           3 43.944968 11.1898730
84    11           4 43.446467 10.8099599
85    11           5 43.420819 11.2696770
86    11           6 42.367318 11.6183869
87    11           7 43.543785 11.1336555
88    11           8 43.750271 12.0311065
89    12           1 46.122429 12.3528733
90    12           2 47.604306 11.4522787
91    12           3 45.568748 13.6906476
92    12           4 48.331177 12.3561907
93    12           5 47.143246 11.7339915
94    12           6 44.461190 13.3898768
95    12           7 46.879044 11.4054972
96    12           8 46.314055 12.3143487
end

eststo, title(mixed): mixed y x || id:, mle // Mixed estimator
eststo, title(GLS): xtreg y x, i(id) re     // GLS RE estimato
eststo, title(MLE): xtreg y x, i(id) mle    // MLE RE estimator 
eststo, title(Within): xtreg y x, i(id) fe  // FE Within estimator 
eststo, title(Between): xtreg y x, i(id) be // Between estimator 

esttab *, b(a5) se(a5) mtitle 
  • +1. Thanks, Dimitriy, it is definitely helpful to see Stata's output on the same toy dataset. I have a question about MLE estimator. I thought that the mixed model approach (mixed in Stata and lmer in R) is also maximum likelihood or sometimes "restricted maximum likelihood" (I can use both in my lmer call by setting REML=T or REML=F and they give almost identical results). However mixed model approach gives a very sensible and correct result, whereas what Stat calls "MLE" gives a nonsense result in this case. What is the difference? What exactly does Stat's "MLE" refer to? – amoeba Oct 4 '16 at 9:04
  • 2
    @amoeba Both mixed, mle and xtreg, mle are MLE estimators, but the likelihood functions are somewhat different. See here for the former, and here for the latter. I don't quite understand why the mixed model is so robust. – Dimitriy V. Masterov Oct 4 '16 at 14:17
  • xtmixed is what mixed was called in older versions of Stata. For your data, that equivalence clearly does not hold, while it holds for my data, as the manual suggests. – Dimitriy V. Masterov Oct 5 '16 at 0:04
  • ssc install estout Although my recollection is that it had different functionality in different versions, and failed to be backward compatible. – StasK Oct 7 '16 at 14:41
  • 1
    @StasK put me in contact with Stata tech support and they said it's likely a bug in xtreg, mle. "In general the result should be the same [...]. These kind of differences normally arise when there are identification problems in the estimation of the parameters of the model. [...] I actually checked the condition number for the variance-covariance matrix resulting from both computations and that number is basically infinite for -xtreg, mle- and over 4000 for -mixed ,mle-. [...] the developers [...] are going to evaluate the problem to determine whether a fixed code is needed." – amoeba Oct 11 '16 at 20:54

Let me confuse things even more:

ECONOMETRICS - FIXED EFFECTS APPROACH
The "fixed effects" approach in econometrics for panel data, is a way to estimate the slope coefficients (the betas), by "by-passing" the existence of the individual effects variable $\alpha_i$, and so by not making any assumption as to whether it is "fixed" or "random". This is what the "First Difference" estimator (using first differences of the data) and the "Within" estimator (using deviations from time-averages) do: they manage to estimate only the betas.

For a more traditional approach that does explicitly treat the individual effects (the "intercepts") as constants, we use the Least Squares Dummy Variable (LSDV) Estimator, which provides also estimates for the $\alpha_i$'s Note: in the linear model the three estimators algebraically coincide as regards the produced estimates for the betas - but only in the linear model.

Discussion (partly excerpted from class notes)

"The main advantage of the fixed effects approach is that we do not need to make any assumptions about the nature of the individual effects. We should apply it whenever we suspect that the latter are correlated with one or more of the regressors since in this case ignoring the presence of such correlation and naively applying OLS on the pooled model produces inconsistent estimators. Despite its appeal on grounds of the minimal assumptions that we need to make concerning the individual effects, the fixed effects approach has certain limitations. First, coefficients of time invariant regressors cannot be estimated since these variables are differenced out along with the unobservable individual effects. Second, the individual effects (in case we use the LSDV estimator) cannot be consistently estimated (except if we let the time dimension go to infinity)."

ECONOMETRICS - RANDOM EFFECTS APPROACH
In the "traditional" econometric Random Effects approach we assume that the individual "intercepts" $\alpha_i$ are "permanent random components" while the "usual" error terms are "transitory" error components.

In an interesting extension, the additional randomness arises from the existence of a random time effect, common to all cross sections but time varying, alongside a fixed(constant) individual effect and the error term. This "time effect" for example may represent an aggregate shock at economy-wide level that affects equally all households. Such aggregate disturbances are indeed observed and so it appears to be a realistic modelling choice.

Here the "Random Effects" Estimator is a Generalized Least Squares (GLS) estimator, for increased efficiency.

Now, one more conceived estimator, the "Between" Estimator, performs OLS on the time-averaged observations. As a matter of algebra it has been shown that the GLS estimator can be obtained as a weighted average of the Within and the Between estimators, where the weights are not arbitrary but relate to the VCV matrices of the two.

...and there is also the variants of "Uncorrelated Random Effects" and "Correlated Random Effects" models.

I hope the above help make the contrast with the "mixed effects" models.

  • +1, thanks Alecos. This is helpful, but the relationship of all of that to the mixed models approach remains unclear to me. I start to suspect that perhaps there is no relationship whatsoever. The between and within estimators (and that the within is equivalent to the class dummies) are clear by the way; my confusion is only about the random effects approach. – amoeba Oct 3 '16 at 18:59

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.