version 0.7

cran · Dec 28, 2013 · 7ab38cd · 7ab38cd
1 parent f13f92d
commit 7ab38cd
Show file tree

Hide file tree

Showing 5 changed files with 45 additions and 231 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -2,9 +2,9 @@ Package: rgabriel
 Type: Package
 Title: Gabriel Multiple Comparison Test and Plot the Confidence
         Interval on Barplot
-Version: 0.4
-Date: 2013-12-15
-Author: Miao YU
+Version: 0.7
+Date: 2013-12-28
+Author: Yihui XIE, Miao YU
 Maintainer: Miao YU <yufreecas@gmail.com>
 Description: This package was created to analyze multi-level one-way
         experimental designs. It is designed to handle vectorized
@@ -17,7 +17,7 @@ Description: This package was created to analyze multi-level one-way
         a confidence interval plot for Multiple Comparison.
 URL: https://github.com/yufree/rgabriel
 License: GPL (>= 2)
-Packaged: 2013-12-15 09:45:25 UTC; yufree
+Packaged: 2013-12-28 07:19:31 UTC; yufree
 NeedsCompilation: no
 Repository: CRAN
-Date/Publication: 2013-12-16 07:43:56
+Date/Publication: 2013-12-28 09:08:33
diff --git a/MD5 b/MD5
@@ -1,7 +1,7 @@
-62f2b9ee9ad5928a4861f144907d41b5 *DESCRIPTION
+b29fc3e280268882adc47e8da5d163f6 *DESCRIPTION
 8b54e5a89fbda3af5e077053d40bec76 *NAMESPACE
 1055bc274b8ce6f57794533c10d808e4 *R/gabriel.plot.R
-a312991190ad1378a477db4896e55548 *R/rgabriel.R
+a39efb6e533d92febef931713ccce5c2 *R/rgabriel.R
 9b08a826b0daaf47b2dbd30b211a1905 *man/gabriel.plot.Rd
-17325b0fdbbfa575677c12b5d5aa7ca7 *man/rgabriel-package.Rd
-e1ea66f9bf12e6e527f09ad405e98f87 *man/rgabriel.Rd
+b03894f41c677dcb0ecf38bb9b309a87 *man/rgabriel-package.Rd
+adde413f774922229e997fc4bf5f204c *man/rgabriel.Rd
diff --git a/R/rgabriel.R b/R/rgabriel.R
@@ -1,6 +1,6 @@
 ###rgabriel
 ###Function for performing the procedure described in Gabriel's 1978 paper
-###13 Dec 2013
+###28 Dec 2013
 
 rgabriel <- function (x, f, a = "alpha level") 
 {
@@ -26,218 +26,27 @@ rgabriel <- function (x, f, a = "alpha level")
     df[i]  <-  n[i] - 1
   }
   s <- tapply(x,f,sd)
-#  if(length(unique(n))==1) {
-#  SR  <-  qtukey(p = a, nmeans = k, df = df, lower.tail = FALSE)
-#  vstar  <-  SR*s/(2*sqrt(n))
-#  return(vstar)
-#  } else if(a == 0.05){
-  if(a == 0.05){
-    kstar <- choose(k,2)
-    dfstar <- sum(df)
-    smmcrit<-function(nuhat,C){
-      #
-      #  Determine the .95 quantile of the C-variate Studentized maximum
-      #  modulus distribution using linear interpolation on inverse
-      #  degrees of freedom
-      #  If C=1, this function returns the .975 quantile of Student's t
-      #  distribution.
-      #
-      if(C-round(C)!=0)stop("The number of contrasts, C, must be an  integer")
-      if(C>=29)stop("C must be less than or equal to 28")
-      if(C<=0)stop("C must be greater than or equal to 1")
-      if(nuhat<2)stop("The degrees of freedom must be greater than or equal to 2")
-      if(C==1)smmcrit<-qt(.975,nuhat)
-      if(C>=2){
-        C<-C-1
-        m1<-matrix(0,20,27)
-        m1[1,]<-c(5.57,6.34,6.89,7.31,7.65,7.93,8.17,8.83,8.57,
-                  8.74,8.89,9.03,9.16,9.28,9.39,9.49,9.59, 9.68,
-                  9.77,9.85,9.92,10.00,10.07,10.13,10.20,10.26,10.32)
-        m1[2,]<-c(3.96,4.43,4.76,5.02,5.23,5.41,5.56,5.69,5.81,
-                  5.92,6.01,6.10,6.18,6.26,6.33,6.39,6.45,6.51,
-                  6.57,6.62,6.67,6.71,6.76,6.80,6.84,6.88, 6.92)
-        m1[3,]<-c(3.38,3.74,4.01,4.20,4.37,4.50,4.62,4.72,4.82,
-                  4.89,4.97,5.04,5.11,5.17,5.22,5.27,5.32, 5.37,
-                  5.41,5.45,5.49,5.52,5.56,5.59,5.63,5.66,5.69)
-        m1[4,]<-c(3.09,3.39,3.62,3.79,3.93,4.04,4.14,4.23,4.31,
-                  4.38,4.45,4.51,4.56,4.61,4.66,4.70,4.74,4.78,
-                  4.82,4.85,4.89,4.92,4.95,4.98,5.00,5.03,5.06)
-        m1[5,]<-c(2.92,3.19,3.39,3.54,3.66,3.77,3.86,3.94,4.01,
-                  4.07,4.13,4.18,4.23,4.28,4.32,4.36,4.39,4.43,
-                  4.46,4.49,4.52,4.55,4.58,4.60,4.63,4.65,4.68)
-        m1[6,]<-c(2.80,3.06,3.24,3.38,3.49,3.59,3.67,3.74,3.80,
-                  3.86,3.92,3.96,4.01,4.05,4.09,4.13,4.16,4.19,
-                  4.22,4.25,4.28,4.31,4.33,4.35,4.38,4.39,4.42)
-        m1[7,]<-c(2.72,2.96,3.13,3.26,3.36,3.45,3.53,3.60,3.66,
-                  3.71,3.76,3.81,3.85,3.89,3.93,3.96,3.99, 4.02,
-                  4.05,4.08,4.10,4.13,4.15,4.18,4.19,4.22,4.24)
-        m1[8,]<-c(2.66,2.89,3.05,3.17,3.27,3.36,3.43,3.49,3.55,
-                  3.60,3.65,3.69,3.73,3.77,3.80,3.84,3.87,3.89,
-                  3.92,3.95,3.97,3.99,4.02,4.04,4.06,4.08,4.09)
-        m1[9,]<-c(2.61,2.83,2.98,3.10,3.19,3.28,3.35,3.41,3.47,
-                  3.52,3.56,3.60,3.64,3.68,3.71,3.74,3.77,3.79,
-                  3.82,3.85,3.87,3.89,3.91,3.94,3.95, 3.97,3.99)
-        m1[10,]<-c(2.57,2.78,2.93,3.05,3.14,3.22,3.29,3.35,3.40,
-                   3.45,3.49,3.53,3.57,3.60,3.63,3.66,3.69,3.72,
-                   3.74,3.77,3.79,3.81,3.83,3.85,3.87,3.89,3.91)
-        m1[11,]<-c(2.54,2.75,2.89,3.01,3.09,3.17,3.24,3.29,3.35,
-                   3.39,3.43,3.47,3.51,3.54,3.57,3.60,3.63,3.65,
-                   3.68,3.70,3.72,3.74,3.76,3.78,3.80,3.82,3.83)
-        m1[12,]<-c(2.49,2.69,2.83,2.94,3.02,3.09,3.16,3.21,3.26,
-                   3.30,3.34,3.38,3.41,3.45,3.48,3.50,3.53,3.55,
-                   3.58,3.59,3.62,3.64,3.66,3.68,3.69,3.71,3.73)
-        m1[13,]<-c(2.46,2.65,2.78,2.89,2.97,3.04,3.09,3.15,3.19,
-                   3.24,3.28,3.31,3.35,3.38,3.40,3.43,3.46,3.48,
-                   3.50,3.52,3.54,3.56,3.58,3.59,3.61,3.63,3.64)
-        m1[14,]<-c(2.43,2.62,2.75,2.85,2.93,2.99,3.05,3.11,3.15,
-                   3.19,3.23,3.26,3.29,3.32,3.35,3.38,3.40,3.42,
-                   3.44,3.46,3.48,3.50,3.52,3.54,3.55,3.57,3.58)
-        m1[15,]<-c(2.41,2.59,2.72,2.82,2.89,2.96,3.02,3.07,3.11,
-                   3.15,3.19,3.22,3.25,3.28,3.31,3.33,3.36,3.38,
-                   3.39,3.42,3.44,3.46,3.47,3.49,3.50,3.52,3.53)
-        m1[16,]<-c(2.38,2.56,2.68,2.77,2.85,2.91,2.97,3.02,3.06,
-                   3.09,3.13,3.16,3.19,3.22,3.25,3.27,3.29,3.31,
-                   3.33,3.35,3.37,3.39,3.40,3.42,3.43,3.45,3.46)
-        m1[17,]<-c(2.35,2.52,2.64,2.73,2.80,2.87,2.92,2.96,3.01,
-                   3.04,3.07,3.11,3.13,3.16,3.18,3.21,3.23,3.25,
-                   3.27,3.29,3.30,3.32,3.33,3.35,3.36,3.37,3.39)
-        m1[18,]<-c(2.32,2.49,2.60,2.69,2.76,2.82,2.87,2.91,2.95,
-                   2.99,3.02,3.05,3.08,3.09,3.12,3.14,3.17, 3.18,
-                   3.20,3.22,3.24,3.25,3.27,3.28,3.29,3.31,3.32)
-        m1[19,]<-c(2.29,2.45,2.56,2.65,2.72,2.77,2.82,2.86,2.90,
-                   2.93,2.96,2.99,3.02,3.04,3.06,3.08,3.10, 3.12,
-                   3.14,3.16,3.17,3.19,3.20,3.21,3.23,3.24,3.25)
-        m1[20,]<-c(2.24,2.39,2.49,2.57,2.63,2.68,2.73,2.77,2.79,
-                   2.83,2.86,2.88,2.91,2.93,2.95,2.97,2.98, 3.01,
-                   3.02,3.03,3.04,3.06,3.07,3.08,3.09,3.11,3.12)
-        if(nuhat>=200)smmcrit<-m1[20,C]
-        if(nuhat<200){
-          nu<-c(2,3,4,5,6,7,8,9,10,11,12,14,16,18,20,24,30,40,60,200)
-          temp<-abs(nu-nuhat)
-          find<-order(temp)
-          if(temp[find[1]]==0)smmcrit<-m1[find[1],C]
-          if(temp[find[1]]!=0){
-            if(nuhat>nu[find[1]]){
-              smmcrit<-m1[find[1],C]-
-                (1/nu[find[1]]-1/nuhat)*(m1[find[1],C]-m1[find[1]+1,C])/
-                (1/nu[find[1]]-1/nu[find[1]+1])
-            }
-            if(nuhat<nu[find[1]]){
-              smmcrit<-m1[find[1]-1,C]-
-                (1/nu[find[1]-1]-1/nuhat)*(m1[find[1]-1,C]-m1[find[1],C])/
-                (1/nu[find[1]-1]-1/nu[find[1]])
-            }
-          }}
-      }
-      smmcrit
-    }
-    SR <- smmcrit(dfstar,kstar)
-    vstar <- SR*s/sqrt(2*n)
-    return(vstar)
-  } else if(a == 0.01){
-    kstar <- choose(k,2)
-    dfstar <- sum(df)
-    smmcrit01<-function(nuhat,C){
-      #
-      #  Determine the .99 quantile of the C-variate Studentized maximum
-      #  modulus distribution using linear interpolation on inverse
-      #  degrees of freedom
-      #  If C=1, this function returns the .995 quantile of Student's t
-      #  distribution.
-      #
-      if(C-round(C)!=0)stop("The number of contrasts, C, must be an  integer")
-      if(C>=29)stop("C must be less than or equal to 28")
-      if(C<=0)stop("C must be greater than or equal to 1")
-      if(nuhat<2)stop("The degrees of freedom must be greater than or equal to 2")
-      if(C==1)smmcrit01<-qt(.995,nuhat)
-      if(C>=2){
-        C<-C-1
-        m1<-matrix(0,20,27)
-        m1[1,]<-c(12.73,14.44,15.65,16.59,17.35,17.99,18.53,19.01,19.43,
-                  19.81,20.15,20.46,20.75,20.99,20.99,20.99,20.99,20.99,
-                  22.11,22.29,22.46,22.63,22.78,22.93,23.08,23.21,23.35)
-        m1[2,]<-c(7.13,7.91,8.48,8.92,9.28,9.58,9.84,10.06,10.27,
-                  10.45,10.61,10.76,10.90,11.03,11.15,11.26,11.37,11.47,
-                  11.56,11.65,11.74,11.82,11.89,11.97,12.07,12.11,12.17)
-        m1[3,]<-c(5.46,5.99,6.36,6.66,6.89,7.09,7.27,7.43,7.57,
-                  7.69,7.80,7.91,8.01,8.09,8.17,8.25,8.32,8.39,
-                  8.45,8.51,8.57,8.63,8.68,8.73,8.78,8.83,8.87)
-        m1[4,]<-c(4.70,5.11,5.39,5.63,5.81,5.97,6.11,6.23,6.33,
-                  6.43,6.52,6.59,6.67,6.74,6.81,6.87,6.93,6.98,
-                  7.03,7.08,7.13,7.17,7.21,7.25,7.29,7.33,7.36)
-        m1[5,]<-c(4.27,4.61,4.85,5.05,5.20,5.33,5.45,5.55,5.64,
-                  5.72,5.79,5.86,5.93,5.99,6.04,6.09,6.14,6.18,
-                  6.23,6.27,6.31,6.34,6.38,6.41,6.45,6.48,6.51)
-        m1[6,]<-c(3.99,4.29,4.51,4.68,4.81,4.93,5.03,5.12,5.19,
-                  5.27,5.33,5.39,5.45,5.50,5.55,5.59,5.64,5.68,
-                  5.72,5.75,5.79,5.82,5.85,5.88,5.91,5.94,5.96)
-        m1[7,]<-c(3.81,4.08,4.27,4.42,4.55,4.65,4.74,4.82,4.89,
-                  4.96,5.02,5.07,5.12,5.17,5.21,5.25,5.29, 5.33,
-                  5.36,5.39,5.43,5.45,5.48,5.51,5.54,5.56,5.59)
-        m1[8,]<-c(3.67,3.92,4.10,4.24,4.35,4.45,4.53,4.61,4.67,
-                  4.73,4.79,4.84,4.88,4.92,4.96,5.01,5.04,5.07,
-                  5.10,5.13,5.16,5.19,5.21,5.24,5.26,5.29,5.31)
-        m1[9,]<-c(3.57,3.80,3.97,4.09,4.20,4.29,4.37,4.44,4.50,
-                  4.56,4.61,4.66,4.69,4.74,4.78,4.81,4.84,4.88,
-                  4.91,4.93,4.96,4.99,5.01,5.03,5.06,5.08,5.09)
-        m1[10,]<-c(3.48,3.71,3.87,3.99,4.09,4.17,4.25,4.31,4.37,
-                   4.42,4.47,4.51,4.55,4.59,4.63,4.66,4.69,4.72,
-                   4.75,4.78,4.80,4.83,4.85,4.87,4.89,4.91,4.93)
-        m1[11,]<-c(3.42,3.63,3.78,3.89,.99,4.08,4.15,4.21,4.26,
-                   4.31,4.36,4.40,4.44,4.48,4.51,4.54,4.57,4.59,
-                   4.62,4.65,4.67,4.69,4.72,4.74,4.76,4.78,4.79)
-        m1[12,]<-c(3.32,3.52,3.66,3.77,3.85,3.93,3.99,.05,4.10,
-                   4.15,4.19,4.23,4.26,4.29,4.33,4.36,4.39,4.41,
-                   4.44,4.46,4.48,4.50,4.52,4.54,4.56,4.58,4.59)
-        m1[13,]<-c(3.25,3.43,3.57,3.67,3.75,3.82,3.88,3.94,3.99,
-                   4.03,4.07,4.11,4.14,4.17,4.19,4.23,4.25,4.28,
-                   4.29,4.32,4.34,4.36,4.38,4.39,4.42,4.43,4.45)
-        m1[14,]<-c(3.19,3.37,3.49,3.59,3.68,3.74,3.80,3.85,3.89,
-                   3.94,3.98,4.01,4.04,4.07,4.10,4.13,4.15,4.18,
-                   4.19,4.22,4.24,4.26,4.28,4.29,4.31,4.33,4.34)
-        m1[15,]<-c(3.15,3.32,3.45,3.54,3.62,3.68,3.74,3.79,3.83,
-                   3.87,3.91,3.94,3.97,3.99,4.03,4.05,4.07,4.09,
-                   4.12,4.14,4.16,4.17,4.19,4.21,4.22,4.24,4.25)
-        m1[16,]<-c(3.09,3.25,3.37,3.46,3.53,3.59,3.64,3.69,3.73,
-                   3.77,3.80,3.83,3.86,3.89,3.91,3.94,3.96,3.98,
-                   4.00,4.02,4.04,4.05,4.07,4.09,4.10,4.12,4.13)
-        m1[17,]<-c(3.03,3.18,3.29,3.38,3.45,3.50,3.55,3.59,3.64,
-                   3.67,3.70,3.73,3.76,3.78,3.81,3.83,3.85,3.87,
-                   3.89,3.91,3.92,3.94,3.95,3.97,3.98,4.00,4.01)
-        m1[18,]<-c(2.97,3.12,3.22,3.30,3.37,3.42,3.47,3.51,3.55,
-                   3.58,3.61,3.64,3.66,3.68,3.71,3.73,3.75,3.76,
-                   3.78,3.80,3.81,3.83,3.84,3.85,3.87,3.88,3.89)
-        m1[19,]<-c(2.91,3.06,3.15,3.23,3.29,3.34,3.38,3.42,3.46,
-                   3.49,3.51,3.54,3.56,3.59,3.61,3.63,3.64,3.66,
-                   3.68,3.69,3.71,3.72,3.73,3.75,3.76,3.77,3.78)
-        m1[20,]<-c(2.81,2.93,3.02,3.09,3.14,3.19,3.23,3.26,3.29,
-                   3.32,3.34,3.36,3.38,3.40,.42,.44,3.45,3.47,
-                   3.48,3.49,3.50,3.52,3.53,3.54,3.55,3.56,3.57)
-        if(nuhat>=200)smmcrit01<-m1[20,C]
-        if(nuhat<200){
-          nu<-c(2,3,4,5,6,7,8,9,10,11,12,14,16,18,20,24,30,40,60,200)
-          temp<-abs(nu-nuhat)
-          find<-order(temp)
-          if(temp[find[1]]==0)smmcrit01<-m1[find[1],C]
-          if(temp[find[1]]!=0){
-            if(nuhat>nu[find[1]]){
-              smmcrit01<-m1[find[1],C]-
-                (1/nu[find[1]]-1/nuhat)*(m1[find[1],C]-m1[find[1]+1,C])/
-                (1/nu[find[1]]-1/nu[find[1]+1])
-            }
-            if(nuhat<nu[find[1]]){
-              smmcrit01<-m1[find[1]-1,C]-
-                (1/nu[find[1]-1]-1/nuhat)*(m1[find[1]-1,C]-m1[find[1],C])/
-                (1/nu[find[1]-1]-1/nu[find[1]])
-            }
-          }}
-      }
-      smmcrit01
-    }
-    SR <- smmcrit01(dfstar,kstar)
-    vstar <- SR*s/sqrt(2*n)
-    return(vstar)
-  } else{
-    print('Sorry, I can not figure it out')
+  dfstar <- sum(df)
+  # the following code is written by Yihui according to Stoline's 1979 paper(equation 2.2). 
+  psmm_x = function(x, c, r, nu) {
+    snu = sqrt(nu)
+    sx  = snu * x  # for the scaled Chi distribution
+    lgx = log(snu) - lgamma(nu/2) + (1 - nu/2) * log(2) + (nu - 1) * log(sx) + (-sx^2/2)
+    exp(r * log(2 * pnorm(c * x) - 1) + lgx)
+  }
+  psmm = function(x, r, nu) {
+    res = integrate(psmm_x, 0, Inf, c = x, r = r, nu = nu)
+    res$value
+  }
+  qsmm = function(q, r, nu) {
+    r = (r * (r - 1)/2)
+    if (!is.finite(nu)) return(qnorm(1 - .5 *(1 - q^(1/r))))
+    res = uniroot(function(c, r, nu, q) {
+      psmm(c, r = r, nu = nu) - q
+    }, c(0, 100), r = r, nu = nu, q = q)
+    res$root
   }
+  SR <- qsmm(1-a, k, dfstar)
+  vstar <- SR*s/sqrt(2*n)
+  return(vstar)
 }
diff --git a/man/rgabriel-package.Rd b/man/rgabriel-package.Rd
@@ -11,25 +11,28 @@ Functions for conducting and plotting Gabriel's (1978) multiple comparison test
 \tabular{ll}{
 Package: \tab rgabriel\cr
 Type: \tab Package\cr
-Version: \tab 0.4\cr
-Date: \tab 2013-12-15\cr
+Version: \tab 0.7\cr
+Date: \tab 2013-12-28\cr
 License: \tab GPL (>= 2)\cr
 }
 This package was created to analyze multi-level one-way experimental designs. It is designed to handle vectorized observation and factor data where there are unequal sample sizes and population variance homogeneity can not be assumed.To conduct the Gabriel test, create two vectors: one for your observations and one for the factor level of each observation. The function, rgabriel, conduct the test and save the output as a vector to input into the gabriel.plot function, which produces a confidence interval plot for Multiple Comparison.  
 }
 \author{
-Miao YU, PhD student in Research Center for Eco-Environmental Sciences(RCEES), Chinese Academy of Sciences
 
-Maintainer: Miao YU <yufreecas@gmail.com>
+Yihui XIE <\url{http://yihui.name}>
+
+Miao YU, PhD student in Research Center for Eco-Environmental Sciences(RCEES), Chinese Academy of Sciences
 
 }
 \note{
-I would like to acknowledge the invaluable help of the author of http://www.gmw.rug.nl/~boomsma/apstatdata/wilcox20.R who show a table of studentized maximum modulus's distribution.
+More details on the simulation of studentized maximum modulus's distribution from http://cos.name/cn/topic/142002. 
+
 }
 
 \references{
 Gabriel, K.R., 1978. A Simple Method of Multiple Comparisons of Means. Journal of the American Statistical Association 73, 724. 
-http://www.gmw.rug.nl/~boomsma/apstatdata/wilcox20.R
+
+Stoline, M.R., Ury, H.K., 1979. Tables of the Studentized Maximum Modulus Distribution and an Application to Multiple Comparisons among Means. Technometrics 21, 87.
 }
 \keyword{ package, plot, Gabriel }
 

diff --git a/man/rgabriel.Rd b/man/rgabriel.Rd
@@ -17,7 +17,7 @@ data vector
 factor vector
 }
   \item{a}{
-alpha level of mutiple comparison. Only alpha = 0.01 or 0.05 could be used as alpha level.
+alpha level of mutiple comparison. 
 }
 }
 \details{
@@ -29,9 +29,11 @@ As shown in Gabriel's paper,use M(alpha,k*,v), the upper alpha point of the Stud
 \references{
 Gabriel, K.R., 1978. A Simple Method of Multiple Comparisons of Means. Journal of the American Statistical Association 73, 724. 
 
-http://www.gmw.rug.nl/~boomsma/apstatdata/wilcox20.R
+Stoline, M.R., Ury, H.K., 1979. Tables of the Studentized Maximum Modulus Distribution and an Application to Multiple Comparisons among Means. Technometrics 21, 87.
 }
 \author{
+Yihui XIE
+
 Miao YU
 }