Разные способы нормализации данных и оценка их успешности

В пакете caret есть функция preProcess, в которой реализован ряд методов препроцессинга данных: “center”, “scale”, “BoxCox”, “YeoJohnson”, “expoTrans”, “range”, “knnImpute”, “bagImpute”, “medianImpute”, “pca”, “ica”, “spatialSign”, “zv”, “nzv”, “conditionalX”.

Из них “bagImpute”, “pca”, “ica”, “spatialSign” у меня не получилось использовать (например, “pca” и “ica” нужны для одновременного преобразования групп переменных, а не по одной), а что с “bagImpute” и “spatialSign” пошло не так я так и не понял…

Кроме того, из книги А.Б. Шипунова и др. “Наглядная статистика. Используем R!” я почерпнул еще несколько методов:

• Логарифмическое: log(data + 1). Если распределение скошено вправо, может дать нормальное распределение. Может также делать более линейными зависимости между переменными и уравнивать дисперсии. «Боится» нулей в данных, поэтому рекомендуется прибавлять единицу.
• Квадратного корня: sqrt(data). Похоже по действию на логарифмическое. «Боится» отрицательных значений. Обратное: 1/(data + 1). Эффективно для стабилизации дисперсии. «Боится» нулей.
• Квадратное: data^2. Если распределение скошено влево, может дать нормальное распределение. Линеаризует зависимости и выравнивает дисперсии.
• Логит: log(p/(1-p)). Чаще всего применяется к пропорциям. Линеаризует так называемую сигмовидную кривую.
• Кроме логит-преобразования, для пропорций часто используют и арксинус-преобразование, asin(sqrt(p))

Ряд перечисленных в книге преобразований не работает с отрицательными значениям. Поэтому я сначала шкалирую переменную (scale(x)), а затем добавляю 3, чтобы > 99% значений стали положительными, а затем применяю вышеозначенные дополнительные преобразования.

После всех трансформаций я снова провожу финальное шкалирование (scale(x)).

Функция строит график типа barplot с нанесенными поверх него точками данных, а также на графике подписывается значение p-level для теста shapiro.test(x), с помощью которого можно определить нормальность распределения трансформированных данных. Нулевая гипотеза этого теста заключается в том, что данные распледелены нормально (подписи красного цвета), а если p-level < 0.1, то значит эта гипотеза отвергается и распределение нельзя считать нормальным (подписи серого цвета).

Как правило, метод BoxCox справляется лучше всего с задачей.

compareTransformations = function(var.nonscaled){

  require(caret)
  require(ipred)
  
  if (is.data.frame(var.nonscaled)) var.nonscaled = var.nonscaled[[1]]
  if (is.matrix(var.nonscaled)) var.nonscaled = as.numeric(var.nonscaled)
  
  var.scaled = scale(var.nonscaled)
  len = length(var.nonscaled)
  
  # different types of transformations
  pp.list = c( "BoxCox", "YeoJohnson", "expoTrans", "range", "knnImpute", "medianImpute", "zv", "nzv", "conditionalX")
  
  dt = data.frame(transform=rep("scaled",len), data=var.scaled)
  
  for (pp in pp.list){
    pp.obj = preProcess(x=data.frame(x=var.nonscaled), method=pp)
    pp.res = predict(pp.obj, data.frame(x=var.nonscaled))
    dt = rbind(dt, data.frame(transform=rep(pp,len),data=scale(pp.res[,1])))
  }
  
  # log
  dt = rbind(dt, data.frame(transform=rep("log",len),data=scale(log(var.scaled+3))))
  # sqrt
  dt = rbind(dt, data.frame(transform=rep("sqrt",len),data=scale(sqrt(var.scaled+3))))
  # 1/(data+1)
  dt = rbind(dt, data.frame(transform=rep("1/(data+1)",len),data=scale(1/(var.scaled+1))))
  # data^2
  dt = rbind(dt, data.frame(transform=rep("data^2",len),data=scale(var.scaled^2)))
  # logit
  dt = rbind(dt, data.frame(transform=rep("logit",len),data=scale((var.scaled+3)/(1-(var.scaled+3)))))
  # arcsin
  dt = rbind(dt, data.frame(transform=rep("arcsin",len),data=scale(asin(scale(sqrt(var.scaled+3))/3))))
  
  dt.sh = c()
  dt.list = list()
  for (pp in unique(dt[,1])){
    dt.list[[pp]] = dt[dt[,1]==pp,2]
    dt.sh = c(dt.sh,sprintf("%.3f",shapiro.test(dt.list[[pp]])$p.value))
  }
  dt.sh.col = numeric(length(dt.sh))
  dt.sh.col[dt.sh>0.1] = "red"
  dt.sh.col[dt.sh<=0.1] = "gray"
  
  
  g = ggplot(dt,aes(x=transform, y=data)) +
    geom_boxplot(fill="cyan", size=.9) +
    geom_jitter(color="blue",alpha=0.3) +
    theme_bw(base_size = 16) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=0.2)) +
    geom_text(data=data.frame(), aes(x=1:length(unique(dt[,1])), y=rep(-4, length(unique(dt[,1]))), label=dt.sh),col=dt.sh.col) +
    coord_cartesian(ylim=c(-5,5))
  print(g)
  
  dt.list
}

Вот пример использования функции. На входе - переменная, имеющая гамма-распределение.

set.seed(2016)
y = rgamma(500, 5)
transformedData = compareTransformations(y)

str(transformedData)

List of 16
 $ scaled      : num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ BoxCox      : num [1:500] -1.008 -1.227 0.255 0.786 1.071 ...
 $ YeoJohnson  : num [1:500] -1.019 -1.236 0.256 0.792 1.078 ...
 $ expoTrans   : num [1:500] -1.045 -1.25 0.253 0.817 1.112 ...
 $ range       : num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ knnImpute   : num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ medianImpute: num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ zv          : num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ nzv         : num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ conditionalX: num [1:500] -0.9695 -1.1105 0.0982 0.6961 1.0574 ...
 $ log         : num [1:500] -1.023 -1.242 0.261 0.798 1.081 ...
 $ sqrt        : num [1:500] -1.008 -1.187 0.182 0.758 1.084 ...
 $ 1/(data+1)  : num [1:500] 1.1121 -0.2703 0.059 0.0483 0.0449 ...
 $ data^2      : num [1:500] -0.0343 0.139 -0.5837 -0.3032 0.0709 ...
 $ logit       : num [1:500] -0.578 -0.885 0.408 0.619 0.707 ...
 $ arcsin      : num [1:500] -0.989 -1.177 0.198 0.771 1.106 ...

А вот другой пример.

y = 1/rnorm(500, 40)*rf(500,40,5)
transformedData = compareTransformations(y)