Select explanatory variables and outcome variables

Define two lists of variables, one for institutional factors and another for outcomes of interest, and create smaller dataset for our purposes

# pick the list of variables we want to use from the QoG
fraser <- c("fi_ftradeint_pd", "fi_legprop_pd", "fi_reg_pd", "fi_sm_pd", "fi_sog_pd")

outcomes <- c(
  # GDP per capita, growth, Gini, life expectancy
  "unna_gdppc", "wdi_gdpcapgr", "wdi_gini", "wdi_lifexp",

  # fractionalization
  "al_ethnic", "al_religion",

  # education female and male ages 35-44
  "gea_ea3544f", "gea_ea3544m",

  # corruption, public sector and political
  "vdem_pubcorr", "vdem_corr")

# prepare smaller work dataset containing institutional factors
institutions <- qog %>% 
  select(cname, ccodealp, fraser, outcomes)

# remove rows with missing data
institutions <- na.omit(institutions)

# use row names to identify the countries
row.names(institutions) <- institutions$cname

qog_summary %>% 
  filter(Variable %in% c(fraser, outcomes)) %>% 
  knitr::kable()

Scaling the data to account for differing units of measure

In general, before doing the cluster analysis, you need to make sure that all your variables are measured in comparable units – otherwise your analysis can be heavily distorted by just one variable that happens to have much larger values than the others.

In the present case, I’m using the Fraser Institute data on economic institutions, and this data is already scaled – all variables are between 1 and 10. But if you add some other instititutional factors you need to consider this possible problem.

The simplest way to standardize the dataset is to use the scale function:

fraser_scaled   <- institutions %>% select(fraser) %>% scale() %>% as.data.frame()
outcomes_scaled <- institutions %>% select(outcomes) %>% scale() %>% as.data.frame()

The basic analysis

The cluster analysis identifies how similar two countries are to each other, either based on Euclidean distance or based on correlations.

If our case, the Euclidean distance between, say, USA and France is:

\[ \begin{align} d_{USA-FR}^2 =& [(Size of Govt)_{USA}-(Size of Govt)_{FR}]^2 + \\ & [(Deregulation)_{USA}-(Deregulation)_{FR}]^2 + \\ & [(Free Trade)_{USA}-(Free Trade)_{FR}]^2 + \\ & [(Sound Money)_{USA}-(Sound Money)_{FR}]^2 + \\ & [(Property Rights)_{USA}-(Property Rights)_{FR}]^2 \end{align} \] The “ward.D” method of clustering is calculating these distances between all pairs of countries and creates the clusters such that countries within a cluster are those that have the smallest distances between each other. This is usually the best method for identifying clearly distinct categories. See this link for more information about clustering methods.

Create a hierarchical cluster tree based on Euclidean distance between cases.

hc <- institutions %>% 
  select(fraser) %>%                # choose only fraser variables
  dist(method = "euclidean") %>%    # calculate distances between pairs of countries
  hclust(method="ward.D")           # build cluster tree

Plot the cluster tree (known as a “dendrogram”)

plot(hc, 
     hang = -1,   # arranges labels at fixed distance
     cex = 0.7,   # size of labels
     main = ""
     )

Plot a circular dendrogram

hc %>% 
  ape::as.phylo() %>% 
  plot(type="fan", cex = 0.8)

Based on the plot we can see that there we can choose to separate the countries into 2, 3 or 4 categories. I’m going to choose 3 here.

Categorize the data based on the clusters

Split the set into 3 groups and add the group number to the dataset.

institutions <- institutions %>% 
  mutate(groups_fraser = cutree(hc, 3)) %>% 
  
  # recode the group numbers to be in a better order
  # (this was discovered by first doing the analysis without recoding)
  mutate(groups_fraser = recode(groups_fraser, `1` = 2, `2` = 1, `3` = 3))

List countries in each group

map_df(1:3, function(g){
  institutions %>% 
    filter(groups_fraser == g) %>% 
    select(cname) %>% 
    rename_at("cname", funs(paste("Group", g)))
  }) %>% 
  
  # shift elements in Group 2 and 3 upwards in the table to replace the NAs
  mutate(`Group 2` = lead(`Group 2`, length(na.omit(`Group 1`)))) %>% 
  mutate(`Group 3` = lead(`Group 3`, length(na.omit(`Group 1`)) + 
                                     length(na.omit(`Group 2`))
                          )
         ) %>% 
  
  # eliminate rowns that have only NAs
  filter(!is.na(`Group 1`) | !is.na(`Group 2`) | !is.na(`Group 3`)) %>% 
  
  # change NAs to ""
  mutate_all(funs(ifelse(is.na(.), "", .))) %>% 
  
  knitr::kable()

Plot some relationships between institutions and outcomes

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = log(unna_gdppc), color = factor(groups_fraser)) +
    geom_point() +
    labs(x = "Property rights", y = "Log of GDP per capita")

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = vdem_pubcorr, color = factor(groups_fraser)) +
    geom_point() +
    labs(x = "Property rights", y = "Corruption")

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = wdi_lifexp, color = factor(groups_fraser)) +
    geom_point() +
    labs(x = "Property rights", y = "Life expectancy")

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = gea_ea3544f, color = factor(groups_fraser)) +
    geom_point() +
    labs(x = "Property rights", y = "Education, female (35-44)")

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = wdi_gdpcapgr, color = factor(groups_fraser)) +
    geom_point() +
    labs(x = "Property rights", y = "Growth rate")

Calculate summary statistics based on the cluster analysis

# calculate means of the outcome variables for each cluster
outcomes_means <- institutions %>%
  group_by(groups_fraser) %>%
  summarise_at(vars(outcomes), funs(mean)) %>%  
  gather(key = "variable", value = "mean", -groups_fraser)

# calculate the standard error 
std.err <- function(x) { sqrt(var(x, na.rm=TRUE)/length(na.omit(x))) }

# calculate standard errors (uses function `std.err` defined above)
outcomes_se <- institutions %>%
  group_by(groups_fraser) %>%
  summarise_at(vars(outcomes), funs(std.err)) %>% 
  gather(key = "variable", value = "se", -groups_fraser)

outcomes_summary <- full_join(outcomes_means, outcomes_se)

# show results as table
knitr::kable(outcomes_summary, digits=2)

Show summary statistics results as bar plots

outcomes_summary %>%
  ggplot() +
    aes(x = groups_fraser, y = mean) +
    geom_bar(stat = "identity") +
    geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.3) +
    facet_wrap(~variable, scale = "free", ncol = 2)

Compare different clustering methods

## calculate hierarchical cluster
## comparison of different clustering methods
hc1 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="single")
hc2 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="complete")
hc3 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="average")
hc4 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="mcquitty")
hc5 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="median")
hc6 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="centroid")
hc7 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="ward.D")
hc8 <- fraser_scaled %>% dist(method = "euclidean") %>% hclust(method="ward.D2")

fraser_dendlist <- dendlist(
  as.dendrogram(hc1),
  as.dendrogram(hc2),
  as.dendrogram(hc3),
  as.dendrogram(hc4),
  as.dendrogram(hc5),
  as.dendrogram(hc6),
  as.dendrogram(hc7),
  as.dendrogram(hc8)
)

hclust_methods <- c("single", "complete", "average", "mcquitty",
                    "median", "centroid", "ward.D", "ward.D2")
names(fraser_dendlist) <- hclust_methods

# calculate correlations
fraser_cor <- cor.dendlist(fraser_dendlist)

# plot correlation table
corrplot.mixed(fraser_cor)

Compare countries based on outcomes

hc_o1 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="single")
hc_o2 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="complete")
hc_o3 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="average")
hc_o4 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="mcquitty")
hc_o5 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="median")
hc_o6 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="centroid")
hc_o7 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="ward.D")
hc_o8 <- outcomes_scaled %>% dist(method = "euclidean") %>% hclust(method="ward.D2")

outcomes_dendlist <- dendlist(
  as.dendrogram(hc_o1),
  as.dendrogram(hc_o2),
  as.dendrogram(hc_o3),
  as.dendrogram(hc_o4),
  as.dendrogram(hc_o5),
  as.dendrogram(hc_o6),
  as.dendrogram(hc_o7),
  as.dendrogram(hc_o8)
)

hclust_methods <- c("single", "complete", "average", "mcquitty",
                    "median", "centroid", "ward.D", "ward.D2")
names(outcomes_dendlist) <- hclust_methods

# calculate correlations
outcomes_cor <- cor.dendlist(outcomes_dendlist)

# plot correlation table
corrplot.mixed(outcomes_cor)

# plot cluster
plot(hc_o7, labels=outcomes_scaled$cname, hang=-1, cex = 0.7)

Correlation between institutions tree and outcomes tree

Which clustering method leads to the highest correlation between the tree based on institutions and the tree based on outcomes?

compar_list <- dendlist(
  as.dendrogram(hc1),
  as.dendrogram(hc2),
  as.dendrogram(hc3),
  as.dendrogram(hc4),
  as.dendrogram(hc5),
  as.dendrogram(hc6),
  as.dendrogram(hc7),
  as.dendrogram(hc8),
  as.dendrogram(hc_o1),
  as.dendrogram(hc_o2),
  as.dendrogram(hc_o3),
  as.dendrogram(hc_o4),
  as.dendrogram(hc_o5),
  as.dendrogram(hc_o6),
  as.dendrogram(hc_o7),
  as.dendrogram(hc_o8)
)

names(compar_list) <- c("fraser1", "fraser2","fraser3","fraser4","fraser5",
                        "fraser6","fraser7","fraser8",
                        "outcomes1", "outcomes2", "outcomes3", "outcomes4",
                        "outcomes5", "outcomes6", "outcomes7", "outcomes8")

# correlations plot
corrplot.mixed(cor.dendlist(compar_list))

Another type of analysis that we could make: Instead of looking at different methods using the same set of variables, we can ask which set of institutional variables leads to the highest correlation to the outcomes tree?

Visual comparison of trees, between institutions and outcomes:

dend_fraser   <- as.dendrogram(hc7, h = 25)
dend_outcomes <- as.dendrogram(hc_o7, h = 25)

tanglegram(dend_fraser, dend_outcomes,
           sort = TRUE,
           common_subtrees_color_lines = FALSE,
           highlight_distinct_edges  = FALSE,
           highlight_branches_lwd = FALSE,
           margin_inner = 0.5,
           main_left = "Tree based on institutions",
           main_right = "Tree based on outcomes")

Visualize the results

Some plots:

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = gea_ea3544f, color = factor(kmeans)) +
    geom_point() +
    labs(x = "Property rights", y = "Education, female (35-44)")

institutions %>% 
  ggplot() +
    aes(x = fi_legprop_pd, y = log(unna_gdppc), color = factor(kmeans)) +
    geom_point() +
    labs(x = "Property rights", y = "Log of GDP per capita")

As you can see, the result of the k-means clustering is very similar to that of hierarchical clustering.

Alternative way to visualize using the first two principal components:

institutions %>% 
  select(fraser) %>% 
  fviz_cluster(km, geom = "point", data = .) + ggtitle("k = 3")

Find the natural number of clusters

How many clusters should we make? See which number maximizes the silhouette:

institutions %>% 
  select(fraser) %>% 
  fviz_nbclust(kmeans, method = "silhouette")

The sihouette is a measure that takes into account both how close the different cases within a cluster are to each other and how distant the different clusters are from one another.

The kmeans with 2 clusters:

institutions %>% 
  select(fraser) %>% 
  kmeans(2) %>% 
  fviz_cluster(geom = "point", 
               data = institutions %>% select(fraser)) + ggtitle("k = 2")

institutions %>% 
  select(fraser) %>% 
  kmeans(8) %>% 
  fviz_cluster(geom = "point", 
               data = institutions %>% select(fraser)) + ggtitle("k = 8")