---
title: "Foreign Direct Investment and Renewable Energy Consumption in Sudan"
subtitle: "Econometrics Exercise 5 | Linear Regression - Basic Tests and Diagnostics"
author: "Ahmed Mustafa <br> (& Google Gemini)"
date: today
lang: en
format:
html:
toc: true
toc-location: left
toc-depth: 3
number-sections: true
theme: cosmo
embed-resources: true
code-fold: show
code-tools: true
execute:
echo: true
message: false
warning: false
---
```{css}
:root {
--brand: #AA272F;
--brand-dark: #7F1D24;
--brand-pale: #F7E9EA;
--ink: #292929;
}
h1, h2, h3 { color: var(--brand-dark); }
h1 { border-bottom: 2px solid var(--brand); padding-bottom: .25rem; }
.callout-note { border-left-color: var(--brand) !important; }
.callout-tip { border-left-color: #477A63 !important; }
.callout-warning { border-left-color: #B7791F !important; }
table caption { color: var(--brand-dark); font-weight: 700; }
table thead th { background: var(--brand-dark) !important; color: white !important; }
```
```{r}
#| label: setup
#| include: false
required_packages <- c(
"dplyr", "ggplot2", "broom", "knitr", "kableExtra",
"tseries", "car", "lmtest", "tidyr"
)
missing_packages <- required_packages[
!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages)) {
stop("Install these packages before rendering: ",
paste(missing_packages, collapse = ", "))
}
library(dplyr)
library(ggplot2)
library(broom)
library(knitr)
library(kableExtra)
brand <- "#AA272F"
brand_dark <- "#7F1D24"
brand_blue <- "#1B4F72"
brand_table <- function(x, caption, digits = 3) {
knitr::kable(
x, format = "html", digits = digits, caption = caption,
row.names = FALSE, escape = TRUE
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE, position = "left"
)
}
fmt_p <- function(p) {
ifelse(is.na(p), "not available",
ifelse(p < 0.001, "< 0.001", sprintf("%.3f", p)))
}
```
# Purpose and learning objectives
This assignment studies the relationship between foreign direct investment (FDI) and renewable energy consumption in Sudan. It uses OLS regression and focuses on the skills specified for Exercise 5:
1. Save and restore an R workspace.
2. Interpret coefficient t-tests, the overall F-test, and \(R^2\).
3. Test residual normality and identify potential outliers.
4. Inspect the four standard diagnostic plots produced by `plot(model)`.
5. Understand leverage and the hat matrix as supporting material.
::: {.callout-note}
## Scope and interpretation
The data contain annual observations. The analysis describes conditional associations in this sample; it does not, by itself, establish that FDI or another regressor causes changes in renewable energy consumption. The sample is small, so diagnostic results and p-values are interpreted cautiously.
:::
# 1. Data preparation
## 1.1 Import and prepare the data
The first chunk looks for a saved workspace. If it exists, `load()` restores the cleaned dataset. If not, the CSV file is imported and prepared.
```{r}
#| label: restore-or-create-workspace
workspace_file <- "sudan_exercise5_workspace.RData"
if (file.exists(workspace_file)) {
loaded_names <- load(workspace_file)
message("Workspace restored: ", paste(loaded_names, collapse = ", "))
} else {
csv_path <- if (file.exists("udaje/LDC-Panel-Data.csv")) {
"udaje/LDC-Panel-Data.csv"
} else if (file.exists("LDC-Panel-Data.csv")) {
"LDC-Panel-Data.csv"
} else {
stop("Cannot find LDC-Panel-Data.csv. Place it in the working directory or in 'udaje/'.")
}
raw_data <- read.csv(
csv_path, na.strings = c("", "NA"),
stringsAsFactors = FALSE
)
required_columns <- c(
"Country_Code", "Year", "REN_ERG_CONS", "FDI_GDP", "GDP_USD",
"IND_VA_GDP", "NODA_RCV", "Trade_GDP"
)
missing_columns <- setdiff(required_columns, names(raw_data))
if (length(missing_columns)) {
stop("Missing columns: ", paste(missing_columns, collapse = ", "))
}
analysis_vars <- c(
"REN_ERG_CONS", "FDI_GDP", "GDP_billions",
"IND_VA_GDP", "NODA_RCV", "Trade_GDP"
)
sudan_clean <- raw_data |>
filter(Country_Code == "SDN") |>
mutate(GDP_billions = GDP_USD / 1e9) |>
arrange(Year) |>
filter(if_all(all_of(analysis_vars), ~ !is.na(.x)))
if (nrow(sudan_clean) <= length(analysis_vars)) {
stop("Too few complete observations for the specified multiple regression.")
}
rownames(sudan_clean) <- as.character(sudan_clean$Year)
rm(raw_data, csv_path, required_columns, missing_columns)
}
```
::: {.callout-tip}
## If the workspace does not yet exist
Run the document once with `LDC-Panel-Data.csv` available. The cleaned dataset will be saved in Section 2.2. On a later run, the opening chunk will restore it with `load()`.
:::
## 2.2 Check the sample and save the cleaned workspace
```{r}
#| label: sample-check
sample_info <- data.frame(
Item = c("Country", "First year", "Last year", "Complete observations"),
Value = c(
"Sudan", min(sudan_clean$Year),
max(sudan_clean$Year), nrow(sudan_clean)
)
)
brand_table(sample_info, "Analysis sample", digits = 0)
save(sudan_clean, analysis_vars, file = workspace_file)
```
Only complete observations for the variables used in the analysis are retained. This is listwise deletion. The actual sample size and years are generated from the data rather than assumed in advance.
## 2.3 Practise removing objects with `rm()`
The following small example follows the classroom demonstration. It removes only the two temporary objects, not the data needed later in this document.
```{r}
#| label: rm-demonstration
temporary_number <- 6
temporary_text <- "temporary object"
rm(list = c("temporary_number", "temporary_text"))
ls()
```
To demonstrate a full workspace reset safely, save the objects first, clear the workspace, and restore them immediately:
```{r}
#| label: workspace-reset-demo
#| eval: false
save.image("sudan_full_workspace.RData")
rm(list = ls())
load("sudan_full_workspace.RData")
ls()
```
::: {.callout-warning}
Do not run `rm(list = ls())` casually in the middle of an analysis. It removes every object from the current environment. Use it only after saving the objects you need and when you are ready to restore them.
:::
# 3. Explore the data
## 3.1 Descriptive statistics
```{r}
#| label: descriptive-statistics
descriptive_stats <- data.frame(
Variable = analysis_vars,
N = vapply(sudan_clean[analysis_vars], length, integer(1)),
Minimum = vapply(sudan_clean[analysis_vars], min, numeric(1)),
Median = vapply(sudan_clean[analysis_vars], median, numeric(1)),
Mean = vapply(sudan_clean[analysis_vars], mean, numeric(1)),
Maximum = vapply(sudan_clean[analysis_vars], max, numeric(1)),
SD = vapply(sudan_clean[analysis_vars], sd, numeric(1))
)
brand_table(descriptive_stats, "Descriptive statistics for Sudan", digits = 2)
```
## 3.2 FDI and renewable energy consumption
```{r}
#| label: fdi-scatterplot
#| fig-cap: "FDI net inflows and renewable energy consumption in Sudan."
ggplot(sudan_clean, aes(x = FDI_GDP, y = REN_ERG_CONS)) +
geom_point(colour = brand_blue, size = 2.6, alpha = .85) +
geom_smooth(
method = "lm", se = TRUE,
colour = brand, fill = "#E8BCC0"
) +
labs(
x = "FDI net inflows (% of GDP)",
y = "Renewable energy consumption (% of total)",
title = "FDI and renewable energy consumption",
caption = "Annual observations; shaded area is the 95% confidence band."
) +
theme_minimal(base_size = 12)
```
::: {.callout-note}
## What the scatterplot can tell us
The plot gives a first impression of the bivariate association. It does not control for other variables, prove causality, or show whether the relationship is stable over time.
:::
# 4. Estimate the regression models
## 4.1 Simple and multiple OLS models
The simple model relates renewable energy consumption to FDI alone. The multiple model adds GDP, industry value added, official development assistance (ODA), and trade openness as controls.
$$
REN_t = \beta_0 + \beta_1 FDI_t + \beta_2 GDP_t + \beta_3 Industry_t + \beta_4 ODA_t + \beta_5 Trade_t + u_t
$$
```{r}
#| label: model-estimation
model_simple <- lm(REN_ERG_CONS ~ FDI_GDP, data = sudan_clean)
model_multiple <- lm(
REN_ERG_CONS ~ FDI_GDP + GDP_billions + IND_VA_GDP + NODA_RCV + Trade_GDP,
data = sudan_clean
)
coef_table <- tidy(model_multiple) |>
mutate(
Significance = case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.10 ~ ".",
TRUE ~ ""
)
)
coef_table |>
kable(
format = "html",
digits = 4,
col.names = c(
"Variable",
"Estimate (β̂)",
"Std. Error (SE)",
"t-statistic",
"p-value",
"Significance"
),
caption = "Multiple Linear Regression Results: Determinants of Renewable Energy Share in Sudan (1995–2021)"
) |>
kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)
```
## 4.2 Coefficient estimates and t-tests
For each coefficient, three diagnostic benchmarks exists:
$|t_j| > 1.96$,
$p_j < 0.05$,
$SE(\hat{\beta}_j) < \frac{1}{2}|\hat{\beta}_j|$.
```{r}
#| label: coefficient-table
coefficient_table <- broom::tidy(model_multiple, conf.int = TRUE) |>
mutate(
Decision_5pct = if_else(
p.value < 0.05, "Reject H0", "Do not reject H0"
)
) |>
select(
Term = term, Estimate = estimate, `Std. error` = std.error,
`t statistic` = statistic, `p-value` = p.value,
`95% CI lower` = conf.low, `95% CI upper` = conf.high,
Decision_5pct
)
brand_table(
coefficient_table,
"OLS coefficient estimates and individual t-tests",
digits = 4
)
```
::: {.callout-note}
## How to interpret a coefficient
A slope coefficient describes the estimated change in the conditional mean of the dependent variable associated with a one-unit increase in that explanatory variable, holding the other included regressors constant. A statistically insignificant coefficient is not proof that the true effect is exactly zero.
:::
## 4.3 Overall F-test and goodness of fit
The overall F-test evaluates whether all slope coefficients in the multiple regression are jointly zero:
\[
H_0:\beta_1=\beta_2=\cdots=\beta_5=0.
\]
The coefficient of determination \(R^2\) is the proportion of sample variation in the dependent variable accounted for by the fitted values. Adjusted \(R^2\) accounts for the number of regressors.
```{r}
#| label: model-fit
fit_table <- bind_rows(
broom::glance(model_simple) |>
mutate(Model = "Simple model"),
broom::glance(model_multiple) |>
mutate(Model = "Multiple model")
) |>
select(
Model, N = nobs, `R-squared` = r.squared,
`Adjusted R-squared` = adj.r.squared,
`Residual standard error` = sigma,
`F statistic` = statistic, `F-test p-value` = p.value
)
brand_table(fit_table, "Model fit and overall F-test", digits = 4)
```
::: {.callout-note}
## Interpreting \(R^2\) and the F-test
- A high \(R^2\) indicates that the fitted model accounts for a large share of the sample variation. It does not prove that the model is correctly specified or causal.
- A small overall F-test p-value indicates that the explanatory variables are jointly statistically significant.
- The F-test is a joint test; individual t-tests answer different questions.
- A higher \(R^2\) in a model with more regressors is not, by itself, proof that the larger model is better.
:::
# 5. Residual diagnostics
Residual diagnostics help assess whether the fitted model is compatible with important OLS assumptions. They do not provide a mechanical pass/fail certificate.
## 5.1 The four standard diagnostic plots
The command `plot(model_multiple)` produces four standard diagnostic plots.
```{r}
#| label: diagnostic-plots
#| fig-cap: "Standard diagnostic plots for the multiple regression model."
#| fig-width: 10
#| fig-height: 8
par(mfrow = c(2, 2))
plot(model_multiple)
par(mfrow = c(1, 1))
```
### 1. Residuals vs Fitted
Use this plot to look for non-linearity and unequal residual spread. Ideally, residuals are scattered around zero without a clear curve or funnel shape.
### 2. Normal Q–Q
This plot compares the residual quantiles with theoretical normal quantiles. Points close to the line are broadly consistent with normality. Systematic departures, especially in the tails, suggest possible non-normality.
### 3. Scale–Location
This plot shows the square root of the absolute standardized residuals against fitted values. A roughly horizontal band with similar spread is consistent with constant variance. A clear trend or funnel shape may indicate heteroscedasticity.
### 4. Residuals vs Leverage
This plot helps identify observations that combine unusual explanatory-variable values with large residuals. Cook's distance contours help flag potentially influential observations. High leverage alone does not imply high influence.
::: {.callout-note}
## How to write a diagnostic interpretation
Describe what is visible in the rendered plots. Do not insert numerical values or name particular years unless those values are produced by your current data and model. A diagnostic plot can suggest a problem, but it does not automatically identify its cause.
:::
## 5.2 Jarque–Bera test of residual normality
The Jarque–Bera test uses residual skewness and kurtosis to test normality.
\[
H_0:\text{the residual distribution is normal}
\qquad
H_1:\text{the residual distribution is not normal}.
\]
```{r}
#| label: jarque-bera
jb_test <- tseries::jarque.bera.test(residuals(model_multiple))
jb_table <- data.frame(
Statistic = unname(jb_test$statistic),
`p-value` = jb_test$p.value,
Test = jb_test$method
)
brand_table(jb_table, "Jarque–Bera test of residual normality", digits = 4)
```
::: {.callout-note}
## Decision rule
At the 5% level, reject the null hypothesis if the p-value is below 0.05.
If the p-value is at least 0.05, do not reject normality; this does not prove that the residuals are normally distributed.
The test is a supplement to the Q–Q plot, not a replacement. With a small sample, the test may have limited power.
:::
## 5.3 Outlier test
The function `car::outlierTest()` examines externally studentized residuals. It reports the most unusual observation and adjusts the p-value using the Bonferroni method to account for testing multiple observations.
```{r}
#| label: outlier-test
outlier_result <- car::outlierTest(model_multiple)
if (is.null(outlier_result)) {
cat("The function did not return an outlier candidate.")
} else {
outlier_table <- data.frame(
Observation = names(outlier_result$rstudent),
`Externally studentized residual` = unname(outlier_result$rstudent),
`Unadjusted p-value` = unname(outlier_result$p),
`Bonferroni adjusted p-value` = unname(outlier_result$bonf.p)
)
brand_table(outlier_table, "Bonferroni outlier test", digits = 4)
}
```
::: {.callout-note}
## Decision rule and caution
Use the Bonferroni-adjusted p-value for the 5% decision. A value below 0.05 flags an observation as statistically unusual relative to the fitted model.
An outlier is not necessarily a data error. Check the original record, the historical context, and the model specification before deciding what to do. Do not delete an observation simply because a test flags it.
:::
# 6. Leverage and the hat matrix
This section provides a short explanation of the concept introduced in the supplementary material. Leverage measures how unusual an observation's explanatory-variable values are relative to the rest of the sample.
For the OLS model,
\[
\widehat{\mathbf y}
= \mathbf X(\mathbf X^\mathsf{T}\mathbf X)^{-1}
\mathbf X^\mathsf{T}\mathbf y
= \mathbf H\mathbf y,
\]
where
\[
\mathbf H =
\mathbf X(\mathbf X^\mathsf{T}\mathbf X)^{-1}\mathbf X^\mathsf{T}
\]
is the hat matrix. Its diagonal elements \(h_{ii}\) are the leverage values. If the design matrix has full column rank and \(p\) columns, then:
\[
\sum_{i=1}^{n}h_{ii}=p,
\qquad
\overline h = \frac{p}{n}.
\]
Here, \(p\) includes the intercept when the model contains one.
```{r}
#| label: leverage-table
leverage_table <- data.frame(
Year = sudan_clean$Year,
Leverage = hatvalues(model_multiple),
Standardized_residual = rstandard(model_multiple),
Studentized_residual = rstudent(model_multiple),
Cooks_distance = cooks.distance(model_multiple)
) |>
arrange(desc(Cooks_distance))
brand_table(
head(leverage_table, 10),
"Ten observations with the largest Cook's distances",
digits = 4
)
p <- ncol(model.matrix(model_multiple))
n <- nobs(model_multiple)
cat("Number of model parameters, including intercept:", p, "\n")
cat("Average leverage p/n:", round(p / n, 4), "\n")
```
::: {.callout-note}
## Leverage is not the same as influence
Leverage depends on the explanatory variables, not on the observed dependent-variable value. An observation can have high leverage but a small residual. Cook's distance combines residual size and leverage to assess how much an observation may affect the fitted model.
Rules such as \(2p/n\), \(3p/n\), \(4/n\), or Cook's distance of 1 are only screening guidelines. They are not universal critical values and should not be used as automatic deletion rules.
:::
# 7. Compare and modify the model — Mini-tasks
We compare specifications estimated on the same observations and with the same dependent variable. This condition matters: AIC/BIC comparisons are not meaningful when models use different samples or different dependent-variable definitions.
## 7.1 Mini-task 1: compare models using AIC (and BIC)
We compare the simple model, the full multiple model, and a reduced model that omits the two variables whose individual significance is assessed in the full model. This reduced model is a teaching exercise, not an automatic variable-selection rule.
```{r}
#| label: reduced-model
model_reduced <- lm(
REN_ERG_CONS ~ FDI_GDP + GDP_billions + Trade_GDP,
data = sudan_clean
)
aic_bic_comparison <- bind_rows(
broom::glance(model_simple) |> mutate(Model = "Simple: FDI only"),
broom::glance(model_multiple) |> mutate(Model = "Full: all controls"),
broom::glance(model_reduced) |> mutate(Model = "Reduced: FDI, GDP, trade")
) |>
select(
Model, N = nobs,
`Log likelihood` = logLik, AIC, BIC,
`Adjusted R-squared` = adj.r.squared
) |>
arrange(AIC)
brand_table(aic_bic_comparison, "AIC and BIC model comparison", digits = 4)
```
```{r}
#| label: aic-bic-winner
aic_winner <- aic_bic_comparison$Model[which.min(aic_bic_comparison$AIC)]
bic_winner <- aic_bic_comparison$Model[which.min(aic_bic_comparison$BIC)]
cat("Lowest AIC: ", aic_winner, ".\n", sep = "")
cat("Lowest BIC: ", bic_winner, ".\n", sep = "")
```
::: {.callout-note}
**Interpretation:** Among the candidate models estimated on the same response and sample, lower AIC or BIC is preferred by that information criterion. AIC and BIC balance fit against model complexity; they do not prove that the selected model is true. BIC generally applies a stronger penalty for additional parameters.
:::
## 7.2 Mini-task 2: modify the model
As a transparent modification, estimate a log-GDP specification while keeping the other regressors. The logarithm is used only if GDP is strictly positive for every observation. The coefficient on log GDP is interpreted approximately as the change in the dependent variable associated with a 1% increase in GDP, holding other regressors constant.
```{r}
#| label: log-gdp-model
if (all(sudan_clean$GDP_billions > 0)) {
model_log_gdp <- lm(
REN_ERG_CONS ~ FDI_GDP + log(GDP_billions) +
IND_VA_GDP + NODA_RCV + Trade_GDP,
data = sudan_clean
)
log_model_available <- TRUE
} else {
model_log_gdp <- NULL
log_model_available <- FALSE
message("Log-GDP model not estimated: GDP contains zero or negative values.")
}
```
For a second, simpler modification, the reduced model above removes `IND_VA_GDP` and `NODA_RCV`. Compare the candidate models' fit and diagnostics rather than choosing a specification from p-values alone.
```{r}
#| label: modified-model-comparison
candidate_models <- list(
"Full model" = model_multiple,
"Reduced model" = model_reduced
)
if (log_model_available) {
candidate_models[["Log-GDP model"]] <- model_log_gdp
}
modified_fit <- bind_rows(lapply(names(candidate_models), function(nm) {
broom::glance(candidate_models[[nm]]) |>
mutate(Model = nm)
})) |>
select(
Model, N = nobs, `R-squared` = r.squared,
`Adjusted R-squared` = adj.r.squared,
`Residual SE` = sigma, `F-test p-value` = p.value,
AIC, BIC
)
brand_table(modified_fit, "Fit statistics for modified specifications", digits = 4)
```
## 7.3 Compare diagnostic results across specifications
```{r}
#| label: diagnostic-comparison
diagnostic_summary <- bind_rows(lapply(names(candidate_models), function(nm) {
m <- candidate_models[[nm]]
jb <- tseries::jarque.bera.test(residuals(m))
data.frame(
Model = nm,
JB_statistic = unname(jb$statistic),
JB_p_value = jb$p.value,
Maximum_abs_studentized_residual = max(abs(rstudent(m))),
Maximum_Cooks_distance = max(cooks.distance(m))
)
}))
brand_table(diagnostic_summary, "Residual diagnostic comparison", digits = 4)
```
```{r}
#| label: modified-diagnostic-plots
#| fig-width: 10
#| fig-height: 8
#| fig-cap: "Standard diagnostic plots for the selected reduced specification."
old_par <- par(no.readonly = TRUE)
par(
mfrow = c(2, 2),
col = brand_dark, col.axis = brand_dark, col.lab = brand_dark,
col.main = brand_dark, fg = brand_dark,
pch = 19, cex = 0.85
)
plot(model_reduced, col = brand, col.smooth = brand_dark)
par(old_par)
```
::: {.callout-note}
A model modification is not automatically an improvement. Consider the AIC/BIC, adjusted \(R^2\), residual standard error, normality test, diagnostic plots and economic interpretability together. Do not select a model solely because one statistic improves.
:::
# 8. Model review and Limitations
::: {.callout-warning}
## Important limitation: annual time-series data
The observations are ordered by year. The standard `lm()` summary and the diagnostic plots above do not test for serial correlation, structural breaks, or non-stationarity. These issues may matter in annual macroeconomic data. They are beyond the core requirements stated for this exercise, but should be considered before making strong economic or causal claims.
:::
# 9. Save the final analysis objects
Save the cleaned data, fitted models, and key diagnostic results so that they can be restored in a later exercise.
```{r}
#| label: save-final-workspace
save(
sudan_clean, analysis_vars,
model_simple, model_multiple,
jb_test, outlier_result, leverage_table,
file = "sudan_exercise5_results.RData"
)
```
To restore them in a future session:
```{r}
#| eval: false
load("sudan_exercise5_results.RData")
ls()
```
# Final checklist
- [ ] The data file and sample selection are documented.
- [ ] The workspace is saved and can be restored.
- [ ] Coefficient t-tests are interpreted using p-values or the correct t critical value.
- [ ] The overall F-test is interpreted as a joint test.
- [ ] \(R^2\) and adjusted \(R^2\) are explained correctly.
- [ ] All four `plot(model)` diagnostic panels are discussed.
- [ ] The Jarque–Bera result is interpreted without claiming that normality has been proved.
- [ ] The outlier test is interpreted using its Bonferroni-adjusted p-value.
- [ ] High leverage is distinguished from influence.
- [ ] No observation is deleted without substantive investigation.
- [ ] Malé úlohy
