Imputação
Introdução
- Introduzir o conceito de imputação de dados;
- Explicar os tipos de dados missing.
Pacote MICE (multiple imputation by chained equations)
- Referência dos dados e metodologia na pagina 9
- Outro estudo com estes dados & imputação
- Artigo: Van Buuren, S., Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Software, 45(3), 1-67. https://www.jstatsoft.org/v45/i03/
Nesta abordagem serão ajustados 3 modelos de regressão linear com efeitos mistos
- Modelo 1: modelo com dados completos;
- Modelo 2: modelo com dados missing ou faltantes;
- Modelo 3: modelo com dados imputados. Serão tiradas conclusões a partir destes ajustes, mostrando as vantagens do método.
Pacotes
Carregar os pacotes - verifique a necessidade de instalar pacotes novos - tirando o comentário (#) para instalar pela primeira vez
require(DT) #para a função datatable
require(tidyverse)
#devtools::install_github("tidymodels/broom")
#install.packages("mice")
require(mice)
#install.packages("broom")
#install.packages("broom.mixed")
require(broom.mixed)
require(lme4)1 Utilização dos dados do pacote MICE
- Pegar os dados de potthoffroy do pacote - os dados consistem das seguintes informações:
- id: identificação do indivíduo - qualitativa codificada em números;
- sex: sexo do indivíduo - qualitativa.
- d8, d10, d12 e d14: distâncias (mm) do centro da glândula pituitária à fissura pteromaxilar, nas idades de 8, 10, 12 e 14 anos respectivamente.
data(potthoffroy)
datatable(potthoffroy,options = list(pageLength = 5))Transformando em long table:
- foi criada uma função para pivotar pois faremos esta mesma tarefa três vezes!
pivotar_long = function(dataset_in){
dataset_in %>%
`colnames<-`(c("id","sex","8","10","12","14")) %>%
pivot_longer(!id:sex, names_to = "age", values_to = "distance") %>%
mutate(age = as.numeric(age))
}
dados_long=pivotar_long(potthoffroy)
datatable(dados_long,cap="Tabela 1: potthoffroy no formato long table",options = list(pageLength = 5))Visualização gráfica
- Boxplot
- o boxplot apresenta uma visão geral dos dados, sem interpretar o efeito dos indivíduos;
- em geral, os indivíduos do sexo masculino apresentam valores superiores de distância.
dados_long %>%
ggplot(aes(x=age,y=distance,fill=sex)) +
geom_boxplot()+
facet_wrap(~sex)
- Gráfico de perfis
- o gráfico de perfis são típicos de uma análise com dados longitudinais;
- temos os perfis dos indivíduos ao longo das idades (age), separados por sexo;
- os perfis exibem comportamento linear;
- Em geral os perfis dos homens são superiores comparados aos das mulheres;
- os perfis dos homens parecem ter intercepto (início) e inclinação superior aos das mulheres;
ggplot(dados_long,aes (x = age, y = distance, color = as.factor(id))) +
geom_point() +
geom_line() +
facet_wrap(~sex)+
theme(legend.position = "none") 
- Gráfico de perfis por indivíduo - para enxergar a idéia do efeito aleatório;
- observamos que os indivíduos exibem perfis com comportamento linear;
- alguns perfis diferem tanto no intercepto quanto na inclinação da reta;
- em geral os perfis do sexo F diferem mais com relação a inclinação da reta.
dados_long %>%
filter(sex=="F") %>%
ggplot(aes(x=age,y=distance))+
geom_point()+
geom_smooth(method="lm")+
facet_wrap(~id)+
ggtitle("Perfis de distância para o sexo F")## `geom_smooth()` using formula 'y ~ x'
dados_long %>%
filter(sex=="M") %>%
ggplot(aes(x=age,y=distance))+
geom_point()+
geom_smooth(method="lm")+
facet_wrap(~id)+
ggtitle("Perfis de distância para o sexo M") ## `geom_smooth()` using formula 'y ~ x'
- Gráfico de interação:
- pelo gráfico de interação parece haver interação entre as variáveis sexo e idade - a partir dos 10 anos os segmentos de reta não são aproximadamente paralelos;
- a idade parece interferir de forma distinta nas distâncias, de acordo com o sexo.
interaction.plot(x.factor = dados_long$age,
trace.factor = dados_long$sex,
response = dados_long$distance,
fun = mean,
type="b",
col=c("black","red"), ### Colors for levels of trace var.
pch=c(19, 17), ### Symbols for levels of trace var.
fixed=TRUE, ### Order by factor order in data
leg.bty = "o",
ylab="distance",
xlab="age")
Modelo 1: Ajustando um modelo linear misto com os dados completos
- o modelo escolhido é o mais simples;
- o efeito aleatório é no intercepto;
- o modelo não considera interação entre sexo e idade;
modelo_1 = lmer(distance ~ age + sex + (1 | id), data = dados_long)
modelo_1## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + sex + (1 | id)
## Data: dados_long
## REML criterion at convergence: 437.5125
## Random effects:
## Groups Name Std.Dev.
## id (Intercept) 1.807
## Residual 1.432
## Number of obs: 108, groups: id, 27
## Fixed Effects:
## (Intercept) age sexM
## 15.3857 0.6602 2.3210broom::tidy(modelo_1, conf.int = TRUE)## # A tibble: 5 x 8
## effect group term estimate std.error statistic conf.low conf.high
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 15.4 0.896 17.2 13.6 17.1
## 2 fixed <NA> age 0.660 0.0616 10.7 0.539 0.781
## 3 fixed <NA> sexM 2.32 0.761 3.05 0.829 3.81
## 4 ran_pars id sd__(Interce~ 1.81 NA NA NA NA
## 5 ran_pars Residu~ sd__Observat~ 1.43 NA NA NA NATransformando alguns dados em “missing data” ou dados faltantes, de forma aleatória
- Utilizando uma semente unica - set.seed
- Colocando os dados missing (NA) de forma aleatória, sendo \(10\%\) dos dados
set.seed(11112020)
dados=potthoffroy
n = nrow(dados)
nmiss=round(0.1*n)
c=2
while(c<=6){
idmis = sample(seq(1,n,1),nmiss,replace=FALSE)
dados[idmis, c] <- NA #na coluna c vai entrar o NA
c=c+1
}
datatable(dados,cap="Tabela 2: potthoffroy com missing",options = list(pageLength = 5))Identificando o padrão de missing
- Auxílio no material do MICE em bookdown material do autor do pacote
- Segundo o autor, é apropriado utilizar o formato wide para identificar os padrões (formato em que as idades dos indivíduos aparecem em várias colunas);
- PENDENTE de analise.
md.pattern(dados)
## id sex d8 d10 d12 d14
## 15 1 1 1 1 1 1 0
## 2 1 1 1 1 1 0 1
## 2 1 1 1 1 0 1 1
## 2 1 1 1 0 1 1 1
## 1 1 1 1 0 0 1 2
## 2 1 1 0 1 1 1 1
## 1 1 0 1 1 1 1 1
## 1 1 0 1 1 1 0 2
## 1 1 0 0 1 1 1 2
## 0 3 3 3 3 3 15p=md.pairs(dados)
p$mr/(p$mr+p$mm)## id sex d8 d10 d12 d14
## id NaN NaN NaN NaN NaN NaN
## sex 1 0.0000000 0.6666667 1.0000000 1.0000000 0.6666667
## d8 1 0.6666667 0.0000000 1.0000000 1.0000000 1.0000000
## d10 1 1.0000000 1.0000000 0.0000000 0.6666667 1.0000000
## d12 1 1.0000000 1.0000000 0.6666667 0.0000000 1.0000000
## d14 1 0.6666667 1.0000000 1.0000000 1.0000000 0.0000000flux(dados)[,1:3]## pobs influx outflux
## id 1.0000000 0.00000000 1.0000000
## sex 0.8888889 0.08843537 0.6666667
## d8 0.8888889 0.09523810 0.7333333
## d10 0.8888889 0.09523810 0.7333333
## d12 0.8888889 0.09523810 0.7333333
## d14 0.8888889 0.09523810 0.7333333Modelo 2: Ajustando um modelo linear misto com os dados faltantes
dados_long_miss = pivotar_long(dados)
modelo_2 = lmer(distance ~ age + sex + (1 | id), data = dados_long_miss)
modelo_2## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + sex + (1 | id)
## Data: dados_long_miss
## REML criterion at convergence: 362.1519
## Random effects:
## Groups Name Std.Dev.
## id (Intercept) 1.901
## Residual 1.547
## Number of obs: 86, groups: id, 24
## Fixed Effects:
## (Intercept) age sexM
## 15.3689 0.6721 2.2762broom::tidy(modelo_2, conf.int = TRUE)## # A tibble: 5 x 8
## effect group term estimate std.error statistic conf.low conf.high
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed <NA> (Intercept) 15.4 1.05 14.6 13.3 17.4
## 2 fixed <NA> age 0.672 0.0749 8.98 0.525 0.819
## 3 fixed <NA> sexM 2.28 0.857 2.66 0.596 3.96
## 4 ran_pars id sd__(Interce~ 1.90 NA NA NA NA
## 5 ran_pars Residu~ sd__Observat~ 1.55 NA NA NA NAImputação pelo pacote mice
- Saiba mais
- a imputação é com os dados na presença de missing e no formato wide (largo) antes de pivotar.
- exibindo o primeiro banco de dados imputados: imp_comp[[1]]
N=100
imp <- mice(dados, m = N, print = FALSE)
imp## Class: mids
## Number of multiple imputations: 100
## Imputation methods:
## id sex d8 d10 d12 d14
## "" "logreg" "pmm" "pmm" "pmm" "pmm"
## PredictorMatrix:
## id sex d8 d10 d12 d14
## id 0 1 1 1 1 1
## sex 1 0 1 1 1 1
## d8 1 1 0 1 1 1
## d10 1 1 1 0 1 1
## d12 1 1 1 1 0 1
## d14 1 1 1 1 1 0## base R approach using lapply
imp_comp <- mice::complete(imp, "all")
imp_comp[[1]]## id sex d8 d10 d12 d14
## 1 1 F 21.5 20.0 21.5 23.0
## 2 2 F 21.0 21.5 24.0 25.5
## 3 3 F 20.5 24.0 24.5 26.0
## 4 4 F 23.5 23.0 25.0 26.5
## 5 5 F 21.5 23.0 22.5 23.5
## 6 6 F 20.0 21.0 21.0 22.5
## 7 7 F 21.5 22.5 23.0 25.0
## 8 8 F 23.0 23.0 23.5 24.0
## 9 9 F 20.0 21.0 22.0 21.5
## 10 10 F 16.5 19.0 19.0 19.5
## 11 11 F 24.5 20.5 28.0 28.0
## 12 12 M 26.0 25.0 25.5 31.0
## 13 13 M 21.5 22.5 23.0 26.5
## 14 14 M 23.0 22.5 24.0 27.5
## 15 15 M 25.5 27.5 26.5 27.0
## 16 16 M 20.0 23.5 22.5 25.5
## 17 17 M 24.0 25.5 27.0 28.5
## 18 18 M 22.0 22.5 24.5 26.5
## 19 19 M 24.0 21.5 24.0 25.5
## 20 20 M 23.0 20.5 31.0 29.5
## 21 21 M 27.5 28.0 31.0 31.5
## 22 22 M 23.0 23.0 23.5 25.0
## 23 23 M 21.5 23.5 24.0 28.0
## 24 24 M 17.0 24.5 26.0 29.5
## 25 25 M 22.5 25.5 25.5 26.5
## 26 26 M 23.0 24.5 26.0 30.0
## 27 27 M 21.5 21.5 23.5 25.0Modelo 3: Aplicar o modelo misto em todos os dados imputados
- imp_comp é uma lista (list do R) de N bancos de dados imputados;
- foi utilizada a função lapply para aplicar pivotar e estimar o modelo;
- exibindo os 2 primeiros conjuntos de dados imputados (no formato long)
imp_comp2 = imp_comp %>%
lapply(pivotar_long)
imp_comp2[[1]]## # A tibble: 108 x 4
## id sex age distance
## <int> <fct> <dbl> <dbl>
## 1 1 F 8 21.5
## 2 1 F 10 20
## 3 1 F 12 21.5
## 4 1 F 14 23
## 5 2 F 8 21
## 6 2 F 10 21.5
## 7 2 F 12 24
## 8 2 F 14 25.5
## 9 3 F 8 20.5
## 10 3 F 10 24
## # ... with 98 more rowsimp_comp2[[2]]## # A tibble: 108 x 4
## id sex age distance
## <int> <fct> <dbl> <dbl>
## 1 1 M 8 22
## 2 1 M 10 20
## 3 1 M 12 21.5
## 4 1 M 14 23
## 5 2 F 8 21
## 6 2 F 10 21.5
## 7 2 F 12 24
## 8 2 F 14 25.5
## 9 3 F 8 20.5
## 10 3 F 10 24
## # ... with 98 more rowsAgora com os dados em formato long table vamos aplicar o modelo linear com efeitos mistos, da mesma forma que os modelos 1 & 2
modelo_3 <- lapply(imp_comp2, lmer, formula = distance ~ age + sex + (1 | id))E com a função pool famos consolidar os resultados obtidos nos 100 modelos ajustados - e fazer a comparação com o primeiro modelo (missing data)
- A função pool combina as estimativas provenientes dos vários conjuntos de dados imputados, utilizando as Rubin’s rules (Rubin, 1987, p. 76);
- Saiba mais sobre pool e as Rubin’s rules (Rubin, 1987, p. 76) ;
summary(modelo_1)## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + sex + (1 | id)
## Data: dados_long
##
## REML criterion at convergence: 437.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7489 -0.5503 -0.0252 0.4534 3.6575
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 3.267 1.807
## Residual 2.049 1.432
## Number of obs: 108, groups: id, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 15.38569 0.89598 17.172
## age 0.66019 0.06161 10.716
## sexM 2.32102 0.76142 3.048
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.756
## sexM -0.504 0.000summary(modelo_2)## Linear mixed model fit by REML ['lmerMod']
## Formula: distance ~ age + sex + (1 | id)
## Data: dados_long_miss
##
## REML criterion at convergence: 362.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4555 -0.5106 0.0284 0.4259 3.1719
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 3.615 1.901
## Residual 2.393 1.547
## Number of obs: 86, groups: id, 24
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 15.36893 1.05412 14.580
## age 0.67210 0.07486 8.978
## sexM 2.27623 0.85724 2.655
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.786
## sexM -0.478 0.010summary(pool(modelo_3), conf.int = TRUE)## term estimate std.error statistic df p.value 2.5 %
## 1 (Intercept) 15.0134789 0.97387427 15.416239 93.00831 0.000000e+00 13.0795599
## 2 age 0.6948333 0.06927964 10.029402 90.84199 2.220446e-16 0.5572146
## 3 sexM 2.3538138 0.78590543 2.995034 96.87147 3.483709e-03 0.7939829
## 97.5 %
## 1 16.9473979
## 2 0.8324521
## 3 3.9136446E quanto ao número adequado de imputações?
Nos dois artigos abaixo, sugere-se um número ótimo m de imputações: sendo a eficiência do método mensurada pela fórmula: \[ \left(1 + \frac{\lambda}{m}\right)^{-1} \] Significa que em um conjunto de dados com \(56\%\) de informações faltantes, seriam necessárias \(m=5\) imputações para se obter \(90\%\) de eficiência no método.
- Rubin, D.B. (1987). Multiple Imputation for Nonresponse in Surveys. New York: John Wiley and Sons.
- Zaninotto, P.; Sacker, A. (2017) Missing Data in Longitudinal Surveys: A Comparison of Performance of Modern Techniques. link do artigo .
Conclusões:
- O modelo ajustado com imputação apresentou resultados satisfatórios, próximo ao modelo completo;
- Os resultados podem ser melhorados a medida que N cresce - número de bancos imputados.
- PENDENTE MELHORAR
Pendências
- Ajustar mais modelos no contexto de imputação - para poder comparar
- efeito aleatório no intercepto
- efeitos aleatórios no intercepto + slope (coeficiente angular considerando a variável age);
- interação de efeitos fixos (age & sex)