set.seed(123) # Allow reproducible random numbers ## Create some data with 5 different categories & Gaussian random noise: mydata <- within(data.frame(x = 1:50), { class <- rep(1:5, each = 10) y <- 0.5 * class * (1:10) + rnorm(50) class <- factor(class) # convert to a factor } ) svg(file="latticePlot.svg",height=5.5,width=6.2) library(lattice) xyplot(y ~ x, data = mydata, group = class, auto.key = list(title = "Class", columns = 5) ) dev.off() ## create a vector of model formulae: forms <- paste("y ~ ", c("x", # same intercept & slope "x + class", # different intercept, same slope "x / class", # same intercept, different slope "x * class", # different intercept & slope "class / x", # different intercept & slope "class / (1 + x) - 1" # different intercept & slope ), sep = "" ) names(forms) <- paste("Model", LETTERS[1:length(forms)]) forms # print out the object ## fit each model to the data: ##install.packages("plyr") # if package not already installed library(plyr) models <- llply(forms, lm, data = mydata) names(models) <- names(forms) # assign model name to model object in list ## Use an anonymous function to add predicted data + AIC fit criterion: predData <- ldply(models, function(m) within(m$model, # combine with fitted x & y data { pred <- predict(m) AIC <- AIC(m) } ) ) head(predData) # print out first few rows of data frame ## add formula column to data frame & create custom panel strip title: predData <- within(predData, { modelName <- .id # preserve model name formula <- rep(as.character(forms), each = nrow(predData) / length(forms)) strip <- factor(paste(formula, ", AIC =", round(AIC, 1))) strip <- reorder(strip, AIC) # sort by AIC value } ) svg(file="ggplot2_data-prediction.svg",height=5.5,width=8.0) ## plot points & predictions for each model: ##install.packages("ggplot2") # if package not already installed library(ggplot2) ggplot(data = predData, aes(x, y, colour = class)) + geom_point(alpha = 0.5) + # use 50% transparent points, to handle overplotting geom_line(aes(y = pred)) + # show model prediction as a line facet_wrap( ~ strip) # use separate panel for each different model dev.off()