My questions is this right approach to do feature selection when data volume is high?
Simply, no.
Basing feature selection on p values is a bad idea, especially when data are large. First, p-values tell you nothing about the effect of the variable. I can always construct a model with a highly significant feature but which performs negligibly different with respect to any classification metric you choose. This is because significant effects can be extremely small.
When data is large, the null is essentially a straw man. You have so much data that you can detect small effects because you have immense power to do so. The effect of any variable is never exactly 0 and you are finding that.
My advice is to use some principled modelling approach. People seem to like AIC (I'm not one of them), you could do forward feature selection (again, not my cup of tea), you could do lasso or ridge regression (I'm more keen on this), or frankly you could do none of them (my preference from what you've said in your post). If you have 12 variables which you know to be important, why aren't you using all of them? That's a rhetorical question.
In short, inference breaks down when you have so much data. The null becomes a straw man, so you reject near everything. People's obsession with p values leads to them using p values for things which they were not intended for (model selection). You should lean on methods which evaluate what you care about via a validation set or lean on your business knowledge.
EDIT:
I claim I can always make a model which performs negligibly better even when the p value is significant. Here is an example using linear regression:
library(tidyverse)
library(Metrics)
set.seed(0)
X = rnorm(1000000)
Z = rnorm(1000000)
y = 2*X + 0.01*Z + rnorm(1000000, 0, 0.3)
d = tibble(X = X, Z = Z, y = y, set =
sample(c('test','train'), replace = T,
size = 1000000))
test = filter(d, set=='test')
train = filter(d, set=='train')
model1 = lm(y~X + Z, data = train)
model2 = lm(y~X, data= train)
rmse(test$y, predict(model1, newdata = test))
#> [1] 0.2996978
rmse(test$y, predict(model2, newdata = test))
#> [1] 0.2998523
Created on 2022-01-06 by the reprex package (v2.0.1)
The rmse for both models agrees up to 3 decimal places. That is good for all intents and purposes in my opinion. Note that the coefficient for Z is highly significant (it gives the smallest p value R can give). The combination of tiny effect size and massive sample is what causes this phenomena.