Analysis of market price GDP per capita across European countries

I compared GDP per capita across European countries in 2014 and modeled it with a few country-level indicators.

# .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

The response is 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 

Only a few countries have very high GDP per capita. The median is about 20,000 EUR, and half of the countries are between roughly 11,000 and 39,000 EUR.

factors affecting GDP

GDP depends on many factors, for example:

economy and trade

• Inflation and interest rates. They affect consumer and investor confidence, and low, stable inflation supports growth.
• Exports and imports. Strong exports indicate a competitive economy, while high imports can indicate dependence on foreign trade.

labor market

• Employment. More people in work means more output and more purchasing power.
• Income levels. Job quality matters as well as employment: many low-wage jobs do not necessarily raise GDP per capita, and people with higher incomes spend more on goods and services.

other factors

• Life satisfaction, which is associated with productivity.

• Technology adoption. The ability to adopt new technology increases productivity and competitiveness.

• Climate. Floods, droughts, and extreme temperatures damage agriculture and infrastructure, while a favorable climate supports tourism and agriculture.

indicators

I picked four indicators as regressors:

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

First, each indicator and its relationship to GDP.

Each indicator is significantly correlated with GDP: at the 5% level, we reject zero correlation for all of them.

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 the last 3 months
• Time frequency: Annual
• Unit of measure: Percentage 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, I merged zones 2 and 3 into one.

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

Next, a test of whether mean GDP differs between climate regions.

\(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 differs between climate regions.

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

I started with a linear model on all four regressors, checked the coefficients, outliers, and model assumptions, and then simplified it. The goal was a model that fits well and is still easy to interpret.

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

All coefficients except “earnings” are insignificant at the 5% level. Because the regressors are correlated, this does not mean they can all be dropped at once, so I simplify the model step by step further below.

Interpretation of the coefficients:

• Each percentage point of employment rate adds €103 to GDP per capita.
• Each €1 of average monthly income adds €15.
• Each percentage point of internet usage subtracts €20.
• For Nordic countries with all regressors at zero, the estimated GDP is €-10,160 per capita. Relative to that baseline, GDP is €917 higher in Mediterranean countries and €3,433 higher in temperate countries.

This model has several problems. The regressors are correlated with each other, which inflates the variance of the individual coefficients and makes their significance unreliable. The negative coefficient for internet usage is implausible and is probably a consequence of this. The fit is also sensitive to outliers at high GDP, which I look at next.

outliers

plot(fit, which = 4)

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

The outliers are Liechtenstein, Luxembourg, and Switzerland, all small countries with large financial sectors.

multicollinearity

condition number 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))

The VIF is below 4 for all regressors, so the collinearity is moderate. The large condition number of \(X^TX\) is partly due to the very different scales of the regressors.

ridge regression

Next, ridge regression, to see how regularization affects the coefficients.

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))

With ridge regression, internet usage has a positive coefficient, which is more plausible than the negative one in the first model.

model assumptions

Before simplifying the model, I checked its assumptions.

Linear regression assumes a linear relationship between the response and the regressors, and independent, identically distributed errors.

ANOVA also assumes independence and equal variances, and the F-tests used for model selection additionally assume normally distributed 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)

normality of residuals

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

I checked 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

Breusch-Pagan test for constant variance:

\(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

The spread of the residuals grows with the fitted values. I took the square root of GDP and removed the outliers identified above.

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

normality of residuals

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

I started with all the regressors and removed terms that did not contribute, using the same transformed response and filtered data as above.

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

R-squared is 0.96, but many coefficients are insignificant, so the model can likely be simplified.

model-submodel test

I used stepwise selection to work between these two models:

• Smallest: GDP' ~ 1 (constant)
• Largest: GDP' ~ 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

Both remaining coefficients are significant.

R-squared during the simplification:

The last model has the highest adjusted R-squared, but the jump from 0.965 to 0.996 is mostly an artifact: without an intercept, R computes R-squared around zero instead of around the mean, so the values are not directly comparable. The F-test below is the fairer comparison.

Model-submodel test:

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

The model has no intercept, so GDP is zero when income and employment are zero. Each percentage point of employment adds 1.33 to the square root of GDP, and each euro of average monthly income adds 0.037.

The final model fits well, shows no systematic errors in the residuals, and has only two terms. Unlike the earlier models, it also gives non-negative predictions at low regressor values.