Генерировать коррелированные переменные в Stata

Question

Я хочу сгенерировать 5 коррелированных переменных в Stata. Четыре нормально распределенных со специальными средствами и стандартными отклонениями, а один - с бернулли с вероятностью 0.60.

Я пытался следовать советам в этой должности: How to generate correlated Uniform[0,1] variables

Мой код выглядит следующим образом:

matrix C =  (1,     ///                 /// 
2*sin(0.05*_pi/6), 1,  /// 
2*sin(-0.45*_pi/6), 2*sin(0.44*_pi/6), 1,   /// 
2*sin(0.22*_pi/6), 2*sin(0.33*_pi/6), 2*sin(-0.54*_pi/6), 1, /// 
2*sin(0.45*_pi/6), 2*sin(0.32*_pi/6), 2*sin(-0.22*_pi/6), 2*sin(-0.13*_pi/6), 1)  

matrix B = (40, 26, 13, 146, 0.35) 
matrix A = (9, 11, 5, 2, 1) 

corr2data var1 var2 var3 var4 var5, n(10000) corr(C) means(B) sds(A) cstorage(lower) 

replace var1 = rnormal(var1) 
replace var2 = rnormal(var2) 
replace var3 = rnormal(var3) 
replace var4 = rnormal(var4) 

replace var5 = normal(var5) 
replace var5 = rbinomial(1,var5) 

Я получил то, что я хотел более или менее в том смысле, что ценности, которые сгенерированные переменные имеют ожидаемое значение.

Однако, мой подход правильный? Если нет, как бы вы изменили код, чтобы правильно дать желаемые результаты, будучи научно обоснованными?

Answer 1

С кодом возникает несколько проблем. Во-первых, преобразование корреляционной матрицы полезно только для специального случая генерирования переменных , но вы хотите, чтобы коррелированные нормали и биномы. Во-вторых, вам не нужно повторно генерировать var1-var4 с rnormal, так как corr2data уже делает это за вас. В-третьих, ваша корреляционная матрица не является положительной (полу) определенной, поэтому код не работает так, как написано для меня. В-четвертых, вам нужно применить обратный CDF распределения Бернулли для имитации ничьей из этого распределения (это шаг 3 в связанном сообщении), а не использовать rbinomial().

Вот упрощенный пример с двумя нормалями и Бернулли:

clear 
local p=0.6 
matrix m = (10,0,0) 
matrix sd = (5,1,1) 

/* I am shooting for corr(n1,b)=0.5 and corr(n2,b)=0.75, so I exaggerate their correlations in the bottom row */ 
matrix c = /// 
(1, /// 
0.5, 1, /// 
0.64, 0.95,1)  

corr2data n1 n2 b, n(10000) corr(c) means(m) sds(sd) cstorage(lower) 

/* Steps 2-3 for the one Bernoulli variable */ 
replace b = cond(normal(b)>=(1-`p'),1,0) 

/* Check that we did things correctly */ 
corr, means 
qnorm n1 
qnorm n2 
prtest b = `p' 

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