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)
{
  # written by Jonathan Harrington and Ulrich Reubold (phsel.emmGrid)
  
  # 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(mmemmGrid, k = 1)
    {
      #written by Ulrich Reubold
      #mmemmGrid needs to be result of pairs(emmeans())
      #from package emmeans with a model fitted by lmer()
      #containing a two-way interaction term, e.g.:
      #stimm.lmer = lmer(vot ~ K * Alter + (K | Vpn), data = stimm)
      #mmphoc = pairs(emmeans(stimm.lmer,~K:Alter))
      ph3 <- data.frame(mmemmGrid)
      rownames(ph3) = sub(pattern = ",", replacement = ":", ph3$contrast)
      rownames(ph3) = sub(pattern = ",",
                          replacement = ":",
                          rownames(ph3))
      rownames(ph3) = sub(pattern = " ",
                          replacement = "",
                          rownames(ph3))
      rownames(ph3) = sub(pattern = " ",
                          replacement = "",
                          rownames(ph3))
      ph3 = ph3[, c(5, 4, 6)]
      colnames(ph3) = c("t", "df", "prob.adj")
      m = rownames(ph3)
      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)
      }
      rphoc <- ph3[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 adding stars and name of adjustment
  
  
  if (any(class(phocresults) == "emmGrid")) {
    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(5, 1, 2, 3, 4)]
    rphoc
  } else{
    rphoc
  }
}


. <-
  function (..., .env = parent.frame())
  {
    structure(as.list(match.call()[-1]), env = .env, class = "quoted")
  }