toxic <- read.table("http://www.stat.ufl.edu/~winner/data/air_int_incap.dat", header=F,col.names=c("matcat","matid","timeinc","Y","CO","HCN","H2S", "HCL","HBr","NO2","SO2")) attach(toxic) library(car) toxic.1 <- lm(Y ~ CO + HCN + H2S + HCL + HBr + NO2 + SO2) summary(toxic.1) anova(toxic.1) drop1(toxic.1,test="F") vif(toxic.1) ######### Perform Backward Elimination, Forward Selection, and Stepwise Regression ######### Based on Model AIC (not individual regression coefficients) ######### fit1 and fit2 represent "extreme" models library(MASS) fit1 <- lm(Y ~ CO + HCN + H2S + HCL + HBr + NO2 + SO2) fit2 <- lm(Y ~ 1) stepAIC(fit1,direction="backward") stepAIC(fit2,direction="forward",scope=list(upper=fit1,lower=fit2)) stepAIC(fit2,direction="both",scope=list(upper=fit1,lower=fit2)) ########## Perform all possible regressions (aka all subset regressions) ########## Prints out best 4 models of each # of predictors install.packages("leaps") library(leaps) alltoxic <- regsubsets(Y ~ CO + HCN + H2S + HCL + HBr + NO2 + SO2, nbest=4,data=toxic) aprout <- summary(alltoxic) n <- length(toxic$Y) pprime <- apply(aprout$which, 1, sum) aprout$aic <- aprout$bic - log(n) * pprime + 2 * pprime with(aprout,round(cbind(which,rsq,adjr2,cp,bic,aic),3)) ## Prints "readable" results plot(alltoxic,scale="bic") plot(alltoxic,scale="adjr2")