Data_Products_Final_Project

Introduction

GGPlot

Barplots

Now lets take a look at some ggplot2 barplots

We’ll start with making a dataframe based on the tooth data.

df <- data.frame(dose = c("D0.5", "D1", "D2"),
                 len = c(4.2, 10, 29.5))
df

##   dose  len
## 1 D0.5  4.2
## 2   D1 10.0
## 3   D2 29.5

And now lets make a second dataframe

df2 <- data.frame(supp=rep(c("VC", "OJ"), each = 3),
                  dose = rep(c("D0.5", "D1", "D2"), 2),
                  len = c(6.8, 15, 33, 4.2, 10, 29.5))

df2

##   supp dose  len
## 1   VC D0.5  6.8
## 2   VC   D1 15.0
## 3   VC   D2 33.0
## 4   OJ D0.5  4.2
## 5   OJ   D1 10.0
## 6   OJ   D2 29.5

Lets load up ggplot2

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.4.3

Lets set our parameters for ggplot

theme_set(
  theme_classic() +
    theme(legend.position = "top")
)

Lets start with some nasic barplots using the tooth data

f <- ggplot(df, aes(x=dose, y=len))

f + geom_col()

Now lets change the fill, and add labels to the top

f + geom_col(fill = "darkblue") + 
  geom_text(aes(label = len), vjust = -0.3)

Now lets add the labels inside the bars

f + geom_col(fill = "darkblue") +
  geom_text(aes(label = len), vjust = 1.6, color = "white")

Now lets change the barplot colors by group

f + geom_col(aes(color = dose), fill = "white") +
  scale_color_manual(values = c("blue", "gold", "red"))

This is kinda hard to see, so lets change the fill.

f + geom_col(aes(fill = dose))+
  scale_fill_manual(values = c("blue", "gold", "red"))

Ok how do we do this with multiple groups

ggplot(df2, aes(x = dose, y=len)) +
  geom_col(aes(color = supp, fill = supp), position = position_stack()) +
  scale_color_manual(values = c("blue", "gold")) +
  scale_fill_manual(values = c("blue", "gold"))

p <- ggplot(df2, aes(x = dose, y = len)) +
  geom_col(aes(color = supp, fill = supp), position = position_dodge(0.8), width = 0.7) +
  scale_color_manual(values = c("blue", "gold")) +
  scale_fill_manual( values = c("blue", "gold"))

p

Now lets add those lables to the dodged barplot

p + geom_text(
  aes(label = len, group = supp),
  position = position_dodge(0.8),
  vjust = -0.3, size = 3.5
)

What if we want to add labels to our stacked barplots? For this we need dplyr

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

df2 <- df2 %>%
  group_by(dose) %>%
  arrange(dose, desc(supp)) %>%
  mutate(lab_ypos = cumsum(len) - 0.5 * len)

df2

## # A tibble: 6 × 4
## # Groups:   dose [3]
##   supp  dose    len lab_ypos
##   <chr> <chr> <dbl>    <dbl>
## 1 VC    D0.5    6.8      3.4
## 2 OJ    D0.5    4.2      8.9
## 3 VC    D1     15        7.5
## 4 OJ    D1     10       20  
## 5 VC    D2     33       16.5
## 6 OJ    D2     29.5     47.8

Now lets recreate our stacked graphs

ggplot(df2, aes(x=dose, y=len))+
  geom_col(aes(fill = supp), width = 0.7) +
  geom_text(aes(y = lab_ypos, label = len, group = supp), color = "white") +
  scale_color_manual(values = c("blue", "gold")) +
  scale_fill_manual(values = c("blue", "gold"))

Box Plots

Lets look at some boxplots

data("ToothGrowth")

Lets change the dose to a factor, and look at the top of the dataframe

ToothGrowth$dose <- as.factor(ToothGrowth$dose)

head(ToothGrowth, 4)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5

Lets load ggplot

library(ggplot2)

Lets set the theme for our plots to classic

theme_set(
  theme_bw() +
    theme(legend.position = "top"))

Lets start with a very basic boxplot with dose vs length

tg <- ggplot(ToothGrowth, aes(x = dose, y = len))
tg + geom_boxplot()

Now lets look at a boxplot with points for the mean

tg + geom_boxplot(notch = TRUE, fill = "lightgrey") +
  stat_summary(fun.y = mean, geom = "point", shape = 18, size = 2.5, color = "indianred")

## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

We can also change the scale number of variables included, and their order

tg + geom_boxplot() +
  scale_x_discrete(limits = c("0.5", "2"))

## Warning: Removed 20 rows containing missing values or values outside the scale range
## (`stat_boxplot()`).

Lets put our x axis in descending order

tg + geom_boxplot() +
  scale_x_discrete(limits = c("2", "1", "0.5"))

We can also change boxplot colors by groups

tg + geom_boxplot(aes(color = dose)) +
  scale_color_manual(values = c("indianred", "blue1", "green2"))

What if we want to display our data subset by oj vs vitamiin c?

tg2 <- tg + geom_boxplot(aes(fill = supp), position = position_dodge(0.9)) +
  scale_fill_manual(values = c("#999999", "#E69F00"))
tg2

We can also arrange this as two plots with facet_wrap

tg2 + facet_wrap(~supp)

Histograms

set.seed(1234)

wdata = data.frame(
  sex = factor(rep(c("F", "M"), each = 200)),
  weight = c(rnorm(200, 56), rnorm(200, 58))
)

head(wdata, 4)

##   sex   weight
## 1   F 54.79293
## 2   F 56.27743
## 3   F 57.08444
## 4   F 53.65430

Now lets load dplyr

library(dplyr)

mu <- wdata %>%
  group_by(sex) %>%
  summarise(grp.mean = mean(weight))

Now lets load the plotting package

library(ggplot2)

theme_set(
  theme_classic() +
    theme(legend.position = "bottom")
  )

Now lets create a ggplot object

a <- ggplot(wdata, aes(x = weight))

a+ geom_histogram(bins = 30, color = "black", fill = "grey") +
  geom_vline(aes(xintercept = mean(weight)),
             linetype = "dashed", size = 0.6)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Now lets change the color by group

a + geom_histogram(aes(color = sex), fill = "white", position = "identity") +
                     scale_color_manual(values = c("pink", "blue"))

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

a + geom_histogram(aes(color = sex, fill = sex), position = "identity") +
                     scale_color_manual(values = c("indianred", "blue")) +
  scale_fill_manual(values = c("pink","lightblue"))

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

What if we want to combine density plots and histograms?

a+ geom_histogram(aes(y= stat(density)),
                  color = "black", fill= "white") +
  geom_density(alpha = 0.2, fill = "#FF6666")

## Warning: `stat(density)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

a + geom_histogram(aes(y = stat(density), color = sex),
                   fill = "white", position = "identity") +
  geom_density(aes(color = sex), size = 1) +
  scale_color_manual(values = c("indianred", "lightblue"))

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Dot Plots

First lets load the required packages

library(ggplot2)

Lets set our theme

theme_set(
  theme_dark() +
    theme(legend.position = "top")
)

First lets initiate a ggplot object called TG

data("ToothGrowth")
ToothGrowth$dose <- as.factor(ToothGrowth$dose)

tg <- ggplot(ToothGrowth, aes(x=dose, y = len))

lets create a dotplot with a summary statistic

tg + geom_dotplot(binaxis = "y", stackdir = "center", fill = "white") +
  stat_summary(fun = mean, fun.args = list(mult=1))

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_segment()`).

Lets add a box plot and a dot plot together

tg + geom_boxplot(width = 0.5) +
  geom_dotplot(binaxis = "y", stackdir = "center", fill = "white")

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

tg + geom_violin(trim = FALSE) +
  geom_dotplot(binaxis = "y", stackdir = "center", fill = "#999999") +
  stat_summary(fun = mean, fun.args = list(mult = 1))

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_segment()`).

Lets create a dotplot with multiple groups

tg + geom_boxplot(width = 0.5) +
  geom_dotplot(aes(fill = supp), binaxis = "y", stackdir = "center") +
  scale_fill_manual(values = c("indianred", "lightblue1"))

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

tg + geom_boxplot(aes(color = supp), width = 0.5, position = position_dodge(0.8)) +
  geom_dotplot(aes(fill = supp, color = supp), binaxis = "y", stackdir = "center", 
               dotsize = 0.8, position = position_dodge(0.8)) +
  scale_fill_manual(values = c("#00AFBB", "#E7B800")) +
  scale_color_manual(values = c("#00AFBB", "#E7B800"))

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

Line Plots

Now lets change it up and look at some line plots

We’ll start by making a custom dataframe kinda like the tooth dataset

This way, we can see the lines and stuff that we’re modifying

df <- data.frame(dose = c("D0.5", "D1", "D2"),
                 len = c(4.2, 10, 29.5))

Now lets create a ssecond dataframe for plotting by groups

df2 <- data.frame(supp = rep(c("VC", "OJ"), each = 3),
                  dose = rep(c("D0.5", "D1", "D2"), 2),
                  len = c(6.8, 15, 33, 4.2, 10, 29.5))

df2

##   supp dose  len
## 1   VC D0.5  6.8
## 2   VC   D1 15.0
## 3   VC   D2 33.0
## 4   OJ D0.5  4.2
## 5   OJ   D1 10.0
## 6   OJ   D2 29.5

Now lets again load ggplot2 fand set a theme

library(ggplot2)

theme_set(
  theme_gray() +
    theme(legend.position = "right")
)

Now lets do some basic line plots. First we will build a function to display all the different line types

generateRLineTypes <- function(){
  oldPar <- par()
  par(font = 2, mar = c(0,0,0,0))
  plot(1, pch = "", ylim = c(0,6), xlim = c(0, 0.7), axes = FALSE, xlab = "", ylab = "")
  
  for(i in 0:6) lines(c(0.3, 0.7), c(i,i), lty = i, lwd = 3)
  text(rep(0.1, 6), 0:6, labels = c("0. 'Blank'", "1.'solid'", "2.'dashed'", "3.'dotted'", "4.'dotdash'", "5.'longdash'", "6.'twodash'"))
  
  par(mar = oldPar$mar, font = oldPar$font)
}

generateRLineTypes()

Now lets build a basic line plot

p <- ggplot(data = df, aes(x = dose, y = len, group = 1))

p + geom_line() + geom_point()

Now lets modify the line type and color

p + geom_line(linetype = "dashed", color= "steelblue") +
  geom_point(color = "lightblue")

Now lets try a step graph, which indicates a threshold type progression

p + geom_step() + geom_point()

Now lets move on to making multiple groups. First we’ll create our ggplot object

p <- ggplot(df2, aes(x = dose, y = len, group = supp))

Now lets change line types and point shapes by group

p + geom_line(aes(linetype = supp, color = supp)) +
  geom_point(aes(shape = supp, color = supp)) +
  scale_color_manual(values = c("red", "blue"))

Now lets look at line plots with a numeric x axis

df3 <- data.frame(supp = rep(c("VC", "OJ"), each = 3),
                  dose = rep(c("0.5", "1", "2"), 2),
                  len = c(6.8, 15, 33, 4.2, 10, 29.5))

df3

##   supp dose  len
## 1   VC  0.5  6.8
## 2   VC    1 15.0
## 3   VC    2 33.0
## 4   OJ  0.5  4.2
## 5   OJ    1 10.0
## 6   OJ    2 29.5

Now lets plot where both axes are treated as continuous labels

df3$dose <- as.numeric(as.vector(df3$dose))
ggplot(data = df3, aes(x = dose, y = len, group = supp, color = supp)) +
  geom_line() + geom_point()

Now lets look at a line graph with having the x axis as dates. We’ll use the built in economics time series for this example.

head(economics)

## # A tibble: 6 × 6
##   date         pce    pop psavert uempmed unemploy
##   <date>     <dbl>  <dbl>   <dbl>   <dbl>    <dbl>
## 1 1967-07-01  507. 198712    12.6     4.5     2944
## 2 1967-08-01  510. 198911    12.6     4.7     2945
## 3 1967-09-01  516. 199113    11.9     4.6     2958
## 4 1967-10-01  512. 199311    12.9     4.9     3143
## 5 1967-11-01  517. 199498    12.8     4.7     3066
## 6 1967-12-01  525. 199657    11.8     4.8     3018

ggplot(data = economics, aes(x = date, y = pop)) +
  geom_line()

Now lets subset the data

ss <- subset(economics, date > as.Date("2006-1-1"))
ggplot(data = ss, aes(x = date, y = pop)) +
  geom_line()

We can also change the line size, for instance, by another variable like unemployment

ggplot(data = economics, aes(x = date, y = pop)) +
  geom_line(aes(size = unemploy / pop))

We can also plot multiple time-series data

ggplot(economics, aes(x = date)) +
  geom_line(aes(y = psavert), color = "darkred") +
  geom_line(aes(y = uempmed), color = "steelblue", linetype = "twodash")

Lastly, lets make this into an are plot

ggplot(economics, aes(x = date)) +
  geom_area(aes(y = psavert), fill = "#999999",
            color = "#999999", alpha = 0.5) +
  geom_area(aes(y = uempmed), fill = "#E69F00",
            color = "#E69F00", alpha = 0.5)

Ridge Plots

First lets load the required packages

library(ggplot2)
library(ggridges)

## Warning: package 'ggridges' was built under R version 4.4.3

#BiocManager::install("ggridges")

Now lets load some sample data

airquality

##     Ozone Solar.R Wind Temp Month Day
## 1      41     190  7.4   67     5   1
## 2      36     118  8.0   72     5   2
## 3      12     149 12.6   74     5   3
## 4      18     313 11.5   62     5   4
## 5      NA      NA 14.3   56     5   5
## 6      28      NA 14.9   66     5   6
## 7      23     299  8.6   65     5   7
## 8      19      99 13.8   59     5   8
## 9       8      19 20.1   61     5   9
## 10     NA     194  8.6   69     5  10
## 11      7      NA  6.9   74     5  11
## 12     16     256  9.7   69     5  12
## 13     11     290  9.2   66     5  13
## 14     14     274 10.9   68     5  14
## 15     18      65 13.2   58     5  15
## 16     14     334 11.5   64     5  16
## 17     34     307 12.0   66     5  17
## 18      6      78 18.4   57     5  18
## 19     30     322 11.5   68     5  19
## 20     11      44  9.7   62     5  20
## 21      1       8  9.7   59     5  21
## 22     11     320 16.6   73     5  22
## 23      4      25  9.7   61     5  23
## 24     32      92 12.0   61     5  24
## 25     NA      66 16.6   57     5  25
## 26     NA     266 14.9   58     5  26
## 27     NA      NA  8.0   57     5  27
## 28     23      13 12.0   67     5  28
## 29     45     252 14.9   81     5  29
## 30    115     223  5.7   79     5  30
## 31     37     279  7.4   76     5  31
## 32     NA     286  8.6   78     6   1
## 33     NA     287  9.7   74     6   2
## 34     NA     242 16.1   67     6   3
## 35     NA     186  9.2   84     6   4
## 36     NA     220  8.6   85     6   5
## 37     NA     264 14.3   79     6   6
## 38     29     127  9.7   82     6   7
## 39     NA     273  6.9   87     6   8
## 40     71     291 13.8   90     6   9
## 41     39     323 11.5   87     6  10
## 42     NA     259 10.9   93     6  11
## 43     NA     250  9.2   92     6  12
## 44     23     148  8.0   82     6  13
## 45     NA     332 13.8   80     6  14
## 46     NA     322 11.5   79     6  15
## 47     21     191 14.9   77     6  16
## 48     37     284 20.7   72     6  17
## 49     20      37  9.2   65     6  18
## 50     12     120 11.5   73     6  19
## 51     13     137 10.3   76     6  20
## 52     NA     150  6.3   77     6  21
## 53     NA      59  1.7   76     6  22
## 54     NA      91  4.6   76     6  23
## 55     NA     250  6.3   76     6  24
## 56     NA     135  8.0   75     6  25
## 57     NA     127  8.0   78     6  26
## 58     NA      47 10.3   73     6  27
## 59     NA      98 11.5   80     6  28
## 60     NA      31 14.9   77     6  29
## 61     NA     138  8.0   83     6  30
## 62    135     269  4.1   84     7   1
## 63     49     248  9.2   85     7   2
## 64     32     236  9.2   81     7   3
## 65     NA     101 10.9   84     7   4
## 66     64     175  4.6   83     7   5
## 67     40     314 10.9   83     7   6
## 68     77     276  5.1   88     7   7
## 69     97     267  6.3   92     7   8
## 70     97     272  5.7   92     7   9
## 71     85     175  7.4   89     7  10
## 72     NA     139  8.6   82     7  11
## 73     10     264 14.3   73     7  12
## 74     27     175 14.9   81     7  13
## 75     NA     291 14.9   91     7  14
## 76      7      48 14.3   80     7  15
## 77     48     260  6.9   81     7  16
## 78     35     274 10.3   82     7  17
## 79     61     285  6.3   84     7  18
## 80     79     187  5.1   87     7  19
## 81     63     220 11.5   85     7  20
## 82     16       7  6.9   74     7  21
## 83     NA     258  9.7   81     7  22
## 84     NA     295 11.5   82     7  23
## 85     80     294  8.6   86     7  24
## 86    108     223  8.0   85     7  25
## 87     20      81  8.6   82     7  26
## 88     52      82 12.0   86     7  27
## 89     82     213  7.4   88     7  28
## 90     50     275  7.4   86     7  29
## 91     64     253  7.4   83     7  30
## 92     59     254  9.2   81     7  31
## 93     39      83  6.9   81     8   1
## 94      9      24 13.8   81     8   2
## 95     16      77  7.4   82     8   3
## 96     78      NA  6.9   86     8   4
## 97     35      NA  7.4   85     8   5
## 98     66      NA  4.6   87     8   6
## 99    122     255  4.0   89     8   7
## 100    89     229 10.3   90     8   8
## 101   110     207  8.0   90     8   9
## 102    NA     222  8.6   92     8  10
## 103    NA     137 11.5   86     8  11
## 104    44     192 11.5   86     8  12
## 105    28     273 11.5   82     8  13
## 106    65     157  9.7   80     8  14
## 107    NA      64 11.5   79     8  15
## 108    22      71 10.3   77     8  16
## 109    59      51  6.3   79     8  17
## 110    23     115  7.4   76     8  18
## 111    31     244 10.9   78     8  19
## 112    44     190 10.3   78     8  20
## 113    21     259 15.5   77     8  21
## 114     9      36 14.3   72     8  22
## 115    NA     255 12.6   75     8  23
## 116    45     212  9.7   79     8  24
## 117   168     238  3.4   81     8  25
## 118    73     215  8.0   86     8  26
## 119    NA     153  5.7   88     8  27
## 120    76     203  9.7   97     8  28
## 121   118     225  2.3   94     8  29
## 122    84     237  6.3   96     8  30
## 123    85     188  6.3   94     8  31
## 124    96     167  6.9   91     9   1
## 125    78     197  5.1   92     9   2
## 126    73     183  2.8   93     9   3
## 127    91     189  4.6   93     9   4
## 128    47      95  7.4   87     9   5
## 129    32      92 15.5   84     9   6
## 130    20     252 10.9   80     9   7
## 131    23     220 10.3   78     9   8
## 132    21     230 10.9   75     9   9
## 133    24     259  9.7   73     9  10
## 134    44     236 14.9   81     9  11
## 135    21     259 15.5   76     9  12
## 136    28     238  6.3   77     9  13
## 137     9      24 10.9   71     9  14
## 138    13     112 11.5   71     9  15
## 139    46     237  6.9   78     9  16
## 140    18     224 13.8   67     9  17
## 141    13      27 10.3   76     9  18
## 142    24     238 10.3   68     9  19
## 143    16     201  8.0   82     9  20
## 144    13     238 12.6   64     9  21
## 145    23      14  9.2   71     9  22
## 146    36     139 10.3   81     9  23
## 147     7      49 10.3   69     9  24
## 148    14      20 16.6   63     9  25
## 149    30     193  6.9   70     9  26
## 150    NA     145 13.2   77     9  27
## 151    14     191 14.3   75     9  28
## 152    18     131  8.0   76     9  29
## 153    20     223 11.5   68     9  30

air <- ggplot(airquality) + aes(Temp, Month, group = Month) + geom_density_ridges()

air

## Picking joint bandwidth of 2.65

Now lets add some pizzaz to our graph

library(viridis)

## Warning: package 'viridis' was built under R version 4.4.3

## Loading required package: viridisLite

ggplot(airquality) + aes(Temp, Month, group = Month, fill = ..x..) +
  geom_density_ridges_gradient() + 
  scale_fill_viridis(option = "C", name = "Temp")

## Warning: The dot-dot notation (`..x..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(x)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Picking joint bandwidth of 2.65

Last thing we will do is create a facet plot for all our data

library(tidyr)

## Warning: package 'tidyr' was built under R version 4.4.3

airquality %>%
gather(key = "Measurement", value = "value", Ozone, Solar.R, Wind, Temp) %>%
  ggplot() + aes(value, Month, group = Month) +
  geom_density_ridges() +
  facet_wrap(~ Measurement, scales = "free")

## Picking joint bandwidth of 11

## Picking joint bandwidth of 40.1

## Picking joint bandwidth of 2.65

## Picking joint bandwidth of 1.44

## Warning: Removed 44 rows containing non-finite outside the scale range
## (`stat_density_ridges()`).

Density Plots

A density plot is a nice alternative to a histogram

set.seed(1234)

wdata = data.frame(
  sex = factor(rep(c("F", "M"), each = 200)),
  weight = c(rnorm(200, 55), rnorm(200, 58))
)

library(dplyr)
mu <- wdata %>%
  group_by(sex) %>%
  summarise(grp.mean = mean(weight))

Now lets load the graphing packages

library(ggplot2)
theme_set(
  theme_classic() +
    theme(legend.position = "right")
)

Now lets do the basic plot function. First we will create a ggplot object

d <- ggplot(wdata, aes(x = weight))

Now lets do a basic density plot

d + geom_density() +
  geom_vline(aes(xintercept = mean(weight)), linetype = "dashed")

Now lets change the y axis to count instead of density

d + geom_density(aes(y = stat(count)), fill = "lightgray") +
  geom_vline(aes(xintercept = mean(weight)), linetype = "dashed")

d + geom_density(aes(color = sex)) +
  scale_color_manual(values = c("pink", "lightblue"))

Lastly, lets fill the density plots

d + geom_density(aes(fill = sex), alpha = 0.4) +
  geom_vline(aes(xintercept = grp.mean, color = sex), data = mu, linetype = "dashed") +
  scale_color_manual(values = c("indianred", "lightblue")) +
  scale_fill_manual(values = c("indianred", "lightblue"))

Plotly

Line Plots

First lets load our required packagef

library(plotly)

## Warning: package 'plotly' was built under R version 4.4.3

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

Orange <- as.data.frame(Orange)

plot_ly(data = Orange, x = ~age, y = ~circumference)

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Now lets add some more infor

plot_ly(data = Orange, x = ~age, y = ~circumference,
        color = ~Tree, size = ~age,
        text = ~paste("Tree ID:", Tree, "<br>Age:", age, "Circ:", circumference))

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.

Now lets create a random distribution and add it to our dataframe

trace_1 <- rnorm(35, mean = 120, sd = 10)
new_data <- data.frame(Orange, trace_1)

We’ll use the random numbers as lines on the graph

plot_ly(data = new_data, x = ~age, y = ~circumference, color = ~Tree, size = ~age, text = ~paste("Tree ID:", Tree, "<br>Age:", age, "Circ:", circumference)) %>%
  add_trace(y = ~trace_1, mode = "lines") %>%
  add_trace(y = ~circumference, mode = "markers")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.

Now lets create a graph with the optioin of showing as a scatter or line, and add labels

plot_ly(data = Orange, x = ~age, y = ~circumference,
        color = ~Tree, size = ~circumference,
        text = ~paste("Tree ID:", Tree, "<br>Age:", age, "Circ:", circumference)) %>%
  add_trace(y = ~circumference, mode = "markers") %>%
  layout(
    title = "Plot of Orange Data with Switchable Trace",
    updatemenus = list(
      list(
        type = "dropdown",
        y = 0.8,
        buttons = list(
          list(method = "restyle",
              args = list("mode", "markers"),
              label = "Marker"),
          list(method = "restyle",
              args = list("mode", "lines"),
              label = "Lines")
        )
      )
    )
  )

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.
## Warning: `line.width` does not currently support multiple values.

Plotly 3D

First lets load our required packages

library(plotly)

Now lets create a random 3d matrix

d <- data.frame(
  x <- seq(1,10, by = 0.5),
  y <- seq(1,10, by = 0.5)
)

z <- matrix(rnorm(length(d$x) * length(d$y)), nrow = length(d$x), ncol = length(d$y))

Now lets plot our 3D data

plot_ly(d, x = ~x, y = ~y, z = ~z) %>%
  add_surface()

Lets add some more aspects to it, such as a topography

plot_ly(d, x = ~x, y = ~y, z = ~z) %>%
  add_surface(
    contours = list(
      z = list(
        show = TRUE,
        usecolormap = TRUE,
        highlightcolor = "#FF0000",
        project = list(z = TRUE)
      )
    )
  )

Now lets look at a 3D scatter plot

plot_ly(longley, x = ~GNP, y = ~Population, z = ~Employed, marker = list(color = ~GNP)) %>%
  add_markers()

Other Graphing Techniques

Error Bars

First lets load our required libraries

library(ggplot2)
library(dplyr)
library(plotrix)

## Warning: package 'plotrix' was built under R version 4.4.3

theme_set(
  theme_classic() +
    theme(legend.position = 'top')
)

Lets again use the tooth data for this exercise

df <- ToothGrowth
df$dose <- as.factor(df$dose)

Now lets use dplyr for manipulation purposes

df.summary <- df %>%
  group_by(dose) %>%
  summarise(
    sd = sd(len, na.rm = TRUE),
    stderr = std.error(len, na.rm = TRUE),
    len = mean(len)
  )

df.summary

## # A tibble: 3 × 4
##   dose     sd stderr   len
##   <fct> <dbl>  <dbl> <dbl>
## 1 0.5    4.50  1.01   10.6
## 2 1      4.42  0.987  19.7
## 3 2      3.77  0.844  26.1

Lets now look at some key functions

• geom_crossbar() for hollow bars with middle indicated by a horizontal line
• geom_errorbar() for error bars
• geom_errorbarh() for horizontal error bars
• geom_linerange() for drawing an interval represented by a vertical line
• geom_pointrange() for creating an interval represented by a vertical line, with a point in the middle

Lets start by creating a ggplot object

tg <- ggplot(
  df.summary,
  aes(x =dose, y = len, ymin = len - sd, ymax = len + sd)
)

Now lets look at the most basic error bars

tg + geom_pointrange()

tg + geom_errorbar(width = 0.2) +
  geom_point(size = 1.5)

Now lets create horizontal error bars by manipulating our graph

ggplot(df.summary, aes(x = len, y = dose, xmin = len - sd, xmax = len + sd)) +
  geom_point() +
  geom_errorbarh(height = 0.2)

## Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `height` was translated to `width`.

This just gives you an idea of error bars on the horizontal axis

Now lets look at adding jitter points (actual measurements) to our data

ggplot(df, aes(dose, len)) +
  geom_jitter(position = position_jitter(0.2), color = "darkgray") +
  geom_pointrange(aes(ymin = len - sd, ymax = len + sd), data = df.summary)

Now lets try error bars oon a violin plot

ggplot(df, aes(dose, len)) +
  geom_violin(color = "darkgrey", trim = FALSE) +
  geom_pointrange(aes(ymin = len - sd, ymax = len + sd), data = df.summary)

Now how about with a line graph?

ggplot(df.summary, aes(dose, len)) +
  geom_line(aes(group = 1)) + # Always specify this when you have 1 line
  geom_errorbar(aes(ymin = len - stderr, ymax = len + stderr), width = 0.2) +
  geom_point(size = 2)

Now lets make a bar graph with half error bars

ggplot(df.summary, aes(dose, len)) +
  geom_col(fill = "lightgray", color = "black") +
  geom_errorbar(aes(ymin = len, ymax = len + stderr), width = 0.2)

You can see that by not specifying ymin = len - stderr, we have in essence cut our error bar in half.

How aboutwe add jitter points to line plots? We need to usse the original dataframe for the jitter plot, and the summary df for the geom layers.

ggplot(df, aes(dose, len)) +
  geom_jitter(position = position_jitter(0.2), color = "darkgray") + 
  geom_line(aes(group = 1), data = df.summary) +
  geom_errorbar(
    aes(ymin = len - stderr, ymax = len + stderr),
    data = df.summary, width = 0.2) +
      geom_point(data = df.summary, size = 0.2)

What about adding jitterpoints to a barplot?

ggplot(df, aes(dose, len)) +
  geom_col(data = df.summary, fill = NA, color = "black") +
  geom_jitter(position = position_jitter(0.2), color = "blue") +
  geom_errorbar(aes(ymin = len - stderr, ymax = len + stderr),
                data = df.summary, width = 0.2)

What if we wanted to have our error bars per group (OJ vs VC)

df.summary2 <- df %>%
  group_by(dose, supp) %>%
  summarise(
    sd = sd(len),
    stderr = std.error(len),
    len = mean(len)
  )

## `summarise()` has grouped output by 'dose'. You can override using the
## `.groups` argument.

df.summary2

## # A tibble: 6 × 5
## # Groups:   dose [3]
##   dose  supp     sd stderr   len
##   <fct> <fct> <dbl>  <dbl> <dbl>
## 1 0.5   OJ     4.46  1.41  13.2 
## 2 0.5   VC     2.75  0.869  7.98
## 3 1     OJ     3.91  1.24  22.7 
## 4 1     VC     2.52  0.795 16.8 
## 5 2     OJ     2.66  0.840 26.1 
## 6 2     VC     4.80  1.52  26.1

Now you can see we have mean and error for each dose and supp

 ggplot(df.summary2, aes(dose, len)) +
  geom_pointrange(
    aes(ymin = len - stderr, ymax = len + stderr, color = supp),
    position = position_dodge(0.3)) +
  scale_color_manual(values = c("indianred", "lightblue"))

How about line plots with multiple error bars

ggplot(df.summary2, aes(dose, len)) +
  geom_line(aes(linetype = supp, group = supp)) +
  geom_point() +
  geom_errorbar(aes(ymin = len - stderr, ymax = len + stderr, group = supp), width = 0.2)

And the same with a bar plot

ggplot(df.summary2, aes(dose, len)) +
  geom_col(aes(fill = supp), position = position_dodge(0.8), width = 0.7) +
  geom_errorbar(
    aes(ymin = len - stderr, ymax = len + stderr, group = supp),
    width = 0.2, position = position_dodge(0.8)) +
  scale_fill_manual(values = c("indianred", "lightblue"))

Now lets add some jitterpoints

ggplot(df, aes(dose, len, color = supp)) +
  geom_jitter(position = position_dodge(0.2)) +
  geom_line(aes(group = supp), data = df.summary2) +
  geom_errorbar(aes(ymin = len - stderr, ymax = len + stderr, group = supp), data = df.summary2, width = 0.2)

ggplot(df, aes(dose, len, color = supp)) +
  geom_col(data = df.summary2, position = position_dodge(0.8), width = 0.7, fill = "white") +
  geom_jitter(
    position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.8)) +
  geom_errorbar(
    aes(ymin = len - stderr, ymax = len + stderr), data = df.summary2,
        width = 0.2, position = position_dodge(0.8)) +
      scale_color_manual(values = c("indianred", "lightblue")) +
      theme(legend.position = "top")

ECDF Plots

Now lets do an empirical cumulative distribution function. This reports any given number percentile of individuals that are above or below that threshold.

set.seed(1234)

wdata = data.frame(
  sex = factor(rep(c("F", "M"), each = 200)),
  weight = c(rnorm(200, 55), rnorm(200, 58))
)

Now lets look at our dataframe

head(wdata, 5)

##   sex   weight
## 1   F 53.79293
## 2   F 55.27743
## 3   F 56.08444
## 4   F 52.65430
## 5   F 55.42912

Now lets load our plotting package

library(ggplot2)

theme_set(
  theme_classic() +
    theme(legend.position = "bottom")
)

Now lets create our ECDF Plot

ggplot(wdata, aes(x = weight)) + 
  stat_ecdf(aes(color = sex, linetype = sex),
            geom = "step", size = 1.5) +
  scale_color_manual(values = c("indianred", "lightblue")) +
  labs(y = "weight")

qq Plots

Now lets take a look at qq plots. These are used to determine if the given data follows a normal distribution.

set.seed(1234)

Now lets randomly generate some data

wdata = data.frame(
  sex = factor(rep(c("F", "M"), each = 200)),
  weight = c(rnorm(200, 55), rnorm(200, 58))
)

Lets set our theme for the graphing with ggplot

library(ggplot2)

theme_set(
  theme_classic() +
    theme(legend.position = "top")
)

create a qq plot of the weight

ggplot(wdata, aes(sample = weight)) +
  stat_qq(aes(color = sex)) +
  scale_color_manual(values = c("indianred", "lightblue")) +
  labs(y = "weight")

library(ggpubr)

## Warning: package 'ggpubr' was built under R version 4.4.3

ggqqplot(wdata, x = "weight",
         color = "sex",
         palettes = c("indianred","lightblue"),
         ggtheme = theme_pubclean())

Now what a non-normal distribution look like?

library(mnonr)

## Warning: package 'mnonr' was built under R version 4.4.3

data2 <- mnonr::mnonr(n = 1000, p = 2, ms = 3, mk = 61, Sigma = matrix(c(1, 0.5, 0.5, 1), 2, 2), initial = NULL)

data2 <- as.data.frame(data2)

Now lets plot the non normal data

ggplot(data2, aes(sample = V1)) +
  stat_qq()

ggqqplot(data2, x = "V1",
         pallets = "indianred",
         ggtheme = theme_pubclean())

Facet Plots

Lets look at how to put multiple plots together into a single figure

library(ggpubr)
library(ggplot2)

theme_set(
  theme_bw() +
    theme(legend.position = "top")
)

First lets create a nice boxplot

Lets load the data

df <- ToothGrowth
df$dose <- as.factor(df$dose)

and create the plot object

p <- ggplot(df, aes(x = dose, y = len)) +
  geom_boxplot(aes(fill = supp), position = position_dodge(0.9)) +
  scale_fill_manual(values = c("#00AFBB", "#E7B800"))

p

Now lets look at the ggplot facet function

p + facet_grid(rows = vars(supp))

Now lets do a facet with multiple variables

p + facet_grid(rows = vars(dose), cols = vars(supp))

p

Now lets look at the facet_wrap function. This allows facets to be placed side-by-side

p + facet_wrap(vars(dose), ncol = 2)

Now how do we combine multiple plots using ggarragnge()

Lets start by making some basic plots. First we will define a color palette and data

my3cols <- c("#e7b800", "#2e9fdf", "#fc4e07")
ToothGrowth$dose <- as.factor(ToothGrowth$dose)

Now lets make some basic plots

p <- ggplot(ToothGrowth, aes(x = dose, y = len))
bxp <- p + geom_boxplot(aes(color = dose)) +
  scale_color_manual(values = my3cols)

Ok, now lets do a dotplot

dp <- p + geom_dotplot(aes(color = dose, fill = dose),
                       binaxis = "y", stackdir = "center") +
  scale_color_manual(values = my3cols) +
  scale_fill_manual(values = my3cols)

Now lastly lets create a lineplot

lp <- ggplot(economics, aes(x = date, y = psavert)) +
  geom_line(color = "indianred")

Now we can make the figure

figure <- ggarrange(bxp, dp, lp, labels = c("A", "B", "C"), ncol = 2, nrow = 2)

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

figure

This looks great, but we can make it look even better

figure2 <- ggarrange(
  lp,
  ggarrange(bxp, dp, ncol = 2, labels = c("B", "C")),
  nrow = 2,
  labels = "A"
)

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

figure2

This looks really good, but you’ll notice that there are two legends that are the same.

ggarrange(
  bxp, dp, labels = c("A", "B"),
  common.legend = TRUE, legend = "bottom"
)

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

Lastly, we should export the plot

ggexport(figure2, filename = "facetfigure.pdf")

## file saved to facetfigure.pdf

We can also export multiple plots to a pdf

ggexport(bxp, dp, lp, filename = "multi.pdf")

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## file saved to multi.pdf

Lastly, we can export to pdf with multiple pages and multiple columns

ggexport(bxp, dp, lp, bxp, filename = "test2.pdf", nrow = 2, ncol = 1)

## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.
## file saved to test2.pdf

Heatmaps

Lets get started with heatmaps

library(heatmap3)

## Warning: package 'heatmap3' was built under R version 4.4.3

Now lets get our data

data <- ldeaths

data2 <- do.call(cbind, split(data, cycle(data)))
dimnames(data2) <- dimnames(.preformat.ts(data))

Now lets generate a heat map

heatmap(data2)

heatmap(data2, Rowv = NA, Colv = NA)

Now lets play with the colors

rc <- rainbow(nrow(data2), start = 0, end = 0.3)
cc <- rainbow(ncol(data2), start = 0, end = 0.3)

Now lets apply our colr selections

heatmap(data2, ColSideColors = cc)

library(RColorBrewer)

heatmap(data2, ColSideColors = cc,
        col = colorRampPalette(brewer.pal(8, "PiYG"))(25))

Theres more that we can customize

library(gplots)

## Warning: package 'gplots' was built under R version 4.4.3

## 
## ---------------------
## gplots 3.3.0 loaded:
##   * Use citation('gplots') for citation info.
##   * Homepage: https://talgalili.github.io/gplots/
##   * Report issues: https://github.com/talgalili/gplots/issues
##   * Ask questions: https://stackoverflow.com/questions/tagged/gplots
##   * Suppress this message with: suppressPackageStartupMessages(library(gplots))
## ---------------------

## 
## Attaching package: 'gplots'

## The following object is masked from 'package:plotrix':
## 
##     plotCI

## The following object is masked from 'package:stats':
## 
##     lowess

heatmap.2(data2, ColSideColors = cc,
          col = colorRampPalette(brewer.pal(8, "PiYG"))(25))

Outlier Detection

Missing Values

Missing Values

If you encounter an unusual value in your dataset, and simply want to move on to the rest of your analysis, you have two options:

Drop the entire row with the strange values:

library(dplyr)
library(ggplot2)

diamonds <- diamonds

diamonds2 <- diamonds %>%
  filter(between(y, 3, 20))

In this instance, y is the width of the diamond, so anything under 3mm or above 20 is excluded

I don’t recommend this option, just because there is one bad measurements doesn’t mean they are all bad

Instead, I recommend replacing the unusual values with missing values

diamonds3 <- diamonds %>%
  mutate(y = ifelse(y < 3 | y > 20, NA, y))

Like R, ggplot2 subscribes to the idea that missing values shouldn’t pass silently into the night.

ggplot(data = diamonds3, mapping = aes(x = x, y = y)) +
  geom_point()

## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).

If you want to suppress that warning you can use na.rm = TRUE

ggplot(data = diamonds3, mapping = aes(x = x, y=y)) +
  geom_point(na.rm = TRUE)

Other times you want to understand what makes observations with missing values different to the observation with recorded values. For example, in the NYCflights13 dataset, missing values in the dep_time variable indicate that the flight was canceled, So you might want to compare the scheduled departure times for canceled and non-canceled times.

library(nycflights13)

## Warning: package 'nycflights13' was built under R version 4.4.3

nycflights13::flights %>%
  mutate(
    canceled = is.na(dep_time),
    sched_hour = sched_dep_time %/% 100,
    sched_min = sched_dep_time %% 100,
    sched_dept_time = sched_hour + sched_min / 60
  ) %>%
  ggplot(mapping = aes(sched_dept_time)) +
  geom_freqpoly(mapping = aes(color = canceled), binwidth = 1/4)

Outliers

See the end of the Exploratory Data chapter

Covariation

CATEGORICAL VARIABLES

library(ggplot2)

ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500)

Its hard to see the difference in distribution because the counts differ so much.

ggplot(diamonds) +
  geom_bar(mapping = aes(x = cut))

To make the comparison easier, we need to swap the display on the y-axis. Instead of displaying count, we’ll display density, which is the count standardized so that the area under the curve is one

ggplot(data = diamonds, mapping = aes(x = price, y = ..density..)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500)

It appears that fair diamonds (the lowest cut quality) have the highest average price, but maybe that’s because frequency polygons are a little hard to interpret.

Another alternative is the boxplot. A boxplot is a type of visual shorthand for a distribution of values.

ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
  geom_boxplot()

We see much less iinformation about the distribution, but the boxplots are much more compact, so we can more easily compare them. It supports the counterintuitive finding that better quality diamonds are cheaper on average!

Lets look at some car data

ggplot(data = mpg, mapping = aes(x = class, y = hwy)) +
  geom_boxplot()

ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median), y = hwy))

If you have long variable names, you can switch the axis and flip it 90 degrees

ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median),y = hwy)) +
  coord_flip()

To visualize the correlation between two continous variables, we can use a scatter plot.

ggplot(data = diamonds) +
  geom_point(mapping = aes(x = carat, y = price))

Scatterplots become less useful as the size of your dataset grows, because we get overplot. We can fix this by using the alpha aesthetic

ggplot(data = diamonds) +
  geom_point(mapping = aes(x = carat, y = price), alpha = 1/100)

Exploratory Data Analysis

Exploratory Data Analysis

First, lets load a required library

library(RCurl)

## Warning: package 'RCurl' was built under R version 4.4.3

## 
## Attaching package: 'RCurl'

## The following object is masked from 'package:tidyr':
## 
##     complete

library(dplyr)

Now lets get our data

site <- "https://raw.githubusercontent.com/nytimes/covid-19-data/master/colleges/colleges.csv"

College_Data <- read.csv(site)

First lets use the str function, this shows the structure of the object

str(College_Data)

## 'data.frame':    1948 obs. of  9 variables:
##  $ date      : chr  "2021-05-26" "2021-05-26" "2021-05-26" "2021-05-26" ...
##  $ state     : chr  "Alabama" "Alabama" "Alabama" "Alabama" ...
##  $ county    : chr  "Madison" "Montgomery" "Limestone" "Lee" ...
##  $ city      : chr  "Huntsville" "Montgomery" "Athens" "Auburn" ...
##  $ ipeds_id  : chr  "100654" "100724" "100812" "100858" ...
##  $ college   : chr  "Alabama A&M University" "Alabama State University" "Athens State University" "Auburn University" ...
##  $ cases     : int  41 2 45 2742 220 4 263 137 49 76 ...
##  $ cases_2021: int  NA NA 10 567 80 NA 49 53 10 35 ...
##  $ notes     : chr  "" "" "" "" ...

What if we want to arrange our data set alphabetically by college?

alphabetical <- College_Data %>%
  arrange(College_Data$college)

The glimpse function is another way to preview data

glimpse(College_Data)

## Rows: 1,948
## Columns: 9
## $ date       <chr> "2021-05-26", "2021-05-26", "2021-05-26", "2021-05-26", "20…
## $ state      <chr> "Alabama", "Alabama", "Alabama", "Alabama", "Alabama", "Ala…
## $ county     <chr> "Madison", "Montgomery", "Limestone", "Lee", "Montgomery", …
## $ city       <chr> "Huntsville", "Montgomery", "Athens", "Auburn", "Montgomery…
## $ ipeds_id   <chr> "100654", "100724", "100812", "100858", "100830", "102429",…
## $ college    <chr> "Alabama A&M University", "Alabama State University", "Athe…
## $ cases      <int> 41, 2, 45, 2742, 220, 4, 263, 137, 49, 76, 67, 0, 229, 19, …
## $ cases_2021 <int> NA, NA, 10, 567, 80, NA, 49, 53, 10, 35, 5, NA, 10, NA, 19,…
## $ notes      <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",…

We can also subset with select()

College_Cases <- select(College_Data, college, cases)

We can also filter or subset with the filter function

Louisiana_Cases <- filter(College_Data, state == "Louisiana")

Lets filter out a smaller amount of states

South_Cases <- filter(College_Data, state == "Louisiana" | state == "Texas" | state == "Arkansas" | state == "Mississippi")

Lets look at some time series data

First we’ll load the required libraries

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(dplyr)
library(ggplot2)
library(gridExtra)

## Warning: package 'gridExtra' was built under R version 4.4.3

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

library(scales)

## Warning: package 'scales' was built under R version 4.4.3

## 
## Attaching package: 'scales'

## The following object is masked from 'package:plotrix':
## 
##     rescale

## The following object is masked from 'package:viridis':
## 
##     viridis_pal

Now lets load some data

state_site <- "https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-states.csv"

State_Data <- read.csv(state_site)

Lets create group_by object using the state column

state_cases <- group_by(State_Data, state)

class(state_cases)

## [1] "grouped_df" "tbl_df"     "tbl"        "data.frame"

How many measurements were made by state? This gives us an idea when states started reporting

Days_since_first_reported <- tally(state_cases)

Lets visualize some data

First lets start off with some definitions

Data - obvious = the stuff we want to visualize

Layer - made of geometric elements and requisite statistical information. Include geometric objects which represent the plot

Scales - used to map values in the data space that is used for creation of values (color, size shape, etc)

Coordinate system - describes how the data coordinates are mapped together in relation to the plan on the graphic

Faceting - How to break up data into subsets to display multiple types or groups of data

Theme - controls the finer points of the display, such as font size and background color

options(repr.plot.width = 6, repr.plot.height = 6)

class(College_Data)

## [1] "data.frame"

head(College_Data)

##         date   state     county       city ipeds_id
## 1 2021-05-26 Alabama    Madison Huntsville   100654
## 2 2021-05-26 Alabama Montgomery Montgomery   100724
## 3 2021-05-26 Alabama  Limestone     Athens   100812
## 4 2021-05-26 Alabama        Lee     Auburn   100858
## 5 2021-05-26 Alabama Montgomery Montgomery   100830
## 6 2021-05-26 Alabama     Walker     Jasper   102429
##                           college cases cases_2021 notes
## 1          Alabama A&M University    41         NA      
## 2        Alabama State University     2         NA      
## 3         Athens State University    45         10      
## 4               Auburn University  2742        567      
## 5 Auburn University at Montgomery   220         80      
## 6  Bevill State Community College     4         NA

summary(College_Data)

##      date              state              county              city          
##  Length:1948        Length:1948        Length:1948        Length:1948       
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##    ipeds_id           college              cases          cases_2021    
##  Length:1948        Length:1948        Min.   :   0.0   Min.   :   0.0  
##  Class :character   Class :character   1st Qu.:  32.0   1st Qu.:  23.0  
##  Mode  :character   Mode  :character   Median : 114.5   Median :  65.0  
##                                        Mean   : 363.5   Mean   : 168.1  
##                                        3rd Qu.: 303.0   3rd Qu.: 159.0  
##                                        Max.   :9914.0   Max.   :3158.0  
##                                                         NA's   :337     
##     notes          
##  Length:1948       
##  Class :character  
##  Mode  :character  
##                    
##                    
##                    
## 

Now lets take a look at a different dataset

iris <- as.data.frame(iris)

class(iris)

## [1] "data.frame"

head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

summary(iris)

##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 

Lets start by creating a scatter plot of the College Data

ggplot(data = College_Data, aes(x = cases, y = cases_2021)) +
  geom_point()+
  theme_minimal()

## Warning: Removed 337 rows containing missing values or values outside the scale range
## (`geom_point()`).

Now lets do the iris data

ggplot(data = iris, aes(x = Sepal.Width, y = Sepal.Length)) +
  geom_point() +
  theme_minimal()

Lets color coordinate our college data

ggplot(data = College_Data, aes(x = cases, y = cases_2021, color = state)) +
  geom_point() +
  theme_minimal()

## Warning: Removed 337 rows containing missing values or values outside the scale range
## (`geom_point()`).

Lets color coordinate the iris data

ggplot(data = iris, aes(x = Sepal.Width, y = Sepal.Length, color = Species)) +
  geom_point() +
  theme_minimal()

Lets run a simple histogram of our Louisiana Case Data

hist(Louisiana_Cases$cases, freq = NULL, density = NULL, breaks = 10, xlab = "Total Cases", ylab = "Frequency",
     main = "Total College Covid-19 Infections (Louisiana)")

Lets run a simple histogram for the Iris data

hist(iris$Sepal.Width, freq = NULL, density = NULL, breaks = 10, xlab = "Sepal Width",
     ylab = "Frequency", main = "Iris Sepal Width")

histogram_college <- ggplot(data = Louisiana_Cases, aes(x = cases))

histogram_college + geom_histogram(binwidth = 100, color = "black", aes(fill = county)) +
  xlab("cases") + ylab("Fequency") + ggtitle("Histogram of Covid-19 Cases in Louisiana")

Lets create a ggplot for the iris data

histogram_iris <- ggplot(data = iris, aes(x = Sepal.Width))

histogram_iris + geom_histogram(binwidth = 0.2, color = "black", aes(fill = Species)) +
  xlab("Sepal Width") + ylab("Frequency") + ggtitle("Histogram of Iris Sepal Width by Species")

Maybe a density plot makes more sense for our college data

ggplot(South_Cases) +
  geom_density(aes(x = cases, fill = state), alpha = 0.50)

Lets do it with the iris data

ggplot(iris) +
  geom_density(aes(x = Sepal.Width, fill = Species, alpha = 0.25))

Lets look at violin plots for iris

ggplot(data = iris, aes(x = Species, y = Sepal.Length, color = Species)) +
  geom_violin() +
  theme_classic() +
  theme(legend.position = "none")

Now lets try the south data

ggplot(data = South_Cases, aes(x = state, y = cases, color = state)) +
  geom_violin() +
  theme_gray() +
  theme(legend.position = "none")

Now lets take a look at risidual plots. This is a graph that displays the residuals on the certical axis, and the independent variable on the horizontal. In the event that the points in a residual plot are dispersed in a random manner around the horizontal axis, it is appropriate to use a linear regression. If they are not randomly dispersed, a non linear model is more appropriate.

Lets start with the iris data

ggplot(lm(Sepal.Length ~Sepal.Width, data = iris)) +
  geom_point(aes(x = .fitted, y = .resid))

## Warning: `fortify(<lm>)` was deprecated in ggplot2 4.0.0.
## ℹ Please use `broom::augment(<lm>)` instead.
## ℹ The deprecated feature was likely used in the ggplot2 package.
##   Please report the issue at <https://github.com/tidyverse/ggplot2/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Now look at the southern states cases

ggplot(lm(cases ~cases_2021, data = South_Cases)) +
  geom_point(aes(x = .fitted, y = .resid))

A linear model is not a good all for the state cases

Now lets do some correlations (Example data was not provided)

#obesity <- read.csv("Obesity_insurance.csv")

Lets look at the structure of the dataset

#str(obesity)

Lets look at the column classes

#class(obesity)

And get a summary of distribution of the variables

#summary(obesity)

Now lets look at the distribution for insurance charges

#hist(obesity$charges)

We can also get an idea of the distribution using a boxplot

#boxplot(obesity$charges)

#boxplot(obesity$bmi)

Now lets look at correlations. The cor() command is used to determine correlations between two vectors, all of the columns of a data frame, or two data frames. The cov() command, on the other hand, examines the covariance. The cor.test() command carries out a test as to the significance of the correlation

#cor(obesity$charges, obesity$bmi)

This test uses a spearman Rho correlation, or you can use Kendall’s tau by specifying it

#cor(obesity$charges, obesity$bmi, method = 'kendall')

This correlation measures strength of a correlation between -1 and 1

Now lets look at the Tietjen-Moore test. This is used for univariate datasets The algorithm depicts the detection of the outliers in a univariate dataset.

TietjenMoore <- function(dataSeries, k){
  n = length(dataSeries)
  # Compute the absolute residuals
  r = abs(dataSeries - mean(dataseries))
  # Sort data according to size of residual
  df = data.frame(dataSeries, r)
  dfs = df[order(df$r),]
  # create a subset of the data without the largest values
  klarge = c((n-k+1):n)
  subdataSeries = dfs$dataSeries[-klarge]
  # compute the sum of squares
  ksub = (subdataSeries = mean(subdataSeries)) ** 2
  all = (df$dataSeries - mean(df$dataSeries)) ** 2
  # compute the test statistic
  sum(ksub)/sum(all)
}

This function helps to compute the absolute residuals and sorts data according to the size of the residuals. Later, we will focus on the computation of sum of squares

FindOutliersTietjenMooreTest <- function(dataSeries, k, alpha = 0.5){
  ek <- TietjenMoore(dataSeries, k)
  # Compute critical values based on simulation.
  test = c(1:10000)
  for(i in 1:10000){
    dataSeriesdataSeries = rnorm(length(dataSeries))
    test[i] = TietgenMoore(dataSeriesdataSeries, k)
  }
  Talpha = quantile(test, alpha)
  list(T = ek, Talpha = Talpha)
}

This function helps us to compute the critical values based on simulation data. Now lets demonstrate these functions with sample data and the obesity dataset for evaluating this algorithm

The critical region for the Tietjen-Moore test is determined by simulation. The simulation is performed by generating a standard normal random sample of size n and compution the Tiethen Moore test statistic. Typically, 10,000 random samples are used. The values of the Tietjen-Moore statistic obtained from the data is compared to this reference distribution. The values of the test statistic is between zero and one. If there are no outliers in the data, the test statistic is close to 1. If there are outliers, the test statistic will be closer to zero. Thus, the test is always a lower, one-tailed test regardless of which test statistic issued, Lk or Ek

First we will look at charges

#boxplot(obesity$charges)

#FindOutliersTietjenMooreTest(obesity$charges, 4)

Lets check out bmi

#boxplot(obesity$bmi)

#FindOutliersTietjenMooreTest(obesity$bmi, 50)

Probability Plots

#library(tigerstats)
# Package is defunct

We will use the probability plot function and their output dnorm: density function of the normal distribution. Using the density, it is possible to determine the probability of events. Or for example, you may wonder “what is the likelihood that a person has a BMI of exactly __?” In this case, you would need to retrieve the density of the BMI distribution at values –. The BMI distribution can be modeled with a mean of 100 and a standard deviation of 15. The corresponding density is:

#bmi.mean <- mean(obesity$bmi)
#bmi.sd <- sd(obesity$bmi)

Lets create a plot of our normal distribution

#bmi.dist <- dnorm(obesity$bmi, mean = bmi.mean, sd = bmi.sd)
#bmi.df <- data.frame("bmi" = obesity$bmi, "Density" = bmi.dist)
#ggplot(bmi.df, aes(x = bmi, y = Density)) +
# geom_point()

This gives us the probability of every single point occuring

Now lets use the pnorm function for more info

#bmi.dist <- pnorm(obesity$bmi, mean = bmi.mean, sd = bmi.sd)
#bmi.df <- data.frame("bmi" = obesity$bmi, "Density" = bmi.dist)

#ggplot(bmi.df, aes(x = bmi, y = Density)) +
# geom_point()

What if we want to find the probability of the bmi being greater than 40 in our distribution?

pp_greater <- function(x){
  paste(round(100 * pnorm(x,, mean = 30.66339, sd = 6.09818, lower.tail =  FALSE), 2), "%")
  
  
}
pp_greater(40)

## [1] "6.29 %"

What is the probability that a bmi is less than 40 in our population?

pp_less <- function(x){
  paste(round(100 * (1 - pnorm(x, mean = 30.66339, sd = 6.09818, lower.tail = FALSE)), 2), "%")
}
pp_less(40)

## [1] "93.71 %"

#pnormGC(40, region = "below", mean = 30.66339, sd = 6.09818, graph = TRUE)

What if we want to find the area in between?

#pnormGC(c(20,40), region = "between", mean = 30.66339, sd = 6.09818, graph = TRUE )

What if we want to know the quantiles? Lets use the pnorm function. We need to assume a normal distribution for this.

What bmi represents the lower 1% of the population?

qnorm(0.01, mean = 30.66339, sd = 6.09818, lower.tail = TRUE)

## [1] 16.4769

What if you want a random sampling of values within your distribution?

subset <- rnorm(50, mean = 30.66339, sd = 6.09818)
hist(subset)

subset2 <- rnorm(5000, mean = 30.66339, sd = 6.09818)

hist(subset2)

Shapiro-Wilk Test

So now we know how to generate a normal distribution, how do we tell if our samples came from a normal distribution?

#shapiro.test(obesity$charges[1:5])

You can see here, with a small sample size, we would reject the null hypothesis that the samples came from a normal distribution. We can increase the power of the test by increasing the sample size

#shapiro.test(obesity$charges[1:1000])

Now lets check out age

#shapiro.test(obesity$age[1:1000])

and lastly BMI

#shapiro.test(obesity$bmi[1:1000])

Time Series Data

First lets load our packages

library(readr)

## 
## Attaching package: 'readr'

## The following object is masked from 'package:scales':
## 
##     col_factor

library(readxl)

Air_data <- read_xlsx("AirQualityUCI.xlsx")

Date - date of measurement Time - time of measurement CO(GT) - average hourly CO2 PT08,s1(CO) - tin oxide hourly average sensor response NMHC - average hourly non-metallic hydrocarbon concentration C6HC - average benzene concentration PT08.s3(NMHC) - titania average hourly sensor response NOx - average hourly NOx concentration NO2 - average hourly NO2 concentration T - temperature RH - relative humidity AH - absolute humidity

str(Air_data)

## tibble [9,357 × 15] (S3: tbl_df/tbl/data.frame)
##  $ Date         : POSIXct[1:9357], format: "2004-03-10" "2004-03-10" ...
##  $ Time         : POSIXct[1:9357], format: "1899-12-31 18:00:00" "1899-12-31 19:00:00" ...
##  $ CO(GT)       : num [1:9357] 2.6 2 2.2 2.2 1.6 1.2 1.2 1 0.9 0.6 ...
##  $ PT08.S1(CO)  : num [1:9357] 1360 1292 1402 1376 1272 ...
##  $ NMHC(GT)     : num [1:9357] 150 112 88 80 51 38 31 31 24 19 ...
##  $ C6H6(GT)     : num [1:9357] 11.88 9.4 9 9.23 6.52 ...
##  $ PT08.S2(NMHC): num [1:9357] 1046 955 939 948 836 ...
##  $ NOx(GT)      : num [1:9357] 166 103 131 172 131 89 62 62 45 -200 ...
##  $ PT08.S3(NOx) : num [1:9357] 1056 1174 1140 1092 1205 ...
##  $ NO2(GT)      : num [1:9357] 113 92 114 122 116 96 77 76 60 -200 ...
##  $ PT08.S4(NO2) : num [1:9357] 1692 1559 1554 1584 1490 ...
##  $ PT08.S5(O3)  : num [1:9357] 1268 972 1074 1203 1110 ...
##  $ T            : num [1:9357] 13.6 13.3 11.9 11 11.2 ...
##  $ RH           : num [1:9357] 48.9 47.7 54 60 59.6 ...
##  $ AH           : num [1:9357] 0.758 0.725 0.75 0.787 0.789 ...

library(lubridate)
library(hms)

## 
## Attaching package: 'hms'

## The following object is masked from 'package:lubridate':
## 
##     hms

Lets get rid of the date in the time column

Air_data$Time <- as_hms(Air_data$Time)
str(Air_data)

## tibble [9,357 × 15] (S3: tbl_df/tbl/data.frame)
##  $ Date         : POSIXct[1:9357], format: "2004-03-10" "2004-03-10" ...
##  $ Time         : 'hms' num [1:9357] 18:00:00 19:00:00 20:00:00 21:00:00 ...
##   ..- attr(*, "units")= chr "secs"
##  $ CO(GT)       : num [1:9357] 2.6 2 2.2 2.2 1.6 1.2 1.2 1 0.9 0.6 ...
##  $ PT08.S1(CO)  : num [1:9357] 1360 1292 1402 1376 1272 ...
##  $ NMHC(GT)     : num [1:9357] 150 112 88 80 51 38 31 31 24 19 ...
##  $ C6H6(GT)     : num [1:9357] 11.88 9.4 9 9.23 6.52 ...
##  $ PT08.S2(NMHC): num [1:9357] 1046 955 939 948 836 ...
##  $ NOx(GT)      : num [1:9357] 166 103 131 172 131 89 62 62 45 -200 ...
##  $ PT08.S3(NOx) : num [1:9357] 1056 1174 1140 1092 1205 ...
##  $ NO2(GT)      : num [1:9357] 113 92 114 122 116 96 77 76 60 -200 ...
##  $ PT08.S4(NO2) : num [1:9357] 1692 1559 1554 1584 1490 ...
##  $ PT08.S5(O3)  : num [1:9357] 1268 972 1074 1203 1110 ...
##  $ T            : num [1:9357] 13.6 13.3 11.9 11 11.2 ...
##  $ RH           : num [1:9357] 48.9 47.7 54 60 59.6 ...
##  $ AH           : num [1:9357] 0.758 0.725 0.75 0.787 0.789 ...

plot(Air_data$AH, Air_data$RH, main = "Humidity Analysis", xlab = "Absolute Humidty", ylab = "Relative Humidity")

Notice we have an outlier in our data

t.test(Air_data$RH, Air_data$AH)

## 
##  Welch Two Sample t-test
## 
## data:  Air_data$RH and Air_data$AH
## t = 69.62, df = 17471, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  45.01707 47.62536
## sample estimates:
## mean of x mean of y 
## 39.483611 -6.837604

What if we want to know what our outliers are?

First we need to load the required libraries

library(outliers)

And reload the dataset because we removed outliers (I have not removed the outliers, so I won’t)

Lets create a function using the grubb test to identify all outliers. The grubbs test identifies outliers in a univariate dataset that is presumed to come from a normal distribution.

grubbs.flag <- function(x) {
  # Lets create a variable called outliers and save nothing in it. We'll add to the variable
  # as we identify them
  outliers <- NULL
  # We'll create a variable called test to identify which univariate we are testing
  test <- x
  # Now using the outliers package, use grubbs.tset to find outliers in our variable
  grubbs.result <- grubbs.test(test)
  # Lets get the p-values of all tested variables
  pv <- grubbs.result$p.value
  # Now lets search through our p-values for ones that are outside of 0.5
  while(pv < 0.05){
    # anything with a pvalue greater than p = 0.05, we add to our empty outliers vector
    outliers <- c(outliers, as.numeric(strsplit(grubbs.result$alternative, " ")[[1]][3]))
    # Now we want to remove those outliers from our test variable
    test <- x[!x %in% outliers]
    # and run the grubbs test again without the outliers
    grubbs.result <- grubbs.test(test)
    # and save the new p values
    pv <- grubbs.result$p.value
  }
  return(data.frame(x = x, Outliers = (x %in% outliers)))
}

identified_outliers <- grubbs.flag(Air_data$AH)

Now we can create a histogram showing where the outliers were

ggplot(grubbs.flag(Air_data$AH), aes(x = Air_data$AH, color = Outliers, fill = Outliers)) +
  geom_histogram(binwidth = diff(range(Air_data$AH))/30) +
  theme_bw()

Text Mining

Text Mining

First we’ll look at the unnest_token function

Lets start by looking at an Emily Dickenson Passage

text <- c("Because I could not stop from Death -",
          "He kindly stopped for me -",
          "The carriage held but just ourselves -",
          "and immortality")

text

## [1] "Because I could not stop from Death -" 
## [2] "He kindly stopped for me -"            
## [3] "The carriage held but just ourselves -"
## [4] "and immortality"

This is a typical character vector that we might want to analyze. In order to turn it into a tidytext dataset, we first need to put it into a dataframe.

library(dplyr)

text_df <- tibble(line = 1:4, text = text)

text_df

## # A tibble: 4 × 2
##    line text                                  
##   <int> <chr>                                 
## 1     1 Because I could not stop from Death - 
## 2     2 He kindly stopped for me -            
## 3     3 The carriage held but just ourselves -
## 4     4 and immortality

Reminder: A tibble is a modern class of data frame within R. It’s available in the dplyr and tibble packages, that has a convenient print method, will not convert string to factors, and does not use row names. Tibbles are great for use with tidy tools

Next, we will use the “unest_tokens()” function

First we have the output column name that will be created as the text is unnested into it

library(tidytext)

## Warning: package 'tidytext' was built under R version 4.4.3

text_df %>%
  unnest_tokens(word, text)

## # A tibble: 20 × 2
##     line word       
##    <int> <chr>      
##  1     1 because    
##  2     1 i          
##  3     1 could      
##  4     1 not        
##  5     1 stop       
##  6     1 from       
##  7     1 death      
##  8     2 he         
##  9     2 kindly     
## 10     2 stopped    
## 11     2 for        
## 12     2 me         
## 13     3 the        
## 14     3 carriage   
## 15     3 held       
## 16     3 but        
## 17     3 just       
## 18     3 ourselves  
## 19     4 and        
## 20     4 immortality

Lets use the janeaustenr package to analyze some Jane Austen texts. There are 6 books in this package

library(janeaustenr)

## Warning: package 'janeaustenr' was built under R version 4.4.3

library(stringr)

original_books <- austen_books() %>%
  group_by(book) %>%
  mutate(linenumber = row_number(),
         chapter = cumsum(str_detect(text, regex("^chapter [\\divxlc]",
                                                 ignore_case = TRUE)))) %>%
  ungroup()

original_books

## # A tibble: 73,422 × 4
##    text                    book                linenumber chapter
##    <chr>                   <fct>                    <int>   <int>
##  1 "SENSE AND SENSIBILITY" Sense & Sensibility          1       0
##  2 ""                      Sense & Sensibility          2       0
##  3 "by Jane Austen"        Sense & Sensibility          3       0
##  4 ""                      Sense & Sensibility          4       0
##  5 "(1811)"                Sense & Sensibility          5       0
##  6 ""                      Sense & Sensibility          6       0
##  7 ""                      Sense & Sensibility          7       0
##  8 ""                      Sense & Sensibility          8       0
##  9 ""                      Sense & Sensibility          9       0
## 10 "CHAPTER 1"             Sense & Sensibility         10       1
## # ℹ 73,412 more rows

To work with this as a tidy dataset, we need to restructure it in the one-token-per-row format, which as we saw earlier is done with the unnest_tokens() function

tidy_books <- original_books %>%
  unnest_tokens(word, text)

tidy_books

## # A tibble: 725,055 × 4
##    book                linenumber chapter word       
##    <fct>                    <int>   <int> <chr>      
##  1 Sense & Sensibility          1       0 sense      
##  2 Sense & Sensibility          1       0 and        
##  3 Sense & Sensibility          1       0 sensibility
##  4 Sense & Sensibility          3       0 by         
##  5 Sense & Sensibility          3       0 jane       
##  6 Sense & Sensibility          3       0 austen     
##  7 Sense & Sensibility          5       0 1811       
##  8 Sense & Sensibility         10       1 chapter    
##  9 Sense & Sensibility         10       1 1          
## 10 Sense & Sensibility         13       1 the        
## # ℹ 725,045 more rows

This function uses the okenizers package to seperate each line of text in the original dataframe into tokens.

The default tokenizing is for words, but other options including characters, n-grams, sentences, lines, or paragraphs can be used.

Now that the data is in a one-word-per-row format, we can manipulate it with tools like dplyr

Often in text analysis, we will want to remove stop words. Stop words are words that are NOT USEFUL for analysis. These include words like the, of, to, and, and so forth

We can remove stop words (kept in the tidytext dataset ‘stop_words’) with an anti_join()

data(stop_words)

tidy_books <- tidy_books %>%
  anti_join(stop_words)

## Joining with `by = join_by(word)`

The stop words dataset in the tidytext package contains stop words from three lexicons. We can use them all together, as we have here, or filter() to only use one set of stop words if that’s more appropriate for your analysis.

tidy_books %>%
  count(word, sort = TRUE)

## # A tibble: 13,914 × 2
##    word       n
##    <chr>  <int>
##  1 miss    1855
##  2 time    1337
##  3 fanny    862
##  4 dear     822
##  5 lady     817
##  6 sir      806
##  7 day      797
##  8 emma     787
##  9 sister   727
## 10 house    699
## # ℹ 13,904 more rows

Because we’ve been using tidy tools, our word counts are stored in a tidy data frame, This allows us to pipe this directly into ggplot2. For example, we can create a visualization of the most common words.

library(ggplot2)

tidy_books %>%
  count(word, sort = TRUE) %>%
  filter(n > 600) %>%
  mutate(word = reorder(word, n)) %>%
  ggplot(aes(n, word)) +
  geom_col() +
  labs(y = NULL, x = "Word Count")

The gutenbergr package

This package provides access to the public domain works from the gutenberg project (www.gutenberg.org) This package provides tools for both downloading books and a complete dataset of project gutenberg metadata that can be used to find works of interest. We will mostly use the function gutenberg_download().

Word frequencies

Lets look at some biology texts, starting with Darwin

The Voyage of the Beagle - 944 On the Origin of Species by the Means of Natural Selection - 1228 The Expression of Emotions in Man and Animals - 1227 The Descent of Man, and Selection in Relation to Sex - 2300

We can access these works using the gutenberg_download() and the ID numbers

(Some gutenberg texts broke when run)

library(gutenbergr)

## Warning: package 'gutenbergr' was built under R version 4.4.3

library(readr)
library(readtext)

## Warning: package 'readtext' was built under R version 4.4.3

#darwin <- gutenberg_download(c(944, 1227, 1228, 2300), mirror = "http://mirrors.xmission.com/gutenberg/")

darwin <- readtext(c("944.txt", "1227.txt", "1228.txt", "2300.txt"))

Lets break these into tokens

tidy_darwin <- darwin %>%
  unnest_tokens(word, text) %>%
  anti_join(stop_words)

## Joining with `by = join_by(word)`

Lets check out what the most common darwin words are.

tidy_darwin %>%
  count(word, sort = TRUE)

## readtext object consisting of 23774 documents and 0 docvars.
## # A data frame: 23,774 × 3
##   word        n text     
##   <chr>   <int> <chr>    
## 1 species  2996 "\"\"..."
## 2 male     1672 "\"\"..."
## 3 males    1337 "\"\"..."
## 4 animals  1314 "\"\"..."
## 5 birds    1291 "\"\"..."
## 6 female   1197 "\"\"..."
## # ℹ 23,768 more rows

Now lets get some work from Thomas Hunt Morgan, who is credited with discovering chromosomes.

Regenerataion - 57198 The Genetic and Operative Evidence Relating to Secondary Sexual Characteristics - 57460 Evolution and Adaptation - 63540

morgan <- gutenberg_download(c(57198, 57460, 63540), mirror = "http://mirrors.xmission.com/gutenberg/")

Lets tokenize THM

tidy_morgan <- morgan %>%
  unnest_tokens(word, text) %>%
  anti_join(stop_words)

## Joining with `by = join_by(word)`

What are TJM’s most common words?

tidy_morgan %>%
  count(word, sort = TRUE)

## # A tibble: 13,855 × 2
##    word             n
##    <chr>        <int>
##  1 species        869
##  2 regeneration   814
##  3 piece          702
##  4 cut            669
##  5 male           668
##  6 forms          631
##  7 selection      604
##  8 cells          576
##  9 found          552
## 10 development    546
## # ℹ 13,845 more rows

Lastly, lets look at Thomas Henry Huxley

Evidence as to Man’s Place in Nature - 2931 On the Reception of the Origin of Species - 2089 Evolution and Ethics, and Other Essays - 2940 Science and Culture, and Other Essays - 52344

huxley <- gutenberg_download(c(2931,2089, 2940, 52344), mirror = "http://mirrors.xmission.com/gutenberg/")

tidy_huxley <- huxley %>%
  unnest_tokens(word, text) %>%
  anti_join(stop_words)

## Joining with `by = join_by(word)`

tidy_huxley %>%
  count(word, sort = TRUE)

## # A tibble: 11,749 × 2
##    word          n
##    <chr>     <int>
##  1 knowledge   194
##  2 de          174
##  3 science     169
##  4 animals     162
##  5 life        160
##  6 time        157
##  7 animal      154
##  8 species     132
##  9 body        127
## 10 nature      126
## # ℹ 11,739 more rows

Now, lets calculate the frequency for each word for the works of Darwin, Morgan, and Huxley by binding the dataframes together

library(tidyr)

frequency <- bind_rows(mutate(tidy_morgan, author = "Thomas Hunt Morgan"),
                       mutate(tidy_darwin, author = "Charles Darwin"),
                       mutate(tidy_huxley, author = "Thomas Henry Huxley")) %>%
  mutate(word = str_extract(word, "[a-z']+")) %>%
  count(author, word) %>%
  group_by(author) %>%
  mutate(proportion = n / sum(n)) %>%
  select(-n) %>%
  pivot_wider(names_from = author, values_from = proportion) %>%
  pivot_longer('Thomas Hunt Morgan': 'Charles Darwin', names_to = "author", values_to = "proportion")

frequency

## # A tibble: 90,042 × 3
##    word    author               proportion
##    <chr>   <chr>                     <dbl>
##  1 a       Thomas Hunt Morgan   0.00206   
##  2 a       Thomas Henry Huxley  0.000151  
##  3 a       Charles Darwin       0.000139  
##  4 ab      Thomas Hunt Morgan   0.000165  
##  5 ab      Thomas Henry Huxley  0.000173  
##  6 ab      Charles Darwin       0.00000631
##  7 abaiss  Thomas Hunt Morgan  NA         
##  8 abaiss  Thomas Henry Huxley NA         
##  9 abaiss  Charles Darwin       0.00000631
## 10 abandon Thomas Hunt Morgan   0.00000752
## # ℹ 90,032 more rows

Now we need to change the table so that each author has its own row

frequency2 <- pivot_wider(frequency, names_from = author, values_from = proportion)

frequency2

## # A tibble: 30,014 × 4
##    word        `Thomas Hunt Morgan` `Thomas Henry Huxley` `Charles Darwin`
##    <chr>                      <dbl>                 <dbl>            <dbl>
##  1 a                     0.00206                0.000151        0.000139  
##  2 ab                    0.000165               0.000173        0.00000631
##  3 abaiss               NA                     NA               0.00000631
##  4 abandon               0.00000752             0.0000216       0.00000315
##  5 abandoned             0.0000150              0.0000216       0.00000315
##  6 abashed              NA                     NA               0.00000315
##  7 abatement            NA                      0.0000216       0.00000315
##  8 abbot                NA                      0.0000432       0.00000315
##  9 abbott               NA                     NA               0.00000631
## 10 abbreviated          NA                     NA               0.0000126 
## # ℹ 30,004 more rows

Now lets plot

library(scales)

ggplot(frequency2, aes(x = `Charles Darwin`, y = `Thomas Hunt Morgan`), color = abs(- "Charles Darwin" - "Thomas Hunt Morgan")) +
  geom_abline(color = "gray40", lty = 2) +
  geom_jitter(alpha = 0.1, size = 2.5, width = 0.3, height = 0.3) +
  geom_text(aes(label = word), check_overlap = TRUE, vjust = 1.5, color = "red") +
  scale_x_log10(labels = percent_format()) +
  scale_y_log10(labels = percent_format()) +
  scale_color_gradient(limits = c(0, 0.001), low = "darkslategray4", high = "gray75") +
  theme(legend.position = "none") +
  labs(y = "Thomas Hunt Morgan", x = "Charles Darwin")

## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • color = abs(-"Charles Darwin" - "Thomas Hunt Morgan")
## ℹ Did you misspell an argument name?

## Warning: Removed 22541 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 22542 rows containing missing values or values outside the scale range
## (`geom_text()`).

ggplot(frequency2, aes(x = `Charles Darwin`, y = `Thomas Henry Huxley`), color = abs(- "Charles Darwin" - "Thomas Henry Huxley")) +
  geom_abline(color = "gray40", lty = 2) +
  geom_jitter(alpha = 0.1, size = 2.5, width = 0.3, height = 0.3) +
  geom_text(aes(label = word), check_overlap = TRUE, vjust = 1.5, color = "red") +
  scale_x_log10(labels = percent_format()) +
  scale_y_log10(labels = percent_format()) +
  scale_color_gradient(limits = c(0, 0.001), low = "darkslategray4", high = "gray75") +
  theme(legend.position = "none") +
  labs(y = "Thomas Henry Huxley", x = "Charles Darwin")

## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • color = abs(-"Charles Darwin" - "Thomas Henry Huxley")
## ℹ Did you misspell an argument name?

## Warning: Removed 23166 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 23167 rows containing missing values or values outside the scale range
## (`geom_text()`).

ggplot(frequency2, aes(x = `Thomas Hunt Morgan`, y = `Thomas Henry Huxley`), color = abs(- "Thomas Hunt Morgan" - "Thomas Henry Huxley")) +
  geom_abline(color = "gray40", lty = 2) +
  geom_jitter(alpha = 0.1, size = 2.5, width = 0.3, height = 0.3) +
  geom_text(aes(label = word), check_overlap = TRUE, vjust = 1.5, color = "red") +
  scale_x_log10(labels = percent_format()) +
  scale_y_log10(labels = percent_format()) +
  scale_color_gradient(limits = c(0, 0.001), low = "darkslategray4", high = "gray75") +
  theme(legend.position = "none") +
  labs(y = "Thomas Henry Huxley", x = "Thomas Hunt Morgan")

## Warning in fortify(data, ...): Arguments in `...` must be used.
## ✖ Problematic argument:
## • color = abs(-"Thomas Hunt Morgan" - "Thomas Henry Huxley")
## ℹ Did you misspell an argument name?

## Warning: Removed 25093 rows containing missing values or values outside the scale range
## (`geom_point()`).

## Warning: Removed 25094 rows containing missing values or values outside the scale range
## (`geom_text()`).

Sentiment Analysis

The Sentiments datasets

There are a variety of methods and dictionaries that exist for evaluating the opinion or emotion of the text.

AFINN bing nrc

bing categorizes words in a binary fashion into positive or negative nrc categorizes into positive, negative, anger, anticipation, disgust, fear, joy, sadness, surprise, and trust AFFIN assigns a score between -5 and 5, with negative indicating negative sentiment, and 5, positive

The function get_sentiments() allows us to get the specific sentiments lexicon with the measures for each one

library(tidytext)
library(textdata)

## Warning: package 'textdata' was built under R version 4.4.3

afinn <- get_sentiments("afinn")

afinn

## # A tibble: 2,477 × 2
##    word       value
##    <chr>      <dbl>
##  1 abandon       -2
##  2 abandoned     -2
##  3 abandons      -2
##  4 abducted      -2
##  5 abduction     -2
##  6 abductions    -2
##  7 abhor         -3
##  8 abhorred      -3
##  9 abhorrent     -3
## 10 abhors        -3
## # ℹ 2,467 more rows

Lets look at bing

bing <- get_sentiments("bing")

bing

## # A tibble: 6,786 × 2
##    word        sentiment
##    <chr>       <chr>    
##  1 2-faces     negative 
##  2 abnormal    negative 
##  3 abolish     negative 
##  4 abominable  negative 
##  5 abominably  negative 
##  6 abominate   negative 
##  7 abomination negative 
##  8 abort       negative 
##  9 aborted     negative 
## 10 aborts      negative 
## # ℹ 6,776 more rows

And lastly nrc

nrc <- get_sentiments("nrc")

nrc

## # A tibble: 13,872 × 2
##    word        sentiment
##    <chr>       <chr>    
##  1 abacus      trust    
##  2 abandon     fear     
##  3 abandon     negative 
##  4 abandon     sadness  
##  5 abandoned   anger    
##  6 abandoned   fear     
##  7 abandoned   negative 
##  8 abandoned   sadness  
##  9 abandonment anger    
## 10 abandonment fear     
## # ℹ 13,862 more rows

These libraries were created either using crowdsourcing or cloud computing/ai like Amazon Mechanical Turk, or by labor of one of the authors, and then validated with crowdsourcing

Lets look at the words with a joy score from NRC

library(gutenbergr)
library(dplyr)
library(stringr)

morgan <- gutenberg_download(c(57198, 57460, 63540), mirror = "http://mirrors.xmission.com/gutenberg/")

tidy_books <- morgan %>%
  group_by(gutenberg_id) %>%
  mutate(linenumber = row_number(), chapter = cumsum(str_detect(text, regex("^chapter [\\divxlc]", ignore_case = TRUE)))) %>%
  ungroup() %>%
  unnest_tokens(word, text)

tidy_books

## # A tibble: 350,473 × 4
##    gutenberg_id linenumber chapter word        
##           <int>      <int>   <int> <chr>       
##  1        57198          1       0 regeneration
##  2        57198          6       0 columbia    
##  3        57198          6       0 university  
##  4        57198          6       0 biological  
##  5        57198          6       0 series      
##  6        57198          8       0 edited      
##  7        57198          8       0 by          
##  8        57198          9       0 henry       
##  9        57198          9       0 fairfield   
## 10        57198          9       0 osborn      
## # ℹ 350,463 more rows

Lets add teh book name instead of GID

colnames(tidy_books)[1] <- "book"

tidy_books$book[tidy_books$book == 57198] <- "Regeneration"
tidy_books$book[tidy_books$book == 57460] <- "The Genetic and Operative Evidence Relating to Secondary Sexual Characteristics"
tidy_books$book[tidy_books$book == 63540] <- "Evolution and Adaptation"

tidy_books

## # A tibble: 350,473 × 4
##    book         linenumber chapter word        
##    <chr>             <int>   <int> <chr>       
##  1 Regeneration          1       0 regeneration
##  2 Regeneration          6       0 columbia    
##  3 Regeneration          6       0 university  
##  4 Regeneration          6       0 biological  
##  5 Regeneration          6       0 series      
##  6 Regeneration          8       0 edited      
##  7 Regeneration          8       0 by          
##  8 Regeneration          9       0 henry       
##  9 Regeneration          9       0 fairfield   
## 10 Regeneration          9       0 osborn      
## # ℹ 350,463 more rows

Now that we have a tidy format with one word per row, we are ready for sentiment analysis. First lets use NRC

nrc_joy <- get_sentiments("nrc") %>%
  filter(sentiment == "joy")

tidy_books %>%
  filter(book == "Regeneration") %>%
  inner_join(nrc_joy) %>%
  count(word, sort = TRUE)

## Joining with `by = join_by(word)`

## # A tibble: 88 × 2
##    word             n
##    <chr>        <int>
##  1 found          271
##  2 present        194
##  3 kind            95
##  4 food            80
##  5 organization    72
##  6 organ           69
##  7 grow            67
##  8 true            52
##  9 special         35
## 10 alive           28
## # ℹ 78 more rows

We can also examine how sentiment changes throughout a work.

library(tidyr)

Thomas_Hunt_Morgan_sentiment <- tidy_books %>%
  inner_join(get_sentiments("bing")) %>%
  count(book, index = linenumber %/% 80, sentiment) %>%
  pivot_wider(names_from = sentiment, values_from = n, values_fill = 0) %>%
  mutate(sentiment = positive - negative)

## Joining with `by = join_by(word)`

Now lets plot it

library(ggplot2)

ggplot(Thomas_Hunt_Morgan_sentiment, aes(index, sentiment, fill = book)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~book, ncol = 2, scales = "free_x")

Lets compare the three sentiment dictionaries

There are several options for sentiment lexicons, you might want some more info on which is appropiate for your purpose, Here we will use all three of our dictionaries and examine how the sentiment changes across the arc of Regeneration

regeneration <- tidy_books %>%
  filter(book == "Regeneration")
regeneration

## # A tibble: 138,523 × 4
##    book         linenumber chapter word        
##    <chr>             <int>   <int> <chr>       
##  1 Regeneration          1       0 regeneration
##  2 Regeneration          6       0 columbia    
##  3 Regeneration          6       0 university  
##  4 Regeneration          6       0 biological  
##  5 Regeneration          6       0 series      
##  6 Regeneration          8       0 edited      
##  7 Regeneration          8       0 by          
##  8 Regeneration          9       0 henry       
##  9 Regeneration          9       0 fairfield   
## 10 Regeneration          9       0 osborn      
## # ℹ 138,513 more rows

Lets again use interger division (%/%) to define larger sectionos of the text that span multiple lines, and we can use the same pattern with count(), pivot_wider(), and mutate() to find the net setement in each of these section of text

afinn <- regeneration %>%
  inner_join(get_sentiments("afinn")) %>%
  group_by(index = linenumber %/% 80) %>%
  summarise(sentiment = sum(value)) %>%
  mutate(method = "AFINN")

## Joining with `by = join_by(word)`

bing_and_nrc <- bind_rows(
  regeneration %>%
    inner_join(get_sentiments("bing")) %>%
    mutate(method = "Bing et al."),
  regeneration %>%
    inner_join(get_sentiments("nrc") %>%
                 filter(sentiment %in% c("positive", "negative"))
               ) %>%
    mutate(method = "NRC")) %>%
  count(method, index = linenumber %/% 80, sentiment) %>%
  pivot_wider(names_from = sentiment,
              values_from = n,
              values_fill = 0) %>%
  mutate(sentiment = positive - negative)

## Joining with `by = join_by(word)`
## Joining with `by = join_by(word)`

## Warning in inner_join(., get_sentiments("nrc") %>% filter(sentiment %in% : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 146 of `x` matches multiple rows in `y`.
## ℹ Row 5299 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

We can now estimate the net sentiment (positive - negative) in each chunk of the novel text for each lexicon (dictionary). Lets bind them all together and visualize with ggplot

bind_rows(afinn, bing_and_nrc) %>%
  ggplot(aes(index, sentiment, fill = method)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~method, ncol = 1, scales = "free_y")

Lets look at the counts based on each dictionary

get_sentiments("nrc") %>%
  filter(sentiment %in% c("positive", "negative")) %>%
  count(sentiment)

## # A tibble: 2 × 2
##   sentiment     n
##   <chr>     <int>
## 1 negative   3316
## 2 positive   2308

get_sentiments("bing") %>%
  count(sentiment)

## # A tibble: 2 × 2
##   sentiment     n
##   <chr>     <int>
## 1 negative   4781
## 2 positive   2005

bing_word_counts <- tidy_books %>%
  inner_join(get_sentiments("bing")) %>%
  count(word, sentiment, sort = TRUE) %>%
  ungroup()

## Joining with `by = join_by(word)`

bing_word_counts

## # A tibble: 1,283 × 3
##    word      sentiment     n
##    <chr>     <chr>     <int>
##  1 like      positive    324
##  2 well      positive    266
##  3 important positive    199
##  4 regard    positive    174
##  5 die       negative    160
##  6 great     positive    159
##  7 better    positive    121
##  8 work      positive    121
##  9 doubt     negative    113
## 10 lost      negative    112
## # ℹ 1,273 more rows

This can be shown visually, and we can pipe straight into ggplot2

bing_word_counts %>%
  group_by(sentiment) %>%
  slice_max(n, n = 10) %>%
  ungroup() %>%
  mutate(word = reorder(word, n)) %>%
  ggplot(aes(n, word, fill = sentiment)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~sentiment, scale = "free_y") +
  labs(x = "Contribution to Sentiment", y = NULL)

Lets spot an anomaly in the dataset.

custom_stop_words <- bind_rows(tibble(word = c("wild", "dark", "great", "like"), lexicon = c("custom")), stop_words)

custom_stop_words

## # A tibble: 1,153 × 2
##    word      lexicon
##    <chr>     <chr>  
##  1 wild      custom 
##  2 dark      custom 
##  3 great     custom 
##  4 like      custom 
##  5 a         SMART  
##  6 a's       SMART  
##  7 able      SMART  
##  8 about     SMART  
##  9 above     SMART  
## 10 according SMART  
## # ℹ 1,143 more rows

Word Clouds!!

We can see that tidy text mining and sentiment analysis works well with ggplot2, but having our data in tidy format leads to other nice graphing techniques

Lets use the wordcloud package

library(wordcloud)

## Warning: package 'wordcloud' was built under R version 4.4.3

## 
## Attaching package: 'wordcloud'

## The following object is masked from 'package:gplots':
## 
##     textplot

tidy_books %>%
  anti_join(custom_stop_words) %>%
  count(word) %>%
  with(wordcloud(word, n, max.words = 100, scale = c(2, 0.05)))

## Joining with `by = join_by(word)`

Lets also look at comparison.cloud(), which may require turning the dataframe into a matrix

We can change to matrix using the acast() function

library(reshape2)

## Warning: package 'reshape2' was built under R version 4.4.3

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

tidy_books %>%
  inner_join(get_sentiments("bing")) %>%
  count(word, sentiment, sort = TRUE) %>%
  acast(word ~sentiment, value.var = "n", fill = 0) %>%
  comparison.cloud(colors = c("gray20", "gray80"), max.words = 100, scale = c(2, 0.05))

## Joining with `by = join_by(word)`

Looking at units beyond words

Lots of useful work can be done by tokenizing at the word level, but sometimes its nice to look at different units of text. For example, we can look beyond just unigrams.

Ex. I am not having a good day.

bingnegative <- get_sentiments("bing") %>%
  filter(sentiment == "negative")

wordcounts <- tidy_books %>%
  group_by(book, chapter) %>%
  summarize(word = n())

## `summarise()` has grouped output by 'book'. You can override using the
## `.groups` argument.

tidy_books %>%
  semi_join(bingnegative) %>%
  group_by(book, chapter) %>%
  summarize(negativewords = n()) %>%
  left_join(wordcounts, by = c("book", "chapter")) %>%
  mutate(ratio = negativewords / word) %>%
  #filter(chapter != 0) %>%
  slice_max(ratio, n = 1) %>%
  ungroup()

## Joining with `by = join_by(word)`
## `summarise()` has grouped output by 'book'. You can override using the
## `.groups` argument.

## # A tibble: 3 × 5
##   book                                       chapter negativewords   word  ratio
##   <chr>                                        <int>         <int>  <int>  <dbl>
## 1 Evolution and Adaptation                         0          2766 158568 0.0174
## 2 Regeneration                                     5           226   8394 0.0269
## 3 The Genetic and Operative Evidence Relati…       0           893  53382 0.0167

N-grams

So far we’ve only looked at single words, but many interesting (more accurate) analyses are based on the relationship between words

Lets look at some methods of tidytext for calculating and visualizing word relationships.

library(dplyr)
library(tidytext)

morgan_books <- gutenberg_download(c(57198, 57460, 63540), mirror = "http://mirrors.xmission.com/gutenberg/")

colnames(morgan_books)[1] <- "book"

morgan_books$book[morgan_books$book == 57198] <- "Regeneration"
morgan_books$book[morgan_books$book == 57460] <- "The Genetic and Operative Evidence Relating to Secondary Sexual Characteristics"
morgan_books$book[morgan_books$book == 63540] <- "Evolution and Adaptation"

morgan_bigrams <- morgan_books %>%
  unnest_tokens(bigram, text, token = "ngrams", n = 2)

morgan_bigrams

## # A tibble: 323,605 × 2
##    book         bigram               
##    <chr>        <chr>                
##  1 Regeneration <NA>                 
##  2 Regeneration <NA>                 
##  3 Regeneration <NA>                 
##  4 Regeneration <NA>                 
##  5 Regeneration <NA>                 
##  6 Regeneration columbia university  
##  7 Regeneration university biological
##  8 Regeneration biological series    
##  9 Regeneration <NA>                 
## 10 Regeneration edited by            
## # ℹ 323,595 more rows

This data is still in tidytext format, and is structured as one-token-per-row. Each token is a bigram

Counting and filtering n-gram

morgan_bigrams %>%
  count(bigram, sort = TRUE)

## # A tibble: 108,762 × 2
##    bigram       n
##    <chr>    <int>
##  1 of the    6514
##  2 <NA>      5317
##  3 in the    3018
##  4 that the  1474
##  5 to the    1452
##  6 it is     1011
##  7 on the     938
##  8 of a       821
##  9 from the   818
## 10 the same   806
## # ℹ 108,752 more rows

Most of the common bigrams are stop-words. This can be a good time to use tidyr’s seperate command, which splits a column into multiples based on a delimiter. This will let us make a column for word one and two

library(tidyr)

bigrams_seperated <- morgan_bigrams %>%
  separate(bigram, c("word1", "word2"), sep = " ")

bigrams_filtered <- bigrams_seperated %>%
  filter(!word1 %in% stop_words$word) %>%
  filter(!word2 %in% stop_words$word)
bigrams_filtered

## # A tibble: 42,464 × 3
##    book         word1      word2     
##    <chr>        <chr>      <chr>     
##  1 Regeneration <NA>       <NA>      
##  2 Regeneration <NA>       <NA>      
##  3 Regeneration <NA>       <NA>      
##  4 Regeneration <NA>       <NA>      
##  5 Regeneration <NA>       <NA>      
##  6 Regeneration columbia   university
##  7 Regeneration university biological
##  8 Regeneration biological series    
##  9 Regeneration <NA>       <NA>      
## 10 Regeneration henry      fairfield 
## # ℹ 42,454 more rows

New bigram counts

bigram_counts <- bigrams_filtered %>%
  unite(bigram, word1, word2, sep = " ")

bigram_counts

## # A tibble: 42,464 × 2
##    book         bigram               
##    <chr>        <chr>                
##  1 Regeneration NA NA                
##  2 Regeneration NA NA                
##  3 Regeneration NA NA                
##  4 Regeneration NA NA                
##  5 Regeneration NA NA                
##  6 Regeneration columbia university  
##  7 Regeneration university biological
##  8 Regeneration biological series    
##  9 Regeneration NA NA                
## 10 Regeneration henry fairfield      
## # ℹ 42,454 more rows

We may also be interested in trigrams, which are three word combos

trigrams <- morgan_books %>%
  unnest_tokens(trigram, text, token = "ngrams", n = 3) %>%
  separate(trigram, c("word1", "word2", "word3"), sep = " ") %>%
  filter(!word1 %in% stop_words$word,
         !word2 %in% stop_words$word,
         !word3 %in% stop_words$word) %>%
  count(word1, word2, word3, sort = TRUE)

trigrams

## # A tibble: 10,841 × 4
##    word1     word2     word3           n
##    <chr>     <chr>     <chr>       <int>
##  1 <NA>      <NA>      <NA>         6499
##  2 secondary sexual    characters     73
##  3 hh        hh        hh             22
##  4 hen       feathered males          18
##  5 black     breasted  game           13
##  6 castrated sebright  male           10
##  7 secondary sexual    character      10
##  8 secondary sexual    differences    10
##  9 entw      mech      ii              9
## 10 plate     1         figure          9
## # ℹ 10,831 more rows

Lets analyze some bigrams

bigrams_filtered %>%
  filter(word2 == "selection") %>%
  count(book, word1, sort = TRUE)

## # A tibble: 32 × 3
##    book                                                              word1     n
##    <chr>                                                             <chr> <int>
##  1 Evolution and Adaptation                                          natu…   130
##  2 Evolution and Adaptation                                          sexu…    77
##  3 Regeneration                                                      natu…    30
##  4 The Genetic and Operative Evidence Relating to Secondary Sexual … sexu…    24
##  5 The Genetic and Operative Evidence Relating to Secondary Sexual … natu…    18
##  6 Evolution and Adaptation                                          arti…    16
##  7 Evolution and Adaptation                                          germ…     6
##  8 Evolution and Adaptation                                          cont…     4
##  9 Evolution and Adaptation                                          indi…     3
## 10 Evolution and Adaptation                                          unco…     3
## # ℹ 22 more rows

Lets again look at tf-idf across bigrams across Morgan’s works

bigram_tf_idf <- bigram_counts %>%
  count(book, bigram) %>%
  bind_tf_idf(bigram, book, n) %>%
  arrange(desc(tf_idf))

bigram_tf_idf

## # A tibble: 26,779 × 6
##    book                                       bigram     n      tf   idf  tf_idf
##    <chr>                                      <chr>  <int>   <dbl> <dbl>   <dbl>
##  1 The Genetic and Operative Evidence Relati… hen f…    73 0.00845 1.10  0.00928
##  2 The Genetic and Operative Evidence Relati… cock …    49 0.00567 1.10  0.00623
##  3 The Genetic and Operative Evidence Relati… hen f…    42 0.00486 1.10  0.00534
##  4 The Genetic and Operative Evidence Relati… cock …    36 0.00417 1.10  0.00458
##  5 The Genetic and Operative Evidence Relati… sex l…    35 0.00405 1.10  0.00445
##  6 Regeneration                               illus…    68 0.00364 1.10  0.00400
##  7 The Genetic and Operative Evidence Relati… hh hh     31 0.00359 1.10  0.00394
##  8 The Genetic and Operative Evidence Relati… secon…    80 0.00926 0.405 0.00375
##  9 The Genetic and Operative Evidence Relati… inter…    27 0.00312 1.10  0.00343
## 10 The Genetic and Operative Evidence Relati… feath…    25 0.00289 1.10  0.00318
## # ℹ 26,769 more rows

bigram_tf_idf %>%
  arrange(desc(tf_idf)) %>%
  group_by(book) %>%
  slice_max(tf_idf, n = 10) %>%
  ungroup() %>%
  mutate(bigram = reorder(bigram, tf_idf)) %>%
  ggplot(aes(tf_idf, bigram, fill = book)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~book, ncol = 2, scales = "free") +
  labs(x = "tf-idf of bigrams", y = NULL)

Using bigrams to provide context in sentiment analysis

bigrams_seperated %>%
  filter(word1 == "not") %>%
  count(word1, word2, sort = TRUE)

## # A tibble: 520 × 3
##    word1 word2     n
##    <chr> <chr> <int>
##  1 not   be      114
##  2 not   only     83
##  3 not   the      73
##  4 not   a        49
##  5 not   in       49
##  6 not   so       41
##  7 not   been     39
##  8 not   have     38
##  9 not   seem     35
## 10 not   know     34
## # ℹ 510 more rows

By doing sentiment analysis of bigrams, we can examine how often sentiment-associated words are preceded by a modifier like “not” or other negating words.

AFINN <- get_sentiments("afinn")

AFINN

## # A tibble: 2,477 × 2
##    word       value
##    <chr>      <dbl>
##  1 abandon       -2
##  2 abandoned     -2
##  3 abandons      -2
##  4 abducted      -2
##  5 abduction     -2
##  6 abductions    -2
##  7 abhor         -3
##  8 abhorred      -3
##  9 abhorrent     -3
## 10 abhors        -3
## # ℹ 2,467 more rows

We can examine the most frequent words that were perceded by “not”, and associate with sentiment.

not_words <- bigrams_seperated %>%
  filter(word1 == "not") %>%
  inner_join(AFINN, by = c(word2 = "word")) %>%
  count(word2, value, sort = TRUE)

not_words

## # A tibble: 69 × 3
##    word2      value     n
##    <chr>      <dbl> <int>
##  1 clear          1     8
##  2 difficult     -1     8
##  3 fail          -2     5
##  4 pretend       -1     5
##  5 accept         1     4
##  6 increase       1     4
##  7 leave         -1     4
##  8 affected      -1     3
##  9 determined     2     3
## 10 forget        -1     3
## # ℹ 59 more rows

Lets visualize

library(ggplot2)

not_words %>%
  mutate(contribution = n * value) %>%
  arrange(desc(abs(contribution))) %>%
  head(20) %>%
  mutate(word2 = reorder(word2, contribution)) %>%
  ggplot(aes(n * value, word2, fill = n * value > 0)) +
  geom_col(show.legend = FALSE) +
  labs(x = "Sentiment Value * Number of Occurences", y = "Words Preceded by \"not\"")

negation_words <- c("not", "no", "never", "non", "without")

negated_words <- bigrams_seperated %>%
  filter(word1 %in% negation_words) %>%
  inner_join(AFINN, by = c(word2 = "word")) %>%
  count(word1, word2, value, sort = TRUE)

negated_words

## # A tibble: 98 × 4
##    word1 word2     value     n
##    <chr> <chr>     <dbl> <int>
##  1 no    doubt        -1    48
##  2 not   clear         1     8
##  3 not   difficult    -1     8
##  4 not   fail         -2     5
##  5 not   pretend      -1     5
##  6 no    advantage     2     4
##  7 no    good          3     4
##  8 no    great         3     4
##  9 no    matter        1     4
## 10 not   accept        1     4
## # ℹ 88 more rows

Lets visualize the negation words

negated_words %>%
  mutate(contribution = n * value,
         word2 = reorder(paste(word2, word1, sep = "_"), contribution)) %>%
  group_by(word1) %>%
  slice_max(abs(contribution), n = 12, with_ties = FALSE) %>%
  ggplot(aes(word2, contribution, fill = n * value > 0)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~word1, scales = "free") +
  scale_x_discrete(labels = function(x) gsub("_.+$", "", x)) +
  xlab("Words Preceded by Negation Term") +
  ylab("Sentiment Value * # of Occurences") +
  coord_flip()

Visualize a network of bigrams with ggraph

library(igraph)

## Warning: package 'igraph' was built under R version 4.4.3

## 
## Attaching package: 'igraph'

## The following objects are masked from 'package:lubridate':
## 
##     %--%, union

## The following object is masked from 'package:plotly':
## 
##     groups

## The following object is masked from 'package:tidyr':
## 
##     crossing

## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union

## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum

## The following object is masked from 'package:base':
## 
##     union

bigrams_filtered <- bigrams_filtered %>%
  filter(!is.na(word1))

bigram_counts <- bigrams_filtered %>%
  count(word1, word2, sort = TRUE)

bigram_graph <- bigram_counts %>%
  filter(n > 20) %>%
  graph_from_data_frame()

bigram_graph

## IGRAPH 9281772 DN-- 78 49 -- 
## + attr: name (v/c), n (e/n)
## + edges from 9281772 (vertex names):
##  [1] natural     ->selection  secondary   ->sexual     sexual      ->selection 
##  [4] sexual      ->characters entw        ->mech       hen         ->feathered 
##  [7] cut         ->surface    illustration->fig        digestive   ->tract     
## [10] über        ->die        cock        ->feathered  germ        ->cells     
## [13] de          ->vries      regeneration->takes      hen         ->feathering
## [16] sea         ->urchin     nervous     ->system     breaking    ->joint     
## [19] external    ->conditions cock        ->feathering fluctuating ->variations
## [22] nerve       ->cord       sex         ->linked     cut         ->edge      
## + ... omitted several edges

library(ggraph)

## Warning: package 'ggraph' was built under R version 4.4.3

set.seed(1234)

ggraph(bigram_graph, layout = "fr") +
  geom_edge_link() +
  geom_node_point() +
  geom_node_text(aes(label = name), vjust = 1, hjust = 1)

We can also add directionality to this network

set.seed(1234)

a <- grid::arrow(type = "closed", length = unit(0.15, "inches"))

ggraph(bigram_graph, layout = "fr") +
  geom_edge_link(aes(edge_alpha = n), show.legend = FALSE, arrow = a, end_cap = circle(0.07, 'inches')) +
  geom_node_point(color = "lightblue", size = 3) +
  geom_node_text(aes(label = name), vjust = 1, hjust = 1) +
  theme_void()

Word Frequencies

A central question in text mining is how to quantify what a document is about. We can do this by looking at words that make up the document, and measuring term frequency.

There are a lot of words that may not be important, these are the stop words.

One way to remedy this is to look at inverse document frequency words, which decreases the weight for commonly used words and increases the weight for words that are not used very much.

Term frequency in Morgan’s works

library(dplyr)
library(tidytext)

book_words <- gutenberg_download(c(57198, 57460, 63540), mirror = "http://mirrors.xmission.com/gutenberg/")

colnames(book_words)[1] <- "book"

book_words$book[book_words$book == 57198] <- "Regeneration"
book_words$book[book_words$book == 57460] <- "The Genetic and Operative Evidence Relating to Secondary Sexual Characteristics"
book_words$book[book_words$book == 63540] <- "Evolution and Adaptation"

Now lets disect

book_words <- book_words %>%
  unnest_tokens(word, text) %>%
  count(book, word, sort = TRUE)

book_words

## # A tibble: 21,787 × 3
##    book                                                              word      n
##    <chr>                                                             <chr> <int>
##  1 Evolution and Adaptation                                          the   13965
##  2 Regeneration                                                      the   12981
##  3 Evolution and Adaptation                                          of     8088
##  4 Regeneration                                                      of     7232
##  5 The Genetic and Operative Evidence Relating to Secondary Sexual … the    5033
##  6 Evolution and Adaptation                                          in     4522
##  7 Regeneration                                                      in     4014
##  8 Evolution and Adaptation                                          to     3975
##  9 Evolution and Adaptation                                          that   3333
## 10 Evolution and Adaptation                                          and    3197
## # ℹ 21,777 more rows

book_words$n <- as.numeric(book_words$n)

total_words <- book_words %>%
  group_by(book) %>%
  summarize(total = sum(n))

book_words

## # A tibble: 21,787 × 3
##    book                                                              word      n
##    <chr>                                                             <chr> <dbl>
##  1 Evolution and Adaptation                                          the   13965
##  2 Regeneration                                                      the   12981
##  3 Evolution and Adaptation                                          of     8088
##  4 Regeneration                                                      of     7232
##  5 The Genetic and Operative Evidence Relating to Secondary Sexual … the    5033
##  6 Evolution and Adaptation                                          in     4522
##  7 Regeneration                                                      in     4014
##  8 Evolution and Adaptation                                          to     3975
##  9 Evolution and Adaptation                                          that   3333
## 10 Evolution and Adaptation                                          and    3197
## # ℹ 21,777 more rows

book_words <- left_join(book_words, total_words)

## Joining with `by = join_by(book)`

book_words

## # A tibble: 21,787 × 4
##    book                                                       word      n  total
##    <chr>                                                      <chr> <dbl>  <dbl>
##  1 Evolution and Adaptation                                   the   13965 158568
##  2 Regeneration                                               the   12981 138523
##  3 Evolution and Adaptation                                   of     8088 158568
##  4 Regeneration                                               of     7232 138523
##  5 The Genetic and Operative Evidence Relating to Secondary … the    5033  53382
##  6 Evolution and Adaptation                                   in     4522 158568
##  7 Regeneration                                               in     4014 138523
##  8 Evolution and Adaptation                                   to     3975 158568
##  9 Evolution and Adaptation                                   that   3333 158568
## 10 Evolution and Adaptation                                   and    3197 158568
## # ℹ 21,777 more rows

You can see that the usual suspects are the most common words, but don’t tell us anything about what the book’s topic is.

library(ggplot2)

ggplot(book_words, aes(n / total, fill = book)) +
  geom_histogram(show.legend = FALSE) +
  xlim(NA, 0.0009) +
  facet_wrap(~book, ncol = 2, scale = "free_y")

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## Warning: Removed 437 rows containing non-finite outside the scale range
## (`stat_bin()`).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_bar()`).

Zipf’s Law

The frequency that a word appears is inversely proportional to its rank when predicting a topic.

Lets apply Zipf’s law to Morgan’s work

freq_by_rank <- book_words %>%
  group_by(book) %>%
  mutate(rank = row_number(),
         "term frequency" = n / total) %>%
  ungroup()
freq_by_rank

## # A tibble: 21,787 × 6
##    book                                word      n  total  rank `term frequency`
##    <chr>                               <chr> <dbl>  <dbl> <int>            <dbl>
##  1 Evolution and Adaptation            the   13965 158568     1           0.0881
##  2 Regeneration                        the   12981 138523     1           0.0937
##  3 Evolution and Adaptation            of     8088 158568     2           0.0510
##  4 Regeneration                        of     7232 138523     2           0.0522
##  5 The Genetic and Operative Evidence… the    5033  53382     1           0.0943
##  6 Evolution and Adaptation            in     4522 158568     3           0.0285
##  7 Regeneration                        in     4014 138523     3           0.0290
##  8 Evolution and Adaptation            to     3975 158568     4           0.0251
##  9 Evolution and Adaptation            that   3333 158568     5           0.0210
## 10 Evolution and Adaptation            and    3197 158568     6           0.0202
## # ℹ 21,777 more rows

freq_by_rank %>%
  ggplot(aes(rank, `term frequency`, color = book)) +
  geom_line(size = 1.1, alpha = 0.8, show.legend = FALSE) +
  scale_x_log10() +
  scale_y_log10()

Lets use TF - IDF to find words for each document by decreasing the weight for commonly used words and increasing the weight for words that are not used very much in a collection of documents

book_tf_idf <- book_words %>%
  bind_tf_idf(word, book, n)

book_tf_idf

## # A tibble: 21,787 × 7
##    book                                   word      n  total     tf   idf tf_idf
##    <chr>                                  <chr> <dbl>  <dbl>  <dbl> <dbl>  <dbl>
##  1 Evolution and Adaptation               the   13965 158568 0.0881     0      0
##  2 Regeneration                           the   12981 138523 0.0937     0      0
##  3 Evolution and Adaptation               of     8088 158568 0.0510     0      0
##  4 Regeneration                           of     7232 138523 0.0522     0      0
##  5 The Genetic and Operative Evidence Re… the    5033  53382 0.0943     0      0
##  6 Evolution and Adaptation               in     4522 158568 0.0285     0      0
##  7 Regeneration                           in     4014 138523 0.0290     0      0
##  8 Evolution and Adaptation               to     3975 158568 0.0251     0      0
##  9 Evolution and Adaptation               that   3333 158568 0.0210     0      0
## 10 Evolution and Adaptation               and    3197 158568 0.0202     0      0
## # ℹ 21,777 more rows

Lets look at terms with high tf-idf in Morgan’s works

book_tf_idf %>%
  select(-total) %>%
  arrange(desc(tf_idf))

## # A tibble: 21,787 × 6
##    book                                        word      n      tf   idf  tf_idf
##    <chr>                                       <chr> <dbl>   <dbl> <dbl>   <dbl>
##  1 The Genetic and Operative Evidence Relatin… feat…    85 1.59e-3 1.10  1.75e-3
##  2 The Genetic and Operative Evidence Relatin… cast…    78 1.46e-3 1.10  1.61e-3
##  3 The Genetic and Operative Evidence Relatin… males   183 3.43e-3 0.405 1.39e-3
##  4 The Genetic and Operative Evidence Relatin… plum…   163 3.05e-3 0.405 1.24e-3
##  5 The Genetic and Operative Evidence Relatin… hh       58 1.09e-3 1.10  1.19e-3
##  6 The Genetic and Operative Evidence Relatin… sebr…   144 2.70e-3 0.405 1.09e-3
##  7 The Genetic and Operative Evidence Relatin… barr…    53 9.93e-4 1.10  1.09e-3
##  8 The Genetic and Operative Evidence Relatin… sex     132 2.47e-3 0.405 1.00e-3
##  9 The Genetic and Operative Evidence Relatin… test…   128 2.40e-3 0.405 9.72e-4
## 10 The Genetic and Operative Evidence Relatin… feat…   125 2.34e-3 0.405 9.49e-4
## # ℹ 21,777 more rows

Lets look at a visualization for these high tf-idf words

library(forcats)

book_tf_idf %>%
  group_by(book) %>%
  slice_max(tf_idf, n = 15) %>%
  ungroup() %>%
  ggplot(aes(tf_idf, fct_reorder(word, tf_idf), fill = book)) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~book, ncol = 2, scales = "free") +
  labs(x = "tf-idf", y = NULL)