Fake data with interaction

[Elef]

Hi,

I would like to generate a fake response using rstan from a model with interaction of Continuous:Categorical variable.

I have made the following model with two variables,

• continuous with two levels and

• categorical of four levels.

  sim_data <- 
   '
  data {
  int N;
  real alpha;
  real Category[N];       // Coefficient for each level
  real Category_sigma[N]; // Uncertainty for the coefficient 
  real x1_b;       //Coefficient for the continues 
  real x1_b_sigma; // Uncertainty for the coefficient
  real x1_Category[N];       // Interaction coefficients for each level
  real x1_Category_sigma[N]; // Uncertainty for the coefficients
  real Category_grid[N];
  real<lower=0> sigma;
  
  real x1_grid[N];
  }
  
  parameters {}
  model {}
  
  generated quantities {
  real y[N];
  real x1_beta;
  real Category_betas[N];
  real x1_Cat_beta[N];
  
  x1_beta = normal_rng(x1_b, x1_b_sigma);
  Category_betas = normal_rng(Category, Category_sigma);
  x1_Cat_beta =  normal_rng(x1_Category, x1_Category_sigma);
  
  for (n in 1:N) {
  y[n] = normal_rng(
    alpha + Category_betas[n] +  // Categorical 
    x1_beta * x1_grid[n] +       // Continues
    x1_Cat_beta[n] * x1_grid[n] * Category_grid[n]  // Interaction
    ,sigma);
    }
  }
  '
  
  compiled_example <- stan_model(model_code = sim_data)
  

and the fake input is

Category <- c('A', 'B','C', 'D')
x1_grid <- seq(0, 1, 1)


data_grid <- expand.grid(Category = Category, 
                        x1_grid  = x1_grid)

Then I sample by setting my own coefficients using the following,

N = nrow(data_grid)
sim_ex <- sampling(
  compiled_example, algorithm = 'Fixed_param',
  data = list(  N = N, 
                alpha = 0,
                Category = rep(c(49.4, 50, 50, 52.9), length.out = N),  # Coefficient for each level
                Category_sigma = rep(c(2, 1.8, 1.7, 2), length.out = N),#Uncertainty for the coefficients
                x1_b = -17,       # Coefficient for the continues 
                x1_b_sigma = 2.3, # Uncertainty for the coefficient
                x1_Category = rep(c(2, -1.5, .3, -.9), length.out = N),         # Interaction Coefficients 
                x1_Category_sigma = rep(c(.4, .4, .4, .4), length.out = N), #Uncertainty for the coefficients
                
                x1_grid = data_grid$x1_grid,
                Category_grid = rep(c(1, 2, 3, 4), length.out = N),

                sigma = 1
                
  ), refresh = 0)

Now, for the simpler case without interaction seems to work. But there is a mistake in the interaction and I don’t know how to solve it.

Also if you have links of examples in simulating fake data with Stan would be great to share them because I have found only a couple (1)(2)

Thanks!

[bbbales2]

@Elef

But there is a mistake in the interaction and I don’t know how to solve it.

Are you getting an error at runtime or is it that the results do not look correct?

[Elef]

Are you getting an error at runtime or is it that the results do not look correct?

Hi,
It runs perfectly fine. My problem is that I don’t know how to model the interaction.

[bbbales2]

I see in your code you have an interaction term:

y[n] = normal_rng(
  alpha + Category_betas[n] +  // Categorical 
  x1_beta * x1_grid[n] +       // Continues
  x1_Cat_beta[n] * x1_grid[n] * Category_grid[n]  // Interaction
  ,sigma);

Do you mean that you’re having a problem fitting the data once you generate it?

[Elef]

Hi,

I have found the mistake, so I post the code that it works as I would expect.

The packages I used

library(dplyr)
library(rstan)
library(tidybayes)
library(reshape2)
library(nlme)
library(AlgDesign)
library(stringr)

The model with one discrete variable, 3 continuous and a 2-way interaction of the discrete with a continuous

sim_data_TA4 ←
’
data {
int N;
real CD_Base[N];
real CD_Base_sigma[N];
real Silicone_b;
real Silicone_b_sigma;
real Silicone_b_sq;
real Silicone_b_sq_sigma;
real Ts_b;
real Ts_b_sigma;
real Gl_b;
real Gl_b_sigma;
real CDb_Ts[N];
real CDb_Ts_sigma[N];
real<lower=0> sigma;

real Silicone_x[N];
real Ts_x[N];
real Gl_x[N];
}

parameters {}
model {}
generated quantities {
real y[N];
real S_beta;
real S_beta_sq;
real CD_betas[N];
real Ts_beta;
real Gl_beta;
real CDb_Ts_beta[N];

S_beta = normal_rng(Silicone_b, Silicone_b_sigma);
S_beta_sq = normal_rng(Silicone_b_sq, Silicone_b_sq_sigma);
for(n in 1:N)
CD_betas[n] = normal_rng(CD_Base[n], CD_Base_sigma[n]);
Ts_beta = normal_rng(Ts_b, Ts_b_sigma);
Gl_beta = normal_rng(Gl_b, Gl_b_sigma);
for(n in 1:N)
CDb_Ts_beta[n] = normal_rng(CDb_Ts[n], CDb_Ts_sigma[n]); // coef for interaction

for (n in 1:N) {
y[n] = normal_rng(
CD_betas[n] +
S_beta * Silicone_x[n] +
S_beta_sq * square(Silicone_x[n]) +
Ts_beta * Ts_x[n] +
Gl_beta * Gl_x[n] +
CDb_Ts_beta[n] * Ts_x[n] // interaction
,sigma);
}
}
’

compiled_TA4 ← stan_model(model_code = sim_data_TA4)

Make a fake DoE

cand_list = expand.grid(CD_levels = c(‘B’, ‘F’),
S_levels = c(0, .5, 1.5, 3),
Ts_levels = c(2.4, 5.6, 7.5),
Gl_levels = c(0, 1.5, 3))

doe_grid ← optFederov( ~ Ts_levels * CD_levels + S_levels + I(S_levels^2) + Gl_levels,
data = cand_list, nTrials = 22, crit=“D”)$design

doe_grid

Set the coefficients as I like

beta_S ← data.frame(beta = -19.7, sd_beta = .07) # beta = NA, sd_beta = 1.7
beta_S_sq ← data.frame(beta = 3.4, sd_beta = .02) # beta = NA, sd_beta = .5
beta_Ts ← data.frame(beta = -.59, sd_beta = .1) # beta = -.38, sd_beta = .48
beta_Gl ← data.frame(beta = -.4, sd_beta = .02) # beta = .4, sd_beta = .2

beta_CD ← data.frame(CD_levels = c(‘B’, ‘F’),
betas = c(51.6, 52.8),
sd_betas = c(.28, .11))
#Here are the coefficients for the the interaction
beta_CD_TS ← data.frame(CD_levels = c(‘B’, ‘F’),
betas_interaction = c(0, 1.7),
sd_betas_interaction = c(.07, .07))

beta_CD1 ← left_join(doe_grid, beta_CD, by = ‘CD_levels’)
beta_CD2 ← left_join(beta_CD1, beta_CD_TS, by = ‘CD_levels’)

Sample

N = beta_CD2 %>% nrow
sim_TA4 ← sampling(
compiled_TA4, algorithm = ‘Fixed_param’,
data = list( N = N,
CD = beta_CD2$betas,
CD_sigma = beta_CD2$sd_betas,

            S_b = beta_S$beta,
            S_b_sigma = beta_S$sd_beta,
            
            S_b_sq = beta_S_sq$beta,
            S_b_sq_sigma = beta_S_sq$sd_beta,
            
            CDb_Ts = beta_CD2$betas_interaction,
            CDb_Ts_sigma = beta_CD2$sd_betas_interaction,
            
            Ts_b = beta_Ts$beta,
            Ts_b_sigma = beta_Ts$sd_beta,
            
            Gl_b = beta_Gl$beta, 
            Gl_b_sigma = beta_Gl$sd_beta,
            
            S_x = doe_grid$S_levels,
            Ts_x = doe_grid$Ts_levels,
            Gl_x = doe_grid$Gl_levels,
            sigma = 6

), refresh = 0)

Test (apologies for using lm function in here :D )

fake_data0 ← rstan::extract(sim_TA4, pars = ‘y’) %>% as.data.frame %>% .[1:4, ]
sample_nr ←
rownames(fake_data0) %>% as.numeric %>% summary %>% .[c(‘Min.’, ‘Max.’)] %>% str_c(collapse = “_”)

fake_data0 ← data.frame(doe_grid, fake_data0 %>% t) %>%
melt(id = c(‘CD_levels’, ‘S_levels’, ‘Ts_levels’, ‘Gl_levels’)) %>%
rename(y_fake = value) %>% select(-variable)

mod_fake0 ← lm(y_fake ~ poly(S_levels, 2, raw = T) + CD_levels * Ts_levels + Gl_levels,
data = fake_data0)

fake_data0 %>% head(10)
mod_fake0 %>% summary

fake_coef ← broom::tidy(mod_fake0) %>%
mutate(
low = estimate - 2 * std.error,
high = estimate + 2 * std.error,
Sample_size = sample_nr
) %>% .[, ]

ggplot(fake_coef, aes(estimate, term, xmin = low, xmax = high, height = 0)) +
geom_point() +
geom_vline(xintercept = 0) +
geom_errorbar(colour = ‘red’, width = 0.25)

[Elef]

The problem was that the fake response didn’t have the desired properties because I made a mess with the part of interaction. Thanks for your time :)

[bbbales2]

Glad it worked out in the end!