R/ICQTNullDist.R

In multtest: Resampling-based multiple hypothesis testing

# No robust correlation test statistics.
# Want to return a 3 by M matrix of observations.
corr.Tn <- function(X,test,alternative,use="pairwise"){
  P <- dim(X)[1]
  M <- P*(P-1)/2
  N <- dim(X)[2]
  VCM <- cov(t(X),use=use)
  Cor <- cov2cor(VCM)
  Cov.v <- VCM[lower.tri(VCM)] # vectorize.
  Cor.v <- Cor[lower.tri(Cor)] # vectorize.
  if(test=="t.cor") num <- sqrt(N-2)*Cor.v/sqrt(1-Cor.v^2)
  if(test=="z.cor") num <- sqrt(N-3)*0.5*log((1+Cor.v)/(1-Cor.v))
  denom <- 1
  if(alternative=="two.sided"){
			snum<-sign(num)
		 	num<-abs(num)
                      }
  else {
    if(alternative=="less"){
      snum<-(-1)
      num<-(-num)
    }
    else snum<-1
    }
  rbind(num,denom,snum)
}

ic.tests <- c("t.onesamp","t.pair","t.twosamp.equalvar","t.twosamp.unequalvar","lm.XvsZ","lm.YvsXZ","t.cor","z.cor")

corr.null <- function(X,W=NULL,Y=NULL,Z=NULL,test="t.twosamp.unequalvar",alternative="two-sided",use="pairwise",B=1000,MVN.method="mvrnorm",penalty=1e-6,ic.quant.trans=FALSE,marg.null=NULL,marg.par=NULL,perm.mat=NULL){
  # Most sanity checks conducted already...
  p <- dim(X)[1]
  m <- dim(X)[1] 
  n <- dim(X)[2] 
  cat("calculating vector influence curve...", "\n")

  if(test=="t.onesamp" | test=="t.pair"){
    #t.pair sanity checks and formatting done in stat.closure section
    #in test.R
    if(is.null(W)) IC.Cor <- cor(t(X),use=use)
    else IC.Cor <- IC.CorXW.NA(X,W,N=n,M=p,output="cor")
  }

  if(test=="t.twosamp.equalvar" | test=="t.twosamp.unequalvar"){
    uY<-sort(unique(Y))
    if(length(uY)!=2) stop("Must have two class labels for this test")
    n1 <- sum(Y==uY[1])
    n2 <- sum(Y==uY[2])
    if(is.null(W)){
      cov1 <- cov(t(X[,Y==uY[1]]),use=use)
      cov2 <- cov(t(X[,Y==uY[2]]),use=use)
    }
    else{
      cov1 <- IC.CorXW.NA(X[,Y==uY[1]],W[,Y==uY[1]],N=n1,M=p,output="cov")
      cov2 <- IC.CorXW.NA(X[,Y==uY[2]],W[,Y==uY[2]],N=n2,M=p,output="cov")
    }
    newcov <- cov1/n1 + cov2/n2
    IC.Cor <- cov2cor(newcov)
  }

  # Regression ICs written to automatically incorporate weights.
  # If W=NULL, then give equal weights.
  if(test=="lm.XvsZ"){
    if(is.null(Z)) Z <- matrix(1,nrow=n,ncol=1)
    else Z <- cbind(Z,1)
    if(is.null(W)) W <- matrix(1/n,nrow=p,ncol=n)
    IC.i <- matrix(0,nrow=m,ncol=n)
    for(i in 1:m){
      drop <- is.na(X[i,]) | is.na(rowSums(Z)) | is.na(W[i,])
      x <- as.numeric(X[i,!drop])
      z <- Z[!drop,]
      w <- W[i,!drop]
      nn <- n-sum(drop)
      EXtWXinv <- solve(t(z)%*%(w*diag(nn))%*%z)*sum(w)
      res.m <- lm.wfit(z,x,w)$res
      if(sum(drop)>0) res.m <- insert.NA(which(drop==TRUE),res.m)
      EXtWXinvXt <- rep(0,n)
      for(j in 1:n){
        EXtWXinvXt[j] <- (EXtWXinv%*%(t(Z)[,j]))[1]
      }
      IC.i[i,] <- res.m * EXtWXinvXt
    }
    IC.Cor <- IC.Cor.NA(IC.i,W,N=n,M=p,output="cor")
  }
  
  if(test=="lm.YvsXZ"){
    if(is.null(Y)) stop("An outcome variable is needed for this test")
    if(length(Y)!=n) stop(paste("Dimension of outcome Y=",length(Y),", not equal dimension of data=",n,sep=""))
    if(is.null(Z)) Z <- matrix(1,n,1)
    else Z <- cbind(Z,1)
    if(is.null(W)) W <- matrix(1,nrow=p,ncol=n)
    IC.i <- matrix(0,nrow=m,ncol=n)
    for(i in 1:m){
      drop <- is.na(X[i,]) | is.na(rowSums(Z)) | is.na(W[i,])
      x <- as.numeric(X[i,!drop])
      z <- Z[!drop,]
      w <- W[i,!drop]
      y <- Y[!drop]
      nn <- n-sum(drop)
      xz <- cbind(x,z)
      XZ <- cbind(X[i,],Z)
      EXtWXinv <- solve(t(xz)%*%(w*diag(nn))%*%xz)*sum(w)
      res.m <- lm.wfit(xz,y,w)$res
      if(sum(drop)>0) res.m <- insert.NA(which(drop==TRUE),res.m)
      EXtWXinvXt <- rep(0,n)
      for(j in 1:n){
        EXtWXinvXt[j] <- (EXtWXinv%*%(t(XZ)[,j]))[1]
      }
      IC.i[i,] <- res.m * EXtWXinvXt
    }
    IC.Cor <- IC.Cor.NA(IC.i,W,N=n,M=p,output="cor")
  }
  
  if(test=="t.cor" | test=="z.cor"){
    if(!is.null(W)) warning("Weights not currently implemented for tests of correlation parameters.  Proceeding with unweighted version")
    # Change of dimension
    P <- dim(X)[1] -> p # Number of variables.
    M <- P*(P-1)/2 -> m # Actual number of pairwise hypotheses.
    N <- dim(X)[2] -> m
    ind <- t(combn(P,2))
    VCM <- cov(t(X),use="pairwise")
    Cor <- cov2cor(VCM)
    Vars <- diag(VCM)
    Cov.v <- VCM[lower.tri(VCM)] # vectorize.
    Cor.v <- Cor[lower.tri(Cor)] # vectorize.
    X2 <- X*X
    EX <- rowMeans(X,na.rm=TRUE)
    E2X <- rowMeans(X2,na.rm=TRUE)
    Var1.v <- Vars[ind[,1]]
    Var2.v <- Vars[ind[,2]]
    EX1.v <- EX[ind[,1]]
    EX2.v <- EX[ind[,2]]
    E2X1.v <- E2X[ind[,1]]
    E2X2.v <- E2X[ind[,2]]
    X.vec1 <- X[ind[,1],]
    X.vec2 <- X[ind[,2],]
    X.vec12 <- X.vec1*X.vec2
    EX1X2.v <- rowMeans(X.vec12,na.rm=TRUE)

    cons <- 1/sqrt(Var1.v*Var2.v)
    gradient <- matrix(1,nrow=M,ncol=5)
    gradient[,1] <- EX1.v*Cov.v/Var1.v - EX2.v
    gradient[,2] <- EX2.v*Cov.v/Var2.v - EX1.v
    gradient[,3] <- -0.5*Cov.v/Var1.v
    gradient[,4] <- -0.5*Cov.v/Var2.v

    IC.i <- matrix(0, nrow=M, ncol=N)
    for(i in 1:N){
      diffs.i <- diffs.1.N(X[ind[,1],i], X[ind[,2],i], EX1.v, EX2.v, E2X1.v, E2X2.v, EX1X2.v)
      IC.M <- rep(0,M)
      for(j in 1:M){
        IC.M[j] <- gradient[j,]%*%diffs.i[,j]
      }
      IC.i[,i] <- IC.M
    }
    IC.i <- cons * IC.i
    IC.Cor <- IC.Cor.NA(IC.i,W=NULL,N=n,M=M,output="cor")
  }

  if(ic.quant.trans==FALSE) cat("sampling null test statistics...", "\n\n")
  else cat("sampling null test statistics...", "\n")
  
  if(MVN.method=="mvrnorm") nulldist <- t(mvrnorm(n=B,mu=rep(0,dim(IC.Cor)[1]),Sigma=IC.Cor))
  if(MVN.method=="Cholesky"){
    IC.chol <- t(chol(IC.Cor+penalty*diag(dim(IC.Cor)[1])))
    norms <- matrix(rnorm(B*dim(IC.Cor)[1]),nrow=dim(IC.Cor)[1],ncol=B)
    nulldist <- IC.chol%*%norms
  }
  if(ic.quant.trans==TRUE){
    cat("applying quantile transform...", "\n\n")
    if(is.null(marg.null)){
	marg.null <- "t"
     	if(test=="t.cor" | test=="z.cor" | test=="t.twosamp.equalvar") marg.par <- matrix(rep(dim(X)[2]-2,dim(IC.Cor)[1]),nrow=dim(IC.Cor)[1],ncol=1)
        if(test=="lm.XvsZ") marg.par <- matrix(rep(dim(X)[2]-dim(Z)[2],dim(IC.Cor)[1]),nrow=dim(IC.Cor)[1],ncol=1)
        if(test=="lm.YvsXZ")  marg.par <- matrix(rep(dim(X)[2]-dim(Z)[2]-1,dim(IC.Cor)[1]),nrow=dim(IC.Cor)[1],ncol=1)
        else marg.par <- matrix(rep(dim(X)[2]-1,dim(IC.Cor)[1]),nrow=dim(IC.Cor)[1],ncol=1)
    }
    if(test=="z.cor" & marg.null=="t") warning("IC nulldist for z.cor already MVN. Transforming to N-2 df t marginal distribution not advised.")
    if(marg.null!="t" & marg.null!="perm") stop("IC nulldists can only be quantile transformed to a marginal t-distribution or user-supplied marginal permutation distribution")
    if(marg.null=="t") nulldist <- tQuantTrans(nulldist,marg.null="t",marg.par,ncp=0,perm.mat=NULL)
    if(marg.null=="perm") nulldist <- tQuantTrans(nulldist,marg.null="perm",marg.par=NULL,ncp=NULL,perm.mat=perm.mat)
  }
  if(alternative=="greater") nulldist <- nulldist
  else if(alternative=="less") nulldist <- -nulldist
  else nulldist <- abs(nulldist)
  nulldist
}

# Function, given ICs for each individual, returns variance covariance
# matrix or corresponding correlation matrix.
IC.Cor.NA <- function(IC,W,N,M,output){
  n <- dim(IC)[2]
  m <- dim(IC)[1]
  if(is.null(W)){
    W <- matrix(1,nrow=dim(IC)[1],ncol=dim(IC)[2])
    Wnew <- W/rowSums(W,na.rm=TRUE) # Equal weight, NA handling.
  }
  else Wnew <- W/rowSums(W,na.rm=TRUE)
  IC.VC <- matrix(0,nrow=m,ncol=m)
  for(i in 1:n){
    temp <- crossprod(t(sqrt(Wnew[,i])*IC[,i]))
    temp[is.na(temp)] <- 0
    IC.VC <- IC.VC + temp
  }
  if(output=="cov") out <- IC.VC
  if(output=="cor") out <- cov2cor(IC.VC)
  out
}
 
# Weighted correlation. Generalizes cov.wt() to account for a matrix
# of weights. Uses IC formulation instead of sweep() and crossprod().
# May be slower/clunkier, but pretty transparent, and allows for NA
# handling much like cor(...,use="pairwise") would.  That is, each
# element of the correlation matrix returned uses the maximum amount
# of information possible in obtaining individual elements of that
# matrix.
IC.CorXW.NA <- function(X,W,N,M,output){
  n <- dim(X)[2]
  m <- dim(X)[1]
  XW <- X*W
  EXW <- rowSums(XW)/rowSums(W)
  ICW.i <- X-EXW
  Wnew <- W/rowSums(W,na.rm=T)
  IC.VC <- matrix(0,nrow=m,ncol=m)
  for(i in 1:n){
    temp <- crossprod(t(sqrt(Wnew[,i])*X[,i]))
    temp[is.na(temp)] <- 0
    IC.VC <- IC.VC + temp
  }
  if(output=="cov") out <- IC.VC
  if(output=="cor") out <- cov2cor(IC.VC)
  out
}

# For regression ICs, a function to insert NAs into appropriate locations
# of a vector of returned residuals.
insert.NA <- function(orig.NA, res.vec){
  for(i in 1:length(orig.NA)){
    res.vec <- append(res.vec, NA, after=orig.NA[i]-1)
  }
  res.vec
}

# For correlation ICS, a function to get diff vectors for all M.
# This is the difference between estimates for
# a sample size of one and a sample of size n.
diffs.1.N <- function(vec1, vec2, e1, e2, e21, e22, e12){
  diff.mat.1.N <- matrix(0,nrow=5,ncol=length(vec1))
  diff.mat.1.N[1,] <- vec1 - e1
  diff.mat.1.N[2,] <- vec2 - e2
  diff.mat.1.N[3,] <- vec1*vec1 - e21
  diff.mat.1.N[4,] <- vec2*vec2 - e22
  diff.mat.1.N[5,] <- vec1*vec2 - e12
  diff.mat.1.N
}

### For quantile transform, take a sample from the marginal null distribution.
marg.samp <- function(marg.null,marg.par,m,B,ncp){
out <- matrix(0,m,B)
for(i in 1:m){
  if(marg.null=="normal") out[i,] <- rnorm(B,mean=marg.par[i,1],sd=marg.par[i,2])
  if(marg.null=="t") out[i,] <- rt(B,df=marg.par[i,1],ncp)
  if(marg.null=="f") out[i,] <- rf(B,df1=marg.par[i,1],df2=marg.par[i,2],ncp)
}
out
}

### Quantile transform streamlined for IC nulldists.
tQuantTrans <- function(rawboot, marg.null, marg.par, ncp, perm.mat=NULL){
  m <- dim(rawboot)[1]
  B <- dim(rawboot)[2] 
  ranks <- t(apply(rawboot,1,rank,ties.method="random"))
  if(marg.null=="t") Z.quant <- marg.samp(marg.null="t",marg.par,m,B,ncp)
  if(marg.null=="perm") Z.quant <- perm.mat
  Z.quant <- t(apply(Z.quant,1,sort))
  if(marg.null!="perm"){                   
      for(i in 1:m){                         
        Z.quant[i,] <- Z.quant[i,][ranks[i,]]
      }
    }
  else{
    Z.quant <- t(apply(Z.quant,1,quantile,probs=seq(0,1,length.out=B),na.rm=TRUE))
      for(i in 1:m){                         
        Z.quant[i,] <- Z.quant[i,][ranks[i,]]
      }
    }
  Z.quant
}

### Effective df for two sample test of means, unequal var.
t.effective.df <- function(X,Y){
  uY<-sort(unique(Y))
  X1 <- X[Y==uY[1]]
  X2 <- X[Y==uY[2]]
  mu <- var(X2)/var(X1)
  n1 <- length(Y[Y==uY[1]])
  n2 <- length(Y[Y==uY[2]])
  df <- (((1/n1)+(mu/n2))^2)/(1/((n1^2)*(n1-1)) + (mu^2)/((n2^2)*(n2-1)))
  df
}