## ----knitr-opts, include=FALSE------------------------------------------------
knitr::opts_chunk$set(comment = "#",collapse = TRUE,
                      fig.align = "center",dpi = 120)

## ----load-pkgs----------------------------------------------------------------
library(mr.mashr)
set.seed(12)

## ----sim-data-1, results="hide"-----------------------------------------------
dat <- simulate_mr_mash_data(n = 150,p = 800,p_causal = 5,r = 5,pve = 0.5,
							 V_cor = 0.25)

## ----sim-data-2---------------------------------------------------------------
ntest  <- 50 
Ytrain <- dat$Y[-(1:ntest),]
Xtrain <- dat$X[-(1:ntest),]
Ytest  <- dat$Y[1:ntest,]
Xtest  <- dat$X[1:ntest,]

## ----choose-prior-1-----------------------------------------------------------
S0 <- compute_canonical_covs(r = 5,singletons = TRUE,
                             hetgrid = seq(0,1,0.25))

## ----choose-prior-2-----------------------------------------------------------
names(S0)

## ----choose-prior-3-----------------------------------------------------------
S0_ind <- compute_canonical_covs(r = 5,singletons = FALSE,
                                 hetgrid = c(0,0.001,0.01))
names(S0_ind)

## ----choose-prior-4-----------------------------------------------------------
univ_sumstats <- compute_univariate_sumstats(Xtrain,Ytrain,standardize = TRUE)
scaling_grid <- autoselect.mixsd(univ_sumstats,mult = sqrt(2))^2
S0 <- expand_covs(S0,scaling_grid)
S0_ind <- expand_covs(S0_ind,scaling_grid)

## ----fit-mr-mash-1, results="hide"--------------------------------------------
null_weight <- 0.99
non_null_weight <- 1-null_weight
w0_init <- c(null_weight, rep(non_null_weight/(length(S0)-1), (length(S0)-1)))
fit <- mr.mash(X=Xtrain,Y=Ytrain,S0=S0, w0=w0_init, update_V = TRUE,
               nthreads = 1)

## ----fit-mr-mash-2, results="hide"--------------------------------------------
w0_init <- c(null_weight, rep(non_null_weight/(length(S0_ind)-1), (length(S0_ind)-1)))
fit_ind <- mr.mash(X=Xtrain,Y=Ytrain,S0=S0_ind,w0=w0_init,update_V = TRUE,
                   nthreads = 1)

## ----plot-coefs, fig.height=3.5, fig.width=6----------------------------------
oldpar <- par(mfrow = c(1,2))
plot(dat$B,fit_ind$mu1,pch = 20,xlab = "true",ylab = "estimated",
     main = sprintf("cor = %0.3f",
	                cor(as.vector(dat$B),as.vector(fit_ind$mu1))))
abline(a = 0,b = 1,col = "royalblue",lty = "dotted")
plot(dat$B,fit$mu1,pch = 20,xlab = "true",ylab = "estimated",
     main = sprintf("cor = %0.3f",
	                cor(as.vector(dat$B),as.vector(fit$mu1))))
abline(a = 0,b = 1,col = "royalblue",lty = "dotted")
par(oldpar)

## ----prior-mixture-weights----------------------------------------------------
head(sort(fit$w0,decreasing = TRUE),n = 10)

## ----resid-var----------------------------------------------------------------
dat$V
fit$V

## ----plot-pred-test, fig.height=3.5, fig.width=6------------------------------
oldpar <- par(mfrow = c(1,2))
Ypred <- predict(fit,Xtest)
Ypred_ind <- predict(fit_ind,Xtest)
plot(Ytest,Ypred_ind,pch = 20,col = "darkblue",xlab = "true",
     ylab = "predicted",
     main = sprintf("cor = %0.3f",cor(as.vector(Ytest),as.vector(Ypred_ind))))
abline(a = 0,b = 1,col = "magenta",lty = "dotted")
plot(Ytest,Ypred,pch = 20,col = "darkblue",xlab = "true",
     ylab = "predicted",
     main = sprintf("cor = %0.3f",cor(as.vector(Ytest),as.vector(Ypred))))
abline(a = 0,b = 1,col = "magenta",lty = "dotted")
par(oldpar)