Jdan Jdan - 9 months ago 73
R Question

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

Original title (vague): How to make circular surface from x and z value only




I have data that relate to an x-axis and a z-axis similar to values of
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")


I want generate a surface / matrix where values of
new.data$fit
represent the z-axis and x- and y-axis are generated from the radius being the
new.data$Dist
.

In other words I want a circle generated from radius
Dist
and cells populated by
z
value of logistic probability curve. I would like to say that I have tried some certain solutions but not even sure where to begin.

Answer Source

So you want to plot surface of revolution, by spinning the logistic curve around vertical line Dist = 0. Statistically I don't know why we need this, but purely from a mathematical aspect and for the sake of 3D visualization, this is sort of useful, hence I decide to answer this.

All we need is a function of the initial 2D curve f(d), where d is distance from a point to spinning centre, and f is some smooth function. As we will use outer to make surface matrix, f must be defined to be vectorized function in R. Now the surface of revolution is generated as f3d(x, y) = f((x ^ 2 + y ^ 2) ^ 0.5).

In your logistic regression setting, the above f is logistic curve, the predicted response of a GLM. It can be obtained from predict.glm, which is a vectorized function. The following code fits a model, and define such function f plus its 3D extension.

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)
  }

Due to symmetry, we only call f3d on the 1st quadrant for a surface matrix X1, while flipping X1 for surface matrices on other quadrants.

## 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]

Finally we combine matrices from different quadrants and make 3D plot. Note the combination order is quadrant 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")

enter image description here


Making a toy routine for plotting surface of revolution

In this part, we define a toy routine for plotting surface of revolution. As noted above, we need for this routine:

  1. a (vectorized) function of 2D curve: f;
  2. evaluation grid on the 1st quadrant {(x, y) | x >= 0, y >= 0} (due to symmetry we take y = x);
  3. possible additional arguments to f, and customized graphical parameters to persp.

The following is a simple implementation:

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))
  }

Now suppose we want to spin a cosine wave on [0, pi] for a surface of revolution, we can do

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

enter image description here

We can also use surfrev to plot your desired logistic curve:

## `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 = "")