Hello,
I have stumbled across a peculiar problem. I am writing code for a
simulation, and would like to package it into a function that completes one
iteration. I'm using covTest sig. tests for LASSO (link is below to info
and the source code). When this function is used within a function (F.1
below) it only works when n=15, and the results do not make sense (Part 1
below). When used outside of the function, it works fine with any n and
the results makes sense (Part 2 below). Can someone please look at this
and see if you can make sense of what is going on. I am dumbfounded as to
why this is acting weird inside the function.
Please note, the package "covTest" link found here should be in your
WD:
http://andrewgelman.com/2013/03/18/tibshirani-announces-new-research-result-a-significance-test-for-the-lasso/
Thank you.
BS
rm(list=ls()) # note: remove all
# Part 1:
F.1 <- function(n){
require("lars")
require("mvtnorm")
source("covTest/R/funcs.R")
# Simulated Data
##############################################################
X <- mvrnorm(n, mu=c(0, 0, 0), Sigma=diag(3))
b1 =2; b2 = 0.0; b3 = - 2
Y <- as.vector(X[,1]*b1 +X[,2]*b2 +X[,3]*b3 + rnorm(n, 0, 0.05) )
fit2 <- lars(as.matrix(X), Y, type="lasso", max.steps=100)
res <- covTest(fitobj=fit2, x=as.matrix(X), y=as.vector(Y)) # package
folder must be in WD
res <- res$results[order(res$results[,1]),3]
lasso <- c(fit2$beta[names(fit2$Cp[which.min(fit2$Cp)]),], res)
names(lasso)[4:6] <- c("X1.pval", "X2.pval",
"X3.pval")
return(lasso)
# for output
#########################################################################
}
# function end
# now use F.1 in a loop...
it = 3
out <- vector("list")
for(i in 1:it){
out[[i]] <- F.1(15) # this will only work for n=15, and p-values cannot be
close to correct
}
######################################################################################
out
rm(list=ls()) # note: remove all
# Part 2:
######################################################################################
require("lars")
require("mvtnorm")
source("covTest/R/funcs.R")
it = 3
n= 30 # here n can be anything, and the results make sense
out <- vector("list")
for(i in 1:it){
X <- mvrnorm(n, mu=c(0, 0, 0), Sigma=diag(3))
b1 =2; b2 = 0.0; b3 = - 2
Y <- as.vector(X[,1]*b1 +X[,2]*b2 +X[,3]*b3 + rnorm(n, 0, 0.05) )
fit2 <- lars(as.matrix(X), Y, type="lasso", max.steps=100)
res <- covTest(fitobj=fit2, x=as.matrix(X), y=as.vector(Y)) # package
folder must be in WD
res <- res$results[order(res$results[,1]),3]
lasso <- c(fit2$beta[names(fit2$Cp[which.min(fit2$Cp)]),], res)
names(lasso)[4:6] <- c("X1.pval", "X2.pval",
"X3.pval")
# for output
#########################################################################
out[[i]] <- lasso
}
######################################################################################
out
[[alternative HTML version deleted]]