Анализ с использованием линейной регрессии на основе подгрупп

Question

Предположим, у меня есть данные (t,y), где я ожидаю линейную зависимость y(t). Кроме того, существуют атрибуты для каждого наблюдения par1, par2, par3. Есть ли алгоритм или метод для решения, если (один или оба или все параметры) релевантны для соответствия или нет? Я попробовал leaps::regsubsets(y ~ t + par1 + par2 + par3, data = mydata, nbest = 10), но не смог получить формулу для наилучшего соответствия.

Окончательный результат должен выглядеть так, как показано на рисунке. Данные см. Ниже.

Таким образом, я хочу информацию

• Добавление par1 и par2 дает наилучшим образом подходит
• Модели y_i = a_i * t_i + b_i с данным a_i и b_i

Возпроизводимо пример:

t <- seq(0,10, length.out = 1000) # large sample of x values 
# Create 3 linear equations of the form y_i = a*t_i + b 
a <- c(1, 0.3, 0.2) # slope 
b <- c(-0.5, 0.5, 0.1) # offset 

# create t_i, y_ti and y_i (including noise) 
d <- list() 
y <- list() 
y_t <- list() 
for (i in 1:3) { 
    set.seed(33*i) 
    d[[i]] <- sort(sample(t, 50, replace = F)) 
    set.seed(33*i) 
    noise <- rnorm(10) 
    y[[i]] <- a[i]*d[[i]] + b[i] + noise 
    y_t[[i]] <- a[i]*d[[i]] + b[i] 
} 
# Final data set 
df1 <- data.frame(t=d[[1]], y=y[[1]], par1=rep(1), par2=rep(10), par3=sample(c(100, 200, 300), length(d[[1]]), replace = T)) 
df2 <- data.frame(t=d[[2]], y=y[[2]], par1=rep(2), par2=rep(20), par3=sample(c(100, 200, 300), length(d[[1]]), replace = T)) 
df3 <- data.frame(t=d[[3]], y=y[[3]], par1=rep(2), par2=rep(30), par3=sample(c(100, 200, 300), length(d[[1]]), replace = T)) 
mydata <- rbind(df1, df2, df3) 
mydata <- mydata[sample(nrow(mydata)), ] 

# That is what the data is looking like: 
plot(mydata$t, mydata$y) 

# This is the result I am looking for (ideally): 
plot(d[[1]], y[[1]], col = "black", xlim = c(0, 10), ylim = c(-2, 10), xlab = "t", ylab = "y", 
    main = "Fit for three different groups") 
points(d[[2]], y[[2]], col = "red") 
points(d[[3]], y[[3]], col = "blue") 
lines(d[[1]], y_t[[1]],col = "black") 
lines(d[[2]], y_t[[2]], col = "red") 
lines(d[[3]], y_t[[3]], col = "blue") 

комментарий и вопрос о @ Роланда ответ:

Я понимаю, что это с заданными тремя параметрами есть 2^3=8 группы с 2*3*3=18 уровнями факторов. Но я бы ожидал, что у нас есть только 8 релевантных групп, поскольку у меня всегда есть выбор между «include parameter x or not». Для меня не имеет смысла только «включать уровень x параметра y».

Я попробовал следующее

g <- 0 
t_lin1 <- mydata$t[mydata$g == g] 
y_lin1 <- mydata$y[mydata$g == g] 
plot(mydata$t, mydata$y) 
points(t_lin1, y_lin1, col = "red") 
abline(lm(y_lin1 ~ t_lin1), col = "red") 
points(pred.1se ~ t, data = mydata, col = as.integer(mydata$g), pch = 16) 

и понял, что подгонка выключена. Оглядываясь назад, это ясно, потому что

• I включают неправильные уровни фактора (скорее всего, параметр 3 не имеет отношения)
• и таким образом получить неверные данные для подгонки

Так что мой последний вопрос :

• Где я могу найти соответствующие группы, включенные в лучшую модель и
• Каковы соответствующие параметры соответствия от регрессии?

К сожалению, если бы это было очевидно, но для меня это тайна

Answer 1

лассо может прийти довольно близко (хотя он идентифицирует все еще слишком много эффектов):

#I assume these are supposed to be factors: 
mydata$par1 <- factor(mydata$par1) 
mydata$par2 <- factor(mydata$par2) 
mydata$par3 <- factor(mydata$par3) 

#create model matrix, remove intercept since glmnet adds it 
x <- model.matrix(y ~ (par1 * par2 * par3) * t, data = mydata)[,-1] 

#cross-validated LASSO 
library(glmnet) 
set.seed(42) 
fit <- cv.glmnet(x, mydata$y, intercept = TRUE, nfolds = 10, alpha = 1) 
plot(fit) 

coef <- as.matrix(coef(fit, s = "lambda.1se")) 
coef[coef != 0,] 
#(Intercept)  par230   t  par12:t par230:t par3300:t 
# 0.47542479 -0.27612966 0.75497711 -0.42493030 -0.15044371 0.03033057 

#The groups: 
mydata$g <- factor((mydata$par2 == 30) + 10 * (mydata$par1 == 2) + 100 * (mydata$par3 == 300)) 



mydata$pred.1se <- predict(fit, newx = x, s = "lambda.1se") 

library(ggplot2) 
ggplot(mydata, aes(x = t, color = g)) + 
    geom_point(aes(y = y)) + 
    geom_line(aes(y = pred.1se)) 

Затем вы можете рассчитать требуемые перехваты и наклоны от коэффициентов.