3-D Plots and Maps
3-D scatter plots using Swiss Fertility and Socioeconomic Indicators (1888) Data
swiss data includes 47 observations on 6 variables. All variables but ‘Fertility’ give proportions of the population.
- Fertility Ig, ‘common standardized fertility measure’
- Agriculture % of males involved in agriculture as occupation
- Examination % draftees receiving highest mark on army examination
- Education % education beyond primary school for draftees.
- Catholic % ‘catholic’ (as opposed to ‘protestant’).
- Infant.Mortality live births who live less than 1 year.
interleave <- function(v1, v2) as.vector(rbind(v1,v2))
plot3d(swiss$Education, swiss$Agriculture, swiss$Fertility,
xlab = "Education", ylab = "Agriculture", zlab = "Fertility",
type = "s", size = 0.75, lit = FALSE)
segments3d(interleave(swiss$Education, swiss$Education),
interleave(swiss$Agriculture, swiss$Agriculture),
interleave(swiss$Fertility, min(swiss$Fertility)),
alpha = 0.7, col = "sienna")
- Change the appearance of the background and the axes
# Make plot without axis ticks or labels
plot3d(swiss$Education, swiss$Agriculture, swiss$Fertility,
xlab = "", ylab = "", zlab = "",
axes = FALSE,
size = .75, type = "s", lit = FALSE)
segments3d(interleave(swiss$Education, swiss$Education),
interleave(swiss$Agriculture, swiss$Agriculture),
interleave(swiss$Fertility, min(swiss$Fertility)),
alpha = 0.7, col = "sienna")
# Draw the box.
rgl.bbox(color = "grey50",
emission = "grey50",
xlen = 0, ylen = 0, zlen = 0)
# Set default color of future objects
rgl.material(color = "black")
axes3d(edges = c("x--", "y+-", "z--"),
ntick = 6,
cex = .75)
mtext3d("Education", edge = "x--", line = 2)
mtext3d("Agriculture", edge = "y+-", line = 3)
mtext3d("Fertility", edge = "z--", line = 3)
- We can add a surface of predicted value to a three-dimensional scatter plot:
# Defaults to range of x and y variables, and a 16x16 grid
predictgrid <- function(model, xvar, yvar, zvar, res = 16, type = NULL) {
xrange <- range(model$model[[xvar]])
yrange <- range(model$model[[yvar]])
newdata <- expand.grid(x = seq(xrange[1], xrange[2], length.out = res),
y = seq(yrange[1], yrange[2], length.out = res))
names(newdata) <- c(xvar, yvar)
newdata[[zvar]] <- predict(model, newdata = newdata, type = type)
newdata
}
# Convert long-style data frame with x, y, and z vars into a list
# with x and y as row/column values, and z as a matrix.
df2mat <- function(p, xvar = NULL, yvar = NULL, zvar = NULL) {
if (is.null(xvar)) xvar <- names(p)[1]
if (is.null(yvar)) yvar <- names(p)[2]
if (is.null(zvar)) zvar <- names(p)[3]
x <- unique(p[[xvar]])
y <- unique(p[[yvar]])
z <- matrix(p[[zvar]], nrow = length(y), ncol = length(x))
m <- list(x, y, z)
names(m) <- c(xvar, yvar, zvar)
m
}
# Function to interleave the elements of two vectors
interleave <- function(v1, v2) as.vector(rbind(v1,v2))
# Make a copy of the data set
m <- swiss
# Get the predicted values
mod <- lm(Fertility ~ Education*Agriculture, data = m)
m$pred_Fertility <- predict(mod)
grid_df <- predictgrid(mod, "Education", "Agriculture", "Fertility")
grid_list <- df2mat(grid_df)
plot3d(m$Education, m$Agriculture, m$pred_Fertility,
xlab = "Education", ylab = "Agriculture", zlab = "Predicted Fertility",
type = "s", size = 0.75, lit = FALSE)
spheres3d(m$Education, m$Agriculture, m$pred_Fertility, alpha = 0.4, type = "s", size = 0.75, lit = FALSE)
#errord
segments3d(interleave(m$Education,m$Education),
interleave(m$Agriculture,m$Agriculture),
interleave(m$Fertility,m$pred_Fertility),
alpha = 0.4, col = "blue")
# Mesh of predicted values
surface3d(grid_list$Education, grid_list$Agriculture, grid_list$Fertility,
alpha = 0.4, front = "lines", back = "lines")
Creating maps
Some Europe countries
# Get map data for world
world_map <- map_data("world")
europe <- map_data("world", region = c('UK', 'France', 'Spain', 'Sweden', 'Belgium', 'Italy', 'Romania', 'Belarus', 'Switzerland', 'Germany', 'Denmark', 'Poland', 'Ukraine', 'Hungary', 'Slovakia', 'Finland'))
ggplot(europe, aes(x = long, y = lat, group = group, fill = region)) +
geom_polygon(colour = "black") +
scale_fill_hue(l=40) +
coord_equal() +
theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
scale_y_continuous(breaks = NULL) +
scale_x_continuous(breaks = NULL)
US maps with 2018 population information
- We can create a map with regions that are colored according to variable values (Choropleth Map)
head(US_pop)# Transform the data set to the correct format
US_pop <- data.frame(region = tolower(US_pop$State), Population_2018=US_pop$Population_2018)
states_map <- map_data("state")
# Merge the data sets together
pop_map <- merge(states_map, US_pop, by = "region")
pal <- wesanderson::wes_palette("BottleRocket1", 20, type = "continuous")
# sort the data.
pop_map <- arrange(pop_map, group, order)
ggplot(pop_map, aes(x = long, y = lat, group = group, fill = Population_2018)) +
geom_polygon(colour = "gray") +
expand_limits(x = states_map$long, y = states_map$lat) +
coord_map("lagrange") +
labs(title = "US State Populations in 2018",fill="") +
scale_fill_gradientn(colours = pal)
# Find the quantile bounds
qa <- quantile(US_pop$Population_2018, c(0, 0.2, 0.4, 0.6, 0.8, 1.0))
qa## 0% 20% 40% 60% 80% 100%
## 573720 1350575 3159345 5684203 9032872 39776830# Add a column of the quantile category
US_pop$Pop_q <- cut(US_pop$Population_2018, qa,
labels = c("0-20%", "20-40%", "40-60%", "60-80%", "80-100%"),
include.lowest = TRUE)
# Generate a discrete color palette with 5 values
pal <- colorRampPalette(c("#FE9929", "grey90", "#D95F0E"))(5)
ggplot(US_pop, aes(map_id = region, fill = Pop_q)) +
geom_map(map = states_map, colour = "black") +
scale_fill_manual(values = pal) +
expand_limits(x = states_map$long, y = states_map$lat) +
coord_map("lagrange") +
labs(fill = "US State\nPopulations\n2018\nPercentile") +
theme_void()