Exploratory Data Analysis 1

R Packages used

• tidyverse: for general data wrangling (includes readr and dplyr)

• ggplot2: to make statistical plots; we use this rather than base R for increased functionality and more aesthetically pleasing plots (also part of tidyverse)

• ggthemes: additional themes for use with ggplot

• Hmisc: contains a LOWESS smoother for ggplot

• gap: to run the chow test

R Commands used

Below follows a list of the commands used in this notebook. For further details and a comprehensive list of options, please consult the R documentation.

• Base R: install.packages, library, head, names, summary, range, var, sd,pdf,dev.off,saveRDS,readRDS, function, lm, str, dim

• tidyverse: read_csv, rename, mutate, if_else, filter

• ggplot2: ggplot, geom_histogram, bins, theme_classic, theme_minimal, xlab, ylab, ggtitle, theme, layer_data, ggsave, geom_boxplot, stat_boxplot, geom_point, coord_fixed, geom_smooth, stat_smooth, labs, scale_color_manual

• ggthemes: theme_tufte

• Hmisc: stat_plsmo

• gap: chow.test

Creating an initial data frame

We use the tidyverse function read_csv to read the data into a data frame nyc.data. We could also have used the base R read.csv, but read_csv is a bit more robust and creates a tibble, a data frame with some additional information. As usual, we check the contents of the data frame with a head command.

nyc.data <- read_csv("nyc/nyc.csv")

## Parsed with column specification:
## cols(
##   .default = col_double(),
##   bor_subb = col_integer(),
##   NAME = col_character(),
##   CODE = col_integer(),
##   SUBBOROUGH = col_character(),
##   RENT2002 = col_integer(),
##   RENT2005 = col_integer(),
##   RENT2008 = col_integer()
## )

## See spec(...) for full column specifications.

head(nyc.data)

## # A tibble: 6 x 34
##   bor_subb NAME   CODE SUBBOROUGH FORHIS06 FORHIS07 FORHIS08 FORHIS09
##      <int> <chr> <int> <chr>         <dbl>    <dbl>    <dbl>    <dbl>
## 1      501 Nort…   501 North Sho…     37.1     34.0    27.4      29.3
## 2      502 Mid-…   502 Mid-Island     28.0     18.1    24.0      31.2
## 3      503 Sout…   503 South Sho…     10.7     12.1     9.69     14.7
## 4      401 Asto…   401 Astoria        52.1     54.0    54.7      47.8
## 5      402 Sunn…   402 Sunnyside…     62.7     69.4    67.1      58.3
## 6      403 Jack…   403 Jackson H…     68.5     68.5    66.5      69.2
## # ... with 26 more variables: FORWH06 <dbl>, FORWH07 <dbl>, FORWH08 <dbl>,
## #   FORWH09 <dbl>, HHSIZ1990 <dbl>, HHSIZ00 <dbl>, HHSIZ02 <dbl>,
## #   HHSIZ05 <dbl>, HHSIZ08 <dbl>, KIDS2000 <dbl>, KIDS2005 <dbl>,
## #   KIDS2006 <dbl>, KIDS2007 <dbl>, KIDS2008 <dbl>, KIDS2009 <dbl>,
## #   RENT2002 <int>, RENT2005 <int>, RENT2008 <int>, RENTPCT02 <dbl>,
## #   RENTPCT05 <dbl>, RENTPCT08 <dbl>, PUBAST90 <dbl>, PUBAST00 <dbl>,
## #   YRHOM02 <dbl>, YRHOM05 <dbl>, YRHOM08 <dbl>

Making the variable names compatible

Note, that in contrast to GeoDa (where the dbf file is read), reading the csv file into a data frame results in almost all the variable names being in caps. We confirm this with a names command:

names(nyc.data)

##  [1] "bor_subb"   "NAME"       "CODE"       "SUBBOROUGH" "FORHIS06"  
##  [6] "FORHIS07"   "FORHIS08"   "FORHIS09"   "FORWH06"    "FORWH07"   
## [11] "FORWH08"    "FORWH09"    "HHSIZ1990"  "HHSIZ00"    "HHSIZ02"   
## [16] "HHSIZ05"    "HHSIZ08"    "KIDS2000"   "KIDS2005"   "KIDS2006"  
## [21] "KIDS2007"   "KIDS2008"   "KIDS2009"   "RENT2002"   "RENT2005"  
## [26] "RENT2008"   "RENTPCT02"  "RENTPCT05"  "RENTPCT08"  "PUBAST90"  
## [31] "PUBAST00"   "YRHOM02"    "YRHOM05"    "YRHOM08"

We now use the tidyverse rename function to turn the all-caps variables into lower case for the examples we will use. As in the GeoDa workbook, we only use three variables, kids2009, kids2000, and pubast00.

nyc.data <- nyc.data %>% rename("kids2009" = "KIDS2009", "kids2000" = "KIDS2000",
                                "pubast00" = "PUBAST00")
names(nyc.data)

##  [1] "bor_subb"   "NAME"       "CODE"       "SUBBOROUGH" "FORHIS06"  
##  [6] "FORHIS07"   "FORHIS08"   "FORHIS09"   "FORWH06"    "FORWH07"   
## [11] "FORWH08"    "FORWH09"    "HHSIZ1990"  "HHSIZ00"    "HHSIZ02"   
## [16] "HHSIZ05"    "HHSIZ08"    "kids2000"   "KIDS2005"   "KIDS2006"  
## [21] "KIDS2007"   "KIDS2008"   "kids2009"   "RENT2002"   "RENT2005"  
## [26] "RENT2008"   "RENTPCT02"  "RENTPCT05"  "RENTPCT08"  "PUBAST90"  
## [31] "pubast00"   "YRHOM02"    "YRHOM05"    "YRHOM08"

Selecting the number of histogram bins

The standard way in ggplot is to adjust the number of bins indirectly, by means of the binwidth option, i.e., the range of values that make up a bin, in the units of the variable under consideration. To keep the parallel with the GeoDa workbook, we instead use the option bins, which sets the number of bins directly. The resulting histogram now matches the one in GeoDa (except for the lack of color, which is immaterial).

ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7)

As in the GeoDa workbook, we can now change the number of bins to 5, which yields the following histogram.

ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=5)

Spiffing up the graph

The graph as shown is just rudimentary. There are many options in ggplot to change the appearance of the graph, too many to cover here. But to illustrate some basic features, below, we add a label for the x and y axes using xlab and ylab, and a title for the graph with ggtitle. An unfortunate aspect of the latter is that it left aligns the text, whereas we would typically want it to be centered over the graph.

We can adjust this using the very powerful theme option. But first the basics. Every graph has a theme, which sets the main parameters for its appearance. The default theme with the grey grids, separated by white lines is theme_grey( ). If we want to change this, we can specify one of the other themes. For example, a classic graph a la base R plot, without background shading or grid lines is theme_classic( ). In order to obtain this specialized look, we set the associated theme command. Our histogram in this theme looks as follows, with a label on the x and y axis, and a title (and back to 7 bins).

ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme_classic()

There are seven built-in themes as well as several contributed ones. Another built-in example is theme_minimal( ), shown next.

ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme_minimal()

In addition, the package ggthemes contains several additional themes that look extremely professional. For example, theme_tufte( ).

ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme_tufte()

Besides selecting a different default theme, we can also override the basic settings associated with the current theme. For example, we adjust the plot.title (of course, you need to know what everything is called). Specifically, we set the element_text property’s horizontal justification (hjust) to 0.5. This centers the title. The number of other refinements is near infinite. Again, using the default theme_grey( ):

ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme(plot.title = element_text(hjust = 0.5))

Histogram summary statistics

The histogram in GeoDa has an option to display the contents of each bin as well as some descriptive statistics.

As anything in R, the plot created by ggplot is nothing but an object. When we enter the commands as above, starting with ggplot, the result is drawn directly to the screen. But we can also assign the plot object to a variable. This variable will then contain all the information needed to draw the graph, which includes the count of observations in each bin, the min and max values for each bin, etc. For example, we can assign our histogram plot to the plot.data object, and then extract the information using the layer_data function. The result is a data frame.

plot.data <- ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme(plot.title = element_text(hjust = 0.5))
layer_data(plot.data)

##    y count       x    xmin    xmax     density     ncount   ndensity PANEL
## 1  1     1  0.0000 -4.0109  4.0109 0.002266551 0.05555556 0.05555556     1
## 2  2     2  8.0218  4.0109 12.0327 0.004533102 0.11111111 0.11111111     1
## 3  4     4 16.0436 12.0327 20.0545 0.009066204 0.22222222 0.22222222     1
## 4  8     8 24.0654 20.0545 28.0763 0.018132407 0.44444444 0.44444444     1
## 5 18    18 32.0872 28.0763 36.0981 0.040797917 1.00000000 1.00000000     1
## 6 17    17 40.1090 36.0981 44.1199 0.038531366 0.94444444 0.94444444     1
## 7  5     5 48.1308 44.1199 52.1417 0.011332755 0.27777778 0.27777778     1
##   group ymin ymax colour   fill size linetype alpha
## 1    -1    0    1     NA grey35  0.5        1    NA
## 2    -1    0    2     NA grey35  0.5        1    NA
## 3    -1    0    4     NA grey35  0.5        1    NA
## 4    -1    0    8     NA grey35  0.5        1    NA
## 5    -1    0   18     NA grey35  0.5        1    NA
## 6    -1    0   17     NA grey35  0.5        1    NA
## 7    -1    0    5     NA grey35  0.5        1    NA

The convention used to create the histogram in ggplot is slightly different from that in GeoDa, hence small differences in the bounds of the bins. The summary statistics give the number of observations in each bin (count), the mid-point of the bin (x), and the lower and upper bound for the bin (xmin and xmax). Unlike GeoDa, where the histogram starts at the minimum value and ends at the maximum value, the histogram in ggplot starts with the minimum value at the mid-point of the lowest bin, and the maximum value at the mid-point of the upper bin.

Assigning (part of) a graph to an object

Any subset of ggplot commands can be assigned to an object, which can save on some typing if the same data set and variables are used for several plots. For example, we assign the main ggplot command with the geom_histogram to the object baseplt. As such, this does not draw anything. Next, we add the different options to the baseplt object, and the graph appears.

baseplt <- ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) 
baseplt +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme(plot.title = element_text(hjust = 0.5))

Other descriptive statistics

The usual descriptive statistics can be displayed by means of the base R summary command. In principle, we could assign these to an object and then add them to the plot using the geom_text geom, but that is beyond the current scope. We can easily obtain the descriptive statistics provided by GeoDa that are not contained in the R summary command, by means of range, var, and sd, for the range, variance and standard deviation.

summary(nyc.data$kids2009)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   26.69   33.53   32.07   39.68   48.13

range(nyc.data$kids2009)

## [1]  0.0000 48.1308

var(nyc.data$kids2009)

## [1] 107.535

sd(nyc.data$kids2009)

## [1] 10.36991

Writing the graph to a file

In our discussion so far, the graphs are drawn to the screen and then disappear. To save a ggplot graph to a file for publication, there are two ways to proceed. One is the classic R approach, in which first a device is opened, e.g., by means of a pdf command, then the plot commands are entered, and finally the device is turned off by means of dev.off().

Note that it is always a good idea to specify the dimension of the graph (in inches). If not, the results can be unexpected.

pdf("hist.pdf",height=3,width=3)
ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme(plot.title = element_text(hjust = 0.5))
dev.off()

## quartz_off_screen 
##                 2

In addiiton to the standard R approach, ggplot also has the ggsave command, which does the same thing. It requires the name for the output file, but derives the proper format from the file extension. For example, an output file with a png file extension will create a png file, and similarly for pdf, etc.

The second argument specifies the plot. It is optional, and when not specified, the last plot is saved. Again, it is a good idea to specify the width and height (in inches). In addition, for raster files, the dots per inch (dpi) can be set as well. The default is 300, which is fine for most use cases, but for high resolution graphs, one can set the dpi to 600, as in the example below.

hist.plot <- ggplot(data=nyc.data,aes(kids2009)) +
  geom_histogram(bins=7) +
  xlab("Percent kids in 2009") +
  ylab("Frequency") +
  ggtitle("Example Histogram") +
  theme(plot.title = element_text(hjust = 0.5))
ggsave("hist.png",plot=hist.plot,width=3,height=3,dpi=600)

Saving the graph object

An alternative approach to keep a plot object is to assign the plot commands to a variable and then save this to disk, using the standard R command saveRDS. This can later be brought back into an R session using readRDS. To save the plot, we need to specify a file name with an .rds file extension.

saveRDS(hist.plot,"hist.rds")

At some later point (or in a different R session), we can then read the object and plot it. Note that we do not need to assign it to the same variable name as before. For example, here we call the graph object newplot.

newplot <- readRDS("hist.rds")
newplot

Default settings

The minimal arguments to create a boxplot are the data set and the x and y variables passed to aes. As mentioned above, the logic behind the graphs in ggplot is two-dimensional, so both x and y need to be specified. The x variable is used to create separate box plots for different subsets of the data. In our simple example, we don’t need this feature, so we set the x variable to empty, i.e., " “. The y variable is the actual variable of interest, kids2009. The resulting graph is shown below.

ggplot(data=nyc.data,aes(x="",y=kids2009)) +
  geom_boxplot()

The box encloses the first and third quartile and shows the median as a thick horizontal bar. The vertical lines show the range of data, not the extent of the fences, as is the case in GeoDa. The single dot at value 0 is the lower outlier.

Box plot statistics

Unlike what is the case in GeoDa, there is no straightforward way to show the descriptive statistics on the graph. However, we can access the statistics, extract them, and then use text labeling techniques to add them to the plot. We will not pursue that here, but we will illustrate the type of statistics associated with the box plot.

As we did for the histogram, we will assign the ggplot object to a variable and then use the layer_data function to extract the information.

box.plt <- ggplot(data=nyc.data,aes(x="",y=kids2009)) +
  geom_boxplot()
box.dta <- layer_data(box.plt)
box.dta

##     ymin    lower  middle   upper    ymax outliers notchupper notchlower x
## 1 8.6623 26.69425 33.5284 39.6773 48.1308        0    36.2944    30.7624 1
##   PANEL group ymin_final ymax_final  xmin  xmax xid newx new_width weight
## 1     1     1          0    48.1308 0.625 1.375   1    1      0.75      1
##   colour  fill size alpha shape linetype
## 1 grey20 white  0.5    NA    19    solid

The result is a data frame that contains all the information needed to create the graph.

The descriptive statistics require a little clarification. The values for lower and upper are, respectively, the values for the first and third quartile, and middle is the median. ymin (8.6623) and ymax (48.1308) are not the smallest and largest values overall, but the smallest and largest values that fall inside the fences.⁷ They are the begin and end points of the vertical lines in the plot. outliers contains a list with the outlier values. In our example, there is just one, the value 0 (compare to Figure 18 in the GeoDa workbook). The fences, also sometimes called whiskers, are not contained among the statistics, but they can be easily computed.

We illustrate a simple function to accomplish this (note that this function does not implement any error checking and is purely illustrative of the concepts involved). As anything else in R, a function is an object that is assigned to a name. It takes arguments and returns the result. For example, the function box.desc given below takes the layer_data object as the argument box.lyr, extracts the quartiles to compute the interquartile range and to calculate the fences. We also pass the multiplier as mult, with a default value of 1.5 (the default value is set by the equal sign).

box.desc <- function(box.lyr,mult=1.5) {
  # function to computer lower and upper fence in a box plot
  # box.lyr: a box plot layer_data object
  # mult: the multiplier for the fence calculation, default = 1.5
  iqr <- box.lyr$upper - box.lyr$lower # inter-quartile range
  upfence <- box.lyr$upper + mult * iqr  # upper fence
  lofence <- box.lyr$lower - mult * iqr  # lower fence
  return(c(lofence,upfence))
}

We can now pass the box.dta results to this function to obtain the lower and upper fences.

box.desc(box.dta)

## [1]  7.219675 59.151875

These results match the values in the GeoDa workbook.

Changing the fence

As we did in the GeoDa workbook, we can change the multiplier value to compute the fences. The default is 1.5, but we can set this to 3.0 by means of the coef option. We again assign the plot to an object to both illustrate the graph and the associated statistics.

box.plt3 <- ggplot(data=nyc.data,aes(x="",y=kids2009)) +
  geom_boxplot(coef=3)
box.plt3

As in the GeoDa example, there are no longer any outliers. In the graph, the lowest point on the vertical line now corresponds with the value of 0.

We extract the statistics as before.

box.dta3 <- layer_data(box.plt3)
box.dta3

##   ymin    lower  middle   upper    ymax outliers notchupper notchlower x
## 1    0 26.69425 33.5284 39.6773 48.1308             36.2944    30.7624 1
##   PANEL group ymin_final ymax_final  xmin  xmax xid newx new_width weight
## 1     1     1          0    48.1308 0.625 1.375   1    1      0.75      1
##   colour  fill size alpha shape linetype
## 1 grey20 white  0.5    NA    19    solid

The statistics are the same, except that the value for ymin is now 0. We double check the results for the fences using our function box.desc, but now pass box.dta3 and set mult=3.0.

box.desc(box.dta3,mult=3.0)

## [1] -12.25490  78.62645

Since the lower fence is negative, the value of 0 is no longer an outlier.

Fancier options

As is, the default box plot is pretty rudimentary. We will illustrate the power of ggplot by adding a number of features to the plot in order to mimic the visual representation given in GeoDa. First, we will remove the label for the x-axis by setting it to “”, and add a title using ggtitle. As we did for the histogram, we will center the title over the graph. To save on some typing, we will assign the ggplot command with its arguments to the variable base.plt and then build the graph by adding layers. First, just the labels.

base.plt <- ggplot(data=nyc.data,aes(x="",y=kids2009))
base.plt + geom_boxplot() +
  xlab("") +
  ggtitle("Example Box Plot") +
  theme(plot.title = element_text(hjust=0.5))

Next, we want to give the box plot a color, as in GeoDa. This is accomplished with the color (for the outlines) and fill (for the inside of the box) options to geom_boxplot. In addition, we can give the outlier point a different color by means of outlier.color. For example, setting the color to black, with the fill to purple (if we set both to the same color, we can no longer distinguish the median), and the outlier.color to red, we obtain:

base.plt + 
  geom_boxplot(color="black",fill="purple",outlier.color="red") +
  xlab("") +
  ggtitle("Example Box Plot") +
  theme(plot.title = element_text(hjust=0.5))

So far, the box plots do not show the fences, the way they do in GeoDa. This can be remedied, but not quite in the same way as in GeoDa. In ggplot, the fences are drawn at the location of the extreme values, the ymin and ymax we saw before, and not at the location of the fence cut-off values, as in GeoDa. The fences are obtained from the stat_boxplot function, by passing the geom as errorbar. The result is as shown below.

base.plt + 
  geom_boxplot(color="black",fill="purple",outlier.color="red") +
  stat_boxplot(geom="errorbar") +
  xlab("") +
  ggtitle("Example Box Plot") +
  theme(plot.title = element_text(hjust=0.5))

One final refinement. In GeoDa, the box plot also shows the locations of the actual observations as points on the central axis. We obtain the same effect by adding geom_point with color blue. We draw the points first, and the box plot on top of it, using the layers logic. However, we want to make sure that the central box doesn’t mask the points, which it does when the transparency is kept as the default. To accomplish this, we set the alpha level for both points and box plot at 0.5.

The result comes as close to the GeoDa visualization as we can get without going overboard.

base.plt + 
  geom_point(color="blue",alpha=0.5) +
  geom_boxplot(color="black",fill="purple",outlier.color="red",alpha=0.5) +
  stat_boxplot(geom="errorbar") +
  xlab("") +
  ggtitle("Example Box Plot") +
  theme(plot.title = element_text(hjust=0.5))

As before, we can write this plot to a file using ggsave, or save the object for future use with saveRDS.

Linear smoother

The linear smoother is added to the plot by including the geom_smooth command after the geom_point call. The method is set as lm. In order to better distinguish the fitted line from the points, we set its color to blue, and add a centered title to specify the smoothing algorithm using ggtitle (also, in what follows, we ignore the aspect ratio issue).

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   geom_point() +
   geom_smooth(method=lm, color="blue") +
   ggtitle("Linear Smoother") +
   theme(plot.title = element_text(hjust=0.5))

By default, the linear smoother includes a 95% confidence interval band (the grey band around the blue line). To turn this off, we need to set the option se to FALSE, as below.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   geom_point() +
   geom_smooth(method=lm, color="blue",se=FALSE) +
   ggtitle("Linear Smoother") +
   theme(plot.title = element_text(hjust=0.5))

Extracting the linear regression results

In GeoDa, a linear fit to the scatter plot also yields the results of the regression, such as the coefficients, their standard errors, p-values and an R² measure of fit. This is not the case in ggplot (by design).

However, we can readily obtain these results from the lm function in base R. In its bare minimum, this function takes as arguments a formula and a data set (actually, the latter is not absolutely necessary, depending on how the variables are referred to). In our example, we specify the regression formula as pubast00 ~ kids2000 and the data set as nyc.data.⁸

Rather than just printing the results of the regression, we assign it to an object, reg1 in our example below. We then apply the summary command to this object, which we assign to yet another object (reg1.sum). When we list the latter, we see a summary of the regression results.

reg1 <- lm(pubast00 ~ kids2000, data=nyc.data)
reg1.sum <- summary(reg1)
reg1.sum

## 
## Call:
## lm(formula = pubast00 ~ kids2000, data = nyc.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.5284 -3.7925 -0.2762  3.6696  9.2054 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.61829    2.13013  -2.638   0.0109 *  
## kids2000     0.39000    0.05645   6.909 6.32e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.682 on 53 degrees of freedom
## Multiple R-squared:  0.4739, Adjusted R-squared:  0.4639 
## F-statistic: 47.73 on 1 and 53 DF,  p-value: 6.322e-09

The reason for taking what may seem like a circuitous route is that the summary object is simply a list of elements that correspond to aspects of the regression result. Of course, one needs to know what these are called, but, in doubt, a str command will reveal the full structure.

str(reg1.sum)

## List of 11
##  $ call         : language lm(formula = pubast00 ~ kids2000, data = nyc.data)
##  $ terms        :Classes 'terms', 'formula'  language pubast00 ~ kids2000
##   .. ..- attr(*, "variables")= language list(pubast00, kids2000)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "pubast00" "kids2000"
##   .. .. .. ..$ : chr "kids2000"
##   .. ..- attr(*, "term.labels")= chr "kids2000"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(pubast00, kids2000)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "pubast00" "kids2000"
##  $ residuals    : Named num [1:55] -3.703 -6.222 -8.528 -0.276 -3.061 ...
##   ..- attr(*, "names")= chr [1:55] "1" "2" "3" "4" ...
##  $ coefficients : num [1:2, 1:4] -5.6183 0.39 2.1301 0.0564 -2.6375 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "(Intercept)" "kids2000"
##   .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
##  $ aliased      : Named logi [1:2] FALSE FALSE
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "kids2000"
##  $ sigma        : num 4.68
##  $ df           : int [1:3] 2 53 2
##  $ r.squared    : num 0.474
##  $ adj.r.squared: num 0.464
##  $ fstatistic   : Named num [1:3] 47.7 1 53
##   ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf"
##  $ cov.unscaled : num [1:2, 1:2] 0.206953 -0.005238 -0.005238 0.000145
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "(Intercept)" "kids2000"
##   .. ..$ : chr [1:2] "(Intercept)" "kids2000"
##  - attr(*, "class")= chr "summary.lm"

For example, we see that the coefficients are in the list element coefficients. They can be extracted by means of the standard $ notation.

reg1.sum$coefficients

##               Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) -5.6182932 2.13013007 -2.637535 1.093595e-02
## kids2000     0.3899956 0.05644857  6.908866 6.322337e-09

Similarly, we can extract the R² and adjusted R².

c(reg1.sum$r.squared,reg1.sum$adj.r.squared)

## [1] 0.4738536 0.4639263

We can now use the text and labeling functionality of ggplot to place these results on the graph. We don’t pursue this any further.

Loess smoother

The default nonlinear smoother in ggplot uses the loess algorithm as a locally weighted regression model. This is similar in spirit to the LOWESS method used in GeoDa, but not the same.⁹ The implementation is along the same lines as the linear smoother, using geom_smooth, with the only difference that the method is now loess, as shown below.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   geom_point() +
   geom_smooth(method=loess, color="blue") +
   ggtitle("Loess Smoother") +
   theme(plot.title = element_text(hjust=0.5))

As in any local regression method, an important parameter is how much of the data is used in the local fit, the so-called span. This is typically set to 2/3 of the data by default. A narrower span will yield a smoother that emphasizes local changes. In order to set the span parameter, we use the stat_smooth command. It is in all respects equivalent to geom_smooth, but allows for somewhat more flexibility.

For example, we see the difference between the default and a smoother with a span = 0.4. We also turn off the confidence interval by setting se = FALSE.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   stat_smooth(method="loess",span=0.4,color="blue",se=FALSE) +
   geom_point() +
   ggtitle("Loess Smoother - Span=0.4") +
   theme(plot.title = element_text(hjust=0.5))

And even more with a span = 0.2.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   stat_smooth(method="loess",span=0.2,color="blue",se=FALSE) +
   geom_point() +
   ggtitle("Loess Smoother - Span=0.2") +
   theme(plot.title = element_text(hjust=0.5))

LOWESS smoother

The LOWESS smoother is not implemented in ggplot, but can be found in the Hmisc package. The approach taken illustrates an alternative way to compute a smoother, similar to the stat_smooth function in ggplot. The latter is equivalent to geom_smooth, but allows for non-standard geoms to visualize the results. We won’t pursue that here, but want to point out that this is the spirit in which the function stat_plsmo is implemented. The default is to use a span of 2/3 of the observations (the stat_plsmo command does not have a confidence interval option, so that we do not need to set se).

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   stat_plsmo(color="blue") +
   geom_point() +
   ggtitle("LOWESS Smoother") +
   theme(plot.title = element_text(hjust=0.5))

When compared to the loess results, we can distinguish minor differences. These become more pronounced when setting the span=0.4.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   stat_plsmo(span=0.4,color="blue") +
   geom_point() +
   ggtitle("LOWESS Smoother - Span=0.4") +
   theme(plot.title = element_text(hjust=0.5))

And even more when setting the span=0.2.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   stat_plsmo(span=0.2,color="blue") +
   geom_point() +
   ggtitle("LOWESS Smoother - Span=0.2") +
   theme(plot.title = element_text(hjust=0.5))

Putting it all together

In GeoDa, it is easy to show both the linear and LOWESS smoother on the same plot. In order to accomplish the same in ggplot, we will use some of the really powerful customization options to create a graph that contains all three smoothing methods. We distinguish between them by setting the color argument in aes to the name of the method, in essence a constant (not a variable).

In other words, the color argument to aes is not set to a variable, where it would take a different color depending on the value of that variable, but to a constant, where the color is fixed. This will yield a different color for each of the methods. In order to highlight the curves themselves, we turn off se for lm and loess. The plot will include a legend by default, showing the color with the name of the method. The default title of the legend will be color, i.e., the argument used to reflect the categories. We override this by means of the labs command (for legend labeling), and set color="Method". The result is as shown below.

ggplot(data=nyc.data,aes(x=kids2000,y=pubast00)) +
   stat_plsmo(aes(color="lowess")) +
   geom_point() +
   geom_smooth(aes(color="lm"),method=lm,se=FALSE) +
   geom_smooth(aes(color="loess"),method=loess,se=FALSE) +
   ggtitle("Comparison of Smoothing Methods") +
   theme(plot.title = element_text(hjust=0.5)) +
   labs(color="Method")