Stat 1221 Research

# Import Dependencies
library(readxl)
require(tidyverse)
library(car)
library(psych)

# Read in data
gini <- read.csv("gini-coefficient-wid-vs-pip.csv", header = TRUE)
gdp_raw <- read_excel("API_NY.GDP.PCAP.PP.KD_DS2_en_excel_v2_269.xls", skip = 3)
tourism_intensity <- read.csv("international-tourist-trips-per-1000-people.csv", header = TRUE)
urban_rural_pop <- read.csv("urban-and-rural-population.csv", header = TRUE)
poverty_rate <- read.csv("worldpoverty.csv")

# poverty rate cleaning

poverty_clean <- poverty_rate %>%
  filter(survey_coverage == "national") %>%
  filter(headcount > 0) %>%
  select(country_code,reporting_year, headcount) %>%
  rename(code = country_code, year = reporting_year) %>%
  group_by(code) %>%
    slice_max(year, n = 1, with_ties = FALSE) %>%
    select(code, headcount)
    

# tourism intensity cleaning

tourism_clean <- tourism_intensity %>%
  select(Code,Year,Arrivals.of.tourists.from.abroad.per.1.000.people) %>%
  rename(code = Code, year = Year, tourism_intensity = Arrivals.of.tourists.from.abroad.per.1.000.people) %>%
  group_by(code) %>%
    slice_max(year, n = 1,with_ties = FALSE )
    

# Urbanization Rate Cleaning
urbanization_clean <- urban_rural_pop %>%
  mutate(population = Urban.population + Rural.population) %>%
  mutate(urbanization_rate = Urban.population/(Urban.population + Rural.population)) %>%
  select(Code,Year,urbanization_rate, population)%>%
  rename(code = Code, year = Year) %>%
  group_by(code) %>%
    slice_max(year, n = 1, with_ties = FALSE) %>%
    select(code, urbanization_rate, population)

# GDP Cleaning
gdp_clean <- gdp_raw %>%
  pivot_longer(cols = -c(`Country Name`, `Country Code`, `Indicator Name`, `Indicator Code`),
             names_to = "year",
             values_to = "gdp_per_capita_ppp") %>%
  drop_na(gdp_per_capita_ppp) %>%
  select(`Country Code`,year,gdp_per_capita_ppp)%>%
  rename(code = `Country Code`) %>%
  mutate(year = as.integer(year)) %>%
  group_by(code) %>%
    slice_max(year, n = 1, with_ties = FALSE) %>%
    select(code, gdp_per_capita_ppp)

# Clean Gini
gini_clean <- gini %>%
  select(Code,Year,Gini.according.to.the.World.Inequality.Database)%>%
  rename(code = Code, year = Year, gini = Gini.according.to.the.World.Inequality.Database) %>%
  drop_na(gini) %>%
    filter(code != "") %>%
    group_by(code) %>%
      slice_max(year, n = 1, with_ties = FALSE) %>%
        select(code, gini)

# Merge the Data
df <- tourism_clean %>%
  inner_join(gdp_clean, by = "code") %>%
  inner_join(gini_clean, by = "code") %>%
  inner_join(urbanization_clean, by = "code") %>%
  inner_join(poverty_clean, by = "code")

nrow(df)

[1] 156

# Work with Cleaned Data
description <- describe(df[, c("tourism_intensity", "gdp_per_capita_ppp", "headcount", "gini", "urbanization_rate", "population")])

description

# Read in Dependencies
require(ggplot2)

# Poverty Rate vs Tourism Intensity
pr_vs_ti_plot <- ggplot(data = df, mapping = aes(x = log(tourism_intensity), y=log(headcount+0.001))) + geom_point() + geom_smooth(method = "lm") + theme_classic() + labs(x = "log(Tourism Intensity)", y = "log(Poverty Rate)")

pr_vs_ti_plot

# GDP per capita vs tourism intensity
gdp_vs_ti_plot <- ggplot(data = df, mapping = aes(x = log(tourism_intensity), y = log(gdp_per_capita_ppp))) + geom_point() + geom_smooth(method = "lm") + theme_classic() + labs(x = "log(Tourism Intensity)", y = "log(GDP Per Capita(PPP)")

gdp_vs_ti_plot

# popluation vs tourism intensity

pop_vs_ti <- ggplot(data = df, mapping = aes(x = log(tourism_intensity), y = log(population))) + geom_point() + geom_smooth(method = "lm") + theme_classic() + labs(x = "log(Tourism Intensity)", y = "log(Population")

pop_vs_ti

# Fit the models
model1 <- lm(log(gdp_per_capita_ppp) ~ log(tourism_intensity) + gini + urbanization_rate + log(population), data = df)

model2 <- lm(log(headcount+ 0.001) ~ log(tourism_intensity) + gini + urbanization_rate + log(population), data = df)

# Model1 Work
summary(model1)

Call:
lm(formula = log(gdp_per_capita_ppp) ~ log(tourism_intensity) + 
    gini + urbanization_rate + log(population), data = df)

Residuals:
     Min       1Q   Median 
-1.54177 -0.35476  0.00422 
      3Q      Max 
 0.35666  2.38459 

Coefficients:
                       Estimate
(Intercept)             7.48426
log(tourism_intensity)  0.26262
gini                   -2.80598
urbanization_rate       1.65190
log(population)         0.07523
                       Std. Error
(Intercept)               0.56849
log(tourism_intensity)    0.02933
gini                      0.51045
urbanization_rate         0.27829
log(population)           0.02632
                       t value
(Intercept)             13.165
log(tourism_intensity)   8.954
gini                    -5.497
urbanization_rate        5.936
log(population)          2.858
                       Pr(>|t|)
(Intercept)             < 2e-16
log(tourism_intensity) 1.18e-15
gini                   1.61e-07
urbanization_rate      1.93e-08
log(population)         0.00487
                          
(Intercept)            ***
log(tourism_intensity) ***
gini                   ***
urbanization_rate      ***
log(population)        ** 
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6134 on 151 degrees of freedom
Multiple R-squared:  0.7152,    Adjusted R-squared:  0.7077 
F-statistic: 94.82 on 4 and 151 DF,  p-value: < 2.2e-16

plot(model1)

vif(model1)

log(tourism_intensity) 
              1.695790 
                  gini 
              1.245817 
     urbanization_rate 
              1.431155 
       log(population) 
              1.206550 

# Model2 Work
summary(model2)

Call:
lm(formula = log(headcount + 0.001) ~ log(tourism_intensity) + 
    gini + urbanization_rate + log(population), data = df)

Residuals:
    Min      1Q  Median      3Q 
-4.7648 -0.7944  0.1379  0.9979 
    Max 
 2.4735 

Coefficients:
                       Estimate
(Intercept)            -3.60577
log(tourism_intensity) -0.34408
gini                    6.59317
urbanization_rate      -2.41382
log(population)        -0.00943
                       Std. Error
(Intercept)               1.21729
log(tourism_intensity)    0.06280
gini                      1.09301
urbanization_rate         0.59589
log(population)           0.05637
                       t value
(Intercept)             -2.962
log(tourism_intensity)  -5.479
gini                     6.032
urbanization_rate       -4.051
log(population)         -0.167
                       Pr(>|t|)
(Intercept)             0.00355
log(tourism_intensity) 1.75e-07
gini                   1.20e-08
urbanization_rate      8.14e-05
log(population)         0.86735
                          
(Intercept)            ** 
log(tourism_intensity) ***
gini                   ***
urbanization_rate      ***
log(population)           
---
Signif. codes:  
  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’
  0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.313 on 151 degrees of freedom
Multiple R-squared:  0.5854,    Adjusted R-squared:  0.5744 
F-statistic:  53.3 on 4 and 151 DF,  p-value: < 2.2e-16

plot(model2)

vif(model2)

log(tourism_intensity) 
              1.695790 
                  gini 
              1.245817 
     urbanization_rate 
              1.431155 
       log(population) 
              1.206550 

# Univariate Descriptive Statistics
hist1 <- ggplot(data = df, mapping = aes(x = gdp_per_capita_ppp)) + geom_histogram() + labs(x = "GDP Per Capita(PPP)", y = "Count") + theme_classic()

hist1

fig2 <- ggplot(data = df, mapping = aes(x = headcount)) + geom_histogram() + labs(x = "Poverty Rate", y = "Count") + theme_classic()

fig2

fig3 <- ggplot(data = df, mapping = aes(x = tourism_intensity)) + geom_histogram() + labs(x = "Tourism Intensity", y = "Count") + theme_classic()

fig3

# Calculate Correlations

cor(log(df$tourism_intensity), log(df$gdp_per_capita_ppp))

[1] 0.7201282

cor(log(df$tourism_intensity), log(df$headcount + 0.001))

[1] -0.637536