Setup: We have 10 connected but distinct wetlands within a study area. The wetlands are not totally independent, there is some exchange of water and organisms between them. Perhaps comparable to states or provinces - sharing borders and with some connections but largely self-governing.

We have collected 21 years of annual data. Each year 65 sites are chosen randomly from within the entire study area, resulting in 2-16 observations within each wetland (panel) per year.

Can I treat each of the 10 wetlands as a "panel" to identify the strongest potential drivers of change over time? I am using the PLM function in R which wont allow varying numbers of observations per time step...at least I can't figure out how to make that happen. ?? Is there a way of using the individual observations within each panel in PLM? Another option?

My main goal is to determine which of 3 predictors is driving a very strong overall decrease in CLARITY over time. This decrease and changes in the predictors differ between panels (see attached figure for an example).

The basic code is:

pSummer <- pdata.frame(Summer, index = c("complex", "Year"))

fixed <- plm (CLARITY ~ VEG + FISH + Inputs, data = pSummer, model = "within")

When I average values within each panel (so each panel contains only 1 "average" individual) I get results that make perfect sense, but I am losing or ignoring a lot of potentially important variation. I know it's not correct, but here are the (partial) results:

Unbalanced Panel: n = 10, T = 17-21, N = 206

Coefficients: Estimate Std. Error t-value Pr(>|t|)
FISH 0.0016 0.0002 6.57 4.398e-10 ***

Inputs 0.6401 0.2793 2.29 0.0230 *

VEG -2.6780 1.3997 -1.91 0.0572 .

If I try to include all observations I get the following error message, which is clear enough.

In pdata.frame(Summer, index = c("Subname", "Year")) : duplicate couples (id-time) in resulting pdata.frame

FYI as an example, here are trends in the VEG predictor within each panel.  Dots are annual average values

Thanks for any ideas.


It is not statistically wrong to estimate the model you have in mind if you think it is realistic. Remember, the "within" estimation for panel data just means you are assuming there is a separate intercept for each group of observations with the same id variable. If it makes sense in your context, there is no reason not to assume the baseline clarity should be the same for different sites within a complex. However, you cannot estimate this model with PLM. Instead, estimate the model using OLS with dummy variables:

fixed <- lm (formula= CLARITY ~ VEG + FISH + Inputs + factor(complex) -1, data = pSummer)

  • $\begingroup$ Thank you. The model I suggest is the most realistic, logical option that I have been able to produce - drivers change over time at different rates in each of the panels. I followed a set of comparisons (youtube.com/watch?v=1pST2lUx6QM) that suggested the fixed effect approach, and it would make sense to me that intercept would vary by panel. It seems to me that slope should also vary by panel because the drivers vary considerably by panel (see figure in the question email). I have some new things to try in any case. Thank you again. $\endgroup$ – DCD Apr 9 '19 at 17:36

This sounds like a case where a mixed effect model, with a random effect for wetland, would be appropriate. Random effect groupings can have multiple observations per grouping level, and are straightforwardly estimated using software like lme4.

A model like:

random_lme4 <-   lmer(CLARITY ~ VEG + FISH + Inputs + (1|complex:year) + (1|complex), 
                      data = Summer)

would estimate a single linear relationship for vegetation, fish, and inputs with clarity across all lakes, and a separate intercept for each wetland (the (1|complex) term) and for each wetland in each year (the (1|complex:year) term). While I'm not an expert on panel models, from what I understand from the plm vignette this would be similar to a model="random" plm model. The model="within" version, where the relationship between each predictor term and the outcome variable is assumed to vary between wetlands, can be modelled with:

fixed_lme4 <- lmList(CLARITY ~ VEG + FISH + Inputs|complex:year, 
                     data = Summer)

Although in that case, I would really recommend a varying slopes random effect model, like:

random_slopes_lme4 <- lmList(CLARITY ~ VEG + FISH + Inputs+ (1|complex:year) + (1+VEG + FISH + Inputs|complex), 
                     data = Summer)

That will tend to somewhat regularize estimates by pulling the estimates for the relationships between each predictor and the outcome in each wetland towards the estimates for the other wetlands.

  • $\begingroup$ Hi Eric- Thank you this gives me a lot to think about. I did a mixed model with complex = random, initially just to convince myself that the panels were distinct and accounted for some of the overall variation in clarity, which they did. However, this approach seems to yield a different overall answer (that makes less sense) than than the panel approach. I will try some permutations of what you suggest here and may contact you through our mutual colleague. Thank you again for your thoughts. $\endgroup$ – DCD Apr 9 '19 at 17:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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