Dear R-sig-mixed-models,

I am doing an analysis of how seedling mortality or growth relates to a few continuous topographic variables, and whether species might respond differently to the topography.  I hope someone might be able to help confirm the correct specification and interpretation of the random effects in a GLMM version of this analysis, which has been suggested by a handling editor on a paper.  I have some experience with these methods, but am unsure in this case whether I am proceeding correctly.  Your advice and/or suggestions for references that would give some guidance would be greatly appreciated. Beyond proper lme4 syntax, my concerns include interpretation of zero or near-zero variances for the random effects, whether to test for spatial autocorrelation in the residuals, and interpreting differences among species when species is a random factor.

I am using glmer() and a logit link for the survival model, including fixed factors of 3 topographic models ("topo1", "topo2", and "topo3" for simplicity), starting height ("ht") I have 130+ species ("sp") found at 200 census stations ("station").  Not all species are found at all stations, and the sample size per species ranges from 10 - 1200 individuals (and I could restrict these further to ones with a sample size greater than some threshold).  I would like to know whether the topographic variables are significant predictors of mortality while including the random factors of census station to account for non-independence of seedlings at the same location (which have the same topo measurements) and species to allow for variation in species' responses.  I expect that both the slope and intercept of species' responses to each variable could be quite different.  To allow for different slopes/intercepts among species, I have centered the continuous variables and specified the model as:

glmer(survival ~ topo1 + topo2 + topo3 + ht + (0 + topo1 | sp) + (0 + topo2 | sp) + (0 + topo3 | sp) + (1 | sp) + (1 | station), data=seedlingdata, family=binomial)

Questions:
When I do this, there is a random intercept for station, a random intercept for species, and then random slopes among species for the relationship with the topographic variables as follows in the model output.  I believe this is allowing for the variation among species that I intend, but would like confirmation of this specification vs. something like (topo1  |  sp) or (1 + topo1 | sp) as someone else has suggested to me.

Any version of these models that I have run results in significant fixed factors and zero or near-zero variances for the random effects.  I interpret this to mean that the topographic variables are important predictors of seedling mortality, but that the relationship does not vary among species groups nor census locations.  Is this your interpretation too or need I worry about model specification or the sample size or variance structure of my variables?

A suggestion was made to confirm a lack of spatial autocorrelation in the residuals of this model, but I am not sure that is appropriate given the inclusion of the random effect of census station and the fixed effects of topography, which are shared by seedlings at the same station.  Can anyone suggest an appropriate reference to support or refute this suggestion?

Finally, if the response to topography DID significantly vary among species, where in this model would I see it?  In a large variance for the species slopes or intercept?  Or would I need to include species as a fixed factor crossed with the topographic variables?  

Thanks very much for your advice.  A sample of the data and the model output follow:

Best,
Margaret Metz


seedlingdata subset:

      sp station   ht survival       topo1      topo2     topo3
  oenoba       1  3.0        D    19.02738  0.2123052  1.962315
  pipeau       1  2.0        D    19.02738  0.2123052  1.962315
  pipeo1       1  8.0        D    19.02738  0.2123052  1.962315
  pipeo1       1  2.5        D    19.02738  0.2123052  1.962315
  pipeo1       1  5.0        D    19.02738  0.2123052  1.962315
  ingaau       1 10.0        A    19.02738  0.2123052  1.962315
  ingaau       2  6.5        D    18.09238 -2.8782072  2.130231
  oenoba       2  6.0        A    18.09238 -2.8782072  2.130231
  pipeo1       3 16.3        D    15.41238  4.8228318  1.043460
  rinoli       3  6.5        A    15.41238  4.8228318  1.043460


Model output

Generalized linear mixed model fit by the Laplace approximation 
Formula: survival ~ topo1 + topo2 + topo3 + ht + (0 + topo1 | sp) + (0 + topo2 |      sp) + (0 + topo3 | sp) + (1 | sp) + (1 | station) 
   Data: seedlingdata 
  AIC  BIC logLik deviance
 9745 9815  -4862     9725
Random effects:
 Groups Name              Variance   Std.Dev.  
 station(Intercept)       1.5394e-14 1.2407e-07
 sp     (Intercept)       1.2477e-14 1.1170e-07
 sp     topo1             5.8720e-16 2.4232e-08
 sp     topo2             4.3709e-13 6.6113e-07
 sp     topo3             0.0000e+00 0.0000e+00
Number of obs: 8109, groups: station, 200; sp, 135

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)       -0.57771    0.02468 -23.406  < 2e-16 ***
topo1		   0.15638    0.05420   2.885  0.00391 ** 
topo2	           0.23070    0.03136   7.356  1.9e-13 ***
topo3	          -0.05046    0.04704  -1.073  0.28339    
ht	           1.21482    0.05085  23.889  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Correlation of Fixed Effects:
            (Intr) topo1     topo2   topo3
topo1	    -0.032                           
topo2       -0.034 -0.571                    
topo3        0.012 -0.839    0.342           
ht           0.003 -0.054    0.032   -0.023  



	[[alternative HTML version deleted]]


