analysis of market price gdp per capita across european countries

comparing GDP per capita across European countries in 2014, and pulling in a few country-level factors to see what lines up with it.

# .libPaths('/tmp')

# install.packages("eurostat")
# install.packages("ggplot2")
# install.packages("GGally")
# install.packages("car")
# install.packages("lmtest")

library(eurostat)
library(ggplot2)
library(GGally)
library(car)
library(MASS)
library(lmtest)

id <- "nama_10_pc"
df <- get_eurostat(id = id)

hdp <- subset(df, TIME_PERIOD == "2014-01-01")
hdp <- subset(hdp, na_item == "B1GQ")
hdp <- subset(hdp, unit == "CP_EUR_HAB")
hdp <- subset(hdp, select = c("geo", "values"))
hdp <- subset(hdp, !grepl("^EA", geo))
hdp <- subset(hdp, !grepl("^EU", geo))
hdp <- label_eurostat(hdp)

Table nama_10_pc cached at /tmp/Rtmp5jJlAc/eurostat/75bde536dea5c4c637c89baee219dd46.rds

Main GDP aggregates per capita

• Code: nama_10_pc
• Time: 2014
• Time frequency: Annual
• National accounts indicator: Gross domestic product at market prices
• Unit of measure: Current prices, euro per capita

hdp$geo

1. 'Albania'
2. 'Austria'
3. 'Belgium'
4. 'Bulgaria'
5. 'Switzerland'
6. 'Cyprus'
7. 'Czechia'
8. 'Germany'
9. 'Denmark'
10. 'Estonia'
11. 'Greece'
12. 'Spain'
13. 'Finland'
14. 'France'
15. 'Croatia'
16. 'Hungary'
17. 'Ireland'
18. 'Iceland'
19. 'Italy'
20. 'Liechtenstein'
21. 'Lithuania'
22. 'Luxembourg'
23. 'Latvia'
24. 'Montenegro'
25. 'North Macedonia'
26. 'Malta'
27. 'Netherlands'
28. 'Norway'
29. 'Poland'
30. 'Portugal'
31. 'Romania'
32. 'Serbia'
33. 'Sweden'
34. 'Slovenia'
35. 'Slovakia'
36. 'Türkiye'
37. 'United Kingdom'

the Nordic countries, plus Ireland, Switzerland, and Luxembourg, top the GDP per capita list. Western Europe sits high too. the lowest values are in Eastern Europe.

data description

we look at the GDP of European countries in 2014, in euros per capita (CP_EUR_HAB) at market prices (B1GQ).

options(repr.plot.width = 10, repr.plot.height = 6)
plot.ecdf(hdp$values, xlab = "value")

options(repr.plot.width = 10, repr.plot.height = 6)
hist(hdp$values, breaks = 9, xlab = "value")

head(hdp[order(hdp$values, decreasing = T), ], 5)
tail(hdp[order(hdp$values, decreasing = T), ], 5)

options(repr.plot.width = 10, repr.plot.height = 3)
boxplot(hdp$values, horizontal = T, xlab = "value")

cat("Count:", length(hdp$values), "\n")
cat("Var:", var(hdp$values), "\n")
cat("Std:", sd(hdp$values), "\n")
summary(hdp$values)

Count: 37 
Var: 730769837 
Std: 27032.75 



   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3450   10770   20170   29245   38990  133180 

we see that only a few countries have very high GDP per capita, most sit in the 1,000 to 4,000 EUR range.

factors affecting GDP

GDP is shaped by a lot of things, here are some. the bolded ones are the regressors used in the next section.

economy and trade

• Inflation and interest rates influence consumer and investor confidence. low inflation boosts confidence and helps growth.
• Exports and imports. high exports point to a strong, self-sufficient economy. high imports hint at dependence on foreign trade and possible imbalance.

labor market

• Employment. higher employment usually means higher GDP, more people with income and purchasing power.
• Income levels. job quality matters too, not just having one. lots of low-wage jobs doesn’t necessarily lift GDP per capita. people earning more tend to spend more on goods and services.

other factors

• Happiness and satisfaction. a happier population usually performs better at work and tends to help growth.

• Adaptation to technological changes. a country’s ability to absorb new tech boosts competitiveness and productivity, and opens up new opportunities in production, services, and trade.

• Climatic conditions. floods, droughts, or extreme temperatures can hurt agriculture and infrastructure. favorable climates support industries like tourism and agriculture.

indicators

the regressors we use are

• Employment rate
• Average monthly income
• Internet usage
• Climate region

next we look at the key properties of these indicators and how they relate to GDP.

all indicators are statistically significant and explain the modeled variable well. at the standard significance level we reject the hypothesis that the correlation between each regressor and GDP is zero.

id <- "lfsa_ergan"
emp <- get_eurostat(id = id)

emp <- subset(emp, TIME_PERIOD == "2014-01-01")
emp <- subset(emp, sex == "T")
emp <- subset(emp, citizen == "TOTAL")
emp <- subset(emp, age == "Y_GE25")
emp <- subset(emp, !grepl("^EU", geo))
emp <- subset(emp, !grepl("^EA", geo))
emp <- subset(emp, select = c("geo", "values"))
emp <- label_eurostat(emp)

# https://data.worldbank.org/indicator/SL.EMP.TOTL.SP.NE.ZS?locations=LI-AL
emp <- rbind(emp, c("Albania", 43.8))
emp <- rbind(emp, c("Liechtenstein", 59.4))

# ===========================================

id <- "earn_ses14_19"
earn <- get_eurostat(id = id)

earn <- subset(earn, cpayagr == "TOTAL")
earn <- subset(earn, currency == "EUR")
earn <- subset(earn, sizeclas == "GE10")
earn <- subset(earn, TIME_PERIOD == "2014-01-01")
earn <- subset(earn, indic_se == "ERN")
earn <- subset(earn, sex == "T")
earn <- subset(earn, nace_r2 == "B-S_X_O")
earn <- subset(earn, !grepl("^EU", geo))
earn <- subset(earn, !grepl("^EA", geo))
earn <- subset(earn, select = c("geo", "values"))
earn <- label_eurostat(earn)

# https://ndiqparate.al/?p=9775&lang=en
earn <- rbind(earn, c("Albania", 480))
# https://archiv.llv.li/files/as/jahrbuch-2017.pdf (p. 121)
earn <- rbind(earn, c("Liechtenstein", 6692))

# ===========================================

id <- "tin00028"
int <- get_eurostat(id = id)

int <- subset(int, indic_is == "I_IU3")
int <- subset(int, TIME_PERIOD == "2014-01-01")
int <- subset(int, !grepl("^EU", geo))
int <- subset(int, !grepl("^EA", geo))
int <- subset(int, select = c("geo", "values"))
int <- label_eurostat(int)

# https://data.worldbank.org/indicator/IT.NET.USER.ZS?locations=AL-LI-ME-RS
int <- rbind(int, c("Albania", 54))
int <- rbind(int, c("Liechtenstein", 95))
int <- rbind(int, c("Montenegro", 61))
int <- rbind(int, c("Serbia", 62))

# ===========================================

countries <- c(
    "Albania", "Greece", "Italy", "Cyprus", "Croatia", "Spain", "North Macedonia", "Malta",
    "Slovenia", "Serbia", "Türkiye", "Montenegro", "Portugal",
    
    "Austria", "Belgium", "Bulgaria", "Switzerland", "Liechtenstein", "Luxembourg", "Netherlands",
    "Slovakia", "United Kingdom", "Poland", "Romania", "Czechia", "Germany", "France", "Hungary",
    "Ireland",
    
    "Denmark", "Estonia", "Finland", "Iceland", "Latvia", "Lithuania", "Norway", "Sweden"
)
climate_zones <- c(rep("Mediterranean climate", 13), rep("Temperate", 16), rep("Cold", 8))
clima <- data.frame(geo = countries, climate = climate_zones)

# ===========================================

colnames(emp)[colnames(emp) == "values"] <- "employment"
colnames(earn)[colnames(earn) == "values"] <- "earnings"
colnames(int)[colnames(int) == "values"] <- "internet"
colnames(clima)[colnames(clima) == "values"] <- "climate"

df <- merge(hdp, emp, by = "geo")
df <- merge(df, earn, by = "geo")
df <- merge(df, int, by = "geo")
df <- merge(df, clima, by = "geo")

df$employment <- as.double(df$employment)
df$earnings <- as.double(df$earnings)
df$internet <- as.double(df$internet)
df$climate <- as.factor(df$climate)

Table lfsa_ergan cached at /tmp/Rtmp5jJlAc/eurostat/6f0b206586a43436eb71f5af3e1aa2e0.rds

Table earn_ses14_19 cached at /tmp/Rtmp5jJlAc/eurostat/499f0f870ce46c24042e585f1a785759.rds

Table tin00028 cached at /tmp/Rtmp5jJlAc/eurostat/8820f6902667cc62164e76efda7ba3ce.rds

employment rate

• Code: lfsa_ergan
• Name: Employment rates by sex, age and citizenship
• Time: 2014
• Sex: Total
• Age: 25 years and over
• Country of citizenship: Total
• Time frequency: Annual
• Unit of measure: Percentage

head(df[order(df$employment, decreasing = T), c("geo", "employment")], 5)
tail(df[order(df$employment, decreasing = T), c("geo", "employment")], 5)

regressor analysis

cat("Var:", var(df$employment), "\n")
cat("Std:", sd(df$employment), "\n")
summary(df$employment)

Var: 52.02444 
Std: 7.212797 



   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   41.4    51.6    55.8    55.5    59.4    79.8 

options(repr.plot.width = 10, repr.plot.height = 6)
plot.ecdf(df$employment, xlab = "employment [%]")

options(repr.plot.width = 10, repr.plot.height = 6)
hist(df$employment, xlab = "employment [%]")

options(repr.plot.width = 10, repr.plot.height = 3)
boxplot(df$employment, horizontal = T, xlab = "employment [%]")

dependence on GDP

cor.test(df$values, df$employment)

 Pearson's product-moment correlation

data:  df$values and df$employment
t = 3.6397, df = 35, p-value = 0.0008728
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2408721 0.7249318
sample estimates:
      cor 
0.5239948 

fit <- lm(df$values ~ df$employment)
summary(fit)

options(repr.plot.width = 10, repr.plot.height = 6)
plot(df$values ~ df$employment)
abline(fit)

Call:
lm(formula = df$values ~ df$employment)

Residuals:
   Min     1Q Median     3Q    Max 
-35047 -14744  -3677   7849  96276 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -79749.9    30191.4  -2.641 0.012255 *  
df$employment   1963.9      539.6   3.640 0.000873 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 23350 on 35 degrees of freedom
Multiple R-squared:  0.2746, Adjusted R-squared:  0.2538 
F-statistic: 13.25 on 1 and 35 DF,  p-value: 0.0008728

average monthly income

• Code: earn_ses14_19
• Name: Mean monthly earnings by sex, economic activity and collective pay agreement
• Time: 2014
• Sex: Total
• collective pay agreement: Total
• Currency: Euro
• Economic activities: Industry, construction and services
• Indicator: Gross earnings
• Time frequency: Annual

head(df[order(df$earnings, decreasing = T), c("geo", "earnings")], 5)
tail(df[order(df$earnings, decreasing = T), c("geo", "earnings")], 5)

regressor analysis

cat("Var:", var(df$earnings), "\n")
cat("Std:", sd(df$earnings), "\n")
summary(df$earnings)

Var: 2644510 
Std: 1626.195 



   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    431     811    1720    2201    3151    6692 

options(repr.plot.width = 10, repr.plot.height = 6)
plot.ecdf(df$earnings, xlab = "earnings [EUR]")

options(repr.plot.width = 10, repr.plot.height = 6)
hist(df$earnings, xlab = "earnings [EUR]")

options(repr.plot.width = 10, repr.plot.height = 3)
boxplot(df$earnings, horizontal = T, xlab = "earnings [EUR]")

dependence on GDP

cor.test(df$values, df$earnings)

 Pearson's product-moment correlation

data:  df$values and df$earnings
t = 14.946, df = 35, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8669865 0.9635352
sample estimates:
      cor 
0.9298047 

fit <- lm(df$values ~ df$earnings)
summary(fit)

options(repr.plot.width = 10, repr.plot.height = 6)
plot(df$values ~ df$earnings)
abline(fit)

Call:
lm(formula = df$values ~ df$earnings)

Residuals:
   Min     1Q Median     3Q    Max 
-21215  -5599   -857   3008  34519 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -4773.653   2816.531  -1.695    0.099 .  
df$earnings    15.456      1.034  14.946   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10090 on 35 degrees of freedom
Multiple R-squared:  0.8645, Adjusted R-squared:  0.8607 
F-statistic: 223.4 on 1 and 35 DF,  p-value: < 2.2e-16

internet use by individuals

• Code: tin00028
• Name: Internet use by individuals
• Time: 2014
• Individual type: All individuals
• Last internet use: in last 3 months
• Time frequency: Annual
• Unit of measure: Percantage of individuals

regressor analysis

cat("Var:", var(df$internet), "\n")
cat("Std:", sd(df$internet), "\n")
summary(df$internet)

Var: 193.9805 
Std: 13.92769 



   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  48.49   66.60   76.19   76.80   89.73   98.16 

options(repr.plot.width = 10, repr.plot.height = 6)
plot.ecdf(df$internet, xlab = "internet usage [%]")

options(repr.plot.width = 10, repr.plot.height = 6)
hist(df$internet, xlab = "internet usage [%]")

options(repr.plot.width = 10, repr.plot.height = 3)
boxplot(df$internet, horizontal = T, xlab = "internet usage [%]")

dependence on GDP

cor.test(df$values, df$internet)

 Pearson's product-moment correlation

data:  df$values and df$internet
t = 6.1225, df = 35, p-value = 5.327e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5151387 0.8460142
sample estimates:
      cor 
0.7191251 

fit <- lm(df$values ~ df$internet)
summary(fit)

options(repr.plot.width = 10, repr.plot.height = 6)
plot(df$values ~ df$internet)
abline(fit)

Call:
lm(formula = df$values ~ df$internet)

Residuals:
   Min     1Q Median     3Q    Max 
-24393 -12278  -4654   6025  78528 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -77947      17786  -4.383 0.000102 ***
df$internet     1396        228   6.122 5.33e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 19050 on 35 degrees of freedom
Multiple R-squared:  0.5171, Adjusted R-squared:  0.5033 
F-statistic: 37.48 on 1 and 35 DF,  p-value: 5.327e-07

climate area

source: https://www.barenbrug.biz/forage/european-climate-zones

there are four climate zones on the map:

1. Nordic climate (cold winters and mild humid summers)
2. Eastern-continental climate (cold winters and hot summers)
3. Oceanic climate (mild winters and humid summers)
4. Mediterranean climate

for simplicity, merge 2 and 3 into one zone.

regressor analysis

tapply(df$values, df$climate, summary)

$Cold
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  11850   14550   39900   35632   45470   73670 

$`Mediterranean climate`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   3450    5560   16270   13776   20170   26980 

$Temperate
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5960   13260   35865   38620   40368  133180 

options(repr.plot.width = 10, repr.plot.height = 6)
barplot(table(df$climate), xlab = "climate", ylab = "Frequency")

dependence on GDP

fit <- lm(df$values ~ df$climate)
summary(fit)

options(repr.plot.width = 10, repr.plot.height = 8)
boxplot(df$values ~ df$climate, ylab = "value", xlab = "climate")

Call:
lm(formula = df$values ~ df$climate)

Residuals:
   Min     1Q Median     3Q    Max 
-32660 -10326  -2470   6394  94560 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        35632       8883   4.011 0.000314 ***
df$climateMediterranean climate   -21856      11290  -1.936 0.061240 .  
df$climateTemperate                 2988      10880   0.275 0.785291    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 25130 on 34 degrees of freedom
Multiple R-squared:  0.1841, Adjusted R-squared:  0.1361 
F-statistic: 3.836 on 2 and 34 DF,  p-value: 0.03146

test whether the mean GDP is the same across climate regions or varies significantly.

\(H_0 : \mu_\text{cold} = \mu_\text{temperate} = \mu_\text{meditteranean}\)

\(H_A : H_0 \text{ does not hold.}\)

aov(df$values ~ df$climate)
anova(aov(df$values ~ df$climate))

Call:
   aov(formula = df$values ~ df$climate)

Terms:
                 df$climate   Residuals
Sum of Squares   4843358867 21464355258
Deg. of Freedom           2          34

Residual standard error: 25125.77
Estimated effects may be unbalanced

at the 5% level we reject the null hypothesis, mean GDP varies by climate region.

dependencies between regressors

d <- subset(df, select = c("employment", "earnings", "internet", "climate"))

cor(d[c(1, 2, 3)])

options(repr.plot.width = 12, repr.plot.height = 8)
ggpairs(d, lower = list(combo = "box"))

the correlations between the numerical indicators are statistically significant. there’s also a significant dependence on the categorical indicator (climate zone).

aov(df$employment ~ df$climate)
anova(aov(df$employment ~ df$climate))

Call:
   aov(formula = df$employment ~ df$climate)

Terms:
                df$climate Residuals
Sum of Squares    939.3983  933.4817
Deg. of Freedom          2        34

Residual standard error: 5.239785
Estimated effects may be unbalanced

aov(df$earnings ~ df$climate)
anova(aov(df$earnings ~ df$climate))

Call:
   aov(formula = df$earnings ~ df$climate)

Terms:
                df$climate Residuals
Sum of Squares    17438893  77763451
Deg. of Freedom          2        34

Residual standard error: 1512.336
Estimated effects may be unbalanced

aov(df$internet ~ df$climate)
anova(aov(df$internet ~ df$climate))

Call:
   aov(formula = df$internet ~ df$climate)

Terms:
                df$climate Residuals
Sum of Squares    3195.714  3787.583
Deg. of Freedom          2        34

Residual standard error: 10.5546
Estimated effects may be unbalanced

statistical modeling

with linear regression (or a variant of it), we’ll model GDP against all the chosen regressors, interpret the coefficients, and check the model’s quality. we’ll spot outliers, look at multicollinearity, and test the model assumptions. if the assumptions fail we’ll try something that compensates or drops them. the goal is a final sub-model that explains GDP well without insignificant pieces.

dependence of GDP on all regressors

fit <- lm(df$values ~ df$employment + df$earnings + df$internet + df$climate)
summary(fit)

options(repr.plot.width = 10, repr.plot.height = 6)
plot(fit, which = 1)
plot(fit, which = 3)

Call:
lm(formula = df$values ~ df$employment + df$earnings + df$internet + 
    df$climate)

Residuals:
   Min     1Q Median     3Q    Max 
-22452  -3813   -697   2572  34164 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     -10160.722  25819.538  -0.394    0.697    
df$employment                      103.315    401.955   0.257    0.799    
df$earnings                         15.163      1.803   8.409 1.69e-09 ***
df$internet                        -19.648    273.609  -0.072    0.943    
df$climateMediterranean climate    917.870   7153.894   0.128    0.899    
df$climateTemperate               3433.730   5092.107   0.674    0.505    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10590 on 31 degrees of freedom
Multiple R-squared:  0.8679, Adjusted R-squared:  0.8466 
F-statistic: 40.72 on 5 and 31 DF,  p-value: 1.004e-12

everything but “earnings” is statistically insignificant at the 5% level, so they can probably be dropped individually. the model could be refined by dropping regressors gradually, clustering them more broadly, or stepping back to simpler ones.

interpretation of the coefficients

• for every 1% increase in the employment rate, GDP goes up by €103.
• for every €1 increase in average monthly income, GDP goes up by €15.
• for every 1% increase in internet usage, GDP goes down by €20.
• in Nordic countries, with all regressors at zero, GDP at market prices is estimated at €-10,160 per capita. relative to that baseline, GDP is €917 higher in Mediterranean countries and €3,433 higher in temperate climate countries.

the current model isn’t great for the reasons above. there’s a multicollinearity problem, the regressors are highly correlated with each other, which distorts the individual coefficients and their significance. the model also has oddities like the negative relationship between internet usage and GDP, which doesn’t match intuition. and it’s sensitive to outliers, especially at high GDP, which we look at next.

outliers

plot(fit, which = 4)

df[c(19, 21, 35), ]

the outliers are Liechtenstein, Luxembourg and Switzerland. relatively small countries with very strong financial sectors.

multicollinearity

conditionality of the matrix \(X^TX\)

Xnum <- df[, c("employment", "earnings", "internet")]
X <- as.matrix(Xnum)
XTX <- t(X) %*% X

XTX
cat("Singular values : ", round(svd(XTX)$d), "\n")
cat("Condition number: ", round(kappa(XTX, norm="2")))

Singular values :  274691304 85341 777 
Condition number:  497733

the matrix \(X^TX\) is very ill-conditioned (497733 \(\gg\) 1).

correlation

round(cor(Xnum), 3)

the regressors are strongly correlated.

variance inflation factor

data.frame(round(vif(lm(df$values ~ df$employment + df$earnings + df$internet)), 3))

VIF isn’t too large for any of the regressors.

ridge regression

to dig further into multicollinearity, apply ridge regression and compare it with standard regression without regularization.

fit_rr <- lm.ridge(df$values ~ df$employment + df$earnings + df$internet + df$climate, lambda = 10)
fit_lm <- lm(df$values ~ df$employment + df$earnings + df$internet + df$climate)

data.frame(fit_lm$coef, coef(fit_rr, scale = FALSE))

interpretation improves. internet usage gets a positive coefficient in the ridge regression model.

model assumptions

before using the diagnostic methods for submodel selection, a few assumptions need to hold.

regression model needs independent and equally distributed errors on top of a linear relationship between the explained and explanatory variables.

analysis of variance, like LR, needs independence and homogeneity of variances, and for the model selection method also normality of errors.

fit <- lm(df$values ~ df$employment + df$earnings + df$internet + df$climate)
residuals <- resid(fit)

options(repr.plot.width = 14, repr.plot.height = 5)
par(mfrow = c(1, 2))
plot(residuals)
hist(residuals, 10)
plot(fit, which = 1)
# plot(fit, which = 2)
plot(fit, which = 3)
# plot(fit, which = 4)

residue normality

the residual plot shows the residuals are distributed fairly symmetrically around zero.

test the normality of the residuals with a Q-Q plot and the Shapiro-Wilk test.

\(H_0:\) The residuals come from a normal distribution.
\(H_A:\) The residuals do not come from a normal distribution.

shapiro.test(residuals)
options(repr.plot.width = 10, repr.plot.height = 6)
plot(fit, which = 2)

 Shapiro-Wilk normality test

data:  residuals
W = 0.79839, p-value = 1.21e-05

at the 5% level we reject the normality of the residuals.

homogeneity of variance of residuals

test the homogeneity of the residual variance with the Breusch-Pagan test.

\(H_0:\) The residuals have identical variances.
\(H_A:\) The residuals do not have identical variances.

bptest(fit)

 studentized Breusch-Pagan test

data:  fit
BP = 18.577, df = 5, p-value = 0.002303

at the 5% level, the homogeneity of variances of the initial model is rejected.

correction

since the residuals grow with the explanatory variable, transform the GDP values. we take the square root and drop the outliers that stand out from the rest.

df2 <- df[-c(19, 21, 35),]
df2$values <- sqrt(df2$values)
fit2 <- lm(df2$values ~ df2$employment + df2$earnings + df2$internet + df2$climate)
summary(fit2)
residuals2 <- resid(fit2)

options(repr.plot.width = 14, repr.plot.height = 5)
par(mfrow = c(1, 2))
plot(residuals2)
hist(residuals2, 10)
plot(fit2, which = 1)
plot(fit2, which = 2)

Call:
lm(formula = df2$values ~ df2$employment + df2$earnings + df2$internet + 
    df2$climate)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.9576  -6.0405   0.1392   7.0079  16.2545 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      -4.54867   25.28852  -0.180   0.8585    
df2$employment                    0.97844    0.39462   2.479   0.0194 *  
df2$earnings                      0.03548    0.00229  15.494 2.89e-15 ***
df2$internet                      0.31952    0.28165   1.134   0.2662    
df2$climateMediterranean climate  3.10285    6.96468   0.446   0.6594    
df2$climateTemperate              3.46487    5.08587   0.681   0.5013    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.15 on 28 degrees of freedom
Multiple R-squared:  0.969, Adjusted R-squared:  0.9634 
F-statistic: 174.8 on 5 and 28 DF,  p-value: < 2.2e-16

residue normality

shapiro.test(residuals2)
options(repr.plot.width = 10, repr.plot.height = 6)
plot(fit2, which = 2)

 Shapiro-Wilk normality test

data:  residuals2
W = 0.95774, p-value = 0.2088

here, at the 5% level, we fail to reject the normality of residuals.

homogeneity of variance of residuals

bptest(fit2)

 studentized Breusch-Pagan test

data:  fit2
BP = 6.0235, df = 5, p-value = 0.3039

here too, at the 5% level, we fail to reject the homogeneity of residual variance.

choosing the final model

picking the final model now. it should be accurate, but we’d rather have one that’s simple and easy to interpret.

start with a model that has all the regressors, then simplify as we go. the final model should explain the variable well and have no insignificant pieces.

we’ll work with an all-flags model, transformed response, no outliers.

fit_overparam <- lm(df2$values ~ df2$employment + df2$earnings + df2$internet + df2$climate)
summary(fit_overparam)

Call:
lm(formula = df2$values ~ df2$employment + df2$earnings + df2$internet + 
    df2$climate)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.9576  -6.0405   0.1392   7.0079  16.2545 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      -4.54867   25.28852  -0.180   0.8585    
df2$employment                    0.97844    0.39462   2.479   0.0194 *  
df2$earnings                      0.03548    0.00229  15.494 2.89e-15 ***
df2$internet                      0.31952    0.28165   1.134   0.2662    
df2$climateMediterranean climate  3.10285    6.96468   0.446   0.6594    
df2$climateTemperate              3.46487    5.08587   0.681   0.5013    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.15 on 28 degrees of freedom
Multiple R-squared:  0.969, Adjusted R-squared:  0.9634 
F-statistic: 174.8 on 5 and 28 DF,  p-value: < 2.2e-16

the R-squared of 0.96 is very good, but a lot of the coefficients are insignificant, so we should try to simplify.

test model-submodel

test whether a richer model can be reduced to a simpler submodel. use a sequential procedure with the smallest and largest sets of parameters below.

• smallest: HDP' ~ 1 (constant)
• largest: HDP' ~ employment + earnings + internet + climate

step(lm(df2$values ~ 1), trace=T, scope = list(
    lower = ~1,
    upper = ~df2$values ~ df2$employment + df2$earnings + df2$internet + df2$climate
))

Start:  AIC=271.08
df2$values ~ 1

                 Df Sum of Sq   RSS    AIC
+ df2$earnings    1     88288  4714 171.69
+ df2$internet    1     64228 28774 233.19
+ df2$employment  1     35915 57087 256.48
+ df2$climate     2     24387 68615 264.74
<none>                        93002 271.08

Step:  AIC=171.69
df2$values ~ df2$earnings

                 Df Sum of Sq   RSS    AIC
+ df2$employment  1      1661  3053 158.92
+ df2$internet    1      1089  3625 164.75
+ df2$climate     2       730  3984 169.97
<none>                         4714 171.69
- df2$earnings    1     88288 93002 271.08

Step:  AIC=158.92
df2$values ~ df2$earnings + df2$employment

                 Df Sum of Sq   RSS    AIC
<none>                         3053 158.92
+ df2$internet    1       118  2935 159.58
+ df2$climate     2        33  3020 162.55
- df2$employment  1      1661  4714 171.69
- df2$earnings    1     54034 57087 256.48




Call:
lm(formula = df2$values ~ df2$earnings + df2$employment)

Coefficients:
   (Intercept)    df2$earnings  df2$employment  
       8.96734         0.03727         1.15468  

fit_reduced_1 <- lm(df2$values ~ df2$earnings + df2$employment)
summary(fit_reduced_1)

Call:
lm(formula = df2$values ~ df2$earnings + df2$employment)

Residuals:
     Min       1Q   Median       3Q      Max 
-18.6950  -7.1481   0.3365   9.1403  16.8266 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)     8.967342  14.205352   0.631 0.532495    
df2$earnings    0.037269   0.001591  23.423  < 2e-16 ***
df2$employment  1.154683   0.281167   4.107 0.000271 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.924 on 31 degrees of freedom
Multiple R-squared:  0.9672, Adjusted R-squared:  0.9651 
F-statistic: 456.7 on 2 and 31 DF,  p-value: < 2.2e-16

the sequential selection stopped at the lowest-AIC model, GDP' ~ earnings + employment, but the intercept is still insignificant at the 5% level.

# -intercept
fit_reduced_2 <- lm(df2$values ~ 0 + df2$earnings + df2$employment)
summary(fit_reduced_2)

Call:
lm(formula = df2$values ~ 0 + df2$earnings + df2$employment)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.1454  -7.1348   0.7599   8.8256  16.7243 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
df2$earnings   0.036914   0.001474   25.04   <2e-16 ***
df2$employment 1.327927   0.060566   21.93   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.83 on 32 degrees of freedom
Multiple R-squared:  0.9961, Adjusted R-squared:  0.9958 
F-statistic:  4068 on 2 and 32 DF,  p-value: < 2.2e-16

this simplification got rid of all the irrelevant components.

quality comparison, the R-squared values during this gradual simplification are interesting.

the last model has the highest adjusted R-squared.

test model-submodel

anova(fit_reduced_2,
      fit_reduced_1,
      fit_overparam,
      test="F")

more complex models can be reduced to the submodel GDP' ~ 0 + earnings + employment without a significant accuracy loss.

final model

m <- lm(df2$values ~ 0 + df2$earnings + df2$employment)
r <- resid(m)
summary(m)

options(repr.plot.width = 14, repr.plot.height = 6)
par(mfrow = c(1, 2))
plot(r)
plot(m, which = 3)

Call:
lm(formula = df2$values ~ 0 + df2$earnings + df2$employment)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.1454  -7.1348   0.7599   8.8256  16.7243 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
df2$earnings   0.036914   0.001474   25.04   <2e-16 ***
df2$employment 1.327927   0.060566   21.93   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.83 on 32 degrees of freedom
Multiple R-squared:  0.9961, Adjusted R-squared:  0.9958 
F-statistic:  4068 on 2 and 32 DF,  p-value: < 2.2e-16

with zero income and zero employment, GDP is zero. for every percent of employment, the square root of GDP increases by 1.33, and for every euro in average income, GDP increases by 0.04.

the final model is fine. it has a high R-squared, doesn’t make systematic errors, and is fairly simple. unlike the earlier models, it also gives valid values at low regressor levels.