I have a binary dataset of infected trees consisting of a lot of independent variables that can be divided into three main groups of drivers. Let's intuitively call them
inherentdriver consist of variables such as:
stand age, etc.
staticdriver consist of variables such as:
dynamicdriver consist of variables such as:
annual average growth rateetc.
I have performed logistic ElasticNet regression for all variables included, and now, I would like to determine the relative or absolute importance of the three groups of drivers in my regression. I have been suggested to run the exact same regression for each of the drivers individually. But in my mind, this is a wrongful test of the drivers as there may be an interaction between variables across drivers.
How would you approach this from a statistical standpoint?
Here are (some) data and code for running ElasticNet regression
df.head() infected age LAI avg_temp avg_moist precip species elevation area <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <dbl> <fct> 0 15 2.46 25.0 6.8 989 4 4.66 1 0 13 1.50 18.3 10.7 631 3 11.12 3 0 21 3.80 10.5 25.8 1207 2 14.73 1 0 56 5.21 24.2 9.2 434 1 3.21 2 0 57 4.31 20.6 10.4 499 1 4.63 2 0 2 0.58 25.3 2.1 801 4 2.58 1 set.seed(123) library(caret) library(tidyverse) library(glmnet) library(ROCR) library(doParallel) registerDoParallel(4, cores = 8) training.samples <- df$infected %>% createDataPartition(p = 0.8, list = FALSE) train <- df[training.samples, ] test <- df[-training.samples, ] x.train <- data.frame(train[, names(train) != "infected"]) x.train <- data.matrix(x.train) y.train <- train$infected x.test <- data.frame(test[, names(test) != "infected"]) x.test <- data.matrix(x.test) y.test <- test$infected model <- cv.glmnet(x.train, y.train, type.measure = "auc", alpha = i/10, family = "binomial", parallel = TRUE) coef <- coef(model, s = model$lambda.min) predicted <- predict(model, s = model$lambda.min, newx = x.test, type = “response”) auc <- max(model$cvm) (Intercept) -1.060656e+01 age -5.851591e-03 LAI 2.715723e-01 avg_temp . avg_moist -6.587949e-04 precip . species . elevation -1.093533e-03 area -1.180528e-01
Previously, people have suggested
plot(anova(model), what='proportion chisq') here, but then I need to run the
model separately for each driver and the interaction is lost, right?