Fixed effects panel models in ecology

Causal inference using fixed effects panel models

Causal inference includes a variety of approaches that seek to identify causal relationships between variables. The idea is to combine theory, logic, and robust statistical approaches to isolate a causal effect (i.e., how a change in x causes a change in y). There are many statistical approaches that can be used to infer causal relationships (see https://mixtape.scunning.com/ for a broader introduction into causal inference). Here I focus on a linear model that has often been used to estimate causal relationships – the fixed effects (FE) panel model.

FE panel models, often referred to as within-estimator models, leverage panel data to estimate causal relationships. Panel data are also referred to as longitudinal or time series data with time period ≥ 2. FE panel models can enable causal inference if certain conditions are met. FE panel models are increasingly applied to conservation and ecological studies, and various packages to estimate these models are being developed and revised in R; it’s a fast-moving area of development both in the literature and the various statistical software we use.

Broadly, the FE panel model approach divides all possible variables (i.e., both the observed and unobserved variables in a study) into time-varying and time-invariant. Simply put, time-varying means that the observed values change through time and time-invariant means that the observed values don’t change through time.

A simple FE panel model can be written as follows: \[ Y_{it} \sim \alpha + \beta X_{it} + c_i + y_t+ \epsilon_{it} \] Where Yit is the dependent variable (observed values from observational unit i at time t; for simplicity, from here on I use “plot” as the observational unit and “year” as the time period), α is the intercept, β is the coefficient of the time-varying explanatory variable (the hypothesized causal relationship of interest), ci is the heterogeneity associated with the unobserved, time-invariant variables unique to each plot, γt represents the heterogeneity associated with the unobserved, time-varying variables unique to each year and shared by all plots, and ɛit is the random error term (see https://doi.org/10.1111/2041-210X.13190).

Whether a variable is time-invariant or time-varying largely depends on the study system and spatial distribution of plots. Topography and soil type are typically considered time-invariant (unless of course a volcano erupts or major landslide occurs). The classification of most other variables depends on 1) the observed variable itself (i.e., species abundance can vary more than presence/absence), 2) the study design, and 3) the unique characteristics of the system of interest.

For example, in longitudinal studies with observations from plots located close together in space, there may be very little variation in rainfall across plots but high variation in rainfall across years. In this case, variation in rainfall is largely captured by year fixed effects. In contrast, if there’s high spatial variation across plots (e.g., the plots were established across large topographic gradients), then rainfall would vary by plot and would be estimated as a time-varying explanatory variable. Furthermore, it’s critical to consider whether there are unobserved, time-varying and time-invariant variables. For example, maybe there was an insect outbreak in a subset of plots that wasn’t captured by the observed variables. If and when these time-varying variables exist, additional tests can be done to determine the likelihood of bias due to their presence. These unobserved, time-varying variables are the major limitation of FE panel models and can lead to biased estimates of the causal relationship of interest.

The major advantage of an FE panel model is that unobserved, time-invariant variables are subtracted from the model, which can allow for more causal estimates of underlying relationships (if unobserved, time-varying variables are not biasing the estimates). Because it’s often impossible to control for all unobserved, time-varying variables in observational settings, coefficient estimates should always be carefully interpreted. Still, FE panel models offer an important step towards estimating causal relationships in ecology.

Centering and standardizing vs. demeaning

The FE panel model eliminates unobserved, time-invariant variables through a simple mathematical process of unit-level demeaning, using the within transformation (hence the alternative name to FE panel model – the within-estimator panel data model). Demeaning enables FE panel models to control for time-invariant measured and unmeasured variables, as well as time-varying variables that are collinear with time period (I’ll explain this further below), that may otherwise be confounding the underlying causal relationship a researcher is seeking to estimate. FE panel models essentially reduce the variation in the observed variables to deviations from plot- and year-level means and use this variation (which has now controlled for time-invariant, unobserved variable bias) to estimate the coefficients.

I’ll explain this demeaning process by comparing it to standardization of non-binary variables, often used prior to estimating statistical models in ecology (thereby increasing the interpretability of coefficients and enabling comparison of coefficient magnitude within and across models). The fundamental difference is that standardizing independent variables in ecology demeans (centers the data) across all of the observations, and often across all time periods, whereas demeaning in the FE panel model context subtracts out the mean from each individual- or plot-level observation. To demonstrate this, I create a simulated dataset with annual temperature and precipitation values measured at each plot across 20 years, as well plot-level elevation and year variables.

## creating a panel dataset with plot-level observations (n = 10 plots) across 20 years

##empty dataframe
paneldata = data.frame()

##dataset for loop
for (i in 1:20){
  
  elevation <- seq(1100,2000, by=100)
  plot <- seq(1:10)
  year <- as.integer(rep(paste(i), 10))
  richness <- sample(x = 1:16, replace = T, size  = 10)
  precip <- sample(x = 200:500, size  = 10)
  temps <- sample(x=3:15, replace = T, size = 10)
  
  mergedata <- data.frame(plot, year, elevation, richness, precip, temps)
  paneldata = rbind(paneldata, mergedata)
}

(Side note: economists often use “cross-sectional data” to describe a dataset with observations from one time period vs. “panel data” to describe a dataset that includes repeated observations from each observed unit over multiple periods of time.)

## A cross-sectional dataset example:

cross_sectional_dataset <- filter(paneldata, year==1)

Standardizing explanatory variables using the z-transformation first centers the data and then divides by the standard deviation.

center <- paneldata %>%
  mutate(temps_centered = temps-mean(temps), 
         temps_standardized = (temps-mean(temps))/sd(temps),
         elevation_centered = elevation-mean(elevation),
         elevation_standardized = (elevation-mean(elevation))/sd(elevation)) %>% 
  select(plot, year, elevation, temps_centered, elevation_centered, elevation_standardized)

head(center)

##   plot year elevation temps_centered elevation_centered elevation_standardized
## 1    1    1      1100          -0.34               -450             -1.5627772
## 2    2    1      1200          -1.34               -350             -1.2154934
## 3    3    1      1300          -5.34               -250             -0.8682096
## 4    4    1      1400          -4.34               -150             -0.5209257
## 5    5    1      1500          -4.34                -50             -0.1736419
## 6    6    1      1600          -6.34                 50              0.1736419

Notice the variation of all explanatory variables is largely maintained. This contrasts to demeaning in the FE panel model context, which seeks to eliminate the variation of time-invariant explanatory variables (elevation is a time-invariant explanatory variable).

To eliminate the variation associated with time-invariant variables, FE panel models can demean the data by centering variables at the plot- and year-level. While demeaning is often implemented using matrix algebra, to demonstrate the underlying idea, I’ll give an overly simplistic example. To demean at the plot-level, we first calculate the average richness, temperature, and precipitation for each plot. Then, we simply subtract out these means from the plot-level observations from each year – resulting in plot-level deviations across 20 years for all variables estimated in the model.

demeaningPlot <- paneldata %>%
  ##first demeaning by plot
  group_by(plot)%>%
  mutate(meanp = mean(precip), meant = mean(temps), 
         meanelev = mean(elevation), meanr = mean(richness)) %>%
  mutate(precip_demean = precip-meanp, temps_demean = temps-meant, 
        richness_demean=richness-meanr, elevation_demean=elevation-meanelev) %>% 
  select(plot, year, elevation, temps_demean, elevation_demean, richness_demean)
head(demeaningPlot)

## # A tibble: 6 × 6
## # Groups:   plot [6]
##    plot  year elevation temps_demean elevation_demean richness_demean
##   <int> <int>     <dbl>        <dbl>            <dbl>           <dbl>
## 1     1     1      1100       -0.100                0            4.9 
## 2     2     1      1200       -2                    0           -5.6 
## 3     3     1      1300       -6.05                 0           -2.5 
## 4     4     1      1400       -3.35                 0           -0.85
## 5     5     1      1500       -3.45                 0            6.05
## 6     6     1      1600       -6.95                 0           -1.55

## Notice that using the within-estimator demeaning, elevation is 0; it's time-invariant.

FE panel models also simultaneously demean variables from each time period. That is, all variables have a mean across all plots in a given year, so the resulting variables in the FE transformation, both the explanatory and outcome variables, are demeaned at the plot and year-levels.

In practice, there are at least three within transformations (for brief descriptions, refer to https://en.wikipedia.org/wiki/Fixed_effects_model). Depending on the package and estimator used, there may not be any demeaning prior to estimating the model, as estimating a plot-level coefficient estimate is mathematically equivalent to demeaning. However, particularly when using large datasets, different packages, like fixest, avoid estimating plot-level coefficients because it’s highly cumbersome for large datasets.

Fixed effects (FE) panel models

Simply speaking, FE panel models are often linear models that can be estimated using a variety of packages and software. I’ll use two estimators here for explanation and comparison: 1) the function feols from fixest package (see https://cran.r-project.org/web/packages/fixest/fixest.pdf), and 2) the function lm from R stats. I like using the fixest package as it’s one of the more developed packages for these models in R and more efficient when using large datasets, though it too is evolving rapidly. Both the lm and feols functions estimate models slightly differently, resulting in the same coefficient estimates with small differences in standard errors. (Side note: these models have very different R^2s, which I’ll ignore for now, since we’re less interested in model prediction and fit and much more interested in the underlying causal relationship of the explanatory variables and the outcome variable).

FE panel models using feols

## The FE panel model using the feols function from the fixest package

fe_mod <-  feols(richness ~ temps + precip | plot + year,
               data = paneldata)

summary(fe_mod)

## OLS estimation, Dep. Var.: richness
## Observations: 200 
## Fixed-effects: plot: 10,  year: 20
## Standard-errors: Clustered (plot) 
##         Estimate Std. Error   t value Pr(>|t|)    
## temps  -0.086097   0.103841 -0.829126 0.428474    
## precip -0.008461   0.003995 -2.117601 0.063284 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## RMSE: 4.42784     Adj. R2: 0.007905
##                 Within R2: 0.028072

The FE panel model I am using here estimates richness (the outcome variable), as a function of temperature and precipitation (i.e., the time-varying explanatory variables) and includes plot and year as the fixed effects (i.e., the dummy variables). Observed time-invariant variables (in this case elevation) are not included because time-invariant variables are dropped during demeaning (they are omitted anyway if included in the model (see side note below)).

## Side note: if you add in elevation, the time-invariant variable, you can see a warning 
## message that it's dropped from the model due to collinearity

fe_mod_elev <-  feols(richness ~ temps + precip + elevation  | plot + year,
               data = paneldata)

## The variable 'elevation' has been removed because of collinearity (see $collin.var).

The feols estimator is based off of Berge (2018). The core of this algorithm is based on a general ML model that “integrates a fixed-point acceleration method to hasten the convergence of the fixed-effect coefficients.” In addition, while there’s ongoing debate about whether and how to cluster standard errors – often used to control for bias associated with spatial correlation and experimental or observational designs that are clustered – in feols, standard errors are clustered by default using the fixed effects in the order they are entered into R. You can also use the vcov argument or other specifications from the package clubSandwich. In ecological settings, data are often nested within plot and/or spatially correlated, so clustering standard errors, by plot for example, should typically be included.

In this specific example, the FE model is essentially estimating deviation from plot- and year-level means of the explanatory variable (richness) as a function of the deviations of the time-varying explanatory variables (temperature and precipitation). The plot fixed effect is controlling for the time-invariant unobserved variables correlated with any variables in the model (i.e., both the explanatory and dependent variables), and the year fixed effect is controlling for variables that changed through time but only if the change was the same across plots (these could include atmospheric C02 fluctuations, solar radiation, etc). Time-varying variables that affect only certain plots, which might include weather variables, fire occurrence, defoliators, are NOT controlled for in the demeaning process.

FE panel models using lm

Below is an example of fixed effects panel model using lm with plot and year as factors (i.e., factoring the variables transforms them to dummy variables). The intercept is not shown to match the feols summary table; the intercept is difficult to estimate in large datasets and therefore often not reported. The ordering of the factors using lm is not important, since lm does not cluster standard errors by default.

## FE panel model using lm

fe_mod_lm <- lm(richness ~ temps + precip + factor(plot) + factor(year) -1,
             data = paneldata)

summary(fe_mod_lm)

## 
## Call:
## lm(formula = richness ~ temps + precip + factor(plot) + factor(year) - 
##     1, data = paneldata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1847  -3.5914   0.4511   3.2809   8.9024 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## temps          -0.086097   0.095325  -0.903   0.3677    
## precip         -0.008461   0.004360  -1.941   0.0540 .  
## factor(plot)1  10.311536   2.518060   4.095 6.53e-05 ***
## factor(plot)2  12.805263   2.517997   5.085 9.66e-07 ***
## factor(plot)3  12.345759   2.415237   5.112 8.57e-07 ***
## factor(plot)4  12.817174   2.443714   5.245 4.63e-07 ***
## factor(plot)5  12.454948   2.312967   5.385 2.40e-07 ***
## factor(plot)6  11.694272   2.499588   4.678 5.89e-06 ***
## factor(plot)7  11.780684   2.496909   4.718 4.97e-06 ***
## factor(plot)8  13.701533   2.367085   5.788 3.38e-08 ***
## factor(plot)9  11.816251   2.494369   4.737 4.57e-06 ***
## factor(plot)10 14.066950   2.443850   5.756 3.96e-08 ***
## factor(year)2   1.422366   2.157508   0.659   0.5106    
## factor(year)3   1.544412   2.157009   0.716   0.4750    
## factor(year)4   1.654800   2.161139   0.766   0.4449    
## factor(year)5   0.945296   2.164300   0.437   0.6628    
## factor(year)6   3.293183   2.155092   1.528   0.1284    
## factor(year)7   1.630927   2.169216   0.752   0.4532    
## factor(year)8   1.807964   2.161372   0.836   0.4041    
## factor(year)9   1.129838   2.186885   0.517   0.6061    
## factor(year)10 -1.950716   2.155674  -0.905   0.3668    
## factor(year)11 -1.821751   2.159952  -0.843   0.4002    
## factor(year)12  0.317525   2.162365   0.147   0.8834    
## factor(year)13  3.791207   2.160238   1.755   0.0811 .  
## factor(year)14  0.486082   2.157050   0.225   0.8220    
## factor(year)15  1.441576   2.157104   0.668   0.5049    
## factor(year)16  0.648251   2.172606   0.298   0.7658    
## factor(year)17  0.096238   2.167887   0.044   0.9646    
## factor(year)18 -0.757660   2.164398  -0.350   0.7267    
## factor(year)19  1.459348   2.168078   0.673   0.5018    
## factor(year)20 -0.104364   2.154683  -0.048   0.9614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.817 on 169 degrees of freedom
## Multiple R-squared:  0.8244, Adjusted R-squared:  0.7922 
## F-statistic: 25.59 on 31 and 169 DF,  p-value: < 2.2e-16

## combining these model coefficients to compare  them
fe_res <- data.frame(coeftable(fe_mod), model="feols summary")
lm_res <- data.frame(summary(fe_mod_lm)$coefficient[c(1:2),c(1:4)], model ="lm summary")

list(fe_res, lm_res)

## [[1]]
##            Estimate Std..Error    t.value   Pr...t..         model
## temps  -0.086097314  0.1038410 -0.8291263 0.42847427 feols summary
## precip -0.008460664  0.0039954 -2.1176014 0.06328369 feols summary
## 
## [[2]]
##            Estimate  Std..Error    t.value   Pr...t..      model
## temps  -0.086097314 0.095324662 -0.9032008 0.36770490 lm summary
## precip -0.008460664 0.004359554 -1.9407179 0.05395646 lm summary

You can see that the lm and foels estimators produce the same coefficient estimates, but with slightly different standard errors and corresponding t-values and p-values. The reason that they are producing the same coefficient estimates is that mathematically, the plot and year coefficients are the same as if these models were demeaned first before being estimated. The reason they produce different standard errors is due to the underlying estimator engine and the default clustering of the standard errors in feols. In theory, both lm (and alternatively glm) and feols estimations are correct; therefore it’s incumbent on the researcher to decide which estimator to use both in R and across other statistical software (often comes down to efficiency and familiarity).

Difference between FE panel models and glmms

Now I will compare fixed effects panel models (using feols) to mixed models, specifically generalized linear mixed models (glmms) that are widely used in ecology. Notice I include elevation as an explanatory variable to demonstrate that glmms estimate both time-invariant and time-varying variables when included as covariates, which contrasts to the FE panel model using feols that drops time-invariant variables through demeaning.

Side note and warning: the following explanation may cause further confusion than good. The differences between mixed models in ecology and fixed effects models in econ can seem confusing due in part to the fact that fixed effects are defined very differently. For example, a fixed effect in a glmm is more equivalent to an explanatory variable in an FE panel model than a fixed effect in a FE panel model, whereas a random effect in a glmm would more likely be estimated as a fixed effect in a panel model (for more detail https://en.wikipedia.org/wiki/Fixed_effects_model).

## GLMM: A glmm using random effects for year and plot
glmm_res_A <- lmer(richness~ temps + precip + elevation + (1|plot) + (1|year),
              data = paneldata)

## boundary (singular) fit: see ?isSingular

summary(glmm_res_A)

## Linear mixed model fit by REML ['lmerMod']
## Formula: richness ~ temps + precip + elevation + (1 | plot) + (1 | year)
##    Data: paneldata
## 
## REML criterion at convergence: 1212.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1658 -0.7560  0.1085  0.8207  1.7819 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  year     (Intercept) 0.000e+00 0.000e+00
##  plot     (Intercept) 8.183e-15 9.046e-08
##  Residual             2.276e+01 4.771e+00
## Number of obs: 200, groups:  year, 20; plot, 10
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 10.347513   2.488097   4.159
## temps       -0.096800   0.088742  -1.091
## precip      -0.007959   0.003964  -2.008
## elevation    0.001809   0.001176   1.538
## 
## Correlation of Fixed Effects:
##           (Intr) temps  precip
## temps     -0.289              
## precip    -0.570 -0.100       
## elevation -0.765  0.018  0.046
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular

I’ll focus on the differences between FE panel models and glmms relevant to causal inference. First, FE panel models control for all unobserved time-invariant variables. Though the random effects in glmms help control for much of the unobserved variable bias (for example, the unobserved variables correlated with plot and year), they cannot control for the possibility that these unobserved variables are also correlated with the explanatory variables (see https://doi.org/10.1111/2041-210X.13190). For instance, unobserved time invariant variables (let’s use soil type) may be correlated both with plot AND the explanatory variables, precip and temperature. Including the random effect of plot only controls for the correlation between soil type and plot, not the correlation between soil type and precip and soil type and temp. Thus, in the random effects glmm model above, correlations between unobserved variables and the explanatory variables are captured in the error term, resulting in biased coefficient estimates of precip and temp.

Second, glmms cannot control for severe multicollinearity among the observed explanatory variables. In contrast, FE panel models can control for severe multicollinearity induced by time-invariant variables. If temperature and precipitation are highly correlated with elevation, for example, elevation would likely have a very high variance inflation factor (VIF) https://www.tandfonline.com/doi/abs/10.1081/QEN-120001878. The higher the VIF, the larger the SEs, resulting in less precise and interpretable coefficient estimates. If there’s a strong correlation between a time-invariant (e.g., elevation) and time-varying variable (e.g., temperature), FE panel models provide a useful solution, as they subtract out the time-invariant variable, thereby estimating a more precise coefficient of the time-varying variable. This is of course much preferable to omitting the problematic variable when the VIF is too high, as it increases the interpretability and precision of the coefficient estimate.

Finally, one of the main advantages of FE panel models, other than they often provide less biased coefficients of explanatory variables than a glmm, is that it can be used to estimate a counterfactual scenario – a world where the treatment didn’t occur. This counterfactual can be equated to a control treatment in an experiment. Because the coefficient estimates in a FE panel model are capturing less biased estimates of the causal relationship of interest, they can be used to estimate a world in which one or more of these explanatory variables didn’t change through time. This counterfactual thus resembles a control in an experimental setting and the difference between control and treatment can be estimated.

FE panel model limitations

There are two major limitations of FE panel models. First, an FE panel model is fairly useless when predicting out-of-sample. If prediction is the main focus of an analysis, using a random effects model like a glmm that estimates the variance associated with plots can facilitate out-of-sample predictions (with the caveat that these predictions are likely biased). Second, FE panel models cannot control for unobserved, time-varying variables, which in some settings, means they are likely less useful than a glm(m).