RS-Interface-code/LinearInteractionsMixing.R at master · willthomson1/RS-Interface-code · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
### linear model with interactions to check robustness to multi-modality

set.seed(42)

require(rstan)

d <- 15
n <- 100

X <- data.frame(matrix(rnorm(d*n),n,d))

designMat <- model.matrix(~ . + .^2 + 0,data = X)

p <- ncol(designMat)

beta0_true <- rnorm(1)

beta_true <- rnorm(p)

beta_true[sample.int(p,p-10)] <- 0

y <- as.numeric(beta0_true + designMat%*%beta_true + 0.1*rnorm(n))

model_code <- "data{
int<lower=1> N;
int<lower=1> ncov;
real y[N];
matrix[N,ncov] x;
real<lower = 0> sigma_indic;
real mu_indic;
real<lower = 0> tau;
}

parameters{
vector[ncov] beta_raw;
vector[ncov] indic_raw;
real beta_0;
real<lower = 0> sigma_noise;
}

transformed parameters{
vector[ncov] beta;
vector<lower=0,upper=1>[ncov] indic;

indic = inv_logit(mu_indic + sigma_indic*indic_raw);

beta = tau * indic .* beta_raw;

}

model{
beta_raw ~ normal(0,1);

indic_raw ~ normal(0,1);

y ~ normal(beta_0 + x*beta,sigma_noise);

}
"

sfit <- stan(model_code = model_code,data = list(N = 100,
                                         ncov = 120,
                                         y = y,
                                         x = designMat,
                                         sigma_indic = 10,
                                         mu_indic = 0,
                                         tau = 5),
              chains = 4,iter = 4000, control = list(adapt_delta = 0.99))

require(shinystan)

launch_shinystan(sfit) ## see parameters beta[13],[28],[70],[87],[106],[113],[115],[119]