1

I need to rank candidate model using the QAIC criterion. This is what I've tried:

library(MuMIn)
model_global <-  glm(vs ~ mpg + disp + wt, family = quasibinomial, mtcars)
model_1 <-  glm(vs ~ mpg, family = quasibinomial, mtcars)
model_2 <-  glm(vs ~ disp, family = quasibinomial, mtcars)
model_3 <-  glm(vs ~ wt, family = quasibinomial, mtcars)
model_null <-  glm(vs ~ 1, family = quasibinomial, mtcars)

mod.sel(model_global, model_1, model_2, model_3, model_null, rank="QAIC", chat=deviance(model_global) / df.residual(model_global))

This returns this error:

Error in formula.default(x) : invalid formula

How can I rank the above models using QAIC?

www
  • 38,575
  • 12
  • 48
  • 84
luciano
  • 13,158
  • 36
  • 90
  • 130

2 Answers2

3

You need to provide your rank.args as a list (see ?model.sel)

library(MuMIn)
model.sel(model_global, model_1, model_2, model_3, model_null,
          rank = QAIC,
          rank.args = list(chat = deviance(model_global) / df.residual(model_global)))

In the (nice little) example you provided, chat happens to be < 1, a warning is generated, and chat is set to 1.

# Model selection table 
#              (Intrc)  disp    mpg    wt     df logLik  QAIC delta weight
# model_global -21.9000 -0.0403 0.6470  5.332 4   -8.844 25.7  0.00 0.569 
# model_2        4.1380 -0.0216               2  -11.348 26.7  1.01 0.344 
# model_1       -8.8330         0.4304        2  -12.767 29.5  3.85 0.083 
# model_3        5.7150                -1.911 2  -15.683 35.4  9.68 0.004 
# model_null    -0.2513                       1  -21.930 45.9 20.17 0.000 
# Warning messages:
#   1: In rank(x, chat = 0.631714762477434) :
#   'chat' given is < 1, increased to 1
# ..snip..
Henrik
  • 65,555
  • 14
  • 143
  • 159
  • Henrik: what version of `MuMIn` are you using? I'm using 1.9.13. When I are run your code, the only variables I get in the model selection table are `(Intrc)` `disp` `mpg` `wt` `df`. I do not get the other 4. – luciano Mar 04 '14 at 09:50
  • @luciano, I use 1.9.13. I am glad you are asking, because at one point I actually got the 'partial' output you describe, for no obvious reason. I can't say if I did something particular to make it work again - it 'suddenly' did. Perhaps `MuMin` is a bit shaky, or clashes with other loaded packages. You may try the old classic: restart R, and run the script in a fresh session. – Henrik Mar 04 '14 at 10:13
  • Exactly same happened to me. Initially I got partial output then after several tries I suddenly got full output. Maybe a bug? – luciano Mar 04 '14 at 14:14
  • 1
    I am actually very surprised that you get log-likelihood with quasi* family that does not provide one (as I explained in my answer earlier and `?QAIC`). So the output that @luciano got is not a bug, but expected behaviour. – Kamil Bartoń Mar 04 '14 at 20:13
3

It's really all in the manual, so please read it first (?model.sel and ?QAIC). Note two issues in your code:

  • arguments to QAIC in model.sel are passed in the rank.args argument, not directly.
  • models of quasi- families do not report a likelihood, needed to calculate QAIC. See Note in ?QAIC and example(QAIC) for a hack to go around this.
Kamil Bartoń
  • 1,482
  • 9
  • 10