Логистическая регрессия 4-го графа поверхности (x, y, z, цвет)

Question

Я пытаюсь построить четырехмерный поверхностный график биномиальных отношений, включающий 4 количественные переменные (3 предикторные переменные и 1 ответ).У меня есть лесной покров (cov), ландшафтная неоднородность (het) и обилие хищников (abund), чтобы объяснить вероятность хищничества (pred).Поэтому я хотел бы поставить на график x = cov, y = het, z = abund и градиент цвета для представления pred.Я бы хотел график, похожий на этот:

Моя модель:

model <- glm(cbind(pred,npred)~ cov + het + abund + abund:cov + cov:het, data=results, family=binomial)

Мои сценарии для создания графика:

het <- seq(0, 2, 0.05288)
cov <- seq(0, 0.8, 0.021600)
abun <- seq(0, 0.9, 0.024)

gg <- expand.grid(het=het, cov=cov, abund=abund)
gg$pred < -predict(model,newdata=gg, type = "response")
head(gg)

with(gg, plot3d(het, cov, abund, type = "n"))
cols <- heat.colors(20)
cuts <- with(gg, cut(pred, breaks = 20))

with(gg, plot3d(het, cov, abund, type = "n"))
with(gg, surface3d(1:20,1:20, matrix(pred, ncol = 20),
                  color = cols[cuts], back = "fill"))

К сожалению, я думаю, что что-то не так и не работает, как я себе представлял.Может ли кто-нибудь помочь мне с сценарием или предложить что-то лучшее?

Мой набор данных:

het cov abun    pred    npred
1.90164 0.636460117 0.730424235 24  14
2.13681 0.409547496 0.735686442 8   42
2.18218 0.631526938 0.744585387 2   48
2.10223 0.568536934 0.746072471 8   36
2.11588 0.390850929 0.747720518 2   50
2.10288 0.465843118 0.754245798 6   42
2.28891 0.272948757 0.757754057 4   50
2.31959 0.281255571 0.776337246 2   32
2.10626 0.426842961 0.778512783 4   46
2.0163  0.546051041 0.780669728 6   50
1.91282 0.433635614 0.809807372 2   48
1.67011 0.594238535 0.827360862 18  26
1.80076 0.719426666 0.830916931 4   26
2.29307 0.471266098 0.839754622 2   44
2.00811 0.574829217 0.839935299 2   40
2.25667 0.416338714 0.849465507 0   54
2.02098 0.306265093 0.860655587 4   42
2.10236 0.263078082 0.860656944 6   36
1.23721 0.236391255 0.879960428 6   44

Answer 1

Ваши прогнозы являются функцией трех переменных-предикторов, поэтому вы не можете представлять их как одну поверхность. Есть два распространенных способа сделать это:

1. Построение 3d контуров с помощью функции misc3d::contour3d.

2. Построить несколько поверхностей для разных значений одной из переменных.

Перед первым шагом необходимо исправить некоторые опечатки в вашем коде: иногда вы используете abun, иногда abund: я всегда буду использовать abun. У вас также есть

gg$pred < -predict(model,newdata=gg, type = "response")

и пробел в середине <- полностью меняет значение. Наконец, использование dput(results) облегчает использование данных другими. Вот все, что вычистили:

results <- structure(list(het = c(1.90164, 2.13681, 2.18218, 2.10223, 2.11588, 
2.10288, 2.28891, 2.31959, 2.10626, 2.0163, 1.91282, 1.67011, 
1.80076, 2.29307, 2.00811, 2.25667, 2.02098, 2.10236, 1.23721
), cov = c(0.636460117, 0.409547496, 0.631526938, 0.568536934, 
0.390850929, 0.465843118, 0.272948757, 0.281255571, 0.426842961, 
0.546051041, 0.433635614, 0.594238535, 0.719426666, 0.471266098, 
0.574829217, 0.416338714, 0.306265093, 0.263078082, 0.236391255
), abun = c(0.730424235, 0.735686442, 0.744585387, 0.746072471, 
0.747720518, 0.754245798, 0.757754057, 0.776337246, 0.778512783, 
0.780669728, 0.809807372, 0.827360862, 0.830916931, 0.839754622, 
0.839935299, 0.849465507, 0.860655587, 0.860656944, 0.879960428
), pred = c(24L, 8L, 2L, 8L, 2L, 6L, 4L, 2L, 4L, 6L, 2L, 18L, 
4L, 2L, 2L, 0L, 4L, 6L, 6L), npred = c(14L, 42L, 48L, 36L, 50L, 
42L, 50L, 32L, 46L, 50L, 48L, 26L, 26L, 44L, 40L, 54L, 42L, 36L, 
44L)), class = "data.frame", row.names = c(NA, -19L))

model <- glm(cbind(pred,npred)~ cov + het + abun + abun:cov + cov:het,
             data=results, family=binomial)

het <- seq(0, 2, 0.05288)
cov <- seq(0, 0.8, 0.021600)
abun <- seq(0, 0.9, 0.024)
gg <- expand.grid(het=het, cov=cov, abun=abun)
gg$pred <- predict(model,newdata=gg, type = "response")
head(gg)
#       het cov abun         pred
# 1 0.00000   0    0 2.220446e-16
# 2 0.05288   0    0 2.220446e-16
# 3 0.10576   0    0 2.220446e-16
# 4 0.15864   0    0 2.220446e-16
# 5 0.21152   0    0 2.220446e-16
# 6 0.26440   0    0 2.220446e-16

Чтобы сделать контурный график, вам нужно выбрать уровни контура. Вы можете увидеть диапазон прогнозов, используя

range(gg$pred)
# [1] 2.220446e-16 1.000000e+00

поэтому ваша модель предсказывает вероятности, по существу, от 0 до 1. Я бы хотел, чтобы 20 контуров были между 0 и 1, но функция contour3d не выходит за пределы допустимого диапазона. контуры, поэтому я буду использовать 1/20, ..., 19/20:

levels <- (1:19)/20
colors <- heat.colors(19)

Функция contour3d хочет значения в массиве, а не в векторе:

pred <- array(gg$pred, c(length(het), length(cov), length(abun)))

Теперь нарисуем несколько осей, затем сюжет:

plot3d(gg, type = "n")
contour3d(pred, levels, het, cov, abun, color = colors, alpha = 0.5, add = TRUE)

Значение alpha = 0.2 делает поверхности частично прозрачными; Вы можете попробовать другие значения от 0 до 1, чтобы увидеть, что вам нравится. Вот что я получил. В R вы можете повернуть график, чтобы увидеть его под разными углами, или позвонить по номеру rglwidget(), чтобы сохранить его во вращающейся форме, которую можно просмотреть в веб-браузере.

Я оставлю это вам, чтобы решить, как сделать второй вид заговора.