#' ---
#' title: "Data Visualization in R"
#' author: "Neil S. Williams^[University of Georgia, snpwill@uga.edu]"
#' date: |
#'   | Data Visualization in R
#'   | September 10, 2021
#' output:
#'   pdf_document:
#'     toc: yes
#'     toc_depth: '4'
#'   html_document:
#'     df_print: paged
#'     highlight: tango
#'     theme: united
#'     toc: yes
#'     toc_depth: 4
#'     toc_float: yes
#'   word_document:
#'     toc: yes
#'     toc_depth: '4'
#' ---
#' 
#' \newpage
#' 
#' \frenchspacing
#' 
## ----setup, include=FALSE------------------------------------------------------------------------------------------
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 

#' 
#' 
#' # Getting started
#' 
#' The purpose of this tutorial is to show how to display and visualize data using R. We will start with plots of means and distributions, then move to regression, and finally some more advanced topics. 
#' 
#' 
#' Please copy code from the R script used to produce this tutorial; [this script can be found here](http://www.neilswilliams.com/uploads/1/2/1/9/121979947/welcome_r.r).
#' 
#' 
#' Let's start by loading some base packages that are always useful when working in RStudio:
#' 
## ------------------------------------------------------------------------------------------------------------------
library(tidyverse)
library(dplyr)
library(ggplot2)
library(knitr)


#' 
#' 
#' -Let's start with frog weight
#' 
#' 
## ---- message=FALSE, warning=FALSE, echo = TRUE--------------------------------------------------------------------
library(rio)
frog_weight.data <- rio::import("frog_weight.dta")

#' 
#' Let's check the structure of these data
#' 
#' 
#' ## Examine data
#' 
## ------------------------------------------------------------------------------------------------------------------
#Check the structure of the data with the str() command
str(frog_weight.data) #note that "Lake" iscoded as 1,2,3 or Erie, Huron, Ontario

#' 
#' Much of data visualization is prepping the data to make sure it looks how you want it to in the plot. To start, I am going to recode values of the ```Lake``` variable so we know which number is associated with which lake.
#' 
## ------------------------------------------------------------------------------------------------------------------

#We can recode this and set the factor level becauese Lake is coded as a factor

frog_weight.data$lake[frog_weight.data$lake== 1] <- "Erie"
frog_weight.data$lake[frog_weight.data$lake== 2] <- "Huron"
frog_weight.data$lake[frog_weight.data$lake== 3] <- "Ontario"

#Save lake as a factor and relevel the variable
frog_weight.data$lake <- as.factor(frog_weight.data$lake)
frog_weight.data$lake <- relevel(frog_weight.data$lake, ref = "Erie")

summary(frog_weight.data$lake) # shows the number of cases of frogs in each lake



#' 
#' 
#' 
#' # Showing this with the ```table1``` package
#' 
#' 
## ---- message=FALSE------------------------------------------------------------------------------------------------
library(table1)
table1::table1(~.| lake, data = frog_weight.data) #the dot means we want to include all variables



#' 
#' \bigskip
#' 
#' We also may just want to select a certain set of variables 
#' 
#' 
#' 
## ---- tidy = TRUE--------------------------------------------------------------------------------------------------
 #we could also just select the ones we want
table1::table1(~frog_weight| lake, data = frog_weight.data)


#' 
#' 
#' 
#' 
#' 
#' # Plotting means 
#' 
#' 
#' We can plot the mean for each group with 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(gplots)
# Plot the mean of teeth length by dose groups
plotmeans(frog_weight ~ lake, data = frog_weight.data, 
          p = .84, # choose confidence level
          connect = FALSE, #whether the points should be connected
          ylab = "Frog Weight", 
          xlab = "Lake")


#' 
#' 
#' 
#' 
#' 
#' # Plotting distributions 
#' 
#' 
#' 
#' ## Boxplots
#' 
#' 
#' We can make a basic boxplot showing the range, quartile range, and mean of a variable using the ```boxplot()``` command
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------

boxplot(frog_weight.data$frog_weight)


#' 
#' 
#' We can spruce up this figure 
## ------------------------------------------------------------------------------------------------------------------
boxplot(frog_weight.data$frog_weight, #variable to use
main = "Boxplot of frog weight", #main title
xlab = "Weight of frog", #label on the x axis
ylab = "Frog weight", #label on the y axis
col = "orange", #fill color
border = "brown", #color of the border of the box
horizontal = TRUE, #make the plot horizontal
notch = TRUE #Creates triangular notch by mean
)



#' 
#' Use boxplot to compare between lakes
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
boxplot(frog_weight~lake, #frog weight plotted for each lake
data=frog_weight.data, #data
main="Boxplots of frog weight by lake",
xlab="Frog weight",
ylab="Lake",
col="orange",
border="brown",
horizontal = TRUE, 
notch = TRUE
)

#' 
#' ## Histograms 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
hist(frog_weight.data$frog_weight)




#' 
## ------------------------------------------------------------------------------------------------------------------
hist(frog_weight.data$frog_weight, breaks=10, col="red")





#' 
## ------------------------------------------------------------------------------------------------------------------
# Add a Normal Curve (Thanks to Peter Dalgaard)
x <- frog_weight.data$frog_weight
h<-hist(x, breaks=10, col="red", xlab="Frog weight",
   main="Histogram with Normal Curve")
xfit<-seq(min(x),max(x),length=40)
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x)) #
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)



#' 
#' 
#' ## Density plots
#' 
#' We can also plot the density as well 
#' 
## ------------------------------------------------------------------------------------------------------------------
plot(density(frog_weight.data$frog_weight))


#' 
## ------------------------------------------------------------------------------------------------------------------
d <- density(frog_weight.data$frog_weight)
plot(d, main="Kernel Density of Frog Weight")
polygon(d, col="red", border="blue") # adds fill for the density plot


#' 
#' 
#' ### Density Ridges
#' 
#' Comparisons using density ridges can be even more intuitive
#' 
## ------------------------------------------------------------------------------------------------------------------
library(ggridges)
library(ggplot2)
ggplot(data=frog_weight.data, aes(x =frog_weight, y = lake)) + #data and base aesthetics 
     theme_bw(12) +
  stat_density_ridges( quantile_lines = TRUE, 
                      quantiles= c(0.025, 0.5, 0.975), alpha= 0.7) +
  geom_vline(xintercept = mean(frog_weight.data$frog_weight), 
             linetype = "dashed", col = "red") +    #mean line
  xlab("Frog weight") +ylab("Lake") 



#' 
#' 
#' 
#' 
#' ## Points and lines
#' 
#' Using ```geom_point``` we can plot points 
#' 
## ------------------------------------------------------------------------------------------------------------------
ggplot(data = frog_weight.data) + 
  geom_point(mapping = aes(x =  num_flies, y = frog_weight)) # aes = aesthetics


#' 
#' We can show a smoothed line here from the data
#' 
## ------------------------------------------------------------------------------------------------------------------
ggplot(data = frog_weight.data) +
  geom_smooth(mapping = aes(x =  num_flies, y = frog_weight)) # aes = aesthetics

#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
ggplot(data = frog_weight.data) +
  geom_smooth(mapping = aes(x =  num_flies, y = frog_weight, group = lake)) # aes = aesthetics 

#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
ggplot(data = frog_weight.data) + #our data
  geom_smooth(
    mapping = aes(x = num_flies, y = frog_weight, color = lake), # color is our grouping variable
    show.legend = TRUE # say we want legend
  ) + 
  xlab("Number of flies") + #x-label
  ylab("Frog weight") #y-label

#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
ggplot(data = frog_weight.data) + #our data
    geom_point(mapping = aes(x =  num_flies, y = frog_weight, color = lake))+ # aes = aesthetics

  geom_smooth(
    mapping = aes(x = num_flies, y = frog_weight, color = lake), # color is our grouping variable
    show.legend = TRUE # say we want legend
  ) + 
  xlab("Number of flies") + #x-label
  ylab("Frog weight") #y-label

#' 
#' 
#' We can use ```ggsave()``` to save ```ggplot``` objects to a particular location
#' 
## ------------------------------------------------------------------------------------------------------------------

ggsave("flies_lakes_workshop.pdf", 
       width = 6.5, height = 4.5, units = "in")



#' 
#' 
#' 
#' # Regression
#' 
#' 
#' ## Estimate model
#' 
## ---- message=FALSE, warning=FALSE---------------------------------------------------------------------------------

frog_base.mod <- lm(frog_weight ~ num_flies + lake, data = frog_weight.data)
summary(frog_base.mod)


#' 
#' If using Rmarkdown, we can use ```tidy()``` and ```kable``` to clean this output up
## ------------------------------------------------------------------------------------------------------------------
library(broom)
kable(tidy(frog_base.mod))


#' 
#' 
#' ## Export results to Latex
#' 
#' To export regression results to latex, the ```stargazer``` package is fantastic!
#' 
#' 
## ---- results = "asis"---------------------------------------------------------------------------------------------
library(stargazer)
stargazer(frog_base.mod)


#' 
#' 
#' ## Plotting coefficients 
#' 
#' 
#' ```plot_summs``` is great for coefficient plots 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(jtools)
plot_summs(frog_base.mod, scale = TRUE)

#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
plot_summs(frog_base.mod, scale = TRUE, inner_ci_level = .9) 

#' 
#' Plotting normal distributions can also be interesting 
## ------------------------------------------------------------------------------------------------------------------
plot_summs(frog_base.mod, 
           scale = TRUE, 
           plot.distributions = TRUE) #plot normal distributions 

#' 
#' 
#' 
#' 
#' ## Effects plots
#' 
#' These are useful for showing the size of the main predicted effect and can also include your points
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(jtools)
jtools::effect_plot(frog_base.mod,
            pred = num_flies,  #main focal predictor
            interval = TRUE, #confidence interval
            plot.points = TRUE) #plot data points



#' 
#' 
#' 
#' 
#' 
#' # Interactions 
#' 
#' 
#' Interactive models are important parts of statistical modelling in social science. However, these models are difficult to interpret by just looking at the coefficients from model output. Given this, we should try to approach results from interactions graphically!
#' 
#' 
#' Let's take a look at a continuous $\times$ categorical interaction with the ```frog_weight.data```
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' ## Estimate interactive model
#' 
#' - For interactions in R all we need to do is use the ```*``` symbol between the two terms of interest. Using  ```*``` will include the constituent variables from the interactions in the model as well.
#' 
## ---- message=FALSE, warning=FALSE---------------------------------------------------------------------------------

frog.mod <- lm(frog_weight ~ num_flies*lake, data = frog_weight.data)
summary(frog.mod)


#' 
#' 
#' We can compare this interactive output with our output from before 
## ------------------------------------------------------------------------------------------------------------------

plot_summs(frog_base.mod, frog.mod, scale = TRUE) 
  

## ------------------------------------------------------------------------------------------------------------------

plot_summs(frog_base.mod, frog.mod, 
           scale = TRUE, 
           plot.distributions = TRUE)

#' 
#' 
#' - Now we can plot the predictions with the ```interact_plot()```. More details on this command can be found \href{https://www.rdocumentation.org/packages/jtools/versions/0.4.5/topics/interact_plot}{here}. Note that we have to set the ```interval``` option to ```TRUE``` if we want to have confidence intervals around our predictions.
#' 
#' 
## ---- message=FALSE, fig.align='center', warning=FALSE-------------------------------------------------------------
library(jtools)
library(interactions)

interact_plot(model = frog.mod, pred = num_flies, modx = lake, #modx refers to the modifier
              interval = TRUE, #Add a confidence interval
              legend.main = "Lake", #Set the name of the legend
              x.label = "Number of Flies", #label of the x-axis
              y.label = "Predicted Frog Weight", #label of the y-axis
              main.title = "Predicted Frog Weight \nby Number of Flies in the Great Lakes") 



#' 
#' 
#' ## Blood Pressure Analysis
#' 
#' 
#'  Now let's move to ```nhanes2```. With this, we can see more complicated interactions
#' 
#' 
## ---- message=FALSE, warning=FALSE, echo = TRUE--------------------------------------------------------------------
nhanes2 <- import("nhanes2.dta")

#' 
#' ### Examine Data
#' 
## ---- message=FALSE, warning=FALSE, echo = TRUE--------------------------------------------------------------------
str(nhanes2)

#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
#note that we need to make race and region factors
nhanes2$race <- as.factor(nhanes2$race)
nhanes2$region <- as.factor(nhanes2$region)

#We would have to recode these as well but let's skip that for now

#We will recode the female variable

nhanes2$female[nhanes2$female== 0] <- "Male"
nhanes2$female[nhanes2$female== 1] <- "Female"

#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
table1::table1(~bpsystol + race + weight + age + region + iron| female, data = nhanes2)

#' 
#' 
#' ### Estimate regression
#' 
#' - Let's start with an interaction between ```bmi``` and ```female```
## ------------------------------------------------------------------------------------------------------------------
blood1.reg <- lm(bpsystol ~ age + race + region + iron + bmi*female, data =nhanes2)
summary(blood1.reg)

#' 
#' If using Rmarkdown, we can use ```tidy()``` and ```kable``` to clean this output up
## ------------------------------------------------------------------------------------------------------------------

kable(tidy(blood1.reg))


#' 
#' ### Create more plots!
#' - Let's create a predicted values plot. Instead of showing the slopes of the two 
#' 
## ---- message=FALSE, fig.align='center', warning=FALSE, fig.width = 5.2, fig.height = 3.8--------------------------
interact_plot(model = blood1.reg, 
              pred = bmi,
              modx = female, #modx refers to the modifier
              interval = TRUE, # confidence interval
              legend.main = "Sex",
              x.label = "Body Mass Index (BMI)", 
              y.label = "Predicted Blood Pressure", #this is the dependent variable
              main.title = "Predicted Blood Pressure \nby BMI for Males and Females")



#' 
#' 
#' 
#' 
#' ### Estimate a continuous$\times$continous interaction
#' 
#' - Now we will try a continuous$\times$continous interaction. These are even more complicated to interpret. We will interact ```bmi``` with ```age```
#' 
## ------------------------------------------------------------------------------------------------------------------
blood2.reg <- lm(bpsystol ~ female + race + region + iron + bmi*age, data =nhanes2)
summary(blood2.reg)


#' 
#' 
#' 
#' 
#' ### Other interaction plots!
#' 
#' 
#' - The```interplot()``` command can be especially useful for continuous interactions
#' 
#' - Note that we set the ```hist``` option to ```TRUE``` to include a histogram of the values of the age data.
#' 
## ---- message=FALSE, fig.align='center', warning=FALSE, fig.width = 5, fig.height = 3------------------------------
library(interplot)
library(ggplot2)

interplot(m = blood2.reg, var1 = "bmi", var2 = "age", hist = TRUE) + 
  ggtitle("Conditional Effect of BMI on Blood Pressure \nat Different Values of Age") 

#' 
#' 
#' 
#' 
#' - Finally let's plot predictions for the interaction
#' 
#' 
## ---- message=FALSE, fig.align='center', warning=FALSE-------------------------------------------------------------
interact_plot(model = blood2.reg,
              pred = bmi, #main predictor on the x-axis
              modx = age, #modx refers to the modifier
              interval = TRUE, #Include confidence intervals around predictions
              legend.main = "Age", #Title of the legend
              x.label = "Body Mass Index (BMI)", #x-axis label
              y.label = "Predicted Blood Pressure", #y-axis label
              main.title = "Predicted Blood Pressure \nby Age")#use the \n to break to next line



#' 
#' 
#' - As we can see, ```interact_plot``` will by default plot lines for the \textit{mean} and $\pm$ \textit{1 standard deviation either side of the mean} of the moderator, which in this case is age. 
#' 
#' 
#' 
#' ### Exporting output to latex
#' 
#' 
#' # Advanced plotting 
#' 
#' 
#' R can do all sorts of plots and visualization that goes beyond basic regression. Here we'll look at census data and mapping. 
#' 
#' 
#' 
#' 
#' # Working with ```tidycensus```
#' 
#' 
#' 
#' ### Download census data
#' 
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
#install.packages("tidycensus")
library(tidycensus)

#' 
#' 
#' -Note that you need to submit to [get a census API key](https://api.census.gov/data/key_signup.html)
#' 
#' - After getting the key you can call it using the below function
#' 
## ------------------------------------------------------------------------------------------------------------------

census_api_key("bbf9b57e8eabfb053b114216983d8348d4e24540", overwrite = FALSE, install = FALSE)


Sys.getenv("CENSUS_API_KEY")

#' 
## ------------------------------------------------------------------------------------------------------------------
v17 <- load_variables(2017, "acs5", cache = TRUE) 


#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
country_2017 <- get_acs(geography = "state", 
              variables = c(medincome = "B19013_001", #selecting each code and renaming it
                            Population = "B01003_001",
                            White_pop = "B02008_001",
                            Black_pop = "B02009_001",
                            Latino_pop = "B03001_003",
                            Asian_pop = "B02015_001",
                            Native_pop = "B02014_001"
                           ), 
              year = 2017)

# Can also use get_decennial to get census data



#Can also get the data for 2018 and any other year

country_2018 <- get_acs(geography = "state", 
              variables = c(medincome = "B19013_001",
                            Population = "B01003_001",
                            White_pop = "B02008_001",
                            Black_pop = "B02009_001",
                            Latino_pop = "B03001_003",
                            Asian_pop = "B02015_001",
                            Native_pop = "B02014_001"
                           ), 
              year = 2018)




#' 
## ------------------------------------------------------------------------------------------------------------------
country_2017 %>%
  filter(variable == "Black_pop")%>%
  mutate(estimate = estimate/1000)%>%
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) +  #aes = aesthetics
  xlab("Black population (in thousands) ") + 
  ylab("State") + 

  geom_point() 
  



#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
country_2018 %>%
  filter(variable == "Latino_pop")%>%
  mutate(estimate = estimate/1000)%>%
  ggplot(aes(x = estimate, y = reorder(NAME, estimate)))  + 
  xlab("Latino population (in thousands) ") + 
  ylab("State") + 

  geom_point() 



#' 
#' 
#' ### Add margin of error
## ------------------------------------------------------------------------------------------------------------------
country_2018 %>%
  filter(variable == "medincome")%>%
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) +
  geom_errorbarh(aes(xmin = estimate - moe, xmax = estimate + moe)) +
  geom_point(color = "red", size = 1) +
  labs(title = "Household income by State",
       y = "State",
       x = "Household income (dollars)")


#' 
#' 
#' 
#' ## Maps with ACS and Census data 
#' 
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------

country_2017_year <- country_2017%>%
  mutate(Year = "2017")


country_2018_year <- country_2018%>%
  mutate(Year = "2018")

country.df <- rbind(country_2017_year, 
      country_2018_year)





#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
country.df %>%
  filter(variable == "medincome")%>%
  ggplot(aes(x = estimate, y = reorder(NAME, estimate))) +
  geom_errorbarh(aes(xmin = estimate - moe, xmax = estimate + moe)) +
      theme_bw(8) + 
  geom_point(color = "red", size = 1) +
  labs(title = "Household income by State",
       y = "State",
       x = "Household income (dollars)") + 
  facet_grid(.~Year)


#' 
#' 
#' ## Create proportions 
#' 
#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(reshape2)
country.wide <- recast(country.df, NAME + Year ~ variable, measure.var=c("estimate"))



#' 
## ------------------------------------------------------------------------------------------------------------------
country.prop <- country.wide%>%
group_by(NAME, Year)%>%
mutate(Black_prop = Black_pop/Population, 
       Latino_prop = Latino_pop/Population, 
       Asian_prop = Asian_pop/Population, 
       Native_prop = Native_pop/Population,
       White_prop = White_pop/Population
       
       )%>%
dplyr::select(State = NAME,
              Year, 
              Black = Black_prop,
              Latino = Latino_prop,
              Asian = Asian_prop,
              Native = Native_prop,
              White =  White_prop
)

#' 
#' 
#' Return back to long 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(tidyr)
country_prop.long <- tidyr::gather(country.prop, #old wide data
                            Race_prop, #variable name for the long column
                            Proportion, #variable name for values 
                            Black:White)#variables to combine



#' 
## ------------------------------------------------------------------------------------------------------------------
country_prop.long %>%
  ggplot(aes(x = Proportion, y = reorder(State, desc(State)))) +
    theme_bw(8) + 
  theme(legend.position="bottom",
        legend.title = element_blank(),
        panel.grid.minor.x = element_blank(), 
        panel.grid.minor.y = element_blank()) +
  geom_point(aes(shape=Race_prop, color=Race_prop),position=position_dodge(width=0)) + 
  labs(title = "Proportion of population by race",
       y = "State",
       x = "Proportion of population ") + 
  facet_grid(.~Year)


#' 
#' 
#' Add conditional text
#' 
#' 
## ---- eval = TRUE--------------------------------------------------------------------------------------------------

  ggplot(country_prop.long, aes(x = Proportion, y = reorder(State, desc(State)))) +
    theme_bw(8) + 
  theme(legend.position="bottom",
        legend.title = element_blank(),
        panel.grid.minor.x = element_blank(), 
        panel.grid.minor.y = element_blank()) +
  geom_point(aes(shape=Race_prop, color=Race_prop),position=position_dodge(width=0)) + 
  
        geom_text(data=country_prop.long[country_prop.long$Race_prop=='White',],
                aes(label = sprintf("%0.3f", Proportion)), 
              colour = "black", 
                 size = 2,
                position = position_nudge(x = .06)) + 
  
  labs(title = "Proportion of population by race",
       y = "State",
       x = "Proportion of population ") + 
  
  facet_grid(.~Year)


#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' 
#' # Maps 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(tmap)
library(sf)
library(raster)
library(dplyr)
library(spData)


#' 
#' 
#' 
#' 
#' ## US map
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
us_states_map = tm_shape(us_states, projection = 2163) + tm_polygons() + 
  tm_layout(frame = FALSE)

#' 
## ------------------------------------------------------------------------------------------------------------------
hawaii_map = tm_shape(hawaii) + tm_polygons() + 
  tm_layout(title = "Hawaii", frame = FALSE, bg.color = NA, 
            title.position = c("LEFT", "BOTTOM"))
alaska_map = tm_shape(alaska) + tm_polygons() + 
  tm_layout(title = "Alaska", frame = FALSE, bg.color = NA)

#' 
## ------------------------------------------------------------------------------------------------------------------
us_states_map
print(hawaii_map, vp = grid::viewport(0.35, 0.1, width = 0.2, height = 0.1)) #setting the location of Hawaii
print(alaska_map, vp = grid::viewport(0.15, 0.15, width = 0.3, height = 0.3)) #setting the location of alaska



#' 
#' 
#' We could put Alaska in the Gulf of Mexico if we wanted!! 
#' 
## ------------------------------------------------------------------------------------------------------------------
us_states_map
print(hawaii_map, vp = grid::viewport(0.35, 0.1, width = 0.2, height = 0.1)) #setting the location of Hawaii
print(alaska_map, vp = grid::viewport(0.6, 0.1, width = 0.3, height = 0.3)) #setting the location of alaska



#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
library(tmap)


urb_1970_2030 = urban_agglomerations %>% 
  filter(year %in% c(1970, 1990, 2010, 2030))

tm_shape(world) +
  tm_polygons() +
  tm_shape(urb_1970_2030) +
  tm_symbols(col = "black", border.col = "white", size = "population_millions") +
  tm_facets(by = "year", nrow = 2, free.coords = FALSE)

#' 
#' 
#' 
#' ## World maps
#' 
## ------------------------------------------------------------------------------------------------------------------
library(tmap)


world_coffee = left_join(world, coffee_data, by = "name_long")
facets = c("coffee_production_2016", "coffee_production_2017")
tm_shape(world_coffee) + tm_polygons(facets) + 
  tm_facets(nrow = 1, sync = TRUE)

#' 
#' 
#' ## Plot our census data
#' 
## ------------------------------------------------------------------------------------------------------------------
library(tmap)

country_prop.map <- country.prop%>% 
  filter(Year == 2018)%>% #only 2018
  mutate(NAME = State) #rename state to match with the map file

#' 
#' 
#' 
## ------------------------------------------------------------------------------------------------------------------

us_race_prop.map <- left_join(us_states, country_prop.map, by = "NAME") #combine the map and census data


facets = c("Black", "Latino", "Asian", "White") # set facets to be included


us.census.plot <- tm_shape(us_race_prop.map) + 
  tm_polygons(facets) + #what should be included-here we have four facets
tm_layout(legend.text.size = 0.5, # adjust layout
          legend.position = c("right","bottom"),
          legend.bg.color = "white",
          legend.bg.alpha = 1) + 
tm_facets(nrow = 2, ncol = 2)  #number of rows and columns 


#' 
#' 
## ------------------------------------------------------------------------------------------------------------------
print(us.census.plot)



#' 
#' 
#' 
#' 
#' - For more on this example of census mapping see here: [Walkthrough of mapping census data with ```tidycensus```](http://zevross.com/blog/2018/10/02/creating-beautiful-demographic-maps-in-r-with-the-tidycensus-and-tmap-packages/#download-census-data)
#' 
#' 
#' 
#' 
#' # Conclusions 
#' 
#' Always think: **"What information helps the reader or viewer interpret my data or results?"** A more complicated or advanced plot is not necessarily always better. In fact, these plots may confuse those who are not familiar with your method, approach, or data.
#' 
#' 
#' 
#' Other tips/notes: 
#' 
#' 
#' - Consider the size of text 
#' 
#' - Be careful with color. Are the colors too similar? Jarring? Are they visible to people with color blindness or other visual impairments?
#' 
#' - For the above point, combining both different colors and line/point types 
#' 
#' - Journals may not want color in your final manuscripts
#' 
#' 
#' 
#' 
#' 
#' 
#' # More Resources
#' 
#' - [Data visualization in R](https://r4ds.had.co.nz/data-visualisation.html)
#' - [Regression output visualization with ```jtools```](https://cran.r-project.org/web/packages/jtools/vignettes/summ.html#effect_plot)
#' - [Vignette and examples of Interplot](https://cran.r-project.org/web/packages/interplot/vignettes/interplot-vignette.html)
#' - [Plotting simple effects in regression models](https://jtools.jacob-long.com/reference/effect_plot.html)
#' - [An introduction to the ```margins``` command](https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html) 
#'     
## ---- echo = FALSE, include=FALSE, eval = FALSE--------------------------------------------------------------------
## library(knitr)
## 
## purl("Data_Visualization.Rmd", output = "Data_Visualization.R", documentation = 2)
## 

#'