Интерпретация подходящих коэффициентов MARS
Я пытаюсь подогнать модель многомерного адаптивного регрессионного сплайна, используя earthRpackage, к моим данным и нахождение интерпретации коэффициентов немного сбивает с толку.
Мои данные представляют собой измерения экспрессии генов с течением времени в 14 временных точках, по 2 выборки от каждого пола на каждый временной момент. Вот пример одного гена:
df1 <- data.frame(y = c(0.145, 0.1, 0.144, -0.002, -0.031, 0.062, 0.089, -0.029, 0.046, -0.008, 0.001, -0.147, -0.009, -0.021, 0.169, 0.21, 0.141, 0.095, -0.125, 0.156, 0.156, 0.252, 0.135, 0.123, -0.037, -0.154, -0.162, -0.112, 0.026, 0.069, -0.068, 0.006, -0.349, -0.07, -0.147, -0.079, -0.146, -0.145, -0.377, -0.273, -0.176, -0.204, -0.013, 0.116, 0.033, 0.049, -0.104, -0.086, -0.068, -0.096, -0.056, -0.016, -0.281, -0.277, -0.29, -0.203),
t = rep(c(-1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39),4),
sex = c(rep("F",28),rep("M",28)))
df1$sex <- factor(df1$sex, levels = c("M","F"))
Который выглядит как:
library(ggplot2)
ggplot(df1, aes(y = y, x = t, color = sex)) + geom_point(size = 3) + theme_minimal()
Я хочу оптимизировать # узлы и их расположение только для t переменная, следовательно, я подхожу:
mars.fit <- earth::earth(y ~ sex + t, data = df1, linpreds = "sex", thresh = 0, nfold = 10, ncross = 30, varmod.method = "lm")
Где используются последние три аргумента, чтобы я мог построить модель, используя: plotmo::plotmo(mars.fit, pt.col=1, level=.95): 
Коэффициенты, которые я получаю:
> summary(lm(y ~ mars.fit$bx, data = df1))$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.20434144 0.04734724 -4.315805 7.504512e-05
mars.fit$bxh(t-3.32) -0.27889610 0.06733161 -4.142127 1.323844e-04
mars.fit$bxh(3.32-t) 0.06078628 0.01666151 3.648306 6.297134e-04
mars.fit$bxsexF 0.15757143 0.02252191 6.996362 6.080110e-09
mars.fit$bxh(t-4.58) 0.15469776 0.04939908 3.131592 2.903537e-03
mars.fit$bxh(t-2) 0.12419721 0.04345528 2.858046 6.198530e-03
И что я хотел бы сделать, это извлечь из этого matrix эффекты t (уклоны) в каждом из 4 сегментов кусочной посадки. Я немного запутался в том, как интерпретировать это резюме matrix чтобы получить эту информацию.
Любая идея?
К тому же, mars.fit$dirs дает:
> mars.fit$dirs
sexF t
(Intercept) 0 0
h(t-3.32) 0 1
h(3.32-t) 0 -1
sexF 2 0
h(t-4.58) 0 1
h(t-2) 0 1
Любая идея, какова ценность 2 за sexF средства?
Мой другой вопрос, как бы я смоделировать interaction между переменной, для которой я хочу оптимизировать узлы, и переменной, которую я не оптимизирую (т. е. для линейного предиктора linpred)?
Вот пример набора данных, который включает в себя данные выше и другой набор из 56 образцов, с тем же дизайном, но из другой группы (группа B):
df2 <- data.frame(y = c(-0.07, -0.05, -0.08, -0.12, -0.02, -0.11, -0.16, 0.02, -0.24, -0.13, -0.2, -0.52, -0.2, -0.37, -0.12, 0.14, -0.02, -0.19, -0.05, -0.09, -0.13, -0.09, -0.44, -0.27, -0.16, -0.38, -0.23, -0.26, -0.09, -0.11, -0.06, -0.08, -0.02, -0.2, -0.21, -0.28, -0.5, -0.34, -0.16, -0.43, -0.16, -0.49, 0.04, 0.09, -0.01, -0.18, 0.16, -0.39, -0.21, -0.04, -0.09, -0.07, -0.13, -0.54, -0.13, -0.57, 0.27, 0.61, 1.57, 1.94, 2.2, 2.4, 2.13, 2.18, 2.12, 1.95, 1.99, 2.16, 2.38, 2.7, 0.4, 0.58, 1.33, 2.04, 1.99, 1.94, 1.59, 1.55, 1.94, 1.67, 2.36, 2.05, 2.25, 2.22, 0.31, 0.75, 1.84, 2.21, 2.51, 2.15, 2.33, 1.91, 2.4, 2.24, 1.02, 2.31, 2.45, 2.79, 0.35, 0.53, 1.32, 2.07, 2.13, 2.22, 2, 2.14, 2.28, 2.11, 1.92, 2.26, 2.1, 2.67),
t = rep(c(-1, 0, 1, 1.58, 2, 2.58, 3, 3.32, 3.58, 4.17, 4.58, 5.58, 6.17, 7.39),8),
group = c(rep("B",56),rep("A",56)),
sex = rep(c(rep("F",28),rep("M",28)),2))
df2$group <- factor(df2$group, levels = c("A","B"))
df2$sex <- factor(df2$sex, levels = c("M","F"))
df2$group_sex <- paste0(df2$group,"_",df2$sex)
Который выглядит как:
library(ggplot2)
ggplot(df2, aes(y = y, x = t, color = group_sex)) + geom_point(size = 3) + theme_minimal()
(group * t взаимодействие здесь довольно понятно)
Это правильно MARS модель:
earth::earth(y ~ sex + group*t, data = df2, thresh = 0, nfold = 10, ncross = 30, varmod.method = "lm")
Составление графика с использованием plotmo дает: 
И резюме этой модели:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.8761333 0.16813559 5.210873 9.704648e-07
mars.fit$bxgroupB -0.5412963 0.09186638 -5.892213 4.844315e-08
mars.fit$bxh(1.58-t) -0.1815588 0.07860784 -2.309678 2.290128e-02
mars.fit$bxh(groupB:t-2) 1.0255966 0.10599440 9.675950 3.890338e-16
mars.fit$bxh(groupB:t-4.58) -0.4184100 0.09712737 -4.307849 3.774816e-05
mars.fit$bxh(t-4.58) 0.4021478 0.06937728 5.796534 7.457303e-08
mars.fit$bxh(t-2) -0.8350190 0.11881038 -7.028165 2.338724e-10
mars.fit$bxh(t-0) 0.6909900 0.09677237 7.140364 1.356748e-10
mars.fit$bxh(groupB:t-0) -0.9231932 0.06616621 -13.952637 1.769394e-25
Здесь то, что я хотел бы извлечь, это group*t наклон на каждом отрезке (т. е. groupB:t-knot).
Спасибо