# SELEÇÃO DE MODELOS - EXEMPLO COM DADOS DE SOLO

mult=read.table("A7_RLM.TXT", h=T);mult
X1 <- mult[,1]
X2 <- mult[,2]
X3 <- mult[,3]
X4 <- mult[,4]
X5 <- mult[,5]
X6 <- mult[,6]
Y <- mult[,7]
cor(mult)
#Entrando com o Modelo Nulo (Y = B0) e com o Modelo Completo (Y = B0+B1X1+...+B6X6)
MN = lm(Y ~ 1, data=mult)
MC = lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6, data=mult)

# Método: Passo atrás - "Backward", utilizando o critério AIC (k=2)
# Para o BIC usar a opção k=log(nrow(data)) no lugar de k=2
step(MC, data = mult, direction="backward",k=2, trace=1)
step(MC, data = mult, direction="backward",k=log(nrow(mult)), trace=1)

# Método: Passo atrás - "Forward", utilizando o critério AIC (k=2)
# Para o BIC usar a opção k=log(nrow(data)) no lugar de k=2
step(MN, scope=list(lower=MN,upper=MC), data = mult, direction="forward", k=2, trace=1)
step(MN, scope=list(lower=MN,upper=MC), data = mult, direction="forward", k=log(nrow(mult)), trace=1)

# Método: Passo a passo - "Stepwise", utilizando o critério AIC (k=2)
# Para o BIC usar a opção k=log(nrow(data)) no lugar de k=2
step(MC, direction="both", k=2, trace=1)

step(MN, scope=list(lower=MN,upper=MC), data = mult, direction="both", k=2, trace=1)

# Selecao de Modelo - método Cp (Menor Cp - melhor o modelo) #
require(leaps)
selcp = leaps(x=cbind(X1,X2,X3,X4,X5,X6), y=Y,method=c("Cp"))
nparcp = selcp$size
Cp = selcp$Cp
cbind(nparcp,Cp,selcp$which)
plot(nparcp,Cp,xlab="No de parametros (p)", ylab="Cp")
abline(0,1)

# Selecao de Modelo - método R2 Ajustado (Maior R2 - melhor o modelo#
#require(leaps), também precisa do pacote mas já foi instalado
selr2 = leaps(x=cbind(X1,X2,X3,X4,X5,X6), y=Y,method=c("adjr2"))
nparr2 = selr2$size
r2.aj = selr2$adjr2
cbind(nparr2,r2.aj,selr2$which)

#############################################################
#### COMPLETE A ANÁLISE COM O MODELO SELECIONADO
mod.sel <- lm(Y~X2+X5+X6)
anova(mod.sel)
summary(mod.sel)
###########################################################