# R crashes when increasing number of studies from 4 to 5, in a regression random effect meta-analysis

**URL:** <https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150>\
**Category:** RStan\
**Created:** [September 27, 2019, 1:00pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150 "2019-09-27T13:00:24Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![anthony\_462](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/anthony_462/32/4456_2.png) [@anthony\_462](https://discourse.mc-stan.org/u/anthony_462)\
**Post date:** [September 27, 2019, 1:00pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/1 "2019-09-27T13:00:24Z")

</div>

**Data**  
There are N studies, each study i has n\_i participants and k\_i successes. For each study, we also have two covariates x\_1 and x\_2.

**Model 1**  
I try to fit a random effect regression model, assuming binomial distribution of the data. We define x as a matrix of dimensions N\times 2, whose columns are the two covariates x\_1 and x\_2.  
k\_i \sim \mathcal{Binomial}(p\_i,n\_i)  
\text{logit}(p\_i)=\theta + \theta\_i + x \cdot \beta  
\theta\_i \sim \mathcal{N}(0,\sigma)

**Prior**  
p\_{init} \sim \mathcal{Beta}(1,10)  
\theta=\text{logit}(p\_{init})  
\sigma \sim \mathcal{Exp}(\sigma\_1)  
\beta \sim \mathcal{N}(0,10)

**Model 2**  
In model 2, I tried to estimate \beta with the help of regularized horseshoe priors as in [J Piironen et al.](https://arxiv.org/pdf/1707.01694.pdf).

**Simulations**  
With the help of simulated data (see R code below), I first try to see whether both models can retrieve true values of \beta, using Rstan.  
When N\<5 (number of studies), both models work, but when N=5 or more with Model 2 (with horseshoe prior) RStudio crashes, while there is no issue with Model 1.

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/e77cb6e17dd3f649fe72cbf35846a433714952af.png)

Same problem happens when assigning `inference=0` (prior predictive check). Any idea why ?  
R version: 3.5.3.  
Rstan version: 2.19.2

> N=4 #work  
> N=5 #does not work  
> x1=rnorm(N,50,20)/100  
> x2=rbeta(N,1.5,1.5)

> x=cbind(x1,x2)  
> x=as.matrix(x,nrow=N)  
> K=dim(x)[2]  
> apply(x,2,mean)  
> theta=logit(0.1)  
> sigma=0.8  
> b=2  
> theta\_i=rnorm(N,0,sigma)  
> n=round(seq(100,200,length.out = N))  
> p=inv.logit(theta+theta\_i+b\*x2)  
> k\<-rbinom(N,n,p)
> 
> stan\_data=list(N=N,  
> K=K,  
> k=k,  
> n=n,  
> x=matrix(as.matrix(x),ncol=K),  
> beta\_1=structure(rep(0,K),dim=K),  
> beta\_2=structure(rep(10,K),dim=K),  
> sigma\_1=1,  
> scale\_global=1/((K-1)\*sqrt(N)),  
> nu\_global=1,  
> nu\_local=1,  
> slab\_scale=10,  
> slab\_df=1,  
> inference=1)
> 
> mod\<-stan\_model(file = “model1.stan”)  
> mod\<-stan\_model(file = “model2.stan”)  
> fit\<-sampling(object=mod7,  
> data=stan\_data,  
> warmup = 1000,  
> iter = 2500,  
> chains = 1,  
> cores = 4,  
> thin = 1,  
> control=list(adapt\_delta=0.8))

**Stan code Model 1 (model1.stan)**

```stan
// Stan model for simple linear regression
data {
  int<lower=0> N; // number of data items
  int<lower=0> K; // number of predictors
  int k[N];
  int n[N];
  matrix[N, K] x; // predictor matrix
  
  int inference;
  
  vector[K] beta_1;
  vector[K] beta_2;
  real sigma_1;
}
parameters {
  vector [K] beta;
  real<lower=0, upper=1> p_init;
  vector[N] theta_i;
  real<lower=0> sigma;
}
transformed parameters {
  real theta = logit(p_init);
  vector[N] xb;
  xb = x * beta;
}
model {
  //Priors
  p_init ~ beta(1,10);
  beta ~ normal(beta_1,beta_2);
  theta_i ~ normal(0,1);
  sigma ~ exponential(sigma_1);
  
  if(inference==1){
    for(i in 1:N){
      target += binomial_lpmf(k[i] | n[i], inv_logit(theta + theta_i[i] * sigma + xb[i])); // likelihood
    }
  }
}

```

**Stan code Model 2 (model2.stan)**

```stan
// Stan model for simple linear regression
data {
  int<lower=0> N; // number of data items
  int<lower=0> K; // number of predictors
  int k[N];
  int n[N];
  matrix[N, K] x; // predictor matrix

  int inference;
  
  real < lower =0 > scale_global; // scale for the half -t prior for tau
  real < lower =1 > nu_global; // degrees of freedom for the half -t priors for tau
  real < lower =1 > nu_local; // degrees of freedom for the half - t priors for lambdas
  real < lower =0 > slab_scale; // slab scale for the regularized horseshoe
  real < lower =0 > slab_df; // slab degrees of freedom for the regularized horseshoe
}
parameters {
  real<lower=0, upper=1> p_init;
  vector[N] theta_i; 
  
  vector [K] z;
  real < lower =0 > tau; // global shrinkage parameter
  vector < lower =0 >[K] lambda; // local shrinkage parameter
  real < lower =0 > caux;
}
transformed parameters {
  real theta = logit(p_init);
  vector[K] beta; // coefficients for predictors
  vector[N] xb;
  
  real < lower =0 > c; // slab scale
  vector < lower =0 >[K] lambda_tilde ; // ’ truncated ’ local shrinkage parameter
   c = slab_scale * sqrt (caux);
  lambda_tilde = sqrt ( c ^2 * square ( lambda ) ./ (c ^2 + tau ^2* square ( lambda )) );

  beta = z .* lambda_tilde * tau;
  xb = x * beta;
}
model {
  //Priors
  z ~ normal(0, 1);
  lambda ~ student_t( nu_local, 0, 1);
  tau ~ student_t( nu_global, 0 , scale_global);
  caux ~ inv_gamma (0.5* slab_df , 0.5* slab_df );
 
  p_init ~ beta(1,10);
  
  theta_i ~ normal(0,1);
  
  if(inference==1){
    for(i in 1:N){
      target += binomial_lpmf(k[i] | n[i], inv_logit(theta + theta_i[i] + xb[i])); // likelihood
    }
  }
}

```

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [September 27, 2019, 1:03pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/2 "2019-09-27T13:03:00Z")

</div>

Couple points. First, I’ve taken the liberty of changing the tag of your post because I think this might be a RStan-specific issue rather than a modelling issue. Secondly, can you clarify what you mean by “crash” so folk like @bgoodri can be able to help?

---

<div class="post-metadata">

**Author:** ![anthony\_462](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/anthony_462/32/4456_2.png) [@anthony\_462](https://discourse.mc-stan.org/u/anthony_462)\
**Post date:** [September 27, 2019, 1:46pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/3 "2019-09-27T13:46:28Z")

</div>

Thanks, I just updated the question. Rstudio crashes.

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [September 27, 2019, 2:42pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/4 "2019-09-27T14:42:55Z")

</div>

Using `-mtune = native` or not?

---

<div class="post-metadata">

**Author:** ![anthony\_462](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/anthony_462/32/4456_2.png) [@anthony\_462](https://discourse.mc-stan.org/u/anthony_462)\
**Post date:** [September 27, 2019, 3:47pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/5 "2019-09-27T15:47:09Z")

</div>

Yes. Here is my Makevars.win file:

> CXX14FLAGS=-O3 -march=native -mtune=native  
> CXX11FLAGS=-O3 -march=native -mtune=native

Should I change that?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [September 27, 2019, 4:09pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/6 "2019-09-27T16:09:06Z")

</div>

Try deleting both `-march=native -mtune=native`.

---

<div class="post-metadata">

**Author:** ![anthony\_462](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/anthony_462/32/4456_2.png) [@anthony\_462](https://discourse.mc-stan.org/u/anthony_462)\
**Post date:** [September 27, 2019, 4:37pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/7 "2019-09-27T16:37:26Z")

</div>

It works when deleting `-mtune=native`, thanks! Can you explain me what does it change exactly and why it did not work before?

---

<div class="post-metadata">

**Author:** ![bgoodri](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bgoodri/32/4451_2.png) [@bgoodri](https://discourse.mc-stan.org/u/bgoodri)\
**Post date:** [September 27, 2019, 5:10pm UTC](https://discourse.mc-stan.org/t/r-crashes-when-increasing-number-of-studies-from-4-to-5-in-a-regression-random-effect-meta-analysis/11150/8 "2019-09-27T17:10:17Z")

</div>

`-mtune=native` enables options in the compiled code that are specific to the processor in your computer and it messes that up on Windows for some processors.