Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

Since regression modeling is often more "art" than science, I often find myself testing many iterations of a regression structure. What are some efficient ways to summarize the information from these multiple model runs in an attempt to find the "best" model? One approach I've used is to put all the models into a list and run summary() across that list, but I imagine there are more efficient ways to compare?

Sample code & models:

ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2,10,20, labels=c("Ctl","Trt"))
weight <- c(ctl, trt)

lm1 <- lm(weight ~ group)
lm2 <- lm(weight ~ group - 1)
lm3 <- lm(log(weight) ~ group - 1)

#Draw comparisions between models 1 - 3?

models <- list(lm1, lm2, lm3)

lapply(models, summary)
share|improve this question
4  
Sounds a bit like data dredging to me. Shouldn't the focus be on what you plausibly think is an appropriate model, what covariates, transformations etc before you start modelling. R doesn't know you did all that model fitting to find a good model. –  Gavin Simpson Feb 3 '11 at 19:08
3  
@Gavin - I can see this getting horrendously off-topic very quickly, but the short answer is no, I am not advocating data dredging or finding spurious relationships between random variables in a dataset. Consider a regression model that includes income. Is it not reasonable to test transformations on income to see their impact on the model? Log of income, log of income in 10s of dollars, log of income in 100s...?Even if this is data dredging - a function / summary tool that can aggregate the output from many model runs would still be very helpful, no? –  Chase Feb 3 '11 at 19:32

2 Answers 2

up vote 17 down vote accepted

Plot them!

enter image description here

Or, if you must, use tables: The apsrtable package or the mtable function in the memisc package.

Using mtable

 mtable123 <- mtable("Model 1"=lm1,"Model 2"=lm2,"Model 3"=lm3,
     summary.stats=c("sigma","R-squared","F","p","N"))

> mtable123

Calls:
Model 1: lm(formula = weight ~ group)
Model 2: lm(formula = weight ~ group - 1)
Model 3: lm(formula = log(weight) ~ group - 1)

=============================================
                 Model 1   Model 2   Model 3 
---------------------------------------------
(Intercept)      5.032***                    
                (0.220)                      
group: Trt/Ctl  -0.371                       
                (0.311)                      
group: Ctl                 5.032***  1.610***
                          (0.220)   (0.045)  
group: Trt                 4.661***  1.527***
                          (0.220)   (0.045)  
---------------------------------------------
sigma             0.696      0.696     0.143 
R-squared         0.073      0.982     0.993 
F                 1.419    485.051  1200.388 
p                 0.249      0.000     0.000 
N                20         20        20     
=============================================

share|improve this answer
1  
@Eduardo, +1, nice graph. It should be used with care though when the different transformation of dependent variable is used in different regressions. –  mpiktas Feb 3 '11 at 19:39
    
mpiktas, that´s true in a table as well. Graphs just make it more compact, at the expense of precision. –  Eduardo Leoni Feb 3 '11 at 19:44
    
@Eduardo can you please share the code for graphs? –  suncoolsu Feb 4 '11 at 1:40
2  
@suncoolsu R code is available on the first link given in @Eduardo's response. He he, it's grid, not lattice :) –  chl Feb 4 '11 at 7:58
    
@Eduardo - Thanks for the detailed answer, I wasn't aware of memisc previously, looks like a very handy package to have in one's quiver! –  Chase Feb 4 '11 at 16:36

The following doesn't answer exactly the question. It may give you some ideas, though. It's something I recently did in order to assess the fit of several regression models using one to four independent variables (the dependent variable was in the first column of the df1 dataframe).

# create the combinations of the 4 independent variables
library(foreach)
xcomb <- foreach(i=1:4, .combine=c) %do% {combn(names(df1)[-1], i, simplify=FALSE) }

# create formulas
formlist <- lapply(xcomb, function(l) formula(paste(names(df1)[1], paste(l, collapse="+"), sep="~")))

The contents of as.character(formlist) was

 [1] "price ~ sqft"                     "price ~ age"                     
 [3] "price ~ feats"                    "price ~ tax"                     
 [5] "price ~ sqft + age"               "price ~ sqft + feats"            
 [7] "price ~ sqft + tax"               "price ~ age + feats"             
 [9] "price ~ age + tax"                "price ~ feats + tax"             
[11] "price ~ sqft + age + feats"       "price ~ sqft + age + tax"        
[13] "price ~ sqft + feats + tax"       "price ~ age + feats + tax"       
[15] "price ~ sqft + age + feats + tax"

Then I collected some useful indices

# R squared
models.r.sq <- sapply(formlist, function(i) summary(lm(i))$r.squared)
# adjusted R squared
models.adj.r.sq <- sapply(formlist, function(i) summary(lm(i))$adj.r.squared)
# MSEp
models.MSEp <- sapply(formlist, function(i) anova(lm(i))['Mean Sq']['Residuals',])

# Full model MSE
MSE <- anova(lm(formlist[[length(formlist)]]))['Mean Sq']['Residuals',]

# Mallow's Cp
models.Cp <- sapply(formlist, function(i) {
SSEp <- anova(lm(i))['Sum Sq']['Residuals',]
mod.mat <- model.matrix(lm(i))
n <- dim(mod.mat)[1]
p <- dim(mod.mat)[2]
c(p,SSEp / MSE - (n - 2*p))
})

df.model.eval <- data.frame(model=as.character(formlist), p=models.Cp[1,],
r.sq=models.r.sq, adj.r.sq=models.adj.r.sq, MSEp=models.MSEp, Cp=models.Cp[2,])

The final dataframe was

                      model p       r.sq   adj.r.sq      MSEp         Cp
1                price~sqft 2 0.71390776 0.71139818  42044.46  49.260620
2                 price~age 2 0.02847477 0.01352823 162541.84 292.462049
3               price~feats 2 0.17858447 0.17137907 120716.21 351.004441
4                 price~tax 2 0.76641940 0.76417343  35035.94  20.591913
5            price~sqft+age 3 0.80348960 0.79734865  33391.05  10.899307
6          price~sqft+feats 3 0.72245824 0.71754599  41148.82  46.441002
7            price~sqft+tax 3 0.79837622 0.79446120  30536.19   5.819766
8           price~age+feats 3 0.16146638 0.13526220 142483.62 245.803026
9             price~age+tax 3 0.77886989 0.77173666  37884.71  20.026075
10          price~feats+tax 3 0.76941242 0.76493500  34922.80  21.021060
11     price~sqft+age+feats 4 0.80454221 0.79523470  33739.36  12.514175
12       price~sqft+age+tax 4 0.82977846 0.82140691  29640.97   3.832692
13     price~sqft+feats+tax 4 0.80068220 0.79481991  30482.90   6.609502
14      price~age+feats+tax 4 0.79186713 0.78163109  36242.54  17.381201
15 price~sqft+age+feats+tax 5 0.83210849 0.82091573  29722.50   5.000000

Finally, a Cp plot (using library wle)

share|improve this answer
2  
Thanks for posting this, very insightful on how to approach issues like this. +1! –  Chase Feb 4 '11 at 16:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.