Imputação múltipla para valores ausentes

Eu gostaria de usar a imputação para substituir valores ausentes no meu conjunto de dados sob certas restrições.

Por exemplo, eu gostaria que a variável imputada x1fosse maior ou igual à soma das minhas outras duas variáveis, digamos x2e x3. Eu também quero x3ser imputado por um 0ou outro >= 14e quero x2ser imputado por um 0ou outro >= 16.

Tentei definir essas restrições no SPSS para imputação múltipla, mas no SPSS só posso definir valores máximos e mínimos. Existe alguma maneira de definir restrições adicionais no SPSS ou você conhece algum pacote R que me permita definir essas restrições para imputação de valores ausentes?

Meus dados são os seguintes:

   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

Uma solução é escrever suas próprias funções de imputação personalizadas para o micepacote. O pacote está preparado para isso e a instalação surpreendentemente indolor.

Primeiro, configuramos os dados conforme sugerido:

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))

Em seguida, carregamos o micepacote e vemos quais métodos ele escolhe por padrão:

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

A pmmsignifica correspondência média preditivo - provavelmente o algoritmo imputação populares mais para a imputação de variáveis contínuas. Ele calcula o valor previsto usando um modelo de regressão e escolhe os 5 elementos mais próximos do valor previsto (por distância euclidiana ). Esses elementos escolhidos são chamados de pool de doadores e o valor final é escolhido aleatoriamente nesse pool de doadores.

A partir da matriz de previsão, descobrimos que os métodos passam as variáveis ​​que são de interesse para as restrições. Observe que a linha é a variável de destino e a coluna os preditores. Se x1 não tivesse 1 na coluna x3, teríamos que adicionar isso na matriz:imp_base$predictorMatrix["x1","x3"] <- 1

Agora, a parte divertida, gerando os métodos de imputação. Eu escolhi um método bastante grosseiro aqui, onde descarto todos os valores se eles não atenderem aos critérios. Isso pode resultar em um longo tempo de loop e pode ser potencialmente mais eficiente manter as imputações válidas e refazer apenas as restantes, exigindo um pouco mais de ajustes.

# 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)
  }

Assim que terminamos de definir os métodos, simplesmente mudamos os métodos anteriores. Se você deseja alterar apenas uma única variável, pode simplesmente usar, imp_base$method["x2"] <- "pmm_x2"mas neste exemplo alteraremos tudo (a nomeação não é necessária):

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)

Agora, vamos dar uma olhada no terceiro conjunto de dados imputados:

> 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

Ok, isso faz o trabalho. Eu gosto dessa solução, pois você pode pegar carona nas funções principais e apenas adicionar as restrições que achar significativas.

Atualizar

Para impor as restrições rigorosas @ t0x1n mencionadas nos comentários, podemos adicionar as seguintes habilidades à função de wrapper:

1. Salve valores válidos durante os loops para que os dados de execuções anteriores, parcialmente bem-sucedidas, não sejam descartados
2. Um mecanismo de escape para evitar loops infinitos
3. Infle o pool de doadores depois de tentar x vezes sem encontrar uma correspondência adequada (isso se aplica principalmente a pmm)

Isso resulta em uma função de invólucro um pouco mais complicada:

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)
}

Observe que isso não funciona tão bem, provavelmente devido ao fato de o conjunto de dados sugerido falhar nas restrições de todos os casos sem faltar. Preciso aumentar o comprimento do loop para 400-500 antes que ele comece a se comportar. Presumo que isso não seja intencional, sua imputação deve imitar como os dados reais são gerados.

Otimização

O argumento rycontém os valores que não faltam e poderíamos acelerar o loop removendo os elementos que encontramos imputações elegíveis, mas como não estou familiarizado com as funções internas que me abstive disso.

Acho que o mais importante quando você tem fortes restrições que levam tempo para preencher é paralelizar suas imputações ( veja minha resposta no CrossValidated ). Atualmente, a maioria tem computadores com 4-8 núcleos e R usa apenas um deles por padrão. O tempo pode ser (quase) cortado ao meio, dobrando o número de núcleos.

Parâmetros ausentes na imputação

Com relação ao problema de x2estar ausente no momento da imputação - os ratos nunca alimentam os valores ausentes no x- data.frame. O método de ratos inclui o preenchimento de algum valor aleatório no início. A parte da cadeia da imputação limita o impacto desse valor inicial. Se você observar a micefunção-, poderá encontrá-lo antes da chamada de imputação (a mice:::sampler-função):

...
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))
  }
}
...

O data.initpode ser fornecido para a micefunção e o mouse.imput.sample é um procedimento básico de amostragem.

Sequência de visitas

Se a sequência de visitas for importante, você poderá especificar a ordem em que a micefunção-executa as imputações. O padrão é de, 1:ncol(data)mas você pode definir o visitSequenceque quiser.

A coisa mais próxima que pude encontrar é a inclusão prévia de informações de Amelia . Veja o capítulo 4.7 da vinheta , especificamente 4.7.2:

Priores a nível de observação

Os pesquisadores geralmente têm informações prévias adicionais sobre a falta de valores de dados com base em pesquisas anteriores, consenso acadêmico ou experiência pessoal. Amelia pode incorporar essas informações para produzir imputações muito melhoradas. O algoritmo Amelia permite que os usuários incluam prévios bayesianos informativos sobre células de dados ausentes individuais em vez dos parâmetros de modelo mais gerais, muitos dos quais com pouco significado direto.

A incorporação de anteriores segue a análise bayesiana básica, em que a imputação acaba sendo uma média ponderada da imputação baseada no modelo e a média anterior, onde os pesos são funções da força relativa dos dados e anteriores: quando o modelo prevê muito bem , a imputação diminuirá o peso do anterior e vice-versa (Honaker e King, 2010).

Os antecedentes sobre observações individuais devem descrever a crença do analista sobre a distribuição da célula de dados ausentes. Isso pode assumir a forma de média e desvio padrão ou intervalo de condicional. Por exemplo, podemos saber que as tarifas de 1986 na Tailândia estão em torno de 40%, mas temos alguma incerteza quanto ao valor exato. Nossa crença anterior sobre a distribuição da célula de dados ausentes, então, centra-se em 40 com um desvio padrão que reflete a quantidade de incerteza que temos sobre nossa crença anterior.

Para inserir anteriores, você deve criar uma matriz de anteriores com quatro ou cinco colunas. Cada linha da matriz representa um prior em uma observação ou em uma variável. Em qualquer linha, a entrada na primeira coluna é a linha da observação e a entrada é a segunda coluna é a coluna da observação. Na matriz de quatro colunas anteriores, a terceira e quarta colunas são a média e o desvio padrão da distribuição anterior do valor ausente.

Portanto, embora você geralmente não possa dizer algo assim x1<x2+x3, pode fazer um loop sobre o conjunto de dados e adicionar um nível de observação antes para cada caso relevante. Limites constantes também podem ser aplicados (como definir x1, x2 e x3 para não serem negativos). Por exemplo:

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);

Provavelmente, é mais fácil implementar restrições na média preditiva correspondente à imputação múltipla. Isso pressupõe que haja um número significativo de observações com variáveis ​​de restrição não ausentes que atendam às restrições. Estou pensando em implementar isso na função de Hmiscpacote R. aregImputeVocê pode verificar novamente em mais ou menos um mês. Será importante especificar a distância máxima do alvo que pode ser uma observação do doador, porque as restrições empurrarão os doadores para mais longe do doador ideal sem restrições.

Acredito que o Ameliapacote (Amelia II) atualmente tenha o suporte mais abrangente para especificar restrições de intervalo de valores de dados. No entanto, o problema é que Ameliaassume que os dados são multivariados normais.

Se, no seu caso, a suposição de normalidade multivariada não se aplicar, convém verificar o micepacote, que implementa a imputação múltipla (MI) por meio de equações em cadeia . Este pacote não tem a suposição de normalidade multivariada . Ele também tem uma função que pode ser suficiente para especificar restrições , mas não tenho certeza em que grau. A função é chamada squeeze(). Você pode ler sobre isso na documentação: http://cran.r-project.org/web/packages/mice/mice.pdf . Um benefício adicional de http://www.ats.ucla.edu/stat/r/faq/R_pmm_mi.htm .mice é sua flexibilidade em termos de permitir a especificação de funções de imputação definidas pelo usuário e uma seleção mais ampla de algoritmos. Aqui está um tutorial sobre como executar MI, usando mice:

Até onde eu entendo, o Hmiscpacote do Dr. Harrell , usando as mesmas abordagens de equações encadeadas ( correspondência preditiva média ), provavelmente suporta dados não normais (com exceção do normpmmmétodo). Talvez ele já tenha implementado a funcionalidade de especificação de restrições conforme sua resposta acima. Eu não usei aregImpute(), então não posso dizer muito mais sobre isso (usei Ameliae mice, mas definitivamente não sou especialista em estatística, apenas tentando aprender o máximo que posso).

Por fim, você pode achar interessante o seguinte, um pouco datado, mas ainda agradável, visão geral de abordagens, métodos e software para imputação múltipla de dados com valores ausentes: http://www.ncbi.nlm.nih.gov/pmc/articles / PMC1839993 . Tenho certeza de que existem documentos de visão geral mais recentes sobre MI, mas é disso que estou ciente no momento. Espero que isso seja um pouco útil.

Se entendi sua pergunta corretamente, parece-me que você já sabe quais valores as variáveis ​​ausentes devem estar sujeitas a algumas restrições. Não sou muito familiarizado com o SPSS, mas no RI acho que você pode escrever uma função para fazer isso (o que não deve ser muito difícil, dependendo da sua experiência, devo dizer). Não conheço nenhum pacote que funcione com essas restrições.