Sparse regression with small effects

[Diogo Melo]

I'm attempting to set up a multivariate sparse regression with a random intercept term and could use a few pointers:

1) How to best set up the random intercepts (family terms)? I currently I have:

    vector[K] mu[N];

    matrix[K,K] L_Sigma_G;

    matrix[K,K] L_Sigma_R;

 

    // Random intercept covariance matrix

    L_Sigma_G = diag_pre_multiply(L_sigma_G, L_Omega_G);

    for (j in 1:n_family)

        beta_family[j] ~ multi_normal_cholesky(zeros, L_Sigma_G);

    // Residual covariance matrix

    L_Sigma_R = diag_pre_multiply(L_sigma_R, L_Omega_R);

    for (n in 1:N)

        mu[n] = w0 + w * ad[n] + beta_family[family[n]];

    y ~ multi_normal_cholesky(mu, L_Sigma_R);

    // priors

    w0 ~ normal(0,5);

    L_Omega_G ~ lkj_corr_cholesky(4);

    L_sigma_G ~ cauchy(0, 2.5);

    L_Omega_R ~ lkj_corr_cholesky(4);

    L_sigma_R ~ cauchy(0, 2.5);

Is there any way to vectorize this further?

2) The w terms are from Piironen & Vehtari 2015 for a horseshoe prior on the regression coefficients, but I think I'm having some problems with over shrinkage of the non-noise components, since they are expected to be small (around 0.1 - 0.2 sd).

So, my question is, how can I modify the priors to take the expected scale of the non-null effects into account?

I currently have: 

transformed parameters{

   

    real<lower=0> tau;

    matrix<lower=0>[K, J] lambda;

    matrix[K, J] w;

    tau = r1_global * sqrt(r2_global);

    lambda = r1_local .* sqrt_vec(r2_local);

    w = beta .* lambda * tau;

}

and in the model section:

 to_vector(beta) ~ normal(0, 1);

 to_vector(r1_local) ~ normal(0.0, 1.0);

 to_vector(r2_local) ~ inv_gamma(0.5, 0.5);

 // half cauchy for tau

 r1_global ~ normal(0.0, 1.0);

 r2_global ~ inv_gamma(0.5, 0.5);

Full model is attached, any input is much appreciated! 

Thanks,

Diogo Melo

[Bob Carpenter]

> On Jul 10, 2016, at 1:55 PM, Diogo Melo <dio...@gmail.com> wrote:
>
> I'm attempting to set up a multivariate sparse regression with a random intercept term and could use a few pointers:
>
> 1) How to best set up the random intercepts (family terms)? I currently I have:
>
> vector[K] mu[N];
> matrix[K,K] L_Sigma_G;
> matrix[K,K] L_Sigma_R;
>
> // Random intercept covariance matrix
> L_Sigma_G = diag_pre_multiply(L_sigma_G, L_Omega_G);
>
> for (j in 1:n_family)
> beta_family[j] ~ multi_normal_cholesky(zeros, L_Sigma_G);
>
> // Residual covariance matrix
> L_Sigma_R = diag_pre_multiply(L_sigma_R, L_Omega_R);
>
> for (n in 1:N)
> mu[n] = w0 + w * ad[n] + beta_family[family[n]];
>
> y ~ multi_normal_cholesky(mu, L_Sigma_R);
>
> // priors
> w0 ~ normal(0,5);
>
> L_Omega_G ~ lkj_corr_cholesky(4);
> L_sigma_G ~ cauchy(0, 2.5);
>
> L_Omega_R ~ lkj_corr_cholesky(4);
> L_sigma_R ~ cauchy(0, 2.5);
>
> Is there any way to vectorize this further?

Not really. The statement

> for (n in 1:N)
> mu[n] = w0 + w * ad[n] + beta_family[family[n]];

looks tempting, but with these types it's not straightforward:

vector[K] mu[N];
matrix[K, J] w;
vector[K] w0;
vector[J] ad[N];
int family[N];
vector[K] beta_family[n_family];

even if you turned the arrays of vectors into matrices, because
you have to broadcast w0. It'd probalby help to compute
(w * ad) using matrix arithmetic and store the result in a temp.
This will be even more advantageous when we get Stan 2.11 out,
which speeds up matrix ops tremendously.

>
> 2) The w terms are from Piironen & Vehtari 2015 for a horseshoe prior on the regression coefficients, but I think I'm having some problems with over shrinkage of the non-noise components, since they are expected to be small (around 0.1 - 0.2 sd).
>
> So, my question is, how can I modify the priors to take the expected scale of the non-null effects into account?

Could you ask this in a separate question with "horseshoe prior"
in the title to get Aki's attention? I don't know much about
these.

- Bob

P.S. Thanks for attaching the full model---it's so much easier
for us that way.