2016-01-06 40 views
1

上下文 理想情况下,我想运行一个模型列表,同时考虑所有的模型一次。由于我一直无法做到这一点,所以我运行了三组分别处理1,2和3个变量的模型。AICcmodavg lm不会运行2变量添加剂模型

问题 虽然图1个3变量模型运行良好,线性回归不会运行添加剂的2变量模型,但将运行一个相互作用项(我不想)。线性回归将在包裹之外运行,所以它不是一个自由度或负值问题,但我很难过。有没有人有任何想法?

数据

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可变=工作

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个变量相加模型=不起作用

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) 

警告消息: 在aictab.AIClm(cand.set = Cand.models,modnames = Modnames,排序= TRUE):仔细检查 模型结构的某些型号可能是多余的

但是......

2变量与相互作用=工作

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个变量中相加模型=工作

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) 

回答

0

该警告有点令人误解。这是报告模型可能是相同的,因为它们都有相同的AICc,在这种情况下是Inf。它是自由度问题,意思是自由度与观测值的数量相符减1。所以,当AICcmodavg:::AICc.lm使用公式:

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

n = 5,样本的数量,和K = 4,参数之一的误差方差的数量(全球拦截,imm1imm2)加时,最后一项是(5/0)(天道酬勤)。幸运的是,在你的情况下,所有潜在的车型具有相同的K,这意味着在AIC的相对差异将匹配那些仿佛AICc工作(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 
second.ord = FALSE

现在(只是返回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 

所以只要你所有的机型都具有相同数量的系数,你可以使用常规的AIC获得相同的(有效)的结果。