Last updated: 2020-03-06
Checks: 2 0
Knit directory: misc/analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.6.0.9000). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 3278a63. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.RData
Ignored: analysis/.Rhistory
Ignored: analysis/ALStruct_cache/
Ignored: data/.Rhistory
Ignored: data/pbmc/
Untracked files:
Untracked: .dropbox
Untracked: Icon
Untracked: analysis/GHstan.Rmd
Untracked: analysis/GTEX-cogaps.Rmd
Untracked: analysis/PACS.Rmd
Untracked: analysis/SPCAvRP.rmd
Untracked: analysis/compare-transformed-models.Rmd
Untracked: analysis/cormotif.Rmd
Untracked: analysis/cp_ash.Rmd
Untracked: analysis/eQTL.perm.rand.pdf
Untracked: analysis/eb_prepilot.Rmd
Untracked: analysis/ebpmf1.Rmd
Untracked: analysis/flash_test_tree.Rmd
Untracked: analysis/ieQTL.perm.rand.pdf
Untracked: analysis/m6amash.Rmd
Untracked: analysis/mash_bhat_z.Rmd
Untracked: analysis/mash_ieqtl_permutations.Rmd
Untracked: analysis/mixsqp.Rmd
Untracked: analysis/mr_ash_modular.Rmd
Untracked: analysis/mr_ash_parameterization.Rmd
Untracked: analysis/nejm.Rmd
Untracked: analysis/normalize.Rmd
Untracked: analysis/pbmc.Rmd
Untracked: analysis/poisson_transform.Rmd
Untracked: analysis/pseudodata.Rmd
Untracked: analysis/qrnotes.txt
Untracked: analysis/ridge_iterative_splitting.Rmd
Untracked: analysis/sc_bimodal.Rmd
Untracked: analysis/shrinkage_comparisons_changepoints.Rmd
Untracked: analysis/susie_en.Rmd
Untracked: analysis/susie_z_investigate.Rmd
Untracked: analysis/svd-timing.Rmd
Untracked: analysis/temp.Rmd
Untracked: analysis/test-figure/
Untracked: analysis/test.Rmd
Untracked: analysis/test.Rpres
Untracked: analysis/test.md
Untracked: analysis/test_qr.R
Untracked: analysis/test_sparse.Rmd
Untracked: analysis/z.txt
Untracked: code/multivariate_testfuncs.R
Untracked: code/rqb.hacked.R
Untracked: data/4matthew/
Untracked: data/4matthew2/
Untracked: data/E-MTAB-2805.processed.1/
Untracked: data/ENSG00000156738.Sim_Y2.RDS
Untracked: data/GDS5363_full.soft.gz
Untracked: data/GSE41265_allGenesTPM.txt
Untracked: data/Muscle_Skeletal.ACTN3.pm1Mb.RDS
Untracked: data/Thyroid.FMO2.pm1Mb.RDS
Untracked: data/bmass.HaemgenRBC2016.MAF01.Vs2.MergedDataSources.200kRanSubset.ChrBPMAFMarkerZScores.vs1.txt.gz
Untracked: data/bmass.HaemgenRBC2016.Vs2.NewSNPs.ZScores.hclust.vs1.txt
Untracked: data/bmass.HaemgenRBC2016.Vs2.PreviousSNPs.ZScores.hclust.vs1.txt
Untracked: data/eb_prepilot/
Untracked: data/finemap_data/fmo2.sim/b.txt
Untracked: data/finemap_data/fmo2.sim/dap_out.txt
Untracked: data/finemap_data/fmo2.sim/dap_out2.txt
Untracked: data/finemap_data/fmo2.sim/dap_out2_snp.txt
Untracked: data/finemap_data/fmo2.sim/dap_out_snp.txt
Untracked: data/finemap_data/fmo2.sim/data
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.config
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.k
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.k4.config
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.k4.snp
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.ld
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.snp
Untracked: data/finemap_data/fmo2.sim/fmo2.sim.z
Untracked: data/finemap_data/fmo2.sim/pos.txt
Untracked: data/logm.csv
Untracked: data/m.cd.RDS
Untracked: data/m.cdu.old.RDS
Untracked: data/m.new.cd.RDS
Untracked: data/m.old.cd.RDS
Untracked: data/mainbib.bib.old
Untracked: data/mat.csv
Untracked: data/mat.txt
Untracked: data/mat_new.csv
Untracked: data/matrix_lik.rds
Untracked: data/paintor_data/
Untracked: data/temp.txt
Untracked: data/y.txt
Untracked: data/y_f.txt
Untracked: data/zscore_jointLCLs_m6AQTLs_susie_eQTLpruned.rds
Untracked: data/zscore_jointLCLs_random.rds
Untracked: explore_udi.R
Untracked: output/fit.k10.rds
Untracked: output/fit.varbvs.RDS
Untracked: output/glmnet.fit.RDS
Untracked: output/test.bv.txt
Untracked: output/test.gamma.txt
Untracked: output/test.hyp.txt
Untracked: output/test.log.txt
Untracked: output/test.param.txt
Untracked: output/test2.bv.txt
Untracked: output/test2.gamma.txt
Untracked: output/test2.hyp.txt
Untracked: output/test2.log.txt
Untracked: output/test2.param.txt
Untracked: output/test3.bv.txt
Untracked: output/test3.gamma.txt
Untracked: output/test3.hyp.txt
Untracked: output/test3.log.txt
Untracked: output/test3.param.txt
Untracked: output/test4.bv.txt
Untracked: output/test4.gamma.txt
Untracked: output/test4.hyp.txt
Untracked: output/test4.log.txt
Untracked: output/test4.param.txt
Untracked: output/test5.bv.txt
Untracked: output/test5.gamma.txt
Untracked: output/test5.hyp.txt
Untracked: output/test5.log.txt
Untracked: output/test5.param.txt
Unstaged changes:
Modified: analysis/ash_delta_operator.Rmd
Modified: analysis/ash_pois_bcell.Rmd
Modified: analysis/minque.Rmd
Modified: analysis/mr_missing_data.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/mr_missing_data2.Rmd) and HTML (docs/mr_missing_data2.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 3278a63 | Matthew Stephens | 2020-03-06 | workflowr::wflow_publish(“../analysis/mr_missing_data2.Rmd”) |
The idea here is to look at multiple regression with missing data, following up on the idea to use a pseudo-data approach.
Let \(\gamma_{ir}\) be an indicator for whether \(Y_{ir}\) is non-missing (so 0 indicates missing), and \(n_{rs}\) denote the number of non-missing observations in both \(r\) and \(s\): \[n_{rs} = \sum_i \gamma_{ir}\gamma_{is}.\]
Let \(W\) denote the anologue of \((1/n)X'Y\) in the full data case, based on the non-missing entries of \(Y\). That is, \[W_{jr} := (1/n_{rr}) \sum_i X_{ij}Y_{ir}\gamma_{ir}\] Let \(R\) denote the empirical covariance matrix of the covariates: \[R:=(1/n)X'X\]
Let \(X'X^{(r,s)}\) denote the matrix \(X'X\) computed using only those individuals \(i\) for which both \(r\) and \(s\) are non-missing (\(\gamma_{ir}\gamma_{is} = 1\)):
\[(X'X)^{(r,s)}_{jk} := \sum_{i} \gamma_{ir}\gamma_{is} X_{ij} X_{ik}\]
Let \(\tilde{V}\) denote the matrix \[\tilde{V}_{rs}:= V_{rs}n_{rs}/(n_{rr}n_{ss}).\]
We propose to make the following approximations:
A1: Treat the matrix \(W\) as approximate sufficient statistic, and do inference based on \(L(b) = p(W|b)\).
A2: \(X'X^{(r,s)} \approx n_{rs} R\).
With these approximations it can be shown that:
\[W \sim MN(Rb, R, \tilde{V})\] And we can use a transformation \(T\) such that \(TRT' = I\) to give: \[TW \sim MN(TRb, I, \tilde{V})\]
So we can solve the problem by fitting a MMR with outcome \(TW\) and covariates \(TR\).
In the special case with no missing data this procedure will be exact.
We have an MMR with \(n\) observations in \(r\) conditions and \(r\) covariates. \[Y_{n \times r} = X_{n \times p} b_{p \times r} + E_{n \times r\] We allow that the covariances of the residuals may be correlated (covariance matrix \(V\)), so the rows of \(E_{i\cdot} \sim N(0,V)\). Assume for now that \(V\) is known.
This model can be rewritten in terms of the matrix normal distribution: \[Y \sim MN( Xb , I, V)\]. The log-likelihood for \(b\) is (up to a constant): \[l(b) = -0.5 \text{tr}[V^{-1}(Y-Xb)'(Y-Xb)]\]
Ignoring terms that don’t depend on \(b\) we get \[l(b) = -0.5 \text{tr}[V^{-1}(b'X'Xb - 2Y'Xb)] + \text{const}\] so we can see that the summary data \(Y'X\) (or equivalently its transpose, \(X'Y\)) is a sufficient statistic for \(b\) (assuming \(X\) or \(X'X\) are known).
Let \(\gamma_{ir}\) be an indicator for whether \(Y_{ir}\) is non-missing (so 0 indicates missing).
Now define \(W\) to be \(X'Y\) computed using only the non-missing entries of \(Y\): \[W_{jr} = \sum_i X_{ij}Y_{ir}\gamma_{ir}\]
A first observation: O1 In the presence of missing data, \(W\) is not sufficient for \(b\) (even if \(\gamma,X\) are known).
[An exception: if every row of \(Y\) contains only one non-missing entry then \(W\) is sufficient]
As an aside, we have Observation O2: the way we compute the single-effect regression Bayes factors in susieR exploits this sufficiency, so they are not “valid” (ie not the actual correct BFs) in the presence of missing data.
The following simple example illustrates why \(W\) is not sufficient. Consider the case with \(r=2\) and \(Y1,Y2\) are completely correlated (\(V=11'\)), and let \(X\) be a column vector of all \(1\)s, so we are just estimating the two means, which we will call \((\mu_1,\mu_2)\).
Suppose the first observation is \((Y_1,Y_2) = (0,0)\). This tells us that \(\mu_1 = \mu_2\) (because the complete correlation implies that \(Y_1-Y_2\) is a constant and equal to \(\mu_1-\mu_2\)).
Now suppose we have a bunch of other observations where sometimes \(Y_1\) is missing, sometimes \(Y_2\) is missing, and sometimes both are observed (in which case we will know they are equal). Since we know \(\mu_1=\mu_2\) from the first observation, and the \(Y\)s are completely correlated, we know \(Y_1=Y_2\) with probability 1, so we can perfectly impute all the missing data. The mles are then easily computed based on these imputed data. But we can’t do this imputation if we only know the marginal sums of the observed \(Y_1\) and observed \(Y_2\).
This leads to a first Approximation A1: do inference for \(b\) based on \(L(b):=p(W|b)\). That is, pretend we have only observed \(W\) and not \(Y\). Because \(W\) is not sufficient, this will lose some efficiency. But it should still provide reasonable inference and will simplify things.
TBD