Блог

Я пишу некоторый код для генерации сбалансированных экспериментальных проектов для исследования рынка, специально для использования в совместном анализе и масштабировании максимальной разницы.

Первым шагом является создание частично сбалансированного проекта с незавершенным блоком (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 различным атрибутам. Например, в наборе 1 респонденту будет предложено выбрать предпочитаемый им вариант из атрибутов 1, 3, 4 и 7.

Порядок элементов в каждом наборе концептуально не важен. Таким образом, порядок (1,4,5,7) идентичен порядку (7,5,4,1).

Однако для получения полностью сбалансированного проекта каждый атрибут будет отображаться равное количество раз в каждом столбце. Этот проект является несбалансированным, поскольку атрибут 1 появляется 4 раза в столбце 1:

 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
• При условии: уже полностью ограничено

Хорошо, я как-то неправильно понял ваш вопрос. Так что пока, Федоров, привет прикладному Федорову.

Следующий алгоритм основан на второй итерации алгоритма Федорова :

1. вычислите все возможные перестановки для каждого набора и сохраните их в списке C0
2. нарисуйте первое возможное решение из пространства C0 (по одной перестановке для каждого набора). Это может быть оригинальным, но поскольку мне нужны индексы, я предпочитаю просто начинать наугад.
3. вычисляйте оценку для каждого нового решения, где первый набор заменяется всеми перестановками.
4. замените первый набор на перестановку, дающую наименьший балл
5. повторите 3-4 для каждого другого набора
6. повторяйте 3-5 раз, пока оценка не достигнет 0 или для n итераций.

При желании вы можете перезапустить процедуру после 10 итераций и начать с другой начальной точки. В вашем тестовом примере оказалось, что несколько начальных точек очень медленно сходились к 0. Приведенная ниже функция обнаружила сбалансированные экспериментальные проекты со счетом 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 amp; restart){
          id <- sample(1:n,nsets)
          Sol <- sapply(1:nsets,function(i)C0[[i]][id[i]])
          IT <- iter
      }
    }
    out
}
  

HTH

РЕДАКТИРОВАТЬ: исправлена небольшая ошибка (она перезапускалась сразу после каждого раунда, так как я забывал ставить на нее условие). При этом он работает еще немного быстрее.