Я пытаюсь смоделировать логистический ответ c на основе 4 непрерывных переменных из эксперимента, который я провел. Первоначально я использовал множественную регрессию и достиг довольно хорошего результата, но недавно мне предложили вместо этого использовать GAM. Я немного растерялся из-за того, как правильно сделать выбор модели для GAM, а также как интерпретировать некоторые предупреждения, которые я получаю из своей множественной регрессии GLM. Я подозреваю, что мои проблемы связаны с переоснащением, но я не знаю, как их обойти.
В основном вопрос заключается в следующем: каков наилучший / самый экономный способ моделирования этих данных?
df:
df <- structure(list(response = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0),
V1 = c(14.2,13.67, 13.05, 14.18, 13.4, 14.12, 14.22, 14.15, 13.35, 13.67,
18.58, 18.27, 18.6, 17.94, 18.38, 18.98, 18.15, 19, 18.55, 18.53,
20.77, 21.65, 21.03, 21.57, 21.25, 21.63, 21.6, 21.09, 21.62,
21.6, 26.23, 26.52, 25.7, 26.57, 26.6, 26.25, 26.48, 26.26, 26.25,
26.4, 28.98, 29.45, 29.2, 29.65, 29.38, 28.6, 28.42, 28.95, 28.85,
28.8), V2 = c(27.2, 37.98, 24.63, 32.97, 30.27, 18.66, 13.77,
33.99, 15.8, 21.32, 14.21, 15.81, 35.83, 21.64, 26.93, 38.62,
34.03, 18.76, 24.12, 29.67, 29.83, 33.22, 27.11, 24.92, 21.72,
39.02, 12.93, 18.44, 36.34, 15.81, 13.29, 21.04, 19.05, 33.62,
30.52, 16.07, 28.43, 24.97, 39.9, 37.05, 19.31, 31.3, 34.08,
13.63, 25.1, 28.93, 22.36, 34.69, 39.5, 16.41),
V3 = c(8.06, 7.87, 7.81, 7.72, 8.04, 7.66, 7.72, 7.87, 7.72, 7.98, 7.59, 7.9,
8.08, 7.64, 8.02, 7.73, 7.77, 7.74, 7.66, 7.71, 8.05, 7.68, 7.63,
7.7, 7.64, 7.8, 7.7, 7.98, 7.86, 7.68, 7.65, 7.74, 7.99, 7.75,
7.91, 7.64, 7.69, 7.78, 7.69, 7.66, 7.72, 7.76, 7.71, 7.88, 7.63,
7.7, 7.99, 7.82, 7.75, 7.93),
V4 = c(362.12, 645.38, 667.54,
957.44, 391.84, 818.34, 732.91, 649.05, 722.02, 406.71, 918.9,
471.32, 363.77, 926.82, 385.4, 1038.91, 850.67, 715.11, 964.79,
890.11, 370.51, 1078.68, 1083.7, 893.76, 1026.1, 887.29, 737.68,
406.76, 690.39, 872.8, 847.26, 738.53, 397.33, 895.3, 563.93,
991.17, 957.28, 734.55, 1140.5, 1199.12, 817.17, 800.5, 992.82,
533.45, 1123.29, 943.25, 411.59, 745.22, 929.75, 460.82)),
row.names = c(NA,-50L), class = "data.frame")
Следует отметить, что при проведении экспериментов и зная о системе, я знаю, что V1 и V2 оказывают наибольшее влияние на реакцию. Вы также можете увидеть это, построив ответ только по этим переменным, поскольку все положительные ответы сгруппированы в этом двумерном пространстве. Кроме того, если взглянуть на некоторые временные сплайны, кажется, что V1 линейно связан с откликом, V2 - квадратично, V3 - совсем нет, а V4 - слабо квадратично c.
Еще одно важное замечание: V3 и V4 - это две разные меры одного и того же, поэтому они сильно коррелированы и не будут используется в любых моделях вместе.
Итак, сначала я попытался смоделировать все это в регрессии mutliple logisti c: мне посоветовали протестировать целую кучу разных моделей в моем выборе модели, поэтому я записал их в список и запустил их все в al oop:
formulas <- list(# single predictors
response ~ V1,
response ~ V2,
response ~ V3,
response ~ V4,
# two predictors
response ~ V1 + V2,
response ~ V1 + V3,
response ~ V1 + V4,
response ~ V2 + V3,
response ~ V2 + V4,
# three predictors
response ~ V1 + V2 + V3,
response ~ V1 + V2 + V4,
# start quadratic models
response ~ V2 + I(V2^2) + V1 + I(V1^2),
response ~ V2 + I(V2^2) + V1 + I(V1^2) + V3,
response ~ V2 + I(V2^2) + V1 + I(V1^2) + V4,
response ~ V1 + V2 + I(V1^2),
response ~ V1 + V2 + I(V1^2) + V3,
response ~ V1 + V2 + I(V1^2) + V4,
response ~ V1 + I(V1^2),
response ~ V1 + V2 + I(V2^2),
response ~ V1 + V2 + I(V2^2) + V3,
response ~ V1 + V2 + I(V2^2) + V4,
response ~ V2 + I(V2^2),
response ~ V2 + I(V2^2) + V1 + I(V1^2),
# add interactions
response ~ V1 + V2 + V1*V2,
response ~ V1 + V2 + V1*V2 + V3,
response ~ V1 + V2 + V1*V2 + V4,
# quadratic with interaction
response ~ V1 + V2 + V1*V2 + V3 + I(V1^2),
response ~ V1 + V2 + V1*V2 + V3 + I(V2^2),
response ~ V1 + V2 + V1*V2 + V4 + I(V1^2),
response ~ V1 + V2 + V1*V2 + V4 + I(V2^2)
)
# run them all in a loop, then order by AIC
selection <- purrr::map_df(formulas, ~{
mod <- glm(.x, data= df, family="binomial")
data.frame(formula = format(.x),
AIC = round(AIC(mod),2),
BIC = round(BIC(mod),2),
R2adj = round(DescTools::PseudoR2(mod,which=c("McFaddenAdj")),3)
)
})
warnings()
warnings()
# this returns a bunch of warnings about coercing the formulas into vectors, ignore those.
# however, this also lists the following for a handful of the models:
# "glm.fit: fitted probabilities numerically 0 or 1 occurred"
# which means perfect separation, but I'm not sure if this is a totally bad thing
# or not, as perfect separation actually pretty much does exist in the real data
# then we arrange by AIC and get our winning model:
selection %>% arrange(desc(AIC))
Таким образом, используя эту технику, мы находим две лучшие модели response ~ V1 + V2 + I(V2^2) + V4 и response ~ V1 + V2 + I(V2^2). Но при запуске их по одному мы получаем ошибку numerically 1 or 0 для обоих, и мы видим, что единственное различие между ними (добавленный V4) не является статистически значимым само по себе в лучшей модели. Итак ... какие мы используем ??
bestmod1 <- glm(response ~ V1 + V2 + I(V2^2) + V4,
family="binomial",
data=df)
summary(bestmod1)$coefficients
bestmod2 <- glm(response ~ V1 + V2 + I(V2^2),
family="binomial",
data=df)
summary(bestmod2)$coefficients
Метод 2: GAMs
Подобная техника здесь, перечислить все формулы и запустить в al oop.
library(mgcv)
gam_formulas <- list( # ind. main effects,
response ~ s(V1),
response ~ s(V2),
response ~ s(V3),
response ~ s(V4),
# two variables
response ~ s(V1) + s(V2),
response ~ s(V1) + s(V3),
response ~ s(V1) + s(V4),
response ~ s(V2) + s(V3),
response ~ s(V2) + s(V4),
# three variables
response ~ s(V1) + s(V2) + s(V3),
response ~ s(V1) + s(V2) + s(V4),
# add interactions
response ~ te(V1, V2),
response ~ te(V1, V2) + s(V3),
response ~ te(V1, V2) + s(V4),
response ~ te(V1, V3),
response ~ te(V1, V3) + s(V2),
response ~ te(V1, V4),
response ~ te(V1, V4) + s(V2),
response ~ te(V2, V3),
response ~ te(V2, V3) + s(V1),
response ~ te(V2, V4),
response ~ te(V2, V4) + s(V1),
response ~ te(V2, by=V1),
response ~ te(V1, by=V2),
response ~ te(V2, by=V3),
# two interactions?
response ~ te(V1, V3) + te(V1, V2),
response ~ te(V1, V4) + te(V1, V2),
response ~ te(V2, V3) + te(V1, V2),
response ~ te(V2, V4) + te(V1, V2)
)
gam_selection <- purrr::map_df(gam_formulas, ~{
gam <- gam(.x,
data= df, # always use same df
family="binomial",
method="REML") # always use REML smoothing method
data.frame(cbind(formula = as.character(list(formula(.x))),
round(broom::glance(gam),2),
R2 = summary(gam)$r.sq
))
})
# similarly, this gives a bunch of warnings about coercing the formulas into characters,
# but also this, for a handful of the models, which I am guessing is an overfitting thing:
# In newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, ... :
# Fitting terminated with step failure - check results carefully
gam_selection %>% arrange(desc(AIC))
но это возвращает кучу дурацких вещей, так как многие модели (даже не обязательно с подобными формулами или значениями AI C) говорят, что R2 = 1,00, и они являются формулами, которые не это имеет смысл, биологически. Почему это происходит? Что мне с этим делать? (Я знаю, что это как-то связано с использованием "REML", потому что некоторые из этих ошибок go исчезли без этой строки). По моему мнению, наиболее точным является третий из лучших, согласно значениям AI C: response ~ te(V2, by = V1), с использованием V2 в качестве сглаженной переменной и V1 в качестве линейной.
Кроме того, если присмотреться к двум лучшим играм, согласно AI C, ни одна из переменных сама по себе не является существенной (p-значения = 1 ... странно), что заставляет меня чувствовать, что я не должен использовать это.
bestgam <- gam(response ~ s(V1) + s(V2) + s(V4),
data= df, # always use same df
family="binomial",
method="REML")
summary(bestgam)
bestgam2 <- gam(response ~ s(V1) + s(V2) + s(V3),
data= df, # always use same df
family="binomial",
method="REML")
summary(bestgam2)
bestgam3 <- gam(response ~ te(V2, by = V1),
data= df, # always use same df
family="binomial",
method="REML")
summary(bestgam3) # this is the one I think I should be using
По сути, я не знаю, почему я бы использовал GAM вместо GLM или наоборот, тогда как выбирать переменные и избегать переобучения в процессе. Любой совет приветствуется.
Спасибо!