SuSiE Z Problem

One effect

Let’s see an example data with one true effect. The PVE is 0.5, the residual variance is 0.47^2.

X = readRDS('data/susie_X.rds') # X is from susieR package, N3finemapping, X is coloumn mean centered.
R = readRDS('data/susie_R.rds')
data = readRDS('data/sim_gaussian_75.rds')
n = data$n
beta = numeric(data$p)
beta[data$beta_idx] = data$beta_val
z = data$ss$effect/data$ss$se
susie_plot(z, y = 'z', b=beta)

Expand here to see past versions of unnamed-chunk-13-1.png:

Version	Author	Date
a6cee3a	zouyuxin	2019-01-17

Fit susie model using z scores,

fit_z = susie_z(z, R=R, max_iter = 20, L=5, track_fit = T)

Warning in susie_ss(XtX = R, Xty = z, n = 2, var_y = 1, L = L,
estimate_prior_variance = TRUE, : IBSS algorithm did not converge in 20
iterations!

The objective value is

susie_get_objective(fit_z)

[1] 409.0831

susie_plot(fit_z, y = 'PIP', b = beta)

Expand here to see past versions of unnamed-chunk-16-1.png:

Version	Author	Date
5028244	zouyuxin	2019-01-17

Now we initialize at the truth (L=1), the model converges to a lower objective value.

s_init = susie_init_coef(data$beta_idx, data$beta_val, data$p)
fit_z_init = susie_z(z, R=R, s_init = s_init)

The objective value is

susie_get_objective(fit_z_init)

[1] 299.4515

susie_plot(fit_z_init, y = 'PIP', b = beta)

Expand here to see past versions of unnamed-chunk-18-1.png:

Version	Author	Date
5028244	zouyuxin	2019-01-17
a6cee3a	zouyuxin	2019-01-17

I mentioned in the write-up, we approximate the BF(\(z_j(\sigma)\)) (eqn 7.61) using BF(\(z_j(\hat{\sigma}_j)\)). When the data contain strong association signals (\(\hat{\sigma}_j > \sigma\)), the z score approximation under-estimates the Bayes Factor.

The 7.61 in the write-up: \[ \begin{align} BF(z_j(\sigma); w) &= \sqrt{\frac{1}{1+w^2}} \exp\left( \frac{1}{2} z_j(\sigma)^2 \frac{w^2}{1+w^2} \right) \quad \quad \text{where }w^2 = (n-1)\frac{\sigma_{0}^{2}}{\sigma^2} \end{align} \]

We check the under-estimate of the BF. We fit the susie model using \(z_j(\sigma)\) and \(z_j(\hat{\sigma}_j)\) separately, with L = 1 and iter = 1, fix prior variance.

fit_z_1 = susie_ss(XtX = R, Xty = z, n = 2, var_y = 1, L = 1,
                   estimate_prior_variance = FALSE,
                   estimate_residual_variance = FALSE, max_iter = 1)

X.s = apply(X, 2, function(x) x/(sd(x)*sqrt(n-1)))
z_true = (t(X.s) %*% data$sim_y) / as.numeric(data$sigma)

fit_z_1_true = susie_ss(XtX = R, Xty = z_true, n = 2, var_y = 1, L = 1,
                        estimate_prior_variance = FALSE,
                        estimate_residual_variance = FALSE, max_iter = 1)

par(mfrow=c(1,3))
{plot(fit_z_1$lbf_mtx, fit_z_1_true$lbf_mtx, xlab = expression(logBF(z[j](hat(sigma)[j]))), ylab = expression(logBF(z[j](sigma))), main='log BF')
abline(0,1)}
{plot(fit_z_1$alpha, fit_z_1_true$alpha, xlab = expression(alpha(hat(sigma)[j])), ylab = expression(alpha(sigma)), main=expression(alpha))
abline(0,1)}
{plot(abs(fit_z_1$mu), abs(fit_z_1_true$mu), xlab = expression(abs(mu(hat(sigma)[j]))), ylab = expression(abs(mu(sigma))), main=expression(abs(mu)))
abline(0,1)}

Expand here to see past versions of unnamed-chunk-21-1.png:

Version	Author	Date
5028244	zouyuxin	2019-01-17
a6cee3a	zouyuxin	2019-01-17

It is clear from the plots that the Bayes Factor is under-estimated. The estimated \(\alpha\), \(\mu\) are smaller.

Now, we fit susie model using \(z_j(\hat{\sigma}_j)\) and \(z_j(\sigma)\), with L = 5 and estimating the prior variance.

fit_z_10_prior = susie_z(z=z, R=R, L=5)

Warning in susie_ss(XtX = R, Xty = z, n = 2, var_y = 1, L = L,
estimate_prior_variance = TRUE, : IBSS algorithm did not converge in 100
iterations!

The estimated prior variances for model using \(z_j(\hat{\sigma}_j)\) are

susie_get_prior_variance(fit_z_10_prior)

[1] 24449.0562  4961.9053  2624.2558  2017.4113   163.6263

fit_z_10_true_prior = susie_z(z = z_true, R=R, L=5)

The estimated prior variances for model using \(z_j(\sigma)\) are

susie_get_prior_variance(fit_z_10_true_prior)

[1] 604.074353   3.315658   0.000000   0.000000   0.000000

The PIP for model using \(z_j(\sigma)\):

susie_plot(fit_z_10_true_prior, y='PIP', b=beta)

Using the true z, \(z_j(\sigma)\), the estimated prior variance for some L becomes zero. Using the estimated z, \(z_j(\hat{\sigma}_j)\), the estimated prior variance for all L are non-zero.

Five effects

Let’s see an example data with 5 true effects. The PVE is 0.8, the residual variance is 0.74^2.

data = readRDS('data/sim_gaussian_475.rds')
n = data$n
beta = numeric(data$p)
beta[data$beta_idx] = data$beta_val
z = data$ss$effect/data$ss$se
susie_plot(z, y = "z", b=beta)

Expand here to see past versions of unnamed-chunk-27-1.png:

Version	Author	Date
5028244	zouyuxin	2019-01-17

Fit susie model using z scores,

fit_z = susie_z(z, R=R, max_iter = 20, L=5, track_fit = T)

Warning in susie_ss(XtX = R, Xty = z, n = 2, var_y = 1, L = L,
estimate_prior_variance = TRUE, : IBSS algorithm did not converge in 20
iterations!

The objective value is

susie_get_objective(fit_z)

[1] 525.9303

susie_plot(fit_z, y = 'PIP', b = beta)

Now we initialize at the truth, the model fails to converge.

s_init = susie_init_coef(data$beta_idx, data$beta_val, data$p)
fit_z_init = susie_z(z, R=R, s_init = s_init)

Warning in susie_ss(XtX = R, Xty = z, n = 2, var_y = 1, L = L,
estimate_prior_variance = TRUE, : IBSS algorithm did not converge in 100
iterations!

susie_get_objective(fit_z_init)

[1] 1089.854

susie_plot(fit_z_init, y = 'PIP', b = beta)

Changing the z scores to \(z_j(\sigma)\), the model converges.

z_true = (t(X.s) %*% data$sim_y) / as.numeric(data$sigma)
fit_z_5_true = susie_z(z = z_true, R=R, L=5)

The objective value is

susie_get_objective(fit_z_5_true)

[1] 1119.521

The estimated prior variances are

susie_get_prior_variance(fit_z_5_true)

[1] 1467.02402  640.34084  185.85349   53.12493    0.00000

susie_plot(fit_z_5_true, y = 'PIP', b = beta)