rm(list=(ls()))
main.directory <- getwd()
data.directory <- paste(main.directory, "LASSO.data", sep = "/")
image.directory <- paste(main.directory, "LASSO.images", sep = "/")
results.directory <- paste(main.directory, "LASSO.results", sep = "/")
get.csv.data.file <- function(pFileName) {
setwd(data.directory)
mydataframe <- read.csv(pFileName)
}
write.csv.data.file <- function(pDf, pFileName) {
setwd(data.directory)
write.csv(pDf, pFileName, row.names=FALSE)
}
write.object.to.text.file.first <- function(pObject, pFile) {
(out <- capture.output(pObject))
setwd(results.directory)
cat(out, file = pFile, sep= "\n", append = FALSE)
}
write.object.to.text.file <- function(pObject, pFile) {
(out <- capture.output(pObject))
setwd(results.directory)
cat(out, file = pFile, sep= "\n", append = TRUE)
}
# install covTest from github
# devtools::install_github("https://github.com/cran/covTest")
suppressMessages(suppressWarnings(library(tidyverse, warn.conflicts=F, quietly=T)))
suppressMessages(suppressWarnings(library(Ecdat, warn.conflicts=F, quietly=T)))
suppressMessages(suppressWarnings(library(glmnet, warn.conflicts=F, quietly=T)))
suppressMessages(suppressWarnings(library(covTest, warn.conflicts=F, quietly=T)))
# ls(package:covTest)
# ls(package:glmnet)
# ls(package:Ecdat)
CAschool <- Ecdat::Caschool
case <- c(1:dim(CAschool)[1])
CAschool <- cbind(case, CAschool)
write.csv.data.file(CAschool,"CAschool.csv")
CAschool <- get.csv.data.file("CAschool.csv")
names(CAschool)[1] <- "case"
write.csv.data.file(CAschool,"CAschool.csv")
rm(CAschool)
# Multivariate Behavioral Research 1036965, 14-089.R2
#
# Using Lasso for Predictor Selection and to Assuage Overfitting: A Method Long Overlooked in Behavioral Sciences
#
# Daniel M. McNeish
# Department of Human Development and Quantitative Methodology,
# University of Maryland, College Park
#
# Appendix
#
# Code for Example Lasso Analysis in R and SAS
#
# R Code
#
# Lasso to select predictors and estimate regression coefficients
ca <- get.csv.data.file("CAschool.csv")
ca <- dplyr::select(ca, enrltot, teachers, calwpct, mealpct, computer, compstu,
expnstu, str, avginc, elpct, everything())
#### Read in the data from CSV###
IV <- (ca[,1:10])
### Create an object of candidate predictors. The example had ten predictors so the value spans from 1 to 10. It is easiest to setup the data such that all candidate predictors are in sequential columns ###
IV1 <- scale(IV, FALSE, FALSE)
### Create a matrix like object from the candidate predictors. the "FALSE" arguments tell R not to center or standardize the variables ###
require(glmnet)
### The glmnet package must be downloaded before using the require command and the following functions ###
attach(ca)
### Attaching the data set does not require using data prefixes, this is just for convenience ###
set.seed(1028)
### set a seed value so that a results are reproducible and do not change every time the program is run##
# Lambda <- cv.glmnet(IV1,testscr);Lambda
Lambda <- cv.glmnet(IV1,testscr)
Lambda
##
## Call: cv.glmnet(x = IV1, y = testscr)
##
## Measure: Mean-Squared Error
##
## Lambda Index Measure SE Nonzero
## min 0.1901 49 72.39 5.616 7
## 1se 1.7728 25 77.76 5.308 3
write.object.to.text.file.first("", "LASSO.CA.school.results.txt")
write.object.to.text.file("***** LASSO lambda *****",
"LASSO.CA.school.results.txt")
write.object.to.text.file("", "LASSO.CA.school.results.txt")
write.object.to.text.file(Lambda, "LASSO.CA.school.results.txt")
### Estimate optimal values for the regularization parameter. This outputs the minimum and the value with one SE. testscr is the dependent variable in this example. Values will change slightly each time this code is run. ###
lambda.results <- coef(Lambda, s=Lambda$lambda.1se)
lambda.results
## 11 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 665.6706595
## enrltot .
## teachers .
## calwpct .
## mealpct -0.3913444
## computer .
## compstu .
## expnstu .
## str .
## avginc 0.5255236
## elpct -0.1311628
write.object.to.text.file("", "LASSO.CA.school.results.txt")
write.object.to.text.file("***** LASSO coefficients *****",
"LASSO.CA.school.results.txt")
write.object.to.text.file("", "LASSO.CA.school.results.txt")
write.object.to.text.file(lambda.results, "LASSO.CA.school.results.txt")
### Obtain regression coefficients for candidate predictors. s = is the value for the regularization parameter ###
RegCoef <- glmnet(IV1,testscr,family = "gaussian",alpha = 1)
### This function is used to plot the the coefficients in a later function ###
plot(Lambda)
### Plot the MSE for values of Lambda. This is what produces Figure 1 ###
### Code for Figure 2 - Run all three at once ###
plot(RegCoef, xvar="lambda", lwd = 1.8)
abline(v=log(Lambda$lambda.1se))
abline(v=log(Lambda$lambda.min))
### This plots the coefficients for different values of lambda without restricting the y-axis scale at all. The "abline" functions impose the minimum and within 1 SE values of Lambda onto the plot. ###
### Code for Figure 3 - Run all three at once ###
plot(RegCoef, xvar="lambda",ylim=c(-1.5,1.5), lwd = 1.8)
abline(v=log(Lambda$lambda.1se))
abline(v=log(Lambda$lambda.min))
### This plot shows the value of the regression coefficients for different values of lambda. The "abline" functions impose the minimum and within 1 SE values of Lambda onto the plot. The ylmin option restricts the range of the y axis to be only between -1.5 and 1.5 because most coefficients are small but one is much larger. Readers may need to change this depending on the scale of their variables ###
### Outputting p-values for Lasso coefficients
ca <- get.csv.data.file("CAschool.csv")
ca <- dplyr::select(ca, enrltot, teachers, calwpct, mealpct, computer,
compstu, expnstu, str, avginc, elpct, everything())
### Same as above, can be omitted is previously code is run first ###
IV <- (ca[,1:10])
### Same as above, can be omitted is previously code is run first ##
df <- nrow(IV)-1
### the number of observations minus 1, used in computations below ###
attach(ca)
### Same as Above ###
require (covTest)
### Must download covTest package before using subsequent commands ###
IV2=scale(IV,TRUE,TRUE)/sqrt(df)
### this command does require that the values are centered and standardized ###
LarsCoef=lars(IV2,testscr)
### the object LarsCoef houses the steps of the Least Angle Regression algorithm for the IV2 candidate predictors and testscr dependent variable ###
covTest.results <- covTest(LarsCoef,IV2,testscr)
covTest.results
## $results
## Predictor_Number Drop_in_covariance P-value
## 4 926.3781 0.0000
## 9 18.4785 0.0000
## 10 31.0445 0.0000
## 6 0.1122 0.8939
## 7 1.5759 0.2081
## 8 0.1130 0.8932
## 3 1.6948 0.1849
## 5 0.0173 0.9829
## 1 0.3509 0.7043
## 2 0.0510 0.9503
##
## $sigma
## [1] 8.4068
##
## $null.dist
## [1] "F(2,410)"
write.object.to.text.file("", "LASSO.CA.school.results.txt")
write.object.to.text.file("***** Significance of LASSO coefficients *****", "LASSO.CA.school.results.txt")
write.object.to.text.file("", "LASSO.CA.school.results.txt")
write.object.to.text.file(covTest.results, "LASSO.CA.school.results.txt")
### CovTest will output a p-value for each candidate predictor and the residual variance ###