version 1.0.0

cran · Aug 31, 2006 · 7154dda · 7154dda
commit 7154dda
Show file tree

Hide file tree

Showing 14 changed files with 1,175 additions and 0 deletions.
diff --git a/CHANGES b/CHANGES
@@ -0,0 +1,12 @@
+
+
+
+		CHANGES IN fdrtool VERSION 1.0.0
+
+
+This is the first standalone release (31 August 2006).
+
+This package implements the "shrinkage t" statistics described
+in Opgen-Rhein and Strimmer (2006).  It also offers a convenient
+interface to a number of other regularized t-type statistics often
+used in high-dimensional case-control studies.    
diff --git a/COPYING b/COPYING
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -0,0 +1,15 @@
+Package: st
+Version: 1.0.0
+Date: 2006-31-08
+Title: Shrinkage t Statistic
+Author: Rainer Opgen-Rhein and Korbinian Strimmer.
+Maintainer: Korbinian Strimmer <korbinian.strimmer@lmu.de>
+Depends: R (>= 2.0.0), corpcor (>= 1.4.3)
+Suggests:
+Description: This package implements the "shrinkage t" statistics 
+   described in Opgen-Rhein and Strimmer (2006).  It also offers
+   a convenient interface to a number of other regularized t-type 
+   statistics often used in high-dimensional case-control studies.    
+License: GPL version 2 or newer
+URL: http://www.strimmerlab.org/software/st/
+Packaged: Thu Aug 31 14:40:24 2006; strimmer
diff --git a/R/efront.R b/R/efront.R
@@ -0,0 +1,60 @@
+### efront.R  (2006-08-30)
+###
+###    Efron t Statistic (2001)
+###
+### Copyright 2006 Rainer Opgen-Rhein and Korbinian Strimmer
+###
+###
+### This file is part of the `st' library for R and related languages.
+### It is made available under the terms of the GNU General Public
+### License, version 2, or at your option, any later version,
+### incorporated herein by reference.
+### 
+### This program is distributed in the hope that it will be
+### useful, but WITHOUT ANY WARRANTY; without even the implied
+### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+### PURPOSE.  See the GNU General Public License for more
+### details.
+### 
+### You should have received a copy of the GNU General Public
+### License along with this program; if not, write to the Free
+### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+### MA 02111-1307, USA
+
+
+efront.stat = function (X, L, verbose=TRUE)
+{
+  FUN = efront.fun(L=L, verbose=verbose)
+  score = FUN(X)
+
+  return( score )
+}
+
+efront.fun <- function (L, verbose=TRUE)
+{
+    if (missing(L))
+      stop("Please specify to which group (1 or 2) each sample is assigned!")
+
+    function(X)
+    {
+      tmp = pvt.group.moments(X, L, variances=TRUE)
+
+      # differences between the two groups
+      diff = tmp$mu1-tmp$mu2
+
+      # standard error of diff
+      n1 = sum(L==1); n2 = sum(L==2)      
+      sd = sqrt( (1/n1 + 1/n2)*tmp$v.pooled )
+
+      # tuning parameter
+      a0 <- quantile(sd, probs=c(0.9))
+
+      if (verbose) cat("Fudge factor a0 =", a0, "\n")
+
+      # t statistic
+      t = diff/(sd+a0)
+
+      return(t)
+    }
+}
+
diff --git a/R/modt.R b/R/modt.R
@@ -0,0 +1,52 @@
+### modt.R  (2006-08-30)
+###
+###    Moderated t Statistic
+###
+### Copyright 2006 Rainer Opgen-Rhein and Korbinian Strimmer
+###
+###
+### This file is part of the `st' library for R and related languages.
+### It is made available under the terms of the GNU General Public
+### License, version 2, or at your option, any later version,
+### incorporated herein by reference.
+### 
+### This program is distributed in the hope that it will be
+### useful, but WITHOUT ANY WARRANTY; without even the implied
+### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+### PURPOSE.  See the GNU General Public License for more
+### details.
+### 
+### You should have received a copy of the GNU General Public
+### License along with this program; if not, write to the Free
+### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+### MA 02111-1307, USA
+
+
+# Note: these function2 require the "limma" library
+
+
+modt.stat = function (X, L)
+{
+  FUN = modt.fun(L=L)
+  score = FUN(X)
+
+  return( score )
+}
+
+modt.fun <- function (L)
+{
+    require("limma")
+
+    if (missing(L))
+      stop("Please specify to which group (1 or 2) each sample is assigned!")
+
+    function(X)
+    {
+      d <- cbind(rep(1, length(L)), L)
+      fit <- lmFit(t(X), design=d)
+      eb.out <- ebayes(fit)
+      modt <- -eb.out$t[,2]
+
+      return(modt) 
+    } 
+}
diff --git a/R/sam.R b/R/sam.R
@@ -0,0 +1,49 @@
+### sam.R  (2006-08-30)
+###
+###    SAM t Statistic
+###
+### Copyright 2006 Rainer Opgen-Rhein and Korbinian Strimmer
+###
+###
+### This file is part of the `st' library for R and related languages.
+### It is made available under the terms of the GNU General Public
+### License, version 2, or at your option, any later version,
+### incorporated herein by reference.
+### 
+### This program is distributed in the hope that it will be
+### useful, but WITHOUT ANY WARRANTY; without even the implied
+### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+### PURPOSE.  See the GNU General Public License for more
+### details.
+### 
+### You should have received a copy of the GNU General Public
+### License along with this program; if not, write to the Free
+### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+### MA 02111-1307, USA
+
+
+# Note: these function2 require the "samr" library
+
+sam.stat = function (X, L)
+{
+  FUN = sam.fun(L=L)
+  score = FUN(X)
+
+  return( score )
+}
+
+sam.fun <- function(L)
+{
+   require("samr")
+
+   if (missing(L))
+      stop("Please specify to which group (1 or 2) each sample is assigned!")
+
+   function(X)
+    {
+      dd = list(x=t(X),y=L, logged2=TRUE)
+      out = samr(dd, resp.type="Two class unpaired", nperms=1)
+
+      return(out$tt) # SAM test statistic
+    } 
+}
diff --git a/R/samL1.R b/R/samL1.R
@@ -0,0 +1,126 @@
+### samL1.R  (2006-08-30)
+###
+###    Wu (2005)Improved SAM Statistic
+###
+### Copyright 2006 Rainer Opgen-Rhein and Korbinian Strimmer
+###
+### This function is based in part on R code provided by Baolin Wu.
+###
+###
+### This file is part of the `st' library for R and related languages.
+### It is made available under the terms of the GNU General Public
+### License, version 2, or at your option, any later version,
+### incorporated herein by reference.
+### 
+### This program is distributed in the hope that it will be
+### useful, but WITHOUT ANY WARRANTY; without even the implied
+### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+### PURPOSE.  See the GNU General Public License for more
+### details.
+### 
+### You should have received a copy of the GNU General Public
+### License along with this program; if not, write to the Free
+### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+### MA 02111-1307, USA
+
+
+#### estimate Wu t-statistic (2005)
+
+# B. Wu. 2005. Differential gene expression detection using penalized 
+# linear regression model: the improved SAM statistics. 
+# Bioinformatics, 21: 1565-1571
+
+samL1.stat = function (X, L, method=c("lowess", "cor"), 
+   plot=FALSE, verbose=TRUE)
+{
+  FUN = samL1.fun(L=L, method=method, verbose=verbose, plot=plot)
+  score = FUN(X)
+
+  return( score )
+}
+
+samL1.fun <- function (L, method=c("lowess", "cor"), 
+   plot=FALSE, verbose=TRUE)
+{
+    method = match.arg(method)
+
+    if (missing(L))
+      stop("Please specify to which group (1 or 2) each sample is assigned!")
+
+    function(X)
+    {
+
+      tmp = pvt.group.moments(X, L, variances=TRUE)
+
+      # differences between the two groups
+      diff = tmp$mu1-tmp$mu2
+
+      # variance of diff
+      n1 = sum(L==1); n2 = sum(L==2)      
+      v.diff = (1/n1 + 1/n2)*tmp$v.pooled
+      sd = sqrt(v.diff)
+
+      lambda = pvt.samL1.get.lambda(diff, sd, method=method, verbose=verbose, plot=plot)
+
+      # penalized t statistic
+      nom = ifelse(abs(diff)>lambda, diff-sign(diff)*lambda, 0)
+      den = sqrt(v.diff + lambda^2/(n1+n2-2))
+
+      tL1 = nom/den
+
+      return(tL1)
+    }
+}
+
+## internal function
+
+pvt.samL1.get.lambda = function(di, si, method=c("cor", "lowess"), verbose=TRUE, plot=FALSE)
+{
+  if (method == "lowess")
+  {
+    if (verbose) cat("Optimizing lambda (lowess)...\n")
+
+    Lambda = seq(5, 100, length=1e2)
+    rL = sapply(Lambda, function(lambda){
+      mu12 = ifelse( abs(di)>lambda, di-sign(di)*lambda, 0)
+      tmp = mu12/sqrt(2*si^2+lambda^2/3)
+      a = lowess(si, tmp, f=2/3)$y; 
+      mean((tmp[order(si)]-a)^2)/var(tmp) ## SSE/SSTO
+    })
+    Lam.Opt = Lambda[order(-rL)[1]] 
+
+    if (plot==TRUE)
+    {
+      plot(Lambda, sqrt(rL), pch=20, xlab=expression(lambda), 
+         ylab=expression(sqrt(SSE/SSTO)))
+      abline(v=Lam.Opt, col=2, lty=2)
+    }
+  }
+
+  if (method == "cor")
+  {
+    if (verbose) cat("Optimizing lambda (cor)...\n")
+
+    Lambda = seq(5, 100, length=1e3)
+    crL = sapply(Lambda, function(lambda){
+      mu12 = ifelse( abs(di)>lambda, di-sign(di)*lambda, 0)
+      tmp = mu12/sqrt(2*si^2+lambda^2/3)
+      idm = order(-si)[1] ## remove the biggest si (outlier)
+      return( cor(tmp[-idm],si[-idm]) )
+    })
+    Lam.Opt = Lambda[order(abs(crL))[1]] ## 27.14
+
+    if (plot==TRUE)
+    {
+      plot(Lambda, abs(crL), pch=20, xlab=expression(lambda), ylab="|R|")
+      abline(v=Lam.Opt, col=2, lty=2)
+    }
+
+  }
+
+  if (verbose) cat("Estimated lamba: ", Lam.Opt, "\n")
+
+  return(Lam.Opt)
+}
+
+