- Other 3D Viewing Options
- Some Calibration Experiments
- Interactive Explorations
- Example: Ethanol Data
- Example: Earth Quake Locations and Magnitudes
- Visualizing 2D and 3D Density Estimates
- Example: Environmental Data
- Medical Imaging Data
- Notes
- Two Classification Problems
- Grand and Guided Tours
- Scagnostics
- Two Sliders Using Shiny

Other options for viewing a 3D scene include:

- Anaglyph 3D using red/cyan glasses.
- Polarized 3D as currently used in many 3D movies.
- Stereograms, stereoscopy

A stereogram

- presents each eye with an image from a slightly different viewing angle;
- the brain fuses the images into a 3D view in a process called
*stereopsis*.

This approach is or was used

- with virtual reality displays like Oculus Rift;
- in the View-Master toy;
- in aerial photography using viewers like this:

You can create your own viewer with paper or cardboard tubes, for example.

A stereo image for the soil data surface:

```
tp <- trellis.par.get()
trellis.par.set(theme=col.whitebg())
ocol <- trellis.par.get("axis.line")$col
oclip <- trellis.par.get("clip")$panel
trellis.par.set(list(axis.line = list(col = "transparent"),
clip = list(panel = FALSE)))
print(wireframe(soi.fit ~ soi.grid$easting * soi.grid$northing,
cuts = 9,
aspect=diff(range(soi.grid$n))/diff(range(soi.grid$e)),
xlab = "Easting (km)",
ylab = "Northing (km)", shade=T,
screen = list(z = 40, x = -70, y = 3)),
split = c(1,1,2,1), more = TRUE)
print(wireframe(soi.fit ~ soi.grid$easting * soi.grid$northing,
cuts = 9,
aspect=diff(range(soi.grid$n))/diff(range(soi.grid$e)),
xlab = "Easting (km)",
ylab = "Northing (km)", shade=T,
screen = list(z = 40, x = -70, y = 0)),
split = c(2,1,2,1))
```

```
trellis.par.set(list(axis.line = list(col = ocol), clip = list(panel = oclip)))
trellis.par.set(tp)
```

Linear structure with normally distributed data:

```
n <- 1000
d1 <- data.frame(x = rnorm(n), y = rnorm(n))
d1 <- mutate(d1, z = x + y)
```

`splom(d1)`

```
xyplot(z ~ x | equal.count(y), data = d1)
```

`xyplot(z ~ x | equal.count(y, 12), data = d1)`

`xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d1)`

`splom(d1, group = abs(x) < 0.1)`

```
clear3d()
points3d(d1)
```

You must enable Javascript to view this page properly.

Linear structure with uniformly distributed data:

```
d2 <- data.frame(x = runif(n, -1, 1), y = runif(n, -1, 1))
d2 <- mutate(d2, z = x + y)
```

`splom(d2)`

`xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d2)`

`splom(d2, group = abs(x) < 0.1)`

```
clear3d()
points3d(d2)
```

You must enable Javascript to view this page properly.

Normal data with non-linearity and noise:

`d3 <- mutate(d1, z = cos(x) + cos(y) + 0.5 * (x + y) + rnorm(x, 0, .2))`

`splom(d3)`

`xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d3)`

`splom(d3, group = abs(x) < 0.1)`

`splom(d3, group = abs(x + 1) < 0.1)`

Uniform data with non-linearity and noise:

`d4 <- mutate(d2, z = cos(x) + cos(y) + 0.5 * (x + y) + rnorm(x, 0, .1))`

`xyplot(z ~ x | equal.count(y, 12, overlap = -0.8), data = d4)`

`splom(d4, group = abs(x) < 0.1)`

Full data in the background for context:

```
xyplot(z ~ x | equal.count(y, 12, overlap = -0.8), data = d4,
panel = function(x, y, ...) {
panel.xyplot(d4$x, d4$z, col = "grey")
panel.xyplot(x, y,...)
})
```

```
clear3d()
points3d(d4)
```

You must enable Javascript to view this page properly.

Two-dimensional views, conditioning on two variables:

```
n <- 10000
d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
d <- mutate(d,
x4 = as.numeric(scale(x1 + cos(x2) + sin(x3) + rnorm(n, 0, 0.1))))
```

`splom(d)`

`splom(~d, data = d, group = abs(x1) < 0.1)`

`splom(~d, data = d, group = abs(x1) < 0.1 & abs(x2) < 0.1)`

```
xyplot(x4 ~ x3 |
equal.count(x1, 4, overlap = -0.8) +
equal.count(x2, 4, overlap = -0.8),
data = d, asp = 1)
```

Three-dimensional views, conditioning on one variable:

`dd <- filter(d, abs(x1) < 0.1)[2:4]`

```
clear3d()
points3d(dd)
```

You must enable Javascript to view this page properly.

```
clear3d()
spheres3d(dd, radius = 0.05)
```

You must enable Javascript to view this page properly.

```
clear3d()
spheres3d(dd, radius = 0.05, col = ifelse(abs(dd$x2) < 0.1, "red", "blue"))
```

You must enable Javascript to view this page properly.

```
dd1 <- filter(d, abs(x1 - 1) < 0.1)[2:4]
dd2 <- filter(d, abs(x1 + 1) < 0.1)[2:4]
clear3d()
spheres3d(dd1, radius = 0.05, color = rep("red", ncol(dd1)))
spheres3d(dd2, radius = 0.05, color = rep("blue", ncol(dd2)))
```

You must enable Javascript to view this page properly.

```
ddx <- mutate(d,
v = ifelse(abs(x1 - 1) < 0.1,
1,
ifelse(abs(x1 + 1) < 0.1, 2, NA)))
ddx <- filter(ddx, ! is.na(v))
cloud(x4 ~ x2 * x3, group = v, data = ddx, auto.key = TRUE)
```

`cloud(x4 ~ x2 * x3 | v, data = ddx)`

`tkrplot`

A simple animation function using the tkrplot` package:

```
library(tkrplot)
## Loading required package: tcltk
tkranim <- function(fun, v, lo, hi, res = 0.1, title, label) {
tt <- tktoplevel()
if (! missing(title))
tkwm.title(tt, title)
draw <- function() fun(v)
img <-tkrplot(tt, draw)
tv <- tclVar(init = v)
f<-function(...) {
vv <- as.numeric(tclvalue(tv))
if (vv != v) {
v <<- vv
tkrreplot(img)
}
}
s <- tkscale(tt, command=f, from=lo, to=hi, variable=tv,
showvalue=FALSE, resolution=res, orient="horiz")
if (! missing(label))
tkpack(img, tklabel(tt, text=label), s)
else tkpack(img,s)
}
```

Interactive conditioning for one of the calibration examples:

```
fd4 <- function(v)
print(splom(~ d4, group = abs(x - v) < 0.1, data = d4))
tkranim(fd4, 0, -1, 1, 0.05)
```

Interactive slicing for the soil resistivity surface:

```
f <- function(i) {
e <- eastseq[i]
idx <- which(sf$easting == e)
r <- range(sf$fit)
print(xyplot(fit ~ northing, data = sf[idx,],
ylim = r, type = "l", main = sprintf("Easting = %f", e)))
}
tkranim(f, 1, 1, 85, 1)
```

Interactive conditioning for the soil resistivity data:

```
binwidth <- 0.1 ## bad idea to use a global variable
g <- function(i) {
e <- eastseq[i]
idx <- which(abs(soil$easting - e) <= binwidth)
ry <- range(soil$resistivity)
rx <- range(soil$northing)
print(xyplot(resistivity ~ northing, data = soil[idx,],
xlim = rx, ylim = ry, main = sprintf("Easting = %f", e)))
}
tkranim(g, 1, 1, 85, 1)
```

A version of the animation function using the `shiny`

package:

```
shinyAnim <- function(fun, v, lo, hi, step = NULL, label = "Frame") {
require(shiny)
shinyApp(
ui = fluidPage(
sidebarLayout(
sidebarPanel(sliderInput("n", label, lo, hi, v, step,
animate = animationOptions(100))),
mainPanel(plotOutput("plot"))
)
),
server = function(input, output, session) {
session$onSessionEnded(function() shiny::stopApp())
output$plot <- renderPlot(fun(input$n))
}
)
}
```

`shinyAnim(fd4, 0, -1, 1, 0.05)`

`shinyAnim(f, 1, 1, 85)`

`shinyAnim(g, 1, 1, 85, 1)`

To show how a plot changes with a single parameter it is possible to pre-compute the images and assemble them into a movie.

This is easy to do in Rmarkdown, but requires the `ffmpeg`

software to be installed.

Currently, `ffmpeg`

is not available on the CLAS Linux systems.

A simple example (Rmd).

Linked plots and *brushing* were pioneered at Bell Laboratories in the early 1980s.

Some videos are available in the ASA Statistics Computing and Statistical Graphics Sections’ video library.

Trellis displays were developed by the same research group

*after*the development of brushing.

The `iplots`

package provides a set of basic plot that are all linked by default.

```
library(iplots)
with(soil, iplot(northing, resistivity))
```

`with(soil, iplot(easting,northing))`

The `rggobi`

package provides another framework for linked plots and brusing within R.

```
library(rggobi)
ggobi(soil)
```

Several other tools with linking/brushing support and good integration with R are under development.

- One of the furthest along seems to be
`rbokeh`

. `plotly`

also provides some support for linked brushing.

The `ethanol`

data frame in the `lattice`

package contains data from an experiment efficiency and emissions in small one-cylinder engines.

The data frame contains 88 observations on three variables:

`NOx`

: Concentration of nitrogen oxides (NO and NO2) in micrograms.`C`

Compression ratio of the engine.`E`

Equivalence ratio, a measure of the richness of the air and ethanol fuel mixture.

A scatterplot matrix:

```
library(lattice)
splom(ethanol)
```

Conditional plots of `NOx`

against `E`

for each level of `C`

:

`xyplot(NOx ~ E | C, data = ethanol)`

```
ggplot(ethanol, aes(x = E, y = NOx)) + geom_point() + facet_wrap(~C)
```

A single plot with color to separate levels of `C`

:

`xyplot(NOx ~ E, group = C, data = ethanol, auto.key = TRUE)`

```
ggplot(ethanol, aes(x = E, y = NOx, color = factor(C))) + geom_point()
```

Augmented with smooths:

```
xyplot(NOx ~ E, group = C, type = c("p", "smooth"), data = ethanol,
auto.key = TRUE)
```

```
ggplot(ethanol, aes(x = E, y = NOx, color = factor(C))) +
geom_point() + geom_smooth(se = FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
```

Conditional plots of `NOx`

against `C`

given levels of `E`

:

```
with(ethanol, {
Equivalence.Ratio <- equal.count(E, number = 9, overlap = 0.25)
xyplot(NOx ~ C | Equivalence.Ratio,
panel = function(x, y) {
panel.xyplot(x, y)
panel.loess(x, y, span = 1)
},
aspect = 2.5,
layout = c(5, 2),
xlab = "Compression Ratio",
ylab = "NOx (micrograms/J)")
})
```

An interactive 3D view:

```
library(rgl)
clear3d()
with(ethanol,points3d(E, C / 5, NOx / 5,
size = 4, col = rep("black", nrow(ethanol))))
```

You must enable Javascript to view this page properly.

The `quakes`

data frame contains data on locations of seismic events of magnitude 4.0 or larger in a region near Fiji.

The time frame is from 1964 to perhaps 2000. More recent data is available from a number of sources on the web.

Some scatter plot matrices:

```
library(lattice)
splom(quakes)
```

`splom(quakes, cex = 0.1)`

`splom(quakes, alpha = 0.1)`

```
splom(quakes, panel = panel.smoothScatter)
## (loaded the KernSmooth namespace)
```

Where the quakes are located:

```
library(maps)
map("world2", c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = "red", cex = 0.1))
```

```
map("world2", c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = heat.colors(nrow(quakes))[rank(depth)],
cex = 0.1))
```

```
map("world2",c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = heat.colors(nrow(quakes))[rank(depth)]))
```

Showing `depth`

with grouping:

```
quakes2 <- mutate(quakes, depth_cut = cut_number(depth, 3))
xyplot(lat ~ long, aspect = "iso", groups = depth_cut,
data = quakes2, auto.key = list(columns = 3, title = "Depth"))
```

```
ggplot(quakes2, aes(x = long, y = lat, color = depth_cut) ) +
geom_point() + theme_bw() + coord_map()
```

```
ggplot(quakes2, aes(x = long, y = lat, colour = depth_cut, size = mag)) +
geom_point() + theme_bw() + coord_map()
```

Showing `depth`

with a trellis display:

`xyplot(lat ~ long | depth_cut, aspect = "iso", data = quakes2)`

```
xyplot(lat ~ long | depth_cut, aspect = "iso", type=c("p","g"),
data = quakes2, auto.key = TRUE)
```

```
ggplot(quakes2, aes(x = long, y = lat)) +
geom_point(color = "blue") + facet_wrap(~depth_cut, nrow = 1) +
theme_bw() + coord_map()
```

Adding the full data for background context:

```
xyplot(lat ~ long | cut(depth, 3), aspect = "iso", type=c("p","g"),
panel = function(x, y, ...) {
panel.xyplot(quakes$long, quakes$lat, col = "grey")
panel.xyplot(x, y, ...)
},
data = quakes)
```

```
## quakes does not contain the depth_cut variable used in the facet
ggplot(quakes2, aes(x = long, y = lat)) +
geom_point(data = quakes, color = "gray") +
geom_point(color = "blue") + facet_wrap(~depth_cut, nrow = 1) +
theme_bw() + coord_map()
```

3D views of `lat`

, `long`

, and `depth`

:

`cloud(-depth ~ long * lat, data = quakes)`

```
library(rgl)
bg3d(col = "black")
clear3d()
with(quakes, points3d(long, lat, -depth / 50, size = 2,
col = rep("white", nrow(quakes))))
```

You must enable Javascript to view this page properly.

```
clear3d()
with(quakes, points3d(long, lat, -depth / 50, size = 2,
col=heat.colors(nrow(quakes))[rank(mag)]))
```

You must enable Javascript to view this page properly.

We have already seen some examples of plotting 2D density contours on a scatter plot using `geom_density2d`

.

Here is a version that explicitly creates the density estimate:

```
data(geyser, package = "MASS")
geyser <- mutate(geyser, duration = lag(duration))
geyser <- filter(geyser, ! is.na(duration))
plot(waiting ~ duration, data = geyser)
```

```
sc <- function(x) as.numeric(scale(x))
de <- with(geyser, MASS::kde2d(sc(duration), sc(waiting), n = 50, h = 0.5))
contour(de)
with(geyser, points(sc(duration), sc(waiting)))
```

The density estimate can also be viewed with a perspective plot or an interactive 3D plot:

`persp(de)`

```
clear3d()
bg3d(col = "white")
surface3d(unique(de$x), unique(de$y), de$z, col = "red")
box3d()
axes3d()
```

You must enable Javascript to view this page properly.

A 3D normal mixture:

```
library(rgl)
library(misc3d)
x <- matrix(rnorm(9000), ncol = 3)
x[1001 : 2000, 2] <- x[1001 : 2000, 2] + 4
x[2001 : 3000, 3] <- x[2001 : 3000, 3] + 4
```

```
clear3d()
bg3d(col = "white")
points3d(x = x[,1], y = x[,2], z = x[,3], size = 2, color = "black")
```