Help making inferences from interactions within nested GLMM?

Question

I want to model survival outcome as a response to stem length for a plant species. All individuals within plots were tagged and sampled in year 1, then re-found and sampled again in year 2, with an interval of one year between measurements. The variables include fate at year 2 (survival = 1, death = 0) and length (normally distributed, continuous variable with min 0 and max 224). This data was collected from a total of 16 plots within two distinct populations (8 plots per population). Each row in the data includes data from both surveys: len is year 1 length, surv indicates that the plant was alive or dead at the time of re-sampling. Additionally, the experimental design included a single treatment variable with two levels: fenced plot ($n=7$) or unfenced plot ($n=9$). I am using length at the time of the first survey len to predict fate surv over the following year.

Data snippet:

 Population Treatment Plot  TagNum   flw   len flw_next len_next  surv  repr
 sc         fenced    nf11a     11     0     5        0     63       1     0
 sc         fenced    nf11a     25     0     5        0      0       0     0
 sc         unfenced  hp2a     549     0     6        0      0       0     0
 fl         fenced    mpna     413     0     7        0      0       0     0
 fl         unfenced  mpsb     920     0     7        0      0       0     0
 fl         unfenced  nnb      315     0     9        0     19.2     1     0

I would like to identify differences in the relationship between length and survival outcome at three levels: (H1) differences between populations, (H2) differences between treatments, and (H3) differences between plots within each population. I have specified the following two GLMMs using lme4:

Model 1:

mod.Surv <- glmer(surv ~ len*Population + (1|Plot:Population), family='binomial', data=plant_data)

Model 2:

mod.Surv.int <- glmer(surv ~ Treatment*Population + (1|Plot:Population), family='binomial', data = plant_data)

Model 2 is included only to make inferences about the effect of the interaction between treatment and population on survival. I don't have the stats background to really interrogate these models, so I'd like to know whether 1) these models are correctly specified and 2) how to interpret their output given my hypotheses 3)how to specify a single, "better" model that includes effects of all three nested variables

Output from summary(mod.Surv):

Generalized linear mixed model fit by maximum likelihood (Laplace  Approximation)
 [glmerMod]
 Family: binomial  ( logit )
Formula: surv ~ len * Population + (len | Plot:Population)
   Data: plant_data

     AIC      BIC   logLik deviance df.resid 
   356.0    388.9   -171.0    342.0      802 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-14.7529   0.0719   0.1636   0.3058   0.8438 

Random effects:
 Groups          Name        Variance Std.Dev. Corr 
 Plot:Population (Intercept) 0.789138 0.88833       
                 len         0.001607 0.04009  -0.87
Number of obs: 809, groups:  Plot:Population, 16

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)      -0.42173    0.69625  -0.606  0.54470    
len               0.10638    0.02905   3.662  0.00025 ***
Populationsc      1.18589    0.85981   1.379  0.16782    
len:Populationsc -0.04913    0.03292  -1.493  0.13555    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) len    Ppltns
len         -0.903              
Populatinsc -0.731  0.637       
len:Ppltnsc  0.692 -0.748 -0.877

Output of summary(mod.Surv.int):

Generalized linear mixed model fit by maximum likelihood (Laplace  Approximation)
 [glmerMod]
 Family: binomial  ( logit )
Formula: surv ~ Treatment * Population + (1 | Plot:Population)
   Data: plant_data

     AIC      BIC   logLik deviance df.resid 
   401.4    424.8   -195.7    391.4      804 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.7187  0.2172  0.2259  0.2941  0.3881 

Random effects:
 Groups          Name        Variance Std.Dev.
 Plot:Population (Intercept) 0.1129   0.336   
Number of obs: 809, groups:  Plot:Population, 16

Fixed effects:
                               Estimate Std. Error z value Pr(>|z|)    
(Intercept)                      3.0219     0.3799   7.955  1.8e-15 ***
Treatmentunfenced               -0.7428     0.4955  -1.499   0.1338    
Populationsc                    -0.6338     0.5219  -1.214   0.2246    
Treatmentunfenced:Populationsc   1.2821     0.6886   1.862   0.0626 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) Trtmnt Ppltns
Trtmntnfncd -0.752              
Populatinsc -0.718  0.544       
Trtmntnfn:P  0.555 -0.724 -0.763

Answer 1

If your understanding of the subject matter suggests that the association of one predictor with outcome depends on the value of a different predictor, then you should include an interaction between those predictors in the model. Consider that for each of the two-way interactions among length, treatment, and population, and for the 3-way interaction among all of them. Include all predictors and interactions in a single model, as binomial regression has an omitted-variable bias if any outcome-associated predictor is omitted.

In a comment, you suggest an additional model:

mod.Srv <- glmer(surv ~ len*Population + Treatment + (0+len|Plot), family='binomial', data=ribes_noNew)

That doesn't allow for the effect of Treatment to differ depending on either length or Population. Consider whether that makes sense. If it does, then the fixed effects are close to what you want.

You need to consider how to model the continuous predictor length. This model assumes that length is linearly associated with the log-odds of one-year survival. Things are seldom that simple. It's usually wise to fit such a continuous variable flexibly, for example with a regression spline.

For specifying the random effects in the model, consult this site's lmer cheat sheet and Ben Bolker's GLMM FAQ page.

The models in the original question have some problems. With only 2 levels of population, there isn't much to be gained by including it as a random effect, and it's already included as a fixed effect. Typically, random effects are most useful when there are on the order of 6 or more levels. Also, I don't think that the way you had specified Plot:Population accomplishes what you might have thought. If you think that the association of population with survival depends on the Plot, allow for a random slope of population among Plots.

In the model from your comment and copied above, the (0 + len|Plot) specification is risky. That doesn't allow for random intercepts among Plots, yet the intercept is the thing that is most likely to vary randomly among Plots: the intercept is the baseline log-odds of survival.

When the model is finally fitted to the data, it can be difficult to interpret the initial summary of a model with lots of interactions. In the default coding in R, each interaction coefficient is an estimated difference from what you would predict based solely on the lower-level coefficients. Furthermore, the value of a lower-level coefficient can depend on the coding of the predictors with which it interacts.

So don't worry about the individual coefficients on their own. Evaluate together all coefficients involving a predictor in a single test, for example with the Anova() function in the R car package. Let post-modeling tools like those in the emmeans package display predictions for particular predictor combinations of interest.