R/mb.1D.R

In timecourse: Statistical Analysis for Developmental Microarray Time Course Data

## (c) 2004-2005 Yu Chuan Tai
## University of California, Berkeley
## Created: August 2004     
## Last Modified: July 16, 2005

mb.1D <- function(object, k, n, nu=NULL, Lambda1=NULL, eta=NULL, k.grp=NULL, n.grp=NULL, 
r=FALSE, vec=FALSE, d=NULL, prop=0.01, T2.only=TRUE)
{

    res <- as.list(NULL)
    M <- NULL
    geneNames <- desc <- NULL
    if (is(object, "MAList") || is(object, "list")) 
          M <- object$M
   
    if (is(object, "marrayNorm"))
    {
          M <- object@maM
          geneNames <- rownames(M)
          desc <- object@maGnames@maLabels
    }
            
    if (is(object, "ExpressionSet")) 
    {
          M <- exprs(object)
          geneNames <- rownames(M)
    }
            
    if (is.null(M)) 
        M <- as.matrix(object)

   G <- nrow(M)
   if(missing(k)) stop("The number of time points is missing.")
   if(missing(n)) stop("The sample sizes are missing.")
   max.n <- max(n)
   if(!is.null(n) & length(n) !=G) stop("The sample sizes are incorrect.")
     
   if(is.null(k.grp)) 
   {
         k.grp <- rep(1:k, max.n)
         cat("Time group assignments are set to default.","\n")
   }
   if(length(unique(k.grp)) != k) 
      stop("The number of time points or the time group assignments is incorrect.")
      
   if(is.null(n.grp)) 
   {
        n.grp <- rep(c(1:max.n), rep(k, max.n))
        cat("Replicate group assignments are set to default.", "\n")
   }
   
   
   mydata <- M
   indx <- order(n.grp, k.grp)
   M <- M[,indx]
  
   mis <- apply(!apply(M, 1, is.na), 2, sum)
   mis <- sum((mis/k-floor(mis/k)) !=0)
   if(mis>0) stop(mis, " genes may have missing values in some of the replicates.")

 
   if((max.n*k) !=ncol(M)) stop("The number of time points or sample sizes are incorrect.")
   if(length(unique(n.grp)) != max.n) stop("Sample sizes or replicate group assignments is incorrect.")   
   c.grp <- rep(1, k*max.n)
   
       S1 <- apply(M, 1, matrix.cov, k, c.grp=c.grp)
       if(is.null(nu)) diagS1 <- apply(S1, 2, function(x) diag(matrix(x,ncol=k-1)))
       S1.avg <- matrix(apply(S1, 1, mean, na.rm=TRUE),ncol=k-1)
     
   

   if(is.null(nu)||is.null(Lambda1))
   {
     nu.lim <- k-1+6
   

    if(!is.null(nu))
    {
     nu0 <- nu
     nu <- max(nu0, nu.lim)
     if(is.infinite(nu) & is.null(Lambda1))  {Lambda1 <- S1.avg} 
     if(is.finite(nu)& is.null(Lambda1))    {Lambda1 <- (nu-k)*S1.avg/nu }
     nu<- nu0
   }

   if(is.null(nu))
   {
     nu0 <- mean(sapply(1:(k-1), function(x) squeezeVar(diagS1[x,], n-1)$df.prior))
     nu <- max(nu0, nu.lim)
     if(is.infinite(nu)& is.null(Lambda1))  {Lambda1 <- S1.avg} 
     if(is.finite(nu)& is.null(Lambda1))    {Lambda1 <- (nu-k)*S1.avg/nu }
     nu<- nu0
   }
   }
   
   max.n <- max(n)
   T1 <- ot.helmert(k)[2:k,]
   xbar <- apply(M, 1, function(x) apply(matrix(x,  byrow=FALSE, ncol=max.n),1,mean,na.rm=TRUE))

   if(r)
   {    
        e <- sapply(1:G, function(x) matrix(M[x,], byrow=FALSE, ncol=max.n)-xbar[,x])
        tune <- sapply(1:G, function(x) d * median(abs(e[, x]),na.rm = TRUE)/0.6745)
        indx <- sapply(1:G, function(x) abs(e[, x]) > tune[x])
        na.indx <- sapply(1:G, function(x) c(1:(k * max.n))[!is.na(indx[,x])])
        wt <- matrix(rep(1, G * max.n * k), ncol = max.n * k)
        wt <- t(sapply(1:G, function(x) {wt[x,][-na.indx[,x]] <- 0 
              wt[x,]}))        
        wt <- t(sapply(1:G, function(x) {wt[x, ][na.indx[,x]][indx[,x][na.indx[,x]]] <- 
        tune[x]/abs(e[,x][na.indx[,x]][indx[,x][na.indx[,x]]])
        wt[x,]}))
        totalw <- sapply(1:G, function(x) apply(matrix(wt[x,], byrow=FALSE, ncol=max.n), 1, sum, 
        na.rm=TRUE))
        w <- sapply(1:G, function(x) matrix(wt[x,], byrow=FALSE, ncol=max.n)/totalw[,x])
        xbar <- sapply(1:G, function(x) 
        apply(matrix(w[, x] * M[x, ], byrow = FALSE, ncol = max.n),1, sum,na.rm=TRUE))

       
            
    }

   tol <- .Machine$double.eps
   if(is.finite(nu) & nu > tol) 
   {
     modS1 <- sapply(1:G, function(x) ((n[x]-1)*matrix(S1[,x],ncol=k-1)+nu*Lambda1)/(n[x]-1+nu))
     if(is.null(eta) & (T2.only==FALSE))
     {
       sqrt.modS1 <- apply(modS1, 2, function(x)
       svd(matrix(x,ncol=k-1))$u%*%diag(sqrt(svd(matrix(x,ncol=k-1))$d))%*%t(svd(matrix(x,ncol=k-1))$v))
       modt1 <- sapply(1:G, function(x) n[x]^(1/2)*solve(matrix(sqrt.modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])
       HotellingT2 <- apply(modt1, 2, function(x) t(x)%*%x)
     }
     if(T2.only)
         HotellingT2 <- sapply(1:G, function(x) 
                n[x]*t(T1%*%xbar[,x])%*%solve(matrix(modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])

      if(!is.null(eta) & (T2.only==FALSE))
         HotellingT2 <- sapply(1:G, function(x)
                n[x]*t(T1%*%xbar[,x])%*%solve(matrix(modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])

   }

   if(is.infinite(nu))
   {
     modS1 <- sapply(1:G, function(x) Lambda1)
     if(is.null(eta) & (T2.only==FALSE))
     {
     sqrt.modS1 <- apply(modS1, 2, function(x)
     svd(matrix(x,ncol=k-1))$u%*%diag(sqrt(svd(matrix(x,ncol=k-1))$d))%*%t(svd(matrix(x,ncol=k-1))$v))
     modt1 <- sapply(1:G, function(x) n[x]^(1/2)*solve(matrix(sqrt.modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])
     HotellingT2 <- apply(modt1, 2, function(x) t(x)%*%x)
     }
     if(T2.only)
     HotellingT2 <- sapply(1:G, function(x)
         n[x]*t(T1%*%xbar[,x])%*%solve(matrix(modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])

     if(!is.null(eta) & (T2.only==FALSE))
         HotellingT2 <- sapply(1:G, function(x)
                n[x]*t(T1%*%xbar[,x])%*%solve(matrix(modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])
   }
   if(nu < tol & nu >=0)
   {
     modS1 <- S1
     if(is.null(eta) & (T2.only==FALSE))
     {
     sqrt.modS1 <- apply(modS1, 2, function(x)
     svd(matrix(x,ncol=k-1))$u%*%diag(sqrt(svd(matrix(x,ncol=k-1))$d))%*%t(svd(matrix(x,ncol=k-1))$v))
     modt1 <- sapply(1:G, function(x) n[x]^(1/2)*ginv(matrix(sqrt.modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])
     HotellingT2 <- apply(modt1, 2, function(x) t(x)%*%x)
     }
     if(T2.only)
     HotellingT2 <- sapply(1:G, function(x)
     n[x]*t(T1%*%xbar[,x])%*%ginv(matrix(modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])

     if(!is.null(eta) & (T2.only==FALSE))
         HotellingT2 <- sapply(1:G, function(x)
                n[x]*t(T1%*%xbar[,x])%*%solve(matrix(modS1[,x],ncol=k-1))%*%T1%*%xbar[,x])
   }

     
   # HotellingT2 does not require eta
   if(is.null(eta) & (T2.only==FALSE)) eta <- mean(sapply(1:(k-1), function(x) tmixture.matrix(modt1[x,], 
   sqrt(1/n), n+nu-1, prop, c(0.1, 4)^2/Lambda1[x,x]))^(-1))

   res$M <- mydata
   res$prop <- prop
   res$nu <- nu
   res$Lambda1 <- Lambda1 
   res$eta <- eta

 
   if(nu < 0) stop("The degrees of moderation <0 !")
   
   tol <- .Machine$double.eps
   if(is.finite(nu) & nu > tol)
   {
      MB <- NULL
      if(T2.only==FALSE)
      {
      MB1 <- log(prop,10)-log(1-prop,10)
      MB2 <- ((k-1)/2)*log((eta/(n+eta)),10)
      MB3 <- ((n+nu)/2)*log(((n-1+nu+HotellingT2)/(n-1+nu+(eta/(n+eta))*HotellingT2)), 10)
      MB <- MB1+MB2+MB3
      }
      unique.n <- unique(n)
      res$percent <- matrix(c(unique.n, round(100*nu/(nu+unique.n-1))), ncol=2)
      res$size <- n
      res$con.group <- rep(1, k*max.n)
      res$rep.group <- n.grp
      res$time.group <- k.grp
      if(r) res$weights <- wt
      if(vec) res$modt1 <- modt1
      res$HotellingT2 <- HotellingT2 
      res$MB <- MB
   }

   if(is.infinite(nu))
   {
      MB <- NULL
      if(T2.only==FALSE)
      {
      MB1 <- log(prop,10)-log(1-prop,10)
      MB2 <- ((k-1)/2)*log((eta/(n+eta)),10)
      MB3 <- 0.5*(n/(n+eta))*HotellingT2-log(10)
      MB <- MB1+MB2+MB3
      }

      unique.n <- unique(n)
      res$percent <- matrix(c(unique.n, 100), ncol=2)
      res$size <- n
      res$con.group <- rep(1, k*max.n)
      res$rep.group <- n.grp
      res$time.group <- k.grp
      if(r) res$weights <- wt
      if(vec) res$modt1 <- modt1
      res$HotellingT2 <- HotellingT2 
      res$MB <- MB
   }

   if(nu < tol & nu >=0)
   {
      MB <- NULL
      if(T2.only==FALSE)
      {
      MB1 <- log(prop,10)-log(1-prop,10)
      MB2 <- ((k-1)/2)*log((eta/(n+eta)),10)
      MB3 <- (n/2)*log(((n-1+HotellingT2)/(n-1+(eta/(n+eta))*HotellingT2)), 10)
      MB <- MB1+MB2+MB3
      }

      unique.n <- unique(n)
      res$percent <- matrix(c(unique.n, 0), ncol=2)
      res$size <- n
      res$con.group <- rep(1, k*max.n)
      res$rep.group <- n.grp
      res$time.group <- k.grp
      if(r) res$weights <- wt
      if(vec) res$modt1 <- modt1
      res$HotellingT2 <- HotellingT2 
      res$MB <- MB
   }
    
   res$pos.HotellingT2 <- order(HotellingT2, decreasing = TRUE)
   if(!is.null(MB)) res$pos.MB <- order(MB, decreasing = TRUE)
   res$geneNames <- geneNames
   res$descriptions <- desc
   new("MArrayTC", res)
   
}