随机平衡实验设计

Question

我正在编写一些代码来为市场研究生成平衡的实验设计，特别是用于联合分析和最大差异缩放。

第一步是生成部分平衡未完成块 (PBIB) 设计。这是直接使用 R 包AlgDesign。

对于大多数类型的研究，这样的设计就足够了。然而，在市场研究中，人们希望控制每个区块中的订单效应。这是我希望得到帮助的地方。

创建测试数据

# The following code is not essential in understanding the problem, 
# but I provide it in case you are curious about the origin of the data itself.
#library(AlgDesign)
#set.seed(12345)
#choices <- 4
#nAttributes <- 7
#blocksize <- 7
#bsize <- rep(choices, blocksize)
#PBIB <- optBlock(~., withinData=factor(1:nAttributes), blocksizes=bsize)
#df <- data.frame(t(array(PBIB$rows, dim=c(choices, blocksize))))
#colnames(df) <- paste("Item", 1:choices, sep="")
#rownames(df) <- paste("Set", 1:nAttributes, sep="")

df <- structure(list(
  Item1 = c(1, 2, 1, 3, 1, 1, 2), 
  Item2 = c(4, 4, 2, 5, 3, 2, 3), 
  Item3 = c(5, 6, 5, 6, 4, 3, 4), 
  Item4 = c(7, 7, 6, 7, 6, 7, 5)), 
  .Names = c("Item1", "Item2", "Item3", "Item4"), 
  row.names = c("Set1", "Set2", "Set3", "Set4", "Set5", "Set6", "Set7"), 
  class = "data.frame")

** 定义两个辅助函数

balanceMatrix计算矩阵的余额：

balanceMatrix <- function(x){
    t(sapply(unique(unlist(x)), function(i)colSums(x==i)))
}

balanceScore 计算“适合”的度量 - 分数越低越好，完美为零：

balanceScore <- function(x){
    sum((1-x)^2)
}

定义一个随机重新采样行的函数

findBalance <- function(x, nrepeat=100){
    df <- x
    minw <- Inf
    for (n in 1:nrepeat){
        for (i in 1:nrow(x)){df[i,] <- sample(df[i, ])}
        w <- balanceMatrix(df)
        sumw <- balanceScore(w)
        if(sumw < minw){
            dfbest <- df
            minw <- sumw
        }
    }
    dfbest
}

主要代码

dataframedf是7组的平衡设计。每组将向受访者展示 4 个项目。 df 中的数值指的是 7 个不同的属性。例如，在 Set1 中，将要求受访者从属性 1、3、4 和 7 中选择他/她的首选选项。

每个集合中项目的顺序在概念上并不重要。因此，(1,4,5,7) 的排序与 (7,5,4,1) 相同。

但是，为了获得完全平衡的设计，每个属性将在每列中出现相同的次数。这种设计是不平衡的，因为属性 1 在第 1 列中出现了 4 次：

df

     Item1 Item2 Item3 Item4
Set1     1     4     5     7
Set2     2     4     6     7
Set3     1     2     5     6
Set4     3     5     6     7
Set5     1     3     4     6
Set6     1     2     3     7
Set7     2     3     4     5

为了尝试找到更平衡的设计，我编写了函数findBalance。这通过在df 的行中随机抽样来随机搜索更好的解决方案。重复 100 次后，它会找到以下 最佳 解决方案：

set.seed(12345)
dfbest <- findBalance(df, nrepeat=100)
dfbest

     Item1 Item2 Item3 Item4
Set1     7     5     1     4
Set2     6     7     4     2
Set3     2     1     5     6
Set4     5     6     7     3
Set5     3     1     6     4
Set6     7     2     3     1
Set7     4     3     2     5

这看起来更平衡，并且计算出的平衡矩阵包含很多。平衡矩阵计算每个属性在每列中出现的次数。例如，下表表明（在左上角的单元格中）属性 1 在第 1 列中根本没有出现两次，在第 2 列中出现两次：

balanceMatrix(dfbest)

     Item1 Item2 Item3 Item4
[1,]     0     2     1     1
[2,]     1     1     1     1
[3,]     1     1     1     1
[4,]     1     0     1     2
[5,]     1     1     1     1
[6,]     1     1     1     1
[7,]     2     1     1     0

此解决方案的平衡分为 6，表示至少有 6 个不等于 1 的单元格：

balanceScore(balanceMatrix(dfbest))
[1] 6

我的问题

感谢您关注这个详细的示例。我的问题是如何将这个搜索功能重写为更系统？我想告诉 R：

• 最小化balanceScore(df)
• 通过更改df 的行顺序
• 受制于：已经完全受限

Answer 1

好的，我不知何故误解了你的问题。所以再见了 Fedorov，你好申请了 Fedorov。

以下算法基于 Fedorov 算法的第二次迭代：

1. 计算每个集合的所有可能排列，并将它们存储在 C0 列表中
2. 从 C0 空间绘制第一个可能的解决方案（每个集合一个排列）。这可以是原始的，但由于我需要索引，我宁愿随机开始。
3. 计算每个新解决方案的分数，其中第一组被所有排列替换。
4. 用给出最低分数的排列替换第一组
5. 每隔一组重复 3-4 次
6. 重复 3-5 直到分数达到 0 或进行 n 次迭代。

或者，您可以在 10 次迭代后重新启动该过程并从另一个起点开始。在您的测试用例中，结果证明很少有起点非常缓慢地收敛到 0。下面的函数在我的计算机上平均 1.5 秒内找到了得分为 0 的平衡实验设计：

> X <- findOptimalDesign(df)
> balanceScore(balanceMatrix(X))
[1] 0
> mean(replicate(20, system.time(X <- findOptimalDesign(df))[3]))
[1] 1.733

这就是现在的函数（给定你原来的 balanceMatrix 和 balanceScore 函数）：

findOptimalDesign <- function(x,iter=4,restart=T){
    stopifnot(require(combinat))
    # transform rows to list
    sets <- unlist(apply(x,1,list),recursive=F)
    nsets <- NROW(x)
    # C0 contains all possible design points
    C0 <- lapply(sets,permn)
    n <- gamma(NCOL(x)+1)

    # starting point
    id <- sample(1:n,nsets)
    Sol <- sapply(1:nsets,function(i)C0[[i]][id[i]])

    IT <- iter
    # other iterations
    while(IT > 0){
      for(i in 1:nsets){
          nn <- 1:n
          scores <- sapply(nn,function(p){
             tmp <- Sol
             tmp[[i]] <- C0[[i]][[p]]
             w <- balanceMatrix(do.call(rbind,tmp))
             balanceScore(w)
          })
          idnew <- nn[which.min(scores)]
          Sol[[i]] <- C0[[i]][[idnew]]

      }
      #Check if score is 0
      out <- as.data.frame(do.call(rbind,Sol))
      score <- balanceScore(balanceMatrix(out))
      if (score==0) {break}
      IT <- IT - 1

      # If asked, restart
      if(IT==0 & restart){
          id <- sample(1:n,nsets)
          Sol <- sapply(1:nsets,function(i)C0[[i]][id[i]])
          IT <- iter
      }
    }
    out
}

HTH

编辑：修复了小错误（因为我忘记了 IT 条件，它在每一轮后立即重新启动）。这样做，它的运行速度还是会快一些。