This is a tilegram of Russian regions with area being distorted according to 2015 population. I relied on Gastner and Newman, 2004 distortion algorithm. Each hexagon is approximately equal to 25,000 people. Hover over the regions to see their names and population counts (in Russian).

# (c) Dmitriy Skougarevskiy, November 2016
#
# Use of this source code is governed by the MIT license
# located at https://opensource.org/licenses/MIT

# Load dependencies
packages <- c("maptools", "raster", "sp", "rgeos", "devtools")
if (length(setdiff(packages, rownames(installed.packages()))) > 0) {
	install.packages(setdiff(packages, rownames(installed.packages())))  
}
lapply(packages, library, character.only = TRUE)
# NB: uncomment the below two packages lines to install
# the packages from Github. They require http://www.fftw.org
# (e.g. brew install fftw)
#devtools::install_github("omegahat/Rcartogram")  
#devtools::install_github("chrisbrunsdon/getcartr", subdir = "getcartr")
library(Rcartogram)
library(getcartr)

# Function to download and open zipped shapefile
# via http://thebiobucket.blogspot.ch/2013/09/batch-downloading-zipped-shapefiles.html
url_shp_to_spdf <- function(URL) {

	require(maptools)

	wd <- getwd()
	td <- tempdir()
	setwd(td)

	temp <- tempfile(fileext = ".zip")
	download.file(URL, temp)
	unzip(temp)

	shp <- dir(tempdir(), "*.shp$")
	y <- readShapePoly(shp)

	unlink(dir(td))
	setwd(wd)
	return(y)
}

# Download Russian regional boundaries shapefile
# that has 2015 regional population embedded therein.
# I aggressively simplified it with mapshaper and by hand
# (removal of the Far North islands, all small islands or lakes
# not exceeding than Sebastopol in area)
rus_bounds <- url_shp_to_spdf("http://donutholes.ch/test/rus_regions_simpl.shp.zip")

# Use Albers Siberia projection
WGS_84 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
Albers_Siberia <- CRS("+proj=aea +lat_1=52 +lat_2=64 +lat_0=0 +lon_0=105 +x_0=18500000 +y_0=0 +ellps=krass +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12 +units=m +no_defs")

proj4string(rus_bounds) <- WGS_84
rus_bounds <- spTransform(rus_bounds, Albers_Siberia)

# Create the cartogram (takes time and CPU, decrease the resolution
# to 512 to enjoy reasonable computation time at the cost of larger error in shape areas)
rus_bounds_transf <- quick.carto(rus_bounds, rus_bounds$region_pop, res = 3000, gapdens = 2, blur = 0)

# A hack to fix possible topology errors
# (via http://gis.stackexchange.com/a/163480)
rus_bounds_transf <- gBuffer(rus_bounds_transf, byid = TRUE, width = -0.0001)

# Export transformed shape now if you feel ready to explore it
#writePolyShape(rus_bounds_transf, "rus_bounds_transf")

proj4string(rus_bounds_transf) <- Albers_Siberia

# Check whether areas are now proportional to population
## Correlation between region population and area now exceeds 0.9:
## cor(rus_bounds_transf@data$region_pop, gArea(rus_bounds_transf, byid = T))
## Difference in area between Crimea and Moscow now almost corresponds to difference in population:
## gArea(rus_bounds_transf, byid = T)[36]/gArea(rus_bounds_transf, byid = T)[54]
## rus_bounds_transf@data$region_pop[36]/rus_bounds_transf@data$region_pop[54]

# Also examine the residuals of the transformation
## discrep <- cart.resids(rus_bounds_transf, rus_bounds$region_pop)
## shs <- auto.shading(c(discrep,-discrep), n = 5, col = brewer.pal(5, "RdBu"))
## choropleth(rus_bounds_transf, discrep, shading = shs)
## choro.legend(600, 600, shs)

# Establish the average number of people per unit of transformed shapefile area
pop_per_area <- mean(rus_bounds_transf@data$region_pop/gArea(rus_bounds_transf, byid = T))
# Set area for hexagons so that each one encompasses 25000 people
hex_area <- 25000/pop_per_area

# Function to produce hexagonal grid
# (via http://strimas.com/spatial/hexagonal-grids/)
make_grid <- function(x, cell_diameter, cell_area, clip = FALSE) {
	if (missing(cell_diameter)) {
		if (missing(cell_area)) {
			stop("Must provide cell_diameter or cell_area")
		} else {
			cell_diameter <- sqrt(2 * cell_area / sqrt(3))
		}
 	}
	ext <- as(extent(x) + cell_diameter, "SpatialPolygons")
	projection(ext) <- projection(x)
	# generate array of hexagon centers
	g <- spsample(ext, type = "hexagonal", cellsize = cell_diameter, 
					offset = c(0.5, 0.5))
	# convert center points to hexagons
	g <- HexPoints2SpatialPolygons(g, dx = cell_diameter)
	# clip to boundary of study area
	if (clip) {
		g <- gIntersection(g, x, byid = T)
	} else {
		g <- g[x, ]
	}
	# clean up feature IDs
	row.names(g) <- as.character(1:length(g))
	return(g)
}

# Create the hexagonal grid
source_hex_grid <- make_grid(rus_bounds_transf, cell_area = hex_area, clip = F)

# Add attributes to the grid to produce the final object
rus_hex <- SpatialPolygonsDataFrame(source_hex_grid, over(source_hex_grid, rus_bounds_transf))

proj4string(rus_hex) <- Albers_Siberia

# Define the color palette
region_colors <- data.frame(REGION = levels(rus_hex$REGION),
							color = brewer.pal(12, name = "Set3")[sample(1:12, 85, replace = T)]
							)

# Add color palette to the grid
rus_hex$color <- as.character(region_colors[match(rus_hex$REGION, region_colors$REGION), "color"])

# Export the resulting shapefile with hexagons
writePolyShape(rus_hex, "rus_hex")

# Convert shapefile to TopoJSON with
# https://github.com/topojson/topojson-1.x-api-reference/blob/master/Command-Line-Reference.md
system(paste0("cd ", getwd(), "; ",
		"topojson --width 960 --height 550 --margin 10 -s .25 --shapefile-encoding utf8 --id-property SP_ID -p oktmo=OKTMO -p region=REGION -p population=region_pop -p color=color -o rus_hex.topojson rus_hex.shp"
		)
	)