приписывать недостающие данные, а форсирование коэффициент корреляции остаются теми же

Question

Рассмотрим следующий (Excel) набор данных:

m | r 
----|------ 
2.0 | 3.3 
0.8 | 
    | 4.0 
1.3 | 
2.1 | 5.2 
    | 2.3 
    | 1.9 
2.5 | 
1.2 | 3.0 
2.0 | 2.6 

Моя цель состоит в том, чтобы заполнить недостающие значения, используя следующее условие:

Обозначим через R парную корреляцию между двумя столбцами (около 0,68). Обозначим через R * корреляцию после пустые ячейки были заполнены. Заполните таблицу так, чтобы (R - R *)^2 = 0. Это, я хочу сохранить корреляционную структуру данных целыми.

До сих пор я сделал это с помощью Matlab:

clear all; 

m = xlsread('data.xlsx','A2:A11') ; 
r = xlsread('data.xlsx','B2:B11') ; 

rho = corr(m,r,'rows','pairwise'); 

x0 = [1,1,1,1,1,1]; 
lb = [0,0,0,0,0,0]; 
f = @(x)my_correl(x,rho); 

SOL = fmincon(f,x0,[],[],[],[],lb) 

где функция my_correl является:

function X = my_correl(x,rho) 

sum_m = (11.9 + x(1) + x(2) + x(3)); 
sum_r = (22.3 + x(1) + x(2) + x(3)); 
avg_m = (11.9 + x(1) + x(2) + x(3))/8; 
avg_r = (22.3 + x(4) + x(5) + x(6))/8; 
rho_num = 8*(26.32 + 4*x(1) + 2.3*x(2) + 1.9*x(3) + 0.8*x(4) + 1.3*x(5) + 2.5*x(6)) - sum_m*sum_r; 
rho_den = sqrt(8*(22.43 + (4*x(1))^2 + (2.3*x(2))^2 + (1.9*x(3))^2) - sum_m^2)*sqrt(8*(78.6 + (0.8*x(4))^2 + (1.3*x(5))^ + (2.5*x(6))^2) - sum_r^2); 

X = (rho - rho_num/rho_den)^2; 

end 

Эта функция вычисляет корреляцию вручную, где каждый отсутствующих данных является переменной x(i) ,

Проблема: мой фактический набор данных имеет более 20 000 наблюдений. Невозможно создать эту формулу вручную вручную.

Как я могу заполнить свой набор данных?

Примечание 1: Я открыт для использования альтернативных языков, таких как Python, Julia или R. Matlab, это мой единственный по умолчанию.

Примечание 2: ответ на 100 баллов будет присужден. Обещай сейчас.

Answer 1

Это, как я хотел бы подойти к нему, с реализацией в R при условии:

Существует не уникальное решения для вменения недостающих точек данных, например, что парная корреляция полных (вмененных) данных равна к парной корреляции неполных данных. Поэтому, чтобы найти «хорошее» решение, а не просто «любое» решение, мы можем ввести дополнительные критерии, согласно которым все вмененные данные также должны иметь одинаковую линейную регрессию с исходными данными. Это приводит нас к довольно простому подходу.

1. расчет линейной регрессионной модели для исходных данных.
2. найти вмененные значения для отсутствующих значений, которые будут точно соответствовать этой линии регрессии.
3. генерировать случайный разброс остатков для вмененных значений вокруг этой линии регрессии
4. шкалы вмененные остатки, чтобы заставить корреляцию полных вмененных данных, чтобы быть равным, что из исходных данных

Раствор как это в 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]