Previously, we described the essentials of R programming and provided quick start guides for importing data into **R** as well as visualizing data using R base graphs.

This chapter describes how to create static and interactive three-dimension (**3D**) **graphs**. We provide also an R package named graph3d to easily build and customize, step by step, 3D graphs in R software.

- Static 3D Scatter Plots
- Interactive 3D Graphs

- Install and load scaterplot3d
- Prepare the data
- The function scatterplot3d()
- Basic 3D scatter plots
- Change the main title and axis labels
- Change the shape and the color of points
- Change point shapes by groups
- Change point colors by groups
- Change the global appearance of the graph
- Remove the box around the plot
- Add grids on scatterplot3d

- Add bars
- Modification of scatterplot3d output
- Add legends
- Add point labels
- Add regression plane and supplementary points

Read more: —> Simple 3D Scatter Plots: scatterplot3d Package.

- Install and load plot3D package
- Prepare the data
- Scatter plots
- Basic scatter plot
- Change the type of the box around the plot
- Change the color by groups
- Change the position of the legend
- 3D viewing direction
- Titles and axis labels
- Tick marks and labels
- Add points and text to an existing plot

- Line plots
- Add confidence interval
- 3D fancy Scatter plot with small dots on basal plane
- Regression plane

- text3D: plot 3-dimensionnal texts
- text3D and scatter3D
- 3D Histogram
- scatter2D: 2D scatter plot
- text2D
- Interactive plot

Read more: —> Advanced 3D Graphs: plot3D Package.

- Install and load required packages
- Prepare the data
- The function scatter3d
- Basic 3D scatter plots
- Plot the points by groups
- Default plot
- Remove the surfaces
- Add concentration ellipsoids
- Change point colors by groups

- Axes
- Change axis labels
- Remove axis scales
- Change axis colors

- Add text labels for the points
- Export images

Read more: —> Interactive 3D Scatter Plots.

- Install the RGL package
- Load the RGL package
- Prepare the data
- Start and close RGL device
- 3D scatter plot
- Basic graph
- Change the background and point colors
- Change the shape of points

- rgl_init(): A custom function to initialize RGL device
- Add a bounding box decoration
- Add axis lines and labels
- Set the aspect ratios of the x, y and z axes
- Change the color of points by groups
- Change the shape of points
- Add an ellipse of concentration
- Regression plane
- Create a movie of RGL scene
- Export images as png or pdf
- Export the plot into an interactive HTML file
- Select a rectangle in an RGL scene
- Identify points in a plot
- R3D Interface
- RGL functions
- Device management
- Shape functions
- Scene management
- Setup the environment
- Appearance setup
- Export screenshot
- Assign focus to an RGL window

Read more —> Guide to RGL 3D Visualization System.

This analysis has been performed using **R statistical software** (ver. 3.2.4).

In my previous articles, I already described how to make **3D graphs** in R using the package below:

- scatterplot3d, non interactive
- scatter3d, interactive
- rgl, interactive

To close the discussion about 3D, in this tutorial I’ll describe the impressive **plot3D** package and its extension **plot3Drgl** package.

**plot3D**, from Karline Soetaert, is an R package containing many functions for 2D and 3D plotting: *scatter3D*, *points3D*, *lines3D*,*text3D*, *ribbon3d*, *hist3D*, etc.

In addition to the x, y (and z) values, an **additional data dimension** can be represented by a color variable (argument *colvar*).

This **“4D”** plot (x, y, z, color) with a color legend is not (easily) possible using the packages mentioned above (scatterplot3d, scatter3d, rgl).

The package **plot3Drgl** allows to plot easily the graph generated with plot3D in openGL, as made available by package rgl. This is described at the end of the present article.

`install.packages("plot3D")`

`library("plot3D")`

We’ll use the *iris* data set in the following examples :

```
data(iris)
head(iris)
```

```
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
```

```
# x, y and z coordinates
x <- sep.l <- iris$Sepal.Length
y <- pet.l <- iris$Petal.Length
z <- sep.w <- iris$Sepal.Width
```

*iris* data set gives the measurements of the variables sepal length and width and petal length and width, respectively, for 50 flowers from each of 3 species of iris. The species are Iris setosa, versicolor, and virginica.

The function below will be used:

```
scatter3D(x, y, z, ..., colvar = z, col = NULL, add = FALSE)
text3D(x, y, z, labels, colvar = NULL, add = FALSE)
points3D(x, y, z, ...)
lines3D(x, y, z, ...)
scatter2D(x, y, colvar = NULL, col = NULL, add = FALSE)
text2D(x, y, labels, colvar = NULL, col = NULL, add = FALSE)
```

**x, y, z**: vectors of point coordinates**colvar**: a variable used for coloring**col**: color palette used for coloring the colvar variable**labels**: the text to be written**add**: logical. If TRUE, then the points will be added to the current plot. If FALSE a new plot is started**…**: additional*persp*arguments including xlim, ylim, zlim, xlab, ylab, zlab, main, sub, r, d, scale, expand, box, axes, nticks, tictype.

Note that:

**points3D**and**lines3D**are shorthand for scatter3D(…, type =“p”) and scatter3D(…, type = “l”), respectively.**points2D**and**lines2D**are shorthand for scatter2D(…, type = “p”) and scatter2D(…, type =“l”), respectively.

`scatter3D(x, y, z, clab = c("Sepal", "Width (cm)"))`

The argument *clab* is used to change the title of the color legend.

By default, the points are colored automatically using the variable *Z*

In the R code below:

**colvar = NULL**: avoids coloring by z variable**col = “blue”**: changes point colors to blue**pch = 19**: changes point shapes**cex = 0.5**: changes the size of points

```
scatter3D(x, y, z, colvar = NULL, col = "blue",
pch = 19, cex = 0.5)
```

The argument **bty** is used. Allowed values are:

- “f”: full box
- “b”: default value. Only the back panels are visible
- “b2”: back panels and grid lines are visible
- “g”: grey background with white grid lines
- “bl”: black background
- “bl2”: black background with grey lines
- “u”: means that the user will specify the arguments col.axis, col.panel, lwd.panel, col.grid, lwd.grid manually
- “n”: no box will be drawn. This is the same as setting box = FALSE

```
# full box
scatter3D(x, y, z, bty = "f", colkey = FALSE, main ="bty= 'f'")
# back panels and grid lines are visible
scatter3D(x, y, z, bty = "b2", colkey = FALSE, main ="bty= 'b2'" )
```

```
# grey background with white grid lines
scatter3D(x, y, z, bty = "g", colkey = FALSE, main ="bty= 'g'")
# User defined
scatter3D(x, y, z, pch = 18, bty = "u", colkey = FALSE,
main ="bty= 'u'", col.panel ="steelblue", expand =0.4,
col.grid = "darkblue")
```

The argument **colkey = FALSE** is used to remove the legend.

Several color palettes are available in **plot3D** package:

**jet.col(n, alpha)**: generates the matlab-type colors. This is the default color palette used in plot3D**jet2.col(n, alpha)**: similar to*jet.col()*but lacks the deep blue colors**gg.col(n, alpha)**and**gg2.col(n, alpha)**generates gg-plot-like colors**ramp.col(col = c(“grey”, “black”), n, alpha)**: creates color schemes by interpolation**alpha.col(col = “grey”, alpha)**: creates transparent colors

**n**: Number of colors to generate. Default value is 100**alpha**: color transparency. Value in the range 0, 1. Default value is 1**col**: Colors to interpolate

```
# gg.col: ggplot2 like color
scatter3D(x, y, z, bty = "g", pch = 18, col = gg.col(100))
# ramp.col: custom palettes
scatter3D(x, y, z, bty = "g", pch = 18,
col = ramp.col(c("blue", "yellow", "red")) )
```

The **colkey** is customized (see **?colkey** for more details):

```
scatter3D(x, y, z, bty = "g", pch = 18,
col.var = as.integer(iris$Species),
col = c("#1B9E77", "#D95F02", "#7570B3"),
pch = 18, ticktype = "detailed",
colkey = list(at = c(2, 3, 4), side = 1,
addlines = TRUE, length = 0.5, width = 0.5,
labels = c("setosa", "versicolor", "virginica")) )
```

```
# Bottom colkey
scatter3D(x, y, z, bty = "g",
colkey = list(side = 1, length = 0.5))
```

The argument *side* is used to specify the colkey position: 1: for bottom, 2: for left, 3: for top, 4: for right.

The arguments *theta* and *phi* can be used to define the angles for the viewing direction. *theta* is the azimuthal direction and *phi* the co-latitude.

`scatter3D(x, y, z, theta = 15, phi = 20)`

`scatter3D(x, y, z, phi = 0, bty ="g")`

The default values for theta and phi are 40.

```
scatter3D(x, y, z, pch = 18, theta = 20, phi = 20,
main = "Iris data", xlab = "Sepal.Length",
ylab ="Petal.Length", zlab = "Sepal.Width")
```

The arguments below can be used:

**ticktype**: Possible values are

*“simple”*draws just an arrow parallel to the axis to indicate direction of increase*“detailed”*draws normal ticks and labels

**nticks**: the number of tick marks to draw on the axes. It has no effect if ticktype =“simple”.

```
scatter3D(x, y, z, phi = 0, bty = "g",
pch = 20, cex = 2, ticktype = "detailed")
```

The functions below can be used:

**scatter3D(x, y, z,…, add = TRUE)**: Adds points**text3D(x, y, z, labels, …, add = TRUE)**: Adds texts

**Add points**to an existing plot:

```
# Create a scatter plot
scatter3D(x, y, z, phi = 0, bty = "g",
pch = 20, cex = 2, ticktype = "detailed")
# Add another point (black color)
scatter3D(x = 7, y = 3, z = 3.5, add = TRUE, colkey = FALSE,
pch = 18, cex = 3, col = "black")
```

**Add texts**to an existing plot:

```
# Create a scatter plot
scatter3D(x, y, z, phi = 0, bty = "g", pch = 20, cex = 0.5)
# Add text
text3D(x, y, z, labels = rownames(iris),
add = TRUE, colkey = FALSE, cex = 0.5)
```

```
# type ="l" for lines only
scatter3D(x, y, z, phi = 0, bty = "g", type = "l",
ticktype = "detailed", lwd = 4)
```

```
# type ="b" for both points and lines
scatter3D(x, y, z, phi = 0, bty = "g", type = "b",
ticktype = "detailed", pch = 20,
cex = c(0.5, 1, 1.5))
```

```
# type ="h" for vertical lines
scatter3D(x, y, z, phi = 0, bty = "g", type = "h",
ticktype = "detailed", pch = 19, cex = 0.5)
```

Vertical lines are useful to see clearly the x-y location of points.

The argument **CI** is used. It’s a list containing the parameters and values for the confidence intervals or NULL.

If *CI* is a list, it should contain at least the item x, y or z (latter for scatter3D).These should be 2-columned matrices, defining the left/right intervals.

Other parameters should be one of: alen = 0.01, lty = par(“lty”), lwd = par(“lwd”), col = NULL, to set the length of the arrow head, the line type and width, and the color.

If```
# Confidence interval
CI <- list(z = matrix(nrow = length(x),
data = rep(0.1, 2*length(x))))
head(CI$z)
```

```
[,1] [,2]
[1,] 0.1 0.1
[2,] 0.1 0.1
[3,] 0.1 0.1
[4,] 0.1 0.1
[5,] 0.1 0.1
[6,] 0.1 0.1
```

```
# 3D Scatter plot with CI
scatter3D(x, y, z, phi = 0, bty = "g", col = gg.col(100),
pch = 18, CI = CI)
```

A helper function **scatter3D_fancy()** is used:

```
# Add small dots on basal plane and on the depth plane
scatter3D_fancy <- function(x, y, z,..., colvar = z)
{
panelfirst <- function(pmat) {
XY <- trans3D(x, y, z = rep(min(z), length(z)), pmat = pmat)
scatter2D(XY$x, XY$y, colvar = colvar, pch = ".",
cex = 2, add = TRUE, colkey = FALSE)
XY <- trans3D(x = rep(min(x), length(x)), y, z, pmat = pmat)
scatter2D(XY$x, XY$y, colvar = colvar, pch = ".",
cex = 2, add = TRUE, colkey = FALSE)
}
scatter3D(x, y, z, ..., colvar = colvar, panel.first=panelfirst,
colkey = list(length = 0.5, width = 0.5, cex.clab = 0.75))
}
```

Fancy scatter plot:

```
scatter3D_fancy(x, y, z, pch = 16,
ticktype = "detailed", theta = 15, d = 2,
main = "Iris data", clab = c("Petal", "Width (cm)") )
```

The mtcars data will be used:

```
data(mtcars)
head(mtcars[, 1:6])
```

```
mpg cyl disp hp drat wt
Mazda RX4 21.0 6 160 110 3.90 2.620
Mazda RX4 Wag 21.0 6 160 110 3.90 2.875
Datsun 710 22.8 4 108 93 3.85 2.320
Hornet 4 Drive 21.4 6 258 110 3.08 3.215
Hornet Sportabout 18.7 8 360 175 3.15 3.440
Valiant 18.1 6 225 105 2.76 3.460
```

- Use the function
*lm()*to compute a linear regression model:*ax + by + cz + d = 0* - Use the argument
*surf*in**scatter3D()**function to add a regression surface.

```
# x, y, z variables
x <- mtcars$wt
y <- mtcars$disp
z <- mtcars$mpg
# Compute the linear regression (z = ax + by + d)
fit <- lm(z ~ x + y)
# predict values on regular xy grid
grid.lines = 26
x.pred <- seq(min(x), max(x), length.out = grid.lines)
y.pred <- seq(min(y), max(y), length.out = grid.lines)
xy <- expand.grid( x = x.pred, y = y.pred)
z.pred <- matrix(predict(fit, newdata = xy),
nrow = grid.lines, ncol = grid.lines)
# fitted points for droplines to surface
fitpoints <- predict(fit)
# scatter plot with regression plane
scatter3D(x, y, z, pch = 18, cex = 2,
theta = 20, phi = 20, ticktype = "detailed",
xlab = "wt", ylab = "disp", zlab = "mpg",
surf = list(x = x.pred, y = y.pred, z = z.pred,
facets = NA, fit = fitpoints), main = "mtcars")
```

**surf** is a list specifying a (fitted) surface to be added on the scatter plot. The list should include at least x, y, z, defining the surface.

Other optional parameters can be specified in the *surf* argument including: colvar, col, NAcol, border, facets, lwd, resfac, clim, ltheta, lphi, shade, lighting, fit. (**see ?surf3D** for more details on these parameters)

Note that, by default

**colvar = z**.- The argument
**fit**should give the fitted*z-values*, in the same order as the z-values of the scatter points, for instance produced by*predict()*. When present, this will produce droplines from points to the fitted surface.

Note that, the function **expand.grid()**, in the R code above, creates a data frame from all combinations of factors

The function **text3D()** is used as follow:

`text3D(x, y, z, labels, ...)`

The USArrests data sets will be used in the example below:

```
data(USArrests)
with(USArrests, text3D(Murder, Assault, Rape,
labels = rownames(USArrests), colvar = UrbanPop,
col = gg.col(100), theta = 60, phi = 20,
xlab = "Murder", ylab = "Assault", zlab = "Rape",
main = "USA arrests", cex = 0.6,
bty = "g", ticktype = "detailed", d = 2,
clab = c("Urban","Pop"), adj = 0.5, font = 2))
```

```
# Plot texts
with(USArrests, text3D(Murder, Assault, Rape,
labels = rownames(USArrests), colvar = UrbanPop,
col = gg.col(100), theta = 60, phi = 20,
xlab = "Murder", ylab = "Assault", zlab = "Rape",
main = "USA arrests", cex = 0.6,
bty = "g", ticktype = "detailed", d = 2,
clab = c("Urban","Pop"), adj = 0.5, font = 2))
# Add points
with(USArrests, scatter3D(Murder, Assault, Rape - 1,
colvar = UrbanPop, col = gg.col(100),
type = "h", pch = ".", add = TRUE))
```

```
# Zoom near origin: choose suitable ranges
plotdev(xlim = c(0, 10), ylim = c(40, 150),
zlim = c(7, 25))
```

Note that, in order to choose suitable ranges for zooming, you can display axis ranges as follow:

```
# display axis ranges
getplist()[c("xlim","ylim","zlim")]
```

```
$xlim
[1] 0.8 17.4
$ylim
[1] 45 337
$zlim
[1] 7.3 46.0
```

The function **hist3D()** is used:

```
hist3D (x, y, z, ..., colvar = z,
col = NULL, add = FALSE)
```

**z**: Matrix containing the values to be plotted**x, y**vectors with x and y values.*x*should be of length equal to nrow(z) and y should be equal to ncol(z)**colvar**: the variable used for coloring. If present, it should have the same dimension as z.**col**: color palette to be used for the colvar variable. By default a*red-yellow-blue*color scheme (?jet.col) is used**add**: Logical. If TRUE, then the surfaces will be added to the current plot. If FALSE a new plot is started.

```
data(VADeaths)
# hist3D and ribbon3D with greyish background, rotated, rescaled,...
hist3D(z = VADeaths, scale = FALSE, expand = 0.01, bty = "g", phi = 20,
col = "#0072B2", border = "black", shade = 0.2, ltheta = 90,
space = 0.3, ticktype = "detailed", d = 2)
```

```
hist3D (x = 1:5, y = 1:4, z = VADeaths,
bty = "g", phi = 20, theta = -60,
xlab = "", ylab = "", zlab = "", main = "VADeaths",
col = "#0072B2", border = "black", shade = 0.8,
ticktype = "detailed", space = 0.15, d = 2, cex.axis = 1e-9)
# Use text3D to label x axis
text3D(x = 1:5, y = rep(0.5, 5), z = rep(3, 5),
labels = rownames(VADeaths),
add = TRUE, adj = 0)
# Use text3D to label y axis
text3D(x = rep(1, 4), y = 1:4, z = rep(0, 4),
labels = colnames(VADeaths),
add = TRUE, adj = 1)
```

**fancy 3D** histograms

```
hist3D_fancy<- function(x, y, break.func = c("Sturges", "scott", "FD"), breaks = NULL,
colvar = NULL, col="white", clab=NULL, phi = 5, theta = 25, ...){
# Compute the number of classes for a histogram
break.func <- break.func [1]
if(is.null(breaks)){
x.breaks <- switch(break.func,
Sturges = nclass.Sturges(x),
scott = nclass.scott(x),
FD = nclass.FD(x))
y.breaks <- switch(break.func,
Sturges = nclass.Sturges(y),
scott = nclass.scott(y),
FD = nclass.FD(y))
} else x.breaks <- y.breaks <- breaks
# Cut x and y variables in bins for counting
x.bin <- seq(min(x), max(x), length.out = x.breaks)
y.bin <- seq(min(y), max(y), length.out = y.breaks)
xy <- table(cut(x, x.bin), cut(y, y.bin))
z <- xy
xmid <- 0.5*(x.bin[-1] + x.bin[-length(x.bin)])
ymid <- 0.5*(y.bin[-1] + y.bin[-length(y.bin)])
oldmar <- par("mar")
par (mar = par("mar") + c(0, 0, 0, 2))
hist3D(x = xmid, y = ymid, z = xy, ...,
zlim = c(-max(z)/2, max(z)), zlab = "counts", bty= "g",
phi = phi, theta = theta,
shade = 0.2, col = col, border = "black",
d = 1, ticktype = "detailed")
scatter3D(x, y,
z = rep(-max(z)/2, length.out = length(x)),
colvar = colvar, col = gg.col(100),
add = TRUE, pch = 18, clab = clab,
colkey = list(length = 0.5, width = 0.5,
dist = 0.05, cex.axis = 0.8, cex.clab = 0.8)
)
par(mar = oldmar)
}
```

```
hist3D_fancy(quakes$long, quakes$lat, colvar=quakes$depth,
breaks =30)
```

```
hist3D_fancy(iris$Sepal.Length, iris$Petal.Width,
colvar=as.numeric(iris$Species))
```

**Create some data**:

```
# x, y coordinates
set.seed(1234)
x <- sort(rnorm(10))
y <- runif(10)
# Variable for coloring points
col.v <- sqrt(x^2 + y^2)
```

**Basic 2D scatter plot**:

```
scatter2D(x, y, colvar = col.v, pch = 16, bty ="n",
type ="b")
```

**type**: plot types. Allowed values are:

**“b”**to draw both points and line**“h”**for vertical line**“l”**for line only**“p”**for points only

**bty**: box type

**2D scatter plot with confidence interval**:

```
# Confidence interval for x variable only
CI <- list()
CI$x <- matrix(nrow = length(x), data = c(rep(0.25, 2*length(x))))
scatter2D(x, y, colvar = col.v, pch = 16, bty ="n", cex = 1.5,
CI = CI, type = "b")
```

```
# Confidence interval for both x and y variables
CI$y <- matrix (nrow = length(y), data = c(rep(0.05, 2*length(y))))
CI$col <- "black"
scatter2D(x, y, colvar = col.v, pch = 16, bty ="n", cex = 1.5,
CI = CI, type ="b")
```

```
CI$y[c(2,4,8,10), ] <- NA # Some points have no CI
CI$x[c(2,4,8,10), ] <- NA # Some points have no CI
CI$alen <- 0.02 # increase arrow head
scatter2D(x, y, colvar = col.v, pch = 16, bty ="n", cex = 1.5,
CI = CI, type ="b")
```

```
# Only text
with(USArrests, text2D(x = Murder, y = Assault + 5, colvar = Rape,
xlab = "Murder", ylab = "Assault", clab = "Rape",
main = "USA arrests", labels = rownames(USArrests), cex = 0.6,
adj = 0.5, font = 2))
```

```
# text with point
with(USArrests, text2D(x = Murder, y = Assault + 5, colvar = Rape,
xlab = "Murder", ylab = "Assault", clab = "Rape",
main = "USA arrests", labels = rownames(USArrests), cex = 0.6,
adj = 0.5, font = 2))
with(USArrests, scatter2D(x = Murder, y = Assault, colvar = Rape,
pch = 16, add = TRUE, colkey = FALSE))
```

It’s also possible to draw arrows, segments and rectangles in a 3D or 2D plot using the functions below:

```
arrows3D(x0, y0, z0, x1, y1, z1, ..., colvar = NULL,
col = NULL, type = "triangle", add = FALSE)
segments3D(x0, y0, z0, x1, y1, z1, ..., colvar = NULL,
col = NULL, add = "FALSE")
rect3D(x0, y0, y0, x1, y1, z1, ..., colvar = NULL,
col = NULL, add = FALSE)
arrows2D(x0, y0, z0, x1, y1, z1, ..., colvar = NULL,
col = NULL, type = "triangle", add = FALSE)
segments2D(x0, y0, z0, x1, y1, z1, ..., colvar = NULL,
col = NULL, add = "FALSE")
rect2D(x0, y0, y0, x1, y1, z1, ..., colvar = NULL,
col = NULL, add = FALSE)
```

**x0, y0, z0**: coordinates of points from which to draw**x1**,**y1**,**z1**: coordinates of points to which to draw. For arrows3D and segments3D, at least one must be supplied. For rect3D exactly one must be NULL.**colvar**: The variable used for coloring.**col**: color palette to be used for coloring. Default is*red-yellow-blue*color scheme.**add**: Logical. If TRUE, then the arrows, segments, … will be added to the current plot. If FALSE a new plot is started.

Prepare the data: we want to plot 4 arrows starting from the point of coordinates c(x0, y0, z0) and ending at c(x1, y1, z1)

```
x0 <- c(0, 0, 0, 0)
y0 <- c(0, 0, 0, 0)
z0 <- c(0, 0, 0, 0)
x1 <- c(0.89, -0.46, 0.99, 0.96)
y1 <- c(0.36, 0.88, 0.02, 0.06)
z1 <- c(-0.28, 0.09, 0.05, 0.24)
cols <- c("#1B9E77", "#D95F02", "#7570B3", "#E7298A")
```

**3D Arrows**:

```
arrows3D(x0, y0, z0, x1, y1, z1, colvar = x1^2, col = cols,
lwd = 2, d = 3, clab = c("Quality", "score"),
main = "Arrows 3D", bty ="g", ticktype = "detailed")
# Add starting point of arrow
points3D(x0, y0, z0, add = TRUE, col="darkred",
colkey = FALSE, pch = 19, cex = 1)
# Add labels to the arrows
text3D(x1, y1, z1, c("Sepal.L", "Sepal.W", "Petal.L", "Petal.W"),
colvar = x1^2, col = cols, add=TRUE, colkey = FALSE)
```

**2D arrows:**

```
arrows2D(x0, y0, x1, y1, colvar = x1^2, col = cols,
lwd = 2, clab = c("Quality", "score"),
bty ="n", xlim = c(-1, 1), ylim = c(-1, 1))
# Add vertical and horizontal lines at c(0,0)
abline(h =0, v = 0, lty = 2)
# Add starting point of arrow
points2D(x0, y0, add = TRUE, col="darkred",
colkey = FALSE, pch = 19, cex = 1)
# Add labels to the arrows
text2D(x1, y1, c("Sepal.L", "Sepal.W", "Petal.L", "Petal.W"),
colvar = x1^2, col = cols, add=TRUE, colkey = FALSE)
```

Note that, **segments3D()** and **segments2D()** are very similar to **arrows3D()** and **arrows2D()** and you can play with them also.

**3D rectangle**: the R code below creates a rectangle with a transparent fill color (alpha = 0.5)

```
rect3D(x0 = 0, y0 = 0.5, z0 = 0, x1 = 1, z1 = 5,
ylim = c(0, 1), bty = "g", facets = TRUE,
border = "red", col ="#7570B3", alpha=0.5,
lwd = 2, phi = 20)
```

In the R code above, *facets = FALSE*, will remove the rectangle fill color.

**2D rectangle**:

```
rect2D(x0 = runif(3), y0 = runif(3),
x1 = runif(3), y1 = runif(3), colvar = 1:3,
alpha = 0.4, lwd = 2, main = "rect2D")
```

To draw an interactive **3D plot** the package **plot3Drgl** can be used.

The package **plot3Drgl** allows to plot the graph generated with **plot3D** in openGL, as made available by package rgl.

The simplest way is to do as follow:

- Create base R-graphics using plot3D package
- Then use the function plotrgl() to draw the same figure in rgl

The package **rgl** allows to interactively rotate, zoom the graphs. However it’s not yet possible to plot a colorkey

```
# Create his3D using plot3D
hist3D_fancy(iris$Sepal.Length, iris$Petal.Width, colvar=as.numeric(iris$Species))
# Make the rgl version
library("plot3Drgl")
plotrgl()
```

Note that, after creating the rgl plot you can use the functions below:

- croprgl(xlim, ylim, zlim, …) to modify the ranges
- cutrgl(…) to zoom in on a selected region of the plot. The current plot will be overwritten
uncutrgl(…) and uncroprgl(…) restore the original plot.

- …: any arguments for par3d, open3d or material3d in rgl package.

This analysis has been performed using **R software** (ver. 3.1.2) and **plot3D** (ver. 1.0-2)

**References**:

- Karline Soetaert. plot3D: Tools for plotting 3-D and 2-D data. http://cran.r-project.org/web/packages/plot3D/vignettes/plot3D.pdf

- Install the RGL package
- Load the RGL package
- Prepare the data
- Start and close RGL device
- 3D scatter plot
- rgl_init(): A custom function to initialize RGL device
- Add a bounding box decoration
- Add axis lines and labels
- Set the aspect ratios of the x, y and z axes
- Change the color of points by groups
- Change the shape of points
- Add an ellipse of concentration
- Regression plane
- Create a movie of RGL scene
- Export images as png or pdf
- Export the plot into an interactive HTML file
- Select a rectangle in an RGL scene
- Identify points in a plot
- R3D Interface
- RGL functions
- Infos

This **R tutorial** describes, step by step, how to build a **3D graphic** using **R software** and the **rgl** package. You’ll learn also how to create a **movie** of your **3D scene** in **R**.

**RGL** is a **3D graphics** package that produces a real-time interactive **3D plot**. It allows to interactively rotate, zoom the graphics and select regions.

The **rgl package** includes also a *generic 3D interface* named **R3D**. **R3D** is a collection of generic 3D objects and functions which are described at the end of this article.

`install.packages("rgl")`

Note that, on Linux operating system, the **rgl** package can be installed as follow:

**sudo apt-get install r-cran-rgl**

`library("rgl")`

We’ll use the *iris* data set in the following examples:

```
data(iris)
head(iris)
```

```
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
```

```
x <- sep.l <- iris$Sepal.Length
y <- pet.l <- iris$Petal.Length
z <- sep.w <- iris$Sepal.Width
```

*iris* data set gives the measurements of the variables sepal length and width, petal length and width, respectively, for 50 flowers from each of 3 species of iris. The species are Iris setosa, versicolor, and virginica.

To make a 3D plot with RGL, you should first start the RGL device in R.

The functions below are used to manage the **RGL** device:

**rgl.open()**: Opens a new device**rgl.close()**: Closes the current device**rgl.clear()**: Clears the current device**rgl.cur()**: Returns the active device ID**rgl.quit()**: Shutdowns the RGL device system

In the first sections of this tutorial, I’ll open a new RGL device for each plot.

Note that, you don’t need to do the same thing.

You can just use the functionThe function **rgl.points()** is used to draw a **3D scatter plot**:

```
rgl.open() # Open a new RGL device
rgl.points(x, y, z, color ="lightgray") # Scatter plot
```

**x, y, z**: Numeric vector specifying the coordinates of points to be drawn. The arguments y and z are optional when:*x*is a matrix or a data frame containing at least 3 columns which will be used as the x, y and z coordinates. Ex:*rgl.points(iris)**x*is a formula of form zvar ~ xvar + yvar (see**?xyz.coords**). Ex:*rgl.points( z ~ x + y)*.**…**: Material properties. See**?rgl.material**for details.

- The function
**rgl.bg(color)**can be used to setup the background environment of the scene - The argument
*color*is used in the function*rgl.points()*to change point colors

Note that, it’s also possible to change the **size of points** using the argument **size**

```
rgl.open()# Open a new RGL device
rgl.bg(color = "white") # Setup the background color
rgl.points(x, y, z, color = "blue", size = 5) # Scatter plot
```

Note that, the equivalent of the functions above for the *3d interface* is:

**open3d()**: Open a new 3D device**bg3d(color)**: Set up the background environment of the scene**points3d(x, y, Z, …)**: plot points of coordinates x, y, z

It’s possible to draw **spheres** using the functions **rgl.spheres()** or **spheres3d()**:

```
spheres3d(x, y = NULL, z = NULL, radius = 1, ...)
rgl.spheres(x, y = NULL, z = NULL, r, ...)
```

**rgl.spheres()** draws spheres with center (x, y, z) and radius r.

**x, y, z**: Numeric vector specifying the coordinates for the center of each sphere. The arguments y and z are optional when:*x*is a matrix or a data frame containing at least 3 columns which will be used as the x, y and z coordinates. Ex:*rgl.spheres(iris, r = 0.2)**x*is a formula of form zvar ~ xvar + yvar (see**?xyz.coords**). Ex:*rgl.spheres( z ~ x + y, r = 0.2)*.**radius**: a vector or single value indicating the radius of spheres**…**: Material properties. See**?rgl.material**for details.

```
rgl.open()# Open a new RGL device
rgl.bg(color = "white") # Setup the background color
rgl.spheres(x, y, z, r = 0.2, color = "grey")
```

Note that, an **rgl** plot can be manually rotated by holding down on the mouse or touchpad. It can be also zoomed using the scroll wheel on a mouse or pressing ctrl + using the touchpad on a PC or two fingers (up or down) on a mac.

The function **rgl_init()** will create a new **RGL** device if requested or if there is no opened device:

```
#' @param new.device a logical value. If TRUE, creates a new device
#' @param bg the background color of the device
#' @param width the width of the device
rgl_init <- function(new.device = FALSE, bg = "white", width = 640) {
if( new.device | rgl.cur() == 0 ) {
rgl.open()
par3d(windowRect = 50 + c( 0, 0, width, width ) )
rgl.bg(color = bg )
}
rgl.clear(type = c("shapes", "bboxdeco"))
rgl.viewpoint(theta = 15, phi = 20, zoom = 0.7)
}
```

**Description of the used RGL functions**:

**rgl.open()**: open a new device**rgl.cur()**: returns active device ID**par3d(windowRect)**: set the window size**rgl.viewpoint(theta, phi, fov, zoom)**: set up viewpoint. The arguments*theta*and*phi*are polar coordinates.**theta**and**phi**are the polar coordinates. Default values are 0 and 15, respectively**fov**is the field-of-view angle in degrees. Default value is 60**zoom**is the zoom factor. Default value is 1**rgl.bg(color)**: define the background color of the device**rgl.clear(type)**: Clears the scene from the specified stack (“shapes”, “lights”, “bboxdeco”, “background”)

In the **R code** above, I used the function **rgl.viewpoint()** to set automatically the viewpoint orientation and the zoom. As you already know, the **RGL** device is interactive and you can adjust the viewpoint and zoom the plot using your mouse.

The function **rgl.bbox()** is used:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow") # Scatter plot
rgl.bbox(color = "#333377") # Add bounding box decoration
```

A simplified format of the function **rgl.bbox()** is:

`rgl.bbox(xlen=5, ylen=5, zlen=5, marklen=15.9, ...)`

**xlen**,**ylen**,**zlen**: values specifying the number of tickmarks on x, y and Z axes, respectively**marklen**: value specifying the length of the tickmarks**…**: other rgl material properties (see*?rgl.material*) including:**color**: a vector of colors. The first color is used for the background color of the bounding box. The second color is used for the tick mark labels.**emission**,**specular**,**shininess**: properties for lighting calculation**alpha**: value specifying the color transparency. The value should be between 0.0 (fully transparent) and 1.0 (opaque)

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
# Add bounding box decoration
rgl.bbox(color=c("#333377","black"), emission="#333377",
specular="#3333FF", shininess=5, alpha=0.8 )
```

- The function
**rgl.lines(x, y = NULL, z = NULL, …)**can be used to add axis lines. - The function
**rgl.texts(x, y = NULL, z = NULL, text)**is used to add axis labels

For the function **rgl.lines()**, the arguments x, y, and z are numeric vectors of length 2 (i.e, : x = c(x1,x2), y = c(y1, y2), z = c(z1, z2) ).

- The values x1, y1 and y3 are the
**3D coordinates**of the line starting point. - The values x2, y2 and y3 corresponds to the
**3D coordinates**of the line ending point.

Note also that, the argument *x* can be a matrix or a data frame containing, at least, 3 columns corresponding to the x, y and z coordinates, respectively. In this case, the argument *y* and *z* can be omitted.

To draw an axis, you should specify the range (minimum and the maximum) of the axis to the function **rgl.lines()**:

```
# Make a scatter plot
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
# Add x, y, and z Axes
rgl.lines(c(min(x), max(x)), c(0, 0), c(0, 0), color = "black")
rgl.lines(c(0, 0), c(min(y),max(y)), c(0, 0), color = "red")
rgl.lines(c(0, 0), c(0, 0), c(min(z),max(z)), color = "green")
```

As you can see, the axes are drawn but the problem is that they don’t cross at the point c(0, 0, 0)

There are two solutions to handle this situation:

**Scale the data**to make things easy. Transform the x, y and z variables so that their min = 0 and their max = 1**Use c(-max, +max)**as the ranges of the axes

```
x1 <- (x - min(x))/(max(x) - min(x))
y1 <- (y - min(y))/(max(y) - min(y))
z1 <- (z - min(z))/(max(z) - min(z))
```

```
# Make a scatter plot
rgl_init()
rgl.spheres(x1, y1, z1, r = 0.02, color = "yellow")
# Add x, y, and z Axes
rgl.lines(c(0, 1), c(0, 0), c(0, 0), color = "black")
rgl.lines(c(0, 0), c(0,1), c(0, 0), color = "red")
rgl.lines(c(0, 0), c(0, 0), c(0,1), color = "green")
```

Let’s define a helper function to calculate the axis limits:

`lim <- function(x){c(-max(abs(x)), max(abs(x))) * 1.1}`

```
# Make a scatter plot
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
# Add x, y, and z Axes
rgl.lines(lim(x), c(0, 0), c(0, 0), color = "black")
rgl.lines(c(0, 0), lim(y), c(0, 0), color = "red")
rgl.lines(c(0, 0), c(0, 0), lim(z), color = "green")
```

```
# x, y, z : numeric vectors corresponding to
# the coordinates of points
# axis.col : axis colors
# xlab, ylab, zlab: axis labels
# show.plane : add axis planes
# show.bbox : add the bounding box decoration
# bbox.col: the bounding box colors. The first color is the
# the background color; the second color is the color of tick marks
rgl_add_axes <- function(x, y, z, axis.col = "grey",
xlab = "", ylab="", zlab="", show.plane = TRUE,
show.bbox = FALSE, bbox.col = c("#333377","black"))
{
lim <- function(x){c(-max(abs(x)), max(abs(x))) * 1.1}
# Add axes
xlim <- lim(x); ylim <- lim(y); zlim <- lim(z)
rgl.lines(xlim, c(0, 0), c(0, 0), color = axis.col)
rgl.lines(c(0, 0), ylim, c(0, 0), color = axis.col)
rgl.lines(c(0, 0), c(0, 0), zlim, color = axis.col)
# Add a point at the end of each axes to specify the direction
axes <- rbind(c(xlim[2], 0, 0), c(0, ylim[2], 0),
c(0, 0, zlim[2]))
rgl.points(axes, color = axis.col, size = 3)
# Add axis labels
rgl.texts(axes, text = c(xlab, ylab, zlab), color = axis.col,
adj = c(0.5, -0.8), size = 2)
# Add plane
if(show.plane)
xlim <- xlim/1.1; zlim <- zlim /1.1
rgl.quads( x = rep(xlim, each = 2), y = c(0, 0, 0, 0),
z = c(zlim[1], zlim[2], zlim[2], zlim[1]))
# Add bounding box decoration
if(show.bbox){
rgl.bbox(color=c(bbox.col[1],bbox.col[2]), alpha = 0.5,
emission=bbox.col[1], specular=bbox.col[1], shininess=5,
xlen = 3, ylen = 3, zlen = 3)
}
}
```

- The function
**rgl.texts(x, y, z, text )**is used to add texts to an**RGL plot**. **rgl.quads(x, y, z)**is used to add planes. x, y and z are numeric vectors of length four specifying the coordinates of the four nodes of the quad.

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
rgl_add_axes(x, y, z)
```

The function **axis3d()** can be used as follow:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
rgl_add_axes(x, y, z)
# Show tick marks
axis3d('x', pos=c( NA, 0, 0 ), col = "darkgrey")
axis3d('y', pos=c( 0, NA, 0 ), col = "darkgrey")
axis3d('z', pos=c( 0, 0, NA ), col = "darkgrey")
```

It’s easier to just add a bounding box decoration:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
rgl_add_axes(x, y, z, show.bbox = TRUE)
```

In the plot above, the bounding box is displayed as rectangle. All the coordinates are displayed at the same scale (*iso-metric*).

The function **aspect3d(x, y = NULL, z = NULL)** can be used to set the apparent ratios of the x, y and z axes for the current plot.

x, y and z are the ratio for x, y and z axes, respectively. x can be a vector of length 3, specifying the ratio for the 3 axes.

If the ratios are (1, 1, 1), the bounding box will be displayed as a cube.

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
rgl_add_axes(x, y, z, show.bbox = TRUE)
aspect3d(1,1,1)
```

Note that, the default display corresponds to the aspect “iso”: *aspect3d(“iso”).*

The values of the ratios can be set larger or smaller to zoom on a given axis:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "yellow")
rgl_add_axes(x, y, z, show.bbox = TRUE)
aspect3d(2,1,1) # zoom on x axis
```

A helper function can be used to select automatically a color for each group:

```
#' Get colors for the different levels of
#' a factor variable
#'
#' @param groups a factor variable containing the groups
#' of observations
#' @param colors a vector containing the names of
# the default colors to be used
get_colors <- function(groups, group.col = palette()){
groups <- as.factor(groups)
ngrps <- length(levels(groups))
if(ngrps > length(group.col))
group.col <- rep(group.col, ngrps)
color <- group.col[as.numeric(groups)]
names(color) <- as.vector(groups)
return(color)
}
```

Change colors by groups :

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2,
color = get_colors(iris$Species))
rgl_add_axes(x, y, z, show.bbox = TRUE)
aspect3d(1,1,1)
```

Use custom colors:

```
cols <- get_colors(iris$Species, c("#999999", "#E69F00", "#56B4E9"))
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = cols)
rgl_add_axes(x, y, z, show.bbox = TRUE)
aspect3d(1,1,1)
```

It’s also possible to use color palettes from the **RColorBrewer** package:

```
library("RColorBrewer")
cols <- get_colors(iris$Species, brewer.pal(n=3, name="Dark2") )
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = cols)
rgl_add_axes(x, y, z, show.bbox = TRUE)
aspect3d(1,1,1)
```

Read more about colors in R: colors in R

6 mesh objects are available in RGL package and can be used as point shapes:

**cube3d()****tetrahedron3d()****octahedron3d()****icosahedron3d()****dodecahedron3d()****cuboctahedron3d()**

To make a plot using the objects above, the function **shapelist3d()** can be used as follow:

`shapelist3d(shapes, x, y, z)`

**shapes**: a single shape3d (ex: shapes = cube3d()) object or a list of them (ex: shapes = list(cube3d(), icosahedron3d()))**x, y, z**: the coordinates of the points to be plotted

```
rgl_init()
shapelist3d(tetrahedron3d(), x, y, z, size = 0.15,
color = get_colors(iris$Species))
rgl_add_axes(x, y, z, show.bbox = TRUE)
aspect3d(1,1,1)
```

The function **ellipse3d()** is used to estimate the ellipse of concentration. A simplified format is:

```
ellipse3d(x, scale = c(1,1,1), centre = c(0,0,0),
level = 0.95, ...)
```

**x**: the correlation or covariance matrix between x, y and z**scale**: If x is a correlation matrix, then the standard deviations of each parameter can be given in the scale parameter. This defaults to c(1, 1, 1), so no rescaling will be done.**centre**: The center of the ellipse will be at this position.**level**: The confidence level of a confidence region. This is used to control the size of the ellipsoid.

The function **ellipse3d()** returns an object of class **mesh3d** which can be drawn using the function **shade3d()** and/or **wired3d()**

- Draw the ellipse using the function
**shade3d()**:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = TRUE)
# Compute and draw the ellipse of concentration
ellips <- ellipse3d(cov(cbind(x,y,z)),
centre=c(mean(x), mean(y), mean(z)), level = 0.95)
shade3d(ellips, col = "#D95F02", alpha = 0.1, lit = FALSE)
aspect3d(1,1,1)
```

- Draw the ellipse using the function
**wired3d()**:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = TRUE)
# Compute and draw the ellipse of concentration
ellips <- ellipse3d(cov(cbind(x,y,z)),
centre=c(mean(x), mean(y), mean(z)), level = 0.95)
wire3d(ellips, col = "#D95F02", lit = FALSE)
aspect3d(1,1,1)
```

- Combine
**shade3d()**and**wired3d()**:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = TRUE)
# Compute and draw the ellipse of concentration
ellips <- ellipse3d(cov(cbind(x,y,z)),
centre=c(mean(x), mean(y), mean(z)), level = 0.95)
shade3d(ellips, col = "#D95F02", alpha = 0.1, lit = FALSE)
wire3d(ellips, col = "#D95F02", lit = FALSE)
aspect3d(1,1,1)
```

**Add ellipse for each group**:

```
# Groups
groups <- iris$Species
levs <- levels(groups)
group.col <- c("red", "green", "blue")
# Plot observations
rgl_init()
rgl.spheres(x, y, z, r = 0.2,
color = group.col[as.numeric(groups)])
rgl_add_axes(x, y, z, show.bbox = FALSE)
# Compute ellipse for each group
for (i in 1:length(levs)) {
group <- levs[i]
selected <- groups == group
xx <- x[selected]; yy <- y[selected]; zz <- z[selected]
ellips <- ellipse3d(cov(cbind(xx,yy,zz)),
centre=c(mean(xx), mean(yy), mean(zz)), level = 0.95)
shade3d(ellips, col = group.col[i], alpha = 0.1, lit = FALSE)
# show group labels
texts3d(mean(xx),mean(yy), mean(zz), text = group,
col= group.col[i], cex = 2)
}
aspect3d(1,1,1)
```

The function **planes3d()** or **rgl.planes()** can be used to add regression plane into 3D rgl plot:

```
rgl.planes(a, b = NULL, c = NULL, d = 0, ...)
planes3d(a, b = NULL, c = NULL, d = 0, ...)
```

**planes3d()** and **rgl.planes()** draw planes using the parameter *ax + by + cz + d = 0*.

**a, b, c**: coordinates of the normal to the plane**d**: coordinates of the offset

Example of usage:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = FALSE)
aspect3d(1,1,1)
# Linear model
fit <- lm(z ~ x + y)
coefs <- coef(fit)
a <- coefs["x"]; b <- coefs["y"]; c <- -1
d <- coefs["(Intercept)"]
rgl.planes(a, b, c, d, alpha=0.2, color = "#D95F02")
```

The regression plane above is very ugly. Let’s try to do a custom one. The steps below are followed:

- Use the function
*lm()*to compute a linear regression model:*ax + by + cz + d = 0* - Use the argument
**rgl.surface()**to add a regression surface.

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = FALSE)
aspect3d(1,1,1)
# Compute the linear regression (y = ax + bz + d)
fit <- lm(y ~ x + z)
# predict values on regular xz grid
grid.lines = 26
x.pred <- seq(min(x), max(x), length.out = grid.lines)
z.pred <- seq(min(z), max(z), length.out = grid.lines)
xz <- expand.grid( x = x.pred, z = z.pred)
y.pred <- matrix(predict(fit, newdata = xz),
nrow = grid.lines, ncol = grid.lines)
# Add regression surface
rgl.surface(x.pred, z.pred, y.pred, color = "steelblue",
alpha = 0.5, lit = FALSE)
# Add grid lines
rgl.surface(x.pred, z.pred, y.pred, color = "black",
alpha = 0.5, lit = FALSE, front = "lines", back = "lines")
```

The function **movie3d()** can be used as follow:

`movie3d(f, duration, dir = tempdir(), convert = TRUE)`

**f**a function created using**spin3d(axis)****axis**: the desired axis of rotation. Default value is c(0, 0, 1).**duration**: the duration of the animation**dir**: A directory in which to create temporary files for each frame of the movie**convert**: If TRUE, tries to convert the frames to a single GIF movie. It uses**ImageMagick**for the image conversion.

You should install *ImageMagick* (http://www.imagemagick.org/) to be able to create a movie from a list of png file.

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = TRUE)
# Compute and draw the ellipse of concentration
ellips <- ellipse3d(cov(cbind(x,y,z)),
centre=c(mean(x), mean(y), mean(z)), level = 0.95)
wire3d(ellips, col = "#D95F02", lit = FALSE)
aspect3d(1,1,1)
# Create a movie
movie3d(spin3d(axis = c(0, 0, 1)), duration = 3,
dir = getwd())
```

The plot can be saved as png or pdf.

- The function
**rgl.snapshot()**is used to save the screenshot as png file:

`rgl.snapshot(filename = "plot.png")`

- The function
**rgl.postscript()**is used to save the screenshot to a file in**ps, eps, tex, pdf, svg or pgf**format:

`rgl.postscript("plot.pdf",fmt="pdf")`

Example of usage:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = T)
aspect3d(1,1,1)
rg.snapshot("plot.png")
```

The function **writeWebGL()** is used to write the current scene to HTML:

`writeWebGL(dir = "webGL", filename = file.path(dir, "index.html"))`

**dir**: Where to write the files**filename**: The file name to use for the main file

The R code below, writes a copy of the scene and then displays it in a browser:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2,
color = get_colors(iris$Species))
rgl_add_axes(x, y, z, show.bbox = FALSE)
# This writes a copy into temporary directory 'webGL',
# and then displays it
browseURL(
paste("file://", writeWebGL(dir=file.path(tempdir(), "webGL"),
width=500), sep="")
)
```

The functions **rgl.select3d()** or **select3d()** can be used to select **3-dimensional** regions.

They return a function f(x, y, z) which tests whether each of the points (x, y, z) is in the selected region.

The R code below, allows the user to select some points, and then redraw them in a different color:

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = F)
aspect3d(1,1,1)
# Select a point
f <- select3d()
sel <- f(x,y,z)
rgl.clear("shapes")
# Redraw the points
rgl.spheres(x[!sel],y[!sel], z[!sel], r = 0.2, color = "#D95F02")
rgl.spheres(x[sel],y[sel], z[sel], r = 0.2, color = "green")
```

The function **identify3d()** is used:

`identify3d(x, y = NULL, z = NULL, labels, n)`

**x, y, z**: Numeric vector specifying the coordinates of points. The arguments y and z are optional when:*x*is a matrix or a data frame containing at least 3 columns which will be used as the x, y and z coordinates.*x*is a formula of form zvar ~ xvar + yvar (see**?xyz.coords**).**labels**an optional character vector giving labels for points*n*the maximum number of points to be identified

The function identify3d(), works similarly to the identify function in base graphics.

The R code below, allow the user to identify 5 points :

```
rgl_init()
rgl.spheres(x, y, z, r = 0.2, color = "#D95F02")
rgl_add_axes(x, y, z, show.bbox = F)
aspect3d(1,1,1)
rgl.material(color = "blue")
identify3d(x, y, z, labels = rownames(iris), n = 5)
```

Use the right button to select, the middle button to quit.

The **rgl** package also includes a higher level interface called **r3d**. This interface is designed to act more like classic 2D R graphics.

The next sections describe how to make 3D graphics using the R3D interface.

The function **plot3d()** is used:

```
## Default method
plot3d(x, y, z, xlab, ylab, zlab, type = "p", col,
size, lwd, radius, add = FALSE, aspect = !add, ...)
## Method for class 'mesh3d'
plot3d(x, xlab = "x", ylab = "y", zlab = "z",
type = c("shade", "wire", "dots"), add = FALSE, ...)
decorate3d(xlim, ylim, zlim,
xlab = "x", ylab = "y", zlab = "z",
box = TRUE, axes = TRUE, main = NULL, sub = NULL,
top = TRUE, aspect = FALSE, expand = 1.03, ...)
```

**x, y, z**: vectors of points to be drawn. Any reasonable way of defining the coordinates is acceptable. See the function**xyz.coords**for details**xlab, yab, zlab**: x, y and z axis labels**type**:**For the default method**: Allowed values are: ‘p’ for points, ‘s’ for spheres, ‘l’ for lines, ‘h’ for line segments from z = 0, and ‘n’ for nothing.**For the mesh3d method**, one of ‘shade’, ‘wire’, or ‘dots’**col**: the color to be used for plotted items**size**: size of points**lwd**: the line width for plotted item**radius**: the radius of sphere**add**: whether to add the points to an existing plot**aspect**: either a logical indicating whether to adjust the aspect ratio, or a new ratio**…**: additional parameters which will be passed to par3d, material3d or decorate3d**box, axes**: whether to draw a box and axes.**main, sub**: main title and subtitle**top**: whether to bring the window to the top when done

Note that, it’s recommended to use the function **open3d()** to initialize the *3d interface. However, in the following R code chunks, I’ll continue to use the custom function **rgl_init()**.

**Draw points**:

```
rgl_init()
plot3d(x, y, z, col="blue", type ="p")
```

**Remove the box and draw spheres**:

```
rgl_init()
plot3d(x, y, z, col="blue", box = FALSE,
type ="s", radius = 0.15)
```

To remove the axes, use the argument **axes = FALSE**.

**Axis labels**:

```
rgl_init()
plot3d(x, y, z, col="blue", box = FALSE,
type ="s", radius = 0.15, xlab ="Sepal.Length",
ylab = "Petal.Length", zlab = "Sepal.Width")
```

**Add ellipse of concentration:**

```
rgl_init()
plot3d(x, y, z, col="blue", box = FALSE,
type ="s", radius = 0.15)
ellips <- ellipse3d(cov(cbind(x,y,z)),
centre=c(mean(x), mean(y), mean(z)), level = 0.95)
plot3d(ellips, col = "blue", alpha = 0.2, add = TRUE, box = FALSE)
```

**Change the ellipse type**: possible values for the argument **type = c(“shade”, “wire”, “dots”)**

```
rgl_init()
plot3d(x, y, z, col="blue", box = FALSE,
type ="s", radius = 0.15)
ellips <- ellipse3d(cov(cbind(x,y,z)),
centre = c(mean(x), mean(y), mean(z)), level = 0.95)
plot3d(ellips, col = "blue", alpha = 0.5, add = TRUE, type = "wire")
```

**bbox3d(): Add bounding box decoration**

```
rgl_init()
plot3d(x, y, z, col="blue", box = FALSE,
type ="s", radius = 0.15)
# Add bounding box decoration
rgl.bbox(color=c("#333377","black"), emission="#333377",
specular="#3333FF", shininess=5, alpha=0.8, nticks = 3 )
# Change the axis aspect ratios
aspect3d(1,1,1)
```

**open3d()**: Open a new device**Add a shape to the current scene**:- points3d(x, y, z, …)
- lines3d(x, y, z, …)
- segments3d(x, y, z, …)
- quads3d(x, y, z, …)
**texts3d(x, y, z, text)**: Adds text**Draw boxes, axes and other text**:**axes3d()**: Add standard axes**title3d(main, sub, xlab, ylab, zlab)**: Add titles**box3d()**: Draws a box around the plot**mtext3d()**: to put a text in the plot margin(s)

Some important functions for managing RGL device are listed below:

Functions | Description |
---|---|

rgl.open() | Opens a new device |

rgl.close() | Closes the current device |

rgl.cur() | Returns the ID of the active device |

rgl.dev.list() | Returns all device IDs |

rgl.set(which) | Set a device as active |

rgl.quit() | Shutdown the RGL device system |

The equivalents in r3d interface:

Adds a shape node to the current scene.

Functions | Description |
---|---|

rgl.points(x, y, z, …) | Draws a point at x, y and z |

rgl.lines(x, y, z, …) | Draws line segments based on pairs of vertices (x = c(x1,x2), y = c(y1,y2), z = c(z1, z2)) |

rgl.triangles(x, y, z, …) | Draws triangles with nodes (xi, yi, zi), i = 1, 2, 3 |

rgl.quads(x, y, z) | Draws quads with nodes (xi, yi, zi), i = 1, 2, 3, 4 |

rgl.spheres(x, y, z, r, …) | Draws spheres with center (x, y, z) and radius r |

rgl.texts(x, y, z, text, …) | Adds text to the scene |

rgl.surface(x, y, z, …) | Adds surface to the current scene. x and y are two vectors defining the grid. z is a matrix defining the height of each grid point |

Functions | Description |
---|---|

rgl.clear(type = “shapes”) | Clears the scene from the specified stack (“shapes”, “lights”, “bboxdeco”, “background”) |

rgl.pop(type = “shapes”) | removes the last-added node from stack |

Functions | Description |
---|---|

rgl.viewpoint(theta, phi, fov, zoom, …) | Set the viewpoint orientation. theta and phi are the polar coordinates. fov is the field-of-view angle in degrees. zoom: zoom factor. |

rgl.light(theta, phi, …) | Adds a light source to the scene |

rgl.bg(…) | Setup the background environment of the scene |

rgl.bbox(…) | Setup the bounding box decoration |

Functions | Description |
---|---|

rgl.material(…) | Set material properties for geometry appearance |

aspect3d(x, y, z) | Set the aspect ratios of the x, y and z axes |

Functions | Description |
---|---|

rgl.snapshot(“plot.png”) | Saves a screenshot as png file |

rgl.postscript(“plot.pdf”,fmt=“pdf”) | Saves the screenshot to a file in ps, eps, tex, pdf, svg or pgf format. |

**Rgl.bringtotop()**: ‘rgl.bringtotop’ brings the current RGL window to the front of the window stack (and gives it focus).

**References**:

- Daniel Alder et al., RGL: A R-library for 3D visualization with OpenGL, http://rgl.neoscientists.org/arc/doc/RGL_INTERFACE03.pdf

This analysis has been performed using **R software** (ver. 3.1.2) and **rgl** (ver. 0.93.1098)

I recently posted an article describing how to make easily a 3D scatter plot in **R** using the package **scatterplot3d**.

This **R tutorial** describes how to perform an **interactive** **3d graphics** using **R software** and the function **scatter3d** from the package **car**.

The function **scatter3d()** uses the **rgl** package to draw and animate **3D scatter plots**.

The packages **rgl** and **car** are required for this tutorial:

`install.packages(c("rgl", "car"))`

Note that, on Linux operating system, the **rgl** package can be installed as follow:

**sudo apt-get install r-cran-rgl**

Load the packages:

`library("car")`

We’ll use the *iris* data set in the following examples :

```
data(iris)
head(iris)
```

```
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
```

```
sep.l <- iris$Sepal.Length
sep.w <- iris$Sepal.Width
pet.l <- iris$Petal.Length
```

*iris* data set gives the measurements of the variables sepal length and width, petal length and width, respectively, for 50 flowers from each of 3 species of iris. The species are Iris setosa, versicolor, and virginica.

The simplified formats are:

```
scatter3d(formula, data)
scatter3d(x, y, z)
```

**x, y, z**are respectively the coordinates of points to be plotted. The arguments*y*and*z*can be optional depending on the structure of*x*.**formula**: a model formula of form**y ~ x + z**. If you want to plot the points by groups, you can use**y ~ x + z | g**where*g*is a factor dividing the data into groups**data**: data frame within which to evaluate the formula

```
library(car)
# 3D plot with the regression plane
scatter3d(x = sep.l, y = pet.l, z = sep.w)
```

Note that, the plot can be manually rotated by holding down on the mouse or touchpad. It can be also zoomed using the scroll wheel on a mouse or pressing ctrl + using the touchpad on a PC or two fingers (up or down) on a mac.

**Change point colors and remove the regression surface**:

```
scatter3d(x = sep.l, y = pet.l, z = sep.w,
point.col = "blue", surface=FALSE)
```

`scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species)`

- To
**remove the grids**only, the argument**grid = FALSE**can be used as follow:

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
grid = FALSE)
```

Note that, the display of the surface(s) can be changed using the argument **fit**. Possible values for *fit* are **“linear”**, **“quadratic”**, **“smooth”** and **“additive”**

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
grid = FALSE, fit = "smooth")
```

**Remove surfaces**. The argument**surface = FALSE**is used.

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
grid = FALSE, surface = FALSE)
```

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
surface=FALSE, ellipsoid = TRUE)
```

Remove the grids from the ellipsoids:

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
surface=FALSE, grid = FALSE, ellipsoid = TRUE)
```

The argument **surface.col** is used. **surface.col** is a vector of colors for the regression planes.

For **multi-group plots**, the colors are used for the regression surfaces and for the points in the several groups.

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
surface=FALSE, grid = FALSE, ellipsoid = TRUE,
surface.col = c("#999999", "#E69F00", "#56B4E9"))
```

Read more about colors in R: colors in R

It’s also possible to use color palettes from the **RColorBrewer** package:

```
library("RColorBrewer")
colors <- brewer.pal(n=3, name="Dark2")
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
surface=FALSE, grid = FALSE, ellipsoid = TRUE,
surface.col = colors)
```

The arguments **xlab**, **ylab** and **zlab** are used:

```
scatter3d(x = sep.l, y = pet.l, z = sep.w,
point.col = "blue", surface=FALSE,
xlab = "Sepal Length (cm)", ylab = "Petal Length (cm)",
zlab = "Sepal Width (cm)")
```

**axis.scales = FALSE**

```
scatter3d(x = sep.l, y = pet.l, z = sep.w,
point.col = "blue", surface=FALSE,
axis.scales = FALSE)
```

By default, different colors are used for the 3 axes. The argument **axis.col** is used to specify colors for the 3 axes:

```
scatter3d(x = sep.l, y = pet.l, z = sep.w, groups = iris$Species,
surface=FALSE, grid = FALSE, ellipsoid = TRUE,
axis.col = c("black", "black", "black"))
```

The arguments below are used:

**labels**: text labels for the points, one for each point**id.n**: Number of relatively extreme points to identify automatically

```
scatter3d(x = sep.l, y = pet.l, z = sep.w,
surface=FALSE, labels = rownames(iris), id.n=nrow(iris))
```

The plot can be saved as png or pdf.

- The function
**rgl.snapshot()**is used to save the screenshot as png file:

`rgl.snapshot(filename = "plot.png")`

- The function
**rgl.postscript()**is used to saves the screenshot to a file in**ps, eps, tex, pdf, svg or pgf**format:

`rgl.postscript("plot.pdf",fmt="pdf")`

The function **Identify3d()**[ *car* package] allows to label points interactively with the mouse.

This analysis has been performed using **R software** (ver. 3.1.2) and **car** (ver. 2.0-25)

- Install and load scaterplot3d
- Prepare the data
- The function scatterplot3d()
- Basic 3D scatter plots
- Change the main title and axis labels
- Change the shape and the color of points
- Change point shapes by groups
- Change point colors by groups
- Change the global appearance of the graph
- Add bars
- Modification of scatterplot3d output
- Infos

There are many packages in R (*RGL*, *car*, *lattice*, *scatterplot3d*, …) for creating **3D graphics**.

This **tutorial** describes how to generate a scatter pot in the **3D space** using **R software** and the package **scatterplot3d**.

**scaterplot3d** is very simple to use and it can be easily extended by adding supplementary points or regression planes into an already generated graphic.

It can be easily installed, as it requires only an installed version of R.

```
install.packages("scatterplot3d") # Install
library("scatterplot3d") # load
```

The *iris* data set will be used:

```
data(iris)
head(iris)
```

```
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
```

*iris* data set gives the measurements of the variables sepal length and width, petal length and width, respectively, for 50 flowers from each of 3 species of iris. The species are Iris setosa, versicolor, and virginica.

A simplified format is:

`scatterplot3d(x, y=NULL, z=NULL)`

x, y, z are the coordinates of points to be plotted. The arguments *y* and *z* can be optional depending on the structure of *x*.

In what cases, *y* and *z* are optional variables?

**Case 1 : x is a formula**of type*zvar ~ xvar + yvar*. xvar, yvar and zvar are used as x, y and z variables**Case 2 : x is a matrix**containing at least 3 columns corresponding to x, y and z variables, respectively

```
# Basic 3d graphics
scatterplot3d(iris[,1:3])
```

```
# Change the angle of point view
scatterplot3d(iris[,1:3], angle = 55)
```

```
scatterplot3d(iris[,1:3],
main="3D Scatter Plot",
xlab = "Sepal Length (cm)",
ylab = "Sepal Width (cm)",
zlab = "Petal Length (cm)")
```

The argument *pch* and *color* can be used:

`scatterplot3d(iris[,1:3], pch = 16, color="steelblue")`

Read more on the different point shapes available in R : Point shapes in R

```
shapes = c(16, 17, 18)
shapes <- shapes[as.numeric(iris$Species)]
scatterplot3d(iris[,1:3], pch = shapes)
```

Read more on the different point shapes available in R : Point shapes in R

```
colors <- c("#999999", "#E69F00", "#56B4E9")
colors <- colors[as.numeric(iris$Species)]
scatterplot3d(iris[,1:3], pch = 16, color=colors)
```

Read more about colors in R: colors in R

The arguments below can be used:

**grid**: a logical value. If TRUE, a grid is drawn on the plot.**box**: a logical value. If TRUE, a box is drawn around the plot

```
scatterplot3d(iris[,1:3], pch = 16, color = colors,
grid=TRUE, box=FALSE)
```

Note that, the argument **grid = TRUE** plots only the grid on the xy plane. In the next section, we’ll see how to add grids on the other facets of the 3D scatter plot.

This section describes how to add *xy-*, *xz-* and *yz-* to **scatterplot3d** graphics.

We’ll use a custom function named **addgrids3d()**. The source code is available here : addgrids3d.r. The function is inspired from the discussion on this forum.

A simplified format of the function is:

```
addgrids3d(x, y=NULL, z=NULL, grid = TRUE,
col.grid = "grey", lty.grid=par("lty"))
```

**x, y, and z**are numeric vectors specifying the x, y, z coordinates of points. x can be a matrix or a data frame containing 3 columns corresponding to the x, y and z coordinates. In this case the arguments y and z are optional**grid**specifies the facet(s) of the plot on which grids should be drawn. Possible values are the combination of “xy”, “xz” or “yz”. Example: grid = c(“xy”, “yz”). The default value is TRUE to add grids only on xy facet.**col.grid, lty.grid**: the color and the line type to be used for grids

**Add grids on the different factes of scatterplot3d graphics**:

```
# 1. Source the function
source('http://www.sthda.com/sthda/RDoc/functions/addgrids3d.r')
# 2. 3D scatter plot
scatterplot3d(iris[, 1:3], pch = 16, grid=FALSE, box=FALSE)
# 3. Add grids
addgrids3d(iris[, 1:3], grid = c("xy", "xz", "yz"))
```

The problem on the above plot is that the grids are drawn over the points.

The **R code** below, we’ll put the points in the foreground using the following steps:

- An empty scatterplot3 graphic is created and the result of
**scatterplot3d()**is assigned to*s3d* - The function
**addgrids3d()**is used to add grids - Finally, the function
**s3d$points3d**is used to add points on the 3D scatter plot

```
# 1. Source the function
source('~/hubiC/Documents/R/function/addgrids3d.r')
# 2. Empty 3D scatter plot using pch=""
s3d <- scatterplot3d(iris[, 1:3], pch = "", grid=FALSE, box=FALSE)
# 3. Add grids
addgrids3d(iris[, 1:3], grid = c("xy", "xz", "yz"))
# 4. Add points
s3d$points3d(iris[, 1:3], pch = 16)
```

The function **points3d()** is described in the next sections.

The argument **type = “h”** is used. This is useful to see very clearly the x-y location of points.

```
scatterplot3d(iris[,1:3], pch = 16, type="h",
color=colors)
```

**scatterplot3d** returns a list of function closures which can be used to add elements on a existing plot.

The returned functions are :

**xyz.convert()**: to convert 3D coordinates to the 2D parallel projection of the existing scatterplot3d. It can be used to add arbitrary elements, such as legend, into the plot.**points3d()**: to add*points*or*lines*into the existing plot**plane3d()**: to add a*plane*into the existing plot**box3d()**: to add or refresh a*box*around the plot

- The result of
**scatterplot3d()**is assigned to*s3d* - The function
**s3d$xyz.convert()**is used to specify the coordinates for legends - the function
**legend()**is used to add legends to plots

```
s3d <- scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend(s3d$xyz.convert(7.5, 3, 4.5), legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16)
```

It’s also possible to specify the **position of legends** using the following keywords: “bottomright”, “bottom”, “bottomleft”, “left”, “topleft”, “top”, “topright”, “right” and “center”.

```
# "right" position
s3d <- scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("right", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16)
```

```
# Use the argument inset
s3d <- scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("right", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16, inset = 0.1)
```

What means the argument **inset** in the R code above?

The argument **inset** is used to inset distance(s) from the margins as a fraction of the plot region when legend is positioned by keyword. ( see ?legend from R). You can play with *inset* argument using negative or positive values.

```
# "bottom" position
s3d <- scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("bottom", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16)
```

Using *keywords* to specify the *legend position* is very simple. However, sometimes, there is an overlap between some points and the legend box or between the axis and legend box.

Is there any solution to avoid this overlap?

Yes, there are several solutions using the combination of the following arguments for the function **legend()**:

**bty = “n”**: to**remove the box around the legend**. In this case the background color of the legend becomes transparent and the overlapping points become visible.**bg = “transparent”**: to change the background color of the legend box to*transparent*color (this is only possible when bty != “n”).**inset**: to modify the distance(s) between plot margins and the legend box.**horiz**: a logical value; if TRUE, set the legend horizontally rather than vertically**xpd**: a logical value; if TRUE, it enables the legend items to be drawn outside the plot.

```
# Custom point shapes
s3d <- scatterplot3d(iris[,1:3], pch = shapes)
legend("bottom", legend = levels(iris$Species),
pch = c(16, 17, 18),
inset = -0.25, xpd = TRUE, horiz = TRUE)
```

```
# Custom colors
s3d <- scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("bottom", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16,
inset = -0.25, xpd = TRUE, horiz = TRUE)
```

```
# Custom shapes/colors
s3d <- scatterplot3d(iris[,1:3], pch = shapes, color=colors)
legend("bottom", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"),
pch = c(16, 17, 18),
inset = -0.25, xpd = TRUE, horiz = TRUE)
```

In the R code above, you can play with the arguments *inset*, *xpd* and *horiz* to see the effects on the appearance of the legend box.

The function **text()** is used as follow:

```
scatterplot3d(iris[,1:3], pch = 16, color=colors)
text(s3d$xyz.convert(iris[, 1:3]), labels = rownames(iris),
cex= 0.7, col = "steelblue")
```

- The result of
**scatterplot3d()**is assigned to*s3d* - A linear model is calculated as follow : lm(zvar ~ xvar + yvar). Assumption : zvar depends on xvar and yvar
- The function
**s3d$plane3d()**is used to add the regression plane - Supplementary points are added using the function
**s3d$points3d()**

The data sets *trees* will be used:

```
data(trees)
head(trees)
```

```
Girth Height Volume
1 8.3 70 10.3
2 8.6 65 10.3
3 8.8 63 10.2
4 10.5 72 16.4
5 10.7 81 18.8
6 10.8 83 19.7
```

This data set provides measurements of the girth, height and volume for black cherry trees.

**3D scatter plot with the regression plane**:

```
# 3D scatter plot
s3d <- scatterplot3d(trees, type = "h", color = "blue",
angle=55, pch = 16)
# Add regression plane
my.lm <- lm(trees$Volume ~ trees$Girth + trees$Height)
s3d$plane3d(my.lm)
# Add supplementary points
s3d$points3d(seq(10, 20, 2), seq(85, 60, -5), seq(60, 10, -10),
col = "red", type = "h", pch = 8)
```

This analysis has been performed using **R software** (ver. 3.1.2) and **scatterplot3d** (ver. 0.3-35)