I am trying to replicate a simulation from an article for a class project. I have created simulated matrix x and a simulated vector y. I want to first reduce the model size from p=1000 to p=n/log(n) using SIS. I want to then further reduce the vector space to p' by using the Primal Dual Dantzig selector.
There is no direct predict function in this package so I am plugging the p' selected predictors and their coefficients into a lm function to get predictions. I am using the predict.lm function to get the predicted values. However, the returned vector is all NA.
The following code may illustrate the problem further. I am not getting any errors after running the following code. I have to follow this procedure of using SIS and then Dantzig Selector to match the simulation in the article.
library(SIS) #import the SIS library
library(fastclime)
errors_DS_small = c() #create empty array to store errors from each simulation
model_sizes_DS_small = c() #creat empty array to store selected model size from each simulation
start_time = Sys.time()
set.seed(123*1) #re-create results at a later time but different seed for each data set
n = 200 #sample size
p = 1000 #variables
x = matrix(rnorm(n*p, mean=0, sd=1), n, p) #creating IID standard Gaussian random predictors
# gaussian response
set.seed(456*1) #re-create results at a later time but different seed for each data set
s = 8
u = rbinom(s, 1, 0.4)
z = rnorm(s, mean=0, sd=1)
a = 4*log(n)/sqrt(n)
b= ((-1)**(u))*(a + abs(z))
y=x[, 1:s]%*%b
#creating SIS-DS model and gaussian response. iter=FALSE means not doing ISIS
modelSIS_small = SIS(x, y, family='gaussian', iter = FALSE, nsis=(n/log(n)))
modelDS = dantzig(x[,modelSIS_small$sis.ix0], y)
selectDS = dantzig.selector(modelDS$lambdalist, modelDS$BETA0, lambda=min(modelDS$lambdalist))
linearMod = lm(y~x[,modelSIS_small$sis.ix0])
linearMod$coefficients = selectDS
newx = x[,modelSIS_small$sis.ix0]
predTest = predict(linearMod, data=newx) #create predictions using test data test
a = modelDS$BETA0[, 9] != 0
mse = mean((y - predTest)^2) #compare predictions to real values of test y
rmse = sqrt(mse) #Square root of MSE (Square-loss function)
errors_DS_small[1] = rmse #store the RMSE
model_sizes_DS_small[1] = table(a)["TRUE"] #Store model size including intercept
print(c(1, errors_DS_small[1], model_sizes_DS_small[1])) #print results
end_time = Sys.time()
DS_total_time_small = end_time - start_time
DS_error_small_median = median(errors_DS_small)
DS_model_Size_small_median = median(model_sizes_DS_small)
print(c("DS", DS_model_Size_small_median, DS_error_small_median, DS_total_time_small))