Variable selection
In each simulation, a sample is drawn randomly. Therefore, we need a seed for the random number generator:
seed <- 340
set.seed(seed)Given the sample, we use SurSearch_s() function to find
the best model. This function will select one best mode for each of the
evaluation methods. We use the results for AIC, SIC, RMSE, and CRPS
measures:
measureOptions <- GetMeasureOptions(
typesIn = c("aic", "sic"),
typesOut = c("rmse", "crps")
)
measureNames <- c("AIC", "SIC", "RMSE", "CRPS")If you are not familiar with the arguments of the
SurSearch_s() function, please read this vignette first. The
measureNames vector is just a display name that we will use
later.
In order to compare the performance of different measures, we use the following scoring schema:
GetScore <- function(indices1, indices2) {
if (length(indices1) == 2 && length(indices2) == 2 &&
all(indices1 == c(6, 12)) && all(indices2 == c(12, 18))) {
return(1)
}
score <- 0
if (any(indices1 == 6)) score <- score + 0.1
if (any(indices1 == 12)) score <- score + 0.1
if (any(indices2 == 12)) score <- score + 0.1
if (any(indices2 == 18)) score <- score + 0.1
return(score)
}This schema uses the indices of the regressors in the (restricted)
best model and calculates the score. If the indices belong to the true
model (i.e., indices1=c(6,12) and
indices2=c(12,18)), the score is 1. otherwise each correct
index (i.e., 6 and 12 for the first equation and 12 and 18 for the
second equation) adds 0.1 points to the score. For example, this means
the maximum score for an incorrect model is 0.4. Such a model has all
the true explanatory variables and some irrelevant ones. The minimum
score is zero. This is the score of a model with no correct explanatory
variable.
The SurSearch_s() function also reports the inclusion
weights. The previous function is used to get the score of such a list.
In this case, indices1 and indices2 are the
index of the first 2 explanatory variables in the inclusion list.
Before we continue, we need a helper function to get the required indices from the estimations:
helper_GetIndices <- function(measureRes, measureSumRes) {
isRest1 <- measureSumRes$target1$model$bests$best1$coefs$isRestricted
isRest2 <- measureSumRes$target2$model$bests$best1$coefs$isRestricted
if (is.null(isRest1)) {
mindices1 <- sort(measureRes$target1$model$bests$best1$exoIndices)
} else {
mindices1 <- sort(measureRes$target1$model$bests$best1$exoIndices[isRest1[, 1] != 1])
}
if (is.null(isRest2)) {
mindices2 <- sort(measureRes$target2$model$bests$best1$exoIndices)
} else {
mindices2 <- sort(measureRes$target2$model$bests$best1$exoIndices[isRest2[, 2] != 1])
}
end <- nrow(measureRes$target1$model$inclusion)
iindices1 <- sort(measureRes$target1$model$inclusion[3:end, 1, drop = F],
index.return = TRUE, decreasing = TRUE
)
iindices2 <- sort(measureRes$target2$model$inclusion[3:end, 1, drop = F],
index.return = TRUE, decreasing = TRUE
)
return(list(
model1 = sort(mindices1), model2 = sort(mindices2),
inclusion1 = sort(iindices1$ix[1:2]), inclusion2 = sort(iindices2$ix[1:2])
))
}The central piece of this experiment is the following function in which we generate a sample, search for the best models for each performance measure, and calculate the score of the best models:
SampleAndEstimate <- function(i, res, n, s, xSizes, xCounts, yGroups,
measureOptions, searchItems, checkItems) {
sample <- GetSample(n, s)
search_res <- SurSearch_s(
y = sample$y, x = sample$z, numTargets = 2,
xSizes = xSizes, counts = xCounts,
yGroups = yGroups,
searchSigMaxIter = 1,
searchSigMaxProb = 0.1,
searchItems = searchItems,
measureOptions = measureOptions,
modelCheckItems = checkItems,
searchOptions = GetSearchOptions(printMsg = FALSE),
savePre = NULL
)
result_sum <- summary(search_res, y = sample$y, x = sample$z, test = FALSE, printMsg = FALSE)
if (i < 0) {
return(list(search_res = search_res, result_sum = result_sum))
}
j <- 0
for (mea in c("aic", "sic", "rmse", "crps")) {
j <- j + 1
h <- helper_GetIndices(search_res[[mea]], result_sum[[mea]])
res$model[[i, j]] <- res$model[[i, j]] + GetScore(h$model1, h$model2)
res$inclusion[[i, j]] <- res$inclusion[[i, j]] + GetScore(h$inclusion1, h$inclusion2)
}
return(res)
}Finally, we run the simulation:
simCount <- 2
sur_sim_res <- list()
for (n in sampleSizes) {
res <- list(model = matrix(0, length(sMultipliers), length(measureNames),
dimnames = list(paste0("sigma=", sMultipliers), measureNames)
))
res$inclusion <- res$model
i <- 0
for (s in sMultipliers) {
i <- i + 1
for (iter in seq(1, simCount)) {
seedi <- seed + iter
measureOptions$seed <- seedi
measureOptions$simFixSize <- 6
measureOptions$trainRatio <- 0.75
searchItems <- GetSearchItems(model = TRUE, type1 = FALSE, bestK = 10, inclusion = TRUE)
res <- SampleAndEstimate(
i, res, n, s, xSizes, xCounts, yGroups, measureOptions,
searchItems, NULL
)
}
}
res$model <- res$model / simCount
res$inclusion <- res$inclusion / simCount
sur_sim_res[[length(sur_sim_res) + 1]] <- list(n = n, sMultipliers = sMultipliers, res = res)
} You can increase the number of simulations (We choose a small number
to decrease the computation time). Note that for each
n\(\in\)sampleSizes and each
s\(\in\)sMultipliers, we sample
and estimate 2 times. This is relatively time-consuming.
Result
Comparing the performance of different measures in selecting the true explanatory variables can be interesting. Therefore, in the following plot we draw the scores of the best models:
Scores in finding the true model (n=30, 50)
There is a similar result for the inclusion weights:
Scores in finding the true inclusion variables (n=30, 100, 100)
