## Seleção de variáveis baseada em outros critérios. str(longley) help(longley, help_type="html") m0 <- lm(Employed ~ ., data=longley) f0 <- formula(m0) class(f0) methods(class="formula") terms(f0) ## Variáveis independentes presentes na formula do modelo. all.vars(f0[-2]) t0 <- labels(terms(f0)) ## Pelos valores de R^2 ajustado. sapply(t0, function(term){ m1 <- lm(Employed ~ ., data=longley[,!names(longley)==term]) summary(m1)$adj.r.squared }) ## Pelos valores dos p-valores do abandono de cada termo. sapply(t0, function(term){ m1 <- lm(Employed ~ ., data=longley[,!names(longley)==term]) anova(m1, m0)[2,6] }) ## Pela soma de quadrados dos erros de previsão. sapply(t0, function(term){ m1 <- lm(Employed ~ ., data=longley[,!names(longley)==term]) sum((residuals(m1)/(1-hatvalues(m0)))^2) }) ## Pela soma de quadrados dos erros de previsão. sapply(t0, function(term){ m1 <- lm(Employed ~ ., data=longley[,!names(longley)==term]) n <- length(residuals) p1 <- length(coef(m1)) p0 <- length(coef(m0)) SQR1 <- deviance(m1) QMR0 <- deviance(m0)/(n-p0) cp <- SQR1/QMR0-n+2*p1 cp }) length(coef(m0)) ##----------------------------------------------------------------------------- da <- read.table("http://www.leg.ufpr.br/~walmes/data/MontgomeryASPE5th/Example12.14.txt", header=TRUE, sep="\t") names(da) <- tolower(names(da)) str(da) m0 <- lm(quality~., data=da) ## step(m0) f0 <- formula(m0) t0 <- labels(terms(f0)) sapply(t0, function(term){ m1 <- lm(quality ~ ., data=da[,!names(da)==term]) n <- length(residuals) p1 <- length(coef(m1)) p0 <- length(coef(m0)) SQR1 <- deviance(m1) QMR0 <- deviance(m0)/(n-p0) cp <- SQR1/QMR0-n+2*p1 cp }) length(coef(m0)) ## Pelos valores de R^2 ajustado. sapply(t0, function(term){ m1 <- lm(quality ~ ., data=da[,!names(da)==term]) summary(m1)$adj.r.squared }) ## Pela soma de quadrados dos erros de previsão. sapply(t0, function(term){ m1 <- lm(quality ~ ., data=da[,!names(da)==term]) sum((residuals(m1)/(1-hatvalues(m0)))^2) }) ##----------------------------------------------------------------------------- da <- read.table("http://www.leg.ufpr.br/~walmes/data/business_economics_dataset/EXAMPLES/EXECSAL.DAT") head(da) ## V1 : int, line ## V2 : num, salary ## V3 : int, experience (years) ## V4 : int, education (years) ## V5 : int, bonus eligibility (0=no, 1=yes) ## V6 : int, number of employees supervised ## V7 : int, corporate assets ## V8 : int, board member (0=no, 1=yes) ## V9 : int, age ## V10: int, campany profits ## V11: int, international resposability (0=no, 1=yes) ## V12: int, company's total sales