g_vec <- rbinom(100000, 2, 0.1)
env_vec <- round(runif(100000,min=0,max=6))
qt_vec <- rnorm(100000) + 0.1 * g_vec + 0.05 * env_vec + 0.05 * g_vec * env_vec
res <- env_interaction.calc(qt_vec, g_vec, env_vec)
Run the code above in your browser using DataLab