R Дизайн экспериментальной модели выбора на основе AICc

Question

У меня довольно простая проблема. Я пытаюсь автоматизировать выбор линейной модели из плана эксперимента на основе исправленного информационного критерия Акаике (=AICc). Для меня важно использовать AICc вместо AIC, в противном случае выбор модели будет иметь тенденцию к переобучению из-за небольшого количества выборок.

В конце я хочу иметь уменьшенную модель.

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

Я предоставляю вам мой код вместе с примерами данных.

Сначала приведем пример данных для дизайна с 5 факторами (A:E) и двумя ответами (y и y2). Я хочу оценить основные эффекты, взаимодействия (только между основными эффектами) и квадратичные эффекты.

#topic:   DoE analysis
#author:  xxxx


#####################################################################################################################################
####preparation######################################################################################################################
#####################################################################################################################################


#packages
library(glmulti)
library(tidyverse)


#####################################################################################################################################
####options##########################################################################################################################
#####################################################################################################################################


#ranges##############################################################################################################################

#range for factors
rf <- 1:5
#range for responses
y <- 6:7


#####################################################################################################################################
####load data########################################################################################################################
#####################################################################################################################################
df <- data.frame(
                 A = c(-1L, -1L, -1L, 1L, -1L, -1L, 1L, 1L, 0L, 1L, -1L, 0L, 1L, -1L,
                       -1L, 1L, -1L, 1L, 1L, -1L, -1L, -1L, 1L, 0L, 0L, -1L,
                       0L, 1L, 1L, 1L),
                 B = c(1L, -1L, 1L, -1L, -1L, 0L, 0L, -1L, 0L, -1L, 0L, -1L, 1L, 1L,
                       1L, -1L, 1L, 1L, 1L, 1L, -1L, -1L, -1L, 1L, 1L, -1L,
                       -1L, 1L, 0L, 1L),
                 C = c(-1L, 0L, 1L, 1L, 1L, -1L, 1L, -1L, -1L, -1L, 1L, 0L, 0L, -1L,
                       -1L, -1L, 1L, -1L, 1L, 0L, -1L, 1L, 1L, -1L, 1L, -1L,
                       1L, 1L, 0L, -1L),
                 D = c(1L, 0L, 1L, -1L, 1L, 1L, 0L, -1L, -1L, 1L, -1L, -1L, -1L, -1L,
                       0L, 0L, -1L, -1L, -1L, 1L, -1L, -1L, 1L, 1L, 0L, 1L, 1L,
                       1L, 1L, 1L),
                 E = c(0L, -1L, -1L, 1L, 1L, 1L, 1L, -1L, -1L, 1L, 0L, 1L, 0L, -1L,
                       1L, 0L, 1L, 1L, -1L, 1L, 1L, -1L, -1L, 1L, -1L, -1L, 0L,
                       1L, -1L, -1L),
                 y = c(97.73, 106.6, 102.5, 100.2, 119.6, 105.9, 103.1, 112.1, 100.7,
                       109.2, 104.9, 97.78, 94.27, 96.65, 102, 109.9, 111.8,
                       87.58, 105.5, 103.8, 92.11, 105.2, 141.5, 98.02, 107, 113,
                       127.1, 96.26, 114.3, 93.15),
                y2 = c(114.9, 116.2, 114.9, 118.3, 113.5, 115.3, 116.9, 116.6, 115.9,
                       118.3, 114.8, 114.9, 116.2, 114.5, 115.5, 116.6, 114.3,
                       116.6, 116.2, 116.1, 114.9, 115.5, 119, 116.2, 115.7,
                       117.3, 115.4, 117.2, 115.7, 118.5)
      )

После этого я извлекаю имена столбцов и создаю векторы для формулы модели.

#####################################################################################################################################
####data preparation#################################################################################################################
#####################################################################################################################################


### Define which columns are factors and which are responses

#extract column names
namesFactors <- colnames(df)[rf]
namesResponses <- colnames(df)[y]

#extract values
factors <- as.matrix(df[,namesFactors])
responses <- df[,namesResponses]

#create formula
mainFactors <- namesFactors
quadFactors <- paste0("I(", namesFactors, "^2)")

quadquadInteraction <- paste0(combn(quadFactors,2)[1,], ":", combn(quadFactors,2)[2,]) 
mainquadInteraction <- paste0(expand.grid(mainFactors, quadFactors)[,2], ":", expand.grid(mainFactors, quadFactors)[,1])

excludeFactors <- c(quadquadInteraction, mainquadInteraction)

Наконец, я пытаюсь найти сокращенную модель на основе AICc:

#reduce model
redModel <- vector("list",length(namesResponses))

for(i in 1:length(namesResponses)){
  redModel[[i]] <- glmulti(y = namesResponses[i], c(mainFactors, quadFactors), data = df[c(namesFactors, namesResponses[i])], crit="aicc", marginality=TRUE, level = 2, fitfunction = "lm", method = "g", deltaM = 10, deltaB = 0, conseq = 3, exclude = excludeFactors)
}

Эта последняя часть немного усложняется. Glmulti допускает только модели с основными эффектами (уровень = 1) или основными эффектами и взаимодействиями (уровень = 2). Чтобы включить квадратичные термины, я должен добавить их как I(Фактор ^2). Таким образом, они будут рассматриваться как основные эффекты. Однако glmulti также будет искать квадратичные взаимодействия и взаимодействия между квадратичными членами и основными эффектами. Но я не заинтересован в этих эффектах! По этой причине я определил вектор с именем exludeFactors, где я рассматриваю эти нежелательные термины. К сожалению, glmulti допускает только исключение некоторых из этих терминов, но не всех, и я не знаю почему!

Если бы кто-то мог помочь мне запустить код, я был бы очень благодарен. Я также могу использовать любой другой пакет и функцию вместо glmulti, при условии, что информационным критерием (=IC) является AICc, а не AIC.

Если кто-то знает, как уменьшить полную линейную модель, основываясь на функции r-base lm(), мы могли бы использовать этот код для соответствия полной модели:

#interactions
inter <- TRUE
quad <- TRUE

#create formula
formula <- vector("list",length(namesResponses))

for(i in 1:length(namesResponses)){
  formula[[i]] <-  paste(namesResponses[i], "~") %>%  #reponse part
    paste(".") %>% #main effects
    ifelse(inter == TRUE, paste(.," + (.)^2"), .) %>% #interactions
    ifelse(quad == TRUE, paste(.,paste0(" + I(", namesFactors, "^2)", collapse="")), .) %>% #quadratic terms
    as.formula
}

#create lm models based on 1t idea
fullModel <- vector("list",length(namesResponses))

for(i in 1:length(namesResponses)){
  fullModel[[i]] <- lm(data = df[c(namesFactors,namesResponses[i])], formula[[i]])
}