I'm building a series of relatively simple random effects models where we repeatedly measure a water quality variable, say conductivity (cond), in different watersheds (ws) and streams (st). Here, streams are nested within watersheds. We're not necessarily interested in the fixed intercept (mean conductivity) produced from these models; rather, we want to use the random coefficients from the models to make statements about how water quality varies across watersheds and streams. What I'm wrestling with is the best way to test the significance and correctly interpret the random watershed (ws) effect, given this sampling design.
We fit the model using lmer from the lme4 package. Note, I'm using ML instead of REML for simplicity.
m.ws.st = lmer(cond~1+(1|ws)+(1|ws:st),REML=F,data=data1)
This produces the following output for the random portion of the model. The values are small because they've been standardized.
Random effects: Groups Name Variance Std.Dev. ws:st (Intercept) 0.04240 0.2059 ws (Intercept) 0.03814 0.1953 Residual 0.25023 0.5002 Number of obs: 406, groups: ws:st, 14; ws, 4
It's pretty straightforward to test the watershed within stream effect by comparing the full model (m.ws.st) to a model with only a random watershed factor
m.ws = lmer(cond~1+(1|ws),REML=F,data=data1) pchisq(-2*logLik(m.ws)+2*logLik(m.ws.st),1,lower.tail=F)
This produces a tiny p-value, p = 8.8878e-08. So conductivity varies significantly between streams within watersheds, and the value s(watershed:stream) = 0.2059 estimates how far we expect a random stream to deviate from the mean conductivity of a random watershed.
For testing the significance of the random watershed effect, however, I imagine we can do this two ways:
- test the full model (m.ws.st) against a model with the just the random watershed:stream interaction (which would simply estimate the variation among streams)
- test the model with just the random watershed factor against a null, intercept only model.
First, we compare the full model against a model with just the random watershed:stream interaction
m.st = lmer(cond~1+(1|ws:st),REML=F,data=data1) pchisq(-2*logLik(m.st)+2*logLik(m.ws.st),1,lower.tail=F)
This model (w.st) produces the following output
Random effects: Groups Name Variance Std.Dev. ws:st (Intercept) 0.08034 0.2834 Residual 0.25023 0.5002 Number of obs: 406, groups: ws:st, 14
The model comparison test produces an insignificant p-value, p=0.1027. Variation among streams (0.08034) is equal to the variation among watersheds plus the variation among streams within watersheds (0.03814 + 0.04240). As part of my confusion, I'm not certain how I should interpret the test. I'd say we're testing whether or not variation among watersheds represents a significant portion of the variation among streams. The value s(watershed) = 0.1953 represents how far we expect a random watershed to deviate from overall mean conductivity.
Second, we compare the model with just the watershed random effect (m.ws, above) to a null, intercept only model
m0 = glm(cond~1,data=data1) pchisq(-2*logLik(m0)+2*logLik(m.ws),1,lower.tail=F)
This also produces a tiny p-value, p = 1.6038e-13. When I think of building the model by adding the watershed factor first and then adding the nested stream factor second, this type of forward step-wise approach makes sense. But the value for s(watershed) in the watershed only model (m.ws) partitions the variability in conductivity differently than the full model (m.ws.st); here s(watershed) is the deviation among watersheds across all streams. This model (m.ws) produces the following output
Random effects: Groups Name Variance Std.Dev. ws (Intercept) 0.05004 0.2237 Residual 0.28084 0.5299 Number of obs: 406, groups: ws, 4
As such, it doesn't exactly seem correct to list the s(watershed)=0.1953 value from the full (m.ws.st) model and then give it a p-value from the simplified (m.ws) model where s(watershed) = 0.2237. If I were to do this, though, I'd probably just report the values from the full model and then make the statement that there exists significant variation in conductivity among watersheds and then also among streams within watersheds. This approach doesn't allow us to compare the watershed and stream effects.
Given this model structure, do you have thoughts on the best way to test and interpret the watershed effect? Are both approaches appropriate depending on that question we're asking? What kind of language would you use to describe the effect?
Many thanks for your help! I know this is a bit confusing and quite a long thread.