为旋转曲面生成 3D 图（GLM 逻辑曲线示例）

Question

原标题（含糊）：如何仅根据x和z值制作圆形表面

我有与 x-axis 和 z-axis 相关的数据，类似于 new.data 的值：

mydata <- structure(list(Dist = c(82, 82, 85, 85, 126, 126, 126, 126, 178, 
178, 178, 178, 178, 236, 236, 236, 236, 236, 312, 368, 368, 368, 
368, 368, 425, 425, 425, 425, 425, 425, 560, 560, 560, 560, 560, 
612, 612, 612, 612), pDet = c(1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 
1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 1, 0, 0)), .Names = c("Dist", "pDet"), row.names = c(NA, 
-39L), class = "data.frame")

model <- glm(pDet ~ Dist, data = mydata, family = binomial(link = "logit"))
new.data <- data.frame(Dist = seq(0, 650, 50))
new.data$fit <- predict(model, newdata = new.data, type="response")

我想生成一个曲面/矩阵，其中 new.data$fit 的值表示 z-axis，x- 和 y-axis 是从 new.data$Dist 的半径生成的。

换句话说，我想要一个由半径 Dist 生成的圆和由逻辑概率曲线的 z 值填充的单元格。我想说我已经尝试了一些特定的解决方案，但甚至不确定从哪里开始。

Answer 1

因此，您想要通过围绕垂直线 Dist = 0 旋转逻辑曲线来绘制 surface of revolution。 统计上我不知道为什么我们需要这个，但纯粹从数学方面和为了 3D 可视化，这是有用的，因此我决定回答这个问题。

我们只需要初始二维曲线的函数 f(d)，其中 d 是点到旋转中心的距离，f 是一些平滑函数。由于我们将使用outer来制作曲面矩阵，因此必须将f定义为R中的矢量化函数。现在旋转曲面生成为f3d(x, y) = f((x ^ 2 + y ^ 2) ^ 0.5)。

在您的逻辑回归设置中，上面的 f 是逻辑曲线，即 GLM 的预测响应。它可以从 predict.glm 中获得，这是一个矢量化函数。下面的代码适合一个模型，并定义这样的函数 f 加上它的 3D 扩展。

mydata <- structure(list(Dist = c(82, 82, 85, 85, 126, 126, 126, 126, 178, 
178, 178, 178, 178, 236, 236, 236, 236, 236, 312, 368, 368, 368, 
368, 368, 425, 425, 425, 425, 425, 425, 560, 560, 560, 560, 560, 
612, 612, 612, 612), pDet = c(1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 
1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 1, 0, 0)), .Names = c("Dist", "pDet"), row.names = c(NA, 
-39L), class = "data.frame")

model <- glm(pDet ~ Dist, data = mydata, family = binomial(link = "logit"))

## original 2D curve
f <- function (d, glmObject) 
  unname(predict.glm(glmObject, newdata = list(Dist = d), type = "response"))

## 3d surface function on `(x, y)`
f3d <- function (x, y, glmObject) {
  d <- sqrt(x ^ 2 + y ^ 2)
  f(d, glmObject)
  }

由于对称性，我们只在第一象限上调用f3d一个表面矩阵X1，而翻转X1其他象限上的表面矩阵

## prediction on the 1st quadrant
x1 <- seq(0, 650, by = 50)
X1 <- outer(x1, x1, FUN = f3d, glmObject = model)

## prediction on the 2nd quadrant
X2 <- X1[nrow(X1):2, ]

## prediction on the 3rd quadrant
X3 <- X2[, ncol(X2):2]

## prediction on the 4th quadrant
X4 <- X1[, ncol(X1):2]

最后我们将来自不同象限的矩阵组合起来制作 3D 图。注意组合顺序是象限 3-4-2-1.

## combined grid
x <- c(-rev(x1), x1[-1])
# [1] -650 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100  -50    0   50
#[16]  100  150  200  250  300  350  400  450  500  550  600  650

## combined matrix
X <- cbind(rbind(X3, X4), rbind(X2, X1))

## make 3D surface plot
persp(x, x, X, col = "lightblue", theta = 35, phi = 40,
      xlab = "", ylab = "", zlab = "pDet")

制作绘制旋转曲面的玩具例程

在这一部分中，我们定义了一个绘制旋转曲面的玩具例程。如上所述，我们需要这个例程：

二维曲线的（矢量化）函数：f；
第一象限的评估网格{(x, y) | x >= 0, y >= 0}（由于对称性我们取y = x）；
f 的可能附加参数，以及 persp 的自定义图形参数。

下面是一个简单的实现：

surfrev <- function (f, x, args.f = list(), ...) {
  ## extend `f` to 3D
  .f3d <- function (x, y) do.call(f, c(list(sqrt(x ^ 2 + y ^ 2)), args.f))
  ## surface evaluation
  X1 <- outer(x, x, FUN = .f3d)
  X2 <- X1[nrow(X1):2, ]
  X3 <- X2[, ncol(X2):2]
  X4 <- X1[, ncol(X1):2]
  xbind <- c(-rev(x), x[-1])
  X <- cbind(rbind(X3, X4), rbind(X2, X1))
  ## surface plot
  persp(xbind, xbind, X, ...)
  ## invisible return
  invisible(list(grid = xbind, z = X))
  }

现在假设我们想在[0, pi]上旋转一个余弦波作为一个旋转面，我们可以这样做

surfrev(cos, seq(0, pi, by = 0.1 * pi), col = "lightblue", theta = 35, phi = 40,
        xlab = "", ylab = "", zlab = "")

我们还可以使用 surfrev 绘制您想要的逻辑曲线：

## `f` and `model` defined at the beginning
surfrev(f, seq(0, 650, by = 50), args.f = list(glmObject = quote(model)),
        col = "lightblue", theta = 35, phi = 40, xlab = "", ylab = "", zlab = "")

为旋转曲面生成 3D 图（GLM 逻辑曲线示例）

Producing 3D plot for surface of revolution (a GLM logistic curve example)

r

plot

3d

glm

surface