Some updates

DESCRIPTION

@@ -4,6 +4,8 @@ Depends:
     methods,
     MASS,
     BMA,
+    mlogitBMA,
+    nnet,
     coloc
 Suggests:
     ggplot2,

R/coloc.mdf.R

@@ -31,8 +31,8 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response),
     n.orig <- length(snps)
     if(n.orig<3)
         return(1)
-    prep <- coloc:::prepare.df(df1, snps, r2.trim=r2.trim, dataset=1, quiet=quiet)
-    df1 <- prep$df
+    prep <- coloc:::prepare.df(df1=df1[,c(response,snps)], r2.trim=r2.trim, dataset=1, quiet=quiet)
+    df1 <- prep$df1[,c(response,prep$snps)]
     snps <- prep$snps
 
     if(!quiet)
@@ -173,7 +173,7 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response),
         snps2<-paste("2:",unlist(modelsnps),sep="")
         if(!quiet)
             cat(".")
-        capture.output(multiglm <- multinom(as.formula(paste(response"Y~",paste(modelsnps,collapse="+"))), data=df1,maxit=1000))
+        capture.output(multiglm <- multinom(as.formula(paste(response,"~",paste(modelsnps,collapse="+"))), data=df1,maxit=1000))
         coef.1[[i]] <- coefficients(multiglm)[1,modelsnps]
         coef.2[[i]] <- coefficients(multiglm)[2,modelsnps]
         var[[i]] <- vcov(multiglm)[c(snps1,snps2),c(snps1,snps2)]
@@ -204,7 +204,8 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response),
     var.2 <- lapply(var.2, diag)
 
     if(plot.coeff) {
-        coloc:::coeff.plot(unlist(coef.1),unlist(coef.2),
+       message("Plots generated directly are deprecated.  Please plot returned object instead.")
+     coloc:::coeff.plot(unlist(coef.1),unlist(coef.2),
                    unlist(var.1),unlist(var.2),
                    eta=stats["eta.hat"],
                    main="Coefficients",
@@ -214,11 +215,15 @@ coloc.var.bma <- function(df1,snps=setdiff(colnames(df1),response),
     if(!bayes) {
         return(new("coloc",
                   result=c(stats["eta.hat"],chisquare=NA,stats["n"],stats["p"]),
-                  method="BMA"))
+                  method="BMA",
+                   plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),
+                     var1=unlist(var.1),var2=unlist(var.2),
+                     model.prob=probs[wh])))
     } else {
         return(new("colocBayes",
                   result=c(stats["eta.hat"],chisquare=NA,stats["n"],stats["p"]),
                   method="BMA",
+                   plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]),
                   ppp=stats["ppp"],
                   credible.interval=ci,
                   bayes.factor=bf))
@@ -412,7 +417,7 @@ coloc.var.test <- function(X,vars.drop=NULL, ...) {
 ##' @param k Theta has a Cauchy(0,k) prior.  The default, k=1, is equivalent to a
 ##' uniform (uninformative) prior.  We have found varying k to have little effect
 ##' on the results.
-##' @param plot.coeff \code{TRUE} if you want to generate a plot showing the
+##' @param plot.coeff DEPRECATED. Please \code{plot()} returned object instead. \code{TRUE} if you want to generate a plot showing the
 ##' coefficients from the two regressions together with confidence regions.
 ##' @param bma parameter set to \code{TRUE} when \code{coloc.test} is called by \code{coloc.bma}.  DO NOT SET THIS WHEN CALLING \code{coloc.test} DIRECTLY!
 ##' @param plots.extra list with 2 named elements, x and y, equal length
@@ -435,7 +440,7 @@ coloc.var.test <- function(X,vars.drop=NULL, ...) {
 ##' approximating the posterior distribution at \code{n.approx} distinct values.
 ##' @param nsnps The SNPs to consider as potential explanatory variables
 ##' @author Mary Fortune
-coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL,bayes=!is.null(bayes.factor),
+coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=FALSE,plots.extra=NULL,bayes=!is.null(bayes.factor),
                                n.approx=1001, level.ci=0.95,
                                bayes.factor=NULL, bma=FALSE, nsnps=2) {
   snpnum <- length(b1)
@@ -541,7 +546,8 @@ coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL,
 
   ## plots
   if(plot.coeff) {
-    coloc:::coeff.plot(b1,b2,diag(V11),diag(V22),eta=eta.hat,
+    message("Plots generated directly are deprecated.  Please plot returned object instead.")
+      coloc:::coeff.plot(b1,b2,diag(V11),diag(V22),eta=eta.hat,
                main="Coefficients",
                #         sub=paste("ppp =",format.pval(ppp$value,digits=2),"p =",format.pval(pchisq(X2,df=nsnps-1,lower.tail=FALSE),digits=2)),
                xlab=expression(b[1]),ylab=expression(b[2]))
@@ -564,19 +570,22 @@ coloc.var.test.summary <- function(b1,b2,V,k=1,plot.coeff=TRUE,plots.extra=NULL,
   if(!bayes) {
     return(new("coloc",
                result=c(eta.hat=eta.hat,chisquare=X2,n=nsnps),
-               method="single"))
+               method="single",
+                   plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh])))
   } else {
     if(!bma) {
       return(new("colocBayes",
                  result=c(eta.hat=eta.hat,chisquare=X2,n=nsnps),
                  method="single",
+                   plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]),
                  ppp=ppp$value,
                  credible.interval=cred.int,
                  bayes.factor=post.bf))
     } else {
       return(new("colocBayesBMA",
                  result=c(eta.hat=eta.hat,chisquare=X2,n=nsnps),
                  method="single",
+                   plot.data=list(coef1=unlist(coef.1),coef2=unlist(coef.2),var1=unlist(var.1),var2=unlist(var.2),model.prob=probs[wh]),
                  ppp=ppp$value,
                  bma=post.bma,
                  bayes.factor=post.bf))
@@ -661,8 +670,8 @@ singlesnp.twotrait<-function(df1,response="Y",snps=setdiff(colnames(df1),respons
     results.t2<-single.snp.tests(df.t2[,response], snp.data=snp.data.t2)
     capture.output(show.results.t2<-show(results.t2))
 
-    cat("The minimum single SNP p value for RA is: ", min(show.results.t1[,4]),"\n")
-    cat("The minimum single SNP p value for T1D is: ", min(show.results.t2[,4]),"\n")
+    message("The minimum single SNP p value for trait 1 is: ", min(show.results.t1[,4]))
+    message("The minimum single SNP p value for trait 2 is: ", min(show.results.t2[,4]))
     return(c(min(show.results.t1[,4]),min(show.results.t2[,4])))
 }
 

README.md

@@ -1,4 +1,7 @@
 colocCommonControl
 =====
 
-Repo for the R package colocCommonControl
+Repo for the R package colocCommonControl
+
+For details on the method, and for citations, please see Fortune et al (2015) [http://www.nature.com/ng/journal/v47/n7/full/ng.3330.html].  Preprint is at bioRxiv [http://biorxiv.org/content/early/2015/06/08/020651].
+