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]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
