R: Advanced Topics
3 February 2017
During this session, we will primarily explore some more advanced topics within R, including alternative syntaxes and languages you can use from within R for various purposes, alternative ways of organizing and analyzing your information, and applications to more advanced analyses like text analysis, Bayesian analysis, or survey techniques.
dplyr
dplyr, a package maintained by Hadley Wickham, provides useful alternative functions for data exploration analogous to what we explored in the last session. It also provides an alternative syntax for performing a series of operations that can increase efficiency.
dplyr improves upon the previous plyr package.
Analysis and Manipulation
Suppose again that you load the QoG data we have used in previous classes.
dplyr provides a variety of functions that can allow us to see our data differently or to change it to suit our needs. Check out this vignette for more details.
For example, you can filter data to only investigate a subset.
# Filter to only include countries with imputed polity scores
# greater than 5
filter <- filter(qog, fh_ipolity2 > 5)
# Compare filtered to unfiltered (just checking some
# pertinent columns)
qog[1:10, c(2, 8:10)]## cname aid_cpsc aid_crnc aid_crnio
## 1 Afghanistan NA 30 3
## 2 Albania NA 20 1
## 3 Algeria NA 20 5
## 4 Andorra NA NA NA
## 5 Angola NA 22 2
## 6 Antigua and Barbuda NA 9 2
## 7 Azerbaijan NA 20 2
## 8 Argentina NA 25 2
## 9 Australia 2606557696 NA NA
## 10 Austria 461456512 NA 1filter[1:10, c(2, 8:10)]## cname aid_cpsc aid_crnc aid_crnio
## 1 Albania NA 20 1
## 2 Andorra NA NA NA
## 3 Antigua and Barbuda NA 9 2
## 4 Argentina NA 25 2
## 5 Australia 2606557696 NA NA
## 6 Austria 461456512 NA 1
## 7 Bahamas NA 1 1
## 8 Bangladesh NA 1 1
## 9 Armenia NA 21 2
## 10 Barbados NA 10 5You can also use the group_by function to group and summarize.
summarise(group_by(qog, ht_region), mean(fh_ipolity2, na.rm = T))## # A tibble: 10 × 2
## ht_region `mean(fh_ipolity2, na.rm = T)`
## <fctr> <dbl>
## 1 1. Eastern Europe and post Soviet Union 6.875867
## 2 2. Latin America 7.487500
## 3 3. North Africa & the Middle East 3.716280
## 4 4. Sub-Saharan Africa 5.147255
## 5 5. Western Europe and North America 9.904270
## 6 6. East Asia 6.500000
## 7 7. South-East Asia 4.892458
## 8 8. South Asia 5.609441
## 9 9. The Pacific 8.516396
## 10 10. The Caribbean 9.288803sample_ allows you to sample random rows.
sample_n(qog, size = 20)count allows you to summarize observations according to a grouping, much like we did with table previously.
count(qog, ht_region)## # A tibble: 10 × 2
## ht_region n
## <fctr> <int>
## 1 1. Eastern Europe and post Soviet Union 28
## 2 2. Latin America 20
## 3 3. North Africa & the Middle East 20
## 4 4. Sub-Saharan Africa 49
## 5 5. Western Europe and North America 27
## 6 6. East Asia 6
## 7 7. South-East Asia 11
## 8 8. South Asia 8
## 9 9. The Pacific 12
## 10 10. The Caribbean 13Pipes
In addition to providing these functions, dplyr inherits piping: a syntax originally found in the magritrr package. Pipes are simply a syntax represented by %>% that allow you to condense code and apply multiple operations at the same time, rather than in separate lines.
For example, suppose you wanted to eliminate observations with high Polity scores from QoG, then group by region, and then find the mean value of fh_ipolity2, as we did previously.
You could construct these in separate lines:
filter <- filter(qog, fh_ipolity2 > 7)
group <- group_by(filter, ht_region)
result <- summarise(group, mean(fh_ipolity2, na.rm = T))
result## # A tibble: 10 × 2
## ht_region `mean(fh_ipolity2, na.rm = T)`
## <fctr> <dbl>
## 1 1. Eastern Europe and post Soviet Union 9.039216
## 2 2. Latin America 8.466667
## 3 3. North Africa & the Middle East 9.083334
## 4 4. Sub-Saharan Africa 8.576987
## 5 5. Western Europe and North America 9.904270
## 6 6. East Asia 9.458334
## 7 7. South-East Asia 7.750000
## 8 8. South Asia 8.500000
## 9 9. The Pacific 9.236342
## 10 10. The Caribbean 9.288803Or, with dplyr and the %>% notation, you could combine them into a single section:
qog %>% filter(fh_ipolity2 > 7) %>% group_by(ht_region) %>% summarise(mean(fh_ipolity2,
na.rm = T))## # A tibble: 10 × 2
## ht_region `mean(fh_ipolity2, na.rm = T)`
## <fctr> <dbl>
## 1 1. Eastern Europe and post Soviet Union 9.039216
## 2 2. Latin America 8.466667
## 3 3. North Africa & the Middle East 9.083334
## 4 4. Sub-Saharan Africa 8.576987
## 5 5. Western Europe and North America 9.904270
## 6 6. East Asia 9.458334
## 7 7. South-East Asia 7.750000
## 8 8. South Asia 8.500000
## 9 9. The Pacific 9.236342
## 10 10. The Caribbean 9.288803You may be thinking: this all sounds great, but what’s a “tibble”…
tidyverse
In addition to many other great R tools, Hadley Wickham has created and maintained a suite of products called the tidyverse that provide an alternative for many of the simple procedures we have seen so far in R. All of these are integrated into a single package (tidyverse) or can be called on their own (e.g., tidyr), but are designed to have a uniform syntax and set of options that makes it easier for users of one set of tools within the tidyverse to learn and transition to another set of tools. The piping as implemented in dplyr for example, might remind you of the + syntax used with ggplot2, which is part of tidyverse.
Loading tidyverse will include several packages, in particular ggplot2, tibble, tidyr, readr, purrr and dplyr. tidyverse also encapsulates stringr and forcats, which we will cover in this segment. (As an aside, the cat jokes are strong here. Check out #Rcatladies on twitter for more.) You can read full descriptions of all the components here.
So to answer that question from above, “tibbles” are the tidyverse alternative to dataframes in base R. Note that while some tidyverse functions will take dataframes as well, if you need a tibble object you can always use as_tibble(). One of the big differences between a tibble and a dataframe is–as previous lectures have cautioned–data types are not automatically transformed. That is, for example, strings are never automatically converted to factors (more in this in a moment). Another difference is in subsetting behavior. You can still use $ as we have been to select a single variable to view or summarize, but you can also use [[:
data(iris)
iris[["Sepal.Length"]]## [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4
## [18] 5.1 5.7 5.1 5.4 5.1 4.6 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5
## [35] 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8 5.1 4.6 5.3 5.0 7.0
## [52] 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8
## [69] 6.2 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4
## [86] 6.0 6.7 6.3 5.6 5.5 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8
## [103] 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8 6.4 6.5 7.7 7.7
## [120] 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7
## [137] 6.3 6.4 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9If you want to go full-on in the native syntax, you can also use a pipe like so:
iris %>% .$Sepal.Length## [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4
## [18] 5.1 5.7 5.1 5.4 5.1 4.6 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5
## [35] 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8 5.1 4.6 5.3 5.0 7.0
## [52] 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8
## [69] 6.2 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4
## [86] 6.0 6.7 6.3 5.6 5.5 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8
## [103] 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8 6.4 6.5 7.7 7.7
## [120] 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7
## [137] 6.3 6.4 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9# or
iris %>% .[["Sepal.Length"]]## [1] 5.1 4.9 4.7 4.6 5.0 5.4 4.6 5.0 4.4 4.9 5.4 4.8 4.8 4.3 5.8 5.7 5.4
## [18] 5.1 5.7 5.1 5.4 5.1 4.6 5.1 4.8 5.0 5.0 5.2 5.2 4.7 4.8 5.4 5.2 5.5
## [35] 4.9 5.0 5.5 4.9 4.4 5.1 5.0 4.5 4.4 5.0 5.1 4.8 5.1 4.6 5.3 5.0 7.0
## [52] 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8
## [69] 6.2 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4
## [86] 6.0 6.7 6.3 5.6 5.5 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7 6.3 5.8
## [103] 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8 6.4 6.5 7.7 7.7
## [120] 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7
## [137] 6.3 6.4 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9For more on tibbles, check out the relevant chapter from the assigned reading in R for Data Science.
forcats
forcats provides some clean methods for establishing factor levels and variables.
haven
haven provides functions for reading in data from other program formats, such as SAS, SPSS, and STATA. Note that haven’s functionality allows you to read in .dta files from STATA 14, whereas the base function in foreign only reads through STATA 12.
readr
readr is particularly useful for evaluating data stored in spreadsheets/tables.
tidyr
tidyr focuses on the structure of your data. Previous materials have covered the reshape package, and reshape2 provides additional, more intuitive functionality. Both are still available for reshaping your data, as we did in previous weeks, from long to wide format and vice versa. tidyr takes over where these packages leave off. The official vignette provides an excellent overview of how its functionality will be helpful in cleaning your data.
stringr
As you may by now be aware, strings are a very important component of data analysis but can be a very cumbersome one if handled improperly. This is where stringr comes in. In general, stringr allows you to manipulate character strings, as well as to extract meaningful information stored in strings. For example, you might want to subset a vector of strings to only include certain values, or you may want to count how many of a certain value are contained in string variables. We have previously covered how to do this kind of manipulation if you are interested in extracting full character strings (for example, every instance of “Afghanistan”). stringr operates alongside R’s built-in capacity to handle regular expressions in order to perform a more advanced version of this function.
Regular expressions
Regular expressions are a series of characters that you can use to format a search. Some of the underlying concepts will be familiar if you have ever used special search terms on, for example, internet searches or library websites, where " " give you items with an exact match, and * serves as a wildcard (so that you can, for example, search for doc* and receive results with both doc, document, documents, documentary, etc.), etc. Regular expressions invoke this type of logic to allow you to more precisely search for terms or character fragments that you need.
The set of regular expressions is fairly long and complex, but some key characters you should know include:
( ) groups characters
[ ] used for “or” conditions
. can match anything
\ is an escape character, so that you can reference punctuation-type characters directly even when they would otherwise be operators
^ searches from the beginning
$ searches from the end
The strings overview in R for Data Science provides a very comprehensive explanation of all of these and more.
You can access regular expressions outside of the stringr package simply using base R’s grep, grepl, or regexp. Use this guide to see more examples of how the arguments for these functions will allow you to incorporate regular expressions into your work. Note that these are great for textual analysis, but are also great for systematically renaming variables, variable values, or even files. For example, suppose you read in a series of .csv files at once that you need to merge into a single dataset. Each file represents data from a different country and is named something like datasource-countryname-year.csv. You can use regular expressions to split this string by -, select the second element to extract the name of the country from the filename, and use it to populate a column named country in the final, merged dataset.
This cheatsheet provides a comprehensive list of different parts of regular expressions you can use, and this website will actually allow you to test your regular expressions as you create them, to see highlighted parts of text that the expression will correctly identify.
lubridate package
While technically not part of the main tidyverse set of packages, you should be aware of the lubridate package, which uses tidy-esque syntax to make your life easier in evaluating and converting date objects (which, like strings, can cause huge problems if handled incorrectly).
As mentioned in previous weeks, R often interprets dates as strings. This makes the creation of time-related variables (for example, duration variables) more challenging. lubridate lets you specify the format of the dates listed as strings or numbers in your data—whether mdy() or ymd() or dmy()—and it automatically interprets them into date objects. These functions are also extensible to include hour and minute timestamps. You can also use the mutate function in the tidyverse to convert date/time information stored across multiple variables into a single date/timestamp variable using make_datetime. Even more explanation is available in the R for Data Science section.
More advanced applications
All of the tools and techniques we have discussed so far in the course can be applied in combination to tackle even the most technical of projects. In this segment we will examine just a few advanced applications that weave together knowledge from earlier parts of this class, and we will end with some tips on how to continue learning on your own.
Introductory text analysis: R vs. tidyr
There are several approaches one can take to text analysis in R, ranging from basic text mining to more advanced models using LDA, for example. In this section, we will cover some basic text analysis techniques across several packages to illustrate the difference in approach between using base R and dplyr syntax.
As a preliminary step before even worrying about packages and functions, you should organize your textual data to read into R. Typically this means that you want each document in your “corpus” to be a separate .txt file within a single directory that you can easily reference from within R. This may require data cleaning done outside of R, for example if you need to OCR .pdf documents or download files before converting to .txt. Your data may also be stored in a .csv format, however, with columns indicating things like “author”, “date”, and then one including the full “text” for each document in your corpus already.
For now assume that you have multiple files, each with text in them. Once you have your text in the right place, you will want to load this as a directory:
# For mac
mydirectory <- file.path("~", "filelocation", "mytexts")
# For PC
mydirectory <- file.path("C:", "filelocation", "mytexts")
# Look at the separate files you loaded
dir(mydirectory)tm
The next steps will illustrate text analysis techniques with standard R syntax using the tm package. The first step here is to convert your text into a “Corpus”-type object.
# If all documents are separate .txt
mydocs <- Corpus(DirSource(mydirectory))
# If text is a column in a dataframe, read in from a .csv for
# example
mydocs <- Corpus(VectorSource(dataframename$text))Note here if you’re reading the tm documentation that Corpus is an alias for VCorpus (a volatile corpus)—either may be used. Also note that when reading in your text data, the default is for documents to be read in English. If your data are in another language, you will have to specify it. For example, if your documents are in French:
mydocs <- Corpus(DirSource(mydirectory), readerControl = list(language = "fre"))The language options are provided with IETF language codes and adopted through the NLP package. You can read more about IETF codes here.
Once you have a corpus object, you are ready to clean your text for analysis. The steps to do this with tm are more fully explained here, but I will highlight the general procedure using a built-in corpus accessible through the tm package, acq, which contains 50 news articles about corporate acquisitions from Reuters.
In general, you will want to:
- Remove characters that interfere with the ability of the program to identify words, including punctuation, numbers and symbols. Here is where you can use your new regex skills to locate specific symbols to exclude as well.
data(acq)
mydocs <- tm_map(acq, removePunctuation)
mydocs <- tm_map(mydocs, removeNumbers)# Exclude the # sign, replace it with nothing
for(i in seq(mydocs))
{
docs[[i]] <- gsub("\#", "", docs[[i]])
}- By convention, change all words to lowercase so that every instance is counted together
mydocs <- tm_map(mydocs, tolower)- Eliminate “stopwords,” which are common words that, while occurring frequently in texts, do not contribute to the meaning (in English, for example “a” and “the” are stopwords). Here, you will need to identify what dictionary
tmshould use to find stopwords. If your data are not in English, use the appropriate stopword dictionary instead, but note that the options for identifying stopwords are available in fewer languages.
mydocs <- tm_map(mydocs, removeWords, stopwords("english"))- Next, stem the words in your text. This again relies on the underlying language of your text, but essentially will remove endings so that when you identify word counts or sentiments, it will treat, for example, “document” and “documents” as instances of the same essential word. This function relies on the
SnowballCpackage.
library(SnowballC)
mydocs <- tm_map(mydocs, stemDocument)- Now fix all the extra white space we created by eliminating text.
mydocs <- tm_map(mydocs, stripWhitespace)- Finally, resave your data as plain text to insure easy analysis.
mydocs <- tm_map(mydocs, PlainTextDocument)Once these preliminary steps are completed, your next general task is to create a DocumentTermMatrix object, which will allow you to “quantify” the text in terms of word counts. For idiosyncratic reasons, you may need to create and immediately transpose this matrix–it depends on what information you want. This matrix will also serve as a basis for evaluating differences across your documents, if you like.
mydocsDTM <- DocumentTermMatrix(mydocs)
# If you want to evaluate differences among documents in
# terms of distance:
mydocsDTM.mat <- as.matrix(mydocsDTM)# You can add back in rownames if you create an object
# listing the names in your text directory
names <- list.files("filelocation")
rownames(mydocsDTM.mat) <- namesIn our case, the original articles in acq just came with numbers, which we can add back in accordingly:
rownames(mydocsDTM.mat) <- names(acq)library(proxy)##
## Attaching package: 'proxy'## The following objects are masked from 'package:stats':
##
## as.dist, dist## The following object is masked from 'package:base':
##
## as.matrix# Note that to evaluate dissimilarity you will need to
# specify the method for calculating distance--check the
# documentation
dissim <- dist(mydocsDTM.mat, method = "cosine")
# You can also use this information to create a cluster plot
# with a built in plot function, again choosing a method. You
# could also add names to this manually with labels=names if
# you are not using built-in data
hcluster <- hclust(dissim, method = "ward.D")
plot(hcluster)
# You may also care about word frequency
freq <- colSums(mydocsDTM.mat)
ordered <- order(freq)
freq[head(ordered)]## ability acceptance accompanied acquiring acres acting
## 1 1 1 1 1 1freq2 <- sort(colSums(mydocsDTM.mat), decreasing = T)
head(freq2, 10)## said dlrs pct mln company inc shares reut stock
## 186 100 70 65 63 53 52 50 46
## will
## 35See the vignette linked above for how to graph word frequency using ggplot2. You can also restrict attention only to words that have frequency above a certain threshhold.
You may also be interested in the association between particular words, which is a beginning step toward topic analysis.
# Fill in words in the second argument, set a lower bound for
# correlation
findAssocs(mydocsDTM, "good", corlimit = 0.75)## $good
## ability accompanied activity along analytical
## 1.00 1.00 1.00 1.00 1.00
## arena believes better big bilion
## 1.00 1.00 1.00 1.00 1.00
## book boost boosting calculated capitalization
## 1.00 1.00 1.00 1.00 1.00
## center command comments components constitute
## 1.00 1.00 1.00 1.00 1.00
## contrary contributing disagreed divestitures dividend
## 1.00 1.00 1.00 1.00 1.00
## elevating endeavors enormous exposure fact
## 1.00 1.00 1.00 1.00 1.00
## fetched fuel globally goes going
## 1.00 1.00 1.00 1.00 1.00
## heavily heavy helped investing jeffery
## 1.00 1.00 1.00 1.00 1.00
## kinds lesser lewis lift lipper
## 1.00 1.00 1.00 1.00 1.00
## mark might others outlined partially
## 1.00 1.00 1.00 1.00 1.00
## percentage performance perrin preclude premium
## 1.00 1.00 1.00 1.00 1.00
## probably questioned realistic reduced road
## 1.00 1.00 1.00 1.00 1.00
## rumor split suggests takes talk
## 1.00 1.00 1.00 1.00 1.00
## thereby theyre theyve travel true
## 1.00 1.00 1.00 1.00 1.00
## tuesday undervalued unhappy valuations variation
## 1.00 1.00 1.00 1.00 1.00
## volume water yesterday reflect market
## 1.00 1.00 1.00 0.95 0.91
## eckenfelder partial place strong operating
## 0.89 0.89 0.89 0.89 0.87
## spinoff higher make since profitable
## 0.85 0.81 0.81 0.81 0.75If you’re interested in more complex topic model designs, check out the topicmodels package. There are other packages like tm that you might want to explore as well, including qdap. See this vignette for more details.
tidytext
Much of the syntax in the above functions will be familiar, even if the operations are not. The sort, head, and plot commands, as well as optional arguments and the regular expressions component weave in techniques that you have already learned. Yet as with many tasks in R, there is more than one way to perform text analysis. Many additional packages and functions build on tm, but another option is to approach text analysis through the language of dplyr using the tidytext package.
Here again we will use a built-in dataset, this time with 2,246 Associated Press articles from the topicmodels package, to walk through the major functions in the package. The same steps apply as above, however, if your data are not already in a corpus format.
library(topicmodels)
data("AssociatedPress", package = "topicmodels")
AssociatedPress## <<DocumentTermMatrix (documents: 2246, terms: 10473)>>
## Non-/sparse entries: 302031/23220327
## Sparsity : 99%
## Maximal term length: 18
## Weighting : term frequency (tf)You can “clean up” these data in a single line of code. Because the AssociatedPress data are already in a DocumentTermMatrix format, you don’t need to separate out by word, but you could do that with unnest_tokens if your data are not structured this way. Note that unnest_tokens also allows you to change the ngram for your analysis–that is, how many words or word segments are evaluated by comparison. For example if your ngram is 2, you are evaluating 2-word phrases.
tidied <- tidy(AssociatedPress)
tidied## # A tibble: 302,031 × 3
## document term count
## <int> <chr> <dbl>
## 1 1 adding 1
## 2 1 adult 2
## 3 1 ago 1
## 4 1 alcohol 1
## 5 1 allegedly 1
## 6 1 allen 1
## 7 1 apparently 2
## 8 1 appeared 1
## 9 1 arrested 1
## 10 1 assault 1
## # ... with 302,021 more rowsNote that using grep and other subsetting techniques we have explored, you can also narrow down what documents you want to look at within your broader set of texts. For example, if you only wanted to include articles that mentioned government actions or actors, you might want to select based on whether the document includes words like “government” or “president” before you tidy your text.
selection <- "government|president"
final.docs <- mydocs[grepl(selection, mydocs$text), ]You can perform similar tasks to the ones we did above with tm using dplyr syntax, such as removing stopwords, stemming, and obtaining word counts:
data("stop_words")
cleaned <- tidied %>% anti_join(stop_words, by = c(term = "word"))
cleaner <- cleaned %>% mutate(word = wordStem(term))
cleaner %>% count(word, sort = TRUE)## # A tibble: 6,722 × 2
## word n
## <chr> <int>
## 1 report 1144
## 2 nation 1061
## 3 offici 935
## 4 time 878
## 5 peopl 872
## 6 presid 870
## 7 dai 833
## 8 govern 773
## 9 call 760
## 10 week 748
## # ... with 6,712 more rowsWith these cleaned data, and with similar join syntax, we can also do other operations, such as sentiment analysis. You have several sentiment lexicons to choose from, but for this example we will simply use bing.
ap_sentiments <- cleaner %>% inner_join(get_sentiments("bing"),
by = c(term = "word"))
ap_sentiments## # A tibble: 27,068 × 5
## document term count word sentiment
## <int> <chr> <dbl> <chr> <chr>
## 1 1689 bug 5 bug negative
## 2 2141 bug 1 bug negative
## 3 1489 disarray 1 disarrai negative
## 4 1507 disarray 1 disarrai negative
## 5 1631 disarray 1 disarrai negative
## 6 2013 disarray 1 disarrai negative
## 7 2101 disarray 2 disarrai negative
## 8 2219 disarray 1 disarrai negative
## 9 1287 uphold 1 uphold positive
## 10 1394 uphold 1 uphold positive
## # ... with 27,058 more rowsCheck out this vignette for more on tidytext techniques, as well as this guide, featuring Jane Austen novels, for a more in-depth look at how to handle different corpus structures and for discussions of plotting sentiment.
More, different syntaxes: R as a front-end
Besides the divide in syntax between base R and the tidyverse, R can also be used in more advanced applications as a front-end to other software and languages.
rPython
The rPython package allows you to run Python code directly through R. As is common with many adaptations of other languages or programs to R, the functions in this program will have syntax from R that runs code from another language or program. This is useful if you need a particular functionality of one program or language, but are most comfortable or familiar with R’s interface and syntax.
RSelenium
RSelenium is a great example of this. If you are not familiar with Python syntax for scraping webpages, for example, or have a more complex scraping project that requires you to manipulate a website at the same time, you can invoke Selenium (a java-based program) through R to do the necessary tasks. Note that to do this, you may want to run a Docker container. Instructions on how to do that are in the linked vignette above.
Rcpp
The Rcpp package allows you to leverage both R and C++ syntaxes and libraries. This can be particularly useful as you optimize computationally intensive model estimation, or for creating your own packages. Read through the explanation in Advanced R for more.
BUGS, JAGS, and Stan
Several different tools exist for Bayesian analysis in R, many requiring external programs. BUGS/WinBUGS/OpenBugs are a set of classic programs for several platforms with its own set of GUIs. If you visit the BUGS site, you’ll notice this warning:
“The programs are reasonably easy to use and come with a wide range of examples. There is, however, a need for caution. A knowledge of Bayesian statistics is assumed, including recognition of the potential importance of prior distributions, and MCMC is inherently less robust than analytic statistical methods. There is no in-built protection against misuse.”
While knowing your data and knowing your method are key in all types of statistical computing practice, it bears mentioning here that relying on packages or programs that automate elements of your analysis also requires you to recognize their functions, limitations, and underlying assumptions.
Note about NIMBLE: NIMBLE is another option for Bayesian analysis that builds on BUGS and allows you to write your own algorithms if you choose. Handle with care.
JAGS works much the same as BUGS and follows similar coding conventions, and at this stage is likely most commonly used. If you are just getting started with JAGS, you will notice that it requires data in a list format, a full model specification, and instructions about the sampling process (burn-in, number of iterations, etc.).
A basic model will take the general form
model{
# where n is number of observations in data
for(i in 1:n){
# each value of the dependent variable will be determined by some density function (e.g., dbern if it is bernoulli) with a given parameter (e.g., for shape or scale of that distribution)
y[i] ~ densityfn(param)
}
# parameter value of the distribution for y has a prior density of its own, with a parameter or set of parameters as well (e.g., dunif(0,1))
param ~ paramdensity(param2)
}You can specify this model within R itself, or save it out as a .bugs file that you later call with JAGS:
modelString = "model{
y[i] ~ densityfn(param)
}
param ~ paramdensity(param2)
}"
writeLines(modelString, con = "mymodel.bug")You can then call this model directly from a .bug file (or if you saved it as a .jag file) in your working directory when running JAGS. Here I’m using the run.jags function from the runjags package. The arguments are effectively the same if you use rjags, as I do below.
post <- run.jags("mymodel.bug",
# what values am I interested in for results
monitor="param",
# data must be in the form of a list object: if you have binary or factor explanatory variables, you will want to use model.matrix and specify which variables are "contr.treatment"
data=list(y=y, n=length(y)),
# how large of a sample do you want? the bigger it is, the longer it will take, but complex models may take longer to converge
sample=10000,
# throw away the first several samples with poor mixing
burnin=1000,
# take only a subset of total samples: saves on memory and smooths results
thin=10,
# you can technically run more than one "chain" of this analysis or spread the sampling over cores (parallelization), but you don't need to
n.chains=1)In terms of preparing your data, as stated above, if you have factor variables for example, you will need to use model.matrix to expand your grid (that is, to create rows and columns for values of each variable for each level of your factors). Factor variables are “contrasts.” For example, if you had a model predicting voting in terms of age, sex, and party, your data preparation might look like this:
y <- vote.variable
x.expand <- model.matrix(~age + sex + party, data = dataname,
contrasts = list(sex = "contr.treatment", party = "contr.treatment"))
# for modeling, you may need to remove the intercept
x.expand <- x.expand[-1]Note that age is not a contrast because it is presumed to be continuous.
Using this setup for your data, you want to estimate a logit (voted y/n), so your JAGS model might look like:
model{
for (i in 1:n){
y[i] ~ dbern(p[i])
y.hat[i] <- b0 + b1*x1[i] + b2*x2[i] + b3*x3[i]
p[i] <- 1/(1+ exp(-y.hat[i]))
}
b0 ~ dnorm(0,1.0e-4)
b1 ~ dnorm(0,1.0e-4)
b2 ~ dnorm(0,1.0e-4)
b3 ~ dnorm(0,1.0e-4)
}In this setup you have one b for each of your covariates, you are looping over each observation in you independent variables, and you have the p[i] allowing you to estimate the logit. Choices of priors are up to you–these ones are all the same, but note that you would specify one for each coefficient.
Alternatively, you could collapse some of this with inprod. Here you will note I have changed to using jags.model and related functions from the rjags package. Either syntax or setup will work.
model{
for (i in 1:n){
y[i] ~ dbern(p[i])
y.hat[i] <- b0 + inprod(b[]*x[i,])
p[i] <- 1/(1+ exp(-y.hat[i]))
}
b0 ~ dnorm(0,1.0e-2)
for(j in 1:k){
b[j] ~ dnorm(0,1.0e-4)
}
}
expanded.model<-jags.model("mymodel.bug", data=list(y=y, x=x.expand, k=3, n=length(y)))
update(expanded.model, 10000)
expanded.model.j <- coda.samples(expanded.model, c("b0","b","p"), n.iter=100000, thin=10)
summary(expanded.model.j)In this case, you have expanded a grid similar to what you did with model.matrix, so that the set of coefficients you have tracks the number of explanatory variables you include. Note that if you do this, you will want to loop over your b’s to provide priors for each, assuming each gets the same prior. Also note that in creating the second loop we have said “loop over each j element for values 1:k,” meaning that you will now need to specify k (the number of explanatory variables in your model) when you provide your data list.
Finally, note that you could also skip the messy business about y.hat and use the logit() function:
model{
for (i in 1:n){
y[i] ~ dbern(p[i])
logit(p[i]) <- b0 + inprod(b[]*x[i,])
}
b0 ~ dnorm(0,1.0e-2)
for(j in 1:k){
b[j] ~ dnorm(0,1.0e-4)
}For more example code and some video lectures and slides, check out these resources:
This Github repo with rjags scripts, simulated data, and .bug files for various models
Besides using BUGS and JAGS, another option available to you is Stan, a project headed by Andrew Gelman. Although you can access it through R with rstan, it features syntax much closer to C++. It can be more efficient for particularly complex hierarchical models and includes the No-U-Turn Sampler (NUTS), an optimization of Hamiltonian Monte Carlo. Stan also has a great set of documentation, including a language manual, examples, and tutorials to help you get started and to compare syntax with BUGS. Stan also has a user community Google group to which you can post issues and questions.
Learning new packages
As you begin to explore R more on your own, you will need to adapt to different coding styles and packages in order to do different types of analyses. Learning how to learn is also a key part of the educational process (and certainly a staple of graduate school). Below are some quick guides on how to teach yourself new techniques in R, with examples from survey and Baysian analyses.
Surveys
Suppose that as a new user of R (useR), you are confronted with a survey analysis project. Taking the tools you have already learned in this class, you feel comfortable reading in survey data and doing basic summary analysis and visualization. When it comes to more complex analyses, however, you need guidance on survey weights and how to handle them, because your results will depend on the weights you choose and the type of sample you have.
- Step 1: Google
If you google “survey weights in R” or “survey analysis R” you will see that a package called survey has documentation in the first few hits. You will also notice a few posts on R-bloggers and similar sites that provide explainers related to the package, situated in broader context about survey analysis.
Having done this initial search, you have several components of the information you need. First, you discover that the survey package can help you in creating a design object that you use in later analysis using svydesign(). With this function, you can directly specify the name of the weights variable if it already exists in the data. You will also specify the stratification of your sample (or samples).
Once you have created a design object, you will use the built-in functions in survey to do your basic analysis. These functions will take you far if you are manipulating data collected by other researchers. But what if the survey data are your own?
Here you may need to calculate weights yourself. To do this, you will need to use the rake() function in survey, which requires both sample and population margins.
Suppose you had voting data for individuals and want to compute weights by party (R, D, I) and income bracket (High, Medium, Low). Take note of the information that blog posts and the package documentation take into account. You can do this primarily by investigating the arguments of functions and the type of information they take.
In this case, both of these variables of interest are factor variables with 3 levels each and you will need relative frequencies for each level. Following the example in this R-bloggers post, you would create dataframes relevant to population characteristics accordingly:
distr.party <- data.frame(party = c("R", "D", "I"), Freq = nrow(mydata) *
c(0.45, 0.5, 0.05))
distr.income <- data.frame(income = c("High", "Medium", "Low"),
Freq = nrow(mydata) * c(0.01, 0.6, 0.39))Now you will use rake() to relate these values to your sample.
raked.leaves <- rake(design = mydata.nowt, sample.margins = list(~party,
~income), population.margins = list(party.dist, income.dist))Unlike the implementation of packages in other programs (cough STATA), you may need to examine and adjust these weights to be within reasonable bounds. This again requires knowing your data.
Note that survey has its own syntax for estimating models using survey design objects, for example svyglm instead of regular glm. Look at an example from the package documentation:
data(api)
dclus2 <- svydesign(id = ~dnum + snum, weights = ~pw, data = apiclus2)
model0 <- svyglm(I(sch.wide == "Yes") ~ ell + meals + mobility,
design = dclus2, family = quasibinomial())
summary(model0)##
## Call:
## svyglm(formula = I(sch.wide == "Yes") ~ ell + meals + mobility,
## design = dclus2, family = quasibinomial())
##
## Survey design:
## svydesign(id = ~dnum + snum, weights = ~pw, data = apiclus2)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.82118 0.76726 1.070 0.2916
## ell -0.05462 0.02220 -2.461 0.0188 *
## meals 0.02873 0.01820 1.579 0.1232
## mobility 0.01460 0.03017 0.484 0.6313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.9856512)
##
## Number of Fisher Scoring iterations: 4You will notice that the syntax looks very similar to what you are familiar with from our discussion of glm: you specify a dependent variable ~ independent variables, the design object where previously you specified data, and the family of model you would like to estimate (in this case a quasibinomial).
You may not be familiar with the I() notation used to specify the dependent variable. Do not panic. Recall our previous discussion that you can modify variables in place within model commands in R, just as we did with log(y). If you search for I() function in R you will easily find documentation to inform you that I() is an inhibiting function that simply “protects” dataframes so that strings are not converted to factors, for example.
- Step 2: More specific searches
Now that you have some basic information about the relevant package and how it functions, you may have more specific questions that pertain to your data. Perhaps you need more information about cluster bootstrapping for resampling. These are questions you can take to Google again, but you might also benefit from taking them directly to specialty sites like Stack Overflow. In doing this, you might find Thomas Lumley’s (the survey package author) page, and will likely find his slides covering more complex applications of the package than are covered in the R-bloggers posts. You will also discover there is a book, “Complex Surveys: A Guide to Analysis Using R” that Lumley wrote that you could reference. [ Note: Because R is open-source, be wary of paying large sums of money for books on its functions or packages—many resources exist for free!]
- Step 3: Iterate
Now that you have completed basic analysis and found a few preliminary answers to your initial tasks and questions, iterate the search process. Try sample code from one site, use examples from a vignette related to a package, and google any errors you receive and read about other users who have had those issues on Stack Overflow. At this point, if not sooner, you may go beyond the lonely internet life to reach out to friends, colleagues, and professors who may have done similar work.
If no one is available locally to help, and answers on the internet are sparse, but you see a technique referenced in a paper, you may also be able to directly find code for your problem on an author’s personal webpage or Github repository. You can also perform a general search of Github repos to look for code directly (try searching cluster robust standard errors).
- Step 4: DIY
Step 4 won’t be for everyone, but in a pinch you can often write your own functions to solve problems that are not well-managed by existing packages. Some examples might include things like cluster robust standard error calculations for complex, rare model types; uncommon types of cross-validation; particular estimators or samplers. For many of these, the internet will still be your best friend as you are (often) converting math to code. Many times, investing in something of your own design is not necessary because current packages have adequate functionality, but if you find yourself doing the same type of analysis routinely and needing something beyond the prefab, coding your own solution can often be efficient in the long-run. This is how packages are introduced in the first place.
Bayesian estimation
The survey example seems relatively straightforward: identify the relevant package for your needs and use internet resources to help you implement its functions. But what if there are multiple such packages: how do you choose?
In addition to complementary programs like JAGS and stan that can be used with an R interface, there are several packages designed to do Bayesian estimation without needing to fully mathematically specify a model. If you are familiar with Bayes, you can imagine why this creates more need than ever to truly understand what kinds of assumptions a package operates under and what functions go on under the hood.
Here again you may begin by simply googling "Bayesian packages R" or a specific task in Bayesian analysis you want to complete. If you are already aware multiple packages exist, however, (or become aware and overwhelmed quickly by googling), another option is to check out the Cran Task View. The Task View does not list many categories, but for each included category it provides an overview of common packages and their relative advantages. In the Task View for Bayes, you will see a host of packages specifically designed for Bayesian analysis, as well as a series of packages that include functions for particular modeling challenges. If you have a specific task in mind, you can easily search this page for mentions; otherwise, everything is in alphabetical order. Often if you are looking to do Bayesian analysis in a canned routine through one of these packages, your problem is already well-defined by the type of model you are trying to estimate or the type of data you have at your dispoal—these will also help you to narrow the set of possible candidates.
While you should always prioritize the package that has functions most suited to your needs, having a mental catalog of general functions that cover your area of interest will help you to identify possible places to start in the future as new analysis tasks arrive. This section will highlight just a few useful packages to keep in mind.
The arm package is a great general package useful for both MLE and Bayes, with code that accompanies Gelman and Hill (2007). Specifically relevant to Bayes,
armincludes thebayesglm()function. When investigating this function, take note of its default behaviors. For example, prior means default to 0 if not set manually. If you specify only one, but have several predictors (independent variables excluding intercept), it is assumed that that single prior mean is the same for all predictors. Otherwise, you need a vector of prior means of length equal to the number of predictors.MCMCpack contains a number of model types and postestimation analysis functions for Bayes. Andrew Martin and Kevin Quinn provide a set of slides explaining its functions and relative advantage versus using something like BUGS. MCMCpack uses
codaas a basis for posterior analysis, so you will need to use extra documentation to understand its full range of functionality. It does include built-in functions for calculating Bayes Factors and posterior model probability. Also note that while functions likeMCMCregressdo not require you to specify priors in quite the same way as JAGS, say, and that might be more comfortable, it does require you to think about how you want to calculate the marginal likelihood (using the argumentmarginal.likelihood=) if you want to deviate from the default.MCMCglmm includes both standard glm’s estimated with MCMC as well as some rarer types that you might be interested in (e.g., zero-inflated and multi-trait models). While similar to lme4 in several ways, MCMCglmm also handles situations where packages like lme4 fail to converge.
lme4 is another general package for glm’s, in particular emphasizing multilevel models, and several Bayes-specific packages require it as a dependency or extend on its functions. This vignette gives a good overview of the differences in lme4’s syntax and approach relative to base glm specifications.
In addition to useful functions for hurdle and zero-inflated models in MLE, pscl also includes Bayesian IRT.
MNP is specifically designed for the implementation of multinomial probit models.
Bolstad provides the
bayes.lmfunction for fitting Bayesian linear models.blme contains the functions
blmerandbglmerfor doing linear or glm mixed-effect models. The syntax here is much the same as you would find inlm,lmer, orglmerbut the functions inblmeallow you to conduct Bayesian or penalized MLE analyses.brms is a newer package building on Stan that uses syntax similar to lme4. It includes less common models, such as zero-inflated or survival models. It also allows you to do posterior prediction and crossvalidation.
A Note about Venue-Shopping
In the coming weeks, you will learn many of the same techniques we have discussed in R for STATA. Learning both will give you the option to select whichever is most comfortable and convenient for you going forward, but note that their functionality on some dimensions may differ. Some methods or routines have more developed packages for STATA than for R and vice versa, often as a function of departmental conventions (where, e.g., sociologists may be more apt to use STATA). If you are comfortable with the syntax of both and discover that there are not many options for the analysis you want to conduct in R, you may also want to explore STATA code for similar functions. For example, code for survey analysis and for time series analysis is very well developed in STATA, whereas the R packages associated with each may require you to do many of your own calculations or much more data transformation than you would prefer.