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

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?

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.

• Lars, thank you for your answer! I have a few questions: 1) why is it "cheaper to put the variability in $\beta$?" 2) Regarding the size of $\sigma^2$ vs $\beta_ix_{i,j}$, how do I compare it properly? Unfortunatelly there is a convergence problem and I have only one model run where estimate of $\sigma$ converged. The estimate is $\sigma = 0.17$, $\beta = -0.47$, $sd(x_{i,j}) = 0.06$ (should I compute sd of $x$? mean is zero), so this gives me $sd(\beta_ix_{i,j}) = 0.03$ which is much less than the $\sigma$! Does it mean that the variability is then consumed by the random effect? Commented Oct 2, 2014 at 14:21
• Note to point 2) above - I am not comparing "variance of $\sigma^2$" as you suggested but instead I am comparing $sd(\beta_ix_{i,j}) = 0.03$ with $\sigma$. I hope this is what you meant? Commented Oct 2, 2014 at 14:25
• I have replaced $\gamma_i$ with $\gamma_j$ in your answer - I guess this is what you meant. Commented Oct 2, 2014 at 14:27
• @Curious: I have updated the answer. You should compare $\sigma$ to the actual values of $\beta_i x_{i,j}$ and not their variability. If a typical value of $\beta_i x_{i,j}$ is 0.00001 and the estimated standard deviation of your random effect is 0.17, it seems that $\beta$ does not explain a lot of variation. Commented Oct 3, 2014 at 7:57
• Lars, but the values of $x_{i,j}$ have approx normal distribution with the mean = zero, so what is the "actual value of $\beta_i x_{i,j}$" then? Mean of this term is actually zero. Commented Oct 3, 2014 at 9:13

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

• Peter, thank you for the answer! I have to go so I will have to look at your answer in more detail later. Now just a small comment - $\epsilon$ is not the only residual variability. $\epsilon$ is only the overdispersion - main residual variability is present within the $N_{i,j} \sim \mbox{Poiss} (\mu_{i,j})$ term. Commented Oct 2, 2014 at 16:43
• And you say $V$ is the "overall variance", don't you mean "overall residual variance" instead? Commented Oct 2, 2014 at 16:44
• Peter, 1) I am using random effects because the residuals of the model where correlated within years. This means there is some part of between-year population growth present over all sites which is not explained by $x_{i,j}$. This is not handled by the time series model. You speak of autocorrelation on the population size, while I speak of the autocorrelation of population growth, this is a different thing. 2) Sometimes it can happen that adding random effect will in fact enlarge the confidence intervals. Commented Oct 3, 2014 at 6:13
• I agree with comment 1). In fact, I should read your index notation more carefully. Regarding your comment 2), when you fit your model, is $\sigma$ significantly larger than zero? Sometimes adding a random/fixed effect not strongly associated with outcome will add more noise rather than explaining variability in outcome. Commented Oct 3, 2014 at 11:51
• "is $\sigma$ significantly larger than zero?" Peter, well, I think $\sigma$ is always significantly larger than zero because it simply can't be negative... Commented Oct 9, 2014 at 9:12