Aleatorización de diseños experimentales equilibrados

Question

Estoy escribiendo algunos códigos para generar diseños experimentales equilibrados para la investigación de mercado, específicamente para su uso en análisis conjunto y escalas de diferencia máxima.

El primer paso es generar un diseño de Bloque Incompleto Parcialmente Equilibrado (PBIB). Esto es sencillo usando el paquete R AlgDesign.

Para la mayoría de los tipos de investigación, un diseño de este tipo sería suficiente. Sin embargo, en la investigación de mercado uno quiere controlar los efectos de orden en cada bloque. Aquí es donde apreciaría algo de ayuda.

Crear datos de prueba

# 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") 

** definir dos funciones de ayuda

balanceMatrix calcula el balance de la matriz:

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

balanceScore calcula una métrica de 'ajuste' - los puntajes más bajos son mejores, con cero perfecto:

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

definir una función que vuelve a muestrear las filas al azar

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 
} 

código principal

La trama de datos df es un diseño equilibrado de 7 conjuntos. Cada conjunto mostrará 4 elementos al encuestado. Los valores numéricos en df se refieren a 7 atributos diferentes. Por ejemplo, en el Conjunto 1 se le pedirá al encuestado que elija su opción preferida de los atributos 1, 3, 4 y 7.

El orden de los elementos en cada conjunto no es conceptualmente importante. Por lo tanto, un orden de (1,4,5,7) es idéntico a (7,5,4,1).

Sin embargo, para obtener un diseño completamente equilibrado, cada atributo aparecerá la misma cantidad de veces en cada columna. Este diseño se había desequilibrado, ya que el atributo 1 aparece 4 veces en la columna 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 

para tratar de encontrar un diseño más equilibrado, que escribí la función findBalance. Esto lleva a cabo una búsqueda aleatoria de mejores soluciones, muestreando aleatoriamente en las filas de df. Con 100 repeticiones que encuentra la siguiente solución mejor:

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 

Esto parece más equilibrado, y la matriz de equilibrio calculado contiene gran cantidad de queridos. La matriz de saldo cuenta el número de veces que aparece cada atributo en cada columna.Por ejemplo, la siguiente tabla indica (en la celda superior izquierda) que atributo 1 no aparece dos veces en absoluto en la columna 1, y dos veces en la columna 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 

El equilibrio puntuación para esta solución es 6 , lo que indica que hay al menos seis células desiguales a 1:

balanceScore(balanceMatrix(dfbest)) 
[1] 6 

Mi pregunta

Gracias por seguir este ejemplo detallado. Mi pregunta es ¿cómo puedo reescribir esta función de búsqueda para que sea más sistemática? Me gustaría decirle a R:

• Minimizar balanceScore(df)
• Cambiando orden de las filas de df
• Sujeto a: ya está totalmente restringida

Answer 1

bien, de alguna manera entendido bien su pregunta. Así que adiós Fedorov, hola aplicó Fedorov.

El siguiente algoritmo se basa en la segunda iteración del algoritmo Fedorov:

1. calcular todas las posibles permutaciones para cada conjunto, y almacenarlas en la lista C0
2. dibujar una primera solución posible del C0 espacio (una permutación para cada conjunto). Este puede ser el original, pero como necesito los índices, prefiero comenzar al azar.
3. calcule el puntaje para cada solución nueva, donde el primer conjunto se reemplaza por todas las permutaciones.
4. reemplazar el primer conjunto con la permutación dando la puntuación más baja
5. repetición 3-4 para cada otro conjunto
6. repetición 3-5 hasta que la puntuación llega a 0 o para n iteraciones.

Opcionalmente, puede reiniciar el procedimiento después de 10 iteraciones y comenzar desde otro punto de inicio. En le caso de prueba, resultó que algunos puntos de partida convergieron muy lentamente a 0. La función más adelante encontró diseños experimentales equilibradas con una puntuación de 0 en un promedio de 1,5 segundos en mi equipo:

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

Así que este es el ahora la función (dado su balanceMatrix original y funciones 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

EDIT: fijo pequeño fallo (se reinicia inmediatamente después de cada ronda que se me olvidó para acondicionar en él). Al hacerlo, corre un poco más rápido todavía.