On Wed, Apr 18, 2018 at 3:04 PM, Karista Hudelson <karistaeh at gmail.com>
wrote:
Hello Mixed Modelers,
I was hoping to get a bit of feedback on the random effects structure of
the model below. I wanted to allow two of the fixed effects (my effects
interest) to have random slopes and also random intercepts for each
watershed (there are 5 watersheds, 27 observations of Spring_MST &
Summer_Rain, and 790 fish (so there are 790 measurements of Hg and
FishLength).
model6<-lmer(Hg~FishLength+Spring_MST+Summer_Rain+(1+Spring_
MST+Summer_Rain|watershed),data=FSV,REML=FALSE)
this model converges, and the output looks "like I want it to look"; ie,
slopes are different for the two effects which logically should have
different slopes:
$watershed
(Intercept) FishLength Spring_MST Summer_Rain
877911 -1.7734565 0.02147013 0.03245297 -0.04841394
3319324 -1.2429440 0.02147013 0.03037712 -0.04472220
11597701 -1.6688121 0.02147013 0.03204351 -0.04768574
31121323 -0.9214729 0.02147013 0.02911923 -0.04248513
97032352 -1.4819562 0.02147013 0.03131236 -0.04638544
Notice that there is not much variability in your coefficients
across watersheds
(zero variability in FishLength effects because you left it out of your
random
effects on purpose, about 10% variability in Spring_MST and Summer_Rain).
Normally it's not recommended to fit a model with three correlated RE terms
(intercept, spring_MST, summer_rain) to only a grouping factor with
only 5 levels -- this is
equivalent to trying to estimate a 3x3 variance-covariance matrix with
5 observations
(I'm a little surprised this worked at all). Are all the values in
getME(fitted_model,"theta") 'reasonable',
i.e. larger than (say) 1e-4?
There is a huge, ongoing debate about what to do in this case,
where the RE model
you want is theoretically identifiable (you do have multiple
measurements of spring_MST
and summer_RAIN in each watershed), but practically very poorly
constrained, suggestions ranging from
"keep it maximal" (Barr et al): fit the full model, reduce when the
model is singular
try to constrain complexity _a priori_, reduce by model selection
(AIC or p-value based)
use Bayesian methods to add priors/regularize the problem
substitute fixed effects instead.
Here is a bit of the summary() for model6:
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.417728 0.144246 6.000000 -9.829 6.18e-05 ***
Length 0.021470 0.001085 783.800000 19.794 < 2e-16 ***
Spring_MST 0.031061 0.007518 312.000000 4.132 4.63e-05 ***
Summer_Rain -0.045938 0.007276 132.100000 -6.314 3.82e-09 ***
It seems to be doing what I intended, but the random effects structure
looks like a bit of a run-on sentence compared to most of the other
I've seen. Does this appear correct?
Also....if I'm allowing for the slopes of the effects of interest to be
different for each watershed and allowing intercepts to be different, in
your opinion(s) would it be easier to understand if I just run separate
models for each lake? Here how that looks for the Lake with watershed
(corresponding to top line of coef() table above):
clim<-Hg~FishLength+Spring_MST+Summer_Rain
climS<-lm(clim,data=Smid)
summary(climS)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.3531180 0.2152754 -6.286 2.19e-09 ***
Length 0.0276015 0.0027818 9.922 < 2e-16 ***
Spring_MST 0.0247300 0.0081347 3.040 0.002699 **
Summer_Rain -0.0025433 0.0007216 -3.525 0.000531 ***
Any thoughts are appreciated. Thanks in advance for your time and brain
power.
Fitting separate models will be almost equivalent to fitting a single
Hg~(FishLength+Spring_MST+Summer_Rain):watershed
model -- the only difference will be whether the residual variance is
pooled/assumed to be
the same across watersheds.
Much Obliged,
Karista Hudelson
[[alternative HTML version deleted]]