01-data-exploration

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(tidyverse)
library(summarytools)
library(xtable)
library(knitr)
library(ExPanDaR)
library(plotly)

1 Use the package “tidyverse” to import the data:

dat <- read_csv("00-raw-data/data.csv")
dat

2 Use the package “summarytools” to summarize the data:

descr(dat[,5:16])

Descriptive Statistics
Data Frame: dat
N: 46

Table continues below

	ceei92	deei92	deei92v2	ceei08	deei08	deei08v2	ceei92re	deei92v2re
Mean	0.51	0.13	0.87	0.58	0.16	0.84	0.38	0.87
Std.Dev	0.20	0.10	0.10	0.21	0.12	0.12	0.15	0.10
Min	0.17	0.00	0.63	0.19	0.00	0.53	0.13	0.63
Q1	0.38	0.06	0.82	0.40	0.08	0.76	0.28	0.82
Median	0.51	0.12	0.88	0.58	0.13	0.87	0.38	0.88
Q3	0.58	0.18	0.94	0.70	0.24	0.92	0.43	0.94
Max	1.34	0.37	1.00	1.33	0.47	1.00	1.00	1.00
MAD	0.14	0.09	0.09	0.17	0.14	0.14	0.11	0.09
IQR	0.20	0.11	0.11	0.28	0.15	0.15	0.15	0.11
CV	0.39	0.79	0.12	0.36	0.76	0.14	0.39	0.12
Skewness	1.43	0.43	-0.43	0.72	0.39	-0.39	1.43	-0.43
SE.Skewness	0.35	0.35	0.35	0.35	0.35	0.35	0.35	0.35
Kurtosis	4.43	-0.73	-0.73	1.96	-0.59	-0.59	4.43	-0.73
N.Valid	46.00	46.00	46.00	46.00	46.00	46.00	46.00	46.00
Pct.Valid	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

	ceei08re	deei08v2re
Mean	0.43	0.84
Std.Dev	0.16	0.12
Min	0.14	0.53
Q1	0.30	0.76
Median	0.44	0.87
Q3	0.52	0.92
Max	1.00	1.00
MAD	0.13	0.14
IQR	0.21	0.15
CV	0.36	0.14
Skewness	0.72	-0.39
SE.Skewness	0.35	0.35
Kurtosis	1.96	-0.59
N.Valid	46.00	46.00
Pct.Valid	100.00	100.00

print(dfSummary(dat[,5:16], graph.magnif = 0.75), method = 'render')

Data Frame Summary

dat

N: 46

No	Variable	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	ceei92 [numeric]	mean (sd) : 0.51 (0.2) min < med < max : 0.17 < 0.51 < 1.34 IQR (CV) : 0.2 (0.39)	41 distinct values		46 (100%)	0 (0%)
2	deei92 [numeric]	mean (sd) : 0.13 (0.1) min < med < max : 0 < 0.12 < 0.37 IQR (CV) : 0.11 (0.79)	36 distinct values		46 (100%)	0 (0%)
3	deei92v2 [numeric]	mean (sd) : 0.87 (0.1) min < med < max : 0.63 < 0.88 < 1 IQR (CV) : 0.11 (0.12)	36 distinct values		46 (100%)	0 (0%)
4	ceei08 [numeric]	mean (sd) : 0.58 (0.21) min < med < max : 0.19 < 0.58 < 1.33 IQR (CV) : 0.28 (0.36)	44 distinct values		46 (100%)	0 (0%)
5	deei08 [numeric]	mean (sd) : 0.16 (0.12) min < med < max : 0 < 0.13 < 0.47 IQR (CV) : 0.15 (0.76)	37 distinct values		46 (100%)	0 (0%)
6	deei08v2 [numeric]	mean (sd) : 0.84 (0.12) min < med < max : 0.53 < 0.87 < 1 IQR (CV) : 0.15 (0.14)	37 distinct values		46 (100%)	0 (0%)
7	ceei92re [numeric]	mean (sd) : 0.38 (0.15) min < med < max : 0.13 < 0.38 < 1 IQR (CV) : 0.15 (0.39)	41 distinct values		46 (100%)	0 (0%)
8	deei92re [logical]		All NA's		0 (0%)	46 (100%)
9	deei92v2re [numeric]	mean (sd) : 0.87 (0.1) min < med < max : 0.63 < 0.88 < 1 IQR (CV) : 0.11 (0.12)	36 distinct values		46 (100%)	0 (0%)
10	ceei08re [numeric]	mean (sd) : 0.43 (0.16) min < med < max : 0.14 < 0.44 < 1 IQR (CV) : 0.21 (0.36)	44 distinct values		46 (100%)	0 (0%)
11	deei08re [logical]		All NA's		0 (0%)	46 (100%)
12	deei08v2re [numeric]	mean (sd) : 0.84 (0.12) min < med < max : 0.53 < 0.87 < 1 IQR (CV) : 0.15 (0.14)	37 distinct values		46 (100%)	0 (0%)

Generated by summarytools 0.8.8 (R version 3.5.1)
2019-02-07

# First save the results
dat_by_area <- by(data = dat[,5:16],
                  INDICES = dat$area, 
                  FUN = descr, stats = c("mean", "med", "sd", "cv" ,"iqr", "min", "max"), 
                  transpose = TRUE)
# Then use view(), like so:
view(dat_by_area, method = "pander")

Descriptive Statistics
Data Frame: dat
Group: area = Chubu
N: 9

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.54	0.52	0.10	0.19	0.09	0.38	0.74
deei92	0.15	0.11	0.12	0.79	0.18	0.00	0.37
deei92v2	0.85	0.89	0.12	0.14	0.18	0.63	1.00
ceei08	0.66	0.61	0.11	0.16	0.18	0.54	0.81
deei08	0.14	0.12	0.11	0.80	0.09	0.00	0.32
deei08v2	0.86	0.88	0.11	0.13	0.08	0.68	1.00
ceei92re	0.40	0.38	0.08	0.19	0.06	0.28	0.55
deei92v2re	0.85	0.89	0.12	0.14	0.18	0.63	1.00
ceei08re	0.49	0.45	0.08	0.16	0.13	0.40	0.61
deei08v2re	0.86	0.88	0.11	0.13	0.08	0.68	1.00

Group: area = Chugoku
N: 5

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.32	0.25	0.15	0.46	0.25	0.20	0.51
deei92	0.08	0.09	0.06	0.66	0.06	0.00	0.14
deei92v2	0.92	0.91	0.06	0.06	0.06	0.86	1.00
ceei08	0.33	0.28	0.15	0.46	0.28	0.19	0.51
deei08	0.13	0.13	0.09	0.70	0.01	0.00	0.25
deei08v2	0.87	0.87	0.09	0.10	0.00	0.75	1.00
ceei92re	0.24	0.19	0.11	0.46	0.19	0.15	0.38
deei92v2re	0.92	0.91	0.06	0.06	0.06	0.86	1.00
ceei08re	0.25	0.21	0.11	0.46	0.21	0.14	0.38
deei08v2re	0.87	0.87	0.09	0.10	0.00	0.75	1.00

Group: area = Hokkaido
N: 1

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.46	0.46	NA	NA	0.00	0.46	0.46
deei92	0.30	0.30	NA	NA	0.00	0.30	0.30
deei92v2	0.70	0.70	NA	NA	0.00	0.70	0.70
ceei08	0.38	0.38	NA	NA	0.00	0.38	0.38
deei08	0.47	0.47	NA	NA	0.00	0.47	0.47
deei08v2	0.53	0.53	NA	NA	0.00	0.53	0.53
ceei92re	0.34	0.34	NA	NA	0.00	0.34	0.34
deei92v2re	0.70	0.70	NA	NA	0.00	0.70	0.70
ceei08re	0.29	0.29	NA	NA	0.00	0.29	0.29
deei08v2re	0.53	0.53	NA	NA	0.00	0.53	0.53

Group: area = Kanto
N: 7

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.64	0.61	0.34	0.53	0.18	0.31	1.34
deei92	0.15	0.18	0.09	0.57	0.08	0.00	0.27
deei92v2	0.85	0.82	0.09	0.10	0.08	0.73	1.00
ceei08	0.66	0.63	0.33	0.49	0.20	0.33	1.33
deei08	0.20	0.22	0.12	0.60	0.14	0.00	0.33
deei08v2	0.80	0.78	0.12	0.15	0.14	0.67	1.00
ceei92re	0.47	0.45	0.25	0.53	0.14	0.23	1.00
deei92v2re	0.85	0.82	0.09	0.10	0.08	0.73	1.00
ceei08re	0.50	0.48	0.24	0.49	0.15	0.24	1.00
deei08v2re	0.80	0.78	0.12	0.15	0.14	0.67	1.00

Group: area = Kinki
N: 7

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.60	0.61	0.24	0.40	0.43	0.31	0.86
deei92	0.06	0.00	0.07	1.28	0.12	0.00	0.15
deei92v2	0.94	1.00	0.07	0.08	0.12	0.85	1.00
ceei08	0.65	0.70	0.24	0.37	0.37	0.32	0.95
deei08	0.11	0.08	0.11	1.06	0.21	0.00	0.25
deei08v2	0.89	0.92	0.11	0.13	0.21	0.75	1.00
ceei92re	0.44	0.45	0.18	0.40	0.32	0.23	0.64
deei92v2re	0.94	1.00	0.07	0.08	0.12	0.85	1.00
ceei08re	0.49	0.52	0.18	0.37	0.28	0.24	0.71
deei08v2re	0.89	0.92	0.11	0.13	0.21	0.75	1.00

Group: area = Kyushu
N: 7

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.47	0.57	0.15	0.33	0.16	0.17	0.58
deei92	0.12	0.14	0.04	0.34	0.07	0.08	0.17
deei92v2	0.88	0.86	0.04	0.05	0.07	0.83	0.92
ceei08	0.57	0.63	0.18	0.31	0.19	0.23	0.72
deei08	0.18	0.17	0.11	0.62	0.09	0.00	0.36
deei08v2	0.82	0.83	0.11	0.13	0.09	0.64	1.00
ceei92re	0.35	0.42	0.11	0.33	0.12	0.13	0.43
deei92v2re	0.88	0.86	0.04	0.05	0.07	0.83	0.92
ceei08re	0.43	0.48	0.13	0.31	0.14	0.17	0.54
deei08v2re	0.82	0.83	0.11	0.13	0.09	0.64	1.00

Group: area = Shikoku
N: 4

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.40	0.39	0.09	0.23	0.09	0.30	0.52
deei92	0.04	0.00	0.07	1.88	0.04	0.00	0.14
deei92v2	0.96	1.00	0.07	0.07	0.04	0.86	1.00
ceei08	0.48	0.50	0.13	0.26	0.20	0.35	0.60
deei08	0.11	0.11	0.10	0.88	0.14	0.01	0.23
deei08v2	0.89	0.89	0.10	0.11	0.14	0.77	0.99
ceei92re	0.30	0.29	0.07	0.23	0.07	0.22	0.39
deei92v2re	0.96	1.00	0.07	0.07	0.04	0.86	1.00
ceei08re	0.36	0.37	0.10	0.26	0.15	0.26	0.45
deei08v2re	0.89	0.89	0.10	0.11	0.14	0.77	0.99

Group: area = Tohoku
N: 6

	Mean	Median	Std.Dev	CV	IQR	Min	Max
ceei92	0.51	0.55	0.08	0.16	0.10	0.37	0.57
deei92	0.22	0.27	0.12	0.55	0.11	0.00	0.30
deei92v2	0.78	0.73	0.12	0.15	0.11	0.70	1.00
ceei08	0.58	0.58	0.11	0.20	0.14	0.40	0.71
deei08	0.21	0.21	0.15	0.74	0.20	0.00	0.40
deei08v2	0.79	0.79	0.15	0.19	0.20	0.61	1.00
ceei92re	0.38	0.41	0.06	0.16	0.07	0.28	0.42
deei92v2re	0.78	0.73	0.12	0.15	0.11	0.70	1.00
ceei08re	0.43	0.43	0.09	0.20	0.11	0.30	0.53
deei08v2re	0.79	0.79	0.15	0.19	0.20	0.61	1.00

3 Use the package “ExPanD” for panel data explorations:

Fist import the dataframe in panel form. Note that I am using the read.csv instead of the read_csv function. This is because the former automatically imports character variables as factor varaibles

pan <- read.csv("00-raw-data/panel.csv")
pan

Let’s use the ExPanD() function. Here is a nice tutorial

Because we will be running analyses over time, we will add a factor variable that indicates the year

pan <- pan %>% 
  mutate(year2 = factor(year, levels = c(1992, 2008)))

Also for a better presentation of tables, let us rename some variables

3.1 Table Displaying Extreme Observations

tab1 <- prepare_ext_obs_table(pan, n = 5,
                           cs_id = c("area", "region"),
                           ts_id = "year", var = "ceeire")
tab1$df

tab1 <- prepare_ext_obs_table(pan, n = 5,
                           cs_id = c("area", "region"),
                           ts_id = "year", var = "deeiv2re")
tab1$df

3.2 Group Bar Graph

graph1 <- prepare_by_group_bar_graph(pan, "area", "ceeire", mean, order_by_stat = TRUE, color = "blue")
p1 <- graph1$plot +
      xlab("") +
      ylab("Conventional Environmental Efficiency Indicator (GDP/CO2)")
ggplotly(p1)

graph2 <- prepare_by_group_bar_graph(pan, "area", "deeiv2re", mean, order_by_stat = TRUE, color = "blue")
p2 <- graph2$plot +
  xlab("") +
  ylab("DEA-Based Environmental Efficiency Indicator")
ggplotly(p2)

graph3 <- prepare_by_group_bar_graph(pan, "year2", "ceeire", mean, order_by_stat = TRUE, color = "blue")
p3 <- graph3$plot +
  xlab("") +
  ylab("Conventional Environmental Efficiency Indicator")
ggplotly(p3)

graph4 <- prepare_by_group_bar_graph(pan, "year2", "deeiv2re", mean, order_by_stat = TRUE, color = "blue")
p4 <- graph4$plot +
  xlab("") +
  ylab("DEA-Based Environmental Efficiency Indicator")
ggplotly(p4)

3.3 Group Violin Graph

graph5 <- prepare_by_group_violin_graph(pan, "area", "ceeire")
p5 <- graph5  +
  xlab("") +
  ylab("Conventional Environmental Efficiency Indicator")
ggplotly(p5)

graph6 <- prepare_by_group_violin_graph(pan, "area", "deeiv2re")
p6 <- graph6  +
  xlab("") +
  ylab("DEA-Based Environmental Efficiency Indicator")
ggplotly(p6)

graph7 <- prepare_by_group_violin_graph(pan, "year2", "ceeire")
p7 <- graph7  +
  xlab("") +
  ylab("Conventional Environmental Efficiency Indicator")
ggplotly(p7)

graph8 <- prepare_by_group_violin_graph(pan, "year2", "deeiv2re")
p8 <- graph8  +
  xlab("") +
  ylab("DEA-Based Environmental Efficiency Indicator")
ggplotly(p8)

3.4 Trend Graph

df90 <- pan %>% 
  select(year, CEEI = ceeire, DEEI = deeiv2re )
graph9 <- prepare_trend_graph(df90[c("year", "CEEI", "DEEI")], "year")
graph9$plot  + 
  xlab("") +
  ylab("Relative Environmental Efficiency")

df90 <- pan %>% 
  select(year, CEEI = ceeire, DEEI = deeiv2re )
graph9 <- prepare_trend_graph(df90[c("year", "CEEI", "DEEI")], "year")
p9 <- graph9$plot  + 
  xlab("") +
  ylab("Relative Environmental Efficiency")
ggplotly(p9)

3.5 Quantile Trend Graph

graph10 <- prepare_quantile_trend_graph(pan[c("year", "ceeire")], "year", c(0.05, 0.25, 0.5, 0.75, 0.95))
p10 <- graph10$plot + 
  ylab("Conventional Environmental Efficiency Indicator") +
  xlab("")
ggplotly(p10)

graph11 <- prepare_quantile_trend_graph(pan[c("year", "deeiv2re")], "year", c(0.05, 0.25, 0.5, 0.75, 0.95))
p11 <- graph11$plot  + 
  ylab("DEA-Based Environmental Efficiency Indicator") +
  xlab("")
ggplotly(p11)

3.6 Correlation Table

df1 <- pan %>% 
  select(CEEI = ceeire, DEEI = deeiv2re)
tab1 <- prepare_correlation_table(df1, digits = 2, bold = 0.05, format = "html")
tab1

$df_corr

$df_prob

$df_n

$kable_ret

	A	B
A: CEEI		0.32
B: DEEI	0.30
This table reports Pearson correlations above and Spearman correlations below the diagonal. Number of observations: 92. Correlations with significance levels below 5% appear in bold print.

NA

graph12 <- prepare_correlation_graph(df1)

graph12

$df_corr
          CEEI      DEEI
CEEI 1.0000000 0.3189966
DEEI 0.2991585 1.0000000

$df_prob
            CEEI        DEEI
CEEI          NA 0.001939981
DEEI 0.003769318          NA

$df_n
     CEEI DEEI
CEEI   92   92
DEEI   92   92

3.7 Scatterplots

graph13 <- prepare_scatter_plot(pan, x="ceei", y="deeiv2", color="area", loess = 1)
p13 <- graph13  + 
  ylab("DEA-Based Environmental Efficiency Indicator") +
  xlab("Conventional Environmental Efficiency Indicator")
ggplotly(p13)

graph14 <- prepare_scatter_plot(pan, x="ceeire", y="deeiv2re", color="area", loess = 1)
p14 <- graph14   + 
  ylab("DEA-Based Environmental Efficiency Indicator") +
  xlab("Conventional Environmental Efficiency Indicator")
ggplotly(p14)

3.8 Regression Tables

First, lets define the dataframe

df100 <- pan %>% 
  select(regionid, year, year2, region, CEEI = ceeire, DEEI = deeiv2re )
df100

Define list of variables and effects

dvs = c("DEEI")
idvs = c("CEEI")
feffects = c("year2")
clusters = c("region")

Run the regression

t1 <- prepare_regression_table(df100, dvs, idvs, feffects, clusters, format = "html")
htmltools::HTML(t1$table)

	Dependent variable:
	DEEI
CEEI	0.256***
	(0.061)
Estimator	ols
Fixed effects	year2
Std. errors clustered	region
Observations	92
R2	0.145
Adjusted R2	0.126
Note:	*p<0.1; **p<0.05; ***p<0.01

For regressions with more variables see the reference guide

t2 <- prepare_regression_table(df100, dvs, idvs, byvar="year2")
htmltools::HTML(t2$table)

	Dependent variable:
	DEEI
	Full Sample	1992	2008
	(1)	(2)	(3)
CEEI	0.231***	0.187*	0.319***
	(0.072)	(0.097)	(0.106)
Constant	0.762***	0.802***	0.701***
	(0.031)	(0.040)	(0.049)
Estimator	ols	ols	ols
Fixed effects	None	None	None
Std. errors clustered	No	No	No
Observations	92	46	46
R2	0.102	0.077	0.172
Adjusted R2	0.092	0.056	0.153
Note:	*p<0.1; **p<0.05; ***p<0.01

t3 <- prepare_regression_table(df100, dvs, idvs, feffects, clusters, byvar="year2")
htmltools::HTML(t3$table)

	Dependent variable:
	DEEI
	Full Sample	1992	2008
	(1)	(2)	(3)
CEEI	0.256***	0.187***	0.319***
	(0.061)	(0.062)	(0.081)
Estimator	ols	ols	ols
Fixed effects	year2	year2	year2
Std. errors clustered	region	region	region
Observations	92	46	46
R2	0.145	0.077	0.172
Adjusted R2	0.126	0.056	0.153
Note:	*p<0.1; **p<0.05; ***p<0.01

feffects2 <-  c( "year2")

t3 <- prepare_regression_table(df100, dvs, idvs, feffects, clusters, byvar="year2")
htmltools::HTML(t3$table)