Блог

Для исследования моделирования я хочу сгенерировать набор случайных величин (как непрерывных, так и двоичных), которые имеют предопределенные ассоциации с уже существующей двоичной переменной, обозначаемой здесь как x .

Для этого сообщения предположим, что x оно сгенерировано в соответствии с приведенным ниже кодом. Но помните: в реальной жизни x это уже существующая переменная.

 set.seed(1245)
x <- rbinom(1000, 1, 0.6)
 

Я хочу сгенерировать как двоичную переменную, так и непрерывную переменную. Я выяснил, как сгенерировать непрерывную переменную (см. Код ниже)

 set.seed(1245)

cor <- 0.8 #Correlation 
y <- rnorm(1000, cor*x, sqrt(1-cor^2))
 

Но я не могу найти способ сгенерировать двоичную переменную, которая коррелирует с уже существующей переменной x . Я нашел несколько пакетов R, например, copula которые могут генерировать случайные переменные с заданной структурой зависимостей. Однако они не предоставляют возможности генерировать переменные с заданной зависимостью от уже существующей переменной.

Кто-нибудь знает, как сделать это эффективным способом?

Спасибо!

Если мы посмотрим на формулу для корреляции:

Для нового вектора y, если мы сохраним среднее значение, проблему будет легче решить. Это означает, что мы копируем вектор x и пытаемся перевернуть равное количество единиц и 0, чтобы получить предполагаемое значение корреляции.

Если мы допустим E(X) = E(Y) = x_bar , и E(XY) = xy_bar , то для данного rho мы упростим вышесказанное до:

 (xy_bar - x_bar^2) / (x_bar - x_bar^2) =  rho
 

Решаем и получаем:

 xy_bar = rho * x_bar   (1-rho)*x_bar^2
 

И мы можем вывести функцию для переворачивания чисел 1 и 0, чтобы получить результат:

 create_vector = function(x,rho){

  n = length(x)
  x_bar = mean(x)
  xy_bar = rho * x_bar   (1-rho)*x_bar^2
  toflip = sum(x == 1) - round(n * xy_bar)

  y = x
  y[sample(which(x==0),toflip)] = 1
  y[sample(which(x==1),toflip)] = 0
  return(y)
}
 

Для вашего примера это работает:

 set.seed(1245)
x <- rbinom(1000, 1, 0.6)
cor(x,create_vector(x,0.8))
[1] 0.7986037
 

Существуют некоторые экстремальные комбинации предполагаемых rho и p, в которых вы можете столкнуться с проблемами, например:

 set.seed(111)

res = lapply(1:1000,function(i){
             
              this_rho = runif(1)
              this_p = runif(1)
              x = rbinom(1000,1,this_p)
              data.frame(
                intended_rho = this_rho,
                p = this_p,
                resulting_cor = cor(x,create_vector(x,this_rho))
              )
           })

res = do.call(rbind,res)

ggplot(res,aes(x=intended_rho,y=resulting_cor,col=p))   geom_point()
 

Вот биномиальная — формула для q зависит только от среднего значения x и желаемой корреляции.

 set.seed(1245)
cor <- 0.8
x <- rbinom(100000, 1, 0.6)
p <- mean(x)
q <- 1/((1-p)/cor^2 p)
y <- rbinom(100000, 1, q)
z <- x*y
cor(x,z)
#> [1] 0.7984781
 

Это не единственный способ сделать это — обратите внимание, что mean(z) это всегда меньше, чем mean(x) в этой конструкции.

Непрерывная переменная еще менее четко определена — вас действительно не волнует ее среднее значение / дисперсия или что-нибудь еще в ее распределении?

Вот еще одна простая версия, в которой переменная переворачивается в обе стороны:

 set.seed(1245)
cor <- 0.8
x <- rbinom(100000, 1, 0.6)
p <- mean(x)
q <- (1 cor/sqrt(1-(2*p-1)^2*(1-cor^2)))/2
y <- rbinom(100000, 1, q)
z <- x*y (1-x)*(1-y)
cor(x,z)
#> [1] 0.8001219
mean(z)
#> [1] 0.57908