Can correlated random effects "steal" the variability (and the significance) from the regression coefficient?

Question

I have time-series count data $N_{i,j}$ (population sizes in site $i$ and year $j$) and I want to correlate year-to-year changes with the environmental conditions $x_{i,j}$. For this, I am fitting this model:

$$\begin{eqnarray} \mbox{log} ( \mu_{i,j+1} ) &=& \mbox{log} ( \mu_{i,j} ) + \alpha + \beta x_{i,j} + \gamma_j + \epsilon_{i,j} \\ \\ N_{i,j} &\sim& \mbox{Poiss} ( \mu_{i,j} ) \\ \gamma_{j} &\sim& \mbox{Norm} (0, \sigma ) \\ \epsilon_{i,j} &\sim& \mbox{Norm} (0, \theta ) \end{eqnarray} $$

So I'm interested in parameter $\beta$, the slope of the relationship. $\gamma_{j}$ is the random effect for year (as the residuals within single year were correlated) and $\epsilon_{i,j}$ is an overdispersion term.

The problem is in the year specific random effect $\gamma_{j}$. I think I have to use it, because the residuals are significantly explained by year. But, also the environmental conditions $x_{i,j}$ (climate conditions) are very much correlated within years! I.e. the climate in the same year is much more similar than accross years (as you would expect).

So the question is whether the introduced year-specific random effect $\gamma_{j}$ cannot "eat out" the variability that would be explained by the yearly variation of the climate (i.e. the term $\beta x_{i,j}$)? In case this problem occurs, the yearly variation of the climate would go to $\gamma_{j}$ instead of $\beta x_{i,j}$ and we could easily miss the significant relationship - i.e. the significant $\beta$ slope! Couldn't this happen? If yes, how to fit this model in such a way that the variability remains in $\beta x_{i,j}$ while at the same time the yearly autocorrelation in residuals is handled?

Answer 1

Since $\gamma_j$ is assumed to follow a zero-mean normal distribution, any deviation of the predicted value of $\gamma_j$ from zero will be penalized in the likelihood function relative to the variance $\sigma^2$. Thus it will be "cheaper" in terms of likelihood to put year-consistent variability into the fixed effect $\beta$.

The reason that it is cheaper to put the variability in a fixed effect rather than a random one is that the fixed effects will be estimated to minimize the residual, and there will not be a penalization if they are too big or too small. In contrast random effects are supposed to vary following a probability distribution (in your case a normal) around zero, and thus if you need an extreme value of the random effect to explain a given observation, the observation will become very unlikely given the model, which will give a lower likelihood value. If you can explain consistent variability in terms of fixed effects rather than random effects, you will increase the likelihood values. Thus, if you use maximum likelihood estimation, you will choose parameters that do that.

For example, consider the simple model

$$y_i=\alpha + x_i + \varepsilon_i,\qquad i = 1,\dots,n$$

where $(x_i)_i=\boldsymbol{x}\sim\mathcal{N}(0, S)$ and $(\varepsilon_i)_i\sim \mathcal{N}(0,\sigma^2\mathbb{I}_n)$. Let $\boldsymbol{y}=(y_i)_i$ and $\boldsymbol{\alpha}=(\alpha)_i$, $i=1,\dots, n$. The log-likelihood function is $$ \ell_{\boldsymbol y}(\alpha, \sigma^2, S) = -\frac{1}{2}\log\mathrm{det}(S+\sigma^2\mathbb{I})- \frac{1}{2}(\boldsymbol{y} - \boldsymbol\alpha)^\top(S+\sigma^2\mathbb{I})^{-1}(\boldsymbol{y} - \boldsymbol\alpha) $$ If we consider the part that depends on $\alpha$, namely the quadratic term, we can use some linear algebra to rewrite it $$ -(\boldsymbol{y} - \boldsymbol\alpha)^\top(S+\sigma^2\mathbb{I})^{-1}(\boldsymbol{y} - \boldsymbol\alpha) = -\frac{1}{\sigma^2}(\boldsymbol{y} - \boldsymbol\alpha)^\top(\boldsymbol{y} - \boldsymbol\alpha- S(S+\sigma^2\mathbb{I})^{-1}(\boldsymbol{y} - \boldsymbol\alpha)). $$ The term $S(S+\sigma^2\mathbb{I})^{-1}(\boldsymbol{y} - \boldsymbol\alpha)$ is the conditional expectation of $\boldsymbol{x}$ given the observation $\boldsymbol{y}$ which we usually denote $\mathrm{E}[\boldsymbol{x}|\boldsymbol{y}]$. This conditional expectation is in fact the best linear unbiased predictor of $\boldsymbol{x}$. With this in mind we can finally rewrite the square as $$ -\frac{1}{\sigma^2}(\boldsymbol{y} - \boldsymbol\alpha - \mathrm{E}[\boldsymbol{x}|\boldsymbol{y}])^\top(\boldsymbol{y} - \boldsymbol\alpha- \mathrm{E}[\boldsymbol{x}|\boldsymbol{y}])-\frac{1}{\sigma^2}\mathrm{E}[\boldsymbol{x}|\boldsymbol{y}]^\top S^{-1}\mathrm{E}[\boldsymbol{x}|\boldsymbol{y}]. $$ From this expression we see that in the first square of residuals, $\alpha$ and the predicted value of $\boldsymbol{x}$ plays a similar role, but the second penalizes deviation of the predicted value from zero. Thus maximum-likelihood estimation will always seek to describe as much variation as possible through the fixed effect $\alpha$. This holds in general.

Answer 2

If I understand your description correctly, I would say you are more likely to see a significant coefficient $\hat{\beta}$ by including a random effect. The reason is that with the introduction of $\gamma$, you now explicitly distinguish between-year variability and within-year variability. The overall variance in your data $V$ now has 2 parts

1. The variability from $\gamma$: $V_\gamma$
2. The variability from $\epsilon$: $V_\epsilon$

Without a random effect, $V_\gamma =0$ and all variability is attributed to $V_\epsilon$. The variance of the fixed effect estimator $\hat{\beta}$ is actually proportional to $V_\epsilon$. If you introduce random effects, $V_\epsilon$ will be smaller meaning $\hat{\beta}$ will have a short confidence interval (recall the confidence interval is inverse proportional to the square root of variance), thus a smaller pvalue.

Having said that, I am not sure why you want to include a random effect at the first place. My understanding is that the purpose of using random effect is to introduce correlation structure in your covariance matrix. In your case, $\gamma_j$ will make sure observations within a year are positively correlated. However, this is already achieved by using a time series model. If you extend the expression of $\log(\mu_{i,j})$ in term of $\epsilon_{i,1... j}$ (assuming no random effect) and calculate covariance between $\log(\mu_{i,j})$ and $\log(\mu_{i,k})$ for any $j \ne k$, you will find they are always positively correlated. That is the purpose of time series models. Random effect models are essentially doing the same thing, except for using some different techniques and being designed for a particular type of data (a few time points but multiple sequences).

Peter