Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Calculating R squared for Poisson GLMM using MuMIn r.squaredGLMM

I am modeling abundance for a species of bird using a Poisson generalized mixed model using glmer in the R package "lme4". An example of my data:

abund      point_id  patch_area vis_per_year year
6      EL_03Plot035   244.69412          C_5 2003
0         RC_BBSM08   101.68909          C_2 2004
0  RP_211021_HSH088  1348.89935          C_3 2011
0         RC_LRSM04   111.74057          C_4 2008
0      RP_225155_p5     1.34007          C_3 2012
0       HO_YORUP105   141.66933          C_3 1998
1         RC_SPSM07   179.16088          C_2 2006
0     BH_MB12_bmh42 16937.30694         <NA> 2002
1         RC_MOSM11   104.43196          C_4 2012
1         RC_YOSM06   141.66933          C_4 2010
0  RP_244006_HMD366 27778.83482          C_3 2012
0      RP_247155_p5  7688.64751          C_3 2012
0      EL_08Plot127          NA          C_5 2008
2        HO_LITRR10   160.81729          C_4 1997
0         RC_BPSM07    38.23207          C_4 2009
0        HO_HARRIV5    10.46441          C_3 1999
1         RC_SPSM16   179.16088          C_4 2009
0         RC_YOSM01   141.66933          C_3 2002
0  RP_222799_HSH360    14.94866          C_3 2012
1         RC_WESM33   381.19813          C_2 2006
0  RP_209841_HSH017  2269.11227          C_3 2011
0         RC_LRSM03   111.74057          C_3 2001
0   RP_26718_HHO097    26.95666          C_3 2012
0     RP_236935_p14  7979.05373          C_3 2012
0 BD_miles_medium_2          NA          C_2 2003

"abund" is count data of max # of birds detected at that point each year, "point_id" is the name of the survey point, "patch_area" is the area of the habitat patch in which the point is located, vis_per_year is an ordered factor observation-level covariate indicating how many times a point was visited in a year, and "year" indicates year of the observation. With a few exceptions, there is only one abundance count (row) for each point per year.

My model specification is:

model=glmer(abund ~ scale(year) + (1|vis_per_year) + (1|patchid/point_id), family = poisson, data)

So far, for diagnostics I have checked my models for overdispersion using gof from the package "aods3", examined QQ plots of the random effects, and compared the models containing my fixed effect (I only have one) to a null model with only the random effect structure using the anova command from "lmerTest". The model is slightly underdispersed, and ranks higher than the null models using either AIC or anova criteria.

I am now trying to calculate R^2 for my final models. I read the blog entry at https://ecologyforacrowdedplanet.wordpress.com/2013/02/26/r-squared-for-mixed-models/ and the subsequent post update as well as the associated manuscript, and installed the MuMIn package to use for calculating marginal and conditional R^2. However, when I try to use r.squared.GLMM(Model), The following error is thrown:

Error in glmer(formula = SALS ~ scale(year) + (1 | vis_per_year) + (1 |  : 
  fitting model with the observation-level random effect term failed. Add the term manually

In addition: Warning message:

In t(mm[!is.na(ff), ]) :
  error in evaluating the argument 'x' in selecting a method for function 't': Error in mm[!is.na(ff), ] : (subscript) logical subscript too long

To isolate the source of the error I have tried running the model without point_id and visit_per year, however the same error is thrown without these covariates included. What exactly does this error mean - how do I add the observation term manually to the model? I have read the documentation for MuMIn, but I'm at a loss as to 1) exactly what the error means and 2) how to fix it. Any help is much appreciated. I don't think I can produce a reproducible example without providing the entire dataset, however knowing the exact meaning of this error will help me on my way.

Update:

At the recommendation of the error message and some explanation (thank you stack overflow!) I included my individual-level effect manually to the model. Now i get a different error:

Warning message:
In r.squaredGLMM.merMod(model) :
  exp(beta0) of 0.2 is too close to zero, estimate may be unreliable 

Does this mean the expected Beta for my fixed effect is too close to zero?

like image 890
Mo Correll Avatar asked Mar 17 '23 20:03

Mo Correll


2 Answers

I had the same problem with the same error message. After removing NA values it worked just fine for me. glmer is ok with NA-values but r.squaredGLMM seems to have a problem with that.

like image 90
Johannes Titz Avatar answered Mar 19 '23 16:03

Johannes Titz


The error means that after adding the individual-level random effect to the original model (see Nakagawa & Schielzeth 2013 paper), glmer gave an error (which has been converted to the warning message). You should add this term to your model by hand. It is of form (1 | IND) where IND is a factor variable with one level per observation (e.g. factor(1:nrow(data))).

However it may happen that after including the term the model won't converge (and it seems to be the case).

like image 33
Kamil Bartoń Avatar answered Mar 19 '23 16:03

Kamil Bartoń