##########
### Handlings of multi-way tables
##########
# Marginal count
# A is an array
#' Calculate marginal counts of array
#'
#' @param A An array or matrix
#' @return A list of marginal counts per variable
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' calc.marg(A)
#' @export
calc.marg <- function(A){
d <- dim(A)
ret <- list()
for(i in 1:length(d)){
ret[[i]] <- apply(A,i,sum)
}
ret
}
# Expected table
#' Make an expected table
#'
#' @param A an array or matrix
#' @return The expected table in a shape of array or matrix
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' make.exp.table(A)
#' @export
make.exp.table <- function(A){
n <- sum(A)
marg <- calc.marg(A)
tmp <- marg[[1]]
for(i in 2:length(marg)){
tmp <- t(marg[[i]]/n) %x% tmp
}
tmp
}
# Difference table
#' Make a table, the given table minus its expected table
#'
#' @param A an array or matrix
#' @return The difference table in a shape of array or matrix
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' make.diff.table(A)
#' @export
make.diff.table <- function(A){
E <- make.exp.table(A)
array(c(A)-c(E),dim(A))
}
# Chisquare value
#' Calculate chi-square value of an array or matrix
#'
#' @param A an array or matrix
#' @return The chi-square value of the table
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' calc.chisq(A)
#' @export
calc.chisq <- function(A){
E <- make.exp.table(A)
D <- make.diff.table(A)
sum(c(D)^2/c(E))
}
# Simplex and rotation of simplex
# (d-1)-simplex is an polytope with d vertices in d-dimensional space
# This function gives coordinates of d vertices in d-1 dimensional space
#' Coordinate of simplex vertices
#'
#' @param d Number of vertices of a simplex
#' @return d times (d-1) matrix, each row of which is a coordinate of a vertex of (d-1)-simplex, whose center is the origin and distance of whose vertices from the origin is 1
#' @examples
#' d <- 3
#' cv <- CategoryVector(d)
#' print(cv)
#' apply(cv^2,1,sum)
#' cv %*% t(cv)
#' @export
CategoryVector<-function(d){
df <- d - 1
diagval <- 1:d
diagval <- sqrt((d)/df) * sqrt((d - diagval)/(d - diagval + 1))
others <- -diagval/(d - (1:d))
m <- matrix(rep(others, d), nrow = d, byrow = TRUE)
diag(m) <- diagval
m[upper.tri(m)] <- 0
as.matrix(m[, 1:df])
}
# This function gives coordinates of d vertices in "d" dimensional space with d-th coordinate of all vertices being equal and distance of all vertices from origin being 1.
# This is a rotation matrix.
# We call the output "simplex-rotation matrix".
#' Make a rotation matrix using simple vertices coordinates (1)
#'
#' @param k A integer
#' @return k x k rotation matrix to transfer k vertices so that all k vertices' k-th coordinates are the same
#' @examples
#' k <- 3
#' M <- make.simplex.0(k)
#' M %*% diag(rep(1,k))
#' @export
make.simplex.0 <- function(k){
cv <- CategoryVector(k)
rbind(t(cv*sqrt(1-1/k)),rep(1/sqrt(k),k))
}
#' Make a rotation matrix using simple vertices coordinates (2)
#'
#' @param k A integer
#' @return k x k rotation matrix to transfer k vertices so that all k vertices' k-th coordinates are the same
#' @examples
#' k <- 3
#' M <- make.simplex(k)
#' M %*% diag(rep(1,k))
#' @export
make.simplex <- function(k){
ret <- matrix(0,k,k)
for(i in 1:(k-1)){
for(j in 1:k){
if(j < i){
}else if(j==i){
ret[i,j] <- sqrt((k-i)/(k-i+1))
}else{
ret[i,j] <- -sqrt(1/((k-i)*(k-i+1)))
}
}
}
ret[k,] <- sqrt(1/k)
ret
}
# Kronecker product of multiple simplex-rotation matrices
#' Kronecker product of simplex-based rotation matrices
#'
#' @param r A vector of positive integers, a vector of number of levels of variables
#' @return Kronecker product of rotation matrices, each of which is a simplex-based rotation matrix of number of levels of each variable
#' @examples
#' r <- c(2,3,4)
#' KM <- make.simplex.multi(r)
#' dim(KM)
#' @export
make.simplex.multi <- function(r){
X <- make.simplex(r[1])
k <- length(r)
if(k > 1){
for(i in 2:k){
X <- make.simplex(r[i]) %x% X
}
}
X
}
# This function's output indicates which elements of table vectors should be zero after the kronecker product rotation.
# r is dimension vector or multi-way table
#' Make a vector indexing rows that should be zero
#'
#' @param r A vector of positive integers, a vector of number of levels of variables
#' @return A vector with 0 or 1; indices with 1 means that elements of table vectors should be zero after the kronecker product rotation.
#' @examples
#' r <- c(2,3,4)
#' v0 <- make.vector0(r)
#' @export
make.vector0 <- function(r){
tmp <- list()
for(i in 1:length(r)){
tmp[[i]] <- c(rep(1,r[i]-1),0)
}
tmp2 <- apply(expand.grid(tmp),1,sum)
as.numeric(tmp2>1)
}
# r is dimension vector or multi-way table
# Kronecker-product rotation matrix: X
# The meaningfull rows of X: X.sub
# The X^(-1) = X^T: X.inv
# The X.sub^(-1): X.inv.sub
# The vector indicating which rows are meaningfull/less
#' Multiway table-related rotation matrix and its features
#'
#' @param r A vector of positive integers, a vector of number of levels of variables
#' @return X A Rotation matrix; output of make.simplex.multi()
#' @return X.sub A matrix only with meaningful rows of X
#' @return X.inv Inverse matrix of X
#' @return X.inv.sub Meaningful part of X.inv
#' @return vector0 A vector indicating meaningful rows of X
#' @examples
#' r <- c(2,3,4)
#' mXv <- make.X.vector0(r)
#' mXv$X %*% mXv$X.inv
#' @export
make.X.vector0 <- function(r){
X <- make.simplex.multi(r)
#X.inv <- solve(X)
X.inv <- t(X)
vector0 <- make.vector0(r)
non.zero <- which(vector0==1)
return(list(X=X,X.sub=X[non.zero,],X.inv=X.inv,X.inv.sub=X.inv[,non.zero],vector0=vector0))
}
###
# multi-way table is vectorized into a R-length vector and it is rotated with map2df matrix to a vector whose R-df elements are zero.
#' Make a matrix that transfer a vectorized table array to a vector in the df-dimensional space
#'
#' @param A An array or matrix of multiway table
#' @return A matrix
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' M <- map2df(A)
#' M %*% c(A)
#' @export
map2df <- function(A){
r <- dim(A)
out <- make.X.vector0(r)
E.A <- make.exp.table(A)
W <- t(out$X.inv.sub) %*% diag(1/c(E.A)) %*% (out$X.inv.sub)
eigen.W <- eigen(W)
tmpdiag <- matrix(0,length(eigen.W[[1]]),length(eigen.W[[1]]))
diag(tmpdiag) <- sqrt(eigen.W[[1]])
tmpdiag %*% t(eigen.W[[2]]) %*% out$X.sub
}
#' Make a matrix that transfer back a vector in the df-dimensional space to a vector with number-of-cell of table elements
#'
#' @param A An array or matrix of multiway table
#' @return A matrix
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' M <- map2df(A)
#' M.inv <- map2full(A)
#' M.inv %*% (M %*% c(A))
#' @export
map2full <- function(A){
r <- dim(A)
out <- make.X.vector0(r)
E.A <- make.exp.table(A)
W <- t(out$X.inv.sub) %*% diag(1/c(E.A)) %*% (out$X.inv.sub)
eigen.W <- eigen(W)
tmpdiag <- matrix(0,length(eigen.W[[1]]),length(eigen.W[[1]]))
diag(tmpdiag) <- 1/sqrt(eigen.W[[1]])
out$X.inv.sub %*% eigen.W[[2]] %*% tmpdiag
}
# A test vector in R-dimensional space is roted to df-dimensional vector with this matrix
#' Make a matrix that transfer test vectors to the direction in the df-dimensional space
#'
#' @param A An array or matrix of multiway table
#' @return A matrix
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' M <- map2dfdir(A)
#' c(A) %*% M
#' @export
map2dfdir <- function(A){
r <- dim(A)
out <- make.X.vector0(r)
E.A <- make.exp.table(A)
W <- t(out$X.inv.sub) %*% diag(1/c(E.A)) %*% (out$X.inv.sub)
eigen.W <- eigen(W)
tmpdiag <- matrix(0,length(eigen.W[[1]]),length(eigen.W[[1]]))
diag(tmpdiag) <- 1/sqrt(eigen.W[[1]])
out$X.inv.sub %*% (eigen.W[[2]]) %*% tmpdiag
}
###
# 3 matrices mentioned above are bundled.
#' Make three matrices that are the output of map2df(), map2full() and map2dfdir
#'
#' @param A An array or matrix of multiway table
#' @return Full2DF A matrix transferring table vectors to df-dimension vectors
#' @return DF2full A matrix transferring back vectors in df-dimensional space back to vectors with number-of-cell elements
#' @return Normalvec2DF A matric transferring test vectors to df-dimension vectors
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' map.matrix(A)
#' @export
map.matrix <- function(A){
r <- dim(A)
out <- make.X.vector0(r)
E.A <- make.exp.table(A)
W <- t(out$X.inv.sub) %*% diag(1/c(E.A)) %*% (out$X.inv.sub)
eigen.W <- eigen(W)
tmpdiag <- matrix(0,length(eigen.W[[1]]),length(eigen.W[[1]]))
diag(tmpdiag) <- sqrt(eigen.W[[1]])
Full2DF <- tmpdiag %*% t(eigen.W[[2]]) %*% out$X.sub
tmpdiag <- matrix(0,length(eigen.W[[1]]),length(eigen.W[[1]]))
diag(tmpdiag) <- 1/sqrt(eigen.W[[1]])
DF2full <- out$X.inv.sub %*% eigen.W[[2]] %*% tmpdiag
NormalVec2DF <- out$X.inv.sub %*% (eigen.W[[2]]) %*% tmpdiag
return(list(Full2DF = Full2DF,DF2full = DF2full, NormalVec2DF = NormalVec2DF))
}
# multi-dimensional random vectors in normal distribution whose length is 1.
#' Generate random unit vectors
#'
#' @param n Number of random vecotrs to be generated
#' @param df An integer indicating dimension
#' @return A n x df matrix, each row of which is a random unit vector
#' @examples
#' n <- 100
#' df <- 2
#' plot(st.mvn(n,df))
#' @export
st.mvn <- function(n,df){
R <- matrix(rnorm(n*df),ncol=df)
L <- sqrt(apply(R^2,1,sum))
R/L
}
# integrate segments with "OR"
#' Union of segment(s) of multiple one-dimensional segments
#'
#' @param s A matrix with two columns each row of which indicating start and end of each segment
#' @return A matrix with two columns each row of which indicating unified segment
#' @examples
#' s <- matrix(c(0,4,1,6,2,3),byrow=TRUE,ncol=2)
#' segment.union(s)
#' @export
segment.union <- function(s){
n.seg <- length(s[,1])
my.seg.1 <- t(apply(s,1,sort))
ord <- order(my.seg.1[,1])
my.seg.2 <- my.seg.1[ord,]
st.end <- rep(0:1,n.seg)
st.end
loc <- c(t(my.seg.2))
loc.ord <- order(loc)
loc.ordered <- loc[loc.ord]
st.end.ordered <- st.end[loc.ord]
new.st <- c()
current <- 0
cnt <- 0
ret <- matrix(0,0,2)
for(i in 1:length(loc)){
if(st.end.ordered[i]==0){
if(length(new.st)==0){
new.st <- loc.ordered[i]
}
current <- current + 1
}else{
if(current==1){
cnt <- cnt+1
#ret[[cnt]] <- c(new.st,loc.ordered[i])
ret <- rbind(ret,c(new.st,loc.ordered[i]))
new.st <- c()
}
current <- current - 1
}
}
ret
}
# integrate segments with "AND"
#' Overlapped segment(s) of multiple one-dimensional segments
#'
#' @param s A matrix with two columns each row of which indicating start and end of each segment
#' @return A matrix with two columns each row of which indicating overlapped segment
#' @examples
#' s <- matrix(c(0,4,1,6,2,3),byrow=TRUE,ncol=2)
#' segment.overlap(s)
#' @export
segment.overlap <- function(s){
n.seg <- length(s[,1])
my.seg <- t(apply(s,1,sort))
max1 <- max(my.seg[,1])
min2 <- min(my.seg[,2])
ret <- c(max1,min2)
if(max1>min2){
ret <- c(0,0)
}
ret
}
# Chisquare distribution for +chi^2 and (-chi)^2 and cumulative function from -inf -> 0 -> +inf
#' Cumulative distribtution of bidirectional chisquare distribution supporting -Inf to +Inf
#'
#' @param a Chi-square value or its negative
#' @param df Degrees of freedom
#' @return Cumulative probability
#' @examples
#' as <- seq(from=-5,to=5,length=101)
#' plot(pchisq.bid(as,1))
#' @export
pchisq.bid <- function(a,df){
tmp <- as.numeric(a>0) - sign(a)*pchisq(a^2,df,lower.tail = FALSE)/2
tmp[which(a==0)] <- 0.5
tmp
}
###
# Marginal counts-dependent matarials
# Arguments
# A multi-way table (array)
# Values
# tables=list(data=A,expected=exp.table,diff=diff.table)
## A: input array, expected: expected array, diff: diff array
# r.vec = r.vec
## Dimensions of array
# df=df
## degree of freedom
# matrices = list(X=out$X,X.sub=out$X.sub,X.inv=out$X.inv,X.inv.sub=out$X.inv.sub)
## matrices to rotate simplex between full and df dimensional spaces
# zero.vec=out$vector0
## Index vector indicating zero elements after rotation to df space
# map.matrices=list(Full2DF=Full2DF,DF2full=DF2full,NormalVec2DF=NormalVec2DF,W=W)))
##
###
#' Spherization of multiway table array
#'
#' @param A An array or matrix of multiway table
#' @return tables A list of three tables: data :Input table A itself, expected : its expected table and diff: table of differecen between A and expected table
#' @return r.vec A integer vector indicating numbers of levels of A
#' @return df Degrees of freedom of A
#' @return matrices A list of four matrices, X,X.sub,X.inv and X.inv.sub, which are output of make.X.vector0, rotating to/from number-of-cell dimension and df dimension with/without meaningless rows/columns
#' @return zero.vec A vector indicating meaningful rows of rotations
#' @return map.matrices A list of four matrices: Full2df is a matrix transforming a table array vector to df-dimensional point, DF2full is a matrix transforming df-dimensional point to a table vector, Normalvec2DF is a matrix transforming test table to df-dimensional vector and W is a key matrix to calculate other matrices
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' t.sphere <- table.sphere(A)
#' @export
table.sphere <- function(A){
n <- sum(A) # sample size
mg.cnt <- calc.marg(A) # marginal counts
exp.table <- make.exp.table(A) # expected array
diff.table <- make.diff.table(A) # original array - expected array
chisq <- calc.chisq(A) # chisq value
r.vec <- dim(A) # dimension vector
out <- make.X.vector0(r.vec) # zero element-indicating vector
df <- sum(out$vector0) # degree of freedom
W <- t(out$X.inv.sub) %*% diag(1/c(exp.table)) %*% (out$X.inv.sub)
eigen.W <- eigen(W)
tmp.diag <- tmp.diag.inv <- matrix(0,df,df)
diag(tmp.diag) <- sqrt(eigen.W[[1]])
diag(tmp.diag.inv) <- 1/sqrt(eigen.W[[1]])
# Rotate an observe table vector to spherical space and its length^2 = chisq
# |Full2DF %*% A|^2 = chisq
Full2DF <- tmp.diag %*% t(eigen.W[[2]]) %*% out$X.sub
# matrix rotates df-vector to full dimensional vector
DF2full <- out$X.inv.sub %*% eigen.W[[2]] %*% tmp.diag.inv
# Rotate test normal vector to spherical space normal vector
NormalVec2DF <- out$X.inv.sub %*% (eigen.W[[2]]) %*% tmp.diag.inv
return(list(tables=list(data=A,expected=exp.table,diff=diff.table),r.vec = r.vec,df=df,matrices = list(X=out$X,X.sub=out$X.sub,X.inv=out$X.inv,X.inv.sub=out$X.inv.sub),zero.vec=out$vector0,map.matrices=list(Full2DF=Full2DF,DF2full=DF2full,NormalVec2DF=NormalVec2DF,W=W)))
}
###
# Marginal counts-dependent matarials
# Difference from table.shepre() is that this function specifies expected table as B.
# Arguments
# A multi-way table (array), observed
# B multi-way table (array), expected
# Values
# tables=list(data=A,expected=exp.table,diff=diff.table)
## A: input array, expected: expected array, diff: diff array
# r.vec = r.vec
## Dimensions of array
# df=df
## degree of freedom
# matrices = list(X=out$X,X.sub=out$X.sub,X.inv=out$X.inv,X.inv.sub=out$X.inv.sub)
## matrices to rotate simplex between full and df dimensional spaces
# zero.vec=out$vector0
## Index vector indicating zero elements after rotation to df space
# map.matrices=list(Full2DF=Full2DF,DF2full=DF2full,NormalVec2DF=NormalVec2DF,W=W)))
##
###
#' Spherization of multiway table array
#'
#' @param A An array or matrix of multiway table
#' @return tables A list of three tables: data :Input table A itself, expected : its expected table and diff: table of differecen between A and expected table
#' @return r.vec A integer vector indicating numbers of levels of A
#' @return df Degrees of freedom of A
#' @return matrices A list of four matrices, X,X.sub,X.inv and X.inv.sub, which are output of make.X.vector0, rotating to/from number-of-cell dimension and df dimension with/without meaningless rows/columns
#' @return zero.vec A vector indicating meaningful rows of rotations
#' @return map.matrices A list of four matrices: Full2df is a matrix transforming a table array vector to df-dimensional point, DF2full is a matrix transforming df-dimensional point to a table vector, Normalvec2DF is a matrix transforming test table to df-dimensional vector and W is a key matrix to calculate other matrices
#' @examples
#' r <- c(2,3,4)
#' A <- array(1:prod(r),r)
#' B <- A
#' B <- array(rep(sum(A)/length(A),length(A)),r)
#' t.sphere <- table.sphere.2(A,B)
#' @export
table.sphere.2 <- function (A,B)
{
n <- sum(A)
mg.cnt <- calc.marg(A)
exp.table <- B
diff.table <- A-B
chisq <- sum((A-B)^2/B)
r.vec <- dim(A)
out <- make.X.vector0(r.vec)
#df <- sum(out$vector0)
rotA <- out$X %*% c(A)
rotB <- out$X %*% c(B)
df.list <- 1:(prod(r.vec)-1)
W <- t(out$X.inv[,df.list]) %*% diag(1/c(exp.table)) %*% (out$X.inv[,df.list])
eigen.W <- eigen(W)
tmp.diag <- tmp.diag.inv <- matrix(0, length(c(A))-1, length(c(A))-1)
diag(tmp.diag) <- sqrt(eigen.W[[1]])
diag(tmp.diag.inv) <- 1/sqrt(eigen.W[[1]])
Full2DF <- tmp.diag %*% t(eigen.W[[2]]) %*% out$X[df.list,]
DF2full <- out$X.inv[,df.list] %*% eigen.W[[2]] %*% tmp.diag.inv
NormalVec2DF <- out$X.inv[,df.list] %*% (eigen.W[[2]]) %*% tmp.diag.inv
return(list(tables = list(data = A, expected = exp.table,
diff = diff.table), r.vec = r.vec, df = df, matrices = list(X = out$X,
X.sub = out$X.sub, X.inv = out$X.inv, X.inv.sub = out$X.inv.sub),
zero.vec = out$vector0, map.matrices = list(Full2DF = Full2DF,
DF2full = DF2full, NormalVec2DF = NormalVec2DF, W = W)))
}
###
# Arguments
## t.sphere: output of table.sphere()
## test.table: indicates alternative hypothesis direction in the shape of array
## odds.ratio.table: is an array that indicates the cells with which odds ratio should be set at the value or. (1,2,3,4) indicate cells as (cells[1] * cells[4])/(cells[2] * cells[3]) --> or
## n and K is arguments when serach the final result
# Values
## alt.table: Vectorized array
## alt.vec: Vector in df space
###
#' Generate a table indicating alternative hypothesis
#'
#' @param t.sphere Output of table.sphere of a table of interest itself or a table sharing marginal counts with the table of interest
#' @param test.table A matrix or array of a table which indicates the direction of alternative hypothesis along with an argument odds.ratio.table
#' @param odds.ratio.table A matrix or array of a table with cell values being 0,1,2,3 or 4: Odds ratio of table is indicated by the (sum(cells of 1)*sum(cells of 4))/(sum(cells of 2)*sum(cells of 3))
#' @param or A number indicating odds ratio defined by odds.ratio.table argument
#' @param n Optional. An integer which makes a one-dimensional grid with n points, among which the table with odds ratio closest to the argument or is selected
#' @param K Optional. A number indicating the maximum value of square of chi-square value to search the table of alternative hypothsis
#' @return alt.table The estimated table of alternative hypothesis
#' @return alt.vec The vector of the estimated table in df-dimensional space
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t.sphere <- table.sphere(A)
#' test.table <- matrix(c(1,0.5,0,0,0,0),byrow=TRUE,2,3)
#' odds.ratio.table <- matrix(c(1,0,2,3,0,4),byrow=TRUE,2,3)
#' or <- 2
#' n <- 1000
#' K <- 10
#' make.alternative.table(t.sphere,test.table,odds.ratio.table,or,n,k)
#' @export
make.alternative.table <- function(t.sphere,test.table,odds.ratio.table,or,n,K){
# Default values
if(missing(n))n<-1000
if(missing(K))K<-10
# alt.hy vector in df-space
tmp <- c(test.table) %*% t.sphere$map.matrices$NormalVec2DF
# dif vector in full-space
dif <- t.sphere$map.matrices$DF2full %*% t(tmp)
e.table <- t.sphere$tables$expected
# Target OR should be calculated with XYZxyz
X <- sum(e.table[which(odds.ratio.table==1)])
Y <- sum(e.table[which(odds.ratio.table==2)])
Z <- sum(e.table[which(odds.ratio.table==3)])
W <- sum(e.table[which(odds.ratio.table==4)])
x <- sum(dif[which(odds.ratio.table==1)])
y <- sum(dif[which(odds.ratio.table==2)])
z <- sum(dif[which(odds.ratio.table==3)])
w <- sum(dif[which(odds.ratio.table==4)])
# or of test.table
tmp.or <- (X+x)*(W+w)/((Y+y)*(Z+z))
# tmp.or and target or are compared and based on the comparison
# search-range is defined
if(tmp.or > or){
ks <- seq(from=0,to=1,length=n)
}else{
ks <- seq(from=1,to=or/tmp.or*K,length=n)
}
# search-range is evenly spaced and their or is calculated
# then closest one is selected.
ors <- rep(0,n)
for(i in 1:n){
difs <- t.sphere$map.matrices$DF2full %*% t(tmp)*ks[i]
x <- sum(difs[which(odds.ratio.table==1)])
y <- sum(difs[which(odds.ratio.table==2)])
z <- sum(difs[which(odds.ratio.table==3)])
w <- sum(difs[which(odds.ratio.table==4)])
ors[i] <- (X+x)*(W+w)/((Y+y)*(Z+z))
}
ret <- mean(ks[which(abs(ors-or)==min(abs(ors-or)))])
ret2 <- t.sphere$map.matrices$DF2full %*% t(tmp)*ret
return(list(alt.table=array(c(e.table)+ret2,t.sphere$r.vec),alt.vec=tmp*ret))
}
# tests is a list of tables representing test model and df-dimentional test vectors are generated using t.sphere information
#' Generate test vectors from tables
#'
#' @param t.sphere table.sphere function's output
#' @param tests A list of arrays/matrix of tables representing tests
#' @return A number-of-test rows x df columns matrix of test vectos, each row of which is a test vector
#' @return alt.vec The vector of the estimated table in df-dimensional space
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t.sphere <- table.sphere(A)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' make.test.vecs(t.sphere,tests)
#' @export
make.test.vecs <- function(t.sphere,tests){
n.test <- length(tests)
test.vecs <- matrix(0,n.test,t.sphere$df)
for(i in 1:n.test){
test.vecs[i,] <- c(c(tests[[i]]) %*% t.sphere$map.matrices$NormalVec2DF)
}
L.test.vecs <- sqrt(apply(test.vecs^2,1,sum))
test.vecs/L.test.vecs
}
# df-dimensioanl unit random vectrors are generated and their inner product with test vectors (standardized statistic valus) area returned.
#' Generate random K statistic values for test vectors under null hypothesis
#'
#' @param test.vecs A number-of-test rows x df columns matrix of test vectos, each row of which is a test vector
#' @param rmvn Optional A matrix of random unit vectors
#' @param n number of random unit vectors
#' @return A number-of-unit-vector rows x number-of-test columns matrix of K statistic values
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t.sphere <- table.sphere(A)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' test.vecs <- make.test.vecs(t.sphere,tests)
#' runitstat(test.vecs)
#' @export
runitstat <- function(test.vecs,rmvn,n){
# Default values
if(missing(n))n <- 1000
df <- length(test.vecs[1,])
if(missing(rmvn)){
rmvn <- st.mvn(n,df)
}
matrix((test.vecs %*% t(rmvn)),nrow=length(rmvn[,1]))
}
# When alternative hx is true, additional term on urunitstat is necessary to calculate statistic values. This function returns the coefficent for the modification.
#' Calculate K statistics of tests for output of make.alternative.table
#'
#' @param alt.tables Output of make.alternative.table, that is a list of table and vector in df-dimensional space of alternative hypothesis
#' @param test.vecs A number-of-test rows x df columns matrix of test vectos, each row of which is a test vector
#' @return A vector of K statisitc values
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t.sphere <- table.sphere(A)
#' test.table <- matrix(c(1,0.5,0,0,0,0),byrow=TRUE,2,3)
#' odds.ratio.table <- matrix(c(1,0,2,3,0,4),byrow=TRUE,2,3)
#' or <- 2
#' n <- 1000
#' K <- 10
#' alt.out <- make.alternative.table(t.sphere,test.table,odds.ratio.table,or,n,k)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' test.vecs <- make.test.vecs(t.sphere,tests)
#' alt.intersect(alt.out,test.vecs)
#' @export
alt.intersect <- function(alt.tables,test.vecs){
c(test.vecs %*% matrix(alt.tables$alt.vec,ncol=1))
}
# Of note, ks argument should be square-root of chisq values to be evaluated.
# Two terms are returned.
# term for null hx, depending on random vectors and test vectors. k/a,-k/a
# term for alternative hx, depending on test vectors. -b/a
#' Generate terms to calculate K statistic values of tests for random unit vectors
#'
#' @param as NUmber-of-random-unit-vector rows x number-of-test matrix
#' @param ks A vector or K values (square of chi-square equivalent values) for which P-values should be estimated
#' @param bs Optional. A Vector of K values of alternative hypothesis table for tests
#' @return A list. k.per.a is 3-way array and is one set of terms to calculate K statistics, which is variable respect to unit vectors. b.per.a is another set of terms which K values for the alternative hypothesis itself
#' @examples
#' # utility function being used in other functions
#' @export
window.ab <- function(as,ks,bs){
zero.a <- which(as==0)
as[zero.a] <- max(abs(as))*2
# term for null hx, depending on random vectors and test vectors and cutoff test statistics. k/a,-k/a
k.per.a <- outer(1/as,ks,"*")
# term for alternative hx, depending on test vectors. -b/a
if(missing(bs)){
b.per.a <- NULL
}else{
#b.per.a <- -outer(1/as,bs,"*")
b.per.a <- -t(t(1/as) * bs)
}
return(list(k.per.a=k.per.a,b.per.a=b.per.a))
}
# For each test, a segment in which |test statisitcs| is less than the cutoff stat value given after consideration of max stat among multiple tests.
# For each random vector, a segment is returned as a list of 1x2 matrix.
#' Calculate a segment to be integrated for individual K-cutoff value for each random unit vector
#'
#' @param window.ab.out An output of window.ab
#' @param one.side Optional. Logical. When TRUE, test vectors indicate one-side test and otherwise two-sided
#' @return A list segments is a matrix with two rows, 1st column of which is the start and 2nd is the end of segments; ids is a vector indicating random unit vector ids to which segments belong
#' @examples
#' # utility function being used in other functions
#' @export
window.select <- function(window.ab.out,one.side){
if(missing(one.side))one.side <- FALSE
dm <- dim(window.ab.out$k.per.a)
n.r <- dm[1]
n.k <- dm[3]
n.t <- dm[2]
if(is.null(window.ab.out$b.per.a)){
ret <- list()
for(i in 1:n.k){
tmp <- apply(abs(matrix(window.ab.out$k.per.a[,,i],ncol=n.t)),1,min)
ret[[i]] <- list(segments=cbind(-tmp,tmp),ids=1:length(tmp))
}
}else{
ret <- list()
for(i in 1:n.k){
ret[[i]] <- list()
tmp.list.seg1 <- tmp.list.seg2 <- tmp.list.id <- list()
for(j in 1:n.r){
s <- cbind(c(-window.ab.out$k.per.a[j,,i]),c(window.ab.out$k.per.a[j,,i])) + window.ab.out$b.per.a[j,]
if(one.side){
s <- cbind(-c(window.ab.out$k.per.a[j,,i]),rep(0,length(window.ab.out$k.per.a[j,,i]))) + window.ab.out$b.per.a[j,]
}
tmp.seg <- segment.overlap(s)
tmp.list.seg1[[j]] <- tmp.seg[1]
tmp.list.seg2[[j]] <- tmp.seg[2]
tmp.list.id[[j]] <- rep(j,length(tmp.list.seg1[[j]]))
}
ret[[i]] <- list(segments=cbind(unlist(tmp.list.seg1),unlist(tmp.list.seg2)),ids=unlist(tmp.list.id))
}
}
ret
}
# for a list of segment, probability of df-chisq distribution within the segment is calculated and their mean is returned.
#' Calculate probability indicated by segments in random unit vector directions
#'
#' @param segs.out Output of window.select
#' @param n NUmber of random unit vectors
#' @param df Degrees of freedom
#' @param one.side Optional. Logical. When TRUE, test vectors indicate one-side test and otherwise two-sided
#' @return Probability estimated
#' @examples
#' # utility function being used in other functions
#' @export
calc.pr <- function(segs.out,n,df,one.side){
if(missing(one.side))one.side <- FALSE
ret <- rep(0,length(segs.out))
for(i in 1:length(ret)){
if(one.side){
segs.out[[i]]$segments[which(segs<0)] <- 0
}
tmp1 <- pchisq.bid(segs.out[[i]]$segments[,1],df)
tmp2 <- pchisq.bid(segs.out[[i]]$segments[,2],df)
ret[i] <- (sum(tmp2)-sum(tmp1))/n
}
ret
}
# Power calculator
#' Power calculator
#'
#' @param A An array or matrix of a table representing alternative hypothesis
#' @param tests A list of array/marix of tables indicating tests
#' @param alpha type 1 error threshold
#' @param n Number of iteration to estimate
#' @return
#' \itemize{
#' \item{"pvalues"}{null.p: estimated cumulative probability for null hypothesis when statistic K being x or higher
#' that is estimated for quantile of 1-alpha. alt.p: estimated value for alternative hypotheis }
#' \item{"chisqs"}{Estimated K^2 value as a quantile of 1-alpha}
#' \item{"alpha.given"}{Argument alpha}
#' \item{"alpha.used"}{Estimated p.value for the estimated K^2 value}
#' \item{"power"}{Power}
#' }
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' alpha <- 0.05
#' n <- 1000
#' power.mway.table(A,tests,alpha,n)
#' @export
power.mway.table <- function(A,tests,alpha,n){
if(missing(alpha))alpha <- 0.05
if(missing(n))n<-1000
tmp <- qmway.table(1-alpha,A,tests)
x <- tmp$x
alpha.estimated <- 1-tmp$p.estimated
tmp.2 <- pmway.table(x,A,tests,nc=TRUE,n=n)
null.p <- tmp.2$p.null
alt.p <- tmp.2$p.alt
#cut.off <- which(abs(null.p-(1-alpha))==min(abs(null.p-(1-alpha))))
return(list(p.values = cbind(null.p,alt.p),chisqs = x,alpha.given=alpha,alpha.used=alpha.estimated,power=1-alt.p))
}
# With alpha-cut off, power is returned using two cumul prob for alt.hx.
#' Power calculator with statistic value being fixed
#'
#' @param x K statistic values
#' @param rmway.out Output of rmway function
#' @param alpha Optional. type 1 error threshold
#' @return
#' \itemize{
#' \item{"pvalues"}{null.p: estimated cumulative probability for null hypothesis when statistic K being x or higher
#' that is estimated for quantile of 1-alpha. alt.p: estimated value for alternative hypotheis }
#' \item{"chisqs"}{Estimated K^2 value as a quantile of 1-alpha}
#' \item{"alpha.given"}{Argument alpha}
#' \item{"alpha.used"}{Estimated p.value for the estimated K^2 value}
#' \item{"power"}{Power}
#' }
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' rmway.out <- rmway(n,A,tests)
#' x <- 8
#' alpha <- 0.05
#' power.mway(x,rmway.out,alpha)
#' @export
power.mway <- function(x,rmway.out,alpha){
if(missing(alpha))alpha <- 0.05
null.p <- pmway(x,rmway.out,lower.tail=TRUE)
alt.p <- pmway(x,rmway.out,alt.intersect=TRUE,lower.tail=TRUE)
cut.off <- which(abs(null.p-(1-alpha))==min(abs(null.p-(1-alpha))))
return(list(p.values = cbind(null.p,alt.p),chisqs = x^2,cut.off = x[cut.off]^2,power=1-alt.p[cut.off]))
}
# alt hx is given as intersect (b/a) instead of alt.hx table and cumul pr of chisq value x is returned.
#' Power calculator with intersect information rather than alternative hypothesis table given
#'
#' @param x K statistic values
#' @param rmway.out Output of rmway function
#' @param intersect alternative hypothesis-dependent term
#' @param lower.tail ogical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @return Power
#' @examples
#' # Utility function for the authors
#' @export
pmway.intersect <- function(x,rmway.out,intersect,lower.tail){
if(missing(lower.tail))lower.tail <- TRUE
w.ab.1 <- window.ab(rmway.out$runitstat.out,x,intersect)
w.select.1 <- window.select(w.ab.1)
ret <- calc.pr(w.select.1,rmway.out$n,rmway.out$t.sphere$df)
if(!lower.tail){
ret <- 1 - ret
}
ret
}
#' Generate althernative hypothesis tables in a direction with multipls K values
#'
#' @param t.sphere Output of table.sphere
#' @param tests A list of tables indicating tests
#' @param test.table An array or matrix of table indicating the direction of alternaitve hypothesis
#' @param odds.ratio.table An array or matrix of table indicating 1,2,3,4 cells for calculation of odds ratio to be controled
#' @param ks Optional A vector of K statistic values
#' @return
#' \itemize{
#' \item{"alt.tables"}{A list of alternative tables}
#' \item{"alt.vecs"}{A matrix of alternative table vectors in df-dimensional space}
#' \item{"alt.intersect.out"}{A matrix of terms used to calculate pertinent information}
#' \item{"stats"}{A vector of max K of alternative tables}
#' \item{"difs"}{A matrix each row of which has difference of cells of alternative hypothesis tables}
#' \item{"ors"}{A vector of odds ratios}
#' }
#' @examples
#' # Utility function for the authors
#' @export
model.table.series <- function(t.sphere,tests,test.table,odds.ratio.table,ks){
if(missing(ks))ks <- seq(from=0,to=20,by=1)
n.k <- length(ks)
tmp <- c(test.table) %*% t.sphere$map.matrices$NormalVec2DF
dif <- t.sphere$map.matrices$DF2full %*% t(tmp)
e.table <- t.sphere$tables$expected
test.vecs <- make.test.vecs(t.sphere,tests)
difs <- matrix(0,n.k,length(e.table))
alt.vecs <- matrix(0,n.k,t.sphere$df)
stats <- ors <- rep(0,n.k)
alt.tables <- list()
alt.intersect.out <- matrix(0,n.k,length(test.vecs[,1]))
X <- sum(e.table[which(odds.ratio.table==1)])
Y <- sum(e.table[which(odds.ratio.table==2)])
Z <- sum(e.table[which(odds.ratio.table==3)])
W <- sum(e.table[which(odds.ratio.table==4)])
for(i in 1:n.k){
alt.vecs[i,] <- tmp * ks[i]
alt.intersect.out[i,] <- c(test.vecs %*% matrix(alt.vecs[i,],ncol=1))
stats[i] <- max((test.vecs %*% matrix(alt.vecs[i,],ncol=1))^2)
alt.tables[[i]] <- array(c(e.table)+ t.sphere$map.matrices$DF2full %*% matrix(alt.vecs[i,],ncol=1),t.sphere$r.vec)
difs[i,] <- t.sphere$map.matrices$DF2full %*% matrix(alt.vecs[i,],ncol=1)
tmp.dif <- difs[i,]
x <- sum(tmp.dif[which(odds.ratio.table==1)])
y <- sum(tmp.dif[which(odds.ratio.table==2)])
z <- sum(tmp.dif[which(odds.ratio.table==3)])
w <- sum(tmp.dif[which(odds.ratio.table==4)])
ors[i] <- (X+x)*(W+w)/((Y+y)*(Z+z))
}
return(list(alt.tables=alt.tables,alt.vecs=alt.vecs,alt.intersect.out=alt.intersect.out,stats=stats,difs=difs,ors=ors))
}
#' Generate random tables under unll or alternative hypothesis
#'
#' @param n Number of random tables to be generated
#' @param A An array or matrix of table that specifies marginal counts and also alternative hypothesis if nc is TRUE
#' @param tests Optional. When TRUE, test statistics and maximum among them are calculated and returned
#' @param nc Optional. When FALSE, random tables are in null hypothesis and when TRUE in alternative hypothesis with its ccenter being table A
#' @param one.side Optional. When TRUE, tests are one-sided, otherwise two-sided
#' @return
#' \itemize{
#' \item{"stats.all"}{A matrix of K statistics of random tables}
#' \item{"stats"}{A matrix of K^2 statistics of random tables}
#' \item{"rtables"}{A list of random tables}
#' \item{"ellipsoid.loc"}{A matrix of coordinates of random tables that are in df-dimensional space but the contours are ellipsoids}
#' }
#' @details This function is one of three major functions in this package. Functions for a distribution are rxxxx,dxxxx, pxxxx and qxxxx
#' in general and this package provides three of them, rmway.table, pmway.table and qmway.table functions.
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' n <- 1000
#' rmway.tables.null <- rmway.table(n,A)
#' rmway.tables.alt <- rmway.table(n,A,nc=TRUE)
#' @export
rmway.table <- function(n,A,tests,nc,one.side){
t.sphere <- table.sphere(A)
if(missing(tests))tests <- NULL
if(missing(nc))nc <- FALSE
if(missing(one.side))one.side <- FALSE
if(!nc){
alt.vec <- rep(0,t.sphere$df)
}else{
alt.vec <- t.sphere$map.matrices$Full2DF %*% c(t.sphere$tables$diff)
}
if(!is.null(tests)){
test.vecs <- make.test.vecs(t.sphere,tests)
}
#rmvn <- st.mvn(n,t.sphere$df)
rmvn <- matrix(rnorm(n*t.sphere$df),nrow=n)
loc.df <- t(rmvn) + c(alt.vec)
#loc.full <- t.sphere$matrices$X.inv.sub %*% loc.df
loc.full <- t.sphere$map.matrices$DF2full %*% loc.df
ellipsoid.loc <- t.sphere$matrices$X.sub %*% loc.full
loc.full.2 <- loc.full + c(t.sphere$tables$expected)
rtables <- list()
for(i in 1:n){
rtables[[i]] <- array(loc.full.2[,i],t.sphere$r.vec)
}
if(is.null(tests)){
stats <- NULL
stats.all <- NULL
}else{
stats <- test.vecs %*% loc.df
if(one.side){
stats[which(stats<0)] <- 0
}
stats.all <- stats^2
stats <- apply(stats,2,max)
}
return(list(stats.all = stats.all,stats=stats^2,rtables=rtables,ellipsoid.loc=t(ellipsoid.loc)))
}
#' Estimate p-values/cumulative probability of maximum statistics on multi-way table under unll or alternative hypothesis
#'
#' @param x A vector of K^2 quantiles
#' @param A An array or matrix of table that specifies marginal counts
#' @param tests A list of tables
#' @param lower.tail logical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @param nc Optional. if FALSE (default), random tables are in null hypothesis and otherwise in alternative hypothesis with its ccenter being table A
#' @param one.side Optional. When TRUE, tests are one-sided, otherwise two-sided
#' @param n number of random unit vectors used for estimation
#' @return
#' \itemize{
#' \item{"p.null"}{cumulative probability under null hypothesis}
#' \item{"p.alt"}{cumulative probability under alternative hypothesis}
#' }
#' @details This function is one of three major functions in this package. Functions for a distribution are rxxxx,dxxxx, pxxxx and qxxxx
#' in general and this package provides three of them, rmway.table, pmway.table and qmway.table functions.
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' x <- 0:10
#' pmway.table(x,A,tests)
#' @export
pmway.table <- function(x,A,tests,lower.tail,nc,one.side,n){
if(missing(lower.tail))lower.tail <- TRUE
if(missing(nc))nc <- FALSE
if(missing(one.side))one.side <- FALSE
if(missing(n))n <- 1000
if(!nc){
ret.null <- pmway.table.null(x,A,tests,lower.tail=lower.tail,one.side=one.side,n=n)
ret.alt <- NULL
return(list(p.null=ret.null,p.alt=ret.alt))
}
t.sphere <- table.sphere(A)
if(!nc){
alt.vec <- rep(0,t.sphere$df)
}else{
alt.vec <- t.sphere$matrices$X.sub %*% c(A)
}
test.vecs <- make.test.vecs(t.sphere,tests)
rmvn <- st.mvn(n,t.sphere$df)
runitstat.out <- runitstat(test.vecs,rmvn)
intersect <- c(test.vecs %*% matrix(alt.vec,ncol=1))
w.ab.1 <- window.ab(runitstat.out,sqrt(x),intersect)
w.select.1 <- window.select(w.ab.1,one.side=one.side)
pr.sum <- calc.pr(w.select.1,n,t.sphere$df)
max.stat <- apply(abs(runitstat.out),1,max)
non.zero <- which(max.stat !=0)
k.per.a <- outer(1/max.stat,sqrt(x),"*")
pr.null <- pchisq(k.per.a^2,t.sphere$df)
pr.sum.null <- apply(pr.null,2,sum)/n
if(lower.tail){
if(one.side){
ret <- 1 - (1 - pr.sum)/2
ret.null <- 1 - (1 - pr.sum.null)/2
}else{
ret <- pr.sum
ret.null <- pr.sum.null
}
}else{
if(one.side){
ret <- (1 - pr.sum)/2
ret.null <- (1 - pr.sum.null)/2
}else{
ret <- 1 - pr.sum
ret.null <- 1 - pr.sum.null
}
}
return(list(p.null=ret.null,p.alt=ret))
}
# This is important.
#' Estimate p-values/cumulative probability of maximum statistics on multi-way table under unll hypothesis
#'
#' @param x A vector of K^2 quantiles
#' @param A An array or matrix of table that specifies marginal counts
#' @param tests A list of tables
#' @param lower.tail logical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @param one.side Optional. When TRUE, tests are one-sided, otherwise two-sided
#' @param n number of random unit vectors used for estimation
#' @return cumulative probability under null hypothesis
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' x <- 0:10
#' pmway.table.null(x,A,tests)
#' @export
pmway.table.null <- function(x,A,tests,lower.tail,one.side,n){
if(missing(lower.tail))lower.tail <- TRUE
if(missing(one.side))one.side <- FALSE
if(missing(n))n <- 1000
t.sphere <- table.sphere(A)
test.vecs <- make.test.vecs(t.sphere,tests)
rmvn <- st.mvn(n,t.sphere$df)
runitstat.out <- runitstat(test.vecs,rmvn)
max.stat <- apply(abs(runitstat.out),1,max)
#non.zero <- which(max.stat !=0)
k.per.a <- outer(1/max.stat,sqrt(x),"*")
pr <- pchisq(k.per.a^2,t.sphere$df)
pr.sum <- apply(pr,2,sum)/n
if(lower.tail){
if(one.side){
ret <- 1 - (1 - pr.sum)/2
}else{
ret <- pr.sum
}
}else{
if(one.side){
ret <- (1 - pr.sum)/2
}else{
ret <- 1 - pr.sum
}
}
ret
}
# This is important to construct qmway.table()
#' Grossly estimate quantile value. Utility function for qmway.table function
#'
#' @param p probability
#' @param A An array or matrix of table that specifies marginal counts
#' @param tests A list of tables
#' @param lower.tail logical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @param x a vector of K statistic values
#' @param nc if FALSE (default) null hypothesis, otherwise alternative hypothesis
#' @param one.side Optional. When TRUE, tests are one-sided, otherwise two-sided
#' @param n number of random unit vectors used for estimation
#' @return two quantile values and two corresponding p values between which p and its quantile exist
#' @examples
#' # Utility function
#' @export
qmway.table.pre <- function(p,A,tests,lower.tail,x,nc,one.side,n){
if(missing(lower.tail))lower.tail <- TRUE
if(missing(x))x <- seq(from=0,to=100,by=10)
if(missing(nc))nc <- FALSE
if(missing(one.side))one.side <- FALSE
if(missing(n))n <- 1000
if(!lower.tail){
p <- 1-p
}
p.out <- pmway.table(x,A,tests,lower.tail=lower.tail,nc=nc,one.side=one.side,n=n)$p.null
if(p-min(p.out) < 0){
print("min of x is too big")
return(NULL)
}else if(max(p.out) - p < 0){
print("max of x is too small")
return(NULL)
}
tmp.1 <- which(p.out <= p)
tmp.1 <- tmp.1[length(tmp.1)]
tmp.2 <- which(p.out >= p)
tmp.2 <- tmp.2[1]
return(c(x[tmp.1],x[tmp.2],p.out[tmp.1],p.out[tmp.2]))
}
# qmway.table returns stat value for p given as an argument. The returned value is selected by interporating between the closest stat values by two step interporations.
#' Estimate quantiles of maximum statistics on multi-way table under unll or alternative hypothesis
#'
#' @param p probability vector
#' @param A An array or matrix of table that specifies marginal counts
#' @param tests A list of tables
#' @param lower.tail logical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @param x a vector of K statistic values
#' @param nc Optional. if FALSE (default), random tables are in null hypothesis and otherwise in alternative hypothesis with its ccenter being table A
#' @param one.side Optional. When TRUE, tests are one-sided, otherwise two-sided
#' @param n number of random unit vectors used for estimation
#' @return
#' \itemize{
#' \item{"x"}{quantile}
#' \item{"p.used"}{Argument p}
#' \item{"p.estimated"}{p value corresponding to the estimated quantile, which is some deviated from p of argument because, estimation process is involeved.}
#' }
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' p <- 0.05
#' qmway.table(p,A,tests)
#' @export
qmway.table <- function(p,A,tests,lower.tail,x,nc,one.side,n){
if(missing(lower.tail))lower.tail <- TRUE
if(missing(x))x <- seq(from=0,to=100,by=10)
if(missing(nc))nc <- FALSE
if(missing(one.side))one.side <- FALSE
if(missing(n))n <- 1000
pre.out <- qmway.table.pre(p,A,tests,lower.tail=lower.tail,x=x,nc=nc,one.side=one.side,n=n)
post.x <- seq(from=pre.out[1],to=pre.out[2],length=n)
post.out <- qmway.table.pre(p,A,tests,lower.tail=lower.tail,x=post.x,nc=nc,one.side=one.side,n=n)
ret <- post.out[1] + (post.out[2]-post.out[1]) * (p-post.out[3])/(post.out[4]-post.out[3])
ret.p <- pmway.table(ret,A,tests)$p.null
return(list(x=ret,p.used=p,p.estimated=ret.p))
}
# Tests
#' Max test
#'
#' @param A An array or matrix of table that specifies marginal counts
#' @param tests A list of tables
#' @param lower.tail logical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @param one.side Optional. When TRUE, tests are one-sided, otherwise two-sided
#' @param n number of random unit vectors used for estimation
#' @return
#' \itemize{
#' \item{"statistic"}{maximum among test statistics; K^2-equivalent}
#' \item{"p.value"}{multiple testing-corrected P value}
#' }
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' mway.test(A,tests,n=10000)
#' @export
mway.test <- function(A,tests,lower.tail,one.side,n){
if(missing(lower.tail))lower.tail <- FALSE
if(missing(one.side))one.side <- one.side <- FALSE
if(missing(n))n <- 1000
t.sphere <- table.sphere(A)
test.vecs <- make.test.vecs(t.sphere,tests)
stats <- test.vecs %*% t.sphere$map.matrices$Full2DF %*% c(t.sphere$tables$diff)
max.stats <- max(abs(stats))
if(one.side){
max.stats <- max(stats)
}
p <- pmway.table(max.stats^2,A,tests,lower.tail=lower.tail,one.side=one.side,n=n)
return(list(statistic=max.stats^2,p.value=p))
}
#' Generate data matrix from table
#'
#' @param m An array or matrix of table that specifies marginal counts
#' @return A matrix of sample size x number of variables, levels of which are integers starting from 1
#' @examples
#' m <- matrix(c(10,20,30,40,50,60),2,3)
#' make.data.from.table(m)
#' @export
make.data.from.table <- function(m){
r <- dim(m)
address <- which(m <= Inf,arr.ind=TRUE)
N <- sum(m)
x <- matrix(0,0,length(r))
for(i in 1:length(address[,1])){
tmp <- matrix(rep(address[i,],m[i]),byrow=TRUE,nrow=m[i])
x <- rbind(x,tmp)
}
x
}
#' Generate an array or matrix from data matrix
#'
#' @param d a matrix of sample size x number of variables, levels of which should be integers starting from 1
#' @return An array or matrix
#' @examples
#' m <- matrix(c(10,20,30,40,50,60),2,3)
#' d <- make.data.from.table(m)
#' make.table.from.data(d)
#' @export
make.table.from.data <- function(d){
rg <- apply(d,2,max)
ret <- array(rep(0,prod(rg)),rg)
prod.rg <- cumprod(rg)
prod.rg <- c(1,prod.rg[-length(rg)])
tmp <- apply(t(d-1) * prod.rg,2,sum) + 1
for(i in 1:length(tmp)){
ret[tmp[i]] <- ret[tmp[i]]+1
}
ret
}
#' Randomize a data matrix
#'
#' @param data a matrix of sample size x number of variables
#' @param v a vector of 0 or 1 with length being the number of variables, which indicates variables to be shuffled.
#' @return data matrix
#' @examples
#' m <- matrix(c(10,20,30,40,50,60),2,3)
#' d <- make.data.from.table(m)
#' make.table.from.data(d)
#' @export
perm.data <- function(data,v){
if(missing(v))v <- rep(1,length(data[1,]))
ret <- data
for(i in 1:length(data[1,])){
if(v[i]==1){
ret[,i] <- sample(ret[,i])
}
}
ret
}
# rmway.out has info of alt hx table and otherts.
# x is chisq value and p values (cumulative prob) is returned
# When alt.intersect is TRUE, cumulative prob of x is returned, otherwise cumulative prob of x when null is true.
#' Randomize a data matrix
#'
#' @param x K statistic value
#' @param rmway.out output of rmway function
#' @param alt.intersect logical, if false (default), null hypotheis is assumued, otherwise alternative hypothesis is assumed
#' @param lower.tail logical; if TRUE (default), probabilities are P[X ?????? x], otherwise, P[X > x].
#' @return probability
#' @examples
#' # Utility function
#' @export
pmway <- function(x,rmway.out,alt.intersect,lower.tail){
if(missing(alt.intersect))alt.intersect <- FALSE
if(missing(lower.tail))lower.tail <- TRUE
if(!alt.intersect){
w.ab.1 <- window.ab(rmway.out$runitstat.out,x)
}else if(alt.intersect){
w.ab.1 <- window.ab(rmway.out$runitstat.out,x,rmway.out$alt.intersect.out)
}
w.select.1 <- window.select(w.ab.1)
ret <- calc.pr(w.select.1,rmway.out$n,rmway.out$t.sphere$df)
if(!lower.tail){
ret <- 1 - ret
}
ret
}
# rmway.out is already calculated and given as an argment and alt.hx is given to calculate alt hx table with other info.
#' Generate random K statistic values with rmway function's output reused
#'
#' @param alt.hx A matrix or array of a table which indicates the direction of alternative hypothesis along with an argument odds.ratio.table
#' @param odds.ratio.table A matrix or array of a table with cell values being 0,1,2,3 or 4: Odds ratio of table is indicated by the (sum(cells of 1)*sum(cells of 4))/(sum(cells of 2)*sum(cells of 3))
#' @param or A number indicating odds ratio defined by odds.ratio.table argument
#' @param rmway.out rmway function's output
#' @param n.alt Optional. An integer which makes a one-dimensional grid with n points, among which the table with odds ratio closest to the argument or is selected
#' @param K.alt Optional. A number indicating the maximum value of square of chi-square value to search the table of alternative hypothsis
#' @return
#' \itemize{
#' \item{"t.sphere"}{ Output of table.sphere(A)}
#' \item{"alt.table"}{Argument alt.table itself}
#' \item{"tests"}{Argument tests}
#' \item{"test.vecs"}{Matrix whose rows are test vectors}
#' \item{"runitstat.out"}{Output of runitstat function}
#' \item{"alt.intersect.out"}{Output of alt.intersect function}
#' \item{"n"}{Argument n}
#' \item{"A"}{Argument A}
#' \item{"alt.hx"}{Argument alt.hx}
#' \item{"odds.ratio.table"}{Argument odds.ratio.table}
#' \item{"or"}{Argument or}
#' \item{"n.alt"}{Argument n.alt}
#' \item{"K.alt"}{Argument K.alt}
#' \item{"rmvn"}{A matrix of random unit vectors, each row of which is a random unit vector}
#' }
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' rmway.out <- rmway(n,A,tests)
#' test.table <- matrix(c(1,0.5,0,0,0,0),byrow=TRUE,2,3)
#' odds.ratio.table <- matrix(c(1,0,2,3,0,4),byrow=TRUE,2,3)
#' or <- 2
#' rmway.alt.only(test.table,odds.ratio.table,or,rmway.out)
#' @export
rmway.alt.only <- function(alt.hx,odds.ratio.table,or,rmway.out,n.alt,K.alt){
if(missing(n.alt))n.alt <- 1000
if(missing(K.alt))K.alt <- 10
t.sphere <- rmway.out$t.sphere
alt.table <- make.alternative.table(t.sphere,alt.hx,odds.ratio.table,or,n.alt,K.alt)
test.vecs <- rmway.out$test.vecs
rmvn <- rmway.out$rmvn
runitstat.out <- rmway.out$runitstat.out
alt.intersect.out <- alt.intersect(alt.table,test.vecs)
return(list(t.sphere=t.sphere,alt.table=alt.table,tests=tests,test.vecs=test.vecs,runitstat.out=runitstat.out,alt.intersect.out=alt.intersect.out,n=n,A=rmway.out$tables$data,tests=tests,alt.hx=alt.hx,odds.ratio.table=odds.ratio.table,or=or,n.alt=n.alt,K.alt=K.alt,rmvn = rmvn))
}
# alt hx is defined with alt.hx being an array and target odds ratio defined by odds ratio table with {1,2,3,4} with target or value. n.alt and K.alt are parameters to select alt hx table.
# alt hx table is returned with other infos.
#' Generate random K statistic values
#'
#' @param n Number of random unit vectors
#' @param A An array or matrix of a table
#' @param tests A list of test tables
#' @param alt.hx A matrix or array of a table which indicates the direction of alternative hypothesis along with an argument odds.ratio.table
#' @param odds.ratio.table A matrix or array of a table with cell values being 0,1,2,3 or 4: Odds ratio of table is indicated by the (sum(cells of 1)*sum(cells of 4))/(sum(cells of 2)*sum(cells of 3))
#' @param or A number indicating odds ratio defined by odds.ratio.table argument
#' @param n.alt Optional. An integer which makes a one-dimensional grid with n points, among which the table with odds ratio closest to the argument or is selected
#' @param K.alt Optional. A number indicating the maximum value of square of chi-square value to search the table of alternative hypothsis
#' @return
#' \itemize{
#' \item{"t.sphere"}{ Output of table.sphere(A)}
#' \item{"alt.table"}{Argument alt.table itself}
#' \item{"tests"}{Argument tests}
#' \item{"test.vecs"}{Matrix whose rows are test vectors}
#' \item{"runitstat.out"}{Output of runitstat function}
#' \item{"alt.intersect.out"}{Output of alt.intersect function}
#' \item{"n"}{Argument n}
#' \item{"A"}{Argument A}
#' \item{"alt.hx"}{Argument alt.hx}
#' \item{"odds.ratio.table"}{Argument odds.ratio.table}
#' \item{"or"}{Argument or}
#' \item{"n.alt"}{Argument n.alt}
#' \item{"K.alt"}{Argument K.alt}
#' \item{"rmvn"}{A matrix of random unit vectors, each row of which is a random unit vector}
#' }
#' @examples
#' A <- matrix(c(10,20,30,40,50,60),2,3)
#' t1 <- t2 <- t3 <- matrix(c(1,0,0,0,0,0),byrow=TRUE,2,3)
#' t2[1,2] <- 0.5
#' t3[1,2] <- 1
#' tests <- list(t1,t2,t3)
#' n <- 1000
#' rmway.out <- rmway(n,A,tests)
#' @export
rmway <- function(n,A,tests,alt.hx,odds.ratio.table,or,n.alt,K.alt){
t.sphere <- table.sphere(A)
if(missing(alt.hx)){
alt.hx <- odds.ratio.table <- or <- n.alt <- K.alt <- NULL
alt.table <- list(alt.table=t.sphere$tables$expected,alt.vec=rep(0,t.sphere$df))
}else{
if(missing(n.alt))n.alt <- 1000
if(missing(K.alt))K.alt <- 10
alt.table <- make.alternative.table(t.sphere,alt.hx,odds.ratio.table,or,n.alt,K.alt)
}
test.vecs <- make.test.vecs(t.sphere,tests)
rmvn <- st.mvn(n,t.sphere$df)
runitstat.out <- runitstat(test.vecs,rmvn)
alt.intersect.out <- alt.intersect(alt.table,test.vecs)
return(list(t.sphere=t.sphere,alt.table=alt.table,tests=tests,test.vecs=test.vecs,runitstat.out=runitstat.out,alt.intersect.out=alt.intersect.out,n=n,A=A,tests=tests,alt.hx=alt.hx,odds.ratio.table=odds.ratio.table,or=or,n.alt=n.alt,K.alt=K.alt,rmvn = rmvn))
}
# A??????????????????0??????????????????????????????????????????(???????????????)
# A???????????????????????????(?????????????????????)
# A.ori???NULL???????????????????????????A.ori??????????????????????????????????????????????????????A??????????????????????????????
#' Utility function to handle tables with cells whose expected value is 0
#'
#' @param A an array or matrix of table, that may be the table with cells whose expected value is zero (0 cells) or
#' that may be the table without such cells
#' @param A.ori if missing(default), A represetnts a table with 0 cell(s) and the function converts A to a table without 0 cell(s).
#' Otherwise, A should be a table without 0 cell(s) and A.ori is the original array/matrix before shrinkage.
#' of table with 0 cells
#' @return An array or matrix of table. if A.ori is missing, shrunken table is returned otherwise a table with 0 cells recovered.
#' with 0 cells.
#' @examples
#' n <- 5
#' M <- list()
#' N <-100
#' for(i in 1:n){
#' M[[i]] <- sample(0:5,sample(2:4,1))
#' M[[i]] <- M[[i]]/sum(M[[i]])*N
#' }
#' tmp <- M[[1]]
#' for(i in 2:length(M)){
#' tmp <- t(M[[i]]/N) %x% tmp
#' }
#' A <- array(tmp,lapply(M,length))
#'
#' A.r <- reshape.table(A)
#' A.2 <- reshape.table(A.r,A)
#' calc.marg(A)
#' calc.marg(A.r)
#' calc.marg(A.2)
#' A-A.2
#' @export
reshape.table <- function(A,A.ori){
if(missing(A.ori)){
A.ori <- NULL
}
inv <- FALSE
if(!is.null(A.ori)){
inv <- TRUE
}
if(!inv){
A.ori <- A
}
M <- calc.marg(A.ori)
S <- lapply(M,sign)
new.dim <- lapply(S,sum)
tmp <- S[[1]]
for(i in 2:length(S)){
tmp <- t(S[[i]]) %x% tmp
}
if(!inv){
return(array(A[which(tmp!=0)],new.dim))
}else{
tmp.A <- rep(0,length(A.ori))
tmp.A[which(tmp!=0)] <- c(A)
return(array(tmp.A,dim(A.ori)))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.