De ce nu ar trebui să utilizați imputația medie pentru datele lipsă

Am întâlnit astăzi întrebarea ce să fac cu valorile lipsă atunci când se efectuează testarea ipotezei nule sau regresia? Am văzut mulți care sugerează să facă imputare răutăcioasă. Adică, înlocuiți pur și simplu orice valoare lipsă cu media variabilei calculată din valorile observate. Eu susțin că imputarea răutăcioasă este mai rea decât a nu face nimic. Să explorăm.

Pentru început, să simulăm un vector, xdin distribuția normală aleatorie.

set.seed(2112)
x <- rnorm(100, mean = 0, sd = 1)
(mean1 <- mean(x))

(sd1 <- sd(x))

Putem vedea că media și abaterea standard sunt destul de apropiate de 0 și, respectiv, 1. În următoarea bucată de cod vom selecta aleatoriu 20% din observații și vom seta valoarea la NA. Putem calcula media și abaterea standard excluzând valorile lipsă (de ex NAs) dar setarea na.rm = TRUE. Media și abaterea standard sunt relativ apropiate.

x(sample(length(x), length(x) * 0.2, replace = FALSE)) <- NA
(mean2 <- mean(x, na.rm = TRUE))

(sd2 <- sd(x, na.rm = TRUE))

Acum vom înlocui NAam introdus mai sus cu media. Putem vedea că abaterea standard este destul de mai mică, reducând astfel varianța estimării noastre. Deoarece multe dintre testele noastre statistice se bazează pe varianță, reducerea varianței poate duce la concluzii false.

x(is.na(x)) <- mean(x, na.rm = TRUE)
(mean3 <- mean(x))

(sd3 <- sd(x))

Pentru a arăta că aceasta nu este o anomalie aleatorie pentru unicul nostru eșantion aleatoriu, să repetăm ​​cele de mai sus de 1.000 de ori.

n_samples <- 1000
percent_missing <- 0.10
sd_diffs <- data.frame(sample = 1:n_samples,
                       sd_drop_miss = numeric(n_samples),
                       sd_impute_miss = numeric(n_samples))
for(i in seq_len(n_samples)) {
    x2 <- x
    x2(sample(length(x), length(x) * percent_missing, replace = FALSE)) <- NA
    sd_diffs(i,)$sd_drop_miss <- sd(x2, na.rm = TRUE)
    x2(is.na(x2)) <- mean(x2, na.rm = TRUE)
    sd_diffs(i,)$sd_impute_miss <- sd(x2)
}

sd_diffs |> 
    reshape2::melt(id.vars="sample", variable.name="calculation_type", value.name="sd") |>
    ggplot(aes(x = sd, color = calculation_type)) +
        geom_vline(xintercept = sd(x)) +
        geom_density() +
        xlab('Standard Deviation') +
        theme_minimal()

După cum arată figura de mai sus, există o diferență semnificativă în estimările abaterii standard atunci când sunt calculate folosind numai valori observate și calculate cu valori lipsă imputate cu media. The t-testul de mai jos confirmă acest lucru.

t.test(sd_diffs$sd_drop_miss, sd_diffs$sd_impute_miss)

    Welch Two Sample t-test

data:  sd_diffs$sd_drop_miss and sd_diffs$sd_impute_miss
t = 54.288, df = 1992.4, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.04782442 0.05140925
sample estimates:
mean of x mean of y 
0.9569447 0.9073278 

Acum să luăm în considerare modul în care imputația medie poate afecta estimarea unei corelații între două variabile. Vom simula două variabile cu o corelație populațională de 0,18.

n <- 100
mean_x <- 0
mean_y <- 0
sd_x <- 1
sd_y <- 1
rho <- 0.18

set.seed(2112)
df <- mvtnorm::rmvnorm(
    n = 100,
    mean = c(mean_x, mean_y),
    sigma = matrix(c(sd_x^2, rho * (sd_x * sd_y),
                     rho * (sd_x * sd_y), sd_y^2), 2, 2)) |>
    as.data.frame() |>
    dplyr::rename(x = V1, y = V2)

cor.test(df$x, df$y)

    Pearson's product-moment correlation

data:  df$x and df$y
t = 1.8314, df = 98, p-value = 0.07008
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.01504323  0.36527878
sample estimates:
      cor 
0.1819124 

Acum vom selecta aleatoriu 20% din x valori de setat NA.

df_miss <- df
df_miss(sample(n, size = 0.2 * n, replace = FALSE),)$x <- NA
cor.test(df_miss$x, df_miss$y)

    Pearson's product-moment correlation

data:  df_miss$x and df_miss$y
t = 1.8392, df = 78, p-value = 0.06969
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.01658176  0.40543327
sample estimates:
      cor 
0.2038779 

Rețineți că p-valoarea atât pentru corelația estimată folosind setul de date complet, cât și estimată doar cu valorile observate este mai mare de 0,05 (adică nu am reuși să respingem nulul că corelația este 0).

Acum vom imputa valorile lipsă cu media și vom calcula corelația.

df_miss(is.na(df_miss$x),) <- mean(df$x, na.rm = TRUE)
cor.test(df_miss$x, df_miss$y)

    Pearson's product-moment correlation

data:  df_miss$x and df_miss$y
t = 2.0582, df = 98, p-value = 0.04223
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.007431517 0.384594022
sample estimates:
      cor 
0.2035525 

Am respinge acum nulul și am concluziona că există o corelație semnificativă statistic între x şi y chiar dacă setul nostru de date inițial din care a fost simulat nu a fost.