生成与现有变量具有预定义相关性的二元变量
Generate a binary variable with a predefined correlation to an already existing variable
对于模拟研究,我想生成一组随机变量(连续变量和二进制变量),这些随机变量与已经存在的 binary 变量具有预定义的关联,此处表示为 x
.
对于这个 post,假设 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
可以生成具有给定依赖结构的随机变量。但是,它们不提供生成对已存在变量具有集合依赖性的变量的可能性。
有谁知道如何有效地做到这一点?
谢谢!
这是一个二项式 - 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)
。
连续变量的定义更不明确 - 你真的不关心它的 mean/variance 或其他任何关于它的分布吗?
这是另一个简单的版本,它双向翻转变量:
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
如果我们看一下相关公式:
对于新向量y,如果我们保留均值,问题更容易解决。这意味着我们复制向量 x
并尝试翻转相等数量的 1 和 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()
对于模拟研究,我想生成一组随机变量(连续变量和二进制变量),这些随机变量与已经存在的 binary 变量具有预定义的关联,此处表示为 x
.
对于这个 post,假设 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
可以生成具有给定依赖结构的随机变量。但是,它们不提供生成对已存在变量具有集合依赖性的变量的可能性。
有谁知道如何有效地做到这一点?
谢谢!
这是一个二项式 - 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)
。
连续变量的定义更不明确 - 你真的不关心它的 mean/variance 或其他任何关于它的分布吗?
这是另一个简单的版本,它双向翻转变量:
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
如果我们看一下相关公式:
对于新向量y,如果我们保留均值,问题更容易解决。这意味着我们复制向量 x
并尝试翻转相等数量的 1 和 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()