R is now a large community. Today we have 6,139 R packages. I am planning here to jot down the code snippets based on the newer packages (by following table of available packages, sorted by date of publication). The chunk of the codes will be taken either from the Vignette or from the main package documentation pdf. I also like to include stack overflow threads if there are interesting questions. These interesting and easy-to-apply codes will be applied to different data sets to develop models and assumptions. Here’s Part 2.
Five Packages used in Part 2. [Why 5? Here’s why.]
BCDating, historydata, USAboundaries, verification, dygraphs
Code chunk for Business Cycle Dating and Plotting Tools.
library(BCDating)
data("Iran.non.Oil.GDP.Cycle")
dat <- BBQ(Iran.non.Oil.GDP.Cycle, name="Dating Business Cycles of Iran")
show(dat)
## Peaks Troughs Duration
## 1 <NA> 1368Q2 <NA>
## 2 1370Q4 1372Q4 8
## 3 1373Q2 1374Q1 3
## 4 1376Q1 1378Q1 8
## 5 1378Q4 1380Q1 5
## 6 1381Q1 1381Q3 2
## 7 1383Q1 1383Q4 3
## 8 1386Q1 1387Q2 5
## 9 1387Q4 1388Q2 2
## 10 1389Q4 <NA> <NA>
summary(dat)
## Phase ]Start ;End] Duration LevStart LevEnd Amplitude
## 1 Recession <NA> 1368Q2 NA NA 0 NA
## 2 Expansion 1368Q2 1370Q4 10 0 0 0.1
## 3 Recession 1370Q4 1372Q4 8 0 0 0.1
## 4 Expansion 1372Q4 1373Q2 2 0 0 0.0
## 5 Recession 1373Q2 1374Q1 3 0 0 0.0
## 6 Expansion 1374Q1 1376Q1 8 0 0 0.1
## 7 Recession 1376Q1 1378Q1 8 0 0 0.0
## 8 Expansion 1378Q1 1378Q4 3 0 0 0.0
## 9 Recession 1378Q4 1380Q1 5 0 0 0.0
## 10 Expansion 1380Q1 1381Q1 4 0 0 0.0
## 11 Recession 1381Q1 1381Q3 2 0 0 0.0
## 12 Expansion 1381Q3 1383Q1 6 0 0 0.0
## 13 Recession 1383Q1 1383Q4 3 0 0 0.0
## 14 Expansion 1383Q4 1386Q1 9 0 0 0.0
## 15 Recession 1386Q1 1387Q2 5 0 0 0.0
## 16 Expansion 1387Q2 1387Q4 2 0 0 0.0
## 17 Recession 1387Q4 1388Q2 2 0 0 0.0
## 18 Expansion 1388Q2 1389Q4 6 0 0 0.0
## 19 Recession 1389Q4 <NA> NA 0 NA NA
##
## Amplitude Duration
## Exp=]T;P] 0 5.6
## Rec=]P;T] 0 4.5
plot(dat)
plot(dat,Iran.non.Oil.GDP.Cycle)
data("MBRI.Iran.Dating")
plot(MBRI.Iran.Dating)
data("Iran.non.Oil.GDP.Quarterly.Growth")
data("MBRI.Iran.Dating")
avggrowth <- avgts(Iran.non.Oil.GDP.Quarterly.Growth,MBRI.Iran.Dating)
## Warning in window.default(x, ...): 'start' value not changed
cbind(avggrowth,Iran.non.Oil.GDP.Quarterly.Growth)
## avggrowth Iran.non.Oil.GDP.Quarterly.Growth
## 1367 Q2 0.1646207 -4.07094167
## 1367 Q3 0.1646207 2.79743640
## 1367 Q4 0.1646207 1.74877726
## 1368 Q1 0.1646207 -1.56021883
## 1368 Q2 0.1646207 1.90805046
## 1368 Q3 3.0723742 2.83623375
## 1368 Q4 3.0723742 10.88055401
## 1369 Q1 3.0723742 -1.33164518
## 1369 Q2 3.0723742 1.91343779
## 1369 Q3 3.0723742 2.41133835
## 1369 Q4 3.0723742 6.33867468
## 1370 Q1 3.0723742 2.22078677
## 1370 Q2 3.0723742 -0.80004894
## 1370 Q3 3.0723742 7.33667195
## 1370 Q4 3.0723742 -1.08226147
## 1371 Q1 -0.1249312 2.01249577
## 1371 Q2 -0.1249312 0.74387567
## 1371 Q3 -0.1249312 -1.03329087
## 1371 Q4 -0.1249312 0.36154040
## 1372 Q1 -0.1249312 4.62375750
## 1372 Q2 -0.1249312 -1.18687451
## 1372 Q3 -0.1249312 -5.21039585
## 1372 Q4 -0.1249312 -1.31055768
## 1373 Q1 4.6276720 3.42239895
## 1373 Q2 4.6276720 5.83294504
## 1373 Q3 -1.6406157 -3.16349076
## 1373 Q4 -1.6406157 -1.28602569
## 1374 Q1 -1.6406157 -0.47233059
## 1374 Q2 2.0312820 4.10462871
## 1374 Q3 2.0312820 2.40468580
## 1374 Q4 2.0312820 2.43551229
## 1375 Q1 2.0312820 0.08929337
## 1375 Q2 2.0312820 1.30855365
## 1375 Q3 2.0312820 2.23719429
## 1375 Q4 2.0312820 3.15036554
## 1376 Q1 2.0312820 0.52002259
## 1376 Q2 0.3255956 0.33111177
## 1376 Q3 0.3255956 -0.08497765
## 1376 Q4 0.3255956 -0.76753897
## 1377 Q1 0.3255956 2.11167719
## 1377 Q2 0.3255956 1.04099629
## 1377 Q3 0.3255956 0.22974797
## 1377 Q4 0.3255956 2.09536052
## 1378 Q1 0.3255956 -2.35161252
## 1378 Q2 2.1903674 -0.22429148
## 1378 Q3 2.1903674 7.38287469
## 1378 Q4 2.1903674 -0.58748106
## 1379 Q1 0.7194218 -1.00336564
## 1379 Q2 0.7194218 0.33011277
## 1379 Q3 0.7194218 3.67630340
## 1379 Q4 0.7194218 1.63624460
## 1380 Q1 0.7194218 -1.04218596
## 1380 Q2 2.8352818 1.96124296
## 1380 Q3 2.8352818 3.26212598
## 1380 Q4 2.8352818 1.17832025
## 1381 Q1 2.8352818 4.93943799
## 1381 Q2 -0.3773534 0.34760895
## 1381 Q3 -0.3773534 -1.10231567
## 1381 Q4 2.1861083 3.44773355
## 1382 Q1 2.1861083 2.87100193
## 1382 Q2 2.1861083 0.55607496
## 1382 Q3 2.1861083 1.35536701
## 1382 Q4 2.1861083 2.11709022
## 1383 Q1 2.1861083 2.76938187
## 1383 Q2 1.1141754 2.72901031
## 1383 Q3 1.1141754 -3.20785553
## 1383 Q4 1.1141754 3.82137128
## 1384 Q1 1.9668678 2.67397255
## 1384 Q2 1.9668678 0.79258952
## 1384 Q3 1.9668678 0.47202507
## 1384 Q4 1.9668678 8.70496531
## 1385 Q1 1.9668678 -1.20380779
## 1385 Q2 1.9668678 -0.18281883
## 1385 Q3 1.9668678 2.12444450
## 1385 Q4 1.9668678 1.64092948
## 1386 Q1 1.9668678 2.67951004
## 1386 Q2 0.5298574 2.28162019
## 1386 Q3 0.5298574 -1.87405660
## 1386 Q4 0.5298574 -1.50965346
## 1387 Q1 0.5298574 1.04873247
## 1387 Q2 0.5298574 -0.27237525
## 1387 Q3 0.5298574 3.34155637
## 1387 Q4 0.5298574 0.72535061
## 1388 Q1 0.5298574 0.00867241
## 1388 Q2 0.5298574 1.01887019
## 1388 Q3 1.3967372 0.30532549
## 1388 Q4 1.3967372 2.50796151
## 1389 Q1 1.3967372 1.72480749
## 1389 Q2 1.3967372 1.31545395
## 1389 Q3 1.3967372 1.58555484
## 1389 Q4 1.3967372 0.94132021
## 1390 Q1 0.3053066 0.08509787
## 1390 Q2 0.3053066 2.86299079
## 1390 Q3 0.3053066 -1.77010274
## 1390 Q4 0.3053066 0.04324058
plot(MBRI.Iran.Dating,avggrowth)
plot(MBRI.Iran.Dating,Iran.non.Oil.GDP.Quarterly.Growth,averages=TRUE)
## Warning in .local(x, y, ...): BCDating: Plotting Averages only in windowed
## mode
## Warning in window.default(x, ...): 'start' value not changed
## Warning in window.default(x, ...): 'start' value not changed
## Warning in window.default(x, ...): 'start' value not changed
data("Iran.non.Oil.GDP.Cycle")
dat <- BBQ(Iran.non.Oil.GDP.Cycle, name="Dating Business Cycles of Iran")
show(dat)
## Peaks Troughs Duration
## 1 <NA> 1368Q2 <NA>
## 2 1370Q4 1372Q4 8
## 3 1373Q2 1374Q1 3
## 4 1376Q1 1378Q1 8
## 5 1378Q4 1380Q1 5
## 6 1381Q1 1381Q3 2
## 7 1383Q1 1383Q4 3
## 8 1386Q1 1387Q2 5
## 9 1387Q4 1388Q2 2
## 10 1389Q4 <NA> <NA>
summary(dat)
## Phase ]Start ;End] Duration LevStart LevEnd Amplitude
## 1 Recession <NA> 1368Q2 NA NA 0 NA
## 2 Expansion 1368Q2 1370Q4 10 0 0 0.1
## 3 Recession 1370Q4 1372Q4 8 0 0 0.1
## 4 Expansion 1372Q4 1373Q2 2 0 0 0.0
## 5 Recession 1373Q2 1374Q1 3 0 0 0.0
## 6 Expansion 1374Q1 1376Q1 8 0 0 0.1
## 7 Recession 1376Q1 1378Q1 8 0 0 0.0
## 8 Expansion 1378Q1 1378Q4 3 0 0 0.0
## 9 Recession 1378Q4 1380Q1 5 0 0 0.0
## 10 Expansion 1380Q1 1381Q1 4 0 0 0.0
## 11 Recession 1381Q1 1381Q3 2 0 0 0.0
## 12 Expansion 1381Q3 1383Q1 6 0 0 0.0
## 13 Recession 1383Q1 1383Q4 3 0 0 0.0
## 14 Expansion 1383Q4 1386Q1 9 0 0 0.0
## 15 Recession 1386Q1 1387Q2 5 0 0 0.0
## 16 Expansion 1387Q2 1387Q4 2 0 0 0.0
## 17 Recession 1387Q4 1388Q2 2 0 0 0.0
## 18 Expansion 1388Q2 1389Q4 6 0 0 0.0
## 19 Recession 1389Q4 <NA> NA 0 NA NA
##
## Amplitude Duration
## Exp=]T;P] 0 5.6
## Rec=]P;T] 0 4.5
plot(dat)
data(MBRI.Iran.Dating)
plot(dat,MBRI.Iran.Dating)
Two stackoverflow threads.
library(historydata)
data.frame(data(package = "historydata")$results)[-c(1, 2)]
## Item
## 1 catholic_dioceses
## 2 early_colleges
## 3 judges_appointments
## 4 judges_people
## 5 naval_promotions
## 6 paulist_missions
## 7 sarna
## 8 tudors
## 9 us_national_population
## 10 us_state_populations
## Title
## 1 Roman Catholic dioceses in the United States, Canada, and Mexico
## 2 Early colleges in the United States
## 3 Federal judges in the United States of America
## 4 Federal judges in the United States of America
## 5 Promotions of U.S. Navy officers, 1798-1849
## 6 Records of missions held by the Paulist Fathers, 1851-1893
## 7 Population estimates of American Jews
## 8 Tudor dynasty
## 9 Population of the United States, 1790-2010
## 10 Populations of US states and territories, 1790-2010
# Roman Catholic dioceses in the United States, Canada, and Mexico
head(catholic_dioceses)
## diocese rite lat long event date
## 1 Baltimore, Maryland Latin 39.29038 -76.61219 erected 1789-04-06
## 2 New Orleans, Louisiana Latin 29.95107 -90.07153 erected 1793-04-25
## 3 Boston, Massachusetts Latin 42.35843 -71.05977 erected 1808-04-08
## 4 Louisville, Kentucky Latin 38.25266 -85.75846 erected 1808-04-08
## 5 New York, New York Latin 40.71435 -74.00597 erected 1808-04-08
## 6 Philadelphia, Pennsylvania Latin 39.95233 -75.16379 erected 1808-04-08
dim(catholic_dioceses)
## [1] 425 6
# Early colleges in the United States
head(early_colleges)
## college original_name city state
## 1 Harvard <NA> Cambridge MA
## 2 William and Mary <NA> Williamsburg VA
## 3 Yale <NA> New Haven CT
## 4 Pennsylvania, Univ. of <NA> Philadelphia PA
## 5 Princeton College of New Jersey Princeton NJ
## 6 Columbia King's College New York NY
## established sponsorship
## 1 1636 Congregational; after 1805 Unitarian
## 2 1693 Anglican
## 3 1701 Congregational
## 4 1740 Nondenominational
## 5 1746 Presbyterian
## 6 1754 Anglican
if(require(ggplot2))
{
ggplot(early_colleges, aes(x = established)) + geom_bar(binwidth = 5) +
ggtitle("Founding Dates of Early American Colleges")
}
## Loading required package: ggplot2
dim(early_colleges)
## [1] 65 6
head(judges_people)
## judge_id name_first name_middle name_last name_suffix birth_date
## 1 3419 Ronnie <NA> Abrams <NA> 1968
## 2 1 Matthew T. Abruzzo <NA> 1889
## 3 2 Marcus Wilson Acheson <NA> 1828
## 4 3 William Marsh Acker Jr. 1927
## 5 4 Harold Arnold Ackerman <NA> 1928
## 6 5 James Waldo Ackerman <NA> 1926
## birthplace_city birthplace_state death_date death_city death_state
## 1 New York NY NA <NA> <NA>
## 2 Brooklyn NY 1971 Potomac MD
## 3 Washington PA 1906 Pittsburgh PA
## 4 Birmingham AL NA <NA> <NA>
## 5 Newark NJ 2009 West Orange NJ
## 6 Jacksonville FL 1984 Springfield IL
## gender race
## 1 F White
## 2 M White
## 3 M White
## 4 M White
## 5 M White
## 6 M White
head(judges_appointments)
## judge_id court_name
## 1 3419 U. S. District Court, Southern District of New York
## 2 1 U. S. District Court, Eastern District of New York
## 3 2 U. S. District Court, Western District of Pennsylvania
## 4 3 U. S. District Court, Northern District of Alabama
## 5 4 U. S. District Court, District of New Jersey
## 6 5 U. S. District Court, Southern District of Illinois
## court_type president_name president_party nomination_date
## 1 USDC Barack Obama Democratic 07/28/2011
## 2 USDC Franklin D. Roosevelt Democratic 02/03/1936
## 3 USDC Rutherford B. Hayes Republican 01/06/1880
## 4 USDC Ronald Reagan Republican 07/22/1982
## 5 USDC Jimmy Carter Democratic 09/28/1979
## 6 USDC Gerald Ford Republican 06/18/1976
## predecessor_last_name predecessor_first_name senate_confirmation_date
## 1 Kaplan Lewis A. 03/22/2012
## 2 new <NA> 02/12/1936
## 3 Ketcham Winthrop 01/14/1880
## 4 McFadden Frank H. 08/18/1982
## 5 Barlow George H. 10/31/1979
## 6 Wood Jr. Harlington 07/02/1976
## commission_date chief_judge_begin chief_judge_end
## 1 03/23/2012 NA NA
## 2 02/15/1936 NA NA
## 3 01/14/1880 NA NA
## 4 08/18/1982 NA NA
## 5 11/02/1979 NA NA
## 6 07/02/1976 NA NA
## retirement_from_active_service termination_date
## 1 <NA> <NA>
## 2 02/15/1966 05/28/1971
## 3 <NA> 02/09/1891
## 4 05/31/1996 <NA>
## 5 02/15/1994 12/02/2009
## 6 <NA> 03/31/1979
## termination_reason
## 1 <NA>
## 2 Death
## 3 Appointment to Another Judicial Position
## 4 <NA>
## 5 Death
## 6 Reassignment
# Population estimates of American Jews
head(sarna)
## year estimate value
## 1 1660 population_low 50
## 2 1700 population_low 200
## 3 1776 population_low 1000
## 4 1790 population_low 1300
## 5 1800 population_low 2500
## 6 1820 population_low 2650
dim(sarna)
## [1] 92 3
head(us_national_population)
## year population
## 1 1790 3929625
## 2 1800 5308483
## 3 1810 7239881
## 4 1820 9638239
## 5 1830 12860702
## 6 1840 17063353
if(require(ggplot2))
{
ggplot(us_national_population,
aes(x = year, y = population)) +
geom_line() +
ggtitle("Population of the United States, 1790-2010")
}
head(us_state_populations)
## GISJOIN year state population
## 1 G090 1790 Connecticut 237655
## 2 G100 1790 Delaware 59096
## 3 G130 1790 Georgia 82548
## 4 G240 1790 Maryland 319728
## 5 G250 1790 Massachusetts 475199
## 6 G330 1790 New Hampshire 141899
dim(us_state_populations)
## [1] 983 4
library(USAboundaries)
data.frame(data(package = "USAboundaries")$results)[-c(1, 2)]
## Item
## 1 hist_us_counties
## 2 hist_us_states
## Title
## 1 Boundaries of counties in the United States of America, 1629-2000
## 2 Boundaries of states in the United States of America, 1783-2000
library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(sp)
library(ggplot2)
us_1844 <- us_boundaries(as.Date("1844-07-04"))
plot(us_1844)
va_1844 <- us_boundaries(as.Date("1844-07-04"), states = "Louisiana", type = "county")
plot(va_1844)
# head(hist_us_counties) # don't run ever
# head(hist_us_states) # don't run ever
library(verification)
## Loading required package: fields
## Loading required package: spam
## Loading required package: grid
## Spam version 1.0-1 (2014-09-09) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
##
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
##
## Loading required package: maps
## Loading required package: boot
## Loading required package: CircStats
## Loading required package: MASS
## Loading required package: dtw
## Loading required package: proxy
##
## Attaching package: 'proxy'
##
## The following objects are masked from 'package:stats':
##
## as.dist, dist
##
## Loaded dtw v1.17-1. See ?dtw for help, citation("dtw") for use in publication.
data.frame(data(package = "verification")$results)[-c(1, 2)]
## Item Title
## 1 disc.dat Discrimination plot dataset.
## 2 pop Probability of precipitation (pop) data.
## 3 precip.ensemble An ensemble of precipitation forecasts
## 4 prob.frcs.dat Probablisitic Forecast Dataset.
## Data from Wilks, table 7.3 page 246.
## Book: Wilks, D. S. (1995) Statistical Methods in the Atmospheric Sciences Chapter 7, San Diego: Aca- demic Press.
y.i <- c(0,0.05, seq(0.1, 1, 0.1))
obar.i <- c(0.006, 0.019, 0.059, 0.15, 0.277, 0.377, 0.511,
0.587, 0.723, 0.779, 0.934, 0.933)
prob.y<- c(0.4112, 0.0671, 0.1833, 0.0986, 0.0616, 0.0366,
0.0303, 0.0275, 0.245, 0.022, 0.017, 0.203)
obar<- 0.162
attribute(y.i, obar.i, prob.y, obar, main = "Sample Attribute Plot")
## Function will work with a prob.bin class objects as well.
## Note this is a random forecast.
obs<- round(runif(100))
pred<- runif(100)
A<- verify(obs, pred, frcst.type = "prob", obs.type = "binary")
## If baseline is not included, baseline values will be calculated from the sample obs.
attribute(A, main = "Alternative plot", xlab = "Alternate x label")
## NULL
## to add a line from another model
obs<- round(runif(100))
pred<- runif(100)
B<- verify(obs, pred, frcst.type = "prob", obs.type = "binary")
## If baseline is not included, baseline values will be calculated from the sample obs.
lines.attrib(B, col = "green")
## Same with confidence intervals
attribute(A, main = "Alternative plot", xlab = "Alternate x label", CI =TRUE)
## NULL
#### add lines to plot
data(pop)
d <- pop.convert()
## internal function used to
## make binary observations for
## the pop figure.
### note the use of bins = FALSE
mod24 <- verify(d$obs_rain, d$p24_rain, bins = FALSE)
## If baseline is not included, baseline values will be calculated from the sample obs.
mod48 <- verify(d$obs_rain, d$p48_rain, bins = FALSE)
## If baseline is not included, baseline values will be calculated from the sample obs.
plot(mod24, freq = FALSE)
## NULL
lines.attrib(mod48, col = "green", lwd = 2, type = "b")
## A simple forecast
head(disc.dat)
## group.id frcst
## 1 0 0.54
## 2 0 0.54
## 3 0 0.42
## 4 0 0.42
## 5 0 0.42
## 6 0 0.42
dim(disc.dat)
## [1] 526 2
discrimination.plot(disc.dat$group.id, disc.dat$frcst, main = "Default Plot")
discrimination.plot(disc.dat$group.id, disc.dat$frcst, main = "New Labels", cex = 1.2,
leg.txt = c("Low", "Med", "High" ) )
discrimination.plot(disc.dat$group.id, disc.dat$frcst, main = "Without Marginal Plots ", marginal = FALSE)
dygraphs is a fast, flexible open source JavaScript charting library.
library(dygraphs)
dygraph(presidents, main = "Presidential Approval") %>%
dyAxis("y", valueRange = c(0, 100)) %>%
dyAnnotation("1950-7-1", text = "A", tooltip = "Korea") %>%
dyAnnotation("1965-1-1", text = "B", tooltip = "Vietnam")
dygraph(presidents, main = "Presidential Approval") %>%
dyAxis("y", valueRange = c(0, 100)) %>%
dyEvent(date = "1950-6-30", "Korea", labelLoc = "bottom") %>%
dyEvent(date = "1965-2-09", "Vietnam", labelLoc = "bottom")
dygraph(nhtemp, main = "New Haven Temperatures") %>%
dyAxis("y", label = "Temp (F)", valueRange = c(40, 60)) %>%
dyOptions(axisLineWidth = 1.5, fillGraph = TRUE, drawGrid = FALSE)
dygraph(nhtemp, main = "New Haven Temperatures") %>%
dySeries("V1", label = "Temperature (F)") %>%
dyLegend(show = "always", hideOnMouseOut = FALSE)
dygraph(nhtemp, main = "New Haven Temperatures") %>%
dyRangeSelector()
dygraph(nhtemp, main = "New Haven Temperatures") %>%
dyRangeSelector(dateWindow = c("1920-01-01", "1960-01-01"))
dygraph(nhtemp, main = "New Haven Temperatures") %>%
dyRangeSelector(height = 20, strokeColor = "")
dygraph(nhtemp, main = "New Haven Temperatures") %>%
dyShading(from = "1920-1-1", to = "1930-1-1") %>%
dyShading(from = "1940-1-1", to = "1950-1-1")
lungDeaths <- cbind(mdeaths, fdeaths)
dygraph(lungDeaths)
lungDeaths <- cbind(ldeaths, mdeaths, fdeaths)
dygraph(lungDeaths, main = "Deaths from Lung Disease (UK)") %>%
dyHighlight(highlightCircleSize = 5, highlightSeriesBackgroundAlpha = 0.2, hideOnMouseOut = FALSE)
lungDeaths <- cbind(ldeaths, mdeaths, fdeaths)
head(lungDeaths)
## ldeaths mdeaths fdeaths
## [1,] 3035 2134 901
## [2,] 2552 1863 689
## [3,] 2704 1877 827
## [4,] 2554 1877 677
## [5,] 2014 1492 522
## [6,] 1655 1249 406
dygraph(lungDeaths, main = "Deaths from Lung Disease (UK)") %>%
dySeries("mdeaths", drawPoints = TRUE, color = "blue") %>%
dySeries("fdeaths", stepPlot = TRUE, color = "red")
dygraph(discoveries, main = "Important Discoveries") %>%
dyRoller(rollPeriod = 5)