来自逻辑矩阵的紧凑型字母显示

Question

我想知道是否有一种方法可以将比较的逻辑矩阵转换为多重比较测试中使用的字母符号。如multcomp::cld。

我的数据是这样的：

test_data <- data.frame(mean=c(1.48, 1.59, 1.81,1.94),CI_lower=c(1.29,1.38,1.54, 1.62),CI_upper=c(1.56,1.84, 2.3, 2.59))

  mean CI_lower CI_upper
1 1.48     1.29     1.56
2 1.59     1.38     1.84
3 1.81     1.54     2.30
4 1.94     1.62     2.59

我感兴趣的是一个表示哪些条目具有重叠 CI 以获得如下所示的最终结果：

final <- data.frame(mean=c(1.48, 1.59, 1.81,1.94),CI_lower=c(1.29, 1.38,1.54, 1.62),CI_upper=c(1.56,1.84, 2.3, 2.59),letters = c("a","ab","ab","b"))

  mean CI_lower CI_upper letters
1 1.48     1.29     1.56       a
2 1.59     1.38     1.84      ab
3 1.81     1.54     2.30      ab
4 1.94     1.62     2.59       b

我做了一个可怜的尝试，结果是这样的：

same <- outer(test_data$CI_lower, test_data$CI_upper,"-")
same <- same<0
same <- lower.tri(same, diag = FALSE) & same

same_ind <- which(same,arr.ind = T)

groups <- as.list(as.numeric(rep(NA,nrow(test_data))))

for(i in 1:nrow(same_ind)){
  group_pos <- as.numeric(same_ind[i,])
  for(i2 in group_pos){
    groups[[i2]] <- c(groups[[i2]],i)
  }
}

letters_notation <- sapply(groups,function(x){
  x <- x[!is.na(x)]
  x <- letters[x]
  x <- paste0(x,collapse="")
  return(x)
}
)

这将给出：

  mean CI_lower CI_upper letters
1 1.48     1.29     1.56      ab
2 1.59     1.38     1.84     acd
3 1.81     1.54     2.30     bce
4 1.94     1.62     2.59      de

关于如何做到这一点的任何想法？

Answer 1

这是使用data.tables 非常有效的foverlaps 函数的可能解决方案。这不完全是您想要的输出（因为我不完全理解它），但您可以从中识别重叠点

library(data.table)
setkey(setDT(test_data), CI_lower, CI_upper)
Overlaps <- foverlaps(test_data, test_data, type = "any", which = TRUE) ## returns overlap indices
test_data[ , overlaps := Overlaps[, paste(letters[yid], collapse = ""), xid]$V1][]
#    mean CI_lower CI_upper overlaps
# 1: 1.48     1.29     1.56      abc <~~ not overlapping with d
# 2: 1.59     1.38     1.84     abcd
# 3: 1.81     1.54     2.30     abcd
# 4: 1.94     1.62     2.59      bcd <~~ not overlapping with a

Answer 2

根据 David Arenburg 的建议和这篇 http://menugget.blogspot.it/2014/05/automated-determination-of-distribution.html 不错的文章，我找到了解决方案。

library(igraph)

test_data <- data.frame(mean=c(1.48, 1.59, 1.81,1.94),CI_lower=c(1.29,1.38,1.54, 1.62),CI_upper=c(1.56,1.84, 2.3, 2.59))

n <- nrow(test_data)


g <- outer(test_data$CI_lower, test_data$CI_upper,"-")
g <- !(g<0)
g <- g + t(g) # not necessary, but make matrix symmetric
g <- g!=1
rownames(g) <- 1:n # change row names
colnames(g) <- 1:n # change column names


# Re-arrange data into an "edge list" for use in igraph (i.e. which groups are "connected") - Solution from "David Eisenstat" ()
same <- which(g==1)
g2 <- data.frame(N1=((same-1) %% n) + 1, N2=((same-1) %/% n) + 1)
g2 <- g2[order(g2[[1]]),] # Get rid of loops and ensure right naming of vertices
g3 <- simplify(graph.data.frame(g2,directed = FALSE))



# Calcuate the maximal cliques - these are groupings where every node is connected to all others
cliq <- maximal.cliques(g3) # Solution from "majom" ()
cliq2 <- lapply(cliq, as.numeric)

# Reorder by level order - Solution from "MrFlick" ()
ml<-max(sapply(cliq, length))
reord <- do.call(order, data.frame(
  do.call(rbind, 
          lapply(cliq2, function(x) c(sort(x), rep.int(0, ml-length(x))))
  )
))
cliq <- cliq[reord]
cliq

# Generate labels to  factor levels
lab.txt <- vector(mode="list", n) # empty list
lab <- letters[seq(cliq)] # clique labels
for(i in seq(cliq)){ # loop to concatenate clique labels
  for(j in cliq[[i]]){
    lab.txt[[j]] <- paste0(lab.txt[[j]], lab[i])
  }
}

 

unlist(lab.txt)
[1] "a"  "ab" "ab" "b"