Fake data with interaction

Modeling
5

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)

1 Like