CODEWITHSUNDEEP
rmumin
agorapotatoes
agorapotatoes

Reputation: 351

dredge in MuMin (R) keeps models with higher order terms without their respective lower order terms

I am using the function dredge from package MuMin for model selection. My model has a polynomial in it. I was under the impression (https://www.rdocumentation.org/packages/MuMIn/versions/1.42.1/topics/dredge, see "Interactions") that MuMin as standard ditches models that include the higher order term without its respective lower order term; however when I use it the top model includes only the higher order term.

Here is the code I am using to make the global model, which includes hour: at both first and second order

sc.mod.env.hour<-clmm(seedcat~hour+I(hour^2)+MeanSpeed+RH+TKE+AirTemp+MeanUpdraft+(1|jdayfac), data=sc.data, na.action = na.fail, Hess =TRUE)

However when I run dredge:

dro.env.hour<-dredge(sc.mod.env.hour)

The output includes models that include I(hour^2) without hour

Model selection table 
    (Intrc)   ArTmp    hour hour^2 MnSpd    MnUpd     RH      TKE df   logLik   AICc  delta weight
46        + -1.2170         -3.340 1.894          -4.219          13 -571.772 1170.5   0.00  0.154
62        + -1.3410         -3.478 1.964  0.30990 -4.269          14 -570.882 1170.9   0.37  0.128
45        +                 -2.882 2.052          -3.214          12 -573.245 1171.3   0.81  0.103
47        +         -0.3595 -2.670 1.942          -3.601          13 -572.374 1171.7   1.20  0.084
63        +         -0.4495 -2.712 2.010  0.33870 -3.645          14 -571.350 1171.8   1.30  0.080

I have tried also using poly(hour,2) to define the global model but this results in only including a single term for hour.

I am using clmm but have tried a simpler model using lm and got the same result.

Any guidance appreciated thankyou. I wasn't sure if this should be on crossvalidated but it is a programming issue rather than a stats issue so thought it should be here.

EDIT: I have solved the issue by using subset:

dro.env.hour<-dredge(sc.mod.env.hour, subset=(dc(hour, I(hour^2)))

Although it is not clear why this is an issue in the first place.

Upvotes: 4

Views: 2497

Answers (2)

Ben Bolker
Ben Bolker

Reputation: 226861

Another possible solution (untested!) is to use poly(hour, 2) rather than hour + I(hour^2) in the model (or, if you really want 'raw' rather than orthogonal polynomials, poly(hour, 2, raw=TRUE)). This way R can recognize that poly() is meant to be a single term. (The problem with hour + I(hour^2) is that it would take some fancy parsing within the model-frame construction machinery to be able to tell that these terms are related ...)

Upvotes: 0

Mart&#237;n A. Colombo
Mart&#237;n A. Colombo

Reputation: 11

In case someone else has this question, I found this solution in MuMin manual:

Make a subset with a name first, in this case I want quadratic effects of X and Y)

msubset <- expression(  dc(X, `I(X^2)`) &
                        dc(`Y`, `I(Y^2)`))

Then the model, for example

mglm <- glm(response~X + I(X^2) + Y + I(Y^2)+ Z,
                  family=binomial,
                  na.action = "na.fail",
                  data=mdata)

Now dredge using the subset you specified

mmodels<-dredge(mglm, subset=msubset)

Check summary table

model.sel(mmodels)

Upvotes: 1

Related Questions