Code
thesis <- read.csv("data.csv")prelims
1 Loading the Data
2 Preliminary Data Diagnostics
Code
# Load necessary libraries
library(dplyr)
library(stargazer)
# Summary of the dataset
stargazer(thesis[, !(names(thesis) %in% c("Country", "Year", "GDPpc2"))],
type = "text", title = "Descriptive Statistics",
summary.stat = c("mean", "sd", "min", "max", "median"), digits = 3)
Descriptive Statistics
=============================================================
Statistic Mean St. Dev. Min Max Median
-------------------------------------------------------------
RDef -0.122 0.692 -2.246 1.388 0.012
GDPpc 11,090.080 6,525.915 3,106.582 29,750.270 9,216.444
PopD 193.390 96.319 73.360 375.897 139.287
AgriLU 34.089 8.153 21.420 45.552 34.373
RENEW 25.257 13.138 2.000 52.200 26.900
Trade 108.709 46.254 32.972 210.374 120.842
CCor -0.420 0.347 -1.137 0.397 -0.483
RegQ -0.041 0.420 -0.866 0.799 0.008
RLaw -0.217 0.401 -0.910 0.573 -0.339
-------------------------------------------------------------Rate of Deforestation (RDef): The mean value of -0.122 percent suggests that, on average, deforestation rates are slightly negative, potentially indicating forest recovery in some areas. However, the variability (standard deviation = 0.692) and the range from -2.246 to 1.388 show that some regions experience significant forest loss while others might be reforesting.
Gross Domestic Product per capita (GDP per capita): The average GDP per capita is around $11,090, but it varies widely across observations, as shown by the high standard deviation ($6,525). The large range (from $3,106 to nearly $29,750) suggests that some countries are wealthier while others remain relatively poor.
Population Density (PopD): With a mean of 193.39 people per unit area, there is a noticeable difference in how densely populated the regions are, ranging from 73.36 to 375.90. Higher population density could exert pressure on natural resources, influencing deforestation.
Agricultural Land Use (AgriLU): The average agricultural land use is about 34.09 percent of the total land area, with moderate variability. This suggests that some areas rely heavily on agriculture, which may contribute to deforestation if expansion occurs at the cost of forests.
Renewable Energy (RENEW): The average share of renewable energy in total consumption is 25.26 percent of the total TFEC, with values ranging from 2.00 to 52.20. This means that some countries have made significant investments in renewable energy, while others still rely heavily on non-renewable sources.
Trade Openness (Trade): With a mean of 108.71 and a wide range (32.97 to 210.37), some countries are highly engaged in international trade while others are more closed.
Institutional Quality (CCor, RegQ, RLaw):
Control of Corruption (CCor): The negative mean (-0.420) suggests that corruption is a prevalent issue in many countries, potentially leading to weak enforcement of environmental policies.
Regulatory Quality (RegQ): The average value close to zero (-0.041) implies a mixed level of regulatory effectiveness, with some countries having strong regulations while others struggle.
Rule of Law (RLaw): The negative mean (-0.217) indicates weaker rule of law in many cases, which could contribute to illegal deforestation and weak environmental governance.
3 Correlation and VIF Checking
Key Observations:
Institutional Variables are highly correlated (CCor, RegQ, RLaw).
GDP per capita has a moderately high correlation with Institutional Variables suggesting that wealthier ASEAN nations will mostly has better institutional regulations.
GDPpc and RENEW has a strong negative correlation which may indicate that wealthier economies tend to rely less on renewable energy, possibly due to continued dependence on fossil fuels.
Trade openness (Trade) is positively correlated with GDP per capita (0.41), suggesting that wealthier economies are more open to trade.
Code
#correlation and VIF check
library(dplyr)
library(ggcorrplot)
numerical_vars <- thesis %>% select(RDef, GDPpc, PopD, AgriLU, RENEW, Trade, CCor, RegQ, RLaw)
cor_matrix <- cor(numerical_vars, use = "complete.obs")
print(cor_matrix) RDef GDPpc PopD AgriLU RENEW Trade
RDef 1.00000000 0.2088263 -0.3921006 -0.2642595 -0.2813786 -0.3290585
GDPpc 0.20882631 1.0000000 -0.6112429 -0.2151636 -0.8319080 0.4082120
PopD -0.39210060 -0.6112429 1.0000000 0.5165452 0.5577752 -0.2603324
AgriLU -0.26425946 -0.2151636 0.5165452 1.0000000 0.2766443 -0.1241431
RENEW -0.28137862 -0.8319080 0.5577752 0.2766443 1.0000000 -0.5395182
Trade -0.32905849 0.4082120 -0.2603324 -0.1241431 -0.5395182 1.0000000
CCor 0.05021369 0.7767199 -0.5055260 -0.3433025 -0.8919754 0.6156180
RegQ 0.33256327 0.7701194 -0.4940333 -0.1290193 -0.8966702 0.3168260
RLaw 0.08677284 0.8262877 -0.5438441 -0.2229963 -0.9153047 0.6826587
CCor RegQ RLaw
RDef 0.05021369 0.3325633 0.08677284
GDPpc 0.77671994 0.7701194 0.82628768
PopD -0.50552602 -0.4940333 -0.54384411
AgriLU -0.34330252 -0.1290193 -0.22299631
RENEW -0.89197535 -0.8966702 -0.91530471
Trade 0.61561795 0.3168260 0.68265866
CCor 1.00000000 0.8322854 0.89616950
RegQ 0.83228544 1.0000000 0.80674692
RLaw 0.89616950 0.8067469 1.00000000Code
ggcorrplot(cor_matrix, method = "circle", type = "upper", lab = TRUE, lab_size = 3, colors = c("red", "white", "green")) +
labs(title = "Correlation Matrix of Numerical Variables",
subtitle = "Panel Dataset 2002-2020") +
theme_minimal()
Code
library(car)
VIFcheck <- lm(RDef ~ GDPpc + PopD + AgriLU + RENEW + Trade + CCor + RegQ + RLaw, data = thesis)
VIFvalues <- vif(VIFcheck)
print(VIFvalues) GDPpc PopD AgriLU RENEW Trade CCor RegQ RLaw
4.531929 2.333785 2.168849 13.376456 3.583514 9.478206 10.207252 12.021071 Code
model2 <- lm(RDef ~ GDPpc + PopD + AgriLU + Trade + RegQ + RLaw + CCor, data = thesis)
summary(model2)$adj.r.squared[1] 0.48324684 Transforming Variables Using PCA
Since the 3 institutional variables are highly correlated, Principal Components Analysis is used to transform these variables into a single policy index names ‘Gov_Index’. Moreover, RENEW was omitted due to having high correlation with ‘Gov_Index’ which may indicate that those countries that has good institutional regulations tend to invest more on sustainable energy resources (subject for verification from related literature).
Code
library(psych)
# Select governance variables
governance <- thesis[, c("CCor", "RegQ", "RLaw")]
# Standardize data before PCA
governance_scaled <- scale(governance)
# Perform PCA (extracting one component)
pca_result <- principal(governance_scaled, nfactors = 1, rotate = "none")
# Create a governance index
thesis$Gov_Index <- pca_result$scores[,1]
# Drop original governance variables (w/ GDPPpc2 and RENEW)
thesis <- thesis[, !(names(thesis) %in% c("GDPpc2", "RENEW", "CCor", "RegQ", "RLaw"))]5 Rechecking of Correlation and VIF
Code
library(car)
library(corrplot)
#check again for VIF
VIFcheck <- lm(RDef ~ GDPpc + PopD +AgriLU + Trade + Gov_Index, data = thesis)
vif(VIFcheck) GDPpc PopD AgriLU Trade Gov_Index
3.929155 2.124442 1.408658 1.518775 4.246384 Code
library(dplyr)
library(ggplot2)
library(corrplot)
library(ggcorrplot)
independent_vars <- thesis %>% select(GDPpc,PopD, AgriLU, Trade, Gov_Index)
cor_matrix <- cor(thesis %>% select(-Country, -Year), use = "complete.obs")
print(cor_matrix) #correlation matrix RDef GDPpc PopD AgriLU Trade Gov_Index
RDef 1.0000000 0.2088263 -0.3921006 -0.2642595 -0.3290585 0.1632786
GDPpc 0.2088263 1.0000000 -0.6112429 -0.2151636 0.4082120 0.8353195
PopD -0.3921006 -0.6112429 1.0000000 0.5165452 -0.2603324 -0.5433278
AgriLU -0.2642595 -0.2151636 0.5165452 1.0000000 -0.1241431 -0.2460225
Trade -0.3290585 0.4082120 -0.2603324 -0.1241431 1.0000000 0.5706921
Gov_Index 0.1632786 0.8353195 -0.5433278 -0.2460225 0.5706921 1.0000000Code
ggcorrplot(cor_matrix, method = "circle", type = "upper", lab = TRUE, lab_size = 3, colors = c("green", "white", "red")) +
labs(title = "Correlation Matrix of Variables",
subtitle = "Panel Dataset from 2002-2020") +
theme_minimal()
Since GDPpc and Gov_Index has a moderately high correlation (0.84), it may be practical to examine both variables’ interaction effect, meaning there might be a possibility that the effect of GDPpc to RDef may be affected by the Gov_Index and vice versa.
6 Converting Data into Panel Data
Code
library(plm)
library(psych)
library(stargazer)
attach(thesis)
Y <- cbind(RDef)
X <- cbind(GDPpc, PopD, AgriLU, Trade, Gov_Index)
paneldata <- pdata.frame(thesis, index = c("Country", "Year"))
pdim(paneldata)Balanced Panel: n = 5, T = 19, N = 95Code
stargazer(paneldata, type = "text", title = "Descriptive Statistics of Panel Data", digits = 3)
Descriptive Statistics of Panel Data
======================================================
Statistic N Mean St. Dev. Min Max
------------------------------------------------------
RDef 95 -0.122 0.692 -2.246 1.388
GDPpc 95 11,090.080 6,525.915 3,106.582 29,750.270
PopD 95 193.390 96.319 73.360 375.897
AgriLU 95 34.089 8.153 21.420 45.552
Trade 95 108.709 46.254 32.972 210.374
Gov_Index 95 -0.000 1.000 -1.927 2.186
------------------------------------------------------7 Stationarity/Unit Root Test
Augmented Dicky-Fuller (ADF) Test is utilized to examine the stationarity of the data, which may result to spurious regression inferences when not treated properly. Based on the ADF Test, all variables are stationary meaning all variables has a constant statistical properties.
Code
library(plm)
library(tseries)
library(urca)
ADFTest_GDPpc <- adf.test(paneldata$GDPpc, alternative = "stationary")
print(ADFTest_GDPpc)
Augmented Dickey-Fuller Test
data: paneldata$GDPpc
Dickey-Fuller = -2.7269, Lag order = 4, p-value = 0.2762
alternative hypothesis: stationaryCode
ADFTest_RDef <- adf.test(paneldata$RDef, alternative = "stationary")
print(ADFTest_RDef)
Augmented Dickey-Fuller Test
data: paneldata$RDef
Dickey-Fuller = -3.4961, Lag order = 4, p-value = 0.04645
alternative hypothesis: stationaryCode
ADFTest_PopD <- adf.test(paneldata$PopD, alternative = "stationary")
print(ADFTest_PopD)
Augmented Dickey-Fuller Test
data: paneldata$PopD
Dickey-Fuller = -2.1902, Lag order = 4, p-value = 0.4978
alternative hypothesis: stationaryCode
ADFTest_Trade <- adf.test(paneldata$Trade, alternative = "stationary")
print(ADFTest_Trade)
Augmented Dickey-Fuller Test
data: paneldata$Trade
Dickey-Fuller = -2.6876, Lag order = 4, p-value = 0.2924
alternative hypothesis: stationaryCode
ADFTest__AgriLU <- adf.test(paneldata$AgriLU, alternative = "stationary")
print(ADFTest__AgriLU)
Augmented Dickey-Fuller Test
data: paneldata$AgriLU
Dickey-Fuller = -1.6466, Lag order = 4, p-value = 0.7224
alternative hypothesis: stationaryCode
ADFTest_Gov <- adf.test(paneldata$Gov_Index, alternative = "stationary")
print(ADFTest_Gov)
Augmented Dickey-Fuller Test
data: paneldata$Gov_Index
Dickey-Fuller = -2.5806, Lag order = 4, p-value = 0.3366
alternative hypothesis: stationaryCode
#<<all are stationary>>#8 Exploratory Descriptive Analysis
Code
library(ggplot2)
#logging GDPpc and RDef
log <- thesis %>% mutate(log_RDef = RDef, log_GDPpc = log(GDPpc))
thesis$RDef[is.na(thesis$RDef)] <- mean(thesis$RDef, na.rm = TRUE)
thesis$GDPpc[is.na(thesis$GDPpc)] <- mean(thesis$GDPpc, na.rm = TRUE)
#variables' visual presentation over time per country
ggplot(log, aes(x = Year, y = log_GDPpc, group = Country, color = Country)) +
geom_line() + theme_minimal() + ggtitle("GDP per Capita Over Time")
Code
ggplot(log, aes(x = Year, y = log_RDef, group = Country, color = Country)) +
geom_line() + theme_minimal() + ggtitle("Rate of Deforestation per Capita Over Time")
Code
ggplot(log, aes(x = Year, y = PopD, group = Country, color = Country)) +
geom_line() + theme_minimal() + ggtitle("Population Density Over Time")
Code
ggplot(log, aes(x = Year, y = Trade, group = Country, color = Country)) +
geom_line() + theme_minimal() + ggtitle("Trade Over Time")
Code
ggplot(log, aes(x = Year, y = AgriLU, group = Country, color = Country)) +
geom_line() + theme_minimal() + ggtitle("Agricultural Land Use Over Time")
Code
ggplot(log, aes(x = Year, y = Gov_Index, group = Country, color = Country)) +
geom_line() + theme_minimal() + ggtitle("Governance Index Over Time")
Code
library(ggplot2)
library(dplyr)
library(scales)
ggplot(log, aes(x = log_GDPpc, y = log_RDef)) +
geom_point() +
geom_smooth(method = "loess", se = FALSE) + # Adding a smooth line to visualize the relationship
facet_wrap(~ Country, scales = "free") + # Create separate scatter plots for each country with free x-axis scale
labs(title = "Scatter Plot of log(RDef) vs. log(GDPpc) by Country",
x = "log(GDP per capita)",
y = "log(RDef)") +
theme_minimal()
GDP per capita over time:
- Malaysia consistently has the highest GDP per capita among the five countries, showing a steady upward trend.
- Thailand follows with a similar pattern but at a lower level.
- Indonesia and Vietnam demonstrate consistent growth, with Vietnam catching up to Indonesia by the end of the period.
- The Philippines also shows a positive trend but remains lower than Thailand and Malaysia.
Rate of Deforestation over time:
Most countries maintained relatively stable deforestation rates before 2010.
There are noticeable abrupt changes in the rate of deforestation for some countries, particularly Indonesia, Malaysia, and Vietnam around 2010-2015.
Vietnam saw a sharp decline in deforestation rates around 2015, followed by a slight recovery.
- Thailand (Blue Line) and the Philippines (Green Line) show more stability. These two countries have relatively consistent deforestation trends compared to the others.
Population Density over time:
- Philippines has the highest population density, with a continuous upward trend over time.
Vietnam also shows a steady increase in population density, ranking second among the five countries.
Indonesia (red line) and Thailand (blue line) have similar population densities, with slow but steady growth.
Malaysia (yellow line) has the lowest population density among the five countries, though it is also increasing over time.
Trade over time:
- Malaysia had the highest trade levels in the early 2000s, peaking around 2005, but experienced a continuous decline after that.
- Vietnam initially had lower trade levels but showed a significant upward trend after 2010, surpassing Thailand and approaching Malaysia by 2020.
Thailand Maintains a Relatively Stable Trade Pattern: Thailand fluctuated but remained at a consistent level until a drop after 2015.
The Philippines shows a downward trend after 2005, with some fluctuations.
Indonesia has the lowest trade levels among the five countries, with a mostly stable but slightly declining trend.
Agricultural Land Use over time
- Thailand (blue line) consistently had the highest agricultural land use, with a steady increase over time, reaching around 45% by 2020.
- The Philippines (green line) had the second-highest agricultural land use, showing a gradual increase from about 38% in 2002 to over 40% by 2020.
- Vietnam (pink line) started with lower agricultural land use but had a significant jump around 2010-2015, stabilizing afterward.
- Indonesia (red line) started at a lower percentage but displayed a gradual upward trend, reaching about 28% by 2020.
- Malaysia (yellow line) had the lowest agricultural land use among the five countries. However, it experienced gradual growth, with a more noticeable increase after 2010.
Governance Index over time:
Higher values indicate better governance.
- Malaysia consistently had the highest governance index, maintaining positive values throughout the period, though it experienced some fluctuations, particularly a decline around 2010, followed by a recovery.
- Thailand maintained relatively stable governance scores but showed a slight downward trend over time.
- Indonesia started with the lowest governance index (around -2) in 2002 but steadily improved, reaching positive values by 2020.
- The Philippines and Vietnam showed a slow but steady improvement, with governance scores gradually increasing over time.
EKC Hypothesis Validation:
Indonesia follows an a not-so clear inverted U-shape relationship between the relationship of rate of deforestation and GDP per capita.
Malaysia shows a U-shaped curve, indicating an initial decline in deforestation as GDP per capita increases, but later, deforestation rises again.
Philippines shows a more evident inverted U-shape relationship between rate of deforestation and GDP per capita.
Thailand shows a clear S-shaped pattern, with deforestation declining sharply as GDP per capita rises beyond a critical point.
Vietnam appears to have a gradual increase in deforestation with GDP per capita, followed by a leveling off.
9 1st Model
Pooled OLS
Code
Pooled <- plm(RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU,
data = paneldata, index = c("Country", "Year"),
model = "pooling")
summary(Pooled) Pooling Model
Call:
plm(formula = RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU,
data = paneldata, model = "pooling", index = c("Country",
"Year"))
Balanced Panel: n = 5, T = 19, N = 95
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-1.53180 -0.33906 -0.14633 0.33222 1.19871
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.8778e+00 3.7594e-01 4.9949 2.959e-06 ***
GDPpc 1.7779e-05 2.1718e-05 0.8186 0.415222
Gov_Index 4.7690e-01 1.5979e-01 2.9845 0.003676 **
PopD -1.8241e-03 9.3040e-04 -1.9605 0.053095 .
Trade -9.0955e-03 1.4834e-03 -6.1315 2.411e-08 ***
AgriLU -2.1187e-02 1.0811e-02 -1.9598 0.053189 .
GDPpc:Gov_Index -2.4740e-05 1.2344e-05 -2.0042 0.048124 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 45.038
Residual Sum of Squares: 25.64
R-Squared: 0.43071
Adj. R-Squared: 0.3919
F-statistic: 11.0965 on 6 and 88 DF, p-value: 3.4947e-09Random Effects
Code
Random <- plm(RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU,
data = paneldata, index = c("Country", "Year"),
model = "random", random.method = "walhus")
summary(Random)Oneway (individual) effect Random Effect Model
(Wallace-Hussain's transformation)
Call:
plm(formula = RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU,
data = paneldata, model = "random", random.method = "walhus",
index = c("Country", "Year"))
Balanced Panel: n = 5, T = 19, N = 95
Effects:
var std.dev share
idiosyncratic 0.24224 0.49218 0.898
individual 0.02765 0.16628 0.102
theta: 0.4382
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-1.430429 -0.307696 -0.065979 0.305308 1.038454
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) 1.3737e+00 5.1211e-01 2.6824 0.007310 **
GDPpc 4.1024e-05 2.0824e-05 1.9700 0.048835 *
Gov_Index 2.7612e-01 1.5427e-01 1.7899 0.073477 .
PopD -2.5275e-03 1.2400e-03 -2.0384 0.041515 *
Trade -4.7358e-03 1.9646e-03 -2.4105 0.015929 *
AgriLU -2.2807e-02 1.4716e-02 -1.5498 0.121182
GDPpc:Gov_Index -3.1489e-05 1.1920e-05 -2.6418 0.008247 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 26.949
Residual Sum of Squares: 21.187
R-Squared: 0.21381
Adj. R-Squared: 0.1602
Chisq: 23.9318 on 6 DF, p-value: 0.00053755Fixed Effects
Code
Fixed <- plm(RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU + Trade,
data = paneldata, index = c("Country", "Year"),
model = "within")
summary(Fixed)Oneway (individual) effect Within Model
Call:
plm(formula = RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU +
Trade, data = paneldata, model = "within", index = c("Country",
"Year"))
Balanced Panel: n = 5, T = 19, N = 95
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-1.3075830 -0.1667352 0.0033175 0.1980956 0.7032613
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
GDPpc 1.3678e-04 2.2351e-05 6.1195 2.875e-08 ***
Gov_Index 3.1644e-01 1.3755e-01 2.3006 0.0238908 *
PopD -2.5118e-02 3.8810e-03 -6.4722 6.139e-09 ***
Trade 8.6006e-03 2.9499e-03 2.9156 0.0045507 **
AgriLU 4.7605e-02 4.5666e-02 1.0425 0.3001900
GDPpc:Gov_Index -4.2134e-05 1.1156e-05 -3.7769 0.0002955 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 18.607
Residual Sum of Squares: 9.2867
R-Squared: 0.50091
Adj. R-Squared: 0.4415
F-statistic: 14.0511 on 6 and 84 DF, p-value: 5.2361e-11Breusch-Pagan LM Test and Hausman Test
Code
# Perform Breusch-Pagan LM Test
bp_test <- plmtest(Pooled, type = "bp")
print(bp_test)
Lagrange Multiplier Test - (Breusch-Pagan)
data: RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU
chisq = 8.9726, df = 1, p-value = 0.002741
alternative hypothesis: significant effectsCode
#p-value = 0.002741, random effects is preferred
#Hausman Test
hausman_test <- phtest(Random, Fixed)
print(hausman_test)
Hausman Test
data: RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU
chisq = 51.662, df = 6, p-value = 2.181e-09
alternative hypothesis: one model is inconsistentCode
#p-value = 2.181e-09, fixed effects is better
#Generate Well-Formatted Table with Proper Labels
stargazer(Pooled, Random, Fixed, type = "text",
title = "Panel Regression Results",
column.labels = c("Pooled", "Random", "Fixed"),
digits = 3, align = TRUE)
Panel Regression Results
=========================================================================
Dependent variable:
---------------------------------------------------------
RDef
Pooled Random Fixed
(1) (2) (3)
-------------------------------------------------------------------------
GDPpc 0.00002 0.00004** 0.0001***
(0.00002) (0.00002) (0.00002)
Gov_Index 0.477*** 0.276* 0.316**
(0.160) (0.154) (0.138)
PopD -0.002* -0.003** -0.025***
(0.001) (0.001) (0.004)
Trade -0.009*** -0.005** 0.009***
(0.001) (0.002) (0.003)
AgriLU -0.021* -0.023 0.048
(0.011) (0.015) (0.046)
GDPpc:Gov_Index -0.00002** -0.00003*** -0.00004***
(0.00001) (0.00001) (0.00001)
Constant 1.878*** 1.374***
(0.376) (0.512)
-------------------------------------------------------------------------
Observations 95 95 95
R2 0.431 0.214 0.501
Adjusted R2 0.392 0.160 0.441
F Statistic 11.097*** (df = 6; 88) 23.932*** 14.051*** (df = 6; 84)
=========================================================================
Note: *p<0.1; **p<0.05; ***p<0.01Post Estimation Analysis
Code
# Load necessary libraries
library(car) # For VIF
library(lmtest) # For heteroscedasticity test
library(sandwich) # For robust standard errors
#Residual Analysis
par(mfrow = c(2, 2)) # Set layout for multiple plots
plot(as.numeric(fitted(Fixed)), as.numeric(residuals(Fixed)),
main = "Residuals vs Fitted", xlab = "Fitted Values", ylab = "Residuals")
hist(as.numeric(residuals(Fixed)), main = "Histogram of Residuals", xlab = "Residuals")
qqnorm(as.numeric(residuals(Fixed))); qqline(as.numeric(residuals(Fixed))) # QQ Plot
par(mfrow = c(1, 1)) # Reset layout
Code
#VIF
vif(lm(RDef ~ GDPpc + PopD + Trade + AgriLU + Gov_Index, data = paneldata)) GDPpc PopD Trade AgriLU Gov_Index
3.929155 2.124442 1.518775 1.408658 4.246384 Code
#Heteroscedasticity Test
bptest(Fixed, studentize = FALSE) # Breusch-Pagan Test for heteroscedasticity
Breusch-Pagan test
data: Fixed
BP = 19.003, df = 6, p-value = 0.004158Code
#Heteroscedasticity was detected (Breusch-Pagan test), corrected it using robust standard errors.
#Serial Correlation Test
pbgtest(Fixed) # Wooldridge test for autocorrelation
Breusch-Godfrey/Wooldridge test for serial correlation in panel models
data: RDef ~ GDPpc * Gov_Index + PopD + Trade + AgriLU + Trade
chisq = 27.893, df = 19, p-value = 0.08552
alternative hypothesis: serial correlation in idiosyncratic errorsCode
#No serial correlation was found (Wooldridge test), so time dependence is not a major issue.
#Robust standard errors (clustered at the country level)
coeftest(Fixed, vcov = vcovHC(Fixed, type = "HC1", cluster = "group"))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
GDPpc 1.3678e-04 1.5974e-05 8.5626 4.505e-13 ***
Gov_Index 3.1644e-01 1.1532e-01 2.7442 0.007416 **
PopD -2.5118e-02 2.7694e-03 -9.0700 4.289e-14 ***
Trade 8.6006e-03 2.0796e-03 4.1357 8.351e-05 ***
AgriLU 4.7605e-02 3.7475e-02 1.2703 0.207475
GDPpc:Gov_Index -4.2134e-05 5.9818e-06 -7.0437 4.780e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Code
# Load necessary libraries
library(plm)
library(lmtest)
library(sandwich)
library(stargazer)
# Compute robust standard errors (clustered at the country level)
robust_se <- vcovHC(Fixed, type = "HC1", cluster = "group")
robust_results <- coeftest(Fixed, vcov = robust_se)
# Extract estimates and standard errors
estimates <- robust_results[, 1]
std_errors <- robust_results[, 2]
t_values <- robust_results[, 3]
p_values <- robust_results[, 4]
# Create a well-formatted table using stargazer
stargazer(Fixed, type = "text",
title = "Fixed Effects Model with Robust Standard Errors",
se = list(std_errors), # Apply robust standard errors
t = list(t_values), # Add t-values
p = list(p_values), # Add p-values
digits = 4, align = TRUE, star.cutoffs = c(0.1, 0.05, 0.01))
Fixed Effects Model with Robust Standard Errors
===========================================
Dependent variable:
---------------------------
RDef
-------------------------------------------
GDPpc 0.0001***
(0.00002)
Gov_Index 0.3164***
(0.1153)
PopD -0.0251***
(0.0028)
Trade 0.0086***
(0.0021)
AgriLU 0.0476
(0.0375)
GDPpc:Gov_Index -0.00004***
(0.00001)
-------------------------------------------
Observations 95
R2 0.5009
Adjusted R2 0.4415
F Statistic 14.0511*** (df = 6; 84)
===========================================
Note: *p<0.1; **p<0.05; ***p<0.0110 2nd Model
11 3rd Model
