SER model with summary STAT bhat shat

library(susieR)
set.seed(1)
n = 800
p = 1000
beta = rep(0,p)
beta[1] = 1
beta[50] = 1
beta[300] = 1
beta[400] = 1
X = matrix(rnorm(n*p),nrow=n,ncol=p)
y = c(X %*% beta + rnorm(n))

X.cs = susieR:::safe_colScale(X, center=TRUE, scale=TRUE)
X.c = susieR:::safe_colScale(X, center=TRUE, scale=FALSE)

# summary stat
bhat = numeric(p)
se = numeric(p)
sigmahat = numeric(p)
for(j in 1:p){
  m = summary(lm(y~X[,j]))
  bhat[j] = m$coefficients[2,1]
  se[j] = m$coefficients[2,2]
  sigmahat[j] = m$sigma
}
R = cor(X)

Know \(bhat\), \(shat\), \(XtX\) (or Cov, or Cor), \(n\), (\(var(y)\))

We convert the \(XtX\) to Cor (R), since the variance from \(XtX\) might be wrong. We can infer the diagonal element of \(XtX\) using bhat and shat.

\[ \hat{z}_{j} = \frac{\hat{\beta}_{j}}{\hat{s}_{j}} = \frac{\mathbf{x}_{cj}^{T}\mathbf{y}_{c}}{\hat{\sigma}_{j}\sqrt{\mathbf{x}_{cj}^{T}\mathbf{x}_{cj}}} = \frac{\mathbf{x}_{csj}^{T}\frac{1}{\alpha_{y}}\mathbf{y}_{c}}{\frac{1}{\alpha_{y}}\hat{\sigma}_{j}\sqrt{\mathbf{x}_{csj}^{T}\mathbf{x}_{csj}}} \]

\(\hat{\Sigma} = diag(\hat{\sigma}_{1}^2, \cdots, \hat{\sigma}_{p}^2)\)

\[ R_{j}^2 = \frac{\hat{z}_{j}^2}{\hat{z}_{j}^2+n-2} \]

Then \[ 1-R_{j}^{2} = \frac{RSS}{SST} = \frac{\hat{\sigma}_{j}^2(n-2)}{\mathbf{y}_{c}^{T}\mathbf{y}_{c}} \quad \hat{\sigma}_{j}^2 = \frac{\mathbf{y}_{c}^{T}\mathbf{y}_{c}(1-R_{j}^2)}{n-2} \] \[ \hat{\sigma}_{j}^2 = var(\mathbf{y}_{c})\frac{(n-1)(1-R_{j}^2)}{n-2} = var(\mathbf{y}_{c})\tilde{\sigma}_{j}^2 \] \[ \Rightarrow \quad \hat{s}_{j}^{2} = \frac{\hat{\sigma}_{j}^2}{\mathbf{x}_{.j}^{T}\mathbf{x}_{.j}} = var(\mathbf{y}_{c})\frac{\tilde{\sigma}_{j}^2}{\mathbf{x}_{.j}^{T}\mathbf{x}_{.j}} \]

\[ \mathbf{x}_{.j}^{T}\mathbf{x}_{.j} = var(\mathbf{y}_{c})\frac{\tilde{\sigma}_{j}^2}{\hat{s}_{j}^{2}} \quad \mathbf{x}_{.j}^{T}\mathbf{y} = var(\mathbf{y}_{c})\frac{\tilde{\sigma}_{j}^2}{\hat{s}_{j}} \hat{z}_{j} \]

Using \(X_{.}^{T}X_{.}\) and \(X_{.}\mathbf{y}_{c}\) in susie_ss is fitting \[ \begin{align*} \mathbf{y}_c &= X_{.} \boldsymbol{\beta} + \boldsymbol{\epsilon} \quad \boldsymbol{\epsilon}\sim N_{n}(0, \frac{1}{\tau}I) \\ \boldsymbol{\beta} &= \boldsymbol{\gamma} \beta \\ \boldsymbol{\gamma} &\sim Multinomial(1,\boldsymbol{\pi}) \\ \beta &\sim N(0, \frac{1}{\tau_{0}}) \end{align*} \]

In susieR, if there is no information about the variance of y, we estimate the coefficients in the standardized X, y scale. Otherwise, we estimate the coefficients in the same scale of bhat, shat.

Fit susie model using individual level data (non standardize):

res0 = susie(X = X, Y = y, 
             estimate_residual_variance = FALSE, estimate_prior_variance = FALSE,
             intercept = TRUE, standardize = FALSE, 
             scaled_prior_variance = 0.2, max_iter = 3)

Fit susie model using summary statistics (non standardize):

res1 = susie_bhat(bhat = bhat, shat = se, R = crossprod(X.c), n = n, var_y = var(y),
                  standardize = FALSE, 
                  estimate_prior_variance = FALSE, 
                  scaled_prior_variance = 0.2, max_iter = 3)

Compare the fitted \(\alpha\) and coefficients:

all.equal(res1$alpha, res0$alpha)

[1] TRUE

plot(coef(res1), coef(res0))

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

Version	Author	Date
a990c1f	zouyuxin	2018-10-08

Fit susie model using individual level data (standardize):

res2 = susie(X=X.cs, Y = (y-mean(y))/sd(y), 
             estimate_residual_variance = FALSE, estimate_prior_variance = FALSE, 
             intercept = FALSE, standardize = FALSE, scaled_prior_variance = 0.2, max_iter = 3)

Fit susie model using summary statistics (standardize):

res3 = susie_bhat(bhat = bhat, shat = se, R = cor(X), n = n, 
                  estimate_prior_variance = FALSE, 
                  scaled_prior_variance = 0.2, max_iter = 3)

Compare the fitted \(\alpha\) and coefficients:

all.equal(res2$alpha, res3$alpha)

[1] TRUE

plot(coef(res2), coef(res3))

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

Version	Author	Date
a990c1f	zouyuxin	2018-10-08

If we know z scores (intead of bhat, shat), with information about sample size n, we can fit susie model as follows

z_scores = bhat/se
res4 = susie_z(z_scores, R = cor(X), n = n, 
              estimate_prior_variance = FALSE, 
              scaled_prior_variance = 0.2, max_iter = 3)

Compare the fitted \(\alpha\) and coefficients:

all.equal(res2$alpha, res4$alpha)

[1] TRUE

plot(coef(res2), coef(res4))

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

Version	Author	Date
50fd59d	zouyuxin	2018-10-12

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] susieR_0.6.0.0364

loaded via a namespace (and not attached):
 [1] workflowr_1.1.1    Rcpp_0.12.19       matrixStats_0.54.0
 [4] lattice_0.20-35    digest_0.6.18      rprojroot_1.3-2   
 [7] R.methodsS3_1.7.1  grid_3.5.1         backports_1.1.2   
[10] magrittr_1.5       git2r_0.23.0       evaluate_0.12     
[13] stringi_1.2.4      whisker_0.3-2      R.oo_1.22.0       
[16] R.utils_2.7.0      Matrix_1.2-14      rmarkdown_1.10    
[19] tools_3.5.1        stringr_1.3.1      yaml_2.2.0        
[22] compiler_3.5.1     htmltools_0.3.6    knitr_1.20        
[25] expm_0.999-3