R/spectral.r

In algstat: Algebraic statistics in R

#' Analyze a Rank Dataset
#'
#' \code{spectral} analyzes a rank dataset for order interactions; see examples for details.  
#' 
#' @param data a vector in the frequency format (see examples)
#' @param n the number of objects to select from
#' @param k the number of objects selected
#' @param levels the names of the outcomes, in the proper order
#' @param iter iterations in metropolis
#' @param burn burn in
#' @param thin thinning
#' @return a list containing named elements
#' \itemize{
#'   \item \code{effects}: the pure ith order effects as a data frame computed by projecting the data onto the isotypic subspaces.  the length of each is the same as the data, choose(n, k).
#'   \item \code{effectsNorms}: the l2 norms of each effect. 
#'   \item \code{statsMatrices}: the lower order statistics calculating matrices, made with Tmaker.
#'   \item \code{moves}: the markov moves for moving around each V, computed with markov on the statsMatrices.  only the positive moves are given.
#'   \item \code{samps}: the samples from each space conditioned on each level of statistics.  this is the output of metropolis.
#'   \item \code{obs}: a list of the observed data and its lower level summaries.
#'   \item \code{exp}: a list of the expected number of samples at each level given the summary statistics at the previous (lower) level.  these are computed from the samples from metropolis by (1) summarizing them with Tmaker and then (2) averaging the samples.  the expected V2 samples (for example), are determined by taking the samples with the same V1 statistics (in V3, say), summarizing them to V2 with Tmaker, and then averaging every cell across the samples.
#'   \item \code{fullExp}: this is the result of taking each of the samples with the same lower-order statistics and averaging them.  see exp, which is a reduction of fullExp.
#'   \item \code{residuals}: obs - exp
#'   \item \code{isotypicBases}: a list of the basis vectors of each isotypic subspace; computed as the eigenvalues of the result of Amaker, and grouped by eigenvalue.
#'   \item \code{sampsEffects}: the effects determined by each of the samples, projected onto the isotypic subspaces.
#'   \item \code{sampsEffectsNorms }: the norms of the effects of the samples.
#'   \item \code{sampsEffectsNormSummary }: a summary of the norms of the effects of the samples.
#'   \item \code{showStages}: a function that prints out the observed, expected, and residuals of sequentially conditioning on the sample size, first order statistics, second order statistics, and so on.
#'   \item \code{showFit}: a function that prints out a summary of the fit of the model.
#'   \item \code{decompose}: a function that takes a vector of the same length of the table given and summarizes it to its lower level statistics.
#'   \item \code{sampsDecomposed}: every sample decomposed.
#'   \item \code{statistic}: the pearson's chi-squared (X2), likelihood ratio (G2), Freeman-Tukey (FT), Cressie-Read (CR), and Neyman modified chi-squared (NM) statistics computed using the observed data (obs) and expected data (exp) at each level.
#'   \item \code{sampsStats}: the statistics computed on each of the samples.
#'   \item \code{p.value}: the exact p-values of individual tests, accurate to Monte-Carlo error.  these are computed as the proportion of samples with statistics equal to or larger than the oberved statistic.
#'   \item \code{p.value.se}: the standard errors of the p-values computed using the standard asymptotic formula of sqrt(p(1-p)/n).  a measure of the Monte-Carlo error.
#' }	

#' @export spectral
#' @examples
#'
#' \dontrun{
#'
#'
#'
#' ## voting statistics at different levels
#' ############################################################
#'
#' # load the cookies dataset:
#' data(cookie)
#' cookie$freq
#' cookie$cookies
#'
#'
#' # performing the spectral analysis
#' (out <- spectral(cookie$freq, 6, 3, cookie$cookies))
#'  
#'
#' out$obs # the original observations, and the summary statistics
#'
#' out$exp # each level is conditional on the previous level's statistics
#'         # (e.g. what you would expect for 1st order effects given sample size)
#'         # these are generated using 10k markov bases based mcmc samples
#'
#' out$p.value # these are approximate exact test p-values using various 
#'             # popular test statistics.  the approximations are good to
#'             # monte carlo error
#'
#' out$p.value.se # these are the standard errors using the sqrt(p*(1-p)/n)
#'                # asymptotic formula, known to have poor performance
#'                # for small/large p; see package binom for better
#'
#' out$statistic # the individual statistics are also available
#'               # the values are not comprable across Vi levels (the rows)
#'               # as they have asymptotic chi-squared distributions with 
#'               # different degrees of freedom
#'
#' out$fullExp # you can also get the expected number of samples at each scale
#'             # for tables with the same ith order statistics, i = 0, ..., k-1
#'
#'
#' # these can be seen to (re)construct an expected picture of the 
#' # complete data given each successive collection of statistics
#' cbind(
#'   obs = cookie$freq, 
#'   as.data.frame(lapply(out$fullExp, function(x) round(x[[4]],1)))
#' )[c(2:4,1)]
#' # notice that the reconstruction given only the first order statistics
#' # (the number of individual cookies selected) is quite good
#'
#' 
#' 
#' # instead of using the reconstructions from the exp coming from 
#' # the samples, you could reconstruct the summaries of the observed
#' # data using bump; it's not quite as good :
#' V0 <- bump(cookie$freq, 6, 3, 3, 0)
#' V1 <- bump(cookie$freq, 6, 3, 3, 1)
#' V2 <- bump(cookie$freq, 6, 3, 3, 2)
#' 
#' cbind(
#'   obs = cookie$freq, 
#'   round(data.frame(
#'     V0 = bump(V0, 6, 3, 0, 3),
#'     V1 = bump(V1, 6, 3, 1, 3),    
#'     V2 = bump(V2, 6, 3, 2, 3)
#'   ), 2)
#' )[c(2:4,1)]
#' 
#' 
#'
#'
#' # you can see the model step-by-step with showStages() :
#' out$showStages()
#' # notice (1) the significant reduction in the residuals after conditioning
#' # on the first order statistics and also (2) the powdery noise after 
#' # conditioning on the second order statistics.
#' # the p-values reflect the same: 
#' #   * the residuals from conditioning on the sample size show the first
#' #     order effects are strongly significant (in out$p.value V1 = 0)
#' #   * the residuals from conditioning on the first order effects suggest
#' #     the second order effects might be significant (V2 ~ .04-.13ish)
#' #   * the residuals from conditioning on the second order effects indicate
#' #     the third order effects are entirely insignificant (V3 > .2)
#' 
#' 
#' # the isotypic subpaces can be used to determine the pure order effects :
#'
#' out$isotypicBases # bases of the isotypic subspaces (here 4)
#'
#' out$effects # pure ith order effects; cookie$freq projected onto the bases
#'             # these are their effects at the data level, so they all have
#'             # the same length as the original dataset: choose(n, k)
#'
#' zapsmall(rowSums(out$effects)) # the effects sum to the data
#'
#'
#' # if the Vk effects are 0, then the conclusion is that Vk is perfectly
#' # predicted with the (k-1)st level statistics.  this may lead to the 
#' # conclusion that the l2 norms (say) of the effects might be used to 
#' # gauge the relative strength of effects :
#' out$effectsNorms # = apply(out$effects, 2, lpnorm)
#' 
#' 
#' # the natural (not full-dimensional) residuals can be seen with the summary
#' out
#' # or with
#' out$residuals
#' # these are the residuals (obs ith level stats) - (exp ith level stats)
#' # given the (i-1)st statistics
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' # bump is a useful function :
#' out$obs
#' bump(cookie$freq, 6, 3, 3, 0) # the 0 level is the number of voters, not votes
#' bump(cookie$freq, 6, 3, 3, 1) 
#' bump(cookie$freq, 6, 3, 3, 2)
#' bump(cookie$freq, 6, 3, 3, 3)
#' 
#' V1 <- out$obs$V1 # = bump(cookie$freq, 6, 3, 3, 1)
#' bump(V1, 6, 3, 1, 0)
#' bump(V1, 6, 3, 1, 1)
#' bump(V1, 6, 3, 1, 2) # cbind(bump(V1, 6, 3, 1, 2), out$exp$V2)
#' bump(V1, 6, 3, 1, 3) # cbind(bump(V1, 6, 3, 1, 3), out$fullExp$V1[[4]])
#' # the differences here are between an observation and an expectation
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' out$obs$V1 - out$exp$V1
#' out$residuals$V1
#' out$decompose(out$effects$V1)$V1
#' 
#' out$obs$V2 - out$exp$V2
#' out$residuals$V2
#' 
#'   out$decompose(out$effects$V0)$V2 + 
#'   out$decompose(out$effects$V1)$V2 +
#'   out$decompose(out$effects$V2)$V2 -
#'   out$exp$V2
#'   
#'
#'
#'
#' 
#' # this is how to reconstruct the observation given the effects
#' # the cols of out$effects are the Vk order effects reconstructed
#' # from the lower level effects
#' out$obs$V0
#' zapsmall(
#'   out$decompose(out$effects$V0)$V0
#' )
#' 
#' out$obs$V1
#' zapsmall(
#'   out$decompose(out$effects$V0)$V1 + 
#'   out$decompose(out$effects$V1)$V1
#' )
#' 
#' out$obs$V2
#' zapsmall(
#'   out$decompose(out$effects$V0)$V2 + 
#'   out$decompose(out$effects$V1)$V2 +
#'   out$decompose(out$effects$V2)$V2
#' )
#' 
#' out$obs$V3
#' zapsmall(
#'   out$decompose(out$effects$V0)$V3 + 
#'   out$decompose(out$effects$V1)$V3 +
#'   out$decompose(out$effects$V2)$V3 +
#'   out$decompose(out$effects$V3)$V3
#' )
#' zapsmall(rowSums(out$effects))
#' 
#' all(cookie$freq == zapsmall(rowSums(out$effects)))
#' 
#' 
#' 
#' out$effects$V0
#' out$effects$V0 + out$effects$V1
#' out$effects$V0 + out$effects$V2
#' out$effects$V0 + out$effects$V3
#' 
#' 
#' 
#' str(out$sampsDecomposed)
#' as.data.frame(lapply(out$sampsDecomposed, function(l) rowMeans(l$V3)))
#' 
#' eff0 <- rowMeans(out$sampsDecomposed$V0$V3)
#' cbind(eff0, out$effects$V0)
#' 
#' eff1 <- rowMeans(out$sampsDecomposed$V1$V3 - eff0)
#' cbind(eff1, out$effects$V1)
#' 
#' eff2 <- rowMeans(out$sampsDecomposed$V2$V3 - eff0 - eff1)
#' cbind(eff2, out$effects$V2)
#' 
#' sum(eff0)
#' sum(eff1)
#' sum(eff2)
#' 
#' 
#' 
#' 
#' 
#'
#'
#' str(out$sampsEffectsNorms)
#' 
#' data <- out$sampsEffectsNorms$V0$V3
#' plot(density(data))
#' curve(dnorm(x, mean(data), sd(data)), col = "red", add = TRUE)
#' 
#' data <- out$sampsEffectsNorms$V0$V2
#' plot(density(data))
#' curve(dnorm(x, mean(data), sd(data)), col = "red", add = TRUE)
#' 
#' data <- out$sampsEffectsNorms$V0$V1
#' plot(density(data))
#' curve(dnorm(x, mean(data), sd(data)), col = "red", add = TRUE)
#' 
#' 
#' data <- out$sampsEffectsNorms$V1$V3
#' plot(density(data))
#' curve(dnorm(x, mean(data), sd(data)), col = "red", add = TRUE)
#' 
#' data <- out$sampsEffectsNorms$V1$V2
#' plot(density(data))
#' curve(dnorm(x, mean(data), sd(data)), col = "red", add = TRUE)
#' 
#' 
#' data <- out$sampsEffectsNorms$V2$V3
#' plot(density(data))
#' curve(dnorm(x, mean(data), sd(data)), col = "red", add = TRUE)
#' 
#' 
#' 
#' 
#' 
#' 
#' ## how to convert data into the right format
#' ############################################################
#' # this essentially just uses some clever indexing tricks
#' # to reorder the data in the way you want
#' 
#' data <- cookie$raw       # an example raw, unordered dataset
#' levels <- cookie$cookies # the order of the objects you want
#' levsNndcs <- 1:length(levels)
#' names(levsNndcs) <- levels
#' 
#' 
#' # arrange selections within rows (order of selection doesn't matter)
#' data <- t(apply(data, 1, function(x) x[order(levsNndcs[x])] ))
#'
#' 
#' # arrange rows (order of selectors doesn't matter)
#' for(k in ncol(data):1) data <- data[order(levsNndcs[data[,k]]),]
#' 
#' 
#' # check that you've done the right thing
#' all( data == cookie$sorted )
#' 
#' # the data frequency order should match that of subsets:
#' subsets(levels, 1)
#'
#' subsets(levels, 2)
#' sapply(subsets(levels, 2), paste, collapse = ", ")
#' 
#' subsets(levels, 3)
#' sapply(subsets(levels, 3), paste, collapse = ", ")
#' 
#' names(cookie$freq)
#' names(cookie$freq) == sapply(subsets(levels, 3), paste, collapse = ", ")
#'
#' 
#' 
#' 
#' 
#'
#' 
#' 
#' 
#' 
#'
#' 
#' 
#' 
#' 
#'
#' ## other examples
#' ############################################################
#' 
#' # rvotes provides uniform samples
#'
#' n <- 4
#' k <- 2
#' 
#' raw <- rvotes(250, n, k)
#' rawTogether <- apply(raw, 1, paste, collapse = " ")
#' levels <- sapply(subsets(n, k), paste, collapse = " ")
#' freq <- table( factor(rawTogether, levels = levels) )
#' (out <- spectral(freq, n, k))
#'
#' out$p.value
#' out$showStages()
#' 
#' out$obs
#' out$exp
#'
#'
#'
#'
#'
#' n <- 6
#' k <- 3
#' raw <- rvotes(250, n, k)
#' rawTogether <- apply(raw, 1, paste, collapse = " ")
#' levels <- sapply(subsets(n, k), paste, collapse = " ")
#' freq <- table( factor(rawTogether, levels = levels) )
#' (out <- spectral(freq, n, k))
#' 
#' 
#' n <- 7
#' k <- 3
#' raw <- rvotes(250, n, k)
#' rawTogether <- apply(raw, 1, paste, collapse = " ")
#' levels <- sapply(subsets(n, k), paste, collapse = " ")
#' freq <- table( factor(rawTogether, levels = levels) )
#' (out <- spectral(freq, n, k))
#'
#'
#'
#' n <- 8
#' k <- 3
#' raw <- rvotes(250, n, k)
#' rawTogether <- apply(raw, 1, paste, collapse = " ")
#' levels <- sapply(subsets(n, k), paste, collapse = " ")
#' freq <- table( factor(rawTogether, levels = levels) )
#' # out <- spectral(freq, n, k) # breaks
#'
#' 
#' 
#' 
#' }
spectral <- function(data, n, k, levels, iter = 1e4, burn = 1e3, thin = 10){
  
  if(missing(levels)) levels <- 1:n     
  if(!missing(levels) && length(levels) != n){
  	stop("if specified, levels must be of length n.", call. = FALSE)
  }
  
  
  # make list of moves
  message("Computing moves... ", appendLF = FALSE)
  listOfMoves <- as.list(rep(NA, k))
  names(listOfMoves) <- paste0("V", 0:(k-1))
  for(d in 0:(k-1)){
    listOfMoves[[d+1]] <- markov(Tmaker(n, k, d))
  }
  message("done.")



  # run metropolis-hastings
  samps <- as.list(rep(NA, k))
  names(samps) <- paste0("V", 0:(k-1))  
  for(d in 0:(k-1)){
    samps[[d+1]] <- metropolis(unname(data), listOfMoves[[d+1]],
      iter = iter, burn = burn, thin = thin
    )
  }
  
  
  

  # compute the A matrix
  A <- Amaker(n, k)
  
  # compute eigen stuff  
  e <- eigen(A)
  
  # clean eigen
  e$values <- zapsmall(e$values)
  e$vectors <- e$vectors
  
  # number of evals
  nEvals <- length(unique(e$values))  
  
  # split evecs based on eval
  listOfBasisVectors <- split(
    split(e$vectors, col(e$vectors)),
    e$values
  )
  
  # list from V0 to Vk
  listOfBasisVectors <- rev(listOfBasisVectors)
  
  # convert to matrices
  listOfBasisVectors <- lapply(listOfBasisVectors, function(x){
    if(is.list(x)) return(unname(Reduce(cbind, x)))
    t(t(x))
  })
  
  # assign names
  names(listOfBasisVectors) <- paste0("V", 0:k)
  
  
 
  
  # compute projections of data  
  qrs <- lapply(listOfBasisVectors, qr)  
  effects <- lapply(qrs, function(x) as.numeric(qr.fitted(x, data)) )
  effects <- as.data.frame(effects)
  row.names(effects) <- sapply(subsets(levels, k), paste, collapse = ", ")

  # norm projections of data
  effectsNorms <- sapply(effects, lpnorm)




  # compute projections of samps and their effectsNorms
  relative2Data <- sampsEffects <- sampsEffectsNorms <- as.list(rep(NA, k))
  names(sampsEffects) <- names(sampsEffectsNorms) <- paste0("V", 0:(k-1))
  ph <- as.list(rep(NA, k+1))  
  names(ph) <- paste0("V", 0:k)
  message("Projecting and norming... ", appendLF = FALSE)       

  for(d in 0:(k-1)){ # cycle through the same ith-order effects

  	sampsEffects[[d+1]] <- sampsEffectsNorms[[d+1]] <- relative2Data[[d+1]] <- ph

    for(dd in 0:k){ # cycle through dimensions to project onto    	

      # project      
      sampsEffects[[d+1]][[dd+1]] <- qr.fitted(qrs[[dd+1]], samps[[d+1]]$steps)        
      
      # norm
      sampsEffectsNorms[[d+1]][[dd+1]] <- apply(sampsEffects[[d+1]][[dd+1]], 2, lpnorm)      
      
      # compare to data
      relative2Data[[d+1]][[dd+1]] <- 
        mean( sampsEffectsNorms[[d+1]][[dd+1]] <= effectsNorms[dd+1] )
    
    }
    
  }
  message("done.")
  
  
  

  # summarize norms
  summarize <- function(v){
    if(length(v) == 1 && is.na(v)) return(NA)
    c(
        #"2.5%" = unname(quantile(v, .025)),
        "Mean" = mean(v),
       #"97.5%" = unname(quantile(v, .975)),
           "Sd" = sd(v)
    )
  }
  sampsEffectsNormsSummary <- lapply(sampsEffectsNorms, lapply, summarize)
  sampsEffectsNormsSummary <- lapply(sampsEffectsNormsSummary, as.data.frame)
  sampsEffectsNormsSummary <- lapply(sampsEffectsNormsSummary, function(df) df[,1:(k+1)])  
  relative2Data <- lapply(relative2Data, as.data.frame)
  for(i in 1:length(sampsEffectsNormsSummary)){
    sampsEffectsNormsSummary[[i]]   <- rbind(
      sampsEffectsNormsSummary[[i]], 
      unname(effectsNorms),      
      relative2Data[[i]]
    )
    row.names(sampsEffectsNormsSummary[[i]]) <- 
      c("Mean", "Sd", "data", "% <= data")
  }
  names(sampsEffectsNormsSummary) <- paste(names(sampsEffectsNormsSummary), "samples")




  # name listOfBasisVectors, make into matrices, zapsmall
  listOfBasisVectors <- lapply(listOfBasisVectors, function(x){
    if(is.list(x)) return(unname(Reduce(cbind, x)))
    t(t(x))
  })
  listOfBasisVectors <- lapply(listOfBasisVectors, zapsmall)
  

  
  
  # generate summary statistics
  observed <- as.list(rep(NA, k+1))
  for(d in 0:k){
  	observed[[d+1]] <- Tmaker(n, k, d) %*% unname(data)
    row.names(observed[[d+1]]) <- 
  	  sapply(subsets(levels, d), paste, collapse = ", ")
    colnames(observed[[d+1]]) <- "N"
  }
  names(observed) <- paste0("V", 0:k)
  
  

  # compute expected number of samples at each level
  # given the observed number at each level
  fullExpected <- as.list(rep(NA, k))
  names(fullExpected) <- paste0("V", 0:(k-1))
  Ts <- lapply(as.list(0:k), function(x) Tmaker(n, k, x)) # list of T's  
  for(d in 0:(k-1)){
  	fullExpected[[d+1]] <- lapply(Ts, function(Tmat){
  	  rowMeans(Tmat %*% samps[[d+1]]$steps)
  	})    
  }
  
  
  
  # summarize the above to the expected number of samples
  # at each level given the observed number of samples
  # at the previous level
  expected <- observed
  for(d in 1:k){
    expected[[d+1]] <- t(t(fullExpected[[d]][[d+1]]))
    row.names(expected[[d+1]]) <- 
      sapply(subsets(levels, d), paste, collapse = ", ")
    colnames(expected[[d+1]]) <- "N"      
  }
  
  
  
  
  # compute residuals
  residuals <- observed
  for(d in 1:k){
    residuals[[d+1]] <- observed[[d+1]] - expected[[d+1]]
  }  
  
  
  
  # showing stages
  showStages <- function(){
  	count <- 0
  	postScript <- c(
  	  "the sample size :\n",
  	  "first order statistics :\n",
  	  "second order statistics :\n",
  	  "third order statistics :\n",
  	  "fourth order statistics :\n",
  	  "fifth order statistics :\n",
  	  "sixth order statistics :\n"
  	)
    lapply(fullExpected, function(l){
      resMat <- obsMat <- expMat <- 
        matrix(numeric(0), nrow = choose(n,k), ncol = k+1)  
      obsMat[2,1] <- expMat[2,1] <- observed[[1]]          
      for(d in 0:k){
        expMat[1:choose(n, d), d+1] <- l[[d+1]]
        obsMat[1:choose(n, d), d+1] <- observed[[d+1]][,1]
        resMat[1:choose(n, d), d+1] <- observed[[d+1]][,1] - l[[d+1]]
      }
      mat <- cbind(obsMat, expMat, resMat)
      mat <- round(mat, 1)
      mat2 <- apply(mat, 2, format)
      mat2[is.na(mat)] <- ""
      mat2[1,0*(k+1)+1] <- "Obs: "
      mat2[1,1*(k+1)+1] <- "Exp: "    
      mat2[1,2*(k+1)+1] <- "Resid:"      
      mat2[3,0*(k+1)+1] <- paste0("(x", k, ")")
      mat2[3,1*(k+1)+1] <- paste0("(x", k, ")")      
      mat2 <- apply(mat2, 2, format)
      cat(paste("Conditioning on", postScript[count+1]))
      cat(apply(mat2, 1, paste, collapse = "  "), sep = "\n")      
      cat("\n")
      count <<- count + 1
    })  
    invisible()
  } 
  
  
  
  
  # showing expected
  showFit <- function(){
    resMat <- obsMat <- expMat <- 
      matrix(numeric(0), nrow = choose(n,k), ncol = k+1)  
    obsMat[2,1] <- expMat[2,1] <- observed[[1]]    
    for(d in 1:k){
      expMat[1:choose(n, d), d+1]  <- fullExpected[[d]][[d+1]]
      obsMat[1:choose(n, d), d+1] <- observed[[d+1]][,1]
      resMat[1:choose(n, d), d+1] <- residuals[[d+1]][,1]
    }
    mat <- cbind(obsMat, expMat, resMat)
    mat <- round(mat, 1)
    mat2 <- apply(mat, 2, format)
    mat2[is.na(mat)] <- ""
    mat2[1,0*(k+1)+1] <- "Obs: "
    mat2[1,1*(k+1)+1] <- "Exp: "    
    mat2[1,2*(k+1)+1] <- "Resid:"
    mat2[3,0*(k+1)+1] <- paste0("(x", k, ")")
    mat2[3,1*(k+1)+1] <- paste0("(x", k, ")")
    mat2 <- apply(mat2, 2, format)    
    cat(apply(mat2, 1, paste, collapse = "  "), sep = "\n")
    cat("\n")
    invisible()  
  } 
  
  
  
  
  # decompose an observation into the various order effects
  decompose <- function(v){
    out <- as.list(rep(NA, k+1))
    names(out) <- paste0("V", 0:k)
    for(d in 0:k){
      out[[d+1]] <- Ts[[d+1]] %*% unname(v)
      row.names(out[[d+1]]) <- 
        sapply(subsets(levels, d), paste, collapse = ", ")
      colnames(out[[d+1]]) <- "N"
    }
    out
  }   
    



  # decompose samps
  sampsDecomposed <- lapply(samps, function(x){
    out <- as.list(rep(NA, k+1))
    names(out) <- paste0("V", 0:k)  	
    for(d in 0:k) out[[d+1]] <- Ts[[d+1]] %*% x$steps
    out <- lapply(out, function(mat){
      matrix(as.integer(mat), nrow = nrow(mat))
    })
    out
  })
  #for(d in 1:k) sampsDecomposed[[d]] <- sampsDecomposed[[d]][[d+1]]
  #names(sampsDecomposed) <- paste0("V", 1:k)
  #lapply(sampsDecomposed, rowMeans) # = expected
  
  


  # compute statistics for each combo of observed/expected  
  X2 <- vector(length = k)
  names(X2) <- paste0("V", 1:k)
  PR <- G2 <- CR <- FT <- NM <- X2
  for(d in 1:k){
    PR[d] <- computeUProbsCpp(observed[[d+1]])  	
    X2[d] <- computeX2sCpp(observed[[d+1]], expected[[d+1]])
    G2[d] <- computeG2sCpp(observed[[d+1]], expected[[d+1]])    
    FT[d] <- computeCRsCpp(observed[[d+1]], expected[[d+1]], -.5)
    CR[d] <- computeCRsCpp(observed[[d+1]], expected[[d+1]], 2/3)    
    NM[d] <- computeNMsCpp(observed[[d+1]], expected[[d+1]])        
  }
  statistic <- rbind(PR, X2, G2, FT, CR, NM)
  
  
  
  
  
  # compute statistics for each combo of samps/expected  
  X2s <- as.list(rep(NA, length = k))
  names(X2s) <- paste0("V", 1:k)
  PRs <- G2s <- CRs <- FTs <- NMs <- X2s
  for(d in 1:k){
    PRs[[d]] <- computeUProbsCpp(sampsDecomposed[[d]][[d+1]])  	
    X2s[[d]] <- computeX2sCpp(sampsDecomposed[[d]][[d+1]], expected[[d+1]])
    G2s[[d]] <- computeG2sCpp(sampsDecomposed[[d]][[d+1]], expected[[d+1]])    
    FTs[[d]] <- computeCRsCpp(sampsDecomposed[[d]][[d+1]], expected[[d+1]], -.5)
    CRs[[d]] <- computeCRsCpp(sampsDecomposed[[d]][[d+1]], expected[[d+1]], 2/3)    
    NMs[[d]] <- computeNMsCpp(sampsDecomposed[[d]][[d+1]], expected[[d+1]])        
  }
  sampsStats <- list(PRs = PRs, X2s = X2s, G2s = G2s, FTs = FTs, CRs = CRs, NMs = NMs)

  
  
  
  # compute p values and their standard errors
  p.value <- matrix(nrow = nrow(statistic), ncol = k)
  colnames(p.value) <- paste0("V", 1:k)
  row.names(p.value) <- c("PR", "X2", "G2" ,"FT", "CR", "NM")
  mid.p.value <- p.value
  for(i in 1:nrow(statistic)){
    for(d in 1:k){
      p.value[i,d] <- mean( sampsStats[[i]][[d]] >= statistic[i,d] )
      mid.p.value[i,d] <- 
        mean( sampsStats[[i]][[d]] == statistic[i,d] )/2 +
        mean( sampsStats[[i]][[d]] > statistic[i,d] )
    }
  }
  p.value[1,] <- 1 - p.value[1,] 
  p.value.se <- sqrt(p.value*(1-p.value)/iter)
  mid.p.value.se <- sqrt(mid.p.value*(1-mid.p.value)/iter)  
  
    



  # return
  out <- list(
    effects = zapsmall(effects), effectsNorms = effectsNorms, 
    statsMatrices = Ts, moves = listOfMoves,
    samps = samps, 
    obs = observed, exp = expected, fullExp = fullExpected,
    residuals = residuals,    
    isotypicBases = listOfBasisVectors,    
    sampsEffects = sampsEffects, sampsEffectsNorms = sampsEffectsNorms,
    sampsEffectsNormsSummary = sampsEffectsNormsSummary,
    showStages = showStages, showFit = showFit,
    decompose = decompose, sampsDecomposed = sampsDecomposed,
    statistic = statistic, sampsStats = sampsStats,
    p.value = p.value, p.value.se = p.value.se,
    mid.p.value = mid.p.value, mid.p.value.se = mid.p.value.se,
    call = match.call()    
  )
  
  class(out) <- "spectral"
  out
}