Skip to content

Commit

Permalink
version 0.4
Browse files Browse the repository at this point in the history
  • Loading branch information
@gaborcsardi
Miao YU authored and gaborcsardi committed Dec 15, 2013
0 parents commit f13f92d
Show file tree
Hide file tree
Showing 8 changed files with 457 additions and 0 deletions.
23 changes: 23 additions & 0 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
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
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
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.
URL: https://github.com/yufree/rgabriel
License: GPL (>= 2)
Packaged: 2013-12-15 09:45:25 UTC; yufree
NeedsCompilation: no
Repository: CRAN
Date/Publication: 2013-12-16 07:43:56
7 changes: 7 additions & 0 deletions MD5
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
62f2b9ee9ad5928a4861f144907d41b5 *DESCRIPTION
8b54e5a89fbda3af5e077053d40bec76 *NAMESPACE
1055bc274b8ce6f57794533c10d808e4 *R/gabriel.plot.R
a312991190ad1378a477db4896e55548 *R/rgabriel.R
9b08a826b0daaf47b2dbd30b211a1905 *man/gabriel.plot.Rd
17325b0fdbbfa575677c12b5d5aa7ca7 *man/rgabriel-package.Rd
e1ea66f9bf12e6e527f09ad405e98f87 *man/rgabriel.Rd
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
exportPattern("^[[:alpha:]]+")
11 changes: 11 additions & 0 deletions R/gabriel.plot.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
gabriel.plot <- function(x, f, upper, lower=upper, length=0.1,...){
if (class(f) != "factor") {
f <- factor(f)
}
input <- cbind(f[!sapply(is.na(x), all)], x[!sapply(is.na(x), all)])
f <- factor(input[, 1])
x <- input[, 2]
meang <- tapply(x,f,mean)
tem <- barplot(t(meang),ylim=c(0,max(meang)+2*max(upper)),...)
arrows(tem,meang+upper, tem, meang-lower, angle=90, code=3, length=length, ...)
}
243 changes: 243 additions & 0 deletions R/rgabriel.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,243 @@
###rgabriel
###Function for performing the procedure described in Gabriel's 1978 paper
###13 Dec 2013

rgabriel <- function (x, f, a = "alpha level")
{
if (a == "alpha level") {
a <- 0.05
}
if (class(f) != "factor") {
f <- factor(f)
}
Level.Name <- levels(f)
k <- length(Level.Name)
input <- cbind(f[!sapply(is.na(x), all)], x[!sapply(is.na(x), all)])
f <- factor(input[, 1])
x <- input[, 2]
n <- numeric(length <- nlevels(f))
for (i in 1:length(n)) {
n[i] <- length(f[f == i])
}
for (i in 1:length(n)) {
if (i == 1) {
df <- numeric(length = length(n))
}
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')
}
}
60 changes: 60 additions & 0 deletions man/gabriel.plot.Rd
Original file line number Diff line number Diff line change
@@ -0,0 +1,60 @@
\name{gabriel.plot}
\alias{gabriel.plot}
\title{
the Gabriel's barplot (or (l-u)-plot)
}
\description{
Make the Gabriel's barplot, if, and only if, their bar intervals are disjoint, they are differ significantly. This function could also be used to plot error bar when the bar vector is imported as upper or lower margin.
}
\usage{
gabriel.plot(x, f, upper, lower = upper, length = 0.1, ...)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{x}{
data vector
}
\item{f}{
factor vector
}
\item{upper}{
the upper margin of error bar
}
\item{lower}{
the upper margin of error bar
}
\item{length}{
the length of error bar
}
\item{...}{
Arguments to be passed to methods, such as graphical parameters.
}
}

\references{
Gabriel, K.R., 1978. A Simple Method of Multiple Comparisons of Means. Journal of the American Statistical Association 73, 724.
}
\author{
Miao YU
}

\seealso{
\code{\link{rgabriel}}, \code{\link{barplot}}
}
\examples{
# equal numbers

g <- c(1:50)
f <- c(rep(1,10),rep(2,10),rep(3,10),rep(4,10),rep(5,10))
gabriel.plot(g,f,rgabriel(g,f))

# unequal numbers

g <- c(1:40)
f <- c(rep(1,3),rep(2,12),rep(3,15),rep(4,5),rep(5,5))
gabriel.plot(g,f,rgabriel(g,f))
}

\keyword{ Gabriel }
\keyword{ plot }
\keyword{ error bar }

0 comments on commit f13f92d

Please sign in to comment.