.packageName <- "metaGEM"
##########################################
### Calculate and combine within study ###
##########################################


effect.sizes <- function(study){

    check.metaGEM.study(study)

    summ <- t( apply( study$expr, 1, getES, g=study$class ) )
    summ <- data.frame( summ, keys=study$keys )
    return(summ)
}


getES <- function(v, g){
    ## function to calculate basic statistics for two-class comparison for a gene

    stopifnot( identical( length(v), length(g) ) )

    x <- cleanNA( v[ which(g==1) ] )
    y <- cleanNA( v[ which(g==0) ] )

    n1 <- length(x); n2 <- length(y)
    if( n1 < 2 | n2 < 2 )
        return( c(n1=NA, m1=NA, sd1=NA,
                  n2=NA, m2=NA, sd2=NA,
                  diff=NA, pooled.sd=NA,
                  g=NA, se.g=NA) )

    m1   <- mean(x); m2 <- mean(y)
    diff <- m1 - m2

    sd1  <- sd(x);  sd2 <- sd(y)
    sp   <- sqrt( ( (n1-1)*sd1^2 + (n2-1)*sd2^2 )/( n1 + n2 - 2 ) )

    cf   <- 1 - 3/( 4*(n1 + n2) - 9 )
    g    <- cf * diff/sp
    se.g <- sqrt( (n1+n2)/(n1*n2) + 0.5*g^2 /(n1+n2-3.94) )

    return( c(n1=n1, m1=m1, sd1=sd1,
              n2=n2, m2=m2, sd2=sd2,
              diff=diff, pooled.sd=sp,
              g=g, se.g=se.g) )

}


## getES2 <- function(Tmat, Nmat){
##
##     mat <- cbind(Tmat, Nmat)
##     grp <- rep( 1:0, c(ncol(Tmat), ncol(Nmat)) )
##     out <- t( apply(mat, 1, getES, g=grp) )
##     return(out)
## }



summ.eff.within <- function(effects, option="fixed.iv"){
    ## summarizing multiple probes that map to the same UniGene within a study

    stopifnot( all( c("g", "se.g", "keys") %in% colnames(effects) ) )
    effects <- effects[ , c("g", "se.g", "keys")]
    effects$keys  <- as.character(effects$keys)

    if( length( grep(",", effects$keys) ) > 0 ) stop("Multiple keys detected. Please expand geneID first")

    if(nrow(effects)==1){ ## hack when only effect size is available per study
        rownames(effects) <- effects$keys
        effects <- effects[ , c("g", "se.g")]
        return(effects)
    }

    ## Deal with singletons first ##
    singles.keys <- names(which(table(effects$keys) == 1))
    singles.ind <- which(effects$keys %in% singles.keys)
    out <- effects[ singles.ind, ]
    rownames(out) <- out$keys;  out$keys <- NULL

    ## Next, deal with multiple keys within a study ##
    multis <- effects[ -singles.ind, ]
    multis$abs.z <- abs( multis$g/multis$se.g ) # used if extreme option is choosen

    if(nrow(multis) > 0){  # skip if no multiple ID found

        tmp <- split(multis, multis$keys)

        if (option == "fixed.iv") {
            out2 <- sapply(tmp, function(m) {
                unlist(meta.summaries(m$g, m$se.g, method = "fixed")[c("summary", "se.summary")])
            })
            out2 <- t(out2)
            colnames(out2) <- c("g", "se.g")
        }

        if (option == "extreme") {
            out2 <- lapply(tmp, function(mat) mat[which.max(mat$abs.z), ])
            out2 <- do.call(rbind, out2)
            out2 <- out2[, c("g", "se.g")]
        }
        out <- rbind(out, out2)
    }

    out <- out[sort(rownames(out)), ]
    return(out)
}



##############################
### Combine across studies ###
##############################

pool.inverseVar <- function( g, se.g, method ){

    stopifnot( identical( rownames(g), rownames(se.g) ) )
    out <- matrix( nr=nrow(g), nc=5,
                  dimnames=list( rownames(g), c("n.studies", "summary", "se.summary", "tau2", "p.value") ) )

    for(j in 1:nrow(g)){

        e  <- cleanNA(    g[j, ] )
        se <- cleanNA( se.g[j, ] )
        n  <- length(e)

        if(n==1){
            summ <- e;   se.summ <- se;   tau2 <- NA
        } else {
            fit <- meta.summaries(e, se, method = method)
            summ <- fit$summary
            se.summ <- fit$se.summary
            tau2 <- ifelse( method=="fixed", NA, fit$tau2 )
            rm(fit)
        }

        pval     <- 2*pnorm( abs(summ/se.summ), lower.tail=FALSE )

        out[j, ] <- c(n, summ, se.summ, tau2, pval)
        rm(e, se, n, summ, se.summ, tau2, pval)
    }
    return(out)
}


combine.effect.sizes <- function (list.of.effects, between.method="random", within.method="fixed.iv", everything=TRUE){

    if( is.null(names(list.of.effects)) )
        names(list.of.effects) <- paste("data", 1:length(list.of.effects), sep="")

    study.effects <- lapply(list.of.effects, function(effects) {

        effects <- data.frame(effects)
        effects$keys <- as.character(effects$keys)

        ## remove probes that cannot be mapped or have insufficient observations to calculate effect size
        bad <- which( is.na(effects$g) | is.na(effects$keys) | effects$keys=="NA" )
        effects <- effects[ setdiff(1:nrow(effects), bad), ]

        ## expand probes that maps to multiple keys
        effects <- expand.df( effects )

        ## summarize multiple probes within a study
        effects <- summ.eff.within(effects, option = within.method)
    })


    tmp <- multimerge(study.effects)
    g    <- tmp[, paste(names(study.effects), "_g", sep = ""), drop=FALSE]
    se.g <- tmp[, paste(names(study.effects), "_se.g", sep = ""), drop=FALSE]

    pooled.estimates <- data.frame( pool.inverseVar(g, se.g, method=between.method ) )

    if (everything) {
        return(list(g=g, se.g=se.g, pooled.estimates=pooled.estimates))
    } else {
        return(pooled.estimates)
    }
}



#####################
### Visualization ###
#####################

plot.sumsum <- function (db, signif.thr=1, K=NULL, labels=NULL, cex=1, lwd=1) {

    ## db <- db$pooled.estimates
    p.adj <- p.adjust(db[, "p.value"], method = "fdr")
    summ <- db[, "summary"]
    se.summ <- db[, "se.summary"]
    if( is.null(labels) ) labels <- rownames(db)
    LCL <- summ - qnorm(0.975) * se.summ
    UCL <- summ + qnorm(0.975) * se.summ
    yrange <- blow(c(LCL, UCL))

    plot(yrange, range(p.adj), type = "n", log = "y", ann = F )
    title(ylab = list("FDR adjusted p-value", font = 4, col = 4),
          xlab = list("Pooled Standardized Mean Difference (log2 scale)",
          font = 4, col = 4), line = 2)
    grid()
    arrows(LCL, p.adj, UCL, p.adj, code=3, length=0.05, angle=45, col=8, lwd=lwd)
    points(summ, p.adj, pch = 20, cex = 0.75, col = 1)
    abline(v = 0, col = 4, lwd = 1)

    if(is.null(K)){
        neg <- which(p.adj < signif.thr & summ < 0)
        pos <- which(p.adj < signif.thr & summ > 0)
        abline(h = signif.thr, col = 4, lwd = 2, lty = 2)
    } else {
        s <- split( p.adj, sign(summ) )
        r <- sapply( s, rank )
        rr <- unsplit( r, sign(summ) )

        pos <- which( summ > 0 & rr < K )
        neg <- which( summ < 0 & rr < K )
    }

    if (length(neg) > 0)
        text(LCL[neg] - 0.125, p.adj[neg], labels = labels[neg], col = "green", cex = cex, font = 4)
    if (length(pos) > 0)
        text(UCL[pos] + 0.125, p.adj[pos], labels = labels[pos], col = "red", cex = cex, font = 4)
}



forest.plot <- function(db, key, main=key, ...){
    stopifnot(identical(rownames(db$g), rownames(db$se.g)))
    stopifnot(identical(rownames(db$g), rownames(db$pooled.estimates)))

    g    <- cleanNA( db$g[ key, ] )
    g    <- g[ sort(names(g)) ]

    se.g <- cleanNA( db$se.g[ key, ] )
    se.g <- se.g[ sort(names(se.g)) ]

    ## stopifnot( identical( names(g), names(se.g) ) )

    pool    <- db$pooled.estimates[ key, "summary" ]
    se.pool <- db$pooled.estimates[ key, "se.summary" ]

    metaplot( g, se.g, labels=names(g),
             summn=pool, sumse=se.pool, sumnn=1/se.pool^2,
             xlab="Standardized Mean Difference (log2 scale)", ylab="", main=main,
             colors=meta.colors(box="blue", lines="lightblue",
             zero="grey", summary="orange", text="red"), ... )
}

##########################################
### Calculate and combine within study ###
##########################################


ttest.Qvalues <- function(study){

    out <- ttest.Pvalues(study)
    out$Q.up   <- p.adjust( out$P.up, method="fdr" )
    out$Q.down <- p.adjust( out$P.down, method="fdr" )

    out <- out[ , c("Q.up", "Q.down", "keys")]
    return(out)
}

ttest.Pvalues <- function(study){

    check.metaGEM.study(study)

    summ <- get.ttest.P( study$expr, study$class )[ , c("P.up", "P.down")]
    summ <- data.frame( summ, keys=study$keys )
    return(summ)
}

get.ttest.P <- function(mat, g){

    ## test statistic and DF calculated using equal variance assumption
    tstat <- mt.teststat( mat, g, test="t.equalvar" )
    df <- length(g) - 2

    P.both <- 2*pt( abs(tstat), df=df, lower=FALSE )
    P.down <- pt( tstat, df=df, lower=TRUE )
    P.up   <- pt( tstat, df=df, lower=FALSE )

    out <- cbind(P.both, P.down, P.up)
    rownames(out) <- rownames(mat)
    return(out)
}


##############################
### Combine across studies ###
##############################

combine.significances <- function(list.of.sigs){

    study.sigs <- lapply(list.of.sigs, function(sigs){
        sigs <- data.frame(sigs)
        sigs$keys <- as.character(sigs$keys)

        ## remove probes that cannot be mapped or have insufficient observations to calculate the p-value
        bad <- which( is.na(sigs[ ,1]) | is.na(sigs$keys) | sigs$keys=="NA" )
        sigs <- sigs[ setdiff(1:nrow(sigs), bad), ]

        ## expand probes that maps to multiple keys
        sigs <- expand.df( sigs )

        ## summarize multiple probes within a study
        out <- summ.sigs.within(sigs)
    })

    tmp <- multimerge(study.sigs)
    sigs.up   <- tmp[ , grep("\\.up$", colnames(tmp), v=TRUE), drop=FALSE]
    sigs.down <- tmp[ , grep("\\.down$", colnames(tmp), v=TRUE), drop=FALSE]

    return(list(sigs.up=sigs.up, sigs.down=sigs.down))
}


count.votes <- function(list.of.sigs, vote.threshold){

    combsigs  <- combine.significances( list.of.sigs )

    sigs.up   <- combsigs$sigs.up
    valid.up <- rowSums( !is.na(sigs.up) )
    yes.up   <- rowSums( sigs.up < vote.threshold, na.rm=TRUE )

    sigs.down <- combsigs$sigs.down
    valid.down <- rowSums( !is.na(sigs.down) )
    yes.down   <- rowSums( sigs.down < vote.threshold, na.rm=TRUE )

    out <- cbind(yes.down, valid.down, yes.up, valid.up)
    return(out)
}

sum.of.logs <- function(list.of.sigs){

    combsigs  <- combine.significances( list.of.sigs )

    sigs.up   <- combsigs$sigs.up
    valid.up  <- rowSums( !is.na(sigs.up) )
    F.stat.up <- -2*rowSums( log(sigs.up), na.rm=TRUE )
    F.pval.up <- pchisq( F.stat.up, 2*valid.up, lower.tail=FALSE )

    sigs.down   <- combsigs$sigs.down
    valid.down  <- rowSums( !is.na(sigs.down) )
    F.stat.down <- -2*rowSums( log(sigs.down), na.rm=TRUE )
    F.pval.down <- pchisq( F.stat.down, 2*valid.down, lower.tail=FALSE )

    out <- cbind(F.stat.up, F.pval.up, F.stat.down, F.pval.down)
    return(out)
}

summ.sigs.within <- function(sigs){

    pval.cols <- setdiff( 1:ncol(sigs), grep("keys", colnames(sigs)) )
    out <- sapply( sigs[ , pval.cols], function(x)
                                           tapply(x, sigs$keys, min) )
    return(out)
}


library(rmeta)
library(multtest)

se <- function(x) sd(x)/length(x)

cleanNA <- function(x) return( x[!is.na(x) & is.finite(x) ] )

gsub.formula <- function(pattern, replacement, x, ...)
                 as.formula( gsub(pattern, replacement, as.expression(x), ...) )


## length.NA <- function(x) length( cleanNA(x) )
## check with length( x, na.rm=TRUE )


install.packages.BioC <- function(x, ...){
    if( !exists("getBioC") ) source("http://www.bioconductor.org/getBioC.R")
    getBioC(x, ...)
}


blow <- function(x, scale=0.10)
    c( mini = ( min(x) - scale*abs(min(x)) ),
      maxi = ( max(x) + scale*abs(max(x)) ) )
## superseeded by extendrange()


expand.df <- function (df, key.name = "keys", keys.sep = ",") {
    ## every row in df may have multiple identities
    ## these identities are defined in keys, with individual keys separated by keys.sep
    ## this function expands every row that has 1:n mapping to n x 1:1 rows

    keys <- as.character(df[ ,key.name])
    skey <- strsplit( keys, split = keys.sep)
    df   <- df[rep(1:nrow(df), sapply(skey, length)), ]
    df[, key.name] <- unlist(skey)
    return(df)
}


multimerge <- function(mylist){
    ## iterative version of merge using all=T option
    ## it assumes rownames is unique and is the common key

    unames <- unique( unlist( lapply( mylist, rownames ) ) )
    n      <- length(unames)

    out <- lapply( mylist, function(df){
        tmp <- matrix( nr=n, nc=ncol(df),
                      dimnames=list( unames, colnames(df) ) )
        tmp[ rownames(df), ] <- as.matrix(df)
        return(tmp)
    })

    bigout <- do.call( cbind, out )
    colnames(bigout) <- paste(rep( names(mylist), sapply(mylist, ncol) ),
                              sapply(mylist, colnames), sep="_")
    return(bigout)
}


sample.stratified <- function( cl , ndg, collapse=TRUE ){

    ## cl  = original classification (can handle any number of classes)
    ## ndg = number of desired groups

    ## Split and then stratify within split
    f   <- function(v) split( sample(v), rep(1:ndg, len=length(v)) )
    tmp <- lapply( split( 1:length(cl), cl ), f )

    ## try to minimize unbalanced groups
    alternate <- function(v) v[ seq(1, length(v), by=2) ]
    toReverse <- alternate( which( (table(cl) %% ndg) != 0 ) )
    for(i in toReverse) tmp[[i]] <- rev( tmp[[i]] )


    ## rearrange everything so that indices are for each allocated group
    out <- vector("list", ndg)
    for(i in 1:ndg)
        out[[i]] <- unlist( lapply(tmp, function(mylist) mylist[[i]]) )

    return(out)

    ## EXAMPLE USAGE
    ## cl <- rep( LETTERS[1:3], c(11, 17, 28) )
    ## sample.stratified( cl, ndg=5 )
}


square.plot <- function(x, y, ...)
    plot( x, y, xlim=range(x,y), ylim=range(x,y), ... )


pairwise.apply <- function(x, FUN, ...){
    ## output of FUN must be scalar

    n <- nrow(x)
    r <- rownames(x)
    output <- matrix(NA, nc=n, nr=n, dimnames=list(r, r))


    for(i in 1:n){
        for(j in 1:n){
            if(i >= j) next()
            output[i, j] <- FUN( x[i,], x[j,], ... )
        }
    }
    return(output)
}



####################################################
## Functions for formatting output and extracting ##
####################################################
se.coefficients <- function(obj){

  cn <- colnames( model.matrix( obj ) )
  x <- rep( NA, length=length(cn) )
  names(x) <- cn

  se.raw <- sqrt( diag( vcov(obj) ) )
  x[ names(se.raw) ] <- se.raw
  return(x)
}

coef.all <- function(obj){

  Estimate  <- coef(obj)
  Std.Error <- se.coefficients(obj)
  z.value   <- Estimate/Std.Error
  p.value   <- 2*pnorm(abs(z.value), lower=F)

  return( data.frame(Estimate, Std.Error, z.value, p.value) )
}


formatOR <- function(f, digits=getOption("digits")){

  if( family(f)$family != "binomial" ) stop("This function only works
on GLM class objects with binomial family")

  est <- coefficients(summary(f))[ ,1]
  se  <- coefficients(summary(f))[ ,2]

  out <- data.frame( odds.ratio = exp( est ),
                     lcl        = exp( est - 1.96*se ),
                     ucl        = exp( est + 1.96*se ) )
  out <- cbind(coefficients(summary(f)), out)
  round(out, digits=digits)
}

resids.all <- function(glmfit){

 r   <- glmfit$residuals
 out <- rep(NA, length(glmfit$na.action) + length(r))
 valid.pts <- setdiff( 1:length(out), glmfit$na.action )
 out[ valid.pts ] <- r
 return(out)
}


#########################
## generic for metaGEM ##
#########################

summarize.probe2geneID <- function(keys){

  bad <- which(is.na(keys) | keys == "NA")
  x <- keys[ setdiff( 1:length(keys), bad ) ]
  x <- strsplit(x, split=",")

  out <- c(n.probes  = length(keys),                   # no. of probeIDs at start
           n.mapped  = length(x),                      # no. of un-mappable probes
           n.one2one = sum( sapply(x, length)==1 ),    # probes with one-to-one mapping
           n.records = length( unlist(x) ),            # no. of records after probe expand
           n.geneIDs = length( unique( unlist(x) ) ) ) # no. of geneIDs after probe summary
  return(out)
}


check.metaGEM.study <- function(study){

    stopifnot( all( c("expr", "class", "keys") %in% names(study) ) )

    if( !all(levels(as.factor(study$class)) == c("0", "1")) )
        stop("study$class must be coded as 0 or 1")

    if( !is.character(study$keys) )
        stop("The keys must be stored as a character vector")

}
plot.missing <- function(mat, sort=FALSE, main="Location of missing
values", ...){

  image2 <- function(m, ...) image( t(m)[ , nrow(m):1 ], ... )

  mat <- 1*is.na(mat)
  if(sort) mat <- mat[ order(rowSums(mat)), order(colSums(mat)) ]

  image2( mat, col=c(0,1), xaxt="n", yaxt="n", main=main, ... )
  box(); grid(col=4)

  if(!sort){ # does not make sense to show the grid position for sorted version
    ticks <- c(0,0.2,0.4,0.6,0.8,1.0)
    axis( side=1, at=ticks,      labels=round(quantile(1:ncol(mat), ticks)) )
    axis( side=2, at=rev(ticks), labels=round(quantile(1:nrow(mat), ticks)) )
  }
}

plotDensity <- function(...){
 inputs <- list(...)
 dens   <- lapply( inputs, density, na.rm=T )

 xranges <- sapply( dens, function(mylist) range(mylist$x) )
 yranges <- sapply( dens, function(mylist) range(mylist$y) )

 plot( dens[[1]], xlim=range(xranges), ylim=range(yranges), main="",
xlab="", type="n" )
 for(i in 1:length(dens)) lines( dens[[i]], col=i )
 grid()
 legend("topright", col=1:length(dens), lty=1, legend=names(inputs), bg="white")
}


targetLines <- function(x, y, approach="bottomleft", col=2, lty=1,
lwd=1, cex=0.85, digits=getOption("digits")){

   if (names(dev.cur()) == "null device")
          stop("No active plot to print to")

   if( !(approach %in% c("bottomleft", "bottomright", "topleft", "topright")) )
          stop("Check the 'approach' argument")

   xmin <- par("usr")[1];  xmax <- par("usr")[2]
   ymin <- par("usr")[3];  ymax <- par("usr")[4]
   par(las=1)

   bottom <- (length( grep("bottom", approach) ) > 0)
   segments( x, ifelse(bottom, ymin, ymax), x, y, col=col, lty=lty, lwd=lwd )
   mtext( round(x, digits=digits), side=ifelse(bottom, 1, 3), at=x,
cex=cex, col=col )

   left <- (length( grep("left$", approach) ) > 0)
   segments( ifelse(left, xmin, xmax), y, x, y, col=col, lty=lty, lwd=lwd )
   mtext( round(y, digits=digits), side=ifelse(left, 2, 4), at=y,
cex=cex, col=col )
}
strata.glm <- function( df, fmla, keep=1:nrow(df), strata="country",
plot=FALSE, family="gaussian" ){

  df <- df[ keep, c(all.vars(fmla), strata) ]
  df <- df[ which(complete.cases(df)), ]
  df[ , strata] <- factor( df[ , strata] )
  cat("Number of individuals used in the analysis:", nrow(df), "\n")

  ##########################################
  ### Finding the strata-level estimates ###
  ##########################################

  dfs  <- split( df, df[ , strata] )
  fits <- lapply( dfs, function(m) glm( fmla, data=m, family=family ) )

  #####################################
  ### Combining the strata-level estimates ###
  #####################################

  e  <- sapply( fits, coefficients )
  se <- sapply( fits, se.coefficients )

  p <- nrow(e)
  meta.e <- meta.se <- numeric( p )
  names(meta.e) <- names(meta.se) <- rownames( e )
  if(plot) par(mfrow=c( ceiling( (p-1) / 3 ), 3 ))

  for(i in 1:nrow(e)){
      m <- meta.summaries( cleanNA(e[i, ]), cleanNA(se[i, ]),
                          method="random", names=names( cleanNA(e[i, ]) ) )
      meta.e[i]  <- m$summary
      meta.se[i] <- m$se.summary
    if(plot){ if(i != 1) plot(m, main=rownames(e)[i] ) }
  }

  return( list(fits=fits, e=e, se=se, meta.e=meta.e, meta.se=meta.se) )
}


show.pooled <- function(obj, digits=2, conf.level=0.95){

    ci.value <- -qnorm((1 - conf.level)/2)

    e   <- cbind( obj$e,  pooled=obj$meta.e  )
    se  <- cbind( obj$se, pooled=obj$meta.se )
    lcl <- round( e - ci.value*se, digits )
    ucl <- round( e + ci.value*se, digits )

    p   <- 2*pnorm( abs(obj$meta.e/obj$meta.se), lower=FALSE )

    out <- matrix( paste( round(e, digits), " (", lcl, ",", ucl, ")",
sep="" ), nc=ncol(e) )
    dimnames(out) <- dimnames(e)
    out <- data.frame(out, pvalue.pooled=p)
    return(out)
}


plot.pooled <- function(fit, what, method="random", logscale=FALSE, ...){

 nm <- paste( colnames(fit$e), " (n=",
             sapply( fit$fits, function(f) nrow(f$data)) , ")", sep="" )

 m <- meta.summaries( fit$e[ what, ], fit$se[ what, ], names=nm,
method=method, logscale=logscale)

 print( summary(m) )
 plot( m, ... )
}

strata.glm.residuals <- function (df, fmla, keep=1:nrow(df), strata="country", plot=FALSE){

   df <- df[keep, c(all.vars(fmla), strata)]
   df <- df[which(complete.cases(df)), ]
   df[, strata] <- factor(df[, strata])
   cat("Number of individuals used in the analysis:", nrow(df), "\n")

   dfs  <- split(df, df[, strata])
   resids <- lapply(dfs, function(m) glm(fmla, data = m)$residuals)
   names(resids) <- NULL
   resids <- unlist(resids)[ rownames(df) ]
   return(resids)
}
