随机化平衡的实验设计

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,你好申请Fedorov.

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

1. 计算每个集合的所有可能的排列,并将它们存储在C0列表中
2. 从C0空间中提取第一个可能的解决方案(每组一个排列).这可能是原始的,但由于我需要索引,我宁愿随意开始.
3. 计算每个新解决方案的得分,其中第一组被所有排列替换.
4. 用排列取代第一组得分最低
5. 每隔一组重复3-4次
6. 重复3-5,直到得分达到0或n次迭代.

(可选)您可以在10次迭代后重新启动该过程,并从另一个起点开始.在你的测试案例中,事实证明很少有起点非常缓慢地收敛到0.下面的函数在我的计算机上找到了平均1.5秒的平衡实验设计:

> 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上的条件).这样做,它仍然运行得更快.