Bivariate distribution dynamics analysis in R

Suggested citation:

Mendez C. (2020). Bivariate distribution dynamics analysis in R. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/tutorial-bivariate-distribution-dynamics

This work is licensed under the Creative Commons Attribution-Non Commercial-Share Alike 4.0 International License.

Acknowledgment:

Material adapted from multiple sources, in particular the dataset is from Magrini (2007).

1 Replication files

• The tutorial is self-contained. No aditional file is needed.

• If you are a member of the QuaRCS lab, you can run this tutorial in R Studio Cloud

2 Libraries

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)
library(skimr)
library(kableExtra)    # html tables 
library(pdfCluster)    # density based clusters
library(hdrcde)        # conditional density estimation 
library(plotly)

library(intoo)
library(barsurf)
library(bivariate)

library(np)

library(basetheme)
basetheme("minimal")

library(viridis)
library(ggpointdensity)
library(isoband)

#library(MASS)
library(KernSmooth)


# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=2, scipen=999)

3 Tutorial objectives

• Study the dynamics of univariate densities

• Compute the bandwidth of a density

• Study mobility plots

• Study bi-variate densities

• Study density-based clustering methods

• Study conditional bi-variate densities

4 Import data

We will use two hypothetical cross-sectional series.

• The first (x) series was produced by drawing a random sample of 1000 observations from a univariate normal distribution.
• The second (y) series was produced by merging and sorting two random samples of 500 observations.

The mean and the standard deviation of these two series respectively matched those of the logarithm of per capita Gross Value Added observed for the Italian Provinces in 1996 and in 2002. For this reason assume that the analysis has been performed over a 6-year time period.

dat <- read_csv("https://raw.githubusercontent.com/ds777/sample-datasets/master/simData.csv")

5 Transform data

Since the data is in log terms, let us rename the variables and add new variables.

dat <- dat %>% 
  rename(log_x = x, log_y = y) %>% 
  mutate(x = exp(log_x), y = exp(log_y)) %>% 
  select(x, y, everything()) 

dat <- dat %>% 
  mutate(
    rel_x = x/mean(x),
    rel_y = y/mean(y),
    rel_log_x = log_x/mean(log_x),
    rel_log_y = log_x/mean(log_y),
  )
dat

6 Descriptive statistics

skim(dat) %>% 
  kable() %>% 
  kable_styling()

skim_type	skim_variable	n_missing	complete_rate	numeric.mean	numeric.sd	numeric.p0	numeric.p25	numeric.p50	numeric.p75	numeric.p100	numeric.hist
numeric	x	0	1	14883.41	3746.62	3344.50	12362.51	14974.61	17321.7	26506.5	▁▃▇▃▁
numeric	y	0	1	16441.00	4011.14	7400.07	12992.51	16244.71	20021.6	25394.5	▂▇▅▇▂
numeric	log_x	0	1	9.57	0.28	8.12	9.42	9.61	9.8	10.2	▁▁▂▇▃
numeric	log_y	0	1	9.68	0.25	8.91	9.47	9.70	9.9	10.1	▁▃▇▇▇
numeric	rel_x	0	1	1.00	0.25	0.22	0.83	1.01	1.2	1.8	▁▃▇▃▁
numeric	rel_y	0	1	1.00	0.24	0.45	0.79	0.99	1.2	1.5	▂▇▅▇▂
numeric	rel_log_x	0	1	1.00	0.03	0.85	0.98	1.00	1.0	1.1	▁▁▂▇▃
numeric	rel_log_y	0	1	0.99	0.03	0.84	0.97	0.99	1.0	1.1	▁▁▂▇▃

7 Univariate dynamics

7.1 Select bandwiths

select bandwidth based on function dpik from the package KernSmooth

h_rel_x <- dpik(dat$rel_x)
h_rel_x

[1] 0.065

h_rel_y <- dpik(dat$rel_y)
h_rel_y

[1] 0.039

7.2 Plot each density

dis_rel_x <- bkde(dat$rel_x, bandwidth = h_rel_x)
dis_rel_x <- as.data.frame(dis_rel_x)
ggplot(dis_rel_x, aes(x, y)) + geom_line() + 
  theme_minimal() 

dis_rel_y <- bkde(dat$rel_y, bandwidth = h_rel_y)
dis_rel_y <- as.data.frame(dis_rel_y)
ggplot(dis_rel_y, aes(x, y)) + geom_line() + 
  theme_minimal() 

7.3 Plot both densities

7.3.1 Method 1

Keep the orignal bandwiths of the package KernSmooth

domain_initial <- dis_rel_x$x
domain_final   <- dis_rel_y$x

density_initial <- dis_rel_x$y
density_final   <- dis_rel_y$y

densities <- cbind(domain_initial, density_initial, domain_final, density_final)

densities <-  as.data.frame(densities)

densities_plot <- densities %>% 
  ggplot()+
  theme_minimal()+
  geom_line(aes(domain_initial, density_initial))+
  geom_line(aes(domain_final,density_final), linetype = "dashed")+
 labs(subtitle = "",
       x = "Relative Variable",
       y = "Density") + 
  geom_label(
    label="Year 2002", 
    x= 1.525,
    y= 0.9,
    label.size = 0.35,
    color = "black",
  ) +
   geom_label(
    label="Year 1996", 
    x= 1.75,
    y= 0.2,
    label.size = 0.35,
    color = "black",
  )
densities_plot

Note that you have adjust the labels manually

• Interactive plotly version

densities_plot2 <- densities %>% 
  ggplot()+
  theme_minimal()+
  geom_line(aes(domain_initial, density_initial))+
  geom_line(aes(domain_final,density_final), linetype = "dashed")+
 labs(subtitle = "",
       x = "Relative Variable",
       y = "Density") 

ggplotly(densities_plot2)

Manual labels are not yet implemented in the ggplotly function

7.3.2 Method 2

using the bandwidth default of ggplot

rel_x <- dat %>% 
  select(rel_x) %>% 
  rename(rel_var = rel_x) %>% 
  mutate(year = 1996)

rel_y <- dat %>% 
  select(rel_y) %>% 
  rename(rel_var = rel_y) %>% 
  mutate(year = 2002)

rel_xy <- bind_rows(rel_x, rel_y)

rel_xy <- rel_xy %>% 
  mutate(year = as.factor(year))
head(rel_xy)

dis_rel_xy <- ggplot(rel_xy, aes(x=rel_var, color=year)) +
  geom_density() + 
  theme_minimal() 
dis_rel_xy

Using plotly

ggplotly(dis_rel_xy)

8 Bivariate density

8.1 Mobility scatterplot

dat %>% 
  ggplot(aes(x = rel_x, y = rel_y)) +
  geom_point(alpha=0.5) +
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative y",
       x = "Relative x",
       y = "") +
  theme(text=element_text(family="Palatino")) 

Fit a non-linear function

dat %>% 
  ggplot(aes(x = rel_x, y = rel_y)) +
  geom_point(alpha=0.5) + 
  geom_smooth() + 
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative y",
       x = "Relative x",
       y = "") +
  theme(text=element_text(family="Palatino")) 

Not that the nonlinear fit crosses the 45-degree line two times from above.

8.2 Using geom_pointdensity

dat %>% 
ggplot(aes(x = rel_x, y = rel_y)) +
  geom_pointdensity() +
  scale_color_viridis() +
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative y") + 
  xlab("Relative x") +
  ylab("")

8.3 Using the KernSmooth package

xy <- cbind(dat$rel_x, dat$rel_y)
xy_dis <- bkde2D(xy, bandwidth = c(h_rel_x, h_rel_y)) 

contour(xy_dis$x1, xy_dis$x2, xy_dis$fhat)
abline(a=0, b=1)

abline(h=1, v=1)

Interactive

plot_ly(x=xy_dis$x1, y=xy_dis$x2, z=xy_dis$fhat, type = "contour", contours = list(showlabels = TRUE))  %>%
  colorbar(title = "Density")

persp(xy_dis$fhat)

Interactive version

plot_ly(x=xy_dis$x1, y=xy_dis$x2, z=xy_dis$fhat) %>% add_surface()

8.4 Using the Bivariate package

bivariate <- kbvpdf(dat$rel_x, dat$rel_y, h_rel_x, h_rel_y) 

plot(bivariate,
      xlab="Relative x", 
      ylab="Relative y")
abline(a=0, b=1)

abline(h=1, v=1)

plot(bivariate,
      TRUE,
      xlab="Relative x", 
      ylab="Relative y")

8.5 Using ggplot (stat_density_2d())

dat %>% 
  ggplot(aes(x = rel_x, y = rel_y)) +
  geom_point(color = "lightgray") + 
  geom_smooth() + 
  #geom_smooth(method=lm, se=FALSE) + 
  stat_density_2d() +
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept = 1, linetype="dashed") + 
  geom_vline(xintercept = 1, linetype="dashed") +
  theme_minimal() +
  labs(subtitle = "Relative y",
       x = "Relative x",
       y = "") +
  theme(text=element_text(family="Palatino")) 

dat %>% 
  ggplot(aes(x = rel_x, y = rel_y)) +
        stat_density_2d(aes(fill = stat(nlevel)), geom = "polygon") + 
  scale_fill_viridis_c() +
        geom_abline(aes(intercept = 0, slope = 1)) +
        geom_hline(yintercept = 1, linetype="dashed") + 
        geom_vline(xintercept = 1, linetype="dashed") + 
  theme_minimal() +
        labs(x = "Relative X",
             y = "Relative Y") +
        theme(text=element_text(size=8, family="Palatino"))

9 Density-based clusters

xy_cluster <- pdfCluster(dat[, 5:6]) 
summary(xy_cluster)

An S4 object of class "pdfCluster"

Call: pdfCluster(x = dat[, 5:6])

Initial groupings: 
 label    1   2  NA 
 count  296 297 407 

Final groupings: 
 label    1   2 
 count  504 496 

Groups tree (here 'h' denotes 'height'):
--[dendrogram w/ 1 branches and 2 members at h = 1]
  `--[dendrogram w/ 2 branches and 2 members at h = 0.593]
     |--leaf "1 " 
     `--leaf "2 " (h= 0.152  )

9.1 Core clusters

plot(xy_cluster, 
      stage = 0,
      which = 3,
      frame = F)
abline(a=0, b=1)

abline(h=1, v=1)

9.2 Full clustering

plot(xy_cluster, 
      which = 3,
      frame = F)
abline(a=0, b=1)

abline(h=1, v=1)

9.3 Cluster tree

plot(xy_cluster, which = 2)

9.4 Mode function

plot(xy_cluster, which = 1)

10 Conditional density analysis

10.1 Using the hdrcde package

xy_cde <- cde(dat$rel_x, dat$rel_y)

Increase the number of intervals to 60

xy_cde2 <- cde(dat$rel_x, dat$rel_y, nxmargin = 60)

plot(xy_cde)

plot(xy_cde2)

High density regions

plot(xy_cde, plot.fn="hdr")
abline(a=0, b=1)

abline(h=1, v=1)

plot(xy_cde2, plot.fn="hdr")
abline(a=0, b=1)

abline(h=1, v=1)

10.2 Using the np package

Compute adaptive bandwith based on cross-validation

bw_c_ad_cv <- npcdensbw(
  formula = dat$rel_y ~ dat$rel_x, 
  bwtype = "adaptive_nn") 

summary(bw_c_ad_cv)

Conditional density data (1000 observations, 2 variable(s))
(1 dependent variable(s), and 1 explanatory variable(s))

Bandwidth Selection Method: Maximum Likelihood Cross-Validation
Formula: dat$rel_y ~ dat$rel_x
Bandwidth Type: Adaptive Nearest Neighbour
Objective Function Value: 5502 (achieved on multistart 1)

Exp. Var. Name: dat$rel_x Bandwidth: 2  

Dep. Var. Name: dat$rel_y Bandwidth: 2 

Continuous Kernel Type: Second-Order Gaussian
No. Continuous Explanatory Vars.: 1
No. Continuous Dependent Vars.: 1
Estimation Time: 68 seconds

Compute conditional density object

cdist_bwCV <- npcdens(bws = bw_c_ad_cv)
summary(cdist_bwCV)

Conditional Density Data: 1000 training points, in 2 variable(s)
(1 dependent variable(s), and 1 explanatory variable(s))

                        dat$rel_y
Dep. Var. Bandwidth(s):         2
                        dat$rel_x
Exp. Var. Bandwidth(s):         2

Bandwidth Type: Adaptive Nearest Neighbour
Log Likelihood: 5948

Continuous Kernel Type: Second-Order Gaussian
No. Continuous Explanatory Vars.: 1
No. Continuous Dependent Vars.: 1

## Better copy and paste it in the console because it is an animation
# plot(bw_c_ad_cv)

plot(cdist_bwCV, 
      xtrim = -0.2, 
      view = "fixed",
      main = "",
      theta=350,
      phi=20)

11 References

• Magrini, S. (2007). Analysing convergence through the distribution dynamics approach: why and how?. University Ca’Foscari of Venice, Dept. of Economics Research Paper Series No, 13.

• Mendez C. (2020). Classical sigma and beta convergence analysis in R: Using the REAT 2.1 Package. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/classical-convergence-reat21

• Mendez C. (2020). Univariate distribution dynamics in R: Using the ggridges package. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/univariate-distribution-dynamics

• Mendez, C. (2020) Regional efficiency convergence and efficiency clusters. Asia-Pacific Journal of Regional Science, 1-21.

• Mendez, C. (2019). Lack of Global Convergence and the Formation of Multiple Welfare Clubs across Countries: An Unsupervised Machine Learning Approach. Economies, 7(3), 74.

• Mendez, C. (2019). Overall efficiency, pure technical efficiency, and scale efficiency across provinces in Indonesia 1990 and 2010. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/efficiency-clusters-indonesia-1990-2010

• Mendez-Guerra, C. (2018). On the distribution dynamics of human development: Evidence from the metropolitan regions of Bolivia’’. Economics Bulletin, 38(4), 2467-2475.

END