phoc <- function(dframe, dv, sid, indfacs)
{
  if (any(class(dframe) == "tbl")) {
    dframe = data.frame(dframe)
  }
  dv = dframe[, names(dframe) == dv]
  sid = dframe[, names(dframe) == sid]
  indfacs = dframe[, names(dframe) %in%  indfacs]
  indfacs = apply(as.matrix(indfacs), 1, paste, collapse = ":")
  
  u = unique(indfacs)
  n = length(u)
  mat = NULL
  for (j in 1:(n - 1)) {
    w = u[j]
    temp = indfacs == w
    wdata = dv[temp]
    wspeak = sid[temp]
    if (any(duplicated(wspeak)))
      stop(
        "Unique values per subject-factor combination required; aggregate over factors and try again"
      )
    
    x = u[(j + 1):n]
    for (k in x) {
      paired = FALSE
      temp = indfacs == k
      kdata = dv[temp]
      kspeak = sid[temp]
      if (any(duplicated(kspeak)))
        stop(
          "Unique values per subject-factor combination required; aggregate over factors and try again"
        )
      
      
      if (length(wspeak) == length(kspeak))
        # if the lengths of the speaker vectors are the same
      {
        # make sure they are all the same
        temp = all(wspeak == kspeak)
        if (temp)
          # if so
        {
          # make sure they're aligned by speaker
          m = match(wspeak, kspeak)
          kdata = kdata[m]
          # and set paired to T
          paired = T
        }
      }
      # do the t-test
      res = t.test(wdata, kdata, paired = paired)
      mat$res = rbind(mat$res, c(res[[1]], res[[2]], res[[3]]))
      mat$name = c(mat$name, paste(u[j], k, sep = "-"))
      mat$paired = c(mat$paired, paired)
    }
  }
  mat$bonf = nrow(mat$res)
  rownames(mat$res) = mat$name
  colnames(mat$res) = c("t", "df", "prob-adj")
  # Bonf. adjust
  mat$res[, 3] = mat$res[, 3] * mat$bonf
  mat$res[, 3][mat$res[, 3] > 1] = 1
  mat
}

###########################################################################


phsel <- function(phocresults, k = 1, digits = 3)
{
    # written by Jonathan Harrington (with slight adaptations to phsel.emmGrid() by Ulrich Reubold)
    
    # choose between whether the output is from
    # phoc() or from glht() or from pairs(emmeans())
    
    
    phsel.glht <-
        function(mmphoc, k = 1)
        {
            mx = names(mmphoc$test$tstat)
            omatx = cbind(mmphoc$test$tstat, mmphoc$test$pvalues)
            colnames(omatx) = c("z value", "Adjusted p values")
            mx.un = matrix(unlist(strsplit(mx, " - ")), ncol = 2, byrow = T)
            # Number of independent variables
            nx = length(unlist(strsplit(mx.un[1, 1], ".", fixed = TRUE)))
            leftx = matrix(unlist(strsplit(mx.un[, 1], ".", fixed = TRUE)), ncol =
                               nx, byrow = T)
            rightx = matrix(unlist(strsplit(mx.un[, 2], ".", fixed = TRUE)), ncol =
                                nx, byrow = T)
            
            # leave out one or more  of the columns
            leftx = as.matrix(leftx[, -k])
            rightx = as.matrix(rightx[, -k])
            matx = NULL
            for (j in 1:nrow(leftx)) {
                vecx = all(leftx[j, ] == rightx[j, ])
                matx = c(matx, vecx)
            }
            (omatx[matx, ])
        }
    
    phsel.anova <- function(tukeyoutput, k = 1)
    {
        # written by Jonathan Harrington
        # the data is modifief from K. Johnson (Pitt_Shoaf2.txt)
        # psh = read.table(file.path(pfad, "psh.txt"))
        # carry out ANOVA
        # psh.aov = aov(rt ~ Overlap * Position, data=psh)
        # summary(psh.aov)
        # Tukey-Test
        # psh.tk = TukeyHSD(psh.aov)
        # Here are the components of the Tukey-Test
        # The interaction term is Overlap:Position
        # and it is is position 3
        # names(psh.tk)
        # Select the results of the Tukey test keeping
        # the first factor, Position constant
        # phsel(psh.tk[[3]])
        # the same
        # phsel(psh.tk[[3]], 1)
        # Select the results of the Tukey test keeping
        # the second factor, Overlap constant
        # phsel(psh.tk[[3]], 2)
        
        m = rownames(tukeyoutput)
        m.un = matrix(unlist(strsplit(m, "-")), ncol = 2, byrow = T)
        # Number of independent variables
        n = length(unlist(strsplit(m.un[1, 1], ":", fixed = TRUE)))
        left = matrix(unlist(strsplit(m.un[, 1], ":", fixed = TRUE)), ncol =
                          n, byrow = T)
        right = matrix(unlist(strsplit(m.un[, 2], ":", fixed = TRUE)), ncol =
                           n, byrow = T)
        
        # leave out one or more  of the columns
        left = as.matrix(left[, -k])
        right = as.matrix(right[, -k])
        mat = NULL
        for (j in 1:nrow(left)) {
            vec = all(left[j, ] == right[j, ])
            mat = c(mat, vec)
        }
        
        data.frame(tukeyoutput[mat, ])
    }
    
    phsel.emmGrid <- function(df.ph,k) 
    {
        # df.ph is the output of 
        # e.g. df.ph = pairs(emmeans(get_model(df.step),  ~ A:B:C:D))
        
        # convert to data-frame
        df.ph = data.frame(df.ph)
        
        # mat for storing the intermediary results
        # rmat final results
        # selects pairs for which all other features are constant
        rmat = mat = NULL
        # the first column has the contrasts
        vec = as.character(df.ph[,1])
        
        # split on the "-" character
        N = length(vec)
        for(j in 1:N){
            vecminus = splitstring(vec[j], "-")
            left = splitstring(vecminus[1], ",")
            right = splitstring(vecminus[2], ",")
            mat$left = rbind(mat$left, left)
            mat$right = rbind(mat$right, right)
        }
        # remove final blank space from the last column of mat$left
        n = nchar(mat$left[,ncol(mat$left)])
        mat$left[,ncol(mat$left)] = substring(mat$left[,ncol(mat$left)], 1, n-1)
        # remove single blank space from 1st column of mat$right
        n = nchar(mat$right[,1])
        mat$right[,1] = substring(mat$right[,1], 2, n)
        
        
        if (k > ncol(mat$right)){
            options(error=NULL)
            stop("k > number of factors")
        }
        else
        {
            # this selects the contrasts for which n-1 features in left and right are equivalent         
            for(j in 1:nrow(mat$right)){
                if (paste(mat$left[j, -k], collapse=".")  ==  paste(mat$right[j, -k], collapse=".") )
                    rmat = rbind(rmat, df.ph[j,])
            }
            rmat
        }
    }  
    
    # also needs:
    splitstring <- function (str, char) 
    {
        if (str == "") 
            mat <- c(str)
        else {
            mat <- NULL
            ind <- 1
            cont <- TRUE
            while (TRUE) {
                ministr <- NULL
                length <- 0
                while (TRUE) {
                    ch <- substring(str, ind, ind)
                    if (ch == char) {
                        ind <- ind + 1
                        break
                    }
                    if (ch == "") {
                        break
                    }
                    ministr <- c(ministr, ch)
                    ind <- ind + 1
                    length <- length + 1
                }
                if (length > 0) 
                    mat <- c(mat, paste(ministr, collapse = ""))
                if (ch == "") 
                    break
            }
        }
        mat
    }
    
    ###########################################################
    
    # function begins here and chooses between ghlt and anova
    if (any(class(phocresults) == "glht")) {
        rphoc = phsel.glht(phocresults, k)
    } else if (any(class(phocresults) == "emmGrid")) {
        rphoc = phsel.emmGrid(phocresults, k)
    } else{
        rphoc = phsel.anova(phocresults, k)
    }
    
    
    # print result after rounding the numbers, adding stars and name of adjustment
    rphoc$sig = ""
    rphoc$sig[rphoc$`prob.adj` < 0.1] = "."
    rphoc$sig[rphoc$`prob.adj` < 0.05] = "*"
    rphoc$sig[rphoc$`prob.adj` < 0.01] = "**"
    rphoc$sig[rphoc$`prob.adj` < 0.001] = "***"  
    
    if (any(class(phocresults) == "emmGrid")) {
        colnames(rphoc)[6] = "prob.adj"
        rownames(rphoc)=rphoc$contrast
        rphoc$sig = ""
        rphoc$sig[rphoc$`prob.adj` < 0.1] = "."
        rphoc$sig[rphoc$`prob.adj` < 0.05] = "*"
        rphoc$sig[rphoc$`prob.adj` < 0.01] = "**"
        rphoc$sig[rphoc$`prob.adj` < 0.001] = "***"  
        #rphoc$adj.method = phocresults@misc$adjust
        #rphoc = rphoc[, c(8, 4:7)]
        rphoc = rphoc[, c(4:7)]
        colnames(rphoc)[2] = substr(colnames(rphoc)[2],1,1)
        rphoc$`t` = round(rphoc$`t`, digits = 1)
    	rphoc$`df` = round(rphoc$`df`, digits = 1)
    	rphoc$`prob.adj` = round(rphoc$`prob.adj`, digits = digits)
        message("P value adjustment method: ",phocresults@misc$adjust)
        print(rphoc[,c(2,1,3,4)])
    } else{
    	rphoc$`t` = round(rphoc$`t`, digits = 1)
    	rphoc$`df` = round(rphoc$`df`, digits = 1)
    	rphoc$`prob.adj` = round(rphoc$`prob.adj`, digits = digits)
        rphoc
    }
}