[R-sig-ME] random slopes model specification

[Ben Bolker]

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 of
> 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:
>
>>   coef(model6)
> $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 models
> 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  877911
> (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]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models