XgBoost Simulation Relative Importance I am using Xgboost to find the relative contribution of features to the dependent variables. I am running a simulation wherein I'm generating a simulated data-set based on a prespecified mean and variance. There are 6 independent variables X1, X2, X3, X4, X5, and X6. The dependent variable dep is an average of these 6 variables. Therefore, XGBOOST should output a rel contribution of 1/6 for each variable. However, this is not what I am getting when I am running the code. I have attached the code for your reference. Can anyone tell where am I going wrong?
#Creating the data set
set.seed(123456789)
cov_sim1 <- matrix(c(857.7941,0,0,0,0,0,0,15065.39,0,0,0,0,0,0,18.38235,0,0,0,0,0,0,985.8595,0,0,0,0,0,0,35719.32,0,0,0,0,0,0,174.683), nrow = 6, ncol =6)
mean_sim <- c(-2.368421, 3.000000, 1.052632, 18.631579, -33.210526, -5.473684) 
dat <- MASS::mvrnorm(19, mu = mean_sim, Sigma = cov_sim1, empirical = TRUE)
simulated_set <- data.frame(dat)
simulated_set$dep <-(simulated_set$X1 + simulated_set$X2 + simulated_set$X3 + simulated_set$X4 + simulated_set$X5 + simulated_set$X6)/6

# Setting parameter and tuning
xgb_trcontrol <- caret::trainControl(
  method = "cv", 
  number = 2,
  allowParallel = TRUE, 
  verboseIter = FALSE, 
  returnData = FALSE
)

obj <- 'reg:squarederror'

xgb_grid <- base::expand.grid(
  list(
    nrounds = c(100, 200),
    max_depth = c(2, 4, 6), # maximum depth of a tree
    colsample_bytree = seq(from=0.1, to=1, by=0.1), # subsample ratio of columns when construction each tree
    eta = seq(from=0.1, to=1, by=0.1),    
    gamma = 0, # minimum loss reduction,
    min_child_weight = 1,  # minimum sum of instance weight (hessian) needed ina child
    subsample = 1    # subsample ratio of the training instances
  ))

#Fitting the model and extracting feature importance

d <- data.matrix(simulated_set[1:6])
set.seed(10302)
xgb_model <- caret::train(
  d, simulated_set$dep,
  trControl = xgb_trcontrol,
  tuneGrid = xgb_grid,
  method = "xgbTree",
  nthread = 10,
  verbose = FALSE,
  metric = "RMSE",
  verbosity = 0
)
bestfit <- xgb_model$bestTune
nrounds_bestfit <- bestfit[1,1]
max_depth_bestfit <- bestfit[1,2]
eta_bestfit <- bestfit[1,3]
colsample_bytree_bestfit <- bestfit[1,4]
min_child_weight_bestfit <- bestfit[1,5]
subsample_bestfit <- bestfit[1,6]
bst <- xgboost(booster = "gbtree", data = d, label = simulated_set$dep, max.depth = max_depth_bestfit, 
           nround = nrounds_bestfit, eta = eta_bestfit,
             min_child_weight <- min_child_weight_bestfit)
importance <- xgb.importance(model = bst)
importance <- importance[order(importance$Feature),]
ac <- as.data.frame(importance)

I have tried to standardize the variables as well.  That does not change the output by much? My initial thinking is that it depends on the variance of the dependent variables. But, not sure.
 A: The differences in variance are most definitely the reason for the observed differences in importance.
Broadly, you would only expect the result to show ~.167 for a lm's set of coefficients.  Because of the way simulated_set[["dep"]] is defined, the mean value is determined by mulitplying all X variables by one-sixth.
> sapply(coef(lm(dep ~ X1 + X2 + X3 + X4 + X5 + X6, simulated_set)), round, digits=5)
(Intercept)          X1          X2          X3          X4          X5 
    0.00000     0.16667     0.16667     0.16667     0.16667     0.16667 
         X6 
    0.16667 

Importance statistics usually involve model fit, and xgboost is no exception.  This idea is easier to see for simpler metrics though like lm-based importance metrics such as Dominance Analysis/DA which is explained well in the linked article by Grömping.
Take the same model and run it through a DA (implemented by domir but equivalent results can be obtained by dominanceanalysis and relaimpo) and you get:
> sapply(domir::domin(dep ~ X1 + X2 + X3 + X4 + X5 + X6, lm, list(summary, "r.squared"), data=simulated_set)$General_Dominance, round, digits=5)
     X1      X2      X3      X4      X5      X6 
0.01624 0.28521 0.00035 0.01866 0.67623 0.00331 

> sapply(simulated_set[-7], sd)
        X1         X2         X3         X4         X5         X6
 29.288122 122.741150   4.287464  31.398400 188.995556  13.216770

Very much dependent on the variances as can be seen.  Whereas this article focuses on random forest, a comparison between linear models and random forest is discussed here; it provides good description of similarities and differences between the approaches that might help to explain the results.
