2

I'm trying to do AICc model selection and model averaging with tweedie (compound Poisson) distributed data in R.

I was working with the AICcmodavg R package with no success, then decided to try out the MuMIn package when I came across the suggestion here (https://stats.stackexchange.com/questions/141806/glm-model-selection-using-aicc-with-tweedie-distribution) that

"You can use AICtweedie directly in MuMIn's functions, just specify it as a rank argument."

I set up my models as follows My response variable (NVIR) is catch-per-unit-effort of eastern newt adults and my explanatory variables are various habitat characteristics of my sampling sites. 

m1<- glm(NVIR~Water_T+cond+DO+ORP+pH+max_depth+type, 
    family = tweedie(link.power=0, var.power=1.3), data = cpue)
m2<- glm(NVIR~Water_T+cond+DO+ORP+pH+littoral_slope+type, 
    family = tweedie(link.power=0, var.power=1.3), data = cpue)
m3<- glm(NVIR~pH+DO+cond+max_depth+type, 
    family = tweedie(link.power=0, var.power=1.3), data = cpue)
m4<- glm(NVIR~pH+DO+cond+littoral_slope+type, 
    family = tweedie(link.power=0, var.power=1.3), data = cpue)
m5<- glm(NVIR~cond+type+pH+max_depth, 
    family = tweedie(link.power=0, var.power=1.3), data = cpue)

and then tried this line

model.sel(m1, m2, m3, m4, m5, rank = AICc, rank.args = AICtweedie)

and received the error

Error in UseMethod("logLik") : 
no applicable method for 'logLik' applied to an object of class "function"
In addition: Warning message:
In model.sel.default(m1, m2, m3, m4, m5, rank = AICc, rank.args =   AICtweedie) :
models are not all fitted to the same data

Alternatively, I also tried this line

model.sel(m1,m2,m3,m4,m5, rank.args=AICtweedie)

and got this error:

Error in get(x) : object 'Tweedie' not found
In addition: Warning message:
In model.sel.default(m1, m2, m3, m4, m5, rank.args = AICtweedie) :
models are not all fitted to the same data

I'm wondering whether the problem is with my code, or if the tweedie family is incompatible with this package.

Thank you for your time.

www
  • 38,575
  • 12
  • 48
  • 84
RachelF
  • 21
  • 2

1 Answers1

0

rank = tweedie::AICtweedie

Just for anyone who comes across this

Darren Norris
  • 1
  • do you have any idea why `rank=AICtweedie` didn't work for the OP (as reported in comments)? Is there an `AICtweedie` object in more than one package? – Ben Bolker Jun 15 '18 at 22:30
  • 1
    @BenBolker `model.sel` assumes `family(model)$family` would be a name of the family function. In this case it is `"Tweedie"` not `"tweedie"` and hence the problem. There are also other problems with handling `tweedie` family - I'll look into it soon. – Kamil Bartoń Jun 16 '18 at 14:36
  • hmm. In `glmmTMB`, `beta_family()` gives rise to an object with `$family` equal to "beta". Is there going to be a way to deal with this more generally ... ? (We might be able to change the `$family` component without breaking anything ... let me know offline if that would be the best way ...) – Ben Bolker Jun 16 '18 at 14:44
  • 1
    @BenBolker Don't worry about it - I see that the code should not expect all `family` objects to conform to the standard structure, so I guess I just remove that bit (which is not crucial in any way). Even the `MASS::negative.binomial` is a pain with it with its `$family` such as `"Negative Binomial(1.23)"`. – Kamil Bartoń Jun 18 '18 at 11:31
  • Thanks all! I appreciate your time and comments. – RachelF Jun 26 '18 at 17:50