Множественное вменение для пропущенных значений

Я хотел бы использовать вменение для замены отсутствующих значений в моем наборе данных при определенных ограничениях.

Например, я бы хотел, чтобы вмененная переменная x1была больше или равна сумме двух других моих переменных, скажем, x2и x3. Я также хочу x3быть вмененным либо 0или, >= 14и я хочу x2быть вмененным либо 0или >= 16.

Я попытался определить эти ограничения в SPSS для множественного вменения, но в SPSS я могу определить только максимальные и минимальные значения. Есть ли способ определить дополнительные ограничения в SPSS или вы знаете какой-нибудь пакет R, который позволил бы мне определить такие ограничения для вменения пропущенных значений?

Мои данные таковы:

   x1 =c(21, 50, 31, 15, 36, 82, 14, 14, 19, 18, 16, 36, 583, NA,NA,NA, 50, 52, 26, 24)
   x2 = c(0, NA, 18,0, 19, 0, NA, 0, 0, 0, 0, 0, 0,NA,NA, NA, 22, NA, 0, 0)
   x3 = c(0, 0, 0, 0, 0, 54, 0 ,0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, 0)
   dat=data.frame(x1=x1, x2=x2, x3=x3)
   > dat
       x1 x2 x3
   1   21  0  0
   2   50 NA  0
   3   31 18  0
   4   15  0  0
   5   36 19  0
   6   82  0 54
   7   14 NA  0
   8   14  0  0
   9   19  0  0
   10  18  0  0
   11  16  0  0
   12  36  0  0
   13 583  0  0
   14  NA NA NA
   15  NA NA NA
   16  NA NA NA
   17  50 22 NA
   18  52 NA  0
   19  26  0  0
   20  24  0  0

Одним из решений является написание ваших собственных функций вменения для miceпакета. Пакет подготовлен для этого, и установка на удивление безболезненна.

Сначала мы настроим данные как предложено:

dat=data.frame(x1=c(21, 50, 31, 15, 36, 82, 14, 14, 19, 18, 16, 36, 583, NA,NA,NA, 50, 52, 26, 24), 
               x2=c(0, NA, 18,0, 19, 0, NA, 0, 0, 0, 0, 0, 0,NA,NA, NA, 22, NA, 0, 0), 
               x3=c(0, 0, 0, 0, 0, 54, 0 ,0, 0, 0, 0, 0, 0, NA, NA, NA, NA, 0, 0, 0))

Далее мы загружаем miceпакет и видим, какие методы он выбирает по умолчанию:

library(mice)
# Do a non-imputation
imp_base <- mice(dat, m=0, maxit = 0)

# Find the methods that mice chooses
imp_base$method
# Returns: "pmm" "pmm" "pmm"

# Look at the imputation matrix
imp_base$predictorMatrix
# Returns:
#   x1 x2 x3
#x1  0  1  1
#x2  1  0  1
#x3  1  1  0

pmmОзначает предсказательной среднего соответствия - вероятно, самый популярный алгоритм вменения для вменения непрерывных переменных. Он рассчитывает прогнозируемое значение с использованием регрессионной модели и выбирает 5 ближайших элементов к прогнозируемому значению (по евклидову расстоянию ). Эти выбранные элементы называются донорским пулом, и окончательное значение выбирается случайным образом из этого донорского пула.

Из матрицы прогнозирования мы находим, что методы получают переданные переменные, которые представляют интерес для ограничений. Обратите внимание, что строка является целевой переменной, а столбец - предикторами. Если x1 не имеет 1 в столбце x3, мы должны добавить это в матрицу:imp_base$predictorMatrix["x1","x3"] <- 1

Теперь самое интересное, генерация методов вменения. Я выбрал довольно грубый метод, где я отбрасываю все значения, если они не соответствуют критериям. Это может привести к продолжительному циклу, и, возможно, будет более эффективно сохранить действительные вменения и только повторить оставшиеся, однако потребуется немного больше настроек.

# Generate our custom methods
mice.impute.pmm_x1 <- 
  function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "", 
            ...) 
  {
    max_sum <- sum(max(x[,"x2"], na.rm=TRUE),
                   max(x[,"x3"], na.rm=TRUE))
    repeat{
      vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
                              version = "", ...)
      if (all(vals < max_sum)){
        break
      }
    }
    return(vals)
  }

mice.impute.pmm_x2 <- 
  function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "", 
            ...) 
  {
    repeat{
      vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
                              version = "", ...)
      if (all(vals == 0 | vals >= 14)){
        break
      }
    }
    return(vals)
  }

mice.impute.pmm_x3 <- 
  function (y, ry, x, donors = 5, type = 1, ridge = 1e-05, version = "", 
            ...) 
  {
    repeat{
      vals <- mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
                              version = "", ...)
      if (all(vals == 0 | vals >= 16)){
        break
      }
    }
    return(vals)
  }

Как только мы закончили с определением методов, мы просто изменили предыдущие методы. Если вы хотите изменить только одну переменную, вы можете просто использовать, imp_base$method["x2"] <- "pmm_x2"но для этого примера мы изменим все (наименование не обязательно):

imp_base$method <- c(x1 = "pmm_x1", x2 = "pmm_x2", x3 = "pmm_x3")

# The predictor matrix is not really necessary for this example
# but I use it just to illustrate in case you would like to 
# modify it
imp_ds <- 
  mice(dat, 
       method = imp_base$method, 
       predictorMatrix = imp_base$predictorMatrix)

Теперь давайте посмотрим на третий вмененный набор данных:

> complete(imp_ds, action = 3)
    x1 x2 x3
1   21  0  0
2   50 19  0
3   31 18  0
4   15  0  0
5   36 19  0
6   82  0 54
7   14  0  0
8   14  0  0
9   19  0  0
10  18  0  0
11  16  0  0
12  36  0  0
13 583  0  0
14  50 22  0
15  52 19  0
16  14  0  0
17  50 22  0
18  52  0  0
19  26  0  0
20  24  0  0

Хорошо, это делает работу. Мне нравится это решение, так как вы можете использовать верхнюю часть основных функций и просто добавить ограничения, которые вы считаете значимыми.

Обновить

Чтобы применить строгие ограничения @ t0x1n, упомянутые в комментариях, мы можем захотеть добавить следующие функции в функцию-обертку:

1. Сохраняйте действительные значения во время циклов, чтобы данные из предыдущих, частично успешных прогонов не отбрасывались
2. Механизм побега, чтобы избежать бесконечных петель
3. Надуйте пул доноров после попытки x раз, не найдя подходящего соответствия (это в первую очередь относится к pmm)

Это приводит к несколько более сложной функции-обертке:

mice.impute.pmm_x1_adv <-   function (y, ry, 
                                      x, donors = 5, 
                                      type = 1, ridge = 1e-05, 
                                      version = "", ...) {
  # The mice:::remove.lindep may remove the parts required for
  # the test - in those cases we should escape the test
  if (!all(c("x2", "x3") %in% colnames(x))){
    warning("Could not enforce pmm_x1 due to missing column(s):",
            c("x2", "x3")[!c("x2", "x3") %in% colnames(x)])
    return(mice.impute.pmm(y, ry, x, donors = 5, type = 1, ridge = 1e-05,
                           version = "", ...))
  }

  # Select those missing
  max_vals <- rowSums(x[!ry, c("x2", "x3")])

  # We will keep saving the valid values in the valid_vals
  valid_vals <- rep(NA, length.out = sum(!ry))
  # We need a counter in order to avoid an eternal loop
  # and for inflating the donor pool if no match is found
  cntr <- 0
  repeat{
    # We should be prepared to increase the donor pool, otherwise
    # the criteria may become imposs
    donor_inflation <- floor(cntr/10)
    vals <- mice.impute.pmm(y, ry, x, 
                            donors = min(5 + donor_inflation, sum(ry)), 
                            type = 1, ridge = 1e-05,
                            version = "", ...)

    # Our criteria check
    correct <- vals < max_vals
    if (all(!is.na(valid_vals) |
              correct)){
      valid_vals[correct] <-
        vals[correct]
      break
    }else if (any(is.na(valid_vals) &
                    correct)){
      # Save the new valid values
      valid_vals[correct] <-
        vals[correct]
    }

    # An emergency exit to avoid endless loop
    cntr <- cntr + 1
    if (cntr > 200){
      warning("Could not completely enforce constraints for ",
              sum(is.na(valid_vals)),
              " out of ",
              length(valid_vals),
              " missing elements")
      if (all(is.na(valid_vals))){
        valid_vals <- vals
      }else{
        valid_vals[is.na(valid_vals)] <- 
          vals[is.na(valid_vals)]
      }
      break
    }
  }
  return(valid_vals)
}

Обратите внимание, что это не очень хорошо работает, скорее всего из-за того, что предлагаемый набор данных не пропускает ограничения для всех случаев. Мне нужно увеличить длину цикла до 400-500, прежде чем он начнет себя вести. Я предполагаю, что это непреднамеренно, ваше вменение должно имитировать, как генерируются фактические данные.

оптимизация

Аргумент ryсодержит не пропущенные значения, и мы могли бы ускорить цикл, удалив элементы, которые мы нашли приемлемые вменения, но, поскольку я не знаком с внутренними функциями, я воздержался от этого.

Я думаю, что наиболее важная вещь, когда у вас есть серьезные ограничения, требующие времени для полного заполнения, - это распараллелить ваши вменения ( см. Мой ответ на CrossValidated ). Большинство из них сегодня имеют компьютеры с 4-8 ядрами, и R по умолчанию использует только одно из них. Время можно (почти) разрезать пополам, удвоив количество ядер.

Отсутствующие параметры при вменении

Что касается проблемы x2отсутствия во время вменения - мыши фактически никогда не подают пропущенные значения в x- data.frame. Метод мышей включает заполнение некоторого случайного значения при запуске. Цепная часть вменения ограничивает влияние этого начального значения. Если вы посмотрите на mice-функцию, вы можете найти ее до вызова вменения ( mice:::sampler-функция):

...
if (method[j] != "") {
  for (i in 1:m) {
    if (nmis[j] < nrow(data)) {
      if (is.null(data.init)) {
        imp[[j]][, i] <- mice.impute.sample(y, 
                                            ry, ...)
      }
      else {
        imp[[j]][, i] <- data.init[!ry, j]
      }
    }
    else imp[[j]][, i] <- rnorm(nrow(data))
  }
}
...

Функция data.initможет быть передана в miceфункцию, и мыши. Imput.sample является основной процедурой выборки.

Последовательность посещений

Если последовательность посещений важна, вы можете указать порядок, в котором mice-функция запускает вычисления. По умолчанию используется значение from, 1:ncol(data)но вы можете установить visitSequenceвсе, что вам нравится.

Самая близкая вещь, которую я мог найти, является предварительным информационным включением Амелии . Смотрите главу 4.7 в виньетке , в частности 4.7.2:

Приоры уровня наблюдения

Исследователи часто имеют дополнительную предварительную информацию о пропущенных значениях данных, основанную на предыдущих исследованиях, академическом консенсусе или личном опыте. Амелия может использовать эту информацию, чтобы получить значительно улучшенную оценку. Алгоритм Амелия позволяет пользователям включать информативные байесовские априоры об отдельных пропущенных ячейках данных вместо более общих параметров модели, многие из которых имеют небольшое прямое значение.

Включение априорных значений следует базовому байесовскому анализу, где вменение оказывается средневзвешенным значением вменения на основе модели и предшествующим средним значением, где веса являются функциями относительной силы данных и априорного значения: когда модель очень хорошо предсказывает , вменение снизит вес предыдущего и наоборот (Honaker and King, 2010).

Приоритеты отдельных наблюдений должны описывать мнение аналитика о распределении недостающей ячейки данных. Это может принимать форму среднего значения и стандартного отклонения или интервал соответствия. Например, мы можем знать, что в 1986 году тарифы в Таиланде составляют около 40%, но у нас есть некоторая неопределенность относительно точного значения. Таким образом, наше предыдущее мнение о распределении отсутствующей ячейки данных основывается на 40 со стандартным отклонением, которое отражает степень неопределенности, которую мы имеем относительно нашего предыдущего убеждения.

Для ввода априоров вы должны построить матрицу априоров с четырьмя или пятью столбцами. Каждая строка матрицы представляет априор либо по одному наблюдению, либо по одной переменной. В любой строке запись в первом столбце является строкой наблюдения, а запись - вторым столбцом столбца наблюдения. В матрице с четырьмя столбцами третьего и четвертого столбцов представлены среднее значение и стандартное отклонение предыдущего распределения отсутствующего значения.

Таким образом, хотя вы не сможете вообще сказать что-то подобное x1<x2+x3, вы можете зациклить свой набор данных и добавить уровень наблюдения до каждого соответствующего случая. Также могут быть применены постоянные границы (например, установка x1, x2 и x3 как неотрицательные). Например:

priors = matrix(NA, nrow=0, ncol=5);
for (i in seq(1, length(data))) 
{
    x1 = data$x1[i];
    x2 = data$x2[i];
    x3 = data$x3[i];

    if (is.na(x1) && !is.na(x2) && !is.na(x3))
    {
        priors = rbind(priors, c(i, 1, 0, x2+x3, 0.999999))
    }
}

amelia(data, m=1, bound = rbind(c(1, 0, Inf), c(2, 0, Inf), c(3, 0, Inf)), pr = priors);

Ограничения, вероятно, легче реализовать в предсказывающем среднем совпадении множественного вменения. Это предполагает, что существует значительное количество наблюдений с не отсутствующими ограничивающими переменными, которые удовлетворяют ограничениям. Я думаю о реализации этого в функции Hmiscпакета R. aregImputeВы можете проверить через месяц или около того. Будет важно указать максимальное расстояние от цели, которым может быть наблюдение донора, потому что ограничения будут отталкивать доноров от идеального неограниченного донора.

Я считаю, что Ameliaпакет (Amelia II) в настоящее время имеет наиболее полную поддержку для определения ограничений диапазона значений данных. Однако проблема состоит в том, что Ameliaпредполагается, что данные многомерны нормально.

Если в вашем случае предположение о многомерной нормальности не применимо, вы можете проверить miceпакет, который реализует множественное вменение (MI) через цепочечные уравнения . Этот пакет не предполагает многомерной нормальности . Он также имеет функцию, которая может быть достаточно для определения ограничений , но я не уверен, в какой степени. Функция называется squeeze(). Вы можете прочитать об этом в документации: http://cran.r-project.org/web/packages/mice/mice.pdf . Дополнительным преимуществом miceявляется его гибкость с точки зрения возможности задания пользовательских функций вменения и более широкого выбора алгоритмов. Вот учебник по выполнению MI, используя mice:http://www.ats.ucla.edu/stat/r/faq/R_pmm_mi.htm .

Насколько я понимаю, Hmiscпакет доктора Харрелла , использующий тот же подход цепных уравнений ( прогнозирующее совпадение среднего ), вероятно, поддерживает ненормальные данные (за исключением normpmmметода). Возможно, он уже реализовал функциональность спецификации ограничений согласно ответу выше. Я не использовал aregImpute(), поэтому не могу сказать больше об этом (я использовал Ameliaи mice, но я определенно не эксперт в области статистики, просто пытаюсь узнать как можно больше).

Наконец, вам может быть интересен следующий, немного устаревший, но все же хороший обзор подходов, методов и программного обеспечения для многократного ввода данных с пропущенными значениями: http://www.ncbi.nlm.nih.gov/pmc/articles / PMC1839993 . Я уверен, что есть новые обзорные статьи по МИ, но это все, что я знаю в настоящее время. Я надеюсь, что это несколько полезно.

Если я правильно понимаю ваш вопрос, мне кажется, что вы уже знаете, какие значения должны принимать отсутствующие переменные с учетом некоторых ограничений. Я не очень разбираюсь в SPSS, но в RI думаю, что вы можете написать функцию для этого (что не должно быть слишком сложным, в зависимости от вашего опыта, я должен сказать). Я не знаю ни одного пакета, который работает с такими ограничениями.