0

I'm using [lme4].

Data:

  tree rep trayid survival
1  Ash   1      1        1
2  Ash   1      1        1
3  Ash   1      1        1
4  Ash   1      1        1
5  Ash   1      1        1
6  Ash   1      1        1

Model:

NPVcorrbinary.glmer <- glmer(survival ~ tree + (1|tree:trayid), data=NPV01Datacorr.data, family="binomial", nAGQ=25)
Generalized linear mixed model fit by maximum likelihood (Adaptive Gauss-Hermite Quadrature, nAGQ = 25) ['glmerMod']
 Family: binomial  ( logit )
Formula: survival ~ tree + (1 | tree:trayid)
   Data: NPV01Datacorr.data

     AIC      BIC   logLik deviance df.resid 
  1081.4   1114.4   -533.7   1067.4      812 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.2942 -0.7744 -0.5773  0.9108  1.7508 

Random effects:
 Groups      Name        Variance Std.Dev.
 tree:trayid (Intercept) 0.01282  0.1132  
Number of obs: 819, groups:  tree:trayid, 36

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.1582     0.1850  -0.855 0.392390    
treeBC3F3     0.6006     0.2638   2.277 0.022779 *  
treeD54      -0.5170     0.2586  -2.000 0.045524 *  
treeD58      -0.9438     0.2799  -3.372 0.000747 ***
treeEllis1   -0.3869     0.2613  -1.481 0.138708    
treeQing      0.3817     0.2509   1.521 0.128177    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Question:

I'm trying to get the probability of mortality for each of the tree types.

> predict(NPVcorrbinary.glmer, data.frame(tree="Ash"), type="response")
Error in eval(expr, envir, enclos) : object 'trayid' not found
RegalPlatypus
  • 43
  • 9

1 Answers1

4

By default the predict() method produces predictions conditional on the random effects, and this is why it asks for the trayid. You can get predictions only using the fixed-effects part of the model using predict(..., re.form = NA).

However, note that these are not marginal/population averaged predictions. I.e., they are not the average survival probabilities of the specific group of trees. If you want such average probabilities, have a look at the GLMMadaptive package, and in the vignette Methods for MixMod Objects that describes how you can obtain coefficients with a marginal intepretation using the marginal_coefs() function, and the corresponding predictions using the predict() method.

Dimitris Rizopoulos
  • 17,519
  • 2
  • 16
  • 37
  • Thanks for the reply! Using re.form = NA works and gives me predictions that make sense for my data. A few more questions: In the interpretation of the coefficients in my model output as opposed to marginal coefficients you mention, would it be correct to say that the coefficients for my model only describe the sample in the data? I.e., only marginal/population averaged coefficients would provide predictions based on the population of trees represented by my trees in my data? If that's so, what does the standard error of the coefficients tell me in my model output? – RegalPlatypus Oct 12 '18 at 21:57
  • No, this is not exactly the issue; for more info check this question, the answer, and the comments to the answer: https://stats.stackexchange.com/questions/365907/interpretation-of-fixed-effects-from-mixed-effect-logistic-regression – Dimitris Rizopoulos Oct 13 '18 at 05:51
  • I took a look and I think I understand the difference. The SS interpretation, in my example, would be the survival of trees (fixed) that were reared in the same rearing tray (random/blocking)? I'm trying out marginal_coefs now. Every time I run it I'm getting slight variation in my std error/z-scores. That makes since since it's a monte carlo simulation. Is there a good way to decrease that variation? I tried doubling the M and K values with little effect (other than processing time) on the magnitude of the variation. As always, thank you so much! – RegalPlatypus Oct 14 '18 at 18:29
  • 1
    For now only increasing the number of Monte Carlo sample. I'm planning to implement better algorithms in the near future for decreasing the Monte Carlo variation. – Dimitris Rizopoulos Oct 15 '18 at 19:05