CODEWITHSUNDEEP
rlm
Bekah
Bekah

Reputation: 11

AICcmodavg lm wont run additive model with 2 variables

Context Ideally I would like to run one model list taking into account all of the models at once. As I have been unable to get this to work, I have run three sets of models dealing with 1, 2 and 3 variables respectively.

Problem Although the 1 and 3 variable models run fine, the linear regression wont run an additive 2 variable model but will run an interaction term (which I don't want). The linear regression will run outside of the package fine so it is not a degrees of freedom or negative values issue but I am stumped. Does anyone have any ideas?

Data

lre<-c(0.398,0,0.9298,1.470,0)
imm1<-c(-0.54,-1.67,0.07.96,-0.862,1.02)
imm2<-c(-0.033,4.3798,0.0358,-1.045,0.592)
met1<-c(-1.689,-1.06,1.156,-1.574,1.632)
met2<-c(-1.980,1.349,1.538,0.6303,-0.310)
phy1<-c(0.202,0.368,-0.643,2.274259,0.847)
phy2<-c(1.079,-0.068,-1.438,-0.716,0.846)

1 variable = works

library(AICcmodavg) 
Cand.models <- list( )
Cand.models[[1]]<-lm(lre~imm1)
Cand.models[[2]]<-lm(lre~imm2)
Cand.models[[3]]<-lm(lre~met1)
Cand.models[[4]]<-lm(lre~met2)
Cand.models[[5]]<-lm(lre~phy1)
Cand.models[[6]]<-lm(lre~phy2)
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)

2 variables in additive model = doesn't work

Cand.models <- list( )
Cand.models[[1]]<-lm(lre~imm1+imm2)
Cand.models[[2]]<-lm(lre~imm1+met2)
Cand.models[[3]]<-lm(lre~imm1+phy2)
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)

Warning message: In aictab.AIClm(cand.set = Cand.models, modnames = Modnames, sort = TRUE) : Check model structure carefully as some models may be redundant

However...

2 variables with interaction = works

Cand.models <- list( )
Cand.models[[1]]<-lm(lre~imm1*imm2)
Cand.models[[2]]<-lm(lre~imm1*met2)
Cand.models[[3]]<-lm(lre~imm1*phy2)
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)

3 variables in additive model = works

Cand.models <- list( )
Cand.models[[1]]<-lm(lre~imm1+imm2+met1)
Cand.models[[2]]<-lm(lre~imm1+imm2+met2)
Cand.models[[3]]<-lm(lre~imm1+imm2+phy1)
Cand.models[[4]]<-lm(lre~imm1+imm2+phy2)
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)

Upvotes: 1

Views: 998

Answers (1)

alexforrence
alexforrence

Reputation: 2744

The warning is a bit misleading. It's reporting the models might be the same because they all have the same AICc, which is Inf in this case. It is a degrees of freedom issue, in the sense that the degrees of freedom matches the number of observations minus one. So when AICcmodavg:::AICc.lm uses the formula:

AICc <- -2 * LL + 2 * K * (n/(n - K - 1))

Where n = 5, the number of samples, and K = 4, the number of parameters (global intercept, imm1 and imm2) plus one for the error variance, the last term is (5/0) (Inf). Luckily, in your case, all of the potential models have the same K, meaning that the relative differences in the AIC will match those as if the AICc worked (Wikipedia):

lre <- c(0.398, 0, 0.9298, 1.470, 0)
imm1 <- c(-0.54, -1.67, 0.0796, -0.862, 1.02) # fixed typo
imm2 <- c(-0.033, 4.3798, 0.0358, -1.045, 0.592)
met1 <- c(-1.689, -1.06, 1.156, -1.574, 1.632)
met2 <- c(-1.980, 1.349, 1.538, 0.6303, -0.310)
phy1 <- c(0.202, 0.368, -0.643, 2.274259, 0.847)
phy2 <- c(1.079, -0.068, -1.438, -0.716, 0.846)

# Put everything in a data frame for cleanliness
dat <- data.frame(lre, imm1, imm2, met1, met2, phy1, phy2)

Cand.models <- list()
Cand.models[[1]]<-lm(lre ~ imm1 * imm2, data = dat)
Cand.models[[2]]<-lm(lre ~ imm1 * met2, data = dat)
Cand.models[[3]]<-lm(lre ~ imm1 * phy2, data = dat)
Modnames <- paste("mod", 1:length(Cand.models), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE)

Model selection based on AICc :

      K   AICc Delta_AICc AICcWt Cum.Wt    LL
mod 1 5 -55.53       0.00   0.99   0.99  2.77
mod 2 5 -45.60       9.94   0.01   1.00 -2.20
mod 3 5 -44.42      11.12   0.00   1.00 -2.79

Now with second.ord = FALSE (just return AIC):

aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE, second.ord = FALSE)

Model selection based on AIC :

      K   AIC Delta_AIC AICWt Cum.Wt    LL
mod 1 5  4.47      0.00  0.99   0.99  2.77
mod 2 5 14.40      9.94  0.01   1.00 -2.20
mod 3 5 15.58     11.12  0.00   1.00 -2.79

So as long as all of your models have the same number of coefficients, you can use the regular AIC to get the same (effective) result.

Upvotes: 1

Related Questions