rmcelreath / glmer2stan

Define Stan models using glmer-style (lme4) formulas
55 stars 11 forks source link

Enhancement: Feed model code directly into glmer2stan #7

Open cgonza12 opened 9 years ago

cgonza12 commented 9 years ago

Enhancement: Feed model code directly into glmer2stan

I noticed it was kind of difficult to re-feed an edited STAN model into the glmer2stan function. So I also included an argument, 'mymodel', which allows the user to specify an R object of text as the STAN model. The idea being you run the the model with sampling turned off first, make your changes save that text to an object and rerun with tweaked priors. my function can be found here.

example3_noSamples <- lmer2stan( Reaction ~ Days + (Days | subject_index), data=sleepstudy,
                               calcWAIC=T,
                               warmup=100, 
                               iter = 500, 
                               chains=4,
                               sample=F) 

cat(example3_noSamples$model) # stan model as text

example3_priors = "data{
  int N;
  real Reaction[N];
real Days[N];
int subject_index[N];
int N_subject_index;
}

transformed data{
vector[2] zeros_subject_index;
for ( i in 1:2 ) zeros_subject_index[i] <- 0;
}

parameters{
real Intercept;
real beta_Days;
real<lower=0> sigma;
vector[2] vary_subject_index[N_subject_index];
cov_matrix[2] Sigma_subject_index;
}

model{
real vary[N];
real glm[N];
// Priors
Intercept ~ normal( 0 , 100 );
beta_Days ~ normal( 10 , 2 ); // Hey! Look! An informed prior!
sigma ~ uniform( 0 , 100 );
// Varying effects
for ( j in 1:N_subject_index ) vary_subject_index[j] ~ multi_normal( zeros_subject_index , 
Sigma_subject_index );
// Fixed effects
for ( i in 1:N ) {
vary[i] <- vary_subject_index[subject_index[i],1]
+ vary_subject_index[subject_index[i],2] * Days[i];
glm[i] <- vary[i] + Intercept
+ beta_Days * Days[i];
}
Reaction ~ normal( glm , sigma );
}

generated quantities{
real dev;
real vary[N];
real glm[N];
dev <- 0;
for ( i in 1:N ) {
vary[i] <- vary_subject_index[subject_index[i],1]
+ vary_subject_index[subject_index[i],2] * Days[i];
glm[i] <- vary[i] + Intercept
+ beta_Days * Days[i];
dev <- dev + (-2) * normal_log( Reaction[i] , glm[i] , sigma );
}
}"

example3_withpriors <- myglmer2stan( Reaction ~ Days + (Days | subject_index), data=sleepstudy,
                        calcWAIC=T,
                        warmup=nwarm, 
                        iter = niter,
                        chains=chains, 
                        mymodel=example3_priors # include object with edited model text
) 

stanmer(example3_withpriors)

glmer2stan model: Reaction ~ Days + (Days | subject_index) [gaussian]

Level 1 estimates:
            Expectation StdDev   2.5%  97.5%
(Intercept)      249.98   7.97 233.28 264.58
Days              10.25   1.34   7.67  13.04
sigma             25.81   1.48  22.97  28.73

Level 2 estimates:
(Std.dev. and correlations)

Group: subject_index (18 groups / imbalance: 0)

                    (1)  (2)
  (1) (Intercept) 32.31   NA
  (2) Days         0.00 7.33

DIC: 1712   pDIC: 32   Deviance: 1648.2

WAIC: 1720   pWAIC: 33.2   -2*lppd: 1653.2