0

I am attempting to do a goodness of fit test for a several competing GLMM models. First, it is unclear to me what the best methods are for this, as there are many different methods with positive and negative aspects. So, I came across the package "rms" which provides both the R-squared and a likelihood Ratio Test for X Square. Problem is, I am not sure if rms, command lrm, is applicable for GLMER and how to interpret the results I have obtained if they are. The X square value seems high to me, but perhaps I am not understanding something.

So, to be clear, my question is if this method for representing goodness of fit is acceptable. If not, what are some recommendations to go about this? If yes, then how am I to understand that chi-square value in my ouput?

If there are comments about the rest of my code, I kindly accept them!

####Sample code### (will provide different ouputs from what I have below)

Factor Patch.Structure PatchGram PatchForbs Microtopography standing.dead
1   4   86.6667 5.4444  2   100
0   3   34.3333 5.4555  0   50
1   3   58.3333 2.7778  2   100
0   3   32.3333 19.6667 1   100
1   2   42.7778 18.8889 0   0
0   4   31.6667 13.8889 1   100
1   4   34.4444 16.3333 1   100
0   4   39.4444 7.4444  1   100
1   4   9.6667  5.6667  2   5
0   7   14.8889 75.5556 0   10
1   4   20.0000 19.7778 1   100
0   4   25.0000 21.6667 1   100
1   4   30.3333 21.8744 0   45
0   4   9.7778  0.5556  1   90



###Multivariate models####
Model1 = glmer(Factor~Patch.Structure+(1|Study.area)+(1|Study.Year),      data=NestLand, family=poisson (link="log"),nAGQ = 0)
Model2 = glmer(Factor~PatchGram+(1|Study.area)+(1|Study.Year), data=NestLand, family=poisson (link="log"),nAGQ = 0)
Model3 = glmer(Factor~PatchGram+standing.dead+(1|Study.area)+(1|Study.Year), data=NestLand, family=poisson (link="log"),nAGQ = 0)
Model4 = glmer(Factor~PatchGram+Microtopography+(1|Study.area)+(1|Study.Year), data=NestLand, family=poisson (link="log"),nAGQ = 0)
Model5 = glmer(Factor~PatchGram+PatchForbs+(1|Study.area)+(1|Study.Year), data=NestLand, family=poisson (link="log"),nAGQ = 0)
Model6 = glmer(Factor~Patch.Structure+PatchGram+(1|Study.area)+(1|Study.Year), data=NestLand, family=poisson (link="log"),nAGQ = 0)

###Outputs##

out.put<-model.sel(Model1, Model2, Model3, Model4, Model5, Model6) 
out.put

####Model fit for top model
Require(rms)
lrm(Model6)

####Output
Logistic Regression Model

 lrm(formula = Model6)

                  Model Likelihood     Discrimination    Rank Discrim.    
                     Ratio Test           Indexes           Indexes       
Obs         238    LR chi2     128.64    R2       0.650    C       0.932    
0           188    d.f.            14    g        6.758    Dxy     0.863    
1            50    Pr(> chi2) <0.0001    gr     861.347    gamma   0.864    
max |deriv| 0.02                         gp       0.290    tau-a   0.288    
                                       Brier    0.078                     

                Coef     S.E.     Wald Z Pr(>|Z|)
 Intercept            7.3794  80.0725  0.09  0.9266  
 Patch.Structure=2   14.4677  50.0347  0.29  0.7725  
 Patch.Structure=3   11.5718  50.0199  0.23  0.8170  
 Patch.Structure=4   12.9510  50.0204  0.26  0.7957  
 Patch.Structure=5   11.9745  50.0347  0.24  0.8109  
 Patch.Structure=6    9.1426  50.0285  0.18  0.8550  
 Patch.Structure=7   11.2092  50.0274  0.22  0.8227  
 Patch.Structure=8    2.6170  74.3649  0.04  0.9719  
 Patch.Structure=9    1.9226  58.8724  0.03  0.9739  
 Patch.Structure=11   4.1922 110.6873  0.04  0.9698  
 PatchGram           -0.0505   0.0118 -4.27  <0.0001 
 Study.area=IR       -1.8942   0.7928 -2.39  0.0169  
 Study.area=NF        0.0446   0.8813  0.05  0.9596  
 Study.Year=2016    -18.8128  62.5319 -0.30  0.7635  
 Study.Year=2017    -17.2469  62.5325 -0.28  0.7827  
  • Yes, of course. I am trying to get the relative goodness of fit. How would I go about migrating this question to CrossValidate? Just reask it in that forum? On a side note Ben, I have often reference your book Ecological Models and Data in R, so, Thank you! – user5613688 Dec 17 '17 at 02:18

1 Answers1

1
  • lrm is a fine substitute for glm(...,family="binomial"). The models it fits are nearly equivalent, with some extra bells and whistles/extensions (penalized likelihood, ordinal models); the output format is slightly different. AFAIK it does not handle mixed models.
  • you may need to be more specific about what your goal is in testing goodness of fit. There are two primary paradigms for assessing the relative goodness of fit of a series of (G)LMMs:
    • likelihood ratio tests (anova(model1,model2,...) in R) perform pairwise tests of a series of nested models, using the null distribution of differences in deviance ("LR chi2" in the lrm output, or "deviance" in summary() output from a glmer fit) to test for statistically significant differences between pairs.
    • information-theoretic approaches (AIC, BIC, etc.) attempt to find/estimate the best predictive model (minimized Kullback-Leibler distance or something else) among a series of not-necessarily-nested models (see here for more information/discussion)

Neither of these approaches tests whether the most complex model is an adequate fit to the data; if not, then none of the results will be sensible. The magnitude of the AIC or residual deviance/$\chi^2$ for a single model usually doesn't tell you very much (it depends most heavily on the size of your data set). Calculating absolute GOF is not always straightforward; there are some tests (e.g. Hosmer-Lemeshow-le Cessie test, via residuals.lrm(fitted,"gof"), or you can do posterior predictive checks, or you can simply examine your predicted values to see that they seem adequate.

Ben Bolker
  • 34,308
  • 2
  • 93
  • 126
  • 1
    To implement `rms` methods for `glmer` models would require you to write these methods. Contributions are welcomed, and `rms` is managed on `github`. – Frank Harrell Dec 17 '17 at 16:20