drizopoulos / GLMMadaptive

GLMMs with adaptive Gaussian quadrature
https://drizopoulos.github.io/GLMMadaptive/
60 stars 14 forks source link

bug in nobs()? #10

Closed florianhartig closed 5 years ago

florianhartig commented 5 years ago

nobs returns for me always the value 10, irrespective of the real size of the data, e.g.

library(GLMMadaptive)
library(DHARMa)
testData = createData(sampleSize = 200, family = poisson())
fittedModel <- mixed_model(observedResponse ~ Environment1, random =  ~ 1 | group , family = "poisson", data = testData)
nobs(fittedModel)

I looked at the implementation and I'm not sure if I understand the sense of this - what's the function of the id slot? In any case, shouldn't this link to the model frame?

https://github.com/drizopoulos/GLMMadaptive/blob/ca5a891d329054f058f9304037cab3df5927e22c/R/methods.R#L1302

drizopoulos commented 5 years ago

The definition of the number of observations is more complex in the case of mixed models. In particular, you have $n$ the number of independent sample units that is determined by the number of levels in the grouping variables for which you specify random effects, and $N$ the total number of observations.

Because nobs() relates to the logLik() method and how the marginal BIC is calculated, it returns $n$. However, to make it more complete, I just committed an update which nobs() has an extra argument level. The default value 0 returns $n$, whereas if you set it to 1, it returns $N$.

florianhartig commented 5 years ago

hmm ... OK, thanks. I think it's definitely good to keep level = 0 as default, as this is how the other mixed model packages implement (try, e.g. lme4)

fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) , family = "poisson", data = testData)
nobs(fittedModel)

I see the point of being able to count RE units, but I wonder if nobs is the ideal function to do this.

Moreover, it's not intuitive for me how the level argument works or should work a mix of crossed / nested REs. I can see that I could go level 1,2,3 for nested REs, but for a mix of crossed / nested, something like re.form would probably work better.

But anyway, my immediate problem is solved.

florianhartig commented 5 years ago

sorry, I just noted the default is just opposite of what I thought. Shouldn't level = 0 go on the data itself (I would understand level = 0 as the data, not group level).

in any case, I would strongly recommend setting the default to data, to conform to the use in other packages.

drizopoulos commented 5 years ago

I changed the default to level = 1 to give the same behavior as other packages.