機械学習や分析時のカテゴリ変数の高度なあれこれ

2*sum(train.dt[, .(importance=abs(sum(y) - .N*train.dt[, (sum(y)/.N)])), by=x1][,importance])/train.dt[,.N]

remain_x1 <- train.dt[, .(importance=abs(sum(y) - .N*train.dt[, (sum(y)/.N)])), by=x1][order(-importance), as.character(x1)][1:100]

train.dt[!is.na(x1) & !as.character(x1) %in% remain_x1, x1 := '(other)']
train.dt[,x1:=droplevels(x1)]

bin_encode <- function(vec, label){
  mat <- matrix(
    as.integer(intToBits(as.integer(vec))),
    ncol = 32,
    nrow = length(vec),
    byrow = TRUE
  )[, 1:ceiling(log(length(unique(vec)) + 1)/log(2))]
  if (!is.matrix(mat)) mat <- as.matrix(mat)
  colnames(mat) <- paste0(label, seq(ncol(mat)))
  return(mat)
}

# Jaccard係数を計算する関数。水準数が多い場合を想定してsparse matrix(dgCMatrix)を使う。
require(Matrix)
jcdSparse <- function(X){
  X <- as(X, "dgCMatrix")
  X@x <- rep(1, length(X@x)) # 非ゼロ要素はすべて1に置換
  vec_sum_col <- colSums(X)
  AB <- crossprod(X)
  sim <- AB / (vec_sum_col - t(AB - vec_sum_col)) # まとめて
  sim <- as(sim, "dgCMatrix")
  return(sim)
}

# Jaccard係数が閾値以上の水準数の全水準数に対する比率を計算する関数。
factors_almost_equal <- function(data, x, y, threshold=0.95) {
  x <- deparse(substitute(x))
  y <- deparse(substitute(y))
  if (is.factor(data[[x]]) == FALSE) {
    #stop(paste(x, 'is not a factor!'))
    data[[x]] <- as.factor(as.character(data[[x]]))
  }
  if (is.factor(data[[y]]) == FALSE) {
    #stop(paste(y, 'is not a factor!'))
    data[[y]] <- as.factor(as.character(data[[y]]))
  }
  x_mat <- fac2sparse(data[[x]], drop=F)
  y_mat <- fac2sparse(data[[y]], drop=F)
  sim <- jcdSparse(t(rbind(x_mat, y_mat)))
  sim <- tril(sim, -1) + triu(sim, 1)
  ratio_gt_threshold <- sum(apply(sim, 2, max) > threshold) / nrow(sim)
  #print(paste('Similarity between', x, 'and', y, ':', ratio_gt_threshold))
  return(ratio_gt_threshold)
}

factor_arrange_all <- function(df1, df2) {
  .col_factor <- intersect(
    names(which(sapply(df1, is.factor))),
    names(which(sapply(df2, is.factor)))
  )
  for (.col in .col_factor) {
    new_level <- union(levels(df1[[.col]]), levels(df2[[.col]]))
    new_level <- sort(new_level, na.last=F)
    df1[[.col]] = factor(as.character(df1[[.col]]), levels=new_level, exclude=NULL)
  }
  return(df1)
}

train.dt[,x1:=factor(x1, levels=c(levels(x1), '(other)'))]
train.df$x1 <- factor(train.df$x1, levels=c(levels(train.df$x1), '(other)'))

factor_arrange_reduced <- function(df1, df2) {
  .col_factor <- intersect(
    names(which(sapply(df1, is.factor))),
    names(which(sapply(df2, is.factor)))
  )
  for (.col in .col_factor) {
    new_level <- levels(df2[[.col]])
    new_level <- sort(new_level)
    df1[[.col]] = factor(as.character(df1[[.col]]), levels=new_level)
  }
  return(df1)
}

converted_predict.df <- factor_arrange_reduced(predict.df, train.df)

saveRDS(train.df[0], file='factor_template.rds')

factor_template.df <- readRDS(file='factor_template.rds')
converted_predict.df <- factor_arrange_reduced(predict.df, factor_template.df)

{1,2,3,NA}
```r

という4水準があるとする。順序尺度として扱うのもありだが（その場合NAが厄介）、単純なfactor型（順序を持たないカテゴリ変数）として扱ってもいい。その場合は**ベース水準をNAにしておく**と便利である。
NAをベース水準にすると、水準1の回帰係数の意味は「無回答と比較した時の1の効果」ということになる。ベース水準は「○○と比較した時の」の「○○」になる。この○○が「無回答」のほうが各選択肢の意味を解釈しやすくなる（無回答の頻度にもよるが）。
以下の関数`factor_add_na()`でNAを水準として追加し、それをベース水準として設定する。`x`は水準化されていないNAを含むfactor型変数。`drop`が`TRUE`の場合、使用されていない水準が取り除かれる。

`factor()`関数の引数`exclude`で水準対象外にするラベルを指定できる。これがデフォルトではNAが採用されているので、NULLを指定することで水準対象外のラベルがなくなる。そして`factor()`関数の引数`levels`の先頭がベース水準になるので、そこにNAを持ってくる。

### カテゴリ変数におけるNAの意味

上記のようにNAを水準として設定することができるが、そもそもfactor型におけるNAは複数の意味を持つことがある。
既存のfactor型変数に対して、水準に含まれない値を追加代入するとNAに変換される。つまりfactor型におけるNAは欠測の場合もあれば**未知の値**を意味することがある。これらは全く意味が異なるわけで、カテゴリ変数におけるNAは欠損（missing）とその他（other）に分けて扱うべきである。
(other)が発生する前のデータを取り込んだ直後の段階で以下の処理を実行し、NAを消しておく。欠損値をすべて(missing)にする。

この後発生したNAはすべて(other)ということになるので、最後に

とすればいい。

### まとめて実装した

学習時 `train.dt`が対象のdata.table

- `x1`はあらかじめGini係数に基づいて水準数を絞り込んでおく
- `x1`を含むすべてのfactor列のNAを(missing)と(other)に分け、水準化する

```r
# Gini係数上位100水準を抜き出しておく
remain_x1 <- train.dt[, .(importance=abs(sum(y) - .N*train.dt[, (sum(y)/.N)])), by=x1][order(-importance), as.character(x1)][1:100]
# 上位100水準に含まれればその値を採用し、含まれなければ(other)にする
train.dt[!is.na(x1) & !as.character(x1) %in% remain_x1, x1 := '(other)']
# 上の処理で不要になった水準を削除
train.dt[,x1:=droplevels(x1)]

# すべてのfactor型の列名を抽出
.col_convert <- which(sapply(train.dt, is.factor))
# factor列におけるNAを(missing)という値に変換
train.dt[, (.col_convert):=lapply(.SD, function(x) as.factor(ifelse(is.na(x), '(missing)', as.character(x)))), .SDcols=.col_convert]

# factor型の水準に(other)と(missing)を入れておく（上で追加されている場合は特になし）。ベース水準に(missing)を設定
train.dt[, (.col_convert):=lapply(.SD, function(x) factor(as.character(x), levels=sort(unique(c('(missing)', '(other)', levels(x)))))), .SDcols=.col_convert]

# テンプレートdata.tableとして保存
saveRDS(train.dt[0], file='factor_template.rds')

# 学習時に抽出した上位100水準に含まれればその値を採用し、含まれなければ(other)にする
predict.dt[!is.na(x1) & !x1 %in% remain_x1, x1 := '(other)']
# 上の処理で不要になった水準を削除
predict.dt[,x1:=droplevels(x1)]

# すべてのfactor型の列名を抽出
.col_convert <- which(sapply(predict.dt, is.factor))
# factor列におけるNAを(missing)という値に変換
predict.dt[, (.col_convert):=lapply(.SD, function(x) as.factor(ifelse(is.na(x), '(missing)', as.character(x)))), .SDcols=.col_convert]

# 学習データと列や水準を合わせるためにテンプレートを読み込む
factor_template.dt <- readRDS(file='factor_template.rds')
# 学習データにない列は削除
predict.dt <- predict.dt[,.SD, .SDcols=colnames(factor_template.dt)]
# factor型の水準を学習データに合わせる（ない場合はNA->(other)と変換する）
converted_predict.dt <- factor_arrange_reduced(predict.dt, factor_template.dt)

# factor型の水準に(other)と(missing)を入れておく（上で追加されている場合は特になし）。ベース水準に(missing)を設定
converted_predict.dt[, (.col_convert):=lapply(.SD, function(x) factor(ifelse(is.na(x),'(other)',as.character(x)), levels=sort(unique(c('(missing)', '(other)', levels(x)))))), .SDcols=.col_convert] # as.factor()だと値のない水準が値のない水準がdropされてしまうので、factor()でlevelsを指定する

import pandas as pd

# 水準数が膨大なカテゴリ変数から重要でない水準を削除
def preprop_reduce_levels(df, key, remain):
  remain = list(map(str, remain))
  categories = sorted(['(missing)', '(other)'] + remain)
  dtype = pd.CategoricalDtype(categories=categories)
  df[key] = (
    df[key].astype('string')
    .fillna('(missing)')
    .apply(lambda x: x if x in remain else '(other)')
    .astype(dtype)
  )
  return df

# 後で予測データと合わせることができるようにするための処理
def format_levels(df):
  _cols = df.select_dtypes(include='category').columns
  for _col in _cols:
    existing_cats = list(df[_col].cat.categories)
    new_cats = sorted(set(existing_cats + ['(missing)', '(other)']))
    dtype = pd.CategoricalDtype(categories=new_cats)
    df[_col] = (
      df[_col].astype('string')
      .fillna('(missing)')
      .str.strip('][>')
      .astype(dtype)
    )
  return df

# 学習データに合わせるための処理
def modify_levels_to_train_data(df_data, df_template):
  # 学習データにない列は削除
  df_data = df_data[df_template.columns]

  # カテゴリの水準を学習データに合わせる
  _cols = df_data.select_dtypes(include='category').columns
  for _col in _cols:
    dtype = pd.CategoricalDtype(
      categories=list(df_template[_col].cat.categories)
    )
    try:
      df_data[_col] = df_data[_col].astype(dtype)
    except Exception as e:
      logger.warn((_col, e))

  df_data[_cols] = df_data[_cols].fillna('(other)')

  return df_data

# remain_x1で指定した上位の水準に含まれればその値を採用し、含まれなければ(other)にする
train_df = preprop_reduce_levels(train_df, 'x1', remain_x1)
# NAを削除して(missing)と(other)を水準に追加
train_df = format_levels(train_df)
# カテゴリ水準テンプレートの保存
train_df.iloc[0:0].to_pickle('category_template.pkl', compression='zip')

# remain_x1で指定した上位の水準に含まれればその値を採用し、含まれなければ(other)にする
predict_df = preprop_reduce_levels(predict_df, 'x1', remain_x1)
# 学習データと列や水準を合わせるためにテンプレートを読み込む
category_template_df = pd.read_pickle('category_template.pkl', compression='zip')
# 列とカテゴリ変数の水準を学習データに合わせる
predict_df = modify_levels_to_train_data(predict_df, category_template_df)