Testing for https://github.com/ModelOriented/ingredients/pull/142

feature_importance_inputData.r

# install.packages(c("abind", "data.table"))
library(keras)

df<- data.frame(id=rep(LETTERS[1:10], each=5), static=rep(1:10, each=5), time=rep(1:5, times=5))
df.cat<- data.frame(id=LETTERS[1:10], cat1=rep(LETTERS[1:5], times=2), cat2=letters[1:10])
df<- merge(df, df.cat)
df$x1<- df$time * df$static
df$x2<- rnorm(nrow(df), mean=df$time * df$static + 10, sd=5)
df$x3<- rnorm(nrow(df), mean=df$time * df$static * 3, sd=2)
df$y<- rnorm(nrow(df), mean=(df$x1 + df$x2) / df$x3, sd=2)

timevar<- "time"
idVars<- "id"
responseVars<- "y"
staticVars<- c("static", "cat1", "cat2")
predTemp<- c("x1", "x2", "x3")
responseTime<- max(df[, timevar], na.rm=TRUE)
regex_time<- "[0-9]+"
perm_dim = 2:3
hidden_shape.RNN<- 32
hidden_shape.static<- 16
hidden_shape.main<- 32
epochs<- 3
batch_size<- 1000
verbose<- 0  # 2

wideTo3Darray.ts<- function(d, vars, idCols){
  d<- as.data.frame(d)
  timevals<- unique(gsub(paste0("^(", paste(vars, collapse="|"), ")_"), "", setdiff(colnames(d), idCols)))

  # Reshape to a 3D array [samples, timesteps, features] Format for RNN layers in NN
  a<- lapply(vars, function(x){
    varTS<- d[, grep(paste0("^", x, "_(", paste(timevals, collapse="|"), ")$"), colnames(d))]
    a<- array(as.matrix(varTS), dim=c(nrow(varTS), ncol(varTS), 1), dimnames=list(case=NULL, t=gsub(paste0("^", x, "_"), "", colnames(varTS)), var=x))
  })
  names(a)<- vars
  a<- abind::abind(a)
  names(dimnames(a))<- c("case", "t", "var")
  dimnames(a)$case<- do.call(paste, c(d[, idCols, drop=FALSE], list(sep="_")))

  return(a)
}

build_modelLTSM<- function(input_shape.ts, input_shape.static=0, output_shape=1,
                           hidden_shape.RNN=32, hidden_shape.static=16, hidden_shape.main=32){
  inputs.ts<- layer_input(shape=input_shape.ts, name="TS_input")
  inputs.static<- layer_input(shape=input_shape.static, name="Static_input")

  predictions.ts<- inputs.ts
  for (i in 1:length(hidden_shape.RNN)){
    predictions.ts<- predictions.ts %>% layer_lstm(units=hidden_shape.RNN[i], name=paste0("LSTM_", i))
  }

  if (input_shape.static > 0){
    predictions.static<- inputs.static
    for (i in 1:length(hidden_shape.static)){
      predictions.static<- predictions.static %>% layer_dense(units=hidden_shape.static[i], name=paste0("Dense_", i))
    }
    output<- layer_concatenate(c(predictions.ts, predictions.static))
  } else {
    output<- predictions.ts
  }

  for (i in 1:length(hidden_shape.main)){
    output<- output %>% layer_dense(units=hidden_shape.main[i], name=paste0("main_dense_", i))
  }
  output<- output %>% layer_dense(units=output_shape, name="main_output")

  if (input_shape.static > 0){
    model<- keras_model(inputs=c(inputs.ts, inputs.static), outputs=output)
  } else {
    model<- keras_model(inputs=inputs.ts, outputs=output)
  }

  compile(model,
          loss="mse",
          optimizer=optimizer_rmsprop(),
          metrics=list("mean_squared_error", "mean_absolute_error", "mean_absolute_percentage_error")
  )

  model
}

train_keras<- function(modelNN, train_data, train_labels, test_data, test_labels, epochs, batch_size, callbacks=NULL, verbose=0){
  validation_data<- list(test_data, test_labels)
  history<- keras::fit(object=modelNN,
                       x=train_data,
                       y=train_labels,
                       batch_size=ifelse(batch_size %in% "all", nrow(train_data), batch_size),
                       epochs=epochs,
                       verbose=verbose,
                       callbacks=callbacks,
                       validation_data=validation_data
  )

  return(modelNN)
}

predVars<- setdiff(colnames(df), c(idVars, timevar))
predVars.cat<- names(which(!sapply(df[, predVars, drop=FALSE], is.numeric)))
predVars.num<- setdiff(predVars, predVars.cat)

df.catBin<- stats::model.matrix(stats::as.formula(paste("~ -1 +", paste(predVars.cat, collapse="+"))), data=df)
predVars.catBin<- colnames(df.catBin)
df<- cbind(df[, setdiff(colnames(df), predVars.cat)], df.catBin)
predVars<- c(predVars.num, predVars.catBin)
staticVars.cat<- staticVars[staticVars %in% predVars.cat]
staticVars<- c(setdiff(staticVars, staticVars.cat), predVars.catBin)

# crossvalidation for timeseries must be done in the wide format data
# scaling must be done on long format data
responseVars.ts<- paste0(responseVars, "_", responseTime)
predVars.tf<- paste0(setdiff(predVars, staticVars), "_", responseTime)

## df to wide format
dt<- data.table::as.data.table(df)
staticCols<- c(idVars, staticVars)
vars<- setdiff(colnames(dt), c(staticCols, timevar))
timevals<- unique(data.table:::`[.data.table`(x=dt, , j=timevar, with=FALSE))[[1]] # without importing data.table functions
LHS<- setdiff(staticCols, timevar)
form<- paste0(paste(LHS, collapse=" + "), " ~ ", timevar)
dt<- data.table::dcast(dt, formula=stats::formula(form), value.var=vars)  # To wide format (var_time columns)
df.wide<- as.data.frame(dt)
predVars.ts<- setdiff(colnames(df.wide), c(idVars, staticVars)) # WARNING: Includes responseVars.ts
timevals<- unique(df[[timevar]])

idxTrain<- sample(1:nrow(df.wide), 6)

train_labels<- df.wide[idxTrain, c(idVars, responseVars.ts), drop=FALSE]
train_data<- df.wide[idxTrain, c(idVars, staticVars, predVars.ts), drop=FALSE]

test_labels<- df.wide[-idxTrain, c(idVars, responseVars.ts), drop=FALSE]
test_data<- df.wide[-idxTrain, c(idVars, staticVars, predVars.ts), drop=FALSE]


# Reshape data to 3D arrays [samples, timesteps, features] as expected by LSTM layer
train_data.3d<- wideTo3Darray.ts(d=train_data, vars=setdiff(predVars, staticVars), idCols=idVars)
test_data.3d<- wideTo3Darray.ts(d=test_data, vars=setdiff(predVars, staticVars), idCols=idVars)

train_data.static<- structure(as.matrix(train_data[, staticVars, drop=FALSE]),
                              dimnames=list(case=do.call(paste, c(train_data[, idVars, drop=FALSE], list(sep="_"))), var=staticVars))
test_data.static<- structure(as.matrix(test_data[, staticVars, drop=FALSE]),
                             dimnames=list(case=do.call(paste, c(test_data[, idVars, drop=FALSE], list(sep="_"))), var=staticVars))

train_labels<- structure(as.matrix(train_labels[, responseVars.ts, drop=FALSE]),
                         dimnames=list(case=do.call(paste, c(train_labels[, idVars, drop=FALSE], list(sep="_"))), var=responseVars.ts))
test_labels<- structure(as.matrix(test_labels[, responseVars.ts, drop=FALSE]),
                        dimnames=list(case=do.call(paste, c(test_labels[, idVars, drop=FALSE], list(sep="_"))), var=responseVars.ts))

train_data.3d<- train_data.3d[, setdiff(dimnames(train_data.3d)[[2]], responseTime), ]
test_data.3d<- test_data.3d[, setdiff(dimnames(train_data.3d)[[2]], responseTime), ]

train_data<- list(TS_input=train_data.3d, Static_input=train_data.static)
test_data<- list(TS_input=test_data.3d, Static_input=test_data.static)

early_stop<- keras::callback_early_stopping(monitor="val_loss", patience=30)

modelNN<- build_modelLTSM(input_shape.ts=dim(train_data.3d)[-1], input_shape.static=length(staticVars), output_shape=length(responseVars),
                                    hidden_shape.RNN=hidden_shape.RNN, hidden_shape.static=hidden_shape.static, hidden_shape.main=hidden_shape.main)
modelNN<- train_keras(modelNN=modelNN, train_data=train_data, train_labels=train_labels,
                                test_data=test_data, test_labels=test_labels,
                                epochs=epochs, batch_size=batch_size, callbacks=early_stop, verbose=verbose)

modelNN.LSTM<- build_modelLTSM(input_shape.ts=dim(train_data.3d)[-1], input_shape.static=0, output_shape=length(responseVars),
                                         hidden_shape.RNN=hidden_shape.RNN, hidden_shape.static=0, hidden_shape.main=hidden_shape.main)
modelNN.LSTM<- train_keras(modelNN=modelNN.LSTM, train_data=train_data.3d, train_labels=train_labels,
                                     test_data=test_data.3d, test_labels=test_labels,
                                     epochs=epochs, batch_size=batch_size, callbacks=early_stop, verbose=verbose)

## Variable importance ----
# explainer<- DALEX::explain(model=modelNN, data=test_data, y=test_labels, predict_function=stats::predict, label="MLP_keras", verbose=FALSE)
# FOR DEBUG
y=test_labels
loss_function = DALEX::loss_root_mean_square
type = c("raw", "ratio", "difference")[1]
predict_function = stats::predict
n_sample = NULL
B = 10
N = n_sample
label = "keras"
perm_dim = NULL
comb_dims = FALSE
variables = NULL
variable_groups = NULL


### Feature_importance: LSTM only ----
## NOTE: REQUIRES array_feature_importance branch
x = modelNN.LSTM
data = test_data.3d

ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, permDim=2:3)


### Feature_importance: LSTM + static ----
## NOTE: REQUIRES multiinput branch
x = modelNN
data = test_data

v_groups.ts <- mapply(function(dimVar, dimNames) {
  v<- lapply(dimVar, function(v) setNames(list(v), dimNames))
  setNames(v, nm = dimVar)
}, dimVar=dimnames(data$TS_input)[-1], dimNames=names(dimnames(data$TS_input))[-1], SIMPLIFY=FALSE)
v_groups.ts <- do.call(c, v_groups.ts)

v_groups.tsCombDim <- expand.grid(dimnames(data$TS_input)[-1], stringsAsFactors=FALSE, KEEP.OUT.ATTRS=FALSE) # All combinations for all dimensions in a dataset
rownames(v_groups.tsCombDim) <- apply(v_groups.tsCombDim, 1, function(v) paste(v, collapse="|"))
v_groups.tsCombDim <- split(v_groups.tsCombDim, rownames(v_groups.tsCombDim))
v_groups.tsCombDim <- lapply(v_groups.tsCombDim, as.list)

v_groups.static<- list(list("static"),
                       list(grep("^cat1", dimnames(data$Static_input)$var, value=TRUE)),
                       list(grep("^cat2", dimnames(data$Static_input)$var, value=TRUE)))
names(v_groups.static)<- c("static", "cat1", "cat2")
variable_groups<- list(TS_input=v_groups.ts, Static_input=v_groups.static)
variable_groups.combDim<- list(TS_input=v_groups.tsCombDim, Static_input=v_groups.static)

variable_groups.noGrNames<-lapply(variable_groups, function(input){
  names(input)<- NULL
  input
})
variable_groups.noVarNames<-lapply(variable_groups, function(input){
  lapply(input, function(gr){
    names(gr)<- NULL
    gr
  })
})
variable_groups.noGrVarNames<-lapply(variable_groups, function(input){
  names(input)<- NULL
  lapply(input, function(gr){
    names(gr)<- NULL
    gr
  })
})

variable_groups.combDim.noGrNames<-lapply(variable_groups.combDim, function(input){
  names(input)<- NULL
  input
})
variable_groups.combDim.noVarNames<-lapply(variable_groups.combDim, function(input){
  lapply(input, function(gr){
    names(gr)<- NULL
    gr
  })
})
variable_groups.combDim.noGrVarNames<-lapply(variable_groups.combDim, function(input){
  names(input)<- NULL
  lapply(input, function(gr){
    names(gr)<- NULL
    gr
  })
})

ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type)
ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, comb_dims=TRUE)
ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, perm_dim=list(TS_input=3, Static_input=2))

ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups)
ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.noGrNames)
# ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.noVarNames)
ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.noGrVarNames)

ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.combDim)
ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.combDim.noGrNames)
# ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.combDim.noVarNames)
ingredients::feature_importance(x=x, data=data, y=y, loss_function=loss_function, type=type, variable_groups=variable_groups.combDim.noGrVarNames)