Jdan - 1 year ago 141
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.

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

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

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 = "")
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download