这就是我的处理方法,并提供了 R 中的实现:
没有唯一的解决方案来估算缺失的数据点,使得完整(估算)数据的成对相关性等于不完整数据的成对相关性。因此,为了找到“好的”解决方案而不仅仅是“任何”解决方案,我们可以引入一个附加标准,即完整的估算数据也应与原始数据共享相同的线性回归。这引导我们采用一种相当简单的方法。
- 计算原始数据的线性回归模型。
- 找到恰好位于该回归线上的缺失值的估算值。
- 为该回归线周围的估算值生成残差的随机散布
- 缩放估算残差以强制完整估算数据的相关性等于原始数据的相关性
R 中的解决方案如下:
library(data.table)
set.seed(123)
rho = cor(dt$m,dt$r,'pairwise')
# calculate linear regression of original data
fit1 = lm(r ~ m, data=dt)
fit2 = lm(m ~ r, data=dt)
# extract the standard errors of regression intercept (in each m & r direction)
# and multiply s.e. by sqrt(n) to get standard deviation
sd1 = summary(fit1)$coefficients[1,2] * sqrt(dt[!is.na(r), .N])
sd2 = summary(fit2)$coefficients[1,2] * sqrt(dt[!is.na(m), .N])
# find where data points with missing values lie on the regression line
dt[is.na(r), r.imp := coefficients(fit1)[1] + coefficients(fit1)[2] * m]
dt[is.na(m), m.imp := coefficients(fit2)[1] + coefficients(fit2)[2] * r]
# generate randomised residuals for the missing data, using the s.d. calculated above
dt[is.na(r), r.ran := rnorm(.N, sd=sd1)]
dt[is.na(m), m.ran := rnorm(.N, sd=sd2)]
# function that scales the residuals by a factor x, then calculates how close correlation of imputed data is to that of original data
obj = function(x, dt, rho) {
dt[, r.comp := r][, m.comp := m]
dt[is.na(r), r.comp := r.imp + r.ran*x]
dt[is.na(m), m.comp := m.imp + m.ran*x]
rho2 = cor(dt$m.comp, dt$r.comp,'pairwise')
(rho-rho2)^2
}
# find the value of x that minimises the discrepencay of imputed versus original correlation
fit = optimize(obj, c(-5,5), dt, rho)
x=fit$minimum
dt[, r.comp := r][, m.comp := m]
dt[is.na(r), r.comp := r.imp + r.ran*x]
dt[is.na(m), m.comp := m.imp + m.ran*x]
rho2 = cor(dt$m.comp, dt$r.comp,'pairwise')
(rho-rho2)^2 # check that rho2 is approximately equal to rho
作为最终检查,计算完整估算数据的线性回归
并绘图以显示回归线与原始数据相同。请注意,下图适用于下面所示的大数据集,以演示此方法在大数据上的使用。
fit.comp = lm(r.comp ~ m.comp, data=dt)
plot(dt$m.comp, dt$r.comp)
points(dt$m, dt$r, col="red")
abline(fit1, col="green")
abline(fit.comp, col="blue")
mtext(paste(" Rho =", round(rho,5)), at=-1)
mtext(paste(" Rho2 =", round(rho2, 5)), at=6)
DATA
OP 示例中的原始玩具数据:
dt=structure(list(m = c(2, 0.8, NA, 1.3, 2.1, NA, NA, 2.5, 1.2, 2),
r = c(3.3, NA, 4, NA, 5.2, 2.3, 1.9, NA, 3, 2.6)),
.Names = c("m", "r"), row.names = c(NA, -10L),
class = c("data.table", "data.frame"))
更大的数据集来演示大数据
dt = data.table(m=rnorm(1e5, 3, 2))[, r:=1.5 + 1.1*m + rnorm(1e5,0,2)]
dt[sample(.N, 3e4), m:=NA]
dt[sample(which(!is.na(m)), 3e4), r := NA]