@@ -14,8 +14,14 @@

 #' @param taxa Taxa list.

 #' @param uniqueNames Number the various taxa.

 #' @param adjustMethod Method to adjust p-values by. Default is "FDR". Options

-#' include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

+#' include "IHW", "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

 #' "none". See \code{\link{p.adjust}} for more details.

+#' @param IHWcov Character value specifying which covariate to use when adjusting pvalues

+#' using IHW. Options include: "nnz" (number of non-zero elements per feature), 

+#' "median" (median abundance value per feature), "Amean" (adjusted mean, used for a 

+#' fitZigResults obj)

+#' @param alpha Value for p-value significance threshold when running IHW. 

+#' The default is set to 0.1 

 #' @param group One of five choices, 0,1,2,3,4. 0: the sort is ordered by a

 #' decreasing absolute value coefficient fit. 1: the sort is ordered by the raw

 #' coefficient fit in decreasing order. 2: the sort is ordered by the raw

@@ -43,21 +49,22 @@

 #' fit = fitFeatureModel(obj = lungTrim,mod=mod)

 #' head(MRcoefs(fit))

 #'

-MRcoefs<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,

-    uniqueNames=FALSE,adjustMethod="fdr",group=0,eff=0,numberEff=FALSE,counts=0,file=NULL){

+MRcoefs<-function(obj,by=2,coef=NULL,number=10,taxa=obj@taxa,

+    uniqueNames=FALSE,adjustMethod="fdr",alpha=0.1,

+    group=0,eff=0,numberEff=FALSE,counts=0,file=NULL){

-    if(length(grep("fitFeatureModel",obj$call))){

-        groups = factor(obj$design[,by])

+    if(length(grep("fitFeatureModel",obj@call))){

+        groups = factor(obj@design[,by])

         by = "logFC"; coef = 1:2;

-        tb = data.frame(logFC=obj$fitZeroLogNormal$logFC,se=obj$fitZeroLogNormal$se)

-        p  = obj$pvalues

+        tb = data.frame(logFC=obj@fitZeroLogNormal$logFC,se=obj@fitZeroLogNormal$se)

+        p  = obj@pvalues

     } else {

-        tb = obj$fit$coefficients

+        tb = obj@fit$coefficients

         if(is.null(coef)){

             coef = 1:ncol(tb)

         }

-        p=obj$eb$p.value[,by]

-        groups = factor(obj$fit$design[,by])

+        p=obj@eb$p.value[,by]

+        groups = factor(obj@fit$design[,by])

         if(eff>0){

             effectiveSamples = calculateEffectiveSamples(obj)

             if(numberEff == FALSE){

@@ -75,8 +82,24 @@ MRcoefs<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,

             tx[ii]=paste(tx[ii],seq_along(ii),sep=":")

         }

     }

-    padj = p.adjust(p,method=adjustMethod)

-

+    

+    # # adding IHW as pvalue adjustment method

+    # if(adjustMethod == "ihw") {

+    #   # run MRihw

+    #   padj = MRihw(obj, p, IHWcov, alpha)

+    # } else {

+    #   padj = p.adjust(p, method = adjustMethod) # use classic pvalue adjusment methods

+    # }

+    

+    # adding 'ihw' as pvalue adjustment method

+    if (adjustMethod == "ihw-ubiquity" | adjustMethod == "ihw-abundance") {

+      # use IHW to adjust pvalues

+      padj = MRihw(obj, p, adjustMethod, alpha)

+    } else {

+      # use classic pvalue adjusment method

+      padj = p.adjust(p, method = adjustMethod)

+    }

+    

     if(group==0){

         srt = order(abs(tb[,by]),decreasing=TRUE)

     } else if(group==1){

@@ -91,7 +114,7 @@ MRcoefs<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,

     valid = 1:length(padj);

     if(counts>0){

-        np=rowSums(obj$counts);

+        np=rowSums(obj@counts);

         valid = intersect(valid,which(np>=counts));

     }

     srt = srt[which(srt%in%valid)][1:min(number,nrow(tb))];

new file mode 100644

@@ -0,0 +1,66 @@

+#' MRihw runs IHW within a MRcoefs() call

+#' 

+#' Function used in MRcoefs() when "IHW" is set as the p value adjustment method

+#' 

+#' @rdname MRihw

+#' @param obj Either a fitFeatureModelResults or fitZigResults object

+#' @param ... other parameters

+#' 

+setGeneric("MRihw", function(obj, ...){standardGeneric("MRihw")})

+

+#' MRihw runs IHW within a MRcoefs() call

+#' 

+#' Function used in MRcoefs() when "IHW" is set as the p value adjustment method

+#' 

+#' @rdname MRihw-fitFeatureModelResults

+#' @param obj Either a fitFeatureModelResults or fitZigResults object

+#' @param p a vector of pvalues extracted from obj

+#' @param IHWcov Value specifying which covariate to use in IHW pvalue adjustment. For obj

+#' of class \code{\link{fitFeatureModelResults}}, options are "median" and "nnz" (number of non-zero 

+#' features per row). For obj of class \code{\link{fitZigResults}}, options are "Amean" and "nnz".

+#' @param alpha pvalue significance level specified for IHW call. Default is 0.1

+#' 

+setMethod("MRihw", signature = "fitFeatureModelResults", function(obj, p, adjustMethod, alpha){

+  if (adjustMethod == "ihw-ubiquity") {

+    # set covariate to be num of non-zero elements per row

+    p <- obj@pvalues

+    covariate <- rowSums(obj@counts != 0)

+    ihwRes <- ihw(p, covariate, alpha)

+    padj <- ihwRes@df$adj_pvalue 

+  }

+  if (adjustMethod == "ihw-abundance"){

+    # use feature median count as covariate

+    covariate <- rowMedians(obj@counts)

+    ihwRes <- ihw(p, covariate, alpha)

+    padj <- ihwRes@df$adj_pvalue

+  }

+  padj

+})

+

+#' MRihw runs IHW within a MRcoefs() call

+#' 

+#' Function used in MRcoefs() when "IHW" is set as the p value adjustment method

+#' 

+#' @rdname MRihw-fitZigResults

+#' @param obj Either a fitFeatureModelResults or fitZigResults object

+#' @param p a vector of pvalues extracted from obj

+#' @param IHWcov Value specifying which covariate to use in IHW pvalue adjustment. For obj

+#' of class \code{\link{fitFeatureModelResults}}, options are "median" and "nnz" (number of non-zero 

+#' features per row). For obj of class \code{\link{fitZigResults}}, options are "Amean" and "nnz".

+#' @param alpha pvalue significance level specified for IHW call. Default is 0.1

+#' 

+setMethod("MRihw", signature = "fitZigResults", function(obj, p, adjustMethod, alpha){

+  if (adjustMethod == "ihw-ubiquity"){

+    #use number of non-zero features per row as the covariate in ihw() call

+    covariate <- rowSums(obj@counts != 0)

+    ihwRes <- ihw(p, covariate, alpha)

+    padj <- ihwRes@df$adj_pvalue

+  }

+  if (adjustMethod == "ihw-abundance"){

+    # use Amean as covariate

+    covariate <- obj@eb$Amean

+    ihwRes <- ihw(p, covariate, alpha)

+    padj <- ihwRes@df$adj_pvalue

+  }

+  padj

+})

\ No newline at end of file

@@ -46,21 +46,21 @@

 #' fit = fitFeatureModel(obj = lungTrim,mod=mod)

 #' head(MRtable(fit))

 #'

-MRtable<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,

+MRtable<-function(obj,by=2,coef=NULL,number=10,taxa=obj@taxa,

     uniqueNames=FALSE,adjustMethod="fdr",group=0,eff=0,numberEff=FALSE,ncounts=0,file=NULL){

-    if(length(grep("fitFeatureModel",obj$call))){

-        groups = factor(obj$design[,by])

+    if(length(grep("fitFeatureModel",obj@call))){

+        groups = factor(obj@design[,by])

         by = "logFC"; coef = 1:2;

-        tb = data.frame(logFC=obj$fitZeroLogNormal$logFC,se=obj$fitZeroLogNormal$se)

-        p  = obj$pvalues

+        tb = data.frame(logFC=obj@fitZeroLogNormal$logFC,se=obj@fitZeroLogNormal$se)

+        p  = obj@pvalues

     } else {

-        tb = obj$fit$coefficients

+        tb = obj@fit$coefficients

         if(is.null(coef)){

             coef = 1:ncol(tb)

         }

-        p=obj$eb$p.value[,by]

-        groups = factor(obj$fit$design[,by])

+        p=obj@eb$p.value[,by]

+        groups = factor(obj@fit$design[,by])

         if(eff>0){

             effectiveSamples = calculateEffectiveSamples(obj)

             if(numberEff == FALSE){

@@ -79,7 +79,7 @@ MRtable<-function(obj,by=2,coef=NULL,number=10,taxa=obj$taxa,

         }

     }

     padj = p.adjust(p,method=adjustMethod)

-    cnts = obj$counts

+    cnts = obj@counts

     posIndices = cnts>0

     np0 = rowSums(posIndices[,groups==0])

@@ -1,32 +1,32 @@

 setClass("MRexperiment", contains=c("eSet"), representation=representation(expSummary = "list"),prototype = prototype( new( "VersionedBiobase",versions = c(classVersion("eSet"),MRexperiment = "1.0.0" ))))

 setMethod("[", "MRexperiment", function (x, i, j, ..., drop = FALSE) {

-        obj= callNextMethod()

-        if(!missing(j)){

-            obj@expSummary = new("list",expSummary=as(expSummary(x)[j,1:2,...,drop=drop],"AnnotatedDataFrame"),cumNormStat=x@expSummary$cumNormStat)

-            if(length(pData(obj))>0){

-              for(i in 1:length(pData(obj))){

-                if(is.factor(pData(obj)[,i])){

-                  pData(obj)[,i] = factor(pData(obj)[,i])

-                } else {

-                  pData(obj)[,i] = pData(obj)[,i]

-                }

-              }

-            }

+  obj= callNextMethod()

+  if(!missing(j)){

+    obj@expSummary = new("list",expSummary=as(expSummary(x)[j,1:2,...,drop=drop],"AnnotatedDataFrame"),cumNormStat=x@expSummary$cumNormStat)

+    if(length(pData(obj))>0){

+      for(i in 1:length(pData(obj))){

+        if(is.factor(pData(obj)[,i])){

+          pData(obj)[,i] = factor(pData(obj)[,i])

+        } else {

+          pData(obj)[,i] = pData(obj)[,i]

         }

-        obj

+      }

+    }

+  }

+  obj

 })

 setMethod("colSums", signature ="MRexperiment", function (x, ...) {

-    callNextMethod(MRcounts(x),...)

+  callNextMethod(MRcounts(x),...)

 })

 setMethod("rowSums", signature="MRexperiment", function (x, ...) {

-    callNextMethod(MRcounts(x),...)

+  callNextMethod(MRcounts(x),...)

 })

 setMethod("rowMeans", signature="MRexperiment", function (x, ...) {

-    callNextMethod(MRcounts(x),...)

+  callNextMethod(MRcounts(x),...)

 })

 setMethod("colMeans", signature="MRexperiment", function (x, ...) {

-    callNextMethod(MRcounts(x),...)

+  callNextMethod(MRcounts(x),...)

 })

 #' Access the normalization factors in a MRexperiment object

@@ -48,11 +48,11 @@ setGeneric("normFactors",function(object){standardGeneric("normFactors")})

 setGeneric("normFactors<-",function(object,value){standardGeneric("normFactors<-")})

 setMethod("normFactors", signature="MRexperiment",function(object) {

-   nf <- expSummary(object)$normFactors

-   nf <- unlist(nf)

-   names(nf) <- sampleNames(object)

-   nf

- })

+  nf <- expSummary(object)$normFactors

+  nf <- unlist(nf)

+  names(nf) <- sampleNames(object)

+  nf

+})

 #' Replace the normalization factors in a MRexperiment object

 #'

@@ -72,11 +72,11 @@ setMethod("normFactors", signature="MRexperiment",function(object) {

 #' head(normFactors(lungData)<- rnorm(1))

 #'

 setReplaceMethod("normFactors", signature=c(object="MRexperiment", value="numeric"),

-  function( object, value ) {

-   pData(object@expSummary$expSummary)$normFactors <- value

-   validObject( object )

-   object

-})

+                 function( object, value ) {

+                   pData(object@expSummary$expSummary)$normFactors <- value

+                   validObject( object )

+                   object

+                 })

 #' Access sample depth of coverage from MRexperiment object

 #'

@@ -98,11 +98,11 @@ setGeneric("libSize",function(object){standardGeneric("libSize")})

 setGeneric("libSize<-",function(object,value){standardGeneric("libSize<-")})

 setMethod("libSize", signature="MRexperiment",function(object) {

-   ls <- expSummary(object)$libSize

-   ls <- unlist(ls)

-   names(ls) <- sampleNames(object)

-   ls

- })

+  ls <- expSummary(object)$libSize

+  ls <- unlist(ls)

+  names(ls) <- sampleNames(object)

+  ls

+})

 #' Replace the library sizes in a MRexperiment object

 #'

@@ -122,11 +122,54 @@ setMethod("libSize", signature="MRexperiment",function(object) {

 #' head(libSize(lungData)<- rnorm(1))

 #'

 setReplaceMethod("libSize", signature=c(object="MRexperiment", value="numeric"),

-  function( object, value ) {

-   pData(object@expSummary$expSummary)$libSize <- value

-   validObject( object )

-   object

-})

+                 function( object, value ) {

+                   pData(object@expSummary$expSummary)$libSize <- value

+                   validObject( object )

+                   object

+                 })

+

+#' Class "fitZigResults" -- a formal class for storing results from a fitZig call

+#' 

+#' This class contains all of the same information expected from a fitZig call, 

+#' but it is defined in the S4 style as opposed to being stored as a list.

+#' 

+#' @slot call the call made to fitZig

+#' @slot fit 'MLArrayLM' Limma object of the weighted fit

+#' @slot countResiduals standardized residuals of the fit

+#' @slot z matrix of the posterior probabilities. It is defined as $z_ij = pr(delta_ij=1 | data)$

+#' @slot zUsed used in \code{\link{getZ}}

+#' @slot eb output of eBayes, moderated t-statistics, moderated F-statistics, etc

+#' @slot taxa vector of the taxa names

+#' @slot counts the original count matrix input

+#' @slot zeroMod the zero model matrix

+#' @slot zeroCoef the zero model fitted results

+#' @slot stillActive convergence

+#' @slot stillActiveNLL nll at convergence

+#' @slot dupcor correlation of duplicates

+#' @exportClass

+#' 

+setClass("fitZigResults", 

+         slots = c(fit = "list", countResiduals = "matrix", z = "matrix", zUsed = "ANY", 

+                   eb = "MArrayLM", zeroMod = "matrix", stillActive = "logical", stillActiveNLL = "numeric", 

+                   zeroCoef = "list", dupcor = "ANY", call = "call", taxa = "character", counts = "matrix"))

+

+#' Class "fitFeatureModelResults" -- a formal class for storing results from a fitFeatureModel call

+#' 

+#' This class contains all of the same information expected from a fitFeatureModel call, 

+#' but it is defined in the S4 style as opposed to being stored as a list.

+#' 

+#' @slot call  the call made to fitFeatureModel

+#' @slot fitZeroLogNormal  list of parameter estimates for the zero-inflated log normal model

+#' @slot design  model matrix

+#' @slot taxa  taxa names

+#' @slot counts  count matrix

+#' @slot pvalues  calculated p-values

+#' @slot permuttedFits  permutted z-score estimates under the null

+#' @exportClass

+#' 

+setClass("fitFeatureModelResults",

+         slots = c(call = "call", fitZeroLogNormal = "list", design = "matrix", taxa = "character",

+                   counts = "matrix", pvalues = "numeric", permuttedFits = "ANY"))

 #' Create a MRexperiment object

 #' 

@@ -147,6 +190,7 @@ setReplaceMethod("libSize", signature=c(object="MRexperiment", value="numeric"),

 #' @param normFactors normFactors, the normalization factors used in either the

 #' model or as scaling factors of sample counts for each particular sample.

 #' @return an object of class MRexperiment

+#' @export

 #' @author Joseph N Paulson

 #' @examples

 #' 

@@ -154,37 +198,37 @@ setReplaceMethod("libSize", signature=c(object="MRexperiment", value="numeric"),

 #' obj <- newMRexperiment(cnts)

 #' 

 newMRexperiment <- function(counts, phenoData=NULL, featureData=NULL,libSize=NULL, normFactors=NULL) {

-    counts= as.matrix(counts)

-

-    if( is.null( featureData ) ){

-      featureData <- annotatedDataFrameFrom(counts, byrow=TRUE)

-    }

-    if( is.null( phenoData ) ){

-      phenoData   <- annotatedDataFrameFrom(counts, byrow=FALSE)

-    }

-    if( is.null( libSize ) ){

-      libSize <- as.matrix(colSums(counts))

-      rownames(libSize) = colnames(counts)

-    }

-    if( is.null( normFactors ) ){

-      normFactors <- as.matrix(rep( NA_real_, length(libSize) ))

-      rownames(normFactors) = rownames(libSize)

-    }

-

-    obj <-new("MRexperiment", assayData = assayDataNew("environment",counts=counts),phenoData = phenoData,featureData = featureData ,expSummary = new("list",expSummary=annotatedDataFrameFrom(counts,byrow=FALSE),cumNormStat=NULL))

-    obj@expSummary$expSummary$libSize = libSize;

-    obj@expSummary$expSummary$normFactors=normFactors;

-    validObject(obj)

-    obj

+  counts= as.matrix(counts)

+  

+  if( is.null( featureData ) ){

+    featureData <- annotatedDataFrameFrom(counts, byrow=TRUE)

+  }

+  if( is.null( phenoData ) ){

+    phenoData   <- annotatedDataFrameFrom(counts, byrow=FALSE)

+  }

+  if( is.null( libSize ) ){

+    libSize <- as.matrix(colSums(counts))

+    rownames(libSize) = colnames(counts)

+  }

+  if( is.null( normFactors ) ){

+    normFactors <- as.matrix(rep( NA_real_, length(libSize) ))

+    rownames(normFactors) = rownames(libSize)

+  }

+  

+  obj <-new("MRexperiment", assayData = assayDataNew("environment",counts=counts),phenoData = phenoData,featureData = featureData ,expSummary = new("list",expSummary=annotatedDataFrameFrom(counts,byrow=FALSE),cumNormStat=NULL))

+  obj@expSummary$expSummary$libSize = libSize;

+  obj@expSummary$expSummary$normFactors=normFactors;

+  validObject(obj)

+  obj

 }

 setValidity( "MRexperiment", function( object ) {

-    if( is.null(assayData(object)$counts))

-        return( "There are no counts!" )

-#    if( ncol(MRcounts(object)) != length(normFactors(object)))

-#        return( "Experiment summary got hacked!" )

-#    if( ncol(MRcounts(object)) != length(libSize(object)))

-#        return( "Experiment summary got hacked!" )

-    TRUE

+  if( is.null(assayData(object)$counts))

+    return( "There are no counts!" )

+  #    if( ncol(MRcounts(object)) != length(normFactors(object)))

+  #        return( "Experiment summary got hacked!" )

+  #    if( ncol(MRcounts(object)) != length(libSize(object)))

+  #        return( "Experiment summary got hacked!" )

+  TRUE

 } )

 #' Accessor for the counts slot of a MRexperiment object

 #' 

@@ -208,22 +252,22 @@ setValidity( "MRexperiment", function( object ) {

 #' head(MRcounts(lungData))

 #' 

 MRcounts <- function(obj,norm=FALSE,log=FALSE,sl=1000) {

-   stopifnot( is( obj, "MRexperiment" ) )

-   if(!norm){

+  stopifnot( is( obj, "MRexperiment" ) )

+  if(!norm){

     x=assayData(obj)[["counts"]]

-   }

-   else{

+  }

+  else{

     if(any(is.na(normFactors(obj)))){

       x=cumNormMat(obj,sl=sl)

     } else{

       x=sweep(assayData(obj)[["counts"]],2,as.vector(unlist(normFactors(obj)))/sl,"/")

     }

-   }

-   if(!log){

+  }

+  if(!log){

     return(x)

-   } else{

+  } else{

     return(log2(x+1))

-   }

+  }

 }

 #' Access the posterior probabilities that results from analysis

 #' 

@@ -248,14 +292,14 @@ MRcounts <- function(obj,norm=FALSE,log=FALSE,sl=1000) {

 #' lungTrim = lungTrim[-k,]

 #' smokingStatus = pData(lungTrim)$SmokingStatus

 #' mod = model.matrix(~smokingStatus)

-#' # The maxit is not meant to be 1 - this is for demonstration/speed

+#' # The maxit is not meant to be 1 -- this is for demonstration/speed

 #' settings = zigControl(maxit=1,verbose=FALSE)

 #' fit = fitZig(obj = lungTrim,mod=mod,control=settings)

 #' head(posteriorProbs(lungTrim))

 #' 

 posteriorProbs <- function( obj ) {

-   stopifnot( is( obj, "MRexperiment" ) )

-   assayData(obj)[["z"]]

+  stopifnot( is( obj, "MRexperiment" ) )

+  assayData(obj)[["z"]]

 }

 #' Access MRexperiment object experiment data

@@ -76,7 +76,13 @@ fitFeatureModel<-function(obj,mod,coef=2,B=1,szero=FALSE,spos=TRUE){

   } else {

     pvals = 2*(1-pnorm(abs(zscore)))

   }

-  res = list(call=match.call(),fitZeroLogNormal=fitzeroln,design=mmCount,

-    taxa=rownames(obj),counts=MRcounts(obj),pvalues=pvals,permuttedFits=permuttedFits)

+  # old way of creating results object

+  # res = list(call=match.call(),fitZeroLogNormal=fitzeroln,design=mmCount,

+  #   taxa=rownames(obj),counts=MRcounts(obj),pvalues=pvals,permuttedFits=permuttedFits)

+  

+  # new way with defined results class

+  res = new("fitFeatureModelResults", call = match.call(), fitZeroLogNormal=fitzeroln, 

+            design = mmCount, taxa = rownames(obj), counts = MRcounts(obj), 

+            pvalues = pvals, permuttedFits = permuttedFits)

   res

 }

\ No newline at end of file

@@ -100,7 +100,16 @@ fitZig <- function(obj,

 	dat$zUsed <- NULL

 	dat <- c(dat, list(call=match.call(),taxa=rownames(obj),counts=y))

-	dat

+	

+	# old way of outputting results with list

+	# dat <- c(dat, list(call=match.call(),taxa=rownames(obj),counts=y))

+	

+	# new output with defined results class

+	dat <- new("fitZigResults", fit=dat$fit, countResiduals=dat$countResiduals,

+	           z=dat$z, zUsed=dat$zUsed, eb=dat$eb, zeroMod=dat$zeroMod, stillActive=dat$stillActive,

+	           stillActiveNLL=dat$stillActiveNLL, zeroCoef=dat$zeroCoef, dupcor=dat$dupcor, call = dat$call,

+	           taxa = rownames(obj), counts = dat$counts)

+ 	dat

 }

 .do_fitZig <- function(y, 

@@ -181,7 +190,7 @@ fitZig <- function(obj,

   eb <- limma::eBayes(fit$fit)

   dat <- list(fit=fit$fit, countResiduals=fit$residuals,

-              z=z, zUsed=zUsed, eb=eb, zeroMod=zero_model_matrix, stillActive=stillActive, 

+              z=z, zUsed=zUsed, eb=eb, zeroMod=zero_model_matrix, stillActive=stillActive,

               stillActiveNLL=stillActiveNLL, zeroCoef=zeroCoef, dupcor=dupcor)

   dat

 }

@@ -4,9 +4,9 @@

 \alias{MRcoefs}

 \title{Table of top-ranked features from fitZig or fitFeatureModel}

 \usage{

-MRcoefs(obj, by = 2, coef = NULL, number = 10, taxa = obj$taxa,

-  uniqueNames = FALSE, adjustMethod = "fdr", group = 0, eff = 0,

-  numberEff = FALSE, counts = 0, file = NULL)

+MRcoefs(obj, by = 2, coef = NULL, number = 10, taxa = obj@taxa,

+  uniqueNames = FALSE, adjustMethod = "fdr", alpha = 0.1,

+  group = 0, eff = 0, numberEff = FALSE, counts = 0, file = NULL)

 }

 \arguments{

 \item{obj}{Output of fitFeatureModel or fitZig.}

@@ -24,9 +24,12 @@ or contrast of the linear model to display.}

 \item{uniqueNames}{Number the various taxa.}

 \item{adjustMethod}{Method to adjust p-values by. Default is "FDR". Options

-include "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

+include "IHW", "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr",

 "none". See \code{\link{p.adjust}} for more details.}

+\item{alpha}{Value for p-value significance threshold when running IHW. 

+The default is set to 0.1}

+

 \item{group}{One of five choices, 0,1,2,3,4. 0: the sort is ordered by a

 decreasing absolute value coefficient fit. 1: the sort is ordered by the raw

 coefficient fit in decreasing order. 2: the sort is ordered by the raw

@@ -40,6 +43,11 @@ of the coefficient fit in increasing order. 4: no sorting.}

 \item{counts}{Filter features to have at least 'counts' counts.}

 \item{file}{Name of output file, including location, to save the table.}

+

+\item{IHWcov}{Character value specifying which covariate to use when adjusting pvalues

+using IHW. Options include: "nnz" (number of non-zero elements per feature), 

+"median" (median abundance value per feature), "Amean" (adjusted mean, used for a 

+fitZigResults obj)}

 }

 \value{

 Table of the top-ranked features determined by the linear fit's

new file mode 100644

@@ -0,0 +1,23 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/MRihw.R

+\docType{methods}

+\name{MRihw,fitFeatureModelResults-method}

+\alias{MRihw,fitFeatureModelResults-method}

+\title{MRihw runs IHW within a MRcoefs() call}

+\usage{

+\S4method{MRihw}{fitFeatureModelResults}(obj, p, adjustMethod, alpha)

+}

+\arguments{

+\item{obj}{Either a fitFeatureModelResults or fitZigResults object}

+

+\item{p}{a vector of pvalues extracted from obj}

+

+\item{alpha}{pvalue significance level specified for IHW call. Default is 0.1}

+

+\item{IHWcov}{Value specifying which covariate to use in IHW pvalue adjustment. For obj

+of class \code{\link{fitFeatureModelResults}}, options are "median" and "nnz" (number of non-zero 

+features per row). For obj of class \code{\link{fitZigResults}}, options are "Amean" and "nnz".}

+}

+\description{

+Function used in MRcoefs() when "IHW" is set as the p value adjustment method

+}

new file mode 100644

@@ -0,0 +1,23 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/MRihw.R

+\docType{methods}

+\name{MRihw,fitZigResults-method}

+\alias{MRihw,fitZigResults-method}

+\title{MRihw runs IHW within a MRcoefs() call}

+\usage{

+\S4method{MRihw}{fitZigResults}(obj, p, adjustMethod, alpha)

+}

+\arguments{

+\item{obj}{Either a fitFeatureModelResults or fitZigResults object}

+

+\item{p}{a vector of pvalues extracted from obj}

+

+\item{alpha}{pvalue significance level specified for IHW call. Default is 0.1}

+

+\item{IHWcov}{Value specifying which covariate to use in IHW pvalue adjustment. For obj

+of class \code{\link{fitFeatureModelResults}}, options are "median" and "nnz" (number of non-zero 

+features per row). For obj of class \code{\link{fitZigResults}}, options are "Amean" and "nnz".}

+}

+\description{

+Function used in MRcoefs() when "IHW" is set as the p value adjustment method

+}

new file mode 100644

@@ -0,0 +1,16 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/MRihw.R

+\name{MRihw}

+\alias{MRihw}

+\title{MRihw runs IHW within a MRcoefs() call}

+\usage{

+MRihw(obj, ...)

+}

+\arguments{

+\item{obj}{Either a fitFeatureModelResults or fitZigResults object}

+

+\item{...}{other parameters}

+}

+\description{

+Function used in MRcoefs() when "IHW" is set as the p value adjustment method

+}

@@ -5,7 +5,7 @@

 \title{Table of top microbial marker gene from linear model fit including sequence

 information}

 \usage{

-MRtable(obj, by = 2, coef = NULL, number = 10, taxa = obj$taxa,

+MRtable(obj, by = 2, coef = NULL, number = 10, taxa = obj@taxa,

   uniqueNames = FALSE, adjustMethod = "fdr", group = 0, eff = 0,

   numberEff = FALSE, ncounts = 0, file = NULL)

 }

new file mode 100644

@@ -0,0 +1,28 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/allClasses.R

+\docType{class}

+\name{fitFeatureModelResults-class}

+\alias{fitFeatureModelResults-class}

+\title{Class "fitFeatureModelResults" -- a formal class for storing results from a fitFeatureModel call}

+\description{

+This class contains all of the same information expected from a fitFeatureModel call, 

+but it is defined in the S4 style as opposed to being stored as a list.

+}

+\section{Slots}{

+

+\describe{

+\item{\code{call}}{the call made to fitFeatureModel}

+

+\item{\code{fitZeroLogNormal}}{list of parameter estimates for the zero-inflated log normal model}

+

+\item{\code{design}}{model matrix}

+

+\item{\code{taxa}}{taxa names}

+

+\item{\code{counts}}{count matrix}

+

+\item{\code{pvalues}}{calculated p-values}

+

+\item{\code{permuttedFits}}{permutted z-score estimates under the null}

+}}

+

new file mode 100644

@@ -0,0 +1,40 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/allClasses.R

+\docType{class}

+\name{fitZigResults-class}

+\alias{fitZigResults-class}

+\title{Class "fitZigResults" -- a formal class for storing results from a fitZig call}

+\description{

+This class contains all of the same information expected from a fitZig call, 

+but it is defined in the S4 style as opposed to being stored as a list.

+}

+\section{Slots}{

+

+\describe{

+\item{\code{call}}{the call made to fitZig}

+

+\item{\code{fit}}{'MLArrayLM' Limma object of the weighted fit}

+

+\item{\code{countResiduals}}{standardized residuals of the fit}

+

+\item{\code{z}}{matrix of the posterior probabilities. It is defined as $z_ij = pr(delta_ij=1 | data)$}

+

+\item{\code{zUsed}}{used in \code{\link{getZ}}}

+

+\item{\code{eb}}{output of eBayes, moderated t-statistics, moderated F-statistics, etc}

+

+\item{\code{taxa}}{vector of the taxa names}

+

+\item{\code{counts}}{the original count matrix input}

+

+\item{\code{zeroMod}}{the zero model matrix}

+

+\item{\code{zeroCoef}}{the zero model fitted results}

+

+\item{\code{stillActive}}{convergence}

+

+\item{\code{stillActiveNLL}}{nll at convergence}

+

+\item{\code{dupcor}}{correlation of duplicates}

+}}

+

@@ -29,7 +29,7 @@ lungTrim = cumNorm(lungTrim)

 lungTrim = lungTrim[-k,]

 smokingStatus = pData(lungTrim)$SmokingStatus

 mod = model.matrix(~smokingStatus)

-# The maxit is not meant to be 1 - this is for demonstration/speed

+# The maxit is not meant to be 1 -- this is for demonstration/speed

 settings = zigControl(maxit=1,verbose=FALSE)

 fit = fitZig(obj = lungTrim,mod=mod,control=settings)

 head(posteriorProbs(lungTrim))