Museums and Nature Centers

There are thousands of museums, aquariums, and zoos across the United States. In this project, we’ll take a look at the distribution of these institutions by geographic region, type, and revenue.

Our data is compiled from administrative records from the Institute of Museum and Library Services, IRS records, and private foundation grantmaking records. This data reflects the status of each institution as of 2013. For each institution, we have information on its name, type, and location. Each institution also has a parent organization – for example, if a museum housed at a university, its parent organization is the university where it resides. Financial data on annual revenue is available at the parent organization level.

library(dplyr)
library(ggplot2)
library(stringr)
library(tidyr)
library(plotrix)

Data Exploration

Let’s start by loading our dataset. We’ve provided a file named museums.csv. Load this file into a data frame named museums_df.

# Load file as data frame
museums_df <- read.csv("museums.csv")

Take a look at the head of this data frame. Make sure to click through using the arrows in the header to see all the available columns.

# Inspect data frame
head(museums_df)

The Museum.Name column represents the name of each individual institution, while the Legal.Name column represents the name of each institution’s parent entity. For example, if “Codecademy University” has two museums on campus, each of those museums would have their own names under Museum.Name and both would share the same Legal.Name i.e. “Codecademy University”.

Museums by Type

In this section, we’ll explore the distribution of institutions in our dataset by type. Our data frame contains a column called Museum.Type describing what kind of museum each location is – a history museum, a zoo, an aquarium, etc. Create and print a bar plot called museum_type that maps Museum.Type to the x axis and counts the frequency of each type on the y axis. Which category is most common?

museum_type <- ggplot(museums_df, aes(x = Museum.Type)) + 
geom_bar()
museum_type

The plot we just created is hard to read because our categories are long. Add a scale_x_discrete() layer to customize our x axis, using the function scales::wrap_format(8) to reformat our labels.

wrap_format() is a function from the scales packages which comes included with ggplot2. By setting the value of wrap_format() to 8, we are telling it that the maximum width per line should be no more than 8 characters.

Now we should be able to see which category is most common.

# Create and print bar plot by type
museum_type <- ggplot(museums_df, aes(x = Museum.Type)) + 
geom_bar() +
scale_x_discrete(labels = scales::wrap_format(8))
museum_type

We’ve included a boolean (TRUE or FALSE) column in our data frame called Is.Museum. The TRUE category includes typical museums like art, history, and science museums. The FALSE category includes zoos, aquariums, nature preserves, and historic sites, which are included in this data but aren’t what most people think of when they hear the word “museum.”

Create a new bar plot called museum_class, mapping Is.Museum to the x axis. Since “TRUE” and “FALSE” aren’t very descriptive, use scale_x_discrete() to rename the x axis labels to more easily understood terms – for example, “Museum” vs “Non-Museum”.

# Create and print bar plot by museum vs non-museum
museum_class <- ggplot(museums_df, aes(x = Is.Museum)) +
geom_bar() +
scale_x_discrete(labels = c( "TRUE" = "Museum", "FALSE" = "Non-Museum"))
museum_class

Instead of looking at the distribution across the entire United States, maybe we’re just interested in a few states. Filter museums_df to include a few states you might be interested in, using the State..Administrative.Location. column; for example, we can choose IL, CA, and NY. Call this filtered data frame museums_states.

After creating museums_states, recreate our bar plot showing the distribution of museums vs non-museums and use facet_grid() to display each state’s distribution in a separate panel. Call this plot museum_facet. How does the distribution of museum vs non-museum vary across the states you chose?

# Filter data frame to select states
museums_states <- museums_df %>%
  filter(State..Administrative.Location. %in% c("IL", "CA", "NY"))


# Create and print bar plot with facets
museums_states <- ggplot(museums_states, aes(x = Is.Museum)) +
geom_bar() +
scale_x_discrete(labels = c( "TRUE" = "Museum", "FALSE" = "Non-Museum")) +
  facet_grid(cols = vars(State..Administrative.Location.))
museums_states

Illinois has less number of Museums than any of the 3 states, and NY the most. California has more museums that “non-museums” as opposed to the other 2 states.

Our data also contains information on each museum’s region, representing groups of states. Create a stacked bar plot using museums_df showing the count of museums by region (Region.Code..AAM.), mapping Is.Museum to the fill aesthetic. Convert Region.Code..AAM. to a factor (e.g. factor(Region.Code..AAM.)) so ggplot2 plots its levels as discrete rather than continuous values. Call this plot museum_stacked.

museum_stacked <- ggplot(museums_df, aes(x = factor(Region.Code..AAM.), fill = Is.Museum)) +
  geom_bar() 
museum_stacked

Our plot is hard to read – right now, we don’t know what the region numbers correspond to. Use scale_x_discrete() to rename the numeric labels to text according to the following table.

Code	Region
1	New England
2	Mid-Atlantic
3	Southeastern
4	Midwest
5	Mountain Plains
6	Western

Similarly, add a scale_fill_discrete() layer to relabel the “TRUE” and “FALSE” labels in our legend to “Museum” and “Non-Museum”.

Based on the plot we created, which region has the most museums?

# Create and print stacked bar plot
museum_stacked <- ggplot(museums_df, aes(x = factor(Region.Code..AAM.), fill = Is.Museum)) +
  geom_bar() +
  scale_x_discrete(labels = c( "1" = "New England",
                               "2" = "Mid-Atlantic",
                               "3" = "Southeastern",
                               "4" = "Midwest",
                               "5" = "Montain Plains",
                               "6" = "Western")) +
  scale_fill_discrete(labels = c( "TRUE" = "Museum", "FALSE" = "Non-Museum"))
museum_stacked

Midwest has the most museums.

Rather than seeing counts, perhaps we’re more interested in the percentage of museums vs non-museums by region. Transform the plot we just created to a stacked bar plot showing values out of 100% by passing position = “fill” to our geom_bar() layer. Apply the scales::percent_format() function to transform our y axis labels into percentage values.

How does the distribution of museum types vary by region?

# Create and print stacked bar plot
museum_stacked <- ggplot(museums_df, aes(x = factor(Region.Code..AAM.), fill = Is.Museum)) +
  geom_bar(position = "fill") +
  scale_x_discrete(labels = c( "1" = "New England",
                               "2" = "Mid-Atlantic",
                               "3" = "Southeastern",
                               "4" = "Midwest",
                               "5" = "Montain Plains",
                               "6" = "Western")) +
  scale_fill_discrete(labels = c( "TRUE" = "Museum", "FALSE" = "Non-Museum"), ) +
  scale_y_continuous(labels = scales::percent_format())
  
museum_stacked

Midwest has less percentage of museums. Southeastern, Montain Plains and Western has a similar distribution. New England and Mid-Atlantic are in between.

Our graph looks pretty good! However, our axes titles are a little non-descript. Using the labs() layer, let’s title this plot “Museum Types by Region”, relabel the x axis title as “Region”, relabel the y axis title as “Percentage of Total”, and relabel the fill legend title as “Type”.

museum_stacked <- ggplot(museums_df, aes(x = factor(Region.Code..AAM.), fill = Is.Museum)) +
  geom_bar(position = "fill") +
  scale_x_discrete(labels = c( "1" = "New England",
                               "2" = "Mid-Atlantic",
                               "3" = "Southeastern",
                               "4" = "Midwest",
                               "5" = "Montain Plains",
                               "6" = "Western")) +
  scale_fill_discrete(labels = c( "TRUE" = "Museum", "FALSE" = "Non-Museum"), ) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title = "Museum Types by Region", x = "Region", y = "Percentage of Total", fill = "Type")
  
museum_stacked

Museums by Revenue

For the next few tasks, we’ll switch to looking at how much money each institution brought in and how that varies across geographies. Because we only have revenue data at the parent organization level, we’ll want to first filter our dataset to omit any duplicates. Next, we’ll create a few data frames from our starting data to look at different groups of museums by how much money they bring in.

Create a new data frame called museums_revenue_df that retains only unique values of Legal.Name in museums_df. Additionally, filter this data frame to include only entities with Annual.Revenue greater than 0.

# Filter data frame
museums_revenue_df <- museums_df %>%
  distinct(Legal.Name, .keep_all = TRUE) %>% 
  filter(Annual.Revenue > 0)

Create a second data frame from museums_revenue_df (the first data frame we created in this task) called museums_small_df that retains only museums with Annual.Revenue less than $1,000,000.

# Filter for only small museums
museums_small_df <- museums_revenue_df %>% 
  filter(Annual.Revenue < 1000000)

Create a third data frame from museums_revenue_df (the first data frame we created in this task) called museums_large_df that retains only museums with Annual.Revenue greater than $1,000,000,000.

# Filter for only large museums
museums_large_df <- museums_revenue_df %>% 
  filter(Annual.Revenue >= 1000000)

Let’s start by visualizing the distribution of annual revenue for our small museums dataset. Create a histogram called revenue_histogram using museums_small_df with Annual.Revenue mapped to the x axis. Experiment with different binwidth values to see what works best for our data, considering that our x axis variable ranges from 0 to $1,000,000.

# Create and print histogram
revenue_histogram <- ggplot(museums_small_df, aes(x = Annual.Revenue)) +
  geom_histogram(binwidth = 25000)
revenue_histogram

Our x axis is a little hard to read. Let’s make it more clear! Add a scale_x_continuous() layer applying the function scales::dollar_format() to our x axis labels. dollar_format() is a function from the scales library included in ggplot2 that adds dollar signs and commas to monetary data.

revenue_histogram +
  scale_x_continuous(labels = scales::dollar_format())

Now, let’s look at the variation in revenue for large museums by region. Create a boxplot called revenue_boxplot using museums_large_df, mapping Region.Code..AAM. to the x axis and Annual.Revenue to the y axis. Remember to convert Region.Code..AAM. to a factor (e.g. factor(Region.Code..AAM.)) so ggplot2 plots its levels as discrete rather than continuous values. Use scale_x_discrete() to rename the numeric region codes to their text equivalents.

revenue_boxplot  <- ggplot(museums_large_df, aes(x = factor(Region.Code..AAM.), y = Annual.Revenue)) +
  geom_boxplot() +
  scale_x_discrete(labels = c( "1" = "New England",
                               "2" = "Mid-Atlantic",
                               "3" = "Southeastern",
                               "4" = "Midwest",
                               "5" = "Montain Plains",
                               "6" = "Western"))
revenue_boxplot

The plot we just created is a little hard to read, since there’s one outlier so far above all the other data points. This one museum made a lot of money in 2013! Let’s zoom in so we can see the rest of our boxes more clearly. Add a coord_cartesian() layer setting ylim to c(1e9, 3e10). This tells our plot to zoom in on the y axis range between $1,000,000,000 and $30,000,000,000. How do the median, 75th, and 25th quantiles vary by region?

revenue_boxplot  <- revenue_boxplot +
  coord_cartesian(ylim = c(1e9, 3e10))
revenue_boxplot

Though we can see the distribution by region more clearly now, our y axis label is hard to understand. Let’s reformat our y axis as billions of dollars. The function defined below will convert values like $1,000,000,000 to $1B, which is much easier to read. Add a scale_y_continuous() layer to map our y axis labels using this

function(x) paste0("$", x/1e9, "B")

We’ve made our box plot much clearer to read.

revenue_boxplot +
  scale_y_continuous(labels = function(x) paste0("$", x/1e9, "B"))

Now, let’s take a look at revenue across all museums in our dataset. Using museums_revenue_df, create a bar plot called revenue_barplot mapping Region.Code..AAM. to the x axis and Annual.Revenue to the y axis. Remember to transform Region.Code..AAM. to a factor (factor(Region.Code..AAM.)) so it shows up as discrete values.

Use stat = “summary” and fun = “mean” to calculate and display the mean revenue by region. Apply the appropriate x and y axis label transformations to make our labels more clear. Which region has the highest and lowest mean revenues?

revenue_barplot <- ggplot(museums_revenue_df, aes(x = factor(Region.Code..AAM.), y = Annual.Revenue)) +
  geom_bar(stat = "summary", fun = "mean") +
  scale_x_discrete(labels = c( "1" = "New England",
                               "2" = "Mid-Atlantic",
                               "3" = "Southeastern",
                               "4" = "Midwest",
                               "5" = "Montain Plains",
                               "6" = "Western")) +
  scale_y_continuous(labels = function(x) paste0("$", x/1e6, "M"))
revenue_barplot

Once again, use the labs layer to make our plot more clear. Title the plot “Mean Annual Revenue by Region”, relabel the y axis title to “Mean Annual Revenue”, and relabel the x axis title to “Region”.

revenue_barplot +
  labs(title = "Mean Annual Revenue by Region", x = "Region", y = "Mean Annual Revenue")

Finally, let’s add some error bars to our means. We’ll need to calculate our standard errors before creating our plot. You can do so using the following code:

museums_error_df <- museums_revenue_df %>%
  group_by(Region.Code..AAM.) %>%
  summarize(
    Mean.Revenue = mean(Annual.Revenue), 
    Mean.SE = std.error(Annual.Revenue)) %>%
  mutate(
    SE.Min = Mean.Revenue - Mean.SE, 
    SE.Max = Mean.Revenue + Mean.SE)

Add error bars to our mean revenue by geography bar plot using the geom_errorbar() layer. Call this new plot revenue_errorbar. Our new y variable is our calculated Mean.Revenue column. Make sure to change the stat being used, since we’re now displaying our calculated means as is rather than calculating them as we create the plot. Which regions have more or less variability around their mean revenues?

museums_error_df <- museums_revenue_df %>%
  group_by(Region.Code..AAM.) %>%
  summarize(
    Mean.Revenue = mean(Annual.Revenue), 
    Mean.SE = std.error(Annual.Revenue)) %>%
  mutate(
    SE.Min = Mean.Revenue - Mean.SE, 
    SE.Max = Mean.Revenue + Mean.SE)

`summarise()` ungrouping output (override with `.groups` argument)

revenue_errorbar <- ggplot(museums_error_df, aes(x = factor(Region.Code..AAM.), y = Mean.Revenue)) +
  geom_bar(stat = "identity") +
  scale_x_discrete(labels = c( "1" = "New England",
                               "2" = "Mid-Atlantic",
                               "3" = "Southeastern",
                               "4" = "Midwest",
                               "5" = "Montain Plains",
                               "6" = "Western")) +
  scale_y_continuous(labels = function(x) paste0("$", x/1e6, "M")) +
  labs(title = "Mean Annual Revenue by Region", x = "Region", y = "Mean Annual Revenue") +
  geom_errorbar(aes(ymin = SE.Min, ymax = SE.Max), width = 0.2)
revenue_errorbar