# List of packages
packages <- c("tidyverse", "fst", "modelsummary", "broom", "sjPlot", "ggplot2", "car", "Lock5Data", "pandoc", "mosaic") # add any you need here

# Install packages if they aren't installed already
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

# Load the packages
lapply(packages, library, character.only = TRUE)
## Warning: package 'ggplot2' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Warning: package 'fst' was built under R version 4.3.3
## Warning: package 'modelsummary' was built under R version 4.3.3
## `modelsummary` 2.0.0 now uses `tinytable` as its default table-drawing
##   backend. Learn more at: https://vincentarelbundock.github.io/tinytable/
## 
## Revert to `kableExtra` for one session:
## 
##   options(modelsummary_factory_default = 'kableExtra')
##   options(modelsummary_factory_latex = 'kableExtra')
##   options(modelsummary_factory_html = 'kableExtra')
## 
## Silence this message forever:
## 
##   config_modelsummary(startup_message = FALSE)
## Warning: package 'sjPlot' was built under R version 4.3.3
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
## Warning: package 'pandoc' was built under R version 4.3.3
## Warning: package 'mosaic' was built under R version 4.3.3
## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## 
## The following object is masked from 'package:Matrix':
## 
##     mean
## 
## The following objects are masked from 'package:car':
## 
##     deltaMethod, logit
## 
## The following object is masked from 'package:modelsummary':
## 
##     msummary
## 
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     stat
## 
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## 
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum
## [[1]]
##  [1] "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
##  [7] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics" 
## [13] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[2]]
##  [1] "fst"       "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"    
##  [7] "readr"     "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"    
## [13] "graphics"  "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[3]]
##  [1] "modelsummary" "fst"          "lubridate"    "forcats"      "stringr"     
##  [6] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [11] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [16] "utils"        "datasets"     "methods"      "base"        
## 
## [[4]]
##  [1] "broom"        "modelsummary" "fst"          "lubridate"    "forcats"     
##  [6] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [11] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [16] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## 
## [[5]]
##  [1] "sjPlot"       "broom"        "modelsummary" "fst"          "lubridate"   
##  [6] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [11] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [16] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [21] "base"        
## 
## [[6]]
##  [1] "sjPlot"       "broom"        "modelsummary" "fst"          "lubridate"   
##  [6] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [11] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [16] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [21] "base"        
## 
## [[7]]
##  [1] "car"          "carData"      "sjPlot"       "broom"        "modelsummary"
##  [6] "fst"          "lubridate"    "forcats"      "stringr"      "dplyr"       
## [11] "purrr"        "readr"        "tidyr"        "tibble"       "ggplot2"     
## [16] "tidyverse"    "stats"        "graphics"     "grDevices"    "utils"       
## [21] "datasets"     "methods"      "base"        
## 
## [[8]]
##  [1] "Lock5Data"    "car"          "carData"      "sjPlot"       "broom"       
##  [6] "modelsummary" "fst"          "lubridate"    "forcats"      "stringr"     
## [11] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [16] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [21] "utils"        "datasets"     "methods"      "base"        
## 
## [[9]]
##  [1] "pandoc"       "Lock5Data"    "car"          "carData"      "sjPlot"      
##  [6] "broom"        "modelsummary" "fst"          "lubridate"    "forcats"     
## [11] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [16] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [21] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## 
## [[10]]
##  [1] "mosaic"       "mosaicData"   "ggformula"    "Matrix"       "lattice"     
##  [6] "pandoc"       "Lock5Data"    "car"          "carData"      "sjPlot"      
## [11] "broom"        "modelsummary" "fst"          "lubridate"    "forcats"     
## [16] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [21] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [26] "grDevices"    "utils"        "datasets"     "methods"      "base"

Task 1

Calculate and interpret the correlation between SAT math scores (sat_math) and freshman GPA (gpa_fy). Next, visualize this relationship with a scatterplot, including a regression line.

# Load the data
sat_gpa <- read_csv("sat_gpa.csv")
## Rows: 1000 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): sex
## dbl (5): sat_verbal, sat_math, sat_total, gpa_hs, gpa_fy
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Calculate the correlation between SAT math scores and freshman GPA
correlation <- cor(sat_gpa$sat_math, sat_gpa$gpa_fy)
print(paste("Correlation between SAT math scores and freshman GPA:", correlation))
## [1] "Correlation between SAT math scores and freshman GPA: 0.387117762067329"

#Interpretation The correlation coefficient is positive, indicating that as SAT math scores increase, freshman GPA tends to increase as well. This suggests a positive linear relationship between the two variables. A correlation coefficient of 0.387 indicates a moderate relationship. While there is a noticeable association between SAT math scores and freshman GPA, it is not very strong.

# Create a scatterplot with a regression line
ggplot(sat_gpa, aes(x = sat_math, y = gpa_fy)) +
  geom_point(size = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = "SAT Math Score", y = "Freshman GPA") +
  ggtitle(paste("Scatterplot of SAT Math Scores vs. Freshman GPA\nCorrelation: ", round(correlation, 2)))
## `geom_smooth()` using formula = 'y ~ x'

Task 2

Create regression tables to analyze the relationship between SAT scores (total and verbal) and freshman GPA, using both the modelsummary and sjPlot packages. Interpret.

# Fit the regression models
model_sat_total <- lm(gpa_fy ~ sat_total, data = sat_gpa)
model_sat_verbal <- lm(gpa_fy ~ sat_verbal, data = sat_gpa)

# Display regression tables using modelsummary
modelsummary(list(
  "SAT Total Score" = model_sat_total,
  "SAT Verbal Score" = model_sat_verbal
))
tinytable_a9z1w55v0i5p8nu032ev
SAT Total Score SAT Verbal Score
(Intercept) 0.002 0.701
(0.152) (0.129)
sat_total 0.024
(0.001)
sat_verbal 0.036
(0.003)
Num.Obs. 1000 1000
R2 0.212 0.161
R2 Adj. 0.211 0.160
AIC 2004.8 2067.2
BIC 2019.5 2081.9
Log.Lik. -999.382 -1030.580
F 268.270 191.673
RMSE 0.66 0.68 
# Display regression tables using sjPlot
tab_model(model_sat_total, model_sat_verbal, 
          title = "Regression Models Predicting Freshman GPA", 
          dv.labels = c("SAT Total Score", "SAT Verbal Score"))
Regression Models Predicting Freshman GPA
  SAT Total Score SAT Verbal Score
Predictors Estimates CI p Estimates CI p
(Intercept) 0.00 -0.30 – 0.30 0.990 0.70 0.45 – 0.95 <0.001
sat total 0.02 0.02 – 0.03 <0.001
sat verbal 0.04 0.03 – 0.04 <0.001
Observations 1000 1000
R2 / R2 adjusted 0.212 / 0.211 0.161 / 0.160

#Interpretation

The regression analysis reveals that both SAT total and verbal scores are significant predictors of freshman GPA. Specifically, for every one-point increase in SAT total score, the freshman GPA increases by an average of 0.02 points (p < 0.001), explaining 21.2% of the variance in GPA. Similarly, a one-point increase in SAT verbal score corresponds to an average GPA increase of 0.04 points (p < 0.001), accounting for 16.1% of the variance. The intercepts indicate the baseline GPA when SAT scores are zero, with the total score model intercept not significant (p = 0.990) and the verbal score model intercept significant (p < 0.001). These findings suggest that while both SAT total and verbal scores positively influence freshman GPA, the SAT total score has a slightly stronger explanatory power.

Task 3

Explore the relationship between Happy Life Years (HLY) and GDP per Capita using the HappyPlanetIndex dataset. Visualize the coefficients using both the modelsummary and Sjplot packages.

data("HappyPlanetIndex")
World <- HappyPlanetIndex
names(HappyPlanetIndex)
##  [1] "Country"        "Region"         "Happiness"      "LifeExpectancy"
##  [5] "Footprint"      "HLY"            "HPI"            "HPIRank"       
##  [9] "GDPperCapita"   "HDI"            "Population"
# Fit the regression model
model_hly_gdp <- lm(HLY ~ GDPperCapita, data = HappyPlanetIndex)

# Display regression tables using modelsummary
modelsummary(model_hly_gdp)
tinytable_z6amdvenm4z225h3gyfu
(1)
(Intercept) 31.182
(1.114)
GDPperCapita 0.001
(0.000)
Num.Obs. 141
R2 0.566
R2 Adj. 0.563
AIC 1043.2
BIC 1052.0
Log.Lik. -518.576
F 181.014
RMSE 9.57 
# Display regression tables using sjPlot
tab_model(model_hly_gdp, 
          title = "Regression Model Predicting Happy Life Years (HLY)", 
          dv.labels = c("GDP per Capita"))
Regression Model Predicting Happy Life Years (HLY)
  GDP per Capita
Predictors Estimates CI p
(Intercept) 31.18 28.98 – 33.38 <0.001
GDPperCapita 0.00 0.00 – 0.00 <0.001
Observations 141
R2 / R2 adjusted 0.566 / 0.563 
# Visualize the coefficients using ggplot2
coef_data <- tidy(model_hly_gdp, conf.int = TRUE)

ggplot(coef_data, aes(x = term, y = estimate, color = term)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2) +
  labs(title = "Coefficients of the Regression Model",
       x = "Predictor",
       y = "Estimate") +
  theme_minimal() +
  scale_color_manual(values = c("(Intercept)" = "blue", "GDP_per_Capita" = "red")) +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5)) +
  guides(color = guide_legend(override.aes = list(size = 5)))

Task 4

Visualize the coefficients of the models specified below using both modelplot from the modelsummary package and plot_model from the sjPlot package. Fit three models to predict happiness (Happiness) using different predictors:

Model 1: Life Expectancy (LifeExpectancy)

Model 2: Ecological Footprint (Footprint)

Model 3: GDP per Capita (GDPperCapita), Human Development Index (HDI), and Population (Population)

# Fit the regression models
m4 <- lm(Happiness ~ LifeExpectancy, data = World)
m5 <- lm(Happiness ~ Footprint, data = World)
m6 <- lm(Happiness ~ GDPperCapita + HDI + Population, data = World)

# Create a list of models
models_world <- list(m4, m5, m6)

# Generate regression tables using modelsummary
modelsummary(models_world)
tinytable_7j3uumia9ciq2gdlugic
(1) (2) (3)
(Intercept) -1.104 4.713 1.528
(0.399) (0.150) (0.357)
LifeExpectancy 0.104
(0.006)
Footprint 0.419
(0.042)
GDPperCapita 0.000
(0.000)
HDI 5.778
(0.577)
Population 0.001
(0.000)
Num.Obs. 143 143 141
R2 0.693 0.414 0.707
R2 Adj. 0.691 0.410 0.700
AIC 332.7 425.0 327.6
BIC 341.6 433.9 342.3
Log.Lik. -163.359 -209.524 -158.779
F 318.074 99.701 109.954
RMSE 0.76 1.05 0.75 
# Generate regression tables using sjPlot
tab_model(m4, m5, m6, 
          title = "Regression Models Predicting Happiness", 
          dv.labels = c("Model 1: Life Expectancy", "Model 2: Ecological Footprint", "Model 3: Multiple Predictors"))
Regression Models Predicting Happiness
  Model 1: Life Expectancy Model 2: Ecological Footprint Model 3: Multiple Predictors
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) -1.10 -1.89 – -0.32 0.006 4.71 4.42 – 5.01 <0.001 1.53 0.82 – 2.23 <0.001
LifeExpectancy 0.10 0.09 – 0.11 <0.001
Footprint 0.42 0.34 – 0.50 <0.001
GDPperCapita 0.00 -0.00 – 0.00 0.100
HDI 5.78 4.64 – 6.92 <0.001
Population 0.00 -0.00 – 0.00 0.247
Observations 143 143 141
R2 / R2 adjusted 0.693 / 0.691 0.414 / 0.410 0.707 / 0.700 
# Display the models summary table with improved presentation 
modelsummary(models_world,
             statistic = "conf.int",
             coef_map = c("LifeExpectancy" = "Life Expectancy",
                          "Footprint" = "Ecological Footprint",
                          "GDPperCapita" = "GDP per Capita",
                          "HDI" = "Human Development Index",
                          "Population" = "Population (millions)"),
             gof_map = c("nobs", "r.squared", "adj.r.squared", "aic", "bic"))
tinytable_hmb3ahcuv3iv8yocfnyw
(1) (2) (3)
Life Expectancy 0.104
[0.092, 0.115]
Ecological Footprint 0.419
[0.336, 0.502]
GDP per Capita 0.000
[0.000, 0.000]
Human Development Index 5.778
[4.636, 6.920]
Population (millions) 0.001
[0.000, 0.001]
Num.Obs. 143 143 141
R2 0.693 0.414 0.707
R2 Adj. 0.691 0.410 0.700
AIC 332.7 425.0 327.6
BIC 341.6 433.9 342.3 
# Visualize 
modelplot(models_world, coef_omit = "Intercept") +
    labs(x = 'Coefficient Estimate',
         y = 'Term',
         title = 'Model Coefficients with Confidence Intervals',
         caption = 'Comparison of Models 4, 5, and 6') +
    theme_minimal()

# Visualize Model 4: Life Expectancy
plot_model(m4, show.values = TRUE, value.offset = .3, ci.lvl = 0.95) +
    labs(title = 'Model Coefficients for Life Expectancy',
         caption = 'Model 1: Life Expectancy') +
    theme_minimal()

# Visualize Model 5: Ecological Footprint
plot_model(m5, show.values = TRUE, value.offset = .3, ci.lvl = 0.95) +
    labs(title = 'Model Coefficients for Ecological Footprint',
         caption = 'Model 2: Ecological Footprint') +
    theme_minimal()

# Visualize Model 6: Multiple Predictors
plot_model(m6, show.values = TRUE, value.offset = .3, ci.lvl = 0.95) +
    labs(title = 'Model Coefficients for Multiple Predictors',
         caption = 'Model 3: Multiple Predictors') +
    theme_minimal()

Task 5

Predict murder rates using both internet penetration and GDP, from the Violence dataset. Create regression tables for the model using both the modelsummary and sjPlot packages. Interpret.

# Load the Violence dataset
load("Violence.Rdata")
names(Violence)
##  [1] "Country"        "Code"           "LandArea"       "Population"    
##  [5] "Energy"         "Rural"          "Military"       "Health"        
##  [9] "HIV"            "Internet"       "Developed"      "BirthRate"     
## [13] "ElderlyPop"     "LifeExpectancy" "CO2"            "GDP"           
## [17] "Cell"           "Electricity"    "MurderRate"
# Check the number of observational units and variables
nrow(Violence)  # Number of observational units
## [1] 8
ncol(Violence)  # Number of variables
## [1] 19
# Check variable names
names(Violence)
##  [1] "Country"        "Code"           "LandArea"       "Population"    
##  [5] "Energy"         "Rural"          "Military"       "Health"        
##  [9] "HIV"            "Internet"       "Developed"      "BirthRate"     
## [13] "ElderlyPop"     "LifeExpectancy" "CO2"            "GDP"           
## [17] "Cell"           "Electricity"    "MurderRate"
# Calculate the average murder rate per 100,000 people
mean_murder_rate <- mean(Violence$MurderRate, na.rm = TRUE)
mean_murder_rate
## [1] 22.55
# Create a scatter between the percentage of people living in a rural area (Rural) and the murder rate per 100,000 people (MurderRate)
plot(Violence$Rural, Violence$MurderRate, 
     xlab = "Percentage in Rural Area", 
     ylab = "Murder Rate per 100,000", 
     pch = 16, 
     main = "Scatterplot of Rural Area Percentage vs Murder Rate")

# Fit a simple linear regression predicting the murder rate using the percentage living in a rural area
model_rural <- lm(MurderRate ~ Rural, data = Violence)
summary(model_rural)
## 
## Call:
## lm(formula = MurderRate ~ Rural, data = Violence)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.172  -5.687   2.190   7.287  19.577 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -13.3276    12.9203  -1.032   0.3421  
## Rural         1.1296     0.3697   3.055   0.0224 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.24 on 6 degrees of freedom
## Multiple R-squared:  0.6087, Adjusted R-squared:  0.5435 
## F-statistic: 9.334 on 1 and 6 DF,  p-value: 0.02236
# Fit a regression model predicting the murder rate using internet penetration and GDP
model_internet_gdp <- lm(MurderRate ~ Internet + GDP, data = Violence)
summary(model_internet_gdp)
## 
## Call:
## lm(formula = MurderRate ~ Internet + GDP, data = Violence)
## 
## Residuals:
##            Austria            Belgium          Guatemala            Jamaica 
##             -2.507             -2.822             14.174             10.453 
## Dominican Republic       South Africa            Ukraine      United States 
##             -6.930             13.338            -29.231              3.526 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 28.9842977 11.9299808   2.430   0.0594 .
## Internet     0.4629043  0.4385000   1.056   0.3394  
## GDP         -0.0013199  0.0006033  -2.188   0.0803 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.84 on 5 degrees of freedom
## Multiple R-squared:  0.602,  Adjusted R-squared:  0.4429 
## F-statistic: 3.782 on 2 and 5 DF,  p-value: 0.09991
# Create a list of models
models_violence <- list("Rural Area Percentage" = model_rural, "Internet Penetration and GDP" = model_internet_gdp)

# Generate regression tables using modelsummary
modelsummary(models_violence)
tinytable_0es3mj2njq6zuqe7bwxr
Rural Area Percentage Internet Penetration and GDP
(Intercept) -13.328 28.984
(12.920) (11.930)
Rural 1.130
(0.370)
Internet 0.463
(0.438)
GDP -0.001
(0.001)
Num.Obs. 8 8
R2 0.609 0.602
R2 Adj. 0.543 0.443
AIC 70.0 72.1
BIC 70.2 72.4
Log.Lik. -31.992 -32.059
F 9.334 3.782
RMSE 13.20 13.31 
# Generate regression tables using sjPlot
tab_model(model_rural, model_internet_gdp, 
          title = "Regression Models Predicting Murder Rate", 
          dv.labels = c("Model 1: Rural Area Percentage", "Model 2: Internet Penetration and GDP"))
Regression Models Predicting Murder Rate
  Model 1: Rural Area Percentage Model 2: Internet Penetration and GDP
Predictors Estimates CI p Estimates CI p
(Intercept) -13.33 -44.94 – 18.29 0.342 28.98 -1.68 – 59.65 0.059
Rural 1.13 0.22 – 2.03 0.022
Internet 0.46 -0.66 – 1.59 0.339
GDP -0.00 -0.00 – 0.00 0.080
Observations 8 8
R2 / R2 adjusted 0.609 / 0.543 0.602 / 0.443 
        # Display the models summary table with improved presentation using modelsummary
modelsummary(models_violence,
             statistic = "conf.int",
             coef_map = c("Rural" = "Rural Area Percentage",
                          "Internet" = "Internet Penetration",
                          "GDP" = "GDP per Capita"),
             gof_map = c("nobs", "r.squared", "adj.r.squared", "aic", "bic"))
tinytable_009slvgf6r69gvzryyz8
Rural Area Percentage Internet Penetration and GDP
Rural Area Percentage 1.130
[0.225, 2.034]
Internet Penetration 0.463
[-0.664, 1.590]
GDP per Capita -0.001
[-0.003, 0.000]
Num.Obs. 8 8
R2 0.609 0.602
R2 Adj. 0.543 0.443
AIC 70.0 72.1
BIC 70.2 72.4 
# Visualize with modelplot from modelsummary
modelplot(models_violence, coef_omit = "Intercept") +
    labs(x = 'Coefficient Estimate',
         y = 'Term',
         title = 'Model Coefficients with Confidence Intervals',
         caption = 'Comparison of Models 1 and 2') +
    theme_minimal()

# Visualize Model 1: Rural Area Percentage
plot_model(model_rural, show.values = TRUE, value.offset = .3, ci.lvl = 0.95) +
    labs(title = 'Model Coefficients for Rural Area Percentage',
         caption = 'Model 1: Rural Area Percentage') +
    theme_minimal()

# Visualize Model 2: Internet Penetration and GDP
plot_model(model_internet_gdp, show.values = TRUE, value.offset = .3, ci.lvl = 0.95) +
    labs(title = 'Model Coefficients for Internet Penetration and GDP',
         caption = 'Model 2: Internet Penetration and GDP') +
    theme_minimal()

Happy Coding