Imputação múltipla para valores ausentes

13

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
rosa
fonte
Eu mudei 0 or 16 or >= 16para 0 or >= 16desde >=16inclui o valor 16. Espero que isso não atrapalhe o seu significado. Mesmo para0 or 14 or >= 14
Alexis

Respostas:

16

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.

Max Gordon
fonte
+1 Isso é ótimo, exatamente o que eu tinha em mente (veja meu comentário na resposta de Frank) e, com certeza, o candidato número 1 à recompensa a partir de agora. No pmm_x1entanto, algumas coisas me incomodam : (1) Tomar a soma máxima de qualquer combinação possível de x2e x3de todo o conjunto de dados é muito mais difícil do que a restrição original. O correto seria testar que para cada linha , x1 < x2 + x3. É claro que quanto mais linhas você tiver, menor será sua chance de cumprir essa restrição (como uma única linha ruim estraga tudo) e mais tempo o loop pode ficar.
t0x1n
(2) se ambos x1e x2estiverem ausentes, você pode atribuir um valor para o x1qual as restrições são mantidas (digamos 50), mas uma vez x2imputadas, elas são quebradas (digamos que sejam imputadas 55). Existe uma maneira de imputar "horizontalmente" e não verticalmente? Dessa forma, poderíamos imputar uma única linha de x1, x2e, x3e simplesmente re-impute-lo até que linha específica cai sob as restrições. Isso deve ser rápido o suficiente e, uma vez feito isso, podemos passar para a próxima linha. Obviamente, se o MI é "vertical" por natureza, estamos sem sorte. Nesse caso, talvez a abordagem que Aleksandr tenha mencionado?
t0x1n
Solução legal, +1! Pode ser especialmente útil, pois atualmente usomice pacote. Obrigado por compartilhar.
Aleksandr Blekh 25/11
1
@ t0x1n Atualizei minha resposta com uma função de invólucro mais avançada, de acordo com seus comentários. Se você quiser mergulhar mais fundo, recomendo que você brinque com o debug()para ver como mice.impute.pmme seus irmãos funcionam sob o capô.
Max Gordon
1
@ t0x1n: Eu acho - inspecione seus valores imputados. Se eles parecerem irrealistas, você poderá escolher minha abordagem para atribuir apenas aqueles que são menos centrais ao modelo. No meu caso, optei por excluir aqueles sem radiografias de acompanhamento, pois eles estão no centro do estudo e as imputações não fornecem valores clinicamente plausíveis (a perna fica mais longa após uma fratura). Não estou totalmente satisfeito com isso, mas parece um compromisso razoável.
Max Gordon
8

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);
t0x1n
fonte
5

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.

Frank Harrell
fonte
Eu adoraria ter isso também. Eu só preciso das restrições inter-variáveis ​​mais básicas, digamos x<y<z.
T1x1n
Perdoe minha ignorância se estou longe, mas fiquei com a impressão de que várias técnicas de imputação envolvem o desenho de valores de uma distribuição apropriada. Não deveria ser uma questão simples usar amostragem por rejeição? por exemplo, continue desenhando até que alguma restrição especificada (como x1<x2) seja atendida?
t0x1n
Isso é o que eu poderia fazer com a aregImputefunção R com correspondência preditiva média. Mas e se nenhuma das observações do doador (quase correspondências de previsões) satisfizer as restrições para a observação do alvo ser imputada, mesmo que elas obviamente tenham que atender às restrições do conjunto de variáveis ​​do doador?
Frank Harrell
Nesse caso, talvez leve diretamente o valor previsto? Isso depende apenas da regressão (sem a fase PMM) para essa amostra?
t0x1n
A imputação de regressão é um pouco mais provável de apresentar valores imputados que são inconsistentes com o restante do registro do sujeito. Portanto, não acho que esse seja um motivo para evitar o PMM.
26414 Frank #: 22340 #:
4

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.

Aleksandr Blekh
fonte
1
Esse bom comentário me faz pensar que a correspondência média preditiva, que substitui as perdas por valores realmente observados, já pode incorporar alguns tipos de restrições se todos os dados observados atenderem a essas restrições. Eu apreciaria alguém pensando nisso. Ainda não implementei nenhuma restrição especial no aregImpute.
Frank Harrell
1
Você está certo. Acabei de perceber que as observações dos doadores fornecem valores consistentes com suas outras variáveis, mas não com as outras variáveis ​​na variável de destino.
Frank Harrell
1
Além das premissas de distribuição feitas por Amelia, você foi capaz de especificar restrições em mais detalhes do que demonstrei em minha resposta? O problema squeezeé que seus limites são constantes, então você não pode especificar nada parecido x1<x2. Além disso, parece ter sido invocado no vetor de resultado imputado, que acredito ser tarde demais. Parece-me que os limites devem ser considerados durante o processo de imputação, para que tenham mais significado do que um ajuste posterior ao fato.
t0x1n
1
@ t0x1n: Infelizmente, não tive a chance de especificar restrições Amelia, porque mudei para mice, assim que meus testes confirmaram que meus dados não são normais multivariados. No entanto, recentemente deparei com este conjunto muito agradável de slides de apresentação sobre o tópico (métodos e software de MI): statistik.lmu.de/~fkreuter/imputation_sose2011/downloads/… . Se entendi corretamente, descreve uma solução em potencial para o problema de restrições (consulte a página 50 do PDF - e não o slide número 50!). Espero que isto ajude.
Aleksandr Blekh
1
@ t0x1n: Na verdade, a solução é descrito nas páginas 50 e 51.
Aleksandr Blekh
0

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.

ThinkStatsme
fonte