| SOCR ≫ | DSPA ≫ | Topics ≫ |
Data Science and Predictive Analytics (UMich HS650)
Specialized Machine Learning Topics
SOCR/MIDAS (Ivo Dinov)
September 2017
In this chapter, we will discuss some technical details about data formats, streaming, optimization of computation, and distributed deployment of optimized learning algorithms. Chapter 21 provides additional optimization details.
The Internet of Things (IoT) leads to a paradigm shift of scientific inference - from static data interrogated in a batch or distributed environment to an on-demand service-based Cloud computing. Here, we will demonstrate how to work with specialized data, data-streams, and SQL databases, as well as develop and assess on-the-fly data modeling, classification, prediction and forecasting methods. Important examples to keep in mind throughout this chapter include high-frequency data delivered real time in hospital ICU’s (e.g., microsecond Electroencephalography signals, EEGs), dynamically changing stock market data (e.g., Dow Jones Industrial Average Index, DJI), and weather patterns.
We will present (1) format conversion and working with XML, SQL, JSON, CSV, SAS and other data objects, (2) visualization of bioinformatics and network data, (3) protocols for managing, classifying and predicting outcomes from data streams, (4) strategies for optimization, improvement of computational performance, parallel (MPI) and graphics (GPU) computing, and (5) processing of very large datasets.
1 Working with specialized data and databases
Unlike the case-studies we saw in the previous chapters, some real world data may not always be nicely formatted, e.g., as CSV files. We must collect, arrange, wrangle, and harmonize scattered information to generate computable data objects that can be further processed by various techniques. Data wrangling and preprocessing may take over 80% of the time researchers spend interrogating complex multi-source data archives. The following procedures will enhance your skills collecting and handling heterogeneous real world data. Multiple examples of handling long-and-wide data, messy and tidy data, and data cleaning strategies can be found in this JSS Tidy Data article by Hadley Wickham.
1.1 Data format conversion
The R package rio imports and exports various types of file formats, e.g., tab-separated (.tsv), comma-separated (.csv), JSON (.json), Stata (.dta), SPSS (.sav and .por), Microsoft Excel (.xls and .xlsx), Weka (.arff), and SAS (.sas7bdat and .xpt).
rio provides three important functions import(), export() and convert(). They are intuitive, easy to understand, and efficient to execute. Take Stata (.dta) files as an example. First we can download 02_Nof1_Data.dta from our datasets folder.
# install.packages("rio")
library(rio)
# Download the SAS .DTA file first locally
# Local data can be loaded by:
#nof1<-import("02_Nof1_Data.dta")
# the data can also be loaded from the server remotely as well:
nof1<-read.csv("https://umich.instructure.com/files/330385/download?download_frd=1")
str(nof1)## 'data.frame': 900 obs. of 10 variables:
## $ ID : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Tx : int 1 1 0 0 1 1 0 0 1 1 ...
## $ SelfEff : int 33 33 33 33 33 33 33 33 33 33 ...
## $ SelfEff25: int 8 8 8 8 8 8 8 8 8 8 ...
## $ WPSS : num 0.97 -0.17 0.81 -0.41 0.59 -1.16 0.3 -0.34 -0.74 -0.38 ...
## $ SocSuppt : num 5 3.87 4.84 3.62 4.62 2.87 4.33 3.69 3.29 3.66 ...
## $ PMss : num 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 ...
## $ PMss3 : num 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 ...
## $ PhyAct : int 53 73 23 36 21 0 21 0 73 114 ...The data is automatically stored as a data frame. Note that rio sets stingAsFactors=FALSE as default.
rio can help us export files into any other format we chose. To do this we have to use the export() function.
#Sys.getenv("R_ZIPCMD", "zip") # Get the C Zip application
Sys.setenv(R_ZIPCMD="E:/Ivo.dir/Ivo_Tools/ZIP/bin/zip.exe")
Sys.getenv("R_ZIPCMD", "zip")## [1] "E:/Ivo.dir/Ivo_Tools/ZIP/bin/zip.exe"export(nof1, "C:/Users/Dinov/Desktop/02_Nof1.xlsx")This line of code exports the Nof1 data in xlsx format located in the R working directory. Mac users may have a problem exporting *.xslx files using rio because of a lack of a zip tool, but still can output other formats such as “.csv”. An alternative strategy to save an xlsx file is to use package xlsx with default row.name=TRUE.
rio also provides a one-step process to convert-and-save data into alternative formats. The following simple code allows us to convert and save the 02_Nof1_Data.dta file we just downloaded into a CSV file.
# convert("02_Nof1_Data.dta", "02_Nof1_Data.csv")
convert("C:/Users/Dinov/Desktop/02_Nof1.xlsx", "C:/Users/Dinov/Desktop/02_Nof1_Data.csv")You can see a new CSV file pop-up in the working directory. Similar transformations are available for other data formats and types.
1.2 Querying data in SQL databases
The CDC Behavioral Risk Factor Surveillance System (BRFSS) Data, 2013-2015. This file for the combined landline and cell phone data set was exported from SAS V9.3 in the XPT transport format. This file contains 330 variables. This format can be imported into SPSS or STATA. Please note: some of the variable labels get truncated in the process of converting to the XPT format
Be careful - this compressed (ZIP) file is over 315MB in size!
# install.packages("Hmisc")
library(Hmisc)## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, round.POSIXt, trunc.POSIXt, unitsmemory.size(max=T)## [1] 111.31pathToZip <- tempfile()
download.file("http://www.socr.umich.edu/data/DSPA/BRFSS_2013_2014_2015.zip", pathToZip)
# let's just pull two of the 3 years of data (2013 and 2015)
brfss_2013 <- sasxport.get(unzip(pathToZip)[1])## Processing SAS dataset LLCP2013 ..brfss_2015 <- sasxport.get(unzip(pathToZip)[3])## Processing SAS dataset LLCP2015 ..dim(brfss_2013); object.size(brfss_2013)## [1] 491773 336## 685581232 bytes# summary(brfss_2013[1:1000, 1:10]) # subsample the data
# report the summaries for
summary(brfss_2013$has_plan)## Length Class Mode
## 0 NULL NULLbrfss_2013$x.race <- as.factor(brfss_2013$x.race)
summary(brfss_2013$x.race)## 1 2 3 4 5 6 7 8 9 NA's
## 376451 39151 7683 9510 1546 2693 9130 37054 8530 25# clean up
unlink(pathToZip)Let’s try to use logistic regression to find out if self-reported race/ethnicity predicts the binary outcome of having a health care plan.
brfss_2013$has_plan <- brfss_2013$hlthpln1 == 1
system.time(
gml1 <- glm(has_plan ~ as.factor(x.race), data=brfss_2013,
family=binomial)
) # report execution time## user system elapsed
## 2.08 0.26 2.35summary(gml1)##
## Call:
## glm(formula = has_plan ~ as.factor(x.race), family = binomial,
## data = brfss_2013)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1862 0.4385 0.4385 0.4385 0.8047
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.293549 0.005649 406.044 <2e-16 ***
## as.factor(x.race)2 -0.721676 0.014536 -49.647 <2e-16 ***
## as.factor(x.race)3 -0.511776 0.032974 -15.520 <2e-16 ***
## as.factor(x.race)4 -0.329489 0.031726 -10.386 <2e-16 ***
## as.factor(x.race)5 -1.119329 0.060153 -18.608 <2e-16 ***
## as.factor(x.race)6 -0.544458 0.054535 -9.984 <2e-16 ***
## as.factor(x.race)7 -0.510452 0.030346 -16.821 <2e-16 ***
## as.factor(x.race)8 -1.332005 0.012915 -103.138 <2e-16 ***
## as.factor(x.race)9 -0.582204 0.030604 -19.024 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 353371 on 491747 degrees of freedom
## Residual deviance: 342497 on 491739 degrees of freedom
## (25 observations deleted due to missingness)
## AIC: 342515
##
## Number of Fisher Scoring iterations: 5Next, we’ll examine the odds (rather the log odds ration, LOR) of having a health care plan (HCP) by race (R). The LORs are calculated for two array dimensions, separately for each race level (presence of health care plan (HCP) is binary, whereas race (R) has 9 levels, \(R1, R2, ..., R9\)). For example, the odds ratio of having a HCP for \(R1:R2\) is:
\[ OR(R1:R2) = \frac{\frac{P \left( HCP \mid R1 \right)}{1 - P \left( HCP \mid R1 \right)}}{\frac{P \left( HCP \mid R2 \right)}{1 - P \left( HCP \mid R2 \right)}} .\]
# install.packages("vcd")
# load the vcd package to compute the LOR
library("vcd")## Loading required package: gridlor_HCP_by_R <- loddsratio(has_plan ~ as.factor(x.race), data = brfss_2013)
lor_HCP_by_R## log odds ratios for has_plan and as.factor(x.race)
##
## 1:2 2:3 3:4 4:5 5:6 6:7
## -0.72167619 0.20990061 0.18228646 -0.78984000 0.57487142 0.03400611
## 7:8 8:9
## -0.82155382 0.74980101Now, let’s see an example of querying a database containing structured relational collection of data records. A query is a machine instruction (typically represented as text) sent by a user to remote database requesting a specific database operation (e.g., search or summary). One database communication protocol relies on SQL (Structured query language). MySQL is an instance of a database management system that supports SQL communication that many web applications utilize, e.g., YouTube, Flickr, Wikipedia, biological databases like GO, ensembl, etc. Below is an example of an SQL query using the package RMySQL. An alternative way to interface an SQL database is using the package RODBC,
# install.packages("DBI"); install.packages("RMySQL")
# install.packages("RODBC"); library(RODBC)
library(DBI)
library(RMySQL)
ucscGenomeConn <- dbConnect(MySQL(),
user='genome',
dbname='hg38',
host='genome-mysql.cse.ucsc.edu')
result <- dbGetQuery(ucscGenomeConn,"show databases;");
# List the DB tables
allTables <- dbListTables(ucscGenomeConn); length(allTables)
# Get dimensions of a table, read and report the head
dbListFields(ucscGenomeConn, "affyU133Plus2")
affyData <- dbReadTable(ucscGenomeConn, "affyU133Plus2"); head(affyData)
# Select a subset, fetch the data, and report the quantiles
subsetQuery <- dbSendQuery(ucscGenomeConn, "select * from affyU133Plus2 where misMatches between 1 and 3")
affySmall <- fetch(subsetQuery); quantile(affySmall$misMatches)
# Get repeat mask
bedFile <- 'repUCSC.bed'
df <- dbSendQuery(ucscGenomeConn,'select genoName,genoStart,genoEnd,repName,swScore, strand,
repClass, repFamily from rmsk') %>%
dbFetch(n=-1) %>%
mutate(genoName = str_replace(genoName,'chr','')) %>%
tbl_df %>%
write_tsv(bedFile,col_names=F)
message('written ', bedFile)
# Once done, close the connection
dbDisconnect(ucscGenomeConn)To complete the above database SQL commands, it requires access to the remote UCSC SQL Genome server and user-specific credentials. The example below can be done by all users, as it relies only on local services.
# install.packages("RSQLite")
library("RSQLite")
# generate an empty DB and store it in RAM
myConnection <- dbConnect(RSQLite::SQLite(), ":memory:")
myConnection## <SQLiteConnection>
## Path: :memory:
## Extensions: TRUEdbListTables(myConnection)## character(0)# Add tables to the local SQL DB
data(USArrests); dbWriteTable(myConnection, "USArrests", USArrests)
dbWriteTable(myConnection, "brfss_2013", brfss_2013)
dbWriteTable(myConnection, "brfss_2015", brfss_2015)
# Check again the DB content
dbListFields(myConnection, "brfss_2013")## [1] "x.state" "fmonth" "idate" "imonth" "iday"
## [6] "iyear" "dispcode" "seqno" "x.psu" "ctelenum"
## [11] "pvtresd1" "colghous" "stateres" "cellfon3" "ladult"
## [16] "numadult" "nummen" "numwomen" "genhlth" "physhlth"
## [21] "menthlth" "poorhlth" "hlthpln1" "persdoc2" "medcost"
## [26] "checkup1" "sleptim1" "bphigh4" "bpmeds" "bloodcho"
## [31] "cholchk" "toldhi2" "cvdinfr4" "cvdcrhd4" "cvdstrk3"
## [36] "asthma3" "asthnow" "chcscncr" "chcocncr" "chccopd1"
## [41] "havarth3" "addepev2" "chckidny" "diabete3" "veteran3"
## [46] "marital" "children" "educa" "employ1" "income2"
## [51] "weight2" "height3" "numhhol2" "numphon2" "cpdemo1"
## [56] "cpdemo4" "internet" "renthom1" "sex" "pregnant"
## [61] "qlactlm2" "useequip" "blind" "decide" "diffwalk"
## [66] "diffdres" "diffalon" "smoke100" "smokday2" "stopsmk2"
## [71] "lastsmk2" "usenow3" "alcday5" "avedrnk2" "drnk3ge5"
## [76] "maxdrnks" "fruitju1" "fruit1" "fvbeans" "fvgreen"
## [81] "fvorang" "vegetab1" "exerany2" "exract11" "exeroft1"
## [86] "exerhmm1" "exract21" "exeroft2" "exerhmm2" "strength"
## [91] "lmtjoin3" "arthdis2" "arthsocl" "joinpain" "seatbelt"
## [96] "flushot6" "flshtmy2" "tetanus" "pneuvac3" "hivtst6"
## [101] "hivtstd3" "whrtst10" "pdiabtst" "prediab1" "diabage2"
## [106] "insulin" "bldsugar" "feetchk2" "doctdiab" "chkhemo3"
## [111] "feetchk" "eyeexam" "diabeye" "diabedu" "painact2"
## [116] "qlmentl2" "qlstres2" "qlhlth2" "medicare" "hlthcvrg"
## [121] "delaymed" "dlyother" "nocov121" "lstcovrg" "drvisits"
## [126] "medscost" "carercvd" "medbills" "ssbsugar" "ssbfrut2"
## [131] "wtchsalt" "longwtch" "dradvise" "asthmage" "asattack"
## [136] "aservist" "asdrvist" "asrchkup" "asactlim" "asymptom"
## [141] "asnoslep" "asthmed3" "asinhalr" "harehab1" "strehab1"
## [146] "cvdasprn" "aspunsaf" "rlivpain" "rduchart" "rducstrk"
## [151] "arttoday" "arthwgt" "arthexer" "arthedu" "imfvplac"
## [156] "hpvadvc2" "hpvadsht" "hadmam" "howlong" "profexam"
## [161] "lengexam" "hadpap2" "lastpap2" "hadhyst2" "bldstool"
## [166] "lstblds3" "hadsigm3" "hadsgco1" "lastsig3" "pcpsaad2"
## [171] "pcpsadi1" "pcpsare1" "psatest1" "psatime" "pcpsars1"
## [176] "pcpsade1" "pcdmdecn" "rrclass2" "rrcognt2" "rratwrk2"
## [181] "rrhcare3" "rrphysm2" "rremtsm2" "misnervs" "mishopls"
## [186] "misrstls" "misdeprd" "miseffrt" "miswtles" "misnowrk"
## [191] "mistmnt" "mistrhlp" "misphlpf" "scntmony" "scntmeal"
## [196] "scntpaid" "scntwrk1" "scntlpad" "scntlwk1" "scntvot1"
## [201] "rcsgendr" "rcsrltn2" "casthdx2" "casthno2" "emtsuprt"
## [206] "lsatisfy" "ctelnum1" "cellfon2" "cadult" "pvtresd2"
## [211] "cclghous" "cstate" "landline" "pctcell" "qstver"
## [216] "qstlang" "mscode" "x.ststr" "x.strwt" "x.rawrake"
## [221] "x.wt2rake" "x.imprace" "x.impnph" "x.chispnc" "x.crace1"
## [226] "x.impcage" "x.impcrac" "x.impcsex" "x.cllcpwt" "x.dualuse"
## [231] "x.dualcor" "x.llcpwt2" "x.llcpwt" "x.rfhlth" "x.hcvu651"
## [236] "x.rfhype5" "x.cholchk" "x.rfchol" "x.ltasth1" "x.casthm1"
## [241] "x.asthms1" "x.drdxar1" "x.prace1" "x.mrace1" "x.hispanc"
## [246] "x.race" "x.raceg21" "x.racegr3" "x.race.g1" "x.ageg5yr"
## [251] "x.age65yr" "x.age.g" "htin4" "htm4" "wtkg3"
## [256] "x.bmi5" "x.bmi5cat" "x.rfbmi5" "x.chldcnt" "x.educag"
## [261] "x.incomg" "x.smoker3" "x.rfsmok3" "drnkany5" "drocdy3."
## [266] "x.rfbing5" "x.drnkdy4" "x.drnkmo4" "x.rfdrhv4" "x.rfdrmn4"
## [271] "x.rfdrwm4" "ftjuda1." "frutda1." "beanday." "grenday."
## [276] "orngday." "vegeda1." "x.misfrtn" "x.misvegn" "x.frtresp"
## [281] "x.vegresp" "x.frutsum" "x.vegesum" "x.frtlt1" "x.veglt1"
## [286] "x.frt16" "x.veg23" "x.fruitex" "x.vegetex" "x.totinda"
## [291] "metvl11." "metvl21." "maxvo2." "fc60." "actin11."
## [296] "actin21." "padur1." "padur2." "pafreq1." "pafreq2."
## [301] "x.minac11" "x.minac21" "strfreq." "pamiss1." "pamin11."
## [306] "pamin21." "pa1min." "pavig11." "pavig21." "pa1vigm."
## [311] "x.pacat1" "x.paindx1" "x.pa150r2" "x.pa300r2" "x.pa30021"
## [316] "x.pastrng" "x.parec1" "x.pastae1" "x.lmtact1" "x.lmtwrk1"
## [321] "x.lmtscl1" "x.rfseat2" "x.rfseat3" "x.flshot6" "x.pneumo2"
## [326] "x.aidtst3" "x.age80" "x.impeduc" "x.impmrtl" "x.imphome"
## [331] "rcsbrac1" "rcsrace1" "rchisla1" "rcsbirth" "typeinds"
## [336] "typework" "has_plan"dbListTables(myConnection); ## [1] "USArrests" "brfss_2013" "brfss_2015"# Retrieve the entire DB table (for the smaller USArrests table)
dbGetQuery(myConnection, "SELECT * FROM USArrests")## Murder Assault UrbanPop Rape
## 1 13.2 236 58 21.2
## 2 10.0 263 48 44.5
## 3 8.1 294 80 31.0
## 4 8.8 190 50 19.5
## 5 9.0 276 91 40.6
## 6 7.9 204 78 38.7
## 7 3.3 110 77 11.1
## 8 5.9 238 72 15.8
## 9 15.4 335 80 31.9
## 10 17.4 211 60 25.8
## 11 5.3 46 83 20.2
## 12 2.6 120 54 14.2
## 13 10.4 249 83 24.0
## 14 7.2 113 65 21.0
## 15 2.2 56 57 11.3
## 16 6.0 115 66 18.0
## 17 9.7 109 52 16.3
## 18 15.4 249 66 22.2
## 19 2.1 83 51 7.8
## 20 11.3 300 67 27.8
## 21 4.4 149 85 16.3
## 22 12.1 255 74 35.1
## 23 2.7 72 66 14.9
## 24 16.1 259 44 17.1
## 25 9.0 178 70 28.2
## 26 6.0 109 53 16.4
## 27 4.3 102 62 16.5
## 28 12.2 252 81 46.0
## 29 2.1 57 56 9.5
## 30 7.4 159 89 18.8
## 31 11.4 285 70 32.1
## 32 11.1 254 86 26.1
## 33 13.0 337 45 16.1
## 34 0.8 45 44 7.3
## 35 7.3 120 75 21.4
## 36 6.6 151 68 20.0
## 37 4.9 159 67 29.3
## 38 6.3 106 72 14.9
## 39 3.4 174 87 8.3
## 40 14.4 279 48 22.5
## 41 3.8 86 45 12.8
## 42 13.2 188 59 26.9
## 43 12.7 201 80 25.5
## 44 3.2 120 80 22.9
## 45 2.2 48 32 11.2
## 46 8.5 156 63 20.7
## 47 4.0 145 73 26.2
## 48 5.7 81 39 9.3
## 49 2.6 53 66 10.8
## 50 6.8 161 60 15.6# Retrieve just the average of one feature
myQuery <- dbGetQuery(myConnection, "SELECT avg(Assault) FROM USArrests"); myQuery## avg(Assault)
## 1 170.76myQuery <- dbGetQuery(myConnection, "SELECT avg(Assault) FROM USArrests GROUP BY UrbanPop"); myQuery## avg(Assault)
## 1 48.00
## 2 81.00
## 3 152.00
## 4 211.50
## 5 271.00
## 6 190.00
## 7 83.00
## 8 109.00
## 9 109.00
## 10 120.00
## 11 57.00
## 12 56.00
## 13 236.00
## 14 188.00
## 15 186.00
## 16 102.00
## 17 156.00
## 18 113.00
## 19 122.25
## 20 229.50
## 21 151.00
## 22 231.50
## 23 172.00
## 24 145.00
## 25 255.00
## 26 120.00
## 27 110.00
## 28 204.00
## 29 237.50
## 30 252.00
## 31 147.50
## 32 149.00
## 33 254.00
## 34 174.00
## 35 159.00
## 36 276.00# Or do it in batches (for the much larger brfss_2013 and brfss_2015 tables)
myQuery <- dbGetQuery(myConnection, "SELECT * FROM brfss_2013")
# extract data in chunks of 2 rows, note: dbGetQuery vs. dbSendQuery
# myQuery <- dbSendQuery(myConnection, "SELECT * FROM brfss_2013")
# fetch2 <- dbFetch(myQuery, n = 2); fetch2
# do we have other cases in the DB remaining?
# extract all remaining data
# fetchRemaining <- dbFetch(myQuery, n = -1);fetchRemaining
# We should have all data in DB now
# dbHasCompleted(myQuery)
# compute the average (poorhlth) grouping by Insurance (hlthpln1)
# Try some alternatives: numadult nummen numwomen genhlth physhlth menthlth poorhlth hlthpln1
myQuery1_13 <- dbGetQuery(myConnection, "SELECT avg(poorhlth) FROM brfss_2013 GROUP BY hlthpln1"); myQuery1_13## avg(poorhlth)
## 1 56.25466
## 2 53.99962
## 3 58.85072
## 4 66.26757# Compare 2013 vs. 2015: Health grouping by Insurance
myQuery1_15 <- dbGetQuery(myConnection, "SELECT avg(poorhlth) FROM brfss_2015 GROUP BY hlthpln1"); myQuery1_15## avg(poorhlth)
## 1 55.75539
## 2 55.49487
## 3 61.35445
## 4 67.62125myQuery1_13 - myQuery1_15## avg(poorhlth)
## 1 0.4992652
## 2 -1.4952515
## 3 -2.5037326
## 4 -1.3536797# reset the DB query
# dbClearResult(myQuery)
# clean up
dbDisconnect(myConnection)1.3 SparQL Queries
The SparQL Protocol and RDF Query Language (SparQL) is a semantic database query language for RDF (Resource Description Framework) data objects. SparQL queries consist of (1) triple patterns, (2) conjunctions, and (3) disjunctions.
The following example uses SparQL to query the prevalence of tuberculosis from the WikiData SparQL server and plot it on a World geographic map.
# install.packages("SPARQL")
# install.packages("rworldmap")
library(SPARQL)## Loading required package: XML## Loading required package: RCurl## Loading required package: bitopslibrary(ggplot2)
library(rworldmap)## Loading required package: sp## ### Welcome to rworldmap ##### For a short introduction type : vignette('rworldmap')# SparQL Formal
# https://www.w3.org/2009/Talks/0615-qbe/
# W3C Turtle - Terse RDF Triple Language:
# https://www.w3.org/TeamSubmission/turtle/#sec-examples
# RDF (Resource Description Framework) is a graphical data model of (subject, predicate, object) triples representing:
# "subject-node to predicate arc to object arc"
# Resources are represented with URIs, which can be abbreviated as prefixed names
# Objects are literals: strings, integers, booleans, etc.
# Syntax
# URIs: <http://example.com/resource> or prefix:name
# Literals:
# "plain string" "13.4""
# xsd:float, or
# "string with language" @en
# Triple: pref:subject other:predicate "object".
wdqs <- "https://query.wikidata.org/bigdata/namespace/wdq/sparql"
query = "PREFIX wd: <http://www.wikidata.org/entity/>
# prefix declarations
PREFIX wdt: <http://www.wikidata.org/prop/direct/>
PREFIX rdfs: <http://www.w3.org/2000/01/rdf-schema#>
PREFIX p: <http://www.wikidata.org/prop/>
PREFIX v: <http://www.wikidata.org/prop/statement/>
PREFIX qualifier: <http://www.wikidata.org/prop/qualifier/>
PREFIX statement: <http://www.wikidata.org/prop/statement/>
# result clause
SELECT DISTINCT ?countryLabel ?ISO3Code ?latlon ?prevalence ?doid ?year
# query pattern against RDF data
# Q36956 Hansen's disease, Leprosy https://www.wikidata.org/wiki/Q36956
# Q15750965 - Alzheimer's disease: https://www.wikidata.org/wiki/Q15750965
# Influenza - Q2840: https://www.wikidata.org/wiki/Q2840
# Q12204 - tuberculosis https://www.wikidata.org/wiki/Q12204
# P699 Alzheimer's Disease ontology ID
# P1193 prevalence: https://www.wikidata.org/wiki/Property:P1193
# P17 country: https://www.wikidata.org/wiki/Property:P17
# Country ISO-3 code: https://www.wikidata.org/wiki/Property:P298
# Location: https://www.wikidata.org/wiki/Property:P625
# Wikidata docs: https://www.mediawiki.org/wiki/Wikidata_query_service/User_Manual
WHERE {
wd:Q12204 wdt:P699 ?doid ; # tuberculosis P699 Disease ontology ID
p:P1193 ?prevalencewithProvenance .
?prevalencewithProvenance qualifier:P17 ?country ;
qualifier:P585 ?year ;
statement:P1193 ?prevalence .
?country wdt:P625 ?latlon ;
rdfs:label ?countryLabel ;
wdt:P298 ?ISO3Code ;
wdt:P297 ?ISOCode .
FILTER (lang(?countryLabel) = \"en\")
# FILTER constraints use boolean conditions to filter out unwanted query results.
# Shortcut: a semicolon (;) can be used to separate two triple patterns that share the same disease (?country is the shared subject above.)
# rdfs:label is a common predicate for giving a human-friendly label to a resource.
}
# query modifiers
ORDER BY DESC(?population)
"
results <- SPARQL(wdqs, query)
resultMatrix <- as.matrix(results$results)
View(resultMatrix)
# join the data to the geo map
sPDF <- joinCountryData2Map(results$results, joinCode = "ISO3", nameJoinColumn = "ISO3Code")## 7 codes from your data successfully matched countries in the map
## 0 codes from your data failed to match with a country code in the map
## 236 codes from the map weren't represented in your data#map the data with no legend
mapParams <- mapCountryData( sPDF
, nameColumnToPlot="prevalence"
# Alternatively , nameColumnToPlot="doid"
, addLegend='FALSE'
)
#add a modified legend using the same initial parameters as mapCountryData
do.call( addMapLegend, c( mapParams
, legendLabels="all"
, legendWidth=0.5
))
text(1, 1, "Partial view of Tuberculosis Prevelance in the World", cex=1)
#do.call( addMapLegendBoxes
# , c(mapParams
# , list(
# legendText=c('Chile', 'US','Brazil','Argentina'),
# x='bottom',title="AD Prevelance",horiz=TRUE)))
# Alternatively: mapCountryData(sPDF, nameColumnToPlot="prevalence", oceanCol="darkblue", missingCountryCol="white")
View(getMap())
# write.csv(file = "C:/Users/Map.csv", getMap())1.4 Real Random Number Generation
We are already familiar with (pseudo) random number generation (e.g., rnorm(100, 10, 4) or runif(100, 10,20)), which generate algorithmically computer values subject to specified distributions. There are also web-services, e.g., random.org, that can provide true random numbers based on atmospheric noise, rather than using a pseudo random number generation protocol. Below is one example of generating a total of 300 numbers arranged in 3 columns, each of 100 rows of random integers (in decimal format) between 100 and 200.
siteURL <- "http://random.org/integers/" # base URL
shortQuery<-"num=300&min=100&max=200&col=3&base=10&format=plain&rnd=new"
completeQuery <- paste(siteURL, shortQuery, sep="?") # concat url and submit query string
rngNumbers <- read.table(file=completeQuery) # and read the data
rngNumbers ## V1 V2 V3
## 1 172 130 170
## 2 163 147 194
## 3 109 129 125
## 4 190 177 166
## 5 122 188 179
## 6 123 142 189
## 7 174 194 150
## 8 154 173 129
## 9 128 195 130
## 10 130 196 174
## 11 165 101 118
## 12 190 136 102
## 13 139 165 151
## 14 117 180 150
## 15 135 168 113
## 16 169 169 124
## 17 152 147 100
## 18 115 183 119
## 19 148 131 159
## 20 117 177 131
## 21 118 186 193
## 22 198 120 173
## 23 147 137 151
## 24 163 102 101
## 25 173 108 194
## 26 124 168 116
## 27 100 109 132
## 28 165 137 104
## 29 161 149 158
## 30 124 176 170
## 31 127 167 144
## 32 104 121 165
## 33 178 168 130
## 34 176 113 152
## 35 120 191 134
## 36 121 189 129
## 37 198 117 186
## 38 191 178 198
## 39 176 163 188
## 40 162 139 119
## 41 163 106 144
## 42 164 132 180
## 43 119 183 118
## 44 119 138 191
## 45 102 123 192
## 46 125 107 156
## 47 123 134 128
## 48 177 147 115
## 49 138 174 162
## 50 176 186 153
## 51 198 123 160
## 52 149 169 106
## 53 117 109 104
## 54 134 106 185
## 55 196 178 135
## 56 164 110 161
## 57 149 129 119
## 58 155 189 180
## 59 196 118 110
## 60 137 131 171
## 61 164 120 138
## 62 118 167 178
## 63 108 197 146
## 64 186 126 199
## 65 133 135 171
## 66 184 180 173
## 67 130 130 152
## 68 103 107 172
## 69 190 143 118
## 70 146 126 133
## 71 197 131 165
## 72 146 130 181
## 73 139 167 148
## 74 142 171 199
## 75 198 165 152
## 76 116 105 107
## 77 188 188 128
## 78 131 114 159
## 79 171 113 159
## 80 125 137 133
## 81 174 107 157
## 82 102 153 108
## 83 138 145 186
## 84 147 133 160
## 85 105 200 188
## 86 127 141 190
## 87 148 111 161
## 88 190 161 153
## 89 185 113 145
## 90 170 120 136
## 91 183 198 126
## 92 191 151 185
## 93 137 179 187
## 94 129 137 138
## 95 191 106 111
## 96 121 199 101
## 97 102 197 138
## 98 181 133 104
## 99 123 160 162
## 100 120 200 1731.5 Downloading the complete text of web pages
RCurl package provides an amazing tool for extracting and scraping information from websites. Let’s install it and extract information from a SOCR website.
# install.packages("RCurl")
library(RCurl)
web<-getURL("http://wiki.socr.umich.edu/index.php/SOCR_Data", followlocation = TRUE)
str(web, nchar.max = 200)## chr "<!DOCTYPE html>\n<html lang=\"en\" dir=\"ltr\" class=\"client-nojs\">\n<head>\n<meta charset=\"UTF-8\" />\n<title>SOCR Data - SOCR</title>\n<meta http-equiv=\"X-UA-Compatible\" conten"| __truncated__The web object looks incomprehensible. This is because most websites are wrapped in XML/HTML hypertext or include JSON formatted meta-data. RCurl deals with special HTML tags and website meta-data.
To deal with the web pages only, httr package would be a better choice than RCurl. It returns a list that makes much more sense.
#install.packages("httr")
library(httr)
web<-GET("http://wiki.socr.umich.edu/index.php/SOCR_Data")
str(web[1:3])## List of 3
## $ url : chr "http://wiki.socr.umich.edu/index.php/SOCR_Data"
## $ status_code: int 200
## $ headers :List of 12
## ..$ date : chr "Sun, 17 Sep 2017 00:10:42 GMT"
## ..$ server : chr "Apache/2.2.15 (Red Hat)"
## ..$ x-powered-by : chr "PHP/5.3.3"
## ..$ x-content-type-options: chr "nosniff"
## ..$ content-language : chr "en"
## ..$ vary : chr "Accept-Encoding,Cookie"
## ..$ expires : chr "Thu, 01 Jan 1970 00:00:00 GMT"
## ..$ cache-control : chr "private, must-revalidate, max-age=0"
## ..$ last-modified : chr "Sat, 22 Oct 2016 21:46:21 GMT"
## ..$ connection : chr "close"
## ..$ transfer-encoding : chr "chunked"
## ..$ content-type : chr "text/html; charset=UTF-8"
## ..- attr(*, "class")= chr [1:2] "insensitive" "list"1.6 Reading and writing XML with the XML package
A combination of the RCurl and the XML packages could help us extract only the plain text in our desired webpages. This would be very helpful to get information from heavy text based websites.
web<-getURL("http://wiki.socr.umich.edu/index.php/SOCR_Data", followlocation = TRUE)
#install.packages("XML")
library(XML)
web.parsed<-htmlParse(web, asText = T)
plain.text<-xpathSApply(web.parsed, "//p", xmlValue)
cat(paste(plain.text, collapse = "\n"))## The links below contain a number of datasets that may be used for demonstration purposes in probability and statistics education. There are two types of data - simulated (computer-generated using random sampling) and observed (research, observationally or experimentally acquired).
##
## The SOCR resources provide a number of mechanisms to simulate data using computer random-number generators. Here are some of the most commonly used SOCR generators of simulated data:
##
## The following collections include a number of real observed datasets from different disciplines, acquired using different techniques and applicable in different situations.
##
## In addition to human interactions with the SOCR Data, we provide several machine interfaces to consume and process these data.
##
## Translate this page:
##
## (default)
##
## Deutsch
##
## Español
##
## Français
##
## Italiano
##
## Português
##
## <U+65E5><U+672C><U+8A9E>
##
## <U+0411><U+044A><U+043B><U+0433><U+0430><U+0440><U+0438><U+044F>
##
## <U+0627><U+0644><U+0627><U+0645><U+0627><U+0631><U+0627><U+062A> <U+0627><U+0644><U+0639><U+0631><U+0628><U+064A><U+0629> <U+0627><U+0644><U+0645><U+062A><U+062D><U+062F><U+0629>
##
## Suomi
##
## <U+0907><U+0938> <U+092D><U+093E><U+0937><U+093E> <U+092E><U+0947><U+0902>
##
## Norge
##
## <U+D55C><U+AD6D><U+C5B4>
##
## <U+4E2D><U+6587>
##
## <U+7E41><U+4F53><U+4E2D><U+6587>
##
## <U+0420><U+0443><U+0441><U+0441><U+043A><U+0438><U+0439>
##
## Nederlands
##
## <U+0395><U+03BB><U+03BB><U+03B7><U+03BD><U+03B9><U+03BA><U+03AC>
##
## Hrvatska
##
## Ceská republika
##
## Danmark
##
## Polska
##
## România
##
## SverigeHere we extracted all plain text between the starting and ending paragraph HTML tags, <p> and </p>.
More information about extracting text from XML/HTML to text via XPath is available here.
1.7 Web-page Data Scraping
The process that extracting data from complete web pages and storing it in structured data format is called scraping. However, before starting a data scrap from a website, we need to understand the underlying HTML structure for that specific website. Also, we have to check the terms of that website to make sure that scraping from this site is allowed.
The R package rvest is a very good place to start “harvesting” data from websites.
To start with, we use read_html() to store SOCR data website into a xmlnode object.
library(rvest)## Loading required package: xml2##
## Attaching package: 'rvest'## The following object is masked from 'package:XML':
##
## xml## The following object is masked from 'package:Hmisc':
##
## htmlSOCR<-read_html("http://wiki.socr.umich.edu/index.php/SOCR_Data")
SOCR## {xml_document}
## <html lang="en" dir="ltr" class="client-nojs">
## [1] <head>\n<meta http-equiv="Content-Type" content="text/html; charset= ...
## [2] <body class="mediawiki ltr sitedir-ltr ns-0 ns-subject page-SOCR_Dat ...From the summary structure of SOCR, we can discover that there are two important hypertext section markups <head> and <body>. Also, notice that the SOCR data website uses <title> and </title> tags to separate title in the <head> section. Let’s use html_node() to extract title information based on this knowledge.
SOCR %>% html_node("head title") %>% html_text()## [1] "SOCR Data - SOCR"Here we used %>% operator, or pipe, to connect two functions. The above line of code creates a chain of functions to operate on the SOCR object. The first function in the chain html_node() extracts the title from head section. Then, html_text() translates HTML formatted hypertext into English. More on R piping can be found in the magrittr package.
Another function, rvest::html_nodes() can be very helpful in scraping. Similar to html_node(), html_nodes() can help us extract multiple nodes in an xmlnode object. Assume that we want to obtain the meta elements (usually page description, keywords, author of the document, last modified, and other metadata) from the SOCR data website. We apply html_nodes() to SOCR object for lines start with <meta in <head> section. It is optional to use html_attrs()(extracts attributes, text and tag name from html) to make texts prettier.
meta<-SOCR %>% html_nodes("head meta") %>% html_attrs()
meta## [[1]]
## http-equiv content
## "Content-Type" "text/html; charset=UTF-8"
##
## [[2]]
## charset
## "UTF-8"
##
## [[3]]
## http-equiv content
## "X-UA-Compatible" "IE=EDGE"
##
## [[4]]
## name content
## "generator" "MediaWiki 1.23.1"
##
## [[5]]
## name content
## "ResourceLoaderDynamicStyles" ""1.8 Parsing JSON from web APIs
Application Programming Interfaces (APIs) allow web-accessible functions to communicate with each other. Today most API is stored in JSON (JavaScript Object Notation) format.
JSON represents a plain text format used for web applications, data structures or objects. Online JSON objects could be retrieved by packages like RCurl and httr. Let’s see a JSON formatted dataset first. We can use 02_Nof1_Data.json in the class file as an example.
library(httr)
nof1<-GET("https://umich.instructure.com/files/1760327/download?download_frd=1")
nof1## Response [https://instructure-uploads.s3.amazonaws.com/account_17700000000000001/attachments/1760327/02_Nof1_Data.json?response-content-disposition=attachment%3B%20filename%3D%2202_Nof1_Data.json%22%3B%20filename%2A%3DUTF-8%27%2702%255FNof1%255FData.json&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Credential=AKIAJFNFXH2V2O7RPCAA%2F20170917%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Date=20170917T001047Z&X-Amz-Expires=86400&X-Amz-SignedHeaders=host&X-Amz-Signature=dc68b4f7c48b732db89cf185356e1228ba2c226b2dd63e04a0f25aa77693d5e0]
## Date: 2017-09-17 00:10
## Status: 200
## Content-Type: application/json
## Size: 109 kB
## [{"ID":1,"Day":1,"Tx":1,"SelfEff":33,"SelfEff25":8,"WPSS":0.97,"SocSuppt...
## {"ID":1,"Day":2,"Tx":1,"SelfEff":33,"SelfEff25":8,"WPSS":-0.17,"SocSuppt...
## {"ID":1,"Day":3,"Tx":0,"SelfEff":33,"SelfEff25":8,"WPSS":0.81,"SocSuppt"...
## {"ID":1,"Day":4,"Tx":0,"SelfEff":33,"SelfEff25":8,"WPSS":-0.41,"SocSuppt...
## {"ID":1,"Day":5,"Tx":1,"SelfEff":33,"SelfEff25":8,"WPSS":0.59,"SocSuppt"...
## {"ID":1,"Day":6,"Tx":1,"SelfEff":33,"SelfEff25":8,"WPSS":-1.16,"SocSuppt...
## {"ID":1,"Day":7,"Tx":0,"SelfEff":33,"SelfEff25":8,"WPSS":0.30,"SocSuppt"...
## {"ID":1,"Day":8,"Tx":0,"SelfEff":33,"SelfEff25":8,"WPSS":-0.34,"SocSuppt...
## {"ID":1,"Day":9,"Tx":1,"SelfEff":33,"SelfEff25":8,"WPSS":-0.74,"SocSuppt...
## {"ID":1,"Day":10,"Tx":1,"SelfEff":33,"SelfEff25":8,"WPSS":-0.38,"SocSupp...
## ...We can see that JSON objects are very simple. The data structure is organized using hierarchies marked by square brackets. Each piece of information is formatted as a {key:value} pair.
The package jsonlite is a very useful tool to import online JSON formatted datasets into data frame directly. Its syntax is very straight forward.
#install.packages("jsonlite")
library(jsonlite)
nof1_lite<-fromJSON("https://umich.instructure.com/files/1760327/download?download_frd=1")
class(nof1_lite)## [1] "data.frame"1.9 Reading and writing Microsoft Excel spreadsheets using XLSX
We can transfer a xlsx dataset into CSV and use read.csv() to load this kind of dataset. However, R provides an alternative read.xlsx() function in package xlsx to simplify this process. Take our 02_Nof1_Data.xls data in the class file as an example. We need to download the file first.
# install.packages("xlsx")
library(xlsx)## Loading required package: rJava##
## Attaching package: 'rJava'## The following object is masked from 'package:RCurl':
##
## clone## Loading required package: xlsxjarsnof1<-read.xlsx("C:/Users/Dinov/Desktop/02_Nof1.xlsx", 1)
str(nof1)## 'data.frame': 900 obs. of 10 variables:
## $ ID : num 1 1 1 1 1 1 1 1 1 1 ...
## $ Day : num 1 2 3 4 5 6 7 8 9 10 ...
## $ Tx : num 1 1 0 0 1 1 0 0 1 1 ...
## $ SelfEff : num 33 33 33 33 33 33 33 33 33 33 ...
## $ SelfEff25: num 8 8 8 8 8 8 8 8 8 8 ...
## $ WPSS : num 0.97 -0.17 0.81 -0.41 0.59 -1.16 0.3 -0.34 -0.74 -0.38 ...
## $ SocSuppt : num 5 3.87 4.84 3.62 4.62 2.87 4.33 3.69 3.29 3.66 ...
## $ PMss : num 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 ...
## $ PMss3 : num 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 ...
## $ PhyAct : num 53 73 23 36 21 0 21 0 73 114 ...The last argument, 1, stands for the first excel sheet, as any excel file may include a large number of tables in it. Also, we can download the xls or xlsx file into our R working directory so that it is easier to find file path.
2 Working with domain-specific data
Powerful Machine-Learning methods have already been applied in many fields. Some of them are very specialized and require unique approaches to address their characteristics.
2.1 Working with bioinformatics data
Genetic data are stored in widely varying formats and usually have more feature variables than observations. They could have 1,000 columns and only 200 rows. One of the commonly used pre-processng steps for such datasets is variable selection. We will talk about this in Chapter 16.
The Bioconductor project created powerful R functionality (packages and tools) for analyzing genomic data, see Bioconductor for more detailed information.
2.2 Visualizing network data
Social network data and graph datasets describe the relations between nodes (vertices) using connections (links or edges) joining the node objects. Assume we have N objects, we can have \(N*(N-1)\) directed links establishing paired associations between the nodes. Let’s use an example with N=4 to demonstrate a simple graph potentially modeling the following linkage table.
| objects | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| 1 | ….. | \(1\rightarrow 2\) | \(1\rightarrow 3\) | \(1\rightarrow 4\) |
| 2 | \(2\rightarrow 1\) | ….. | \(2\rightarrow 3\) | \(2\rightarrow 4\) |
| 3 | \(3\rightarrow 1\) | \(3\rightarrow 2\) | ….. | \(3\rightarrow 4\) |
| 4 | \(4\rightarrow 1\) | \(4\rightarrow 2\) | \(4\rightarrow 3\) | ….. |
If we change the \(a\rightarrow b\) to an indicator variable (0 or 1) capturing whether we have an edge connecting a pair of nodes, then we get the graph adjacency matrix.
Edge lists provide an alternative way to represent network connections. Every line in the list contains a connection between two nodes (objects).
| Vertex | Vertex |
|---|---|
| 1 | 2 |
| 1 | 3 |
| 2 | 3 |
The above edge list is listing three network connections: object 1 is linked to object 2; object 1 is linked to object 3; and object 2 is linked to object 3. Note that edge lists can represent both directed as well as undirected networks or graphs.
We can imagine that if N is very large, e.g., social networks, the data representation and analysis may be resource intense (memory or computation). In R, we have multiple packages that can deal with social network data. One user-friendly example is provided using the igraph package. First, let’s build a toy example and visualize it using this package.
#install.packages("igraph")
library(igraph)##
## Attaching package: 'igraph'## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## uniong<-graph(c(1, 2, 1, 3, 2, 3, 3, 4), n=10)
plot(g)
Here c(1, 2, 1, 3, 2, 3, 3, 4) is an edge list with 4 rows and n=10 means we have 10 nodes (objects) in total. The small arrows in the graph shows us directed network connections. We might notice that 5-10 nodes are scattered out in the graph. This is because they are not included in the edge list, so there are no network connections between them and the rest of the network.
Now let’s examine the co-appearance network of Facebook circles. The data contains anonymized circles (friends lists) from Facebook collected from survey participants using a Facebook app. The dataset only includes edges (circles, 88, 234) connecting pairs of nodes (users, 4, 039) in the ego networks.
The values on the connections represent the number of links/edges within a circle. We have a huge edge-list made of scrambled Facebook user IDs. Let’s load this dataset into R first. The data is stored in a text file. Unlike CSV files, text files in table format need to be imported using read.table(). We are using header=F option to let R know that we don’t have a header in the text file that contains only tab-separated node pairs (indicating the social connections, edges, between Facebook users).
soc.net.data<-read.table("https://umich.instructure.com/files/2854431/download?download_frd=1", sep=" ", header=F)
head(soc.net.data)## V1 V2
## 1 0 1
## 2 0 2
## 3 0 3
## 4 0 4
## 5 0 5
## 6 0 6Now the data is stored in a data frame. To make this dataset ready for igraph processing and visualization, we need to convert soc.net.data into a matrix object.
soc.net.data.mat <- as.matrix(soc.net.data, ncol=2)By using ncol=2, we made a matrix with two columns. The data is now ready and we can apply graph.edgelist().
# remove the first 347 edges (to wipe out the degenerate "0" node)
graph_m<-graph.edgelist(soc.net.data.mat[-c(0:347), ], directed = F)Before we display the social network graph we may want to examine our model first.
summary(graph_m)## IGRAPH aacbcc1 U--- 4038 87887 --This is an extremely brief yet informative summary. The first line U--- 4038 87887 includes potentially four letters and two numbers. The first letter could be U or D indicating undirected or directed edges. A second letter N would mean that the objects set has a “name” attribute. A third letter is for weighted (W) graph. Since we didn’t add weight in our analysis the third letter is empty (“-”). A forth character is an indicator for bipartite graphs (whose vertices can be divided into two disjoint sets and (that is, represent independent sets where each vertex from one set connects to one vertex in the other set). The two numbers following the 4 letters represent the number of nodes and the number of edges, respectively. Now let’s render the graph.
plot(graph_m)
This graph is very complicated. We can still see that some words are surrounded by more nodes than others. To obtain such information we can use degree() function which list the number of edges for each node.
degree(graph_m)Skimming the table we can find that the 107-th user has as many as 1,044 connections, which makes the user a highly-connected hub. Likely, this node may have higher social relevance.
Some edges might be more important than other edges because they serve as a bridge to link a cloud of nodes. To compare their importance, we can use the betweenness centrality measurement. Betweenness centrality measures centrality in a network. High centrality for a specific node indicates influence. betweenness() can help us to calculate this measurement.
betweenness(graph_m)Again, the 107-th node has the highest betweenness centrality (\(3.556221e+06\)).
We can try another example using SOCR hierarchical data, which is also available for dynamic exploration as a tree graph. Let’s read its JSON data source using the jsonlite package.
tree.json<-fromJSON("http://socr.ucla.edu/SOCR_HyperTree.json", simplifyDataFrame = FALSE)This generates a list object representing the hierarchical structure of the network. Note that this is quite different from edge list. There is one root node, its sub nodes are called children nodes, and the terminal notes are call leaf nodes. Instead of presenting the relationship between nodes in pairs, this hierarchical structure captures the level for each node. To draw the social network graph, we need to convert it as a Node object. We can utilize as.Node() function in data.tree package to do so.
# install.packages("data.tree")
library(data.tree)
tree.graph<-as.Node(tree.json, mode = "explicit")Here we use mode="explicit" option to allow “children” nodes to have their own “children” nodes. Now, the tree.json object has been separated into four different node structures - "About SOCR", "SOCR Resources", "Get Started", and "SOCR Wiki". Let’s plot the first one using igraph package.
plot(as.igraph(tree.graph$`About SOCR`), edge.arrow.size=5, edge.label.font=0.05)
In this graph, "About SOCR" that located at the center is the root.
3 Data Streaming
The proliferation of Cloud services and the emergence of modern technology in all aspects of human experiences leads to a tsunami of data much of which is steamed real-time. The interrogation of such voluminous data is an increasingly important area of research. Data streams are ordered, often unbounded sequences of data points created continuously by a data generator. All of the data mining, interrogation and forecasting methods we discuss here are also applicable to data streams.
3.1 Definition
Mathematically, a data stream in an ordered sequence of data points \[Y = \{y_1, y_2, y_3, \cdots, y_t, \cdots \},\] where the (time) index, \(t\), reflects the order of the observation/record, which may be single numbers, simple vectors in multidimensional space, or objects, e.g., structured Ann Arbor Weather (JSON) and its corresponding structured form. Some streaming data is streamed because it’s too large to be downloaded shotgun style and some is streamed because it’s continually generated and serviced. This presents the potential problem of dealing with data streams that may be unlimited.
Notes:
- Data sources: Real or synthetic stream data can be used. Random simulation streams may be created by
rstream. Real stream data may be piped from financial data providers, the WHO, World Bank, NCAR and other sources. - Inference Techniques: Many of the data interrogation techniques we have seem can be employed for dynamic stream data, e.g.,
factas, for PCA,rEMMandbirchfor clustering, etc. Clustering and classification methods capable of processing data streams have been developed, e.g., Very Fast Decision Trees (VFDT), time window-based Online Information Network (OLIN), On-demand Classification, and the APRIORI streaming algorithm. - Cloud distributed computing: Hadoop2/HadoopStreaming, SPARK, Storm3/RStorm provide an environments to expand batch/script-based R tools to the Cloud.
3.2 The stream package
The R stream package provides data stream mining algorithms using fpc, clue, cluster, clusterGeneration, MASS, and proxy packages. In addition, the package streamMOA provides an rJava interface to the Java-based data stream clustering algorithms available in the Massive Online Analysis (MOA) framework for stream classification, regression and clustering.
If you need a deeper exposure to data streaming in R, we recommend you go over the stream vignettes.
3.3 Synthetic example - random Gaussian stream
This example shows the creation and loading of a mixture of 5 random 2D Gaussians, centers at (x_coords, y_coords) with paired correlations rho_corr, representing a simulated data stream.
- Generate the stream:
# install.packages("stream")
library("stream") ## Loading required package: 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##
## Attaching package: 'stream'## The following object is masked from 'package:httr':
##
## write_streamx_coords <- c(0.2,0.3, 0.5, 0.8, 0.9)
y_coords <- c(0.8,0.3, 0.7, 0.1, 0.5)
p_weight <- c(0.1, 0.9, 0.5, 0.4, 0.3) # A vector of probabilities that determines the likelihood of generated a data point from a particular cluster
set.seed(12345)
stream_5G <- DSD_Gaussians(k = 5, d = 2, mu=cbind(x_coords, y_coords), p=p_weight)- K-Means clustering
We will now try k-means and density-based data stream clustering algorithm, D-Stream, where micro-clusters are formed by grid cells of size gridsize with density of a grid cell (Cm) is least 1.2 times the average cell density. The model is updated with the next 500 data points from the stream.
dstream <- DSC_DStream(gridsize = .1, Cm = 1.2)
update(dstream, stream_5G, n = 500)First, let’s run the k-means clustering with \(k=5\) clusters and plot the resulting micro and macro clusters
kmc <- DSC_Kmeans(k = 5)
recluster(kmc, dstream)
plot(kmc, stream_5G, type = "both", xlab="X-axis", ylab="Y-axis")
In this clustering plot, micro-clusters are shown as circles and macro-clusters are shown as crosses and their sizes represent the corresponding cluster weight estimates.
Next try the density-based data stream clustering algorithm D-Stream. Prior to updating the model with the next 1,000 data points from the stream, we specify the grid cells as micro-clusters, grid cell size (gridsize=0.1), and a micro-cluster (Cm=1.2) that specifies the density of a grid cell as a multiple of the average cell density.
dstream <- DSC_DStream(gridsize = 0.1, Cm = 1.2)
update(dstream, stream_5G, n=1000)We can re-cluster the data using k-means with 5 clusters and plot the resulting micro and macro clusters.
km_G5 <- DSC_Kmeans(k = 5)
recluster(km_G5, dstream)
plot(km_G5, stream_5G, type = "both")
Note the subtle changes in the clustering results between kmc and km_G5.
3.4 Sources of Data Streams
3.4.1 Static structure streams
- DSD_BarsAndGaussians generates two uniformly filled rectangular and two Gaussians clusters with different density.
- DSD_Gaussians generates randomly placed static clusters with random multivariate Gaussian distributions.
- DSD_mlbenchData provides streaming access to machine learning benchmark data sets found in the
mlbenchpackage. - DSD_mlbenchGenerator interfaces the generators for artificial data sets defined in the mlbench package.
- DSD_Target generates a ball in circle data set.
- DSD_UniformNoise generates uniform noise in a d-dimensional (hyper) cube.
3.4.2 Concept drift streams
- DSD_Benchmark provides a collection of simple benchmark problems including splitting and joining clusters, and changes in density or size, which can be used as a comprehensive benchmark set for algorithm comparison.
- DSD_MG is a generator to specify complex data streams with concept drift. The shape as well as the behavior of each cluster over time can be specified using keyframes.
- DSD_RandomRBFGeneratorEvents generates streams using radial base functions with noise. Clusters move, merge and split.
3.4.3 Real data streams
- DSD_Memory provides a streaming interface to static, matrix-like data (e.g., a data frame, a matrix) in memory which represent a fixed portion of a data stream. Matrix-like objects also include large objects potentially stored on disk like
ff::ffdf. - DSD_ReadCSV reads data line by line in text format from a file or an open connection and makes it available in a streaming fashion. This way data that is larger than the available main memory can be processed.
- DSD_ReadDB provides an interface to an open result set from a SQL query to a relational database.
3.5 Printing, plotting and saving streams
For DSD objects, some basic stream functions include print(), plot() and write_stream(). These can save part of a data stream to disk. DSD_Memory and DSD_ReadCSV objects also include member functions like reset_stream() to reset the position in the stream to its beginning.
to request a new batch of data points from the stream we use get_points(). This chooses a random cluster (based on the probability weights in p\_weight) and a point is drawn from the multivariate Gaussian distribution (\(mean=mu, covariance\ matrix=\Sigma\)) of that cluster. Below, we pull \(n = 10\) new data points from the stream.
new_p <- get_points(stream_5G, n = 10)
new_p## X1 X2
## 1 0.4017803 0.2999017
## 2 0.4606262 0.5797737
## 3 0.4611642 0.6617809
## 4 0.3369141 0.2840991
## 5 0.8928082 0.5687830
## 6 0.8706420 0.4282589
## 7 0.2539396 0.2783683
## 8 0.5594320 0.7019670
## 9 0.5030676 0.7560124
## 10 0.7930719 0.0937701new_p <- get_points(stream_5G, n = 100, class = TRUE)
head(new_p, n = 20)## X1 X2 class
## 1 0.7915730 0.09533001 4
## 2 0.4305147 0.36953997 2
## 3 0.4914093 0.82120395 3
## 4 0.7837102 0.06771246 4
## 5 0.9233074 0.48164544 5
## 6 0.8606862 0.49399269 5
## 7 0.3191884 0.27607324 2
## 8 0.2528981 0.27596700 2
## 9 0.6627604 0.68988585 3
## 10 0.7902887 0.09402659 4
## 11 0.7926677 0.09030248 4
## 12 0.9393515 0.50259344 5
## 13 0.9333770 0.62817482 5
## 14 0.7906710 0.10125432 4
## 15 0.1798662 0.24967850 2
## 16 0.7985790 0.08324688 4
## 17 0.5247573 0.57527380 3
## 18 0.2358468 0.23087585 2
## 19 0.8818853 0.49668824 5
## 20 0.4255094 0.81789418 3plot(stream_5G, n = 700, method = "pc")
Note that if you add noise to your stream, e.g., stream_Noise <- DSD_Gaussians(k = 5, d = 4, noise = .1, p = c(0.1, 0.5, 0.3, 0.9, 0.1)), then the noise points won’t be part of any clusters will have an NA class label.
3.6 Stream animation
Clusters can be animated over time by animate_data(). Use reset_stream() to start the animation at the beginning of the stream and note that this method is not implemented for streams of class DSD_Gaussians, DSD_R, DSD_data.frame, and DSD. We’ll create a new DSD_Benchmark data stream.
set.seed(12345)
stream_Bench <- DSD_Benchmark(1)
stream_Bench## Benchmark 1: Two clusters moving diagonally from left to right, meeting in
## the center (5% noise).
## Class: DSD_MG, DSD_R, DSD_data.frame, DSD
## With 2 clusters in 2 dimensions. Time is 1library("animation")
reset_stream(stream_Bench)
animate_data(stream_Bench, n=10000, horizon=100, xlim = c(0, 1), ylim = c(0, 1))
This benchmark generator creates two 2D clusters moving in 2D. One moves from top-left to bottom-right, the other from bottom-left to top-right. Then they meet at the center of the domain, the 2 clusters overlap and then split again.
Concept drift in the stream can be depicted by requesting (\(10\)) times \(300\) data points from the stream and animating the plot. Fast-forwarding the stream can be accomplished by requesting, but ignoring, (\(2000\)) points in between the (\(10\)) plots.
for(i in 1:10) {
plot(stream_Bench, 300, xlim = c(0, 1), ylim = c(0, 1))
tmp <- get_points(stream_Bench, n = 2000)
}
reset_stream(stream_Bench)
animate_data(stream_Bench, n=8000, horizon = 120, xlim=c(0, 1), ylim=c(0, 1))
# Animations can be saved as HTML or GIF
#saveHTML(ani.replay(), htmlfile = "stream_Bench_Animation.html")
#saveGIF(ani.replay())Streams can also be saved locally by write_stream(stream_Bench, "dataStreamSaved.csv", n = 100, sep=",") and loaded back in R by DSD_ReadCSV().
3.7 Case-Study: SOCR Knee Pain Data
These data represent the \(X\) and \(Y\) spatial knee-pain locations for over \(8,000\) patients, along with labels about the knee \(F\)ront, \(B\)ack, \(L\)eft and \(R\)ight. Let’s try to read the SOCR Knee Pain Datasest as a stream.
library("XML"); library("xml2"); library("rvest")
wiki_url <- read_html("http://wiki.socr.umich.edu/index.php/SOCR_Data_KneePainData_041409")
html_nodes(wiki_url, "#content")## {xml_nodeset (1)}
## [1] <div id="content" class="mw-body-primary" role="main">\n\t<a id="top ...kneeRawData <- html_table(html_nodes(wiki_url, "table")[[2]])
normalize<-function(x){
return((x-min(x))/(max(x)-min(x)))
}
kneeRawData_df <- as.data.frame(cbind(normalize(kneeRawData$x), normalize(kneeRawData$Y), as.factor(kneeRawData$View)))
colnames(kneeRawData_df) <- c("X", "Y", "Label")
# randomize the rows of the DF as RF, RB, LF and LB labels of classes are sequential
set.seed(1234)
kneeRawData_df <- kneeRawData_df[sample(nrow(kneeRawData_df)), ]
summary(kneeRawData_df)## X Y Label
## Min. :0.0000 Min. :0.0000 Min. :1.000
## 1st Qu.:0.1331 1st Qu.:0.4566 1st Qu.:2.000
## Median :0.2995 Median :0.5087 Median :3.000
## Mean :0.3382 Mean :0.5091 Mean :2.801
## 3rd Qu.:0.3645 3rd Qu.:0.5549 3rd Qu.:4.000
## Max. :1.0000 Max. :1.0000 Max. :4.000# View(kneeRawData_df)We can use the DSD::DSD_Memory class to get a stream interface for matrix or data frame objects, like the Knee pain location dataset. The number of true clusters \(k=4\) in this dataset.
# use data.frame to create a stream (3rd column contains the label assignment)
kneeDF <- data.frame(x=kneeRawData_df[,1], y=kneeRawData_df[,2],
class=as.factor(kneeRawData_df[,3]))
head(kneeDF) ## x y class
## 1 0.1188590 0.5057803 4
## 2 0.3248811 0.6040462 2
## 3 0.3153724 0.4971098 2
## 4 0.3248811 0.4161850 2
## 5 0.6941363 0.5289017 1
## 6 0.3217116 0.4595376 2streamKnee <- DSD_Memory(kneeDF[,c("x", "y")], class=kneeDF[,"class"], loop=T)
streamKnee## Memory Stream Interface
## Class: DSD_Memory, DSD_R, DSD_data.frame, DSD
## With NA clusters in 2 dimensions
## Contains 8666 data points - currently at position 1 - loop is TRUE# Each time we get a point from *streamKnee*, the stream pointer moves to the next position (row) in the data.
get_points(streamKnee, n=10)## x y
## 1 0.11885895 0.5057803
## 2 0.32488114 0.6040462
## 3 0.31537242 0.4971098
## 4 0.32488114 0.4161850
## 5 0.69413629 0.5289017
## 6 0.32171157 0.4595376
## 7 0.06497623 0.4913295
## 8 0.12519810 0.4682081
## 9 0.32329635 0.4942197
## 10 0.30744849 0.5086705streamKnee## Memory Stream Interface
## Class: DSD_Memory, DSD_R, DSD_data.frame, DSD
## With NA clusters in 2 dimensions
## Contains 8666 data points - currently at position 11 - loop is TRUE# Stream pointer is in position 11 now
# We can redirect the current position of the stream pointer by:
reset_stream(streamKnee, pos = 200)
get_points(streamKnee, n=10)## x y
## 200 0.9413629 0.5606936
## 201 0.3217116 0.5664740
## 202 0.3122029 0.6416185
## 203 0.1553090 0.6040462
## 204 0.3645008 0.5346821
## 205 0.3122029 0.5000000
## 206 0.3549921 0.5404624
## 207 0.1473851 0.5260116
## 208 0.1870048 0.6329480
## 209 0.1220285 0.4132948streamKnee## Memory Stream Interface
## Class: DSD_Memory, DSD_R, DSD_data.frame, DSD
## With NA clusters in 2 dimensions
## Contains 8666 data points - currently at position 210 - loop is TRUE3.8 Data Stream clustering and classification (DSC)
Let’s demonstrate clustering using DSC_DStream, which assigns points to cells in a grid. First, initialize the clustering, as an empty cluster and then use the update() function to implicitly alter the mutable DSC object.
dsc_streamKnee <- DSC_DStream(gridsize = 0.1, Cm = 0.4, attraction=T)
dsc_streamKnee## DStream
## Class: DSC_DStream, DSC_Micro, DSC_R, DSC
## Number of micro-clusters: 0
## Number of macro-clusters: 0# stream::update
reset_stream(streamKnee, pos = 1)
update(dsc_streamKnee, streamKnee, n = 500)
dsc_streamKnee## DStream
## Class: DSC_DStream, DSC_Micro, DSC_R, DSC
## Number of micro-clusters: 16
## Number of macro-clusters: 11head(get_centers(dsc_streamKnee))## [,1] [,2]
## [1,] 0.05 0.45
## [2,] 0.05 0.55
## [3,] 0.15 0.35
## [4,] 0.15 0.45
## [5,] 0.15 0.55
## [6,] 0.15 0.65plot(dsc_streamKnee, streamKnee, xlim=c(0,1), ylim=c(0,1))
# plot(dsc_streamKnee, streamKnee, grid = TRUE)
# Micro-clusters are plotted in red on top of gray stream data points
# The size of the micro-clusters indicates their weight - it's proportional to the number of data points represented by each micro-cluster.
# Micro-clusters are shown as dense grid cells (density is coded with gray values).The purity metric represent an external evaluation criterion of cluster quality, which is the proportion of the total number of points that were correctly classified: $ 0Purity = _{i=1}^k { _j a|c_i t_j |} 1$, where \(N\)=number of observed data points, \(k\) = number of clusters, \(c_i\) is the \(i\)th cluster, and \(t_j\) is the classification that has the maximum number of points with \(c_i\) class labels. High purity suggests that we correctly label points.
Next, we can use K-means clustering.
kMeans_Knee <- DSC_Kmeans(k = 5) # choose 4-5 clusters as we have 4 knee labels
recluster(kMeans_Knee, dsc_streamKnee)
plot(kMeans_Knee, streamKnee, type = "both")
animate_data(streamKnee, n=1000, horizon=100, xlim = c(0, 1), ylim = c(0, 1))
# purity <- animate_cluster(kMeans_Knee, streamKnee, n=2500, type="both", xlim=c(0,1), ylim=c(-,1), evaluationMeasure="purity", horizon=10)
animate_cluster(kMeans_Knee, streamKnee, horizon = 100, n = 5000, measure = "purity", plot.args = list(xlim = c(0, 1), ylim = c(0, 1)))## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## Warning in `[<-.factor`(`*tmp*`, is.na(actual), value = 0L): invalid factor
## level, NA generated
## points purity
## 1 1 0.9600000
## 2 101 0.9043478
## 3 201 0.9500000
## 4 301 0.9500000
## 5 401 0.9333333
## 6 501 0.8982872
## 7 601 0.9916667
## 8 701 0.8570115
## 9 801 0.9037037
## 10 901 0.9047619
## 11 1001 1.0000000
## 12 1101 0.9210526
## 13 1201 0.8130435
## 14 1301 0.8815789
## 15 1401 0.9000000
## 16 1501 0.8596491
## 17 1601 0.9738095
## 18 1701 0.9400000
## 19 1801 0.9428571
## 20 1901 0.9259259
## 21 2001 0.9066667
## 22 2101 0.9230769
## 23 2201 0.9111111
## 24 2301 0.8093651
## 25 2401 0.9090909
## 26 2501 0.9685484
## 27 2601 0.9000000
## 28 2701 0.9111111
## 29 2801 0.9153846
## 30 2901 0.8411765
## 31 3001 0.9052632
## 32 3101 0.9200000
## 33 3201 0.9176471
## 34 3301 0.8945055
## 35 3401 0.9333333
## 36 3501 0.9166667
## 37 3601 0.9263158
## 38 3701 0.9100000
## 39 3801 1.0000000
## 40 3901 0.9100000
## 41 4001 0.9181818
## 42 4101 0.8700000
## 43 4201 0.8928571
## 44 4301 0.9217391
## 45 4401 0.9708333
## 46 4501 0.8977273
## 47 4601 0.9500000
## 48 4701 0.9500000
## 49 4801 0.9047619
## 50 4901 0.88500003.9 Evaluation of data stream clustering
# Synthetic Gaussian example
# stream <- DSD_Gaussians(k = 3, d = 2, noise = .05)
# dstream <- DSC_DStream(gridsize = .1)
# update(dstream, stream, n = 2000)
# evaluate(dstream, stream, n = 100)
evaluate(dsc_streamKnee, streamKnee, measure = c("crand", "SSQ", "silhouette"), n = 100, type = c("auto", "micro", "macro"), assign = "micro", assignmentMethod = c("auto", "model", "nn"), noise = c("class", "exclude"))## Evaluation results for micro-clusters.
## Points were assigned to micro-clusters.
##
## cRand SSQ silhouette
## 0.3473634 0.3382900 0.1373143clusterEval <- evaluate_cluster(dsc_streamKnee, streamKnee, measure = c("numMicroClusters", "purity"), n = 5000, horizon = 100)
head(clusterEval)## points numMicroClusters purity
## 1 0 16 0.9555556
## 2 100 17 0.9733333
## 3 200 18 0.9671053
## 4 300 21 0.9687500
## 5 400 21 0.9880952
## 6 500 22 0.9750000plot(clusterEval[ , "points"], clusterEval[ , "purity"], type = "l", ylab = "Avg Purity", xlab = "Points")
animate_cluster(dsc_streamKnee, streamKnee, horizon = 100, n = 5000, measure = "purity", plot.args = list(xlim = c(0, 1), ylim = c(0, 1)))
## points purity
## 1 1 0.9714286
## 2 101 0.9833333
## 3 201 0.9722222
## 4 301 0.9699248
## 5 401 0.9722222
## 6 501 0.9642857
## 7 601 0.9652778
## 8 701 0.9736842
## 9 801 0.9761905
## 10 901 0.9727273
## 11 1001 0.9714286
## 12 1101 0.9736842
## 13 1201 0.9833333
## 14 1301 0.9736842
## 15 1401 0.9857143
## 16 1501 0.9684211
## 17 1601 0.9736842
## 18 1701 0.9629630
## 19 1801 0.9833333
## 20 1901 0.9632353
## 21 2001 0.9647059
## 22 2101 0.9888889
## 23 2201 0.9736842
## 24 2301 0.9809524
## 25 2401 0.9855072
## 26 2501 0.9785714
## 27 2601 0.9562500
## 28 2701 0.9809524
## 29 2801 0.9824561
## 30 2901 0.9750000
## 31 3001 0.9750000
## 32 3101 0.9750000
## 33 3201 0.9714286
## 34 3301 1.0000000
## 35 3401 0.9809524
## 36 3501 0.9702381
## 37 3601 1.0000000
## 38 3701 0.9740260
## 39 3801 0.9855072
## 40 3901 0.9761905
## 41 4001 0.9642857
## 42 4101 0.9761905
## 43 4201 0.9722222
## 44 4301 0.9727273
## 45 4401 0.9814815
## 46 4501 0.9609375
## 47 4601 0.9714286
## 48 4701 0.9722222
## 49 4801 0.9772727
## 50 4901 0.9777778The dsc_streamKnee is the evaluated clustering, where \(n\) is the data points are taken from streamKnee. The evaluation measure can be specified as a vector of character strings. Points are assigned to clusters in dsc_streamKnee using get_assignment() and can be used to assess the quality of the classification. By default, points are assigned to micro-clusters, or can be assigned to macro-cluster centers by assign = "macro". ALso, new points can be assigned to clusters by the rule used in the clustering algorithm by assignmentMethod = "model" or using nearest-neighbor assignment (nn).
4 Optimization and improving the computational performance
Here and in previous chapters, e.g., Chapter 14, we notice that R may sometimes be slow and memory-inefficient. These problems may be severe, especially for datasets with millions of records or when using complex functions. There are packages for processing large datasets and memory optimization – bigmemory, biganalytics, bigtabulate, etc.
4.1 Generalizing tabular data structures with dplyr
We have also seen long execution times when running processes that ingest, store or manipulate huge data.frame objects. The dplyr package, created by Hadley Wickham and Romain Francoi, provides a faster route to manage such large datasets in R. It creates an object called tbl, similar to data.frame, which has an in-memory column-like structure. R reads these objects a lot faster than data frames.
To make a tbl object we can either convert an existing data frame to tbl or connect to an external database. Converting from data frame to tbl is quite easy. All we need to do is call the function as.tbl().
#install.packages("dplyr")
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:igraph':
##
## as_data_frame, groups, union## The following objects are masked from 'package:Hmisc':
##
## combine, src, summarize## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnof1_tbl<-as.tbl(nof1)
nof1_tbl## # A tibble: 900 x 10
## ID Day Tx SelfEff SelfEff25 WPSS SocSuppt PMss PMss3 PhyAct
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 1 33 8 0.97 5.00 4.03 1.03 53
## 2 1 2 1 33 8 -0.17 3.87 4.03 1.03 73
## 3 1 3 0 33 8 0.81 4.84 4.03 1.03 23
## 4 1 4 0 33 8 -0.41 3.62 4.03 1.03 36
## 5 1 5 1 33 8 0.59 4.62 4.03 1.03 21
## 6 1 6 1 33 8 -1.16 2.87 4.03 1.03 0
## 7 1 7 0 33 8 0.30 4.33 4.03 1.03 21
## 8 1 8 0 33 8 -0.34 3.69 4.03 1.03 0
## 9 1 9 1 33 8 -0.74 3.29 4.03 1.03 73
## 10 1 10 1 33 8 -0.38 3.66 4.03 1.03 114
## # ... with 890 more rowsThis looks like a normal data frame. If you are using R Studio by viewing the nof1_tbl you can see the same output as nof1.
4.2 Making data frames faster with data.table
Similar to tbl, the data.table package provides another alternative to data frame object representation. data.table objects are processed in R much faster compared to standard data frames. Also, all of the functions that can accept data frame could be applied to data.table objects as well. The function fread() is able to read a local CSV file directly into a data.table.
#install.packages("data.table")
library(data.table)##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, lastnof1<-fread("C:/Users/Dinov/Desktop/02_Nof1_Data.csv")Another amazing property of data.table is that we can use subscripts to access a specific location in the dataset just like dataset[row, column]. It also allows the selection of rows with Boolean expression and direct application of functions to those selected rows. Note that column names can be used to call the specific column in data.table, whereas with data frames, we have to use the dataset$columnName syntax).
nof1[ID==1, mean(PhyAct)]## [1] 52.66667This useful functionality can also help us run complex operations with only a few lines of code. One of the drawbacks of using data.table objects is that they are still limited by the available system memory.
4.3 Creating disk-based data frames with ff
The ff (fast-files) package allows us to overcome the RAM limitations of finite system memory. For example, it helps with operating datasets with billion rows. ff creates objects in ffdf formats, which is like a map that points to a location of the data on a disk. However, this makes ffdf objects inapplicable for most R functions. The only way to address this problem is to break the huge dataset into small chunks. After processing a batch of these small chunks, we have to combine the results to reconstruct the complete output. This strategy is relevant in parallel computing, which will be discussed in detail in the next section. First, let’s download one of the large datasets in our datasets archive, UQ_VitalSignsData_Case04.csv.
# install.packages("ff")
library(ff)## Loading required package: bit## Attaching package bit## package:bit (c) 2008-2012 Jens Oehlschlaegel (GPL-2)## creators: bit bitwhich## coercion: as.logical as.integer as.bit as.bitwhich which## operator: ! & | xor != ==## querying: print length any all min max range sum summary## bit access: length<- [ [<- [[ [[<-## for more help type ?bit##
## Attaching package: 'bit'## The following object is masked from 'package:data.table':
##
## setattr## The following object is masked from 'package:rJava':
##
## clone## The following object is masked from 'package:RCurl':
##
## clone## The following object is masked from 'package:base':
##
## xor## Attaching package ff## - getOption("fftempdir")=="C:/Users/Dinov/AppData/Local/Temp/Rtmp08lk9C"## - getOption("ffextension")=="ff"## - getOption("ffdrop")==TRUE## - getOption("fffinonexit")==TRUE## - getOption("ffpagesize")==65536## - getOption("ffcaching")=="mmnoflush" -- consider "ffeachflush" if your system stalls on large writes## - getOption("ffbatchbytes")==168034304 -- consider a different value for tuning your system## - getOption("ffmaxbytes")==8401715200 -- consider a different value for tuning your system##
## Attaching package: 'ff'## The following objects are masked from 'package:bit':
##
## clone, clone.default, clone.list## The following object is masked from 'package:rJava':
##
## clone## The following object is masked from 'package:RCurl':
##
## clone## The following objects are masked from 'package:utils':
##
## write.csv, write.csv2## The following objects are masked from 'package:base':
##
## is.factor, is.ordered# vitalsigns<-read.csv.ffdf(file="UQ_VitalSignsData_Case04.csv", header=T)
vitalsigns<-read.csv.ffdf(file="https://umich.instructure.com/files/366335/download?download_frd=1", header=T)As mentioned earlier, we cannot apply functions directly on this object.
mean(vitalsigns$Pulse)## Warning in mean.default(vitalsigns$Pulse): argument is not numeric or
## logical: returning NA## [1] NAFor basic calculations in datasets we can download another package ffbase. This allows operations on ffdf objects using simple tasks like: mathematical operations, query functions, summary statistics and bigger regression models using packages like biglm, which will be mentioned later in this chapter.
# install.packages("ffbase")
library(ffbase)##
## Attaching package: 'ffbase'## The following objects are masked from 'package:ff':
##
## [.ff, [.ffdf, [<-.ff, [<-.ffdf## The following objects are masked from 'package:base':
##
## %in%, tablemean(vitalsigns$Pulse)## [1] 108.71854.4 Using massive matrices with bigmemory
The previously introduced packages include alternatives to data.frames. FOr instance, the bigmemory package creates alternative objects to 2D matrices (second-order tensors). It can store huge datasets and can be divided into small chunks that can be converted to data frames. However, we cannot directly apply machine learning methods on this types of objects. More detailed information about the bigmemory package is available online.
5 Parallel computing
In previous chapters, we saw various machine learning techniques applied as serial computing tasks. The traditional protocol involves: First, applying function 1 to our raw data. Then, using the output from function 1 as an input to function 2. This process in iterated for a series of functions. Finally, we have the terminal output generated by the last function. This serial or linear computing method is straight forward but time consuming and perhaps sub-optimal.
Now we introduce a more efficient way of computing - parallel computing, which provides a mechanism to deal with different tasks at the same time and combine the outputs for all of processes to get the final answer faster. However, parallel algorithms may require special conditions and cannot be applied to all problems. If two tasks have to be run in a specific order, this problem cannot be parallelized.
5.1 Measuring execution time
To measure how much time can be saved for different methods, we can use function system.time().
system.time(mean(vitalsigns$Pulse))## user system elapsed
## 0.01 0.00 0.01This means calculating the mean of Pulse column in the vitalsigns dataset takes 0.001 seconds. These values will vary between computers, operating systems, and states of operations.
5.2 Parallel processing with multiple cores
We will introduce two packages for parallel computing multicore and snow (their core components are included in the package parallel). They both have a different way of multitasking. However, to run these packages, you need to have a relatively modern multicore computer. Let’s check how many cores your computer has. This function parallel::detectCores() provides this functionality. parallel is a base package, so there is no need to install it prior to using it.
library(parallel)
detectCores()## [1] 8So, there are eight (8) cores in my computer. I will be able to run up to 6-8 parallel jobs on this computer.
The multicore package simply uses the multitasking capabilities of the kernel, the computer’s operating system, to “fork” additional R sessions that share the same memory. Imagine that we open several R sessions in parallel and let each of them does part of the work. Now, let’s examine how this can save time when running complex protocols or dealing with large datasets. To start with, we can use the mclapply() function, which is similar to lapply(), which applies functions to a vector and returns a vector of lists. Instead of applying functions to vectors mcapply() divides the complete computational task and delegates portions of it to each available core. We will apply a simple, yet time consuming, task-generating random numbers for demonstrating this procedure. Also, we can use the system.time() to track the time differences.
set.seed(123)
system.time(c1<-rnorm(10000000))## user system elapsed
## 0.60 0.01 0.62# Note the multi core calls may not work on Windows, but will work on Linux/Mac.
#This shows a 2-core and 4-vore invocations
# system.time(c2<-unlist(mclapply(1:2, function(x){rnorm(5000000)}, mc.cores = 2)))
# system.time(c4<-unlist(mclapply(1:4, function(x){rnorm(2500000)}, mc.cores = 4)))
# And here is a Windows (single core invocation)
system.time(c2<-unlist(mclapply(1:2, function(x){rnorm(5000000)}, mc.cores = 1)))## user system elapsed
## 0.64 0.03 0.67The unlist() is used at the end to combine results from different cores into a single vector. Each line of code creates 10,000,000 random numbers. c1 is a regular R command, which used longest time. c2 used two cores to finish the task (each core handle 5,000,000 numbers) and used less time than the first one. c4 used all four cores to finish the task and successfully reduce the time again. We can see that when we use more cores the time is significantly reduced.
The snow package allows parallel computing on multicore multiprocessor machines or a network of multiple machines. It might be more difficult to use but it’s also certainly more flexible. First we can set how many cores we want to use via makeCluster() function.
# install.packages("snow")
library(snow)##
## Attaching package: 'snow'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, clusterSplit, makeCluster,
## parApply, parCapply, parLapply, parRapply, parSapply,
## splitIndices, stopClustercl<-makeCluster(2)This call might cause your computer to pop up a message warning about access though the firewall. To do the same task we can use parLapply() function in the snow package. Note that we have to call the object we created with the previous makeCluster() function.
system.time(c2<-unlist(parLapply(cl, c(5000000, 5000000), function(x) {rnorm(x)})))## user system elapsed
## 0.07 0.13 0.61While using parLapply(), we have to specify the matrix and the function that will be applied to this matrix. Remember to stop the cluster we made after completing the task, to release back the system resources.
stopCluster(cl)5.3 Parallelization using foreach and doParallel
The foreach package provides another option of parallel computing. It relies on a loop-like process basically applying a specified function for each item in the set, which again is somewhat similar to apply(), lapply() and other regular functions. The interesting thing is that these loops can be computed in parallel saving substantial amounts of time. The foreach package alone cannot provide parallel computing. We have to combine it with other packages like doParallel. Let’s reexamine the task of creating a vector of 10,000,000 random numbers. First, register the 4 compute cores using registerDoParallel().
# install.packages("doParallel")
library(doParallel)## Loading required package: foreach## Loading required package: iteratorscl<-makeCluster(4)
registerDoParallel(cl)Then we can examine the time saving foreach command.
#install.packages("foreach")
library(foreach)
system.time(c4<-foreach(i=1:4, .combine = 'c')
%dopar% rnorm(2500000))## user system elapsed
## 0.22 0.11 0.63Here we used four items (each item runs on a separate core), .combine=c allows foreach to combine the results with the parameter c() generating the aggregate result vector.
Also, don’t forget to close the doParallel by registering the sequential backend.
unregister<-registerDoSEQ()5.4 GPU computing
Modern computers have graphics cards, GPU (Graphics Processing Unit), that consists of thousands of cores, however they are very specialized, unlike the standard CPU chip. If we can use this feature for parallel computing, we may reach amazing performance improvements, at the cost of complicating the processing algorithms and increasing the constraints on the data format. Specific disadvantages of GPU computing include relying on a proprietary manufacturer (e.g., NVidia) frameworks and Complete Unified Device Architecture (CUDA) programming language. CUDA allows programming of GPU instructions into a common computing language. This paper provides one example of using GPU computation to improve significantly the performance of advanced neuroimaging and brain mapping processing of multidimensional data.
The R package gputools is created for parallel computing using NVidia CUDA. Detailed GPU computing in R information is available online.
6 Deploying optimized learning algorithms
As we mentioned earlier, some tasks can be parallelized easier than others. In real word situations, we can pick the algorithms that lend themselves well to parallelization. Some of the R packages that allow parallel computing using ML algorithms are listed below.
6.1 Building bigger regression models with biglm
biglm allows training regression models with data from SQL databases or large data chunks obtained from the ff package. The output is similar to teh standard lm() function that builds linear models. However, biglm operates efficiently on massive datasets.
6.2 Growing bigger and faster random forests with bigrf
The bigrf package can be used to train random forests combining the foreach and doParallel packages. In Chapter 14, we presented random forests as machine learners ensembling multiple tree learners. With parallel computing, we can split the task of creating thousands of trees into smaller tasks that can be outsourced to each available compute core. We only need to combine the results at the end. Then, we will obtain the exact same output in a relatively shorter amount of time.
6.3 Training and evaluation models in parallel with caret
Combining the caret package with foreach and we can obtain a powerful method to deal with time-consuming tasks like building a random forest learner. Utilizing the same example we presented in [Chapter 14(http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/14_ImprovingModelPerformance.html), we can see the time difference of utilizing the foreach package.
#library(caret)
system.time(m_rf <- train(CHARLSONSCORE ~ ., data = qol, method = "rf",
metric = "Kappa", trControl = ctrl, tuneGrid = grid_rf))## Loading required package: randomForest## randomForest 4.6-12## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:Hmisc':
##
## combine## The following object is masked from 'package:ggplot2':
##
## margin## user system elapsed
## 128.40 1.38 130.05It took more than a minute to finish this task in standard execution model purely relying on the regular caret function. Below, this same model training completes much faster using parallelization (less than half the time) compared to the standard call above.
set.seed(123)
cl<-makeCluster(4)
registerDoParallel(cl)
getDoParWorkers()## [1] 4system.time(m_rf <- train(CHARLSONSCORE ~ ., data = qol, method = "rf",
metric = "Kappa", trControl = ctrl, tuneGrid = grid_rf))## Warning in .Internal(gc(verbose, reset)): closing unused connection 8 (<-
## nur-cnd4194cfw:11310)## Warning in .Internal(gc(verbose, reset)): closing unused connection 7 (<-
## nur-cnd4194cfw:11310)## Warning in .Internal(gc(verbose, reset)): closing unused connection 6 (<-
## nur-cnd4194cfw:11310)## Warning in .Internal(gc(verbose, reset)): closing unused connection 5 (<-
## nur-cnd4194cfw:11310)## user system elapsed
## 5.14 0.05 49.11unregister<-registerDoSEQ()7 Practice problem
Try to analyze the co-appearance network in the novel “Les Miserablese”. The data contains the weighted network of co-appearances of characters in Victor Hugo’s novel “Les Miserables”. Nodes represent characters as indicated by the labels and edges connect any pair of characters that appear in the same chapter of the book. The values on the edges are the number of such co-appearances.
miserablese<-read.table(“https://umich.instructure.com/files/330389/download?download_frd=1”, sep=“”, header=F) head(miserablese)
Also, try to interrogate some of the larger datasets we have using alternative parallel computing and big data analytics.
