Ограниченная рандомизация в трехмерном массиве в R для экспериментального дизайна

# Создание пустого трехмерного массива
tube_array <- array(dim = c(4, 7, 14))
dimnames(tube_array) <- list(
 replicate = paste0("R", 1:4),
 block = paste0("B", 1:7),
 depth = paste0("D", 1:14)
)

# Пример генерации candidate set
generate_candidate_set <- function(n_farms, n_tubes, n_depths, treatments) {
 candidates <- list()
 max_attempts <- 1000
 attempt <- 0
 
 while(length(candidates) < 100 && attempt < max_attempts) {
 # Генерация одной схемы
 scheme <- generate_random_scheme(n_farms, n_tubes, n_depths, treatments)
 
 # Проверка ограничений
 if(check_constraints(scheme)) {
 candidates[[length(candidates) + 1]] <- scheme
 }
 
 attempt <- attempt + 1
 }
 
 return(candidates)
}

# Пример использования designit
library(designit)

# Определение контейнера пакетов
container <- list(
 tubes = rep(1:28, each = 14),
 depths = rep(1:14, 28)
)

# Функция оценки
evaluation_function <- function(design) {
 # Расчет баланса по глубинам и трубкам
 balance_score <- calculate_balance(design)
 return(balance_score)
}

# Генерация дизайна
balanced_design <- generate_balanced_design(
 container = container,
 evaluation_function = evaluation_function,
 n_schemes = 1000
)

# Пример использования cvcrand
library(cvcrand)

# Ограничение пространства рандомизации
constrained_design <- cvcrand(
 treatment_vector = treatment_assignments,
 covariates = farm_covariates,
 strata = strata_definition,
 q = 0.1 # Квантиль для ограничения
)

# Претод отклонения для предотвращения повторов ферм в одной трубке
farm_rejection_sampling <- function(tube_assignments, farm_ids) {
 valid_assignments <- FALSE
 max_attempts <- 1000
 attempt <- 0
 
 while(!valid_assignments && attempt < max_attempts) {
 # Случайное назначение ферм
 random_assignment <- sample(farm_ids)
 
 # Проверка отсутствия повторов в трубках
 valid_assignments <- check_farm_duplicates(tube_assignments, random_assignment)
 
 attempt <- attempt + 1
 }
 
 return(random_assignment)
}

# Определение требуемого распределения по глубинам
required_depth_distribution <- function(depths, treatments) {
 # Рассчитаем требуемое количество образцов каждого типа на каждой глубине
 depth_dist <- data.frame(
 depth = depths,
 treatment_A = rep(10.8, length(depths)), # 152/14 = 10.857
 treatment_B = rep(10.8, length(depths)),
 treatment_C = rep(6.4, length(depths)) # 88/14 = 6.285
 )
 
 # Корректировка для целочисленности
 depth_dist$treatment_A[1] <- 11 # 152 - 13*10.8 = 11
 depth_dist$treatment_B[1] <- 11 # 152 - 13*10.8 = 11
 depth_dist$treatment_C[1] <- 6 # 88 - 13*6.4 = 6
 
 return(depth_dist)
}

improved_constrained_randomization <- function(n_farms, n_tubes, n_depths, 
 treatments, farm_ids) {
 # Инициализация трехмерного массива
 result_array <- array(dim = c(n_tubes, 1, n_depths))
 
 # Шаг 1: Определение распределения по глубинам
 depth_dist <- required_depth_distribution(1:n_depths, treatments)
 
 # Шаг 2: Распределение обработок по глубинам с соблюдением баланса
 for(depth in 1:n_depths) {
 # Количество образцов каждого типа на текущей глубине
 n_A <- depth_dist$treatment_A[depth]
 n_B <- depth_dist$treatment_B[depth]
 n_C <- depth_dist$treatment_C[depth]
 
 # Создание вектора обработок для текущей глубины
 depth_treatments <- c(
 rep("A", n_A),
 rep("B", n_B),
 rep("C", n_C)
 )
 
 # Случайное перемешивание обработок
 depth_treatments <- sample(depth_treatments)
 
 # Назначение обработок на текущей глубине
 for(tube in 1:n_tubes) {
 result_array[tube, 1, depth] <- depth_treatments[tube]
 }
 }
 
 # Шаг 3: Назначение ферм с учетом ограничения на повторы
 farm_assignments <- array(dim = c(n_tubes, 1, n_depths))
 used_farms <- list()
 
 for(depth in 1:n_depths) {
 for(tube in 1:n_tubes) {
 treatment <- result_array[tube, 1, depth]
 
 # Получаем доступные фермы для текущей обработки
 available_farms <- get_available_farms(treatment, farm_ids, used_farms, tube, depth)
 
 if(length(available_farms) > 0) {
 # Случайный выбор фермы
 selected_farm <- sample(available_farms, 1)
 farm_assignments[tube, 1, depth] <- selected_farm
 
 # Обновление списка использованных ферм
 if(!tube %in% names(used_farms)) {
 used_farms[[as.character(tube)]] <- list()
 }
 used_farms[[as.character(tube)]][[as.character(depth)]] <- selected_farm
 } else {
 # Восстановление баланса при невозможности назначения
 restore_balance(result_array, farm_assignments, tube, depth, used_farms)
 }
 }
 }
 
 return(list(
 treatments = result_array,
 farms = farm_assignments
 ))
}

get_available_farms <- function(treatment, farm_ids, used_farms, tube, depth) {
 # Получаем фермы для текущего типа обработки
 treatment_farms <- farm_ids[farm_ids$treatment == treatment]
 
 # Исключаем уже использованные фермы в текущей трубке
 if(as.character(tube) %in% names(used_farms)) {
 used_tube_farms <- unlist(used_farms[[as.character(tube)]])
 available_farms <- setdiff(treatment_farms, used_tube_farms)
 } else {
 available_farms <- treatment_farms
 }
 
 return(available_farms)
}

restore_balance <- function(result_array, farm_assignments, tube, depth, used_farms) {
 # В случае невозможности назначения фермы, перестраиваем распределение
 # для текущей глубины
 
 # Получаем текущее распределение обработок на глубине
 depth_treatments <- result_array[, 1, depth]
 
 # Пересчитываем требуемое распределение
 current_counts <- table(depth_treatments)
 required_counts <- c(A = 11, B = 11, C = 6) # Для первой глубины
 
 # Корректируем распределение обработок
 for(treatment in names(required_counts)) {
 if(current_counts[treatment] < required_counts[treatment]) {
 # Добавляем образцы этого типа обработки
 deficit <- required_counts[treatment] - current_counts[treatment]
 positions <- which(depth_treatments == treatment)
 
 # Заменяем другие типы обработки на нужный
 replace_positions <- sample(which(depth_treatments != treatment), deficit)
 depth_treatments[replace_positions] <- treatment
 }
 }
 
 # Обновляем массив обработок
 result_array[, 1, depth] <- depth_treatments
 
 # Повторно назначаем фермы
 for(t in 1:length(depth_treatments)) {
 treatment <- depth_treatments[t]
 available_farms <- get_available_farms(treatment, farm_ids, used_farms, t, depth)
 
 if(length(available_farms) > 0) {
 selected_farm <- sample(available_farms, 1)
 farm_assignments[t, 1, depth] <- selected_farm
 
 if(!as.character(t) %in% names(used_farms)) {
 used_farms[[as.character(t)]] <- list()
 }
 used_farms[[as.character(t)]][[as.character(depth)]] <- selected_farm
 }
 }
}

# Параметры эксперимента
n_tubes <- 28 # 4×7
n_depths <- 14
n_farms <- 152
treatments <- c(rep("A", 152), rep("B", 152), rep("C", 88))

# Создание идентификаторов ферм
farm_data <- data.frame(
 farm_id = 1:n_farms,
 treatment = treatments
)

# Запуск улучшенной рандомизации
randomization_result <- improved_constrained_randomization(
 n_farms = n_farms,
 n_tubes = n_tubes,
 n_depths = n_depths,
 treatments = treatments,
 farm_ids = farm_data
)

# Проверка результатов
validate_randomization(randomization_result, farm_data)

# Проверка баланса по трубкам
validate_tube_balance <- function(result_array, expected_counts) {
 # Агрегация по трубкам
 tube_counts <- apply(result_array, 1, function(x) {
 table(factor(x, levels = c("A", "B", "C")))
 })
 
 # Преобразование в data frame
 tube_balance <- data.frame(
 tube = 1:nrow(tube_counts),
 treatment_A = sapply(tube_counts, function(x) x["A"]),
 treatment_B = sapply(tube_counts, function(x) x["B"]),
 treatment_C = sapply(tube_counts, function(x) x["C"])
 )
 
 # Расчет отклонений от ожидаемых значений
 tube_balance$deviation_A <- abs(tube_balance$treatment_A - expected_counts$A)
 tube_balance$deviation_B <- abs(tube_balance$treatment_B - expected_counts$B)
 tube_balance$deviation_C <- abs(tube_balance$treatment_C - expected_counts$C)
 
 # Общая оценка баланса
 tube_balance$balance_score <- rowMeans(tube_balance[, c("deviation_A", "deviation_B", "deviation_C")])
 
 return(tube_balance)
}

# Проверка баланса по глубинам
validate_depth_balance <- function(result_array, expected_counts) {
 # Агрегация по глубинам
 depth_counts <- apply(result_array, 3, function(x) {
 table(factor(x, levels = c("A", "B", "C")))
 })
 
 # Преобразование в data frame
 depth_balance <- data.frame(
 depth = 1:ncol(depth_counts),
 treatment_A = sapply(depth_counts, function(x) x["A"]),
 treatment_B = sapply(depth_counts, function(x) x["B"]),
 treatment_C = sapply(depth_counts, function(x) x["C"])
 )
 
 # Расчет отклонений от ожидаемых значений
 depth_balance$deviation_A <- abs(depth_balance$treatment_A - expected_counts$A)
 depth_balance$deviation_B <- abs(depth_balance$treatment_B - expected_counts$B)
 depth_balance$deviation_C <- abs(depth_balance$treatment_C - expected_counts$C)
 
 # Общая оценка баланса
 depth_balance$balance_score <- rowMeans(depth_balance[, c("deviation_A", "deviation_B", "deviation_C")])
 
 return(depth_balance)
}

# Проверка отсутствия повторных ферм в одной трубке
validate_farm_uniqueness <- function(farm_assignments) {
 n_tubes <- dim(farm_assignments)[1]
 n_depths <- dim(farm_assignments)[3]
 violations <- 0
 
 for(tube in 1:n_tubes) {
 # Получаем все фермы в текущей трубке
 tube_farms <- farm_assignments[tube, 1, ]
 
 # Проверяем дубликаты
 unique_farms <- unique(tube_farms)
 if(length(unique_farms) != n_depths) {
 violations <- violations + (n_depths - length(unique_farms))
 }
 }
 
 return(violations)
}

# Общая функция валидации
validate_randomization <- function(result_array, farm_assignments, expected_counts) {
 # Проверка баланса по трубкам
 tube_validation <- validate_tube_balance(result_array, expected_counts)
 
 # Проверка баланса по глубинам
 depth_validation <- validate_depth_balance(result_array, expected_counts)
 
 # Проверка уникальности ферм
 farm_violations <- validate_farm_uniqueness(farm_assignments)
 
 # Общая оценка
 overall_score <- list(
 tube_balance = mean(tube_validation$balance_score),
 depth_balance = mean(depth_validation$balance_score),
 farm_violations = farm_violations,
 total_score = mean(c(mean(tube_validation$balance_score), 
 mean(depth_validation$balance_score), 
 farm_violations))
 )
 
 return(overall_score)
}