随机平衡实验设计

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
}

主代码

dataframe 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，表示至少有6个单元格等于1：

balanceScore(balanceMatrix(dfbest))
[1] 6

我的问题

谢谢你跟随这个详细的例子。我的问题是，如何重写这个搜索函数，使其更加系统化？我想告诉R：

最小化balanceScore(df)

• By将df

• Subject的行顺序更改为:已经完全受限的

Answer 1

好吧，我有点误解了你的问题。再见，费多罗夫，你好，费多罗夫。

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

1. 计算每个集合的所有可能排列，并将它们存储在C0 list
2. 中，从C0空间绘制出第一个可能的解决方案(每组一个置换)。这可以是最初的，但由于我需要索引，我宁愿从random.
3. calculate开始，为每一个新的解决方案评分，其中第一组被替换为所有排列。
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

编辑:修正了小错误(它在每一轮后立即重新启动，因为我忘记了对它的条件)。这样做，它运行得更快一些。