随机平衡实验设计

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
}

主码

数据框 df是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，表示至少有六个不等于 1 的单元格:

balanceScore(balanceMatrix(dfbest))
[1] 6

我的问题

感谢您关注此详细示例。我的问题是如何重写此搜索功能以使其更系统？我想告诉 R:

最小化 balanceScore(df)

通过更改 df 的行顺序

受制于:已完全受约束

Answer 1

好的，我以某种方式误解了你的问题。所以再见费多罗夫，你好申请费多罗夫。

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

计算每个集合的所有可能排列，并将它们存储在 C0 列表中

从 C0 空间中得出第一个可能的解(每组一个排列)。这可以是原始的，但由于我需要索引，我宁愿随机开始。

计算每个新解决方案的分数，其中第一组被所有排列替换。

用得分最低的排列替换第一组

每隔一组重复 3-4

重复 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 的条件)。这样做，它仍然运行得更快一些。