Testing the Li Keqiang Index: Analysis and Discussion

To Recap¶

The previous blog post introduced the project and performed an EDA on data related to the Li Keqiang index which raised some important concerns. Because of our EDA, we chose to perform our analysis using the logged forms of our main variables of interest: GDP, rail cargo volume, and electricity consumption. (Bank loan volume had been omitted from the analysis because no reliable data exists.) We also noticed that there was autocorrelation within countries, which would have to be dealt with. Finally, we noticed problems with our data set in general, which will be discussed at the end of this analysis.

Model Building¶

The goal of this analysis was to see if rail cargo volume and electricity consumption are good estimators of GDP. I also wanted to see if this would be impacted by economic variables such as the percentage of the population in urban areas and the percentage of the economy made up of the service sector. Additionally, I wanted to see if this would be impacted by the strength of democracy in the country of interest. I expected countries with lower levels of democracy to be more likely to overreport their GDP. To this end, I built two models: one with just the variables of interest and one that included interactions between my control variables and the economic indicators. Using the statsmodels formula API, the first model is:

log_GDP ~ log_energy_cons+log_rail_cargo

And the second model is:

log_GDP ~ (log_energy_cons+log_rail_cargo)*(urban_per+service_per+democracy)-(urban_per+service_per+democracy)

urban_per, service_per and democracy are removed from the model because I am only interested in how they change the effect of the two economic indicators.

Model Fitting¶

Before we begin, you can follow along by downloading the code for this part of the project here.

Preparing Data and Packages¶

First, I imported the packages we'll need for the analysis:

In [1]:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import math

import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.compat import lzip
import statsmodels.stats.api as sms

from statsmodels.stats.diagnostic import het_breuschpagan

import seaborn as sns

Then, I imported the data we produced in the previous blog post:

In [2]:

df = pd.read_csv('keqiang.csv')
# Some of our names kinda suck so let's fix them
df.rename(columns={'Democracy Index': 'democracy'},inplace = True)
# this makes all of our quantitative vars lowercase, our qualitative ones capitalized, and y var in ALL CAPS

# as mentioned in EDA post, we also need to logarithm energy_cons, rail_cargo, and GDP.
df['GDP'] = [math.log(item) for item in df['GDP']]
df['rail_cargo'] = [math.log(item) for item in df['rail_cargo']]
df['energy_cons'] = [math.log(item) for item in df['energy_cons']]
df.rename(columns={'energy_cons': 'log_energy_cons',
                   'rail_cargo': 'log_rail_cargo',
                   'GDP': 'log_GDP'},inplace = True)
df.head()

Out[2]:

	Country	Year	Region	democracy	log_energy_cons	log_rail_cargo	urban_per	service_per	log_GDP
0	Canada	2019	North America	9.22	6.304906	13.008256	81.482	67.671163	28.187045
1	United States	2019	North America	7.96	8.291391	14.672950	82.459	77.181114	30.693523
2	Austria	2019	Western Europe	8.29	4.228670	9.986725	58.515	63.129879	26.820488
3	Finland	2019	Western Europe	9.25	4.431378	9.236982	85.446	60.194782	26.316172
4	France	2019	Western Europe	8.12	6.104842	10.430967	80.709	70.043587	28.634909

Ordinary Least Squares¶

First, I wanted to try an ordinary least squares model to see if there was anything I needed to fix. I started with a model that just looks at the variables of interest:

In [3]:

#fit regression model
fit = smf.ols('log_GDP ~ log_energy_cons+log_rail_cargo', data=df).fit()

#view model summary
print(fit.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                log_GDP   R-squared:                       0.882
Model:                            OLS   Adj. R-squared:                  0.881
Method:                 Least Squares   F-statistic:                     3046.
Date:                Sun, 17 Dec 2023   Prob (F-statistic):               0.00
Time:                        18:33:04   Log-Likelihood:                -829.49
No. Observations:                 820   AIC:                             1665.
Df Residuals:                     817   BIC:                             1679.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
===================================================================================
                      coef    std err          t      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------
Intercept          22.0179      0.090    244.318      0.000      21.841      22.195
log_energy_cons     1.0960      0.020     55.939      0.000       1.058       1.134
log_rail_cargo     -0.0542      0.014     -3.808      0.000      -0.082      -0.026
==============================================================================
Omnibus:                       21.372   Durbin-Watson:                   1.409
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               22.832
Skew:                          -0.365   Prob(JB):                     1.10e-05
Kurtosis:                       3.367   Cond. No.                         39.4
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

One thing I notice here is that while both energy consumption and rail cargo volume are significant here, rail cargo actually has a negative effect on GDP rather than a positive one. In the pair plots from the previous post, both rail cargo and energy consumption were positively correlated with GDP. This means that once energy consumption has been accounted for, higher rail cargo volume is actually a very weak indicator of economic success and is actually a negative indicator at that. This makes some sense if we reason through it. In a modern service-oriented economy, increased rail cargo volume indicates that the economy is still largely industrial, and thus not as productive.

Next, I tried a model that uses the control variables as interactions.

In [4]:

#fit regression model
fit = smf.ols("log_GDP ~ (log_energy_cons+log_rail_cargo)*(urban_per+service_per+democracy)-(urban_per+service_per+democracy)", data=df).fit()

#view model summary
print(fit.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                log_GDP   R-squared:                       0.923
Model:                            OLS   Adj. R-squared:                  0.922
Method:                 Least Squares   F-statistic:                     1217.
Date:                Sun, 17 Dec 2023   Prob (F-statistic):               0.00
Time:                        18:33:04   Log-Likelihood:                -652.94
No. Observations:                 820   AIC:                             1324.
Df Residuals:                     811   BIC:                             1366.
Df Model:                           8                                         
Covariance Type:            nonrobust                                         
===============================================================================================
                                  coef    std err          t      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------------------
Intercept                      22.1749      0.076    293.437      0.000      22.027      22.323
log_energy_cons                 0.8845      0.091      9.768      0.000       0.707       1.062
log_rail_cargo                 -0.1811      0.046     -3.937      0.000      -0.271      -0.091
log_energy_cons:urban_per      -0.0066      0.001     -6.163      0.000      -0.009      -0.005
log_energy_cons:service_per     0.0072      0.002      3.762      0.000       0.003       0.011
log_energy_cons:democracy       0.0292      0.012      2.531      0.012       0.007       0.052
log_rail_cargo:urban_per        0.0039      0.001      7.638      0.000       0.003       0.005
log_rail_cargo:service_per     -0.0015      0.001     -1.691      0.091      -0.003       0.000
log_rail_cargo:democracy       -0.0057      0.006     -1.035      0.301      -0.017       0.005
==============================================================================
Omnibus:                       16.319   Durbin-Watson:                   1.711
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               22.052
Skew:                          -0.212   Prob(JB):                     1.63e-05
Kurtosis:                       3.682   Cond. No.                     4.92e+03
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 4.92e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

In this model, the pattern remains, and in fact rail cargo volume is an even larger detriment to economy size than before. The only variables not significant at a p <0.05 level are the interactions that rail cargo volume has with percent of GDP from services and level of democracy. This is counterintuitive. We would expect the GDP in economies with higher contributions from services to be noticably harder to predict using rail cargo volume; that is, we would expect the interaction to make the effect of rail cargo volume smaller. However, this is not the case; in fact, the coefficient suggests that if anything, it is making the negative effect of rail cargo volume greater. It seems that in large service industries, higher rail cargo volumne could indicate that more money is being spent on transportation of goods (a low value-add) which is a detriment to the GDP.

Also unusually, the results seem to suggest that weaker democracies and authoritarian states are not inflating their GDP. If they were, we would expect the cumulative effect of the democracy interactions to be negative---that is, in democracies, the same level of energy consumption would predict a smaller economy than it would in an authoritarian country. However, the opposite is true: a 1-point increase in Democracy Index results in a 0.0235 increase in log(GDP). This suggests that democracies themselves are either inflating their GDP (not likely, because parties outside power have an incentive to point out false economic statistics to delegitimize their opponents) or they have intrinsically stronger economies.

Our model accounts for 92.2% of the variability in our data, which is very good. The lower AIC of this model suggests that the interaction variables are a significant improvement over the previous model which just had the two variables of interest.

Let's look at the a plot of residuals and fitted values to see if there are likely to be any problems we need to resolve.

In [5]:

df['fitted'] = fit.fittedvalues
df['resid'] = fit.resid
sns.lmplot(data = df, x = 'fitted', y = 'resid')
plt.show()

The fit looks linear, and our correlation is very, very high, which is good. However, there is also heteroskedasticity---that is, the variance of our residuals is not constant across our fitted values. This is a problem because it means that our standard errors will be incorrect. Since we've already done the transforms that I would have tried for resolving this, we need to fit a model that accounts for this. Let's also check for normality.

In [6]:

# Histogram of residuals
sns.histplot(x = fit.resid)
plt.show()

In [7]:

# QQ plot
sm.qqplot(fit.resid, line='s')
plt.show()

I'm comfortable calling this good enough to say that the residuals are normally distributed.

Mixed Linear Model¶

As mentioned before, I wanted to account for the correlation within countries. A mixed linear model allows us to do a longitudinal multiple linear regression on the data. It takes into account clusters that form when observations are taken multiple times from the same subject. We can tell that different groups formed different clusters from the pairplot in the EDA, which shows the distribution of variable values differed by region. Because they differed by region, it is reasonable to assume that they differed by country as well, and this is the group that we are going to use here. Year also needs to be included as a variable to account for autocorrelation over time here.

First, I created a model with just the variables of interest:

In [8]:

lmm_inter = smf.mixedlm("log_GDP ~ Year +(log_energy_cons+log_rail_cargo)", df, groups=df["Country"],
                        re_formula="~1").fit()

print(lmm_inter.summary())

           Mixed Linear Model Regression Results
===========================================================
Model:              MixedLM   Dependent Variable:   log_GDP
No. Observations:   820       Method:               REML   
No. Groups:         93        Scale:                0.0299 
Min. group size:    1         Log-Likelihood:       45.2410
Max. group size:    12        Converged:            Yes    
Mean group size:    8.8                                    
-----------------------------------------------------------
                Coef.  Std.Err.   z    P>|z|  [0.025 0.975]
-----------------------------------------------------------
Intercept       -8.966    3.638 -2.465 0.014 -16.097 -1.836
Year             0.016    0.002  8.578 0.000   0.012  0.019
log_energy_cons  0.821    0.038 21.786 0.000   0.747  0.895
log_rail_cargo   0.032    0.020  1.600 0.110  -0.007  0.070
Group Var        0.484    0.479                            
===========================================================

Once we account for the differences between countries, increased rail cargo appears to have a positive correlation with GDP. However, its effect is also no longer significant. Let's look at a model that has the interactions as well:

In [9]:

lmm_inter = smf.mixedlm("log_GDP ~ Year + (log_energy_cons+log_rail_cargo)*(urban_per+service_per+democracy)-(urban_per+service_per+democracy)", df, groups=df["Country"],
                        re_formula="~1").fit()

print(lmm_inter.summary())

                 Mixed Linear Model Regression Results
=======================================================================
Model:                  MixedLM       Dependent Variable:       log_GDP
No. Observations:       820           Method:                   REML   
No. Groups:             93            Scale:                    0.0287 
Min. group size:        1             Log-Likelihood:           39.3498
Max. group size:        12            Converged:                Yes    
Mean group size:        8.8                                            
-----------------------------------------------------------------------
                            Coef.  Std.Err.   z    P>|z|  [0.025 0.975]
-----------------------------------------------------------------------
Intercept                   -6.317    4.014 -1.574 0.116 -14.185  1.551
Year                         0.014    0.002  7.164 0.000   0.010  0.018
log_energy_cons              0.623    0.103  6.037 0.000   0.421  0.826
log_rail_cargo              -0.013    0.050 -0.254 0.800  -0.110  0.085
log_energy_cons:urban_per   -0.006    0.002 -3.693 0.000  -0.009 -0.003
log_energy_cons:service_per  0.003    0.002  2.010 0.044   0.000  0.007
log_energy_cons:democracy    0.066    0.013  4.947 0.000   0.040  0.093
log_rail_cargo:urban_per     0.004    0.001  5.141 0.000   0.003  0.006
log_rail_cargo:service_per  -0.001    0.001 -2.089 0.037  -0.003 -0.000
log_rail_cargo:democracy    -0.025    0.006 -4.167 0.000  -0.037 -0.013
Group Var                    0.393    0.411                            
=======================================================================

Once country-level effects are accounted for, the effect of rail cargo is found to be insignificant. Democracy still has a cumulative positive correlation with the GDP, but it is much smaller. Once again, we'll go through the regression assumptions.

The data appear to be normal enough based on a histogram of the residuals and the Q-Q plot:

In [10]:

sns.histplot(lmm_inter.resid)
plt.show()

In [11]:

# QQ plot
sm.qqplot(lmm_inter.resid, line='s')
plt.show()

There is a somewhat long left tail, but we have enough observations that I don't think it should really matter that much.

More concerning is that, according to a fitted values versus residual plot, our data is still heteroskedastic.

In [12]:

df['fitted'] = lmm_inter.fittedvalues
df['resid'] = lmm_inter.resid
sns.lmplot(data = df, x = 'fitted', y = 'resid')
plt.show()

We can confirm this with a Breusch-Pagan test:

In [13]:

# Conduct the Breusch-Pagan test
names = ['Lagrange multiplier statistic', 'p-value',
         'f-value', 'f p-value']
 
# Get the test result
test_result = sms.het_breuschpagan(lmm_inter.resid, lmm_inter.model.exog)
 
lzip(names, test_result)

Out[13]:

[('Lagrange multiplier statistic', 96.27687606428948),
 ('p-value', 8.889448381522433e-17),
 ('f-value', 11.972698618036382),
 ('f p-value', 7.474758802932374e-18)]

With p <0.05, we can reject the null hypothesis and conclude that there is still heteroskedasticity in the data. Unfortunately, the mixedlm() function in statsmodels doesn't have a way to deal with this, despite the fact that statsmodels has provided the ols() function with a way to do so. While our estimators will not be biased by this, this means our standard errors will be incorrect.

Conclusion¶

While there remain some problems with our data, we are still able to draw some conclusions. The first is that while energy consumption has a clear correlation with GDP, rail cargo volume does not seem to add much information to a model that already includes energy consumption and may in fact have a negative correlation with GDP. This suggests that the Li Keqiang index is a flawed measure for the modern economy. While rail cargo volume seems to indicate that more goods are being produced and transported, it also shows that the economy relies on industrial production rather than a strong service or financial industry that produces more value-add.

Additionally, in countries with higher levels of democracy, energy consumption had an even higher effect on GDP than in countries with low levesl of democracy. Since I had expected non-democracies to be overreporting their GDP, I had expected energy consumption to have a lower effect on GDP in countries with higher levels of democracy. As mentioned above, this may suggest that democracies have intrinsically better economies or that they use their resources (such as electricity) for projects that result in stronger development. While percentage of GDP from services and percentage of population in urban areas both had statistically significant interactions with the variables of interest, their effect size was overall very small compared to that of the interactions with levels of democracy. This suggests that they had little impact on the efficacy of the index.

Discussion¶

Throughout this project, I have raised some issues with the data and with the models being used.

Data Issues¶

In the end, most of the data used was from Europe alone. This biases the results of the project toward the behavior of Western economies. This is problematic because European economies are generally more developed and ship things shorter distances. The GDP of less developed and larger countries, such as those of Asia and Africa, would a priori, be believed to be better predicted by rail cargo volume. This means that the effect of rail cargo volume may not be detectable, simply because the kinds of economies that rely on it have not been examined adequately. Additionally, many of the non-European countries in the analysis had missing values. Consider the following graph from the previous blog post:

In [14]:

vc = df['Country'].value_counts() < 12
vc = vc[vc]

dropped_countries = df.loc[df['Country'].isin(vc.index)][['Country','Region']].drop_duplicates()
dropped_countries['Region'].value_counts()
ax = sns.countplot(data = dropped_countries, y = 'Region')
ax.set_title('Countries with missing values by region')

Out[14]:

Text(0.5, 1.0, 'Countries with missing values by region')

Non-European countries, especially in Africa, were the countries which were most likely to have some data missing. This means that the countries whose economies were a priori believed to be best represented by measures of rail cargo volume were also the most likely to be themselves inadequately represented. This is especially problematic when you look at countries that have fewer than four data points:

In [15]:

vc = df['Country'].value_counts() < 4
vc = vc[vc]

dropped_countries = df.loc[df['Country'].isin(vc.index)][['Country','Region']].drop_duplicates()
dropped_countries['Region'].value_counts()
ax = sns.countplot(data = dropped_countries, y = 'Region')
ax.set_title('Countries with less than four values by region')

Out[15]:

Text(0.5, 1.0, 'Countries with less than four values by region')

Africa and the Middle East are extremely overrepresented among these countries, and in fact most of the African and Middle Eastern countries with missing values have fewer than four total values. Because three is the minimum number of data points needed to create a linear regression, this suggests that the data simply does not exist to make an adequate conclusion based on these countries. Meanwhile, Europe and North America are represented strongly in the data and have the most influence on the data results show. This has data ethics implications: is it really ethical to make conclusions about economic indicators like the Li Keqiang Index based on data that prioritizes wealthy, Western economies? It's possible that the woefully inadequate data being gathered on the Global South is leading to misconceptions about those economies.

Model Issues¶

As I mentioned above, the mixedlm() function in the statsmodels package for Python doesn't appear to have any way to deal with heteroskedasticity, despite this presumably being a relatively common problem for statisticians using mixed linear models to face. This means that the standard errors of this analysis, thus decisions about whether a particular variable in a model is statistically significant or not, are skewed. Additionally, one glaring omission from this model is the third element of the Li Keqiang index: bank loan volume. As far as I could tell, there is no public data on this for any country in the world, and it would be extremely difficult to get it in most countries. Li Keqiang could do so because he had developed his economic indicators while he was the effective governor of his province in a country where the banks are directly controlled by the government. Most other researchers, however, do not have such a luxury.

All of this, combined with the data issues above, limits the usability of the mixed linear model for evaluating the Li Keqiang index.