Libraries for Python & R
study the sample() function
totalRows = 15Introduction to R
Introduction to Rfree books for R
Cookbook for R ♦RCookbook bookdown R books bookdown all books bookdown r-programming books R Programming for Data ScienceData Frame
Data Frame is a list of vectors of equal length to create a dataframe: n = c(2,3,5) s = c('a','b','c') b = c(TRUE, FALSE, FALSE) df = data.frame(n,s,b) Components of dataframe: header, column names, data row, name of the row cell single square bracket "[]", comma Functions: nrow(), ncol(), head() Inport Data: read.table("mydata.txt") read.csv("mydata.csv") retrieve the column vector by the double square bracket or the "$" operator mtcars[[9]] mtcars[["am"]] mtcars$am mtcars[,"am"] retrieve a column slice with the single square bracket "[]" mtcars[1] mtcars["mpg"] mtcars[c("mpg", "hp")] Data frame Row Slice mtcars[24,] mtcars[c(3,24),] mtcars["camaro z28",] mtcars[c("datsun 710","camaro z28"),]# MLFundStat and Hangseng Fund Stat
#================= MLFundStat.html the computation is long, it is possible to cut time by adjusting the cutdate variable. this should be modified to new version using r chart.# Start Of R
#================= Sys.setlocale(category = 'LC_ALL', 'Chinese') use the .Rprofile.site file to run R commands for all users when their R session starts. D:\R-3.5.1\etc\Rprofile.site See: Initialization at startup. #c:\R-4.2.1\etc\Rprofile.site #loadhistory("C:\Users\User\Desktop\.Rhistory") check environment: Sys.getenv() This command could be an environment set: Sys.setenv(FAME="/opt/fame") Start Of R Initialization set rstudio locale to check locale: Sys.getlocale() LC_COLLATE=Chinese (Simplified)_China.936;LC_CTYPE=Chinese (Simplified)_China.936;LC_MONETARY=Chinese (Simplified)_China.936;LC_NUMERIC=C;LC_TIME=Chinese (Simplified)_China.936" to permanently change to english: open C:\R-4.0.3 file Rprofile.site add last line: Sys.setlocale("LC_ALL","English")Startup : Initialization at Start of an R Session
Description
In R, the startup mechanism is as follows. Unless --no-environ was given on the command line, R searches for site and user files to process for setting environment variables. The name of the site file is the one pointed to by the environment variableR_ENVIRON; if this is unset, R_HOME/etc/Renviron.site is used (if it exists, which it does not in a "factory-fresh" installation).
The name of the user file can be specified by the R_ENVIRON_USER environment variable; if this is unset, the files searched for are .Renviron in the current or in the user's home directory (in that order).
See "Details" for how the files are read.
Then R searches for the site-wide startup profile file of R code unless the command line option --no-site-file was given.
The path of this file is taken from the value of the R_PROFILE environment variable (after tilde expansion).
If this variable is unset, the default is R_HOME/etc/Rprofile.site, which is used if it exists (which it does not in a "factory-fresh" installation). (it contains settings from the installer in a "factory-fresh" installation). This code is sourced into the base package.
Users need to be careful not to unintentionally overwrite objects in base, and it is normally advisable to use local if code needs to be executed: see the examples.
Then, unless --no-init-file was given, R searches for a user profile, a file of R code.
The path of this file can be specified by the R_PROFILE_USER environment variable (and tilde expansion will be performed).
If this is unset, a file called .Rprofile is searched for in the current directory or in the user's home directory (in that order).
The user profile file is sourced into the workspace.
Note that when the site and user profile files are sourced only the base package is loaded, so objects in other packages need to be referred to by e.g.utils::dump.frames or after explicitly loading the package concerned.
R then loads a saved image of the user workspace from .RData in the current directory if there is one (unless --no-restore-data or --no-restore was specified on the command line).
Next, if a function .First is found on the search path, it is executed as .First().
Finally, function .First.sys() in the base package is run.
This calls require to attach the default packages specified by options("defaultPackages").
If the methods package is included, this will have been attached earlier (by function .OptRequireMethods()) so that namespace initializations such as those from the user workspace will proceed correctly.
A function .First (and .Last) can be defined in appropriate .Rprofile or Rprofile.site files or have been saved in .RData.
If you want a different set of packages than the default ones when you start, insert a call to options in the .Rprofile or Rprofile.site file.
For example, options(defaultPackages = character()) will attach no extra packages on startup (only the base package) (or set R_DEFAULT_PACKAGES=NULL as an environment variable before running R).
Using options(defaultPackages = "") or R_DEFAULT_PACKAGES="" enforces the R system default.
On front-ends which support it, the commands history is read from the file specified by the environment variable R_HISTFILE (default .Rhistory in the current directory) unless --no-restore-history or --no-restore was specified.
The command-line option --vanilla implies --no-site-file, --no-init-file, --no-environ and (except for R CMD) --no-restore Under Windows, it also implies --no-Rconsole, which prevents loading the Rconsole file.
Arguments
Details
Note that there are two sorts of files used in startup: environment files which contain lists of environment variables to be set, and profile files which contain R code. Lines in a site or user environment file should be either comment lines starting with#, or lines of the form name=value. The latter sets the environmental variable name to value, overriding an existing value.
If value contains an expression of the form ${foo-bar}, the value is that of the environmental variable foo if that exists and is set to a non-empty value, otherwise bar.
(If it is of the form ${foo}, the default is "".) This construction can be nested, so bar can be of the same form (as in ${foo-${bar-blah}}).
Note that the braces are essential: for example $HOME will not be interpreted.
Leading and trailing white space in value are stripped. value is then processed in a similar way to a Unix shell: in particular the outermost level of (single or double) quotes is stripped, and backslashes are removed except inside quotes.
On systems with sub-architectures (mainly Windows), the files Renviron.site and Rprofile.site are looked for first in architecture-specific directories, e.g.R_HOME/etc/i386/Renviron.site. And e.g..Renviron.i386 will be used in preference to .Renviron.
See Also
For the definition of the "home" directory on Windows see therw-FAQ Q2.14.
It can be found from a running R by Sys.getenv("R_USER").
.Last for final actions at the close of an R session. commandArgs for accessing the command line arguments.
There are examples of using startup files to set defaults for graphics devices in the help for windows.options. X11 and quartz.
An Introduction to R for more command-line options: those affecting memory management are covered in the help file for Memory.
readRenviron to read .Renviron files.
For profiling code, see Rprof.
Examples
# NOT RUN { ## Example ~/.Renviron on Unix R_LIBS=~/R/library PAGER=/usr/local/bin/less ## Example .Renviron on Windows R_LIBS=C:/R/library MY_TCLTK="c:/Program Files/Tcl/bin" ## Example of setting R_DEFAULT_PACKAGES (from R CMD check) R_DEFAULT_PACKAGES='utils,grDevices,graphics,stats' # this loads the packages in the order given, so they appear on # the search path in reverse order. ## Example of .Rprofile options(width=65, digits=5) options(show.signif.stars=FALSE) setHook(packageEvent("grDevices", "onLoad"), function(...) grDevices::ps.options(horizontal=FALSE)) set.seed(1234) .First = function() cat("\n Welcome to R!\n\n") .Last = function() cat("\n Goodbye!\n\n") ## Example of Rprofile.site local({ # add MASS to the default packages, set a CRAN mirror old = getOption("defaultPackages"); r = getOption("repos") r["CRAN"] = "http://my.local.cran" options(defaultPackages = c(old, "MASS"), repos = r) ## (for Unix terminal users) set the width from COLUMNS if set cols = Sys.getenv("COLUMNS") if(nzchar(cols)) options(width = as.integer(cols)) # interactive sessions get a fortune cookie (needs fortunes package) if (interactive()) fortunes::fortune() }) ## if .Renviron contains FOOBAR="coo\bar"doh\ex"abc\"def'" ## then we get # > cat(Sys.getenv("FOOBAR"), "\n") # coo\bardoh\exabc"def' # }# Encoding Problems
To write text UTF8 encoding on Windows
Firstly, set encoding
options(encoding = "UTF-8")
To write text UTF8 encoding on Windows one has to use the useBytes=T options in functions like writeLines or readLines:
txt = "在"
writeLines(txt, "test.txt", useBytes=T)
readLines("test.txt", encoding="UTF-8")
[1] "在"
writeLines(wholePage, theFilename, useBytes=T)
The UTF-8 BOM is a sequence of bytes at the start of a text stream
(0xEF, 0xBB, 0xBF) that allows the reader to more reliably guess a file as being encoded in UTF-8.
Normally, the BOM is used to signal the endianness of an encoding, but since endianness is irrelevant to UTF-8, the BOM is unnecessary.
According to the Unicode standard, the BOM for UTF-8 files is not recommended
#=================
# Encoding Problems
Sys.getlocale()
getOption("encoding")
options(encoding = "UTF-8")
Encoding(txtstring) = "UTF-8"
Encoding(txtstring)
txtstring
Sys.setlocale
Sys.setlocale(category = 'LC_ALL', 'Chinese')
Sys.setlocale(category = "LC_ALL", locale = "chs")
Sys.setlocale(category = "LC_ALL", locale = "cht") # fanti
Note:
default: options("encoding" = "native.enc")
statTxtFile = "test.txt"
write("建设银行", statTxtFile, append=TRUE)
result file is ansi
add:
options("encoding" = "UTF-8")
write("建设银行", statTxtFile, append=TRUE)
result file is utf-8
mytext = "this is my text"
Encoding(mytext)
options(encoding = "UTF-8")
getOption("encoding")
options(encoding='native.enc')
getOption("encoding")
iconvlist()
theHeader = "http://qt.gtimg.cn/r=2&q=r_hk"
onecode = "02009"
con = url(paste0(theHeader,onecode), encoding = "GB2312")
thepage=readLines(con)
close(con)
Info=unlist(strsplit(thepage,"~"))
codename=Info[2]
codename
Encoding(codename)
==================
readLines(textConnection("Z\u00FCrich", encoding="UTF-8"), encoding="UTF-8")
readLines(filename, encoding="UTF-8")
readLines(con = stdin(), n = -1L, ok = TRUE, warn = TRUE, encoding = "unknown", skipNul = FALSE)
# note! the chiname encoding is ok inside R, but will be wrong when write to file by local pc locale, to solve the problem, set Sys.setlocale(category = 'LC_ALL', 'Chinese')
readLines(con = file("Unicode.txt", encoding = "UCS-2LE"))
close(con)
unique(Encoding(A)) # will most likely be UTF-8
==================
guess_encoding(pageHeader)
pageHeader = repair_encoding(pageHeader, from="utf-8")
pageHeader = repair_encoding(pageHeader, "UTF-8")
iconv(pageHeader, to="UTF-8")
Encoding(pageHeader) = "UTF-8"
Sys.getlocale("LC_ALL")
https://rpubs.com/mauriciocramos/encoding
==================
Read text as UTF-8 encoding
the following reads in encoding twice and works but reasons unknown
readLines(textConnection("Z\u00FCrich", encoding="UTF-8"), encoding="UTF-8")
[1] "Zürich"
==================
the page source claim to be using UTF-8 encoding:
meta http-equiv="Content-Type" content="text/html; charset=utf-8"
So, the question is, are they really using a different enough encoding,
or can we just convert to utf-8, guessing that any errors will be negligible?
A quick and dirty approach just force utf-8 using iconv:
TV_Audio_Video = read_html(iconv(page_source[[1]], to = "UTF-8"), encoding = "utf8")
In general, this is a bad idea - better to specify the encoding it's from.
In this case, maybe the error is theirs, so this quick and dirty approach might be ok.
to remove leading zeros
substr(t,regexpr("[^0]",t),nchar(t))Pop up message in windows 8.1
use the tcl/tk package in R to create a messageBox. Here is a very simple example: require(tcltk) tkmessageBox(title = "Title of message box", message = "Hello, world!", icon = "info", type = "ok") library(tcltk) tk_messageBox(type='ok',message='I am a tkMessageBox!') different types of messagebox (yesno, okcancel, etc). See ?tk_messageBox. or use cmd system('CMD /C "ECHO The R process has finished running && PAUSE"', or use hta in one line: mshta "about:<script>alert('Hello, world!');close()</script>" or mshta "javascript:alert('message');close()" or mshta.exe vbscript:Execute("msgbox ""message"",0,""title"":close") mshta "about:<script src='file://%~f0'></script><script>close()</script>" %* msg = paste0( 'mshta ', "\"about:<script>alert('Hello, world!');close()</script>\"" ) to show web page, use script to create #================= Pop up message in windows 8.1 c.bat: start MessageBox.vbs "This will be shown in a popup." MessageBox.vbs : Set objArgs = WScript.Arguments messageText = objArgs(0) MsgBox messageText in fact, save a file named test.vbs with content: MsgBox "some message" double click the file will run directly # options("scipen"=999) # format(xx, scientific=F) # options("scipen"=100, "digits"=4) # getOption("scipen") # or as.integer(functionResult); df = data.frame(matrix(ncol = 10000, nrow = 0)) colnames(df) = c("a", "b," "c") rm(list=ls()) Extracting a Single, Simple Table The first step is to load the ¡§XML¡¨ package, then use the htmlParse() function to read the html document into an R object, and readHTMLTable() to read the table(s) in the document. The length() function indicates there is a single table in the document, simplifying our work. The plot3d() function in the rgl package library(rgl) open3d() attach(mtcars) plot3d(disp,wt,mpg, col = rainbow(10))library(stringr)
#============ library(stringr) library(htmltools) library(threejs) data(mtcars) data = mtcars[order(mtcars$cyl),] uv = tabulate(mtcars$cyl) col = c(rep("red",uv[4]),rep("yellow",uv[6]),rep("blue",uv[8])) row.names(mtcars) scatterplot3d(data[,c(3,6,1)], labels=row.names(mtcars), size=mtcars$hp/100, flip.y=TRUE, color=col,renderer="canvas") tabulate(bin, nbins = max(1, bin, na.rm = TRUE)) tabulate takes the integer-valued vector bin and counts the number of times each integer occurs in it. tabulate(c(2,3,3,5), nbins = 10) [1] 0 1 2 0 1 0 0 0 0 0 table(c(2,3,3,5)) 2 3 5 1 2 1 tabulate(c(-2,0,2,3,3,5)) # -2 and 0 are ignored [1] 0 1 2 0 1 tabulate(c(-2,0,2,3,3,5), nbins = 3) [1] 0 1 2 tabulate(factor(letters[1:10]) [1] 1 1 1 1 1 1 1 1 1 1Scatterplot3d: 3D graphics - R software and data visualization
1 Install and load scaterplot3d
2 Prepare the data
3 The function scatterplot3d()
4 Basic 3D scatter plots
5 Change the main title and axis labels
6 Change the shape and the color of points
7 Change point shapes by groups
8 Change point colors by groups
9 Change the global appearance of the graph
10 Remove the box around the plot
11 Add grids on scatterplot3d
12 Add bars
13 Modification of scatterplot3d output
14 Add legends
15 Specify the legend position using xyz.convert()
16 Specify the legend position using keywords
17 Customize the legend position
18 Add point labels
19 Add regression plane and supplementary points
There are many packages in R (RGL, car, lattice, scatterplot3d, …) for creating 2 Prepare the data
3 The function scatterplot3d()
4 Basic 3D scatter plots
5 Change the main title and axis labels
6 Change the shape and the color of points
7 Change point shapes by groups
8 Change point colors by groups
9 Change the global appearance of the graph
10 Remove the box around the plot
11 Add grids on scatterplot3d
12 Add bars
13 Modification of scatterplot3d output
14 Add legends
15 Specify the legend position using xyz.convert()
16 Specify the legend position using keywords
17 Customize the legend position
18 Add point labels
19 Add regression plane and supplementary points
Install and load scaterplot3d
install.packages("scatterplot3d") # Install
library("scatterplot3d") # load
Prepare the data
The iris data set will be used:data(iris)
head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
iris data set gives the measurements of the variables sepal length and width, petal length and width, respectively, for 50 flowers from each of 3 species of iris.
The species are Iris setosa, versicolor, and virginica.
The function scatterplot3d()
A simplified format is:scatterplot3d(x, y=NULL, z=NULL)
x, y, z are the coordinates of points to be plotted.
The arguments y and z can be optional depending on the structure of x.
In what cases, y and z are optional variables?
Basic 3D scatter plots
# Basic 3d graphics
scatterplot3d(iris[,1:3])
# Change the angle of point view
scatterplot3d(iris[,1:3], angle = 55)
Change the main title and axis labels
scatterplot3d(iris[,1:3],
main="3D Scatter Plot",
xlab = "Sepal Length (cm)",
ylab = "Sepal Width (cm)",
zlab = "Petal Length (cm)")
Change the shape and the color of points
The argument pch and color can be used:scatterplot3d(iris[,1:3], pch = 16, color="steelblue")
Read more on the different point shapes available in R : Point shapes in R
Change point shapes by groups
shapes = c(16, 17, 18)
shapes = shapes[as.numeric(iris$Species)]
scatterplot3d(iris[,1:3], pch = shapes)
Read more on the different point shapes available in R : Point shapes in R
Change point colors by groups
colors = c("#999999", "#E69F00", "#56B4E9")
colors = colors[as.numeric(iris$Species)]
scatterplot3d(iris[,1:3], pch = 16, color=colors)
Read more about colors in R: colors in R
Change the global appearance of the graph
The arguments below can be used:Remove the box around the plot
scatterplot3d(iris[,1:3], pch = 16, color = colors,
grid=TRUE, box=FALSE)
Note that, the argument Add grids on scatterplot3d
This section describes how to add xy-, xz- and yz- toaddgrids3d(x, y=NULL, z=NULL, grid = TRUE,
col.grid = "grey", lty.grid=par("lty"))
# 1.
Source the function
source('http://www.sthda.com/sthda/RDoc/functions/addgrids3d.r')
# 2.
3D scatter plot
scatterplot3d(iris[, 1:3], pch = 16, grid=FALSE, box=FALSE)
# 3.
Add grids
addgrids3d(iris[, 1:3], grid = c("xy", "xz", "yz"))
The problem on the above plot is that the grids are drawn over the points.
The # 1.
Source the function
source('~/hubiC/Documents/R/function/addgrids3d.r')
# 2.
Empty 3D scatter plot using pch=""
s3d = scatterplot3d(iris[, 1:3], pch = "", grid=FALSE, box=FALSE)
# 3.
Add grids
addgrids3d(iris[, 1:3], grid = c("xy", "xz", "yz"))
# 4.
Add points
s3d$points3d(iris[, 1:3], pch = 16)
The function Add bars
The argumentscatterplot3d(iris[,1:3], pch = 16, type="h",
color=colors)
Modification of scatterplot3d output
Add legends
Specify the legend position using xyz.convert()
The result ofs3d = scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend(s3d$xyz.convert(7.5, 3, 4.5), legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16)
It’s also possible to specify the
Specify the legend position using keywords
# "right" position
s3d = scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("right", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16)
# Use the argument inset
s3d = scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("right", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16, inset = 0.1)
What means the argument # "bottom" position
s3d = scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("bottom", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16)
Using keywords to specify the legend position is very simple.
However, sometimes, there is an overlap between some points and the legend box or between the axis and legend box.
Is there any solution to avoid this overlap?
Yes, there are several solutions using the combination of the following arguments for the function Customize the legend position
# Custom point shapes
s3d = scatterplot3d(iris[,1:3], pch = shapes)
legend("bottom", legend = levels(iris$Species),
pch = c(16, 17, 18),
inset = -0.25, xpd = TRUE, horiz = TRUE)
# Custom colors
s3d = scatterplot3d(iris[,1:3], pch = 16, color=colors)
legend("bottom", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"), pch = 16,
inset = -0.25, xpd = TRUE, horiz = TRUE)
# Custom shapes/colors
s3d = scatterplot3d(iris[,1:3], pch = shapes, color=colors)
legend("bottom", legend = levels(iris$Species),
col = c("#999999", "#E69F00", "#56B4E9"),
pch = c(16, 17, 18),
inset = -0.25, xpd = TRUE, horiz = TRUE)
In the R code above, you can play with the arguments inset, xpd and horiz to see the effects on the appearance of the legend box.
Add point labels
The functionscatterplot3d(iris[,1:3], pch = 16, color=colors)
text(s3d$xyz.convert(iris[, 1:3]), labels = rownames(iris),
cex= 0.7, col = "steelblue")
Add regression plane and supplementary points
The result ofdata(trees)
head(trees)
Girth Height Volume
1 8.3 70 10.3
2 8.6 65 10.3
3 8.8 63 10.2
4 10.5 72 16.4
5 10.7 81 18.8
6 10.8 83 19.7
This data set provides measurements of the girth, height and volume for black cherry trees.
# 3D scatter plot
s3d = scatterplot3d(trees, type = "h", color = "blue",
angle=55, pch = 16)
# Add regression plane
my.lm = lm(trees$Volume ~ trees$Girth + trees$Height)
s3d$plane3d(my.lm)
# Add supplementary points
s3d$points3d(seq(10, 20, 2), seq(85, 60, -5), seq(60, 10, -10),
col = "red", type = "h", pch = 8)
scatterplot3d(data[,c(3,6,1)],
scatterplot3d interactive 3d scatterplots Interactive 3D Scatterplots complete guide to 3D visualization Data Visualization 3D and 4D graph Three.js Fundamentals #============ scatterplot3d(data[,c(3,6,1)], labels=row.names(mtcars), size=mtcars$hp/100, flip.y=TRUE, color=col,renderer="canvas") # Gumball machine N = 100 i = sample(3, N, replace=TRUE) x = matrix(rnorm(N*3),ncol=3) lab = c("small", "bigger", "biggest") scatterplot3d(x, color=rainbow(N), labels=lab[i], size=i, renderer="canvas") # Example 1 from the scatterplot3d package (cf.) z = seq(-10, 10, 0.1) x = cos(z) y = sin(z) scatterplot3d(x,y,z, color=rainbow(length(z)), labels=sprintf("x=%.2f, y=%.2f, z=%.2f", x, y, z)) # Interesting 100,000 point cloud example, should run this with WebGL! N1 = 10000 N2 = 90000 x = c(rnorm(N1, sd=0.5), rnorm(N2, sd=2)) y = c(rnorm(N1, sd=0.5), rnorm(N2, sd=2)) z = c(rnorm(N1, sd=0.5), rpois(N2, lambda=20)-20) col = c(rep("#ffff00",N1),rep("#0000ff",N2)) scatterplot3d(x,y,z, color=col, size=0.25) cat("\014") CLS Screen # match returns a vector of the positions v1 = c("a","b","c","d") v2 = c("g","x","d","e","f","a","c") x = match(v1,v2) 6 NA 7 3 v1 %in% v2 TRUE FALSE TRUE TRUE x = match(v1,v2,nomatch=-1) 6 -1 7 3 %in% returns a logical vector indicating if there is a match or notthis check whether an element is inside a group
#============= this check whether an element is inside a group v = c('a','b','c','e') 'b' %in% vcheck vector includes in 31:37 %in% 0:36
#============= 31:37 %in% 0:36 if(all(31:36 %in% 0:36)){cat("good")} # dmInfo=data.matrix(Info) # convert dataframe to matrix, but the row and column is exchanged # bob = data.frame(lapply(bob, as.character), stringsAsFactors=FALSE) #Change numeric to characters # write.csv(Info,quote=FALSE, row.names = FALSE) # write csv is the proper way to write the datafile # attach an excel file in R: 1: Install packages XLConnect and foreign and run both libraries 2: abcd = readWorksheet(loadWorkbook('file extension'),sheet=1) # allocate vector of size 1.7 Gb Try memory.limit() for the current memory limit Use memory.limit (size=50000) to increase memory limit. Try using a cloud based environment, try using package slam use factors Concatenate and Split Strings in R ================================== use the paste() function to concatenate strsplit() function to split pangram = "The quick brown fox jumps over the lazy dog" strsplit(pangram, " ") "The" "quick" "brown" "fox" "jumps" "over" "the" "lazy" "dog" the unique elements unique() function unique(tolower(words)) "the" "quick" "brown" "fox" "jumps" "over" "lazy" "dog" # find duplicates # the intersect function is used for different set, not in inside a vector # instead, use the duplicated function will be OK. words = unlist(strsplit(pangram, " ")) words = tolower(words) duplicated(words) words[duplicated(words)] arr = sample(1:36,6,replace=TRUE) cat(arr, "\n") arr[duplicated(arr)] # test run remove duplicate items from a vector originalArr = c(1,1,3,4,5,5,6,7,8,8,8,8,9,9) cat(originalArr, "\n") # find out duplicates removeItems = unique(originalArr[duplicated(originalArr)]) # use unique to remove repeated duplicates cat(removeItems, "\n") finalArr = originalArr for(item in removeItems){ cat("remove this:", item," ") cat("they are:", which(finalArr == item)," ") finalArr = finalArr[-(which(finalArr == item))] cat("result vec:", finalArr, "\n") } # unique will not remove duplicates originalArr = unique(originalArr) # rmItems(fmList, itemList) remove itemList from fmList rmItems <- function(fmList, itemList){ commons = unique(fmList[fmList %in% itemList]) for(item in commons){fmList = fmList[-(which(fmList == item))]} return(fmList) } rmItems(originalArr, removeItems) # R base functions duplicated(): for identifying duplicated elements and unique(): for extracting unique elements, distinct() [dplyr package] to remove duplicate rows in a data frame. R split Function ================ split() function divides the data in a vector. unsplit() funtion do the reverse. split(x, f, drop = FALSE, ...) split(x, f, drop = FALSE, ...) = value unsplit(value, f, drop = FALSE) x: vector, data frame f: indices drop: discard non existing levels or not # R not recognizing Chinese characters # I have this saved as a script in RStudio: # this works without problem in windows 8.1 a <- "中文" cat("这是中文", a) aaa = readline(prompt="输入汉字：") cat(aaa) This seems to be a Windows/UTF-8 encoding problem. It works if you use eval(parse('test.R', encoding = 'UTF-8')) instead of source(). I try to use read_csv to read my csv file and the source code as follow: ch4sample <- "D:/Rcode/最近一年內.csv" ch4sample.exp1 <-read_csv(ch4sample, col_names = TRUE) Unfortunately, the R console was showing the error message You might use list.files() function to find out how R names these files, and refer to them that way. For example > list.files() [1] "community_test" "community-sandbox.Rproj" [3] "poobär.r" To source .R file saved using UTF-8 encoding first of all: Sys.setlocale(category = 'LC_ALL', 'Chinese') and then source(filename, encoding = 'UTF-8') but remember to save output file in utflist objects in the working environment
ls() data() will give you a list of the datasets of all loaded packages help(package = "datasets") show structure of datasets dataStr = function(package="datasets", ...) { d = data(package=package, envir=new.env(), ...)$results[,"Item"] d = sapply(strsplit(d, split=" ", fixed=TRUE), "[", 1) d = d[order(tolower(d))] for(x in d){ message(x, ": ", class(get(x))); message(str(get(x)))} } dataStr()x = read.csv("anova.csv",header=T,sep=",")
#============= x = read.csv("anova.csv",header=T,sep=",") Subtype,Gender,Expression A,m,-0.54 A,m,-0.8 Split the "Expression" values into two groups based on "Gender" variable, "f" for female group, and "m" for male group: >g = split(x$Expression, x$Gender) >g $f [1] -0.66 -1.15 -0.30 -0.40 -0.24 -0.92 0.48 -1.68 -0.80 -0.55 -0.11 -1.26 $m [1] -0.54 -0.80 -1.03 -0.41 -1.31 -0.43 1.01 0.14 1.42 -0.16 0.15 -0.62 Calculate the length, mean value of each group: sapply(g,length) f m 135 146 sapply(g,mean) f m -0.3946667 -0.2227397 You may use lapply, return is a list: lapply(g,mean) unsplit() function combines the groups: unsplit(g,x$Gender)Apply
===== m = matrix(data=cbind(rnorm(30, 0), rnorm(30, 2), rnorm(30, 5)), nrow=30, ncol=3) apply(m, 1, mean) a 1 in the second argument, giving the mean of each row. apply(m, 2, mean)giving the mean of each column. apply(m, 2, function(x) length(x[x<0])) # count -ve values apply(m, 2, function(x) is.matrix(x)) apply(m, 2, is.vector) apply(m, 2, function(x) mean(x[x>0])) #========= ma = matrix(c(1:4, 1, 6:8), nrow = 2) apply(ma, 1, table) apply(ma, 1, stats::quantile) apply(ma, 2, mean) apply(m, 2, function(x) length(x[x<0])) sapply lapply rollapply sapply(1:3, function(x) x^2) lapply return a list: lapply(1:3, function(x) x^2) use unlist with lapply to get a vector sapply(1:3, function(x, y) mean(y[,x]), y=m) A=matrix(1:9, 3,3) B=matrix(4:15, 4,3) C=matrix(8:10, 3,2) MyList=list(A,B,C) Z=sapply(MyList,"[", 1,1 ) #========== te=matrix(1:20,nrow=2) sapply(te,mean) # this is a vector, order arrange in matrix direction matrix(sapply(te,mean),nrow=2) # this is changed to matrix subset() apply() sapply() lapply() tapply() aggregate() apply apply a function to the rows or columns of a matrix M = matrix(seq(1,16), 4, 4) apply(M, 1, min) lapply apply a function to each element of a list in turn and get a list back x = list(a = 1, b = 1:3, c = 10:100) lapply(x, FUN = length) sapply apply a function to each element of a list in turn, but you want a vector back x = list(a = 1, b = 1:3, c = 10:100) sapply(x, FUN = length) vapply squeeze some more speed out of sapply x = list(a = 1, b = 1:3, c = 10:100) vapply(x, FUN = length, FUN.VALUE = 0L) mapply apply a function to the 1st elements of each, and then the 2nd elements of each, etc., coercing the result to a vector/array as in sapply Note: mApply(X, INDEX, FUN, …, simplify=TRUE, keepmatrix=FALSE) from Hmisc package is different from mapply(FUN, ..., MoreArgs = NULL, SIMPLIFY = TRUE, USE.NAMES = TRUE) Examples #Sums the 1st elements, the 2nd elements, etc. mapply(sum, 1:5, 1:5, 1:5) [1] 3 6 9 12 15 mapply(rep, 1:4, 4:1) mapply(rep, times = 1:4, x = 4:1) mapply(rep, times = 1:4, MoreArgs = list(x = 42)) mapply(function(x, y) seq_len(x) + y, c(a = 1, b = 2, c = 3), # names from first c(A = 10, B = 0, C = -10)) word = function(C, k) paste(rep.int(C, k), collapse = "") utils::str(mapply(word, LETTERS[1:6], 6:1, SIMPLIFY = FALSE)) mapply(function(x,y){x^y},x=c(2,3),y=c(3,4)) 8 81 values1 = list(a = c(1, 2, 3), b = c(4, 5, 6), c = c(7, 8, 9)) values2 = list(a = c(10, 11, 12), b = c(13, 14, 15), c = c(16, 17, 18)) mapply(function(num1, num2) max(c(num1, num2)), values1, values2) a b c 12 15 18 Map A wrapper to mapply with SIMPLIFY = FALSE, so it is guaranteed to return a list rapply For when you want to apply a function to each element of a nested list structure, recursively tapply For when you want to apply a function to subsets of a vector and the subsets are defined by some other vector, usually a factor lapply is a list apply which acts on a list or vector and returns a list. sapply is a simple lapply (function defaults to returning a vector or matrix when possible) vapply is a verified apply (allows the return object type to be prespecified) rapply is a recursive apply for nested lists, i.e. lists within lists tapply is a tagged apply where the tags identify the subsets apply is generic: applies a function to a matrix's rows or columns by a "wrapper" for tapply. The power of by arises when we want to compute a task that tapply can't handle aggregate can be seen as another a different way of use tapply if we use it in such a way xx = c(1,3,5,7,9,8,6,4,2,1,5) duplicated(xx) xx[duplicated(xx)] Accessing dataframe by names: mtcars["mpg"] QueueNo = 12 mtcars[QueueNo,"mpg"] some functions to remember charToRaw(key) as.raw(key) A motion chart is a dynamic chart to explore several indicators over time. subset(airquality, Temp > 80, select = c(Ozone, Temp)) subset(airquality, Day == 1, select = -Temp) subset(airquality, select = Ozone:Wind) with(airquality, subset(Ozone, Temp > 80)) ## sometimes requiring a logical 'subset' argument is a nuisance nm = rownames(state.x77) start_with_M = nm %in% grep("^M", nm, value = TRUE) subset(state.x77, start_with_M, Illiteracy:Murder) # but in recent versions of R this can simply be subset(state.x77, grepl("^M", nm), Illiteracy:Murder) join 3 dataframes library("plyr") join() function names(gdp)[3] = "GDP" names(life_expectancy)[3] = "LifeExpectancy" names(population)[3] = "Population" gdp_life_exp = join(gdp, life_expectancy) development = join(gdp_life_exp, population) subset() function dev_2005 = subset(development, Year == 2005) dev_2005_big = subset(dev_2005, GDP >= 30000) development_motion = subset(development_complete, Country %in% selection) library(googleVis) gvisMotionChart() function motion_graph = gvisMotionChart(development_motion, idvar = "Country", timevar = "Year") plot(motion_graph) motion_graph = gvisMotionChart(development_motion, idvar = "Country", timevar = "Year", xvar = "GDP", yvar = "LifeExpectancy", sizevar = "Population") development_motion$logGDP = log(development_motion$GDP) motion_graph = gvisMotionChart(development_motion, idvar = "Country", timevar = "Year", xvar = "logGDP", yvar = "LifeExpectancy", sizevar = "Population") my_list[[1]] extracts the first element of the list my_list, and my_list[["name"]] extracts the element in my_list that is called name. If the list is nested you can travel down the heirarchy by recursive subsetting. mylist[[1]][["name"]] is the element called name inside the first element of my_list. A data frame is just a special kind of list, so you can use double bracket subsetting on data frames too. my_df[[1]] will extract the first column of a data frame and my_df[["name"]] will extract the column named name from the data frame. names() and str() is a great way to explore the structure of a list. i in 1:ncol(df) This is a pretty common model for a sequence: a sequence of consecutive integers designed to index over one dimension of our data. What might surprise you is that this isn't the best way to generate such a sequence, especially when you are using for loops inside your own functions. Let's look at an example where df is an empty data frame: df = data.frame() 1:ncol(df) for (i in 1:ncol(df)) { print(median(df[[i]])) } Our sequence is now the somewhat non-sensical: 1, 0. You might think you wouldn't be silly enough to use a for loop with an empty data frame, but once you start writing your own functions, there's no telling what the input will be. A better method is to use the seq_along() function. if you grow the for loop at each iteration (e.g. using c()), your for loop will be very slow. A general way of creating an empty vector of given length is the vector() function. It has two arguments: the type of the vector ("logical", "integer", "double", "character", etc.) and the length of the vector. Then, at each iteration of the loop you must store the output in the corresponding entry of the output vector, i.e. assign the result to output[[i]]. (You might ask why we are using double brackets here when output is a vector. It's primarily for generalizability: this subsetting will work whether output is a vector or a list.) A time series can be thought of as a vector or matrix of numbers, along with some information about what times those numbers were recorded. This information is stored in a ts object in R. read in some time series data from an xlsx file using read_excel(), a function from the readxl package, and store the data as a ts object. Use the read_excel() function to read the data from "exercise1.xlsx" into mydata. mydata = read_excel("exercise1.xlsx") Create a ts object called myts using the ts() function. myts = ts(mydata[,2:4], start = c(1981, 1), frequency = 4) The first step in any data analysis task is to plot the data. Graphs enable you to visualize many features of the data, including patterns, unusual observations, changes over time, and relationships between variables. The features that you see in the plots must then be incorporated into the forecasting methods that you use. Just as the type of data determines which forecasting method to use, it also determines which graphs are appropriate. You will use the autoplot() function to produce time plots of the data. In each plot, look out for outliers, seasonal patterns, and other interesting features. Use which.max() to spot the outlier in the gold series. library("fpp2") autoplot(a10) ggseasonplot(a10) An interesting variant of a season plot uses polar coordinates, where the time axis is circular rather than horizontal. ggseasonplot(a10, polar = TRUE) beer = window(a10, start=1992) autoplot(beer) ggseasonplot(beer) Use the window() function to consider only the ausbeer data from 1992 and save this to beer. Set a keyword start to the appropriate year. x = tryCatch( readLines("wx.qq.com/"), warning=function(w){ return(paste( "Warning:", conditionMessage(w)));}, error = function(e) { return(paste( "this is Error:", conditionMessage(e)));}, finally={print("This is try-catch test. check the output.")});x = c(sort(sample(1:20, 9)), NA)
#=================== x = c(sort(sample(1:20, 9)), NA) y = c(sort(sample(3:23, 7)), NA) union(x, y) intersect(x, y) setdiff(x, y) setdiff(y, x) setequal(x, y) alist = readLines("alist.txt") blist = readLines("blist.txt") out = setdiff(blist, alist) writeClipboard(out) use of sample command: newData = sample[sample$x > 0 & sample$y > 0.4, ]# To skip 3rd iteration and go to next iteration
#=================== # To skip 3rd iteration and go to next iteration for(n in 1:5) { if(n==3) next cat(n) }googleVis chart
#=================== googleVis chart =============== library(googleVis) Line chart ========== df=data.frame(country=c("US", "GB", "BR"), val1=c(10,13,14), val2=c(23,12,32)) Line = gvisLineChart(df) plot(Line) Scatter chart ======================= # example 1 dat = data.frame(x=c(1,2,3,4,5), y1=c(0,3,7,5,2), y2=c(1,NA,0,3,2)) plot(gvisScatterChart(dat, options=list(lineWidth=2, pointSize=2, width=900, height=600))) # example 2, women Scatter = gvisScatterChart(women, options=list( legend="none", lineWidth=1, pointSize=2, title="Women", vAxis="{title:'weight (lbs)'}", hAxis="{title:'height (in)'}", width=900, height=600) ) plot(Scatter) # example 3 ex3dat = data.frame(x=c(1,2,3,4,5,6,7,8), y1=c(0,3,7,5,2,0,8,6), y2=c(1,NA,0,3,2,6,4,2)) ex3 = gvisScatterChart(ex3dat, options=list( legend="none", lineWidth=1, pointSize=2, title="ex3", vAxis="{title:'weight (lbs)'}", hAxis="{title:'height (in)'}", width=900, height=600) ) plot(ex3) # Note: to plot timeline chart, arrange the time in x axis, beginning with -ve and the last is 1 to show the sequencecat to a file using file(filename, open = "a")
cat("TITLE extra line", "2 3 5 7", "11 13 17", file = "data.txt", sep = "\n")cat append to a file, open file in "a" mode
#=================== textVector = c("First thing","Second thing","c") catObj = file("theappend.txt", open = "a") cat(textVector, file = catObj, sep="\n") close(catObj)install.packages("readr")
#=================== install.packages("readr") library(readr) to read rectangular data (like csv, tsv, and fwf) readr is part of the core tidyverse library(tidyverse) readr supports seven file formats with seven read_ functions: read_csv(): comma separated (CSV) files read_tsv(): tab separated files read_delim(): general delimited files read_fwf(): fixed width files read_table(): tabular files where columns are separated by white-space. read_log(): web log filesiconv(keyword, "unknown", "GB2312")
#=================== iconv(keyword, "unknown", "GB2312")Grabbing HTML Tags
#========== Grabbing HTML Tags \b[^>]*>(.*?) matches the opening and closing pair of a specific HTML tag. Anything between the tags is captured into the first backreference. The question mark in the regex makes the star lazy, to make sure it stops before the first closing tag rather than before the last, like a greedy star would do. This regex will not properly match tags nested inside themselves, like in one two one. <([A-Z][A-Z0-9]*)\b[^>]*>(.*?) will match the opening and closing pair of any HTML tag. Be sure to turn off case sensitivity. The key in this solution is the use of the backreference \1 in the regex. Anything between the tags is captured into the second backreference. This solution will also not match tags nested in themselvesfind the new item
#========== find the new item theList = c("00700","02318","02007") newList=c("03333","01398","02007") newList[!(newList %in% theList)]formating numbers
#========== formating numbers a = seq(1,101,25) sprintf("%03d", a) format(round(a, 2), nsmall = 2)the match function:
#========== the match function: match(x, table, nomatch = NA_integer_, incomparables = NULL) %in% match returns a vector of the positions of (first) matches of its first argument in its second. Corpus= c('animalada', 'fe', 'fernandez', 'ladrillo') Lexicon= c('animal', 'animalada', 'fe', 'fernandez', 'ladr', 'ladrillo') Lexicon %in% Corpus Lexicon[Lexicon %in% Corpus]Machine Learning:
Machine Learning with R and TensorFlow machine-learning-in-r-step-by-step An Introduction to Machine Learning with R mxnet image classification Image Recognition & Classification with Keras #========== Machine Learning: The caret package Caret contains wrapper functions that allow you to use the exact same functions for training and predicting with dozens of different algorithms. On top of that, it includes sophisticated built-in methods for evaluating the effectiveness of the predictions you get from the model. Use The Titanic dataset Training a model training a bunch of different decision trees and having them vote Random forests work pretty well in *lots* of different situations, so I often try them first. Evaluating the model Cross-validation is a way to evaluate the performance of a model without needing any other data than the training data. Making predictions on the test set Improving the modelto handle error 404 when scraping: use tryCatch()
#========== to handle error 404 when scraping: use tryCatch() for (i in urls) { tmp = tryCatch(readLines(url(i), warn=F), error = function (e) NULL) if (is.null(tmp)) { next() # skip to the next url. } } #========== try(readLines(url), silent = TRUE) tryCatch(readLines(url), error = function (e) conditionMessage(e))write.table
write.table(matrixname, file = "outputname", append = FALSE, quote = FALSE, sep = "\t", #========== write.table(matrixname, file = "outputname", append = FALSE, quote = FALSE, sep = "\t", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"), fileEncoding = "") write.table(finalTableList, theOutputname, row.names=FALSE, col.names=FALSE, quote = FALSE, sep = "\t" )Add Gaussian noise to vector
# Create a vector x x <- 1:10 # Add Gaussian noise with mean 0 and standard deviation 0.1 noise <- rnorm(length(x), mean = 0, sd = 0.1) # Add noise to the vector x x_noisy <- x + noise # Print the original and noisy vectors print(x) print(x_noisy)Generate Random Numbers
Method 1: Generate One Random Number #generate one random number between 1 and 20 runif(n=1, min=1, max=20) Method 2: Generate Multiple Random Numbers #generate five random numbers between 1 and 20 runif(n=5, min=1, max=20) Method 3: Generate One Random Integer in sample pool sample(1:20, 1) Method 4: Generate Multiple Random Integers in sample pool #generate five random integers between 1 and 20 (sample with replacement) sample(1:20, 5, replace=TRUE) #generate five random integers between 1 and 20 (sample without replacement) sample(1:20, 5, replace=FALSE) # Generate Random Number From Uniform Distribution > runif(1) # generates 1 random number [1] 0.3984754 > runif(3) # generates 3 random number [1] 0.8090284 0.1797232 0.6803607 > runif(3, min=5, max=10) # define the range between 5 and 10 [1] 7.099781 8.355461 5.173133 # Generate Random Number From Normal Distribution > rnorm(1) # generates 1 random number [1] 1.072712 > rnorm(3) # generates 3 random number [1] -1.1383656 0.2016713 -0.4602043 > rnorm(3, mean=10, sd=2) # provide our own mean and standard deviation [1] 9.856933 9.024286 10.822507Four normal distribution functions:
#========== Four normal distribution functions: R - Normal Distribution Distribution of data is normal means, on plotting a graph with the value of the variable in the horizontal axis and the count of the values in the vertical axis we get a bell shape curve. The center of the curve represents the mean of the data set. In the graph, fifty percent of values lie to the left of the mean and the other fifty percent lie to the right of the graph. This is referred as normal distribution in statistics. R has four in built functions to generate normal distribution. They are: dnorm(x, mean, sd) pnorm(x, mean, sd) qnorm(p, mean, sd) rnorm(n, mean, sd)dnorm()
This function gives height of the probability distribution at each point for a given mean and standard deviation. # Create a sequence of numbers between -10 and 10 incrementing by 0.1. x <- seq(-10, 10, by = .1) # Choose the mean as 2.5 and standard deviation as 0.5. y <-
pnorm()
This function gives the probability of a normally distributed random number to be less that the value of a given number. It is also called "
qnorm()
This function takes the probability value and gives a number whose cumulative value matches the probability value. # Create a sequence of probability values incrementing by 0.02. x <- seq(0, 1, by = 0.02) # Choose the mean as 2 and standard deviation as 3. y <-
rnorm()
This function is used to generate random numbers whose distribution is normal. It takes the sample size as input and generates that many random numbers. We draw a histogram to show the distribution of the generated numbers. # Create a sample of 50 numbers which are normally distributed. y <-
RNORM Generates random numbers from normal distribution
rnorm(n, mean, sd)
rnorm(1000, 3, .25) Generates 1000 numbers from a normal with mean 3 and sd=.25
DNORM Probability Density Function(PDF)
dnorm(x, mean, sd)
dnorm(0, 0, .5) Gives the density (height of the PDF) of the normal with mean=0 and sd=.5.
dnorm returns the value of the normal distribution given parameters for x, μ, and σ.
# x = 0, mu = 0 and sigma = 0
dnorm(0, mean = 0, sd = 1)
dnorm(1, mean = 1.2, sd = 0.5) # result: 0.7365403
change x to dataset
dataset = seq(-3, 3, by = .1)
dvalues = dnorm(dataset)
plot(dvalues, # y = values and x = index
xaxt = "n", # Don't label the x-axis
type = "l", # Make it a line plot
main = "pdf of the Standard Normal",
xlab= "Data Set")
compare the data with dnorm:
dataset = c( 5, 1,2,5,3,5,6,4,7,4,5,4,8,6,3,3,6,5,4,3,4,3,4,3)
plot(dvalues, # y = values and x = index
xaxt = "n", # Don't label the x-axis
type = "l", # Make it a line plot
main = "pdf of the Standard Normal",
xlab= "Data Set")
to create a dnorm of a dataset to compare with current dataset
make a cut index
cutindex = seq(min(dataset),max(dataset),length = 10)
yfit = dnorm(cutindex, mean=mean(dataset), sd=sd(dataset))
lines(cutindex, yfit)
# Kernel Density Plot
d = density(mtcars$mpg) # returns the density data
plot(d) # plots the results
# Filled Density Plot
d = density(mtcars$mpg)
plot(d, main="Kernel Density of Miles Per Gallon")
polygon(d, col="red", border="blue")
Kernel density estimation is a technique that let's you create a smooth curve given a set of data.
PNORM Cumulative Distribution Function
(CDF) pnorm(q, mean, sd)
pnorm(1.96, 0, 1) Gives the area under the standard normal curve to the left of 1.96, i.e. ~0.975
QNORM Quantile Function – inverse of
pnorm qnorm(p, mean, sd)
qnorm(0.975, 0, 1) Gives the value at which the CDF of the standard normal is .975, i.e. ~1.96
Note that for all functions, leaving out the mean and standard deviation would result in default values of mean=0 and sd=1, a standard normal distribution.
pnorm students scoring higher than 84
#========== pnorm students scoring higher than 84 > pnorm(84, mean=72, sd=15.2, lower.tail=FALSE) [1] 0.21492 Answer The percentage of students scoring 84 or more in the college entrance exam is 21.5%.plot a histogram of 1000
draws from a normal distribution with mean 10, standard deviation 2. #========== plot a histogram of 1000 draws from a normal distribution with mean 10, standard deviation 2. set.seed(seed) x = rnorm(1000, 10, 2) plot(x) hist(x) Using a QQ plot. Assess the normality: qqnorm(x) qqline(x) In statistics, a Q–Q (quantile-quantile) plot is a probability plot, which is a graphical method for comparing two probability distributions by plotting their quantiles against each other. First, the set of intervals for the quantiles is chosen. A point (x, y) on the plot corresponds to one of the quantiles of the second distribution (y-coordinate) plotted against the same quantile of the first distribution (x-coordinate). Thus the line is a parametric curve with the parameter which is the number of the interval for the quantile.format leading zeros
#========== format leading zeros formatC(1, width = 2, format = "d", flag = "0") "01" formatC(125, width = 5, format = "d", flag = "0") "00125"library(pdftools)
#========== setwd("C:/Users/User/Desktop") library(pdftools) txt = pdf_text("a.pdf") str(txt) # 361 pages writeClipboard(txt[1]) Sys.setlocale(category = 'LC_ALL', 'Chinese') options("encoding" = "UTF-8") sink("war.txt") for(i in txt){ cat(i, sep="\n")} sink() shell("war.txt") txt1 = gsub(".*ORIGINATOR", "", txt) txt1 = gsub(" ", "", txt1) list = c(13:16, 19:22, 25:28, 31:34, 37:42, 45:48, 52:58, 62:68, 71:75, 78:85, 88:95, 98:105, 108:115, 118:124, 127:133, 136:142, 145:156, 159:169, 173:202, 206:221, 225:240, 244:258, 261:274, 277:290, 294:298, 302:308, 312:318, 323:331, 334:345, 348:359) txt1 = txt1[list] writeClipboard(txt1) pdf_info("a.pdf") pdf_text("a.pdf") pdf_fonts("a.pdf") pdf_attachments("a.pdf") pdf_toc("a.pdf") toc = pdf_toc("a.pdf") sink("test.txt") print(toc) sink() #========== library(pdftools) txt = pdf_text("a.pdf") str(txt) txtList = unlist(strsplit(txt, "\\s{2,}")) writeClipboard(txtList) pdftools.pdf pdftools Usage pdf_text(pdf)pdfimages
https://stackoverflow.com/questions/47133072/how-to-extract-images-from-a-scanned-pdf http://www.xpdfreader.com/pdfimages-man.html http://www.xpdfreader.com/download.html https://rdrr.io/cran/metagear/src/R/PDF_extractImages.R pdfimages a.pdf -j Quote a string to be passed to an operating system shell. Usage: shQuote(string, type = c("sh", "csh", "cmd", "cmd2")) #("PDF to PPM") files <- list.files(path = dest, pattern = "pdf", full.names = TRUE) lapply(files, function(i){ shell(shQuote(paste0("pdftoppm -f 1 -l 10 -r 300 ", i,".pdf", " ",i))) }) You could also just use the CMD prompt and type pdftoppm -f 1 -l 10 -r 300 stuff.pdf stuff.ppmOCR Extract Text from Images
download Using the Tesseract OCR engine in R library(tesseract) i = "https://upload.wikimedia.org/wikipedia/commons/thumb/9/97/Chineselanguage.svg/1200px-Chineselanguage.svg.png" chi <- tesseract("chi_sim") text <- ocr(i, engine = chi) cat(text) # In love # text <- ocr(i) # for english, default engine library(tesseract) eng <- tesseract("eng") text <- tesseract::ocr("http://jeroen.github.io/images/testocr.png", engine = eng) cat(text) results <- tesseract::ocr_data("http://jeroen.github.io/images/testocr.png", engine = eng) # list the languages have installed. tesseract_info() $datapath [1] "/Users/jeroen/Library/Application Support/tesseract4/tessdata/" $available [1] "chi_sim" "eng" "osd" chinese character recognition using Tesseract OCR download chinese trained data (it will be a file like chi_sim.traineddata) and add it to your tessdata folder. C:/Users/User/AppData/Local/tesseract4/tesseract4/tessdata/ To download the file https://github.com/tesseract-ocr/tessdata/raw/master/chi_sim.traineddata library(tesseract) chi <- tesseract("chi_sim") datapath = "C:/Users/User/Desktop/testReact/" setwd(datapath) shell(shQuote("D:/XpdfReader-win64/xpdf-tools-win-4.03/bin64/pdfimages a.pdf -j")) allFiles <- list.files(path = datapath, pattern = "jpg", full.names = TRUE) allText = character() # for(i in allFiles){ for(file in 1:5){ i = allFiles[file] cat(i, "\n") text <- tesseract::ocr(i, engine = chi) allText = c(allText, text) } setwd(datapath) Sys.setlocale(category = 'LC_ALL', 'Chinese') options("encoding" = "UTF-8") sink("result.txt") cat(allText, sep="\n") sink() options("encoding" = "native.enc") thepage = readLines("result.txt", encoding="UTF-8") thepage = gsub(" ","", thepage) sink("resultNew.txt") cat(thepage, sep="\n") sink() thepage = readLines("resultNew.txt", encoding="UTF-8") thepage = gsub("。","。\n", thepage) sink("resultNew.txt") cat(thepage, sep="\n") sink() Tesseract-OCR 實用心得 Xpdf language support packages with XpdfViewer, XpdfPrint, XpdfText First, download whichever language support package(s) you need and unpack them. You can unpack them anywhere you like – in step 3, you'll set up the config file with the path to wherever you unpacked them. Create an xpdfrc configuration file (if you haven't done this already). All of the Glyph & Cog tools read an (optional) configuration file with various global settings. To use this config file with the Windows DLLs and COM components, simply create a text file called "xpdfrc" in the same directory as the DLL, COM component, or ActiveX control. This must be a plain text file (not Word or RTF) with no file name extension (correct: xpdfrc; incorrect: xpdfrc.txt). Documentation on the configuration settings, i.e., available commands for the xpdfrc file, can be found in the documentation for the DLL or COM component. Each language support package comes with a file called "add-to-xpdfrc". You need to insert the contents of that file into your own xpdfrc file (created in step 2). This information includes pointers to the various files installed when you unpacked the language support package – make sure you modify these paths to match your install directory. The GPG/PGP key used to sign the packages is available here, or from the PGP keyservers (search for xpdf@xpdfreader.com). https://cran.r-project.org/web/packages/tesseract/vignettes/intro.html tesseract_info() to show environment remember to copy the train data to: C:/Users/william/AppData/Local/tesseract4/tesseract4/tessdata/ High Quality OCR in R Using the Tesseract OCR engine in R 實用心得 Tesseract-OCRtrain tessdata library
creating training data improve characters recognition lots of tessdata traindata tessdata langdata Traineddata Files Tesseract ocr train tessdata library on batch with lots of single character image If they are of same font, put them in a multi-page TIFF and conduct training on it. jTessBoxEditor can help you with the TIFF merging and box editing. jTessBoxEditor Here is a summary: 3. The more data, the better the OCR result, so repeat (1) and (2) until you have at least 4 pages. Limit is 32 4. Execute tesseract command to obtain the box files 5. Edit the box file using the bbTesseract editing tool 6. Execute tesseract command to generate the data files (clustering) 7. Rename files with "vie." prefix and copy the files to tessdata directory, overriding the existing data 8. Run OCR on the original images to validate your work. The accuracy rate should be in the high 90% So that the community can benefit from your work, please submit your data files. They will be posted in the VietOCR's Download page. Be sure to indicate the names of the fonts that you have trained for, so users can know which data set they should load into tessdata directory when OCRing their document. training tesseract models from scratchtesseract extra spaces in result when ocr chinese
# workaround to remove extra spaces in OCR result # https://github.com/tesseract-ocr/tesseract/issues/991, 988 and 1009 This fix can be applied via adding the following to the config file and then running combine_tessdata. preserve_interword_spaces 1 SetVariable("preserve_interword_spaces", false); these files need to be fixed: tessdata/chi_sim/chi_sim.config tessdata/chi_tra/chi_tra.config tessdata/jpn/jpn.config tessdata/tha/tha.config tessdata_best/chi_sim/chi_sim.config tessdata_best/chi_tra/chi_tra.config tessdata_best/jpn/jpn.config tessdata_best/tha/tha.config fixed tessdata_best/jpn_vert/jpn_vert.config which is included by tessdata_best/jpn/jpn.configThe name of the site environment variable R_ENVIRON
#========== The name of the site environment variable R_ENVIRON "R_HOME/etc/Renviron.site" the default is "R_HOME/etc/Rprofile.site" Sys.getenv("R_USER") Examples ## Example ~/.Renviron on Unix R_LIBS=~/R/library PAGER=/usr/local/bin/less ## Example .Renviron on Windows R_LIBS=C:/R/library MY_TCLTK="c:/Program Files/Tcl/bin" ## Example of setting R_DEFAULT_PACKAGES (from R CMD check) R_DEFAULT_PACKAGES='utils,grDevices,graphics,stats' # this loads the packages in the order given, so they appear on # the search path in reverse order. ## Example of .Rprofile options(width=65, digits=5) options(show.signif.stars=FALSE) setHook(packageEvent("grDevices", "onLoad"), function(...) grDevices::ps.options(horizontal=FALSE)) set.seed(1234) .First = function() cat("\n Welcome to R!\n\n") .Last = function() cat("\n Goodbye!\n\n") ## Example of Rprofile.site local({ # add MASS to the default packages, set a CRAN mirror old = getOption("defaultPackages"); r = getOption("repos") r["CRAN"] = "http://my.local.cran" options(defaultPackages = c(old, "MASS"), repos = r) ## (for Unix terminal users) set the width from COLUMNS if set cols = Sys.getenv("COLUMNS") if(nzchar(cols)) options(width = as.integer(cols)) # interactive sessions get a fortune cookie (needs fortunes package) if (interactive()) fortunes::fortune() }) ## if .Renviron contains FOOBAR="coo\bar"doh\ex"abc\"def'" ## then we get # > cat(Sys.getenv("FOOBAR"), "\n") # coo\bardoh\exabc"def'How to Convert Factor into Numerical?
#========== How to Convert Factor into Numerical? When you convert factors to numeric, first you should convert it into characters and then convert into numeric. as.numeric(as.character(X)) Df$column=as.numeric(as.factor(df$column) as.integer(as.factor(region))options(error=recover)
#========== options(error=recover) recover {utils} Browsing after an Error This function allows the user to browse directly on any of the currently active function calls, and is suitable as an error option. The expression options(error = recover) will make this the error option. Usage recover() When called, recover prints the list of current calls, and prompts the user to select one of them. The standard R browser is then invoked from the corresponding environment; the user can type ordinary R language expressions to be evaluated in that environment. Turning off the options() debugging mode in R options(error=NULL)Extract hyperlink from Excel file in R
#========== library(XML) # rename file to .zip my.zip.file = sub("xlsx", "zip", my.excel.file) file.copy(from = my.excel.file, to = my.zip.file) # unzip the file unzip(my.zip.file) # unzipping produces a bunch of files which we can read using the XML package # assume sheet1 has our data xml = xmlParse("xl/worksheets/sheet1.xml") # finally grab the hyperlinks hyperlinks = xpathApply(xml, "//x:hyperlink/@display", namespaces="x") To repair Hyperlink address corrupted: copy file to desk top and rename to zip file open zip file and locate: \xl\worksheets\_rels open the sheet1.xml.rels with editor remove all text: D:\Users\Lawht\AppData\Roaming\Microsoft\Excel\Extract part of a string
#========== x = c("75 to 79", "80 to 84", "85 to 89") substr(x, start = 1, stop = 2) substr(x, start, stop) x = "1234567890" substr(x, 5, 7) "567"alter grades
#========== alter grades locate the word get the line location alter the score table #========== locate the word v = c('a','b','c','e') 'b' %in% v ## returns TRUE match('b',v) ## returns the first location of 'b', in this case: 2 subv = c('a', 'f') subv %in% v ## returns a vector TRUE FALSE is.element(subv, v) ## returns a vector TRUE FALSE which() which('a' == v) #[1] 2 4 For finding all occurances as vector of indices grep() returns a vector of integers, which indicate where matches are. yo = c("a", "a", "b", "b", "c", "c") grep("b", yo) # [1] 3 4 ROC="中華民國 – 維基百科，自由的百科全書" grep("中華民國",ROC) Partial String Matching pmatch("med", c("mean", "median", "mode")) # returns 2table, cut and barplot
atab=c(1,2,3,2,1,2,3,4,5,4) table(atab)non-paste answer to concatenate two strings
capture.output(cat(counts, sep = ","))V8 is an R interface JavaScript engine.
This package helps us execute javascript code in R #Loading both the required libraries library(rvest) library(V8) #URL with js-rendered content to be scraped link = 'https://food.list.co.uk/place/22191-brewhemia-edinburgh/' #Read the html page content and extract all javascript codes that are inside a list emailjs = read_html(link) %>% html_nodes('li') %>% html_nodes('script') %>% html_text() # Create a new v8 context ct = v8() #parse the html content from the js output and print it as text read_html(ct$eval(gsub('document.write','',emailjs))) %>% html_text() info@brewhemia.co.uk Thus we have used rvest to extract the javascript code snippet from the desired location (that is coded in place of email ID) and used V8 to execute the javascript snippet (with slight code formatting) and output the actual email (that is hidden behind the javascript code). #################### Getting email address through rvest You need a javascript engine here to process the js code. R has got V8. Modify your code after installing V8 package: library(rvest) library(V8) link = 'https://food.list.co.uk/place/22191-brewhemia-edinburgh/' page = read_html(link) name_html = html_nodes(page,'.placeHeading') business_adr = html_text(adr_html) tel_html = html_nodes(page,'.value') business_tel = html_text(tel_html) emailjs = page %>% html_nodes('li') %>% html_nodes('script') %>% html_text() ct = v8() read_html(ct$eval(gsub('document.write','',emailjs))) %>% html_text()extract protected pdf document
library(pdftools) setwd("C:/Users/User/Desktop") txt = pdf_text("a.pdf") str(txt) # 361 pages # copy page 1 writeClipboard(txt[1]) # copy page 2 writeClipboard(txt[2]) # copy page 3 writeClipboard(txt[3]) Convert unicode character to string format: remove "\u" theStr = "\u9999\u6e2f\u98df\u54c1\u6295\u8d44" # "香港食品投资" ============================= Sys.setlocale(category = 'LC_ALL', 'Chinese') library(pdftools) setwd("C:/Users/User/Desktop") txt = pdf_text("45.pdf") str(txt) chi1 = gsub('\\u' , '&#x', txt[1]) chi2 = gsub('\\u' , '&#x', txt[2]) chi3 = gsub('\\u' , '&#x', txt[3]) sink("aaa.txt") cat(chi1) cat(chi2) cat(chi3) sink()Writing an R package
Develop Packages with RStudio rpackage_instructions.pdf Writing an R package from scratch Writing an R packagetable, cut and breaks
table(cut(as.numeric(resultTable[,3]), 10)) cut(as.numeric(resultTable[,3]),10) breaks = c(seq(lower, 0, by = 5), 0, seq(0, upper, by = 5)) tableA = c(1,3,5,7,9) tableB = c(1,3,5,7,2,4,6,8) tableA = c(tableA, tableB) tableA = sort(tableA) table(tableA) table(cut(tableA, 3)) breaks = c(seq(1, 3, by = 1), 4, seq(5, 9, by = 2)) table(cut(tableA, breaks))List the Files in a Directory
List the Files in a Directory/Folder list.files() list.dirs(R.home("doc")) list.dirs()best way is to run dos command from R
shell("dir /s >thisdir.txt") this will show all file details instead of only filenames in list.files() commandtest url exist
library(httr) http_error(theUrl) Load image from website download.file("url", destfile="tmp.png", mode="wb") url.exists {RCurl} return true of false With httr use url_success()download.file
This function can be used to download a file from the Internet. download.file(url, destfile, method, quiet = FALSE, mode = "w", cacheOK = TRUE, extra = getOption("download.file.extra"), headers = NULL, ...) example: destfile <- "C:/Users/User/Desktop/aaaa.jpg" url <- "https://i.pinimg.com/originals/22/2d/b8/222db84256aecf2a7532dcb1a3bab9af.jpg" download.file(url, destfile, mode = "w", method='curl') method Method to be used for downloading files. Current download methods are "internal", "wininet" (Windows only) "libcurl", "wget" and "curl", and there is a value "auto": see ‘Details’ and ‘Note’. The method can also be set through the option "download.file.method": see options(). quiet If TRUE, suppress status messages (if any), and the progress bar. mode character. The mode with which to write the file. Useful values are "w", "wb" (binary), "a" (append) and "ab". Not used for methods "wget" and "curl". See also ‘Details’, notably about using "wb" for Windows. cacheOK logical. Is a server-side cached value acceptable? extra character vector of additional command-line arguments for the "wget" and "curl" methods. headers named character vector of HTTP headers to use in HTTP requests. It is ignored for non-HTTP URLs. The User-Agent header, coming from the HTTPUserAgent option (see options) is used as the first header, automatically. ... allow additional arguments to be passed, unused.Passing arguments to R script
Passing arguments to R script Rscript --vanilla testargument.R iris.txt newname To avoid Rscript.exe loop forever for keyboard input: use this: cat("a string please: "); a = readLines("stdin",n=1);School Revision Papers
http://schoolsnetkenya.com/form-1-revision-papers-for-term-1-2018/ http://schoolsnetkenya.com/form-1-revision-papers-for-term-1-2017/ https://curriculum.gov.mt/en/Examination-Papers/Pages/list_secondary_papers.aspx http://www2.hkedcity.net/sch_files/a/hf1/hf1-lin/visitor_cabinet/67726/F1-2ndTest-Eng.pdf http://www2.hkedcity.net/sch_files/a/hf1/hf1-lin/visitor_cabinet/67726/F2-2ndTest-Eng.pdf http://www.sttss.edu.hk/parents_corner/pastpaper.phpdifference between 1L and 1
L specifies an integer type, rather than a double, it uses only 4 bytes per element the function as.integer is simplified yb "L " suffix > str(1) num 1 > str(1L) int 1Datatable
data.table FAQ ♦Data.Table Tutorial R Data.Table Tutorial Datatable Cheat Sheetsetkey does two things:
reorders the rows of the data.table DT by the column(s) provided (a, b) by reference, always in increasing order. marks those columns as key columns by setting an attribute called sorted to DT. The reordering is both fast (due to data.table's internal radix sorting) and memory efficient (only one extra column of type double is allocated). When is setkey() required? For grouping operations, setkey() was never an absolute requirement. That is, we can perform a cold-by or adhoc-by. A key is basically an index into a dataset, which allows for very fast and efficient sort, filter, and join operations. These are probably the best reasons to use data tables instead of data frames (the syntax for using data tables is also much more user friendly, but that has nothing to do with keys). library(data.table) dt=data.table(read.table("wAveTable.txt", header=TRUE, colClasses=c('character', 'numeric', 'numeric'))) colnames(dt) "Code" "WAve5" "WAve10" dt[WAve5 > 5, ] summary(dt[WAve5 = 5, ]) summary(dt[WAve5 %between% c(7,9), ]) data.table dt subset rows using i, and manipulate columns with j, grouped according to by dt[i, j, by] Create a data.table data.table(a = c(1, 2), b = c("a", "b")) convert a data frame or a list to a data.table setDT(df) or as.data.table(df) Subset data.table rows using i dt[1:2, ] subset data.table rows based on values in one or more columns dt[a > 5, ] data.table Logical Operators To Use In i >,<,<=,>=, |, !,&, is.na(),!is.na(), %in%, %like%, %between% data.table extract column(s) by number. Prefix column numbers with “-” to drop dt[, c(2)] data.table extract column(s) by name dt[, .(b, c)] create a data.table with new columns based on the summarized values of rows dt[, .(x = sum(a))] compute a data.table column based on an expression dt[, c := 1 + 2] compute a data.table column based on an expression but only for a subset of rows dt[a == 1, c := 1 + 2] compute a data.table multiple columns based on separate expressions dt[, `:=`(c = 1 , d = 2)] delete a data.table column dt[, c := NULL] convert the type of a data.table column using as.integer(), as.numeric(), as.character(), as.Date(), etc.. dt[, b := as.integer(b)] group data.table rows by values in specified column(s) dt[, j, by = .(a)] group data.table and simultaneously sort rows according to values in specified column(s) dt[, j, keyby = .(a)] summarize data.table rows within groups dt[, .(c = sum(b)), by = a] create a new data.table column and compute rows within groups dt[, c := sum(b), by = a] extract first data.table row of groups dt[, .SD[1], by = a] extract last data.table row of groups dt[, .SD[.N], by = a] perform a sequence of data.table operations by chaining multiple “[]” dt[…][…] reorder a data.table according to specified columns setorder(dt, a, -b), “-” for descending data.table’s functions prefixed with “set” and the operator “:=” work without “=” to alter data without making copies in memory df = as.data.table(df) setDT(df) extract unique data.table rows based on columns specified in “by”. Leave out “by” to use all columns unique(dt, by = c("a", "b")) return the number of unique data.table rows based on columns specified in “by” uniqueN(dt, by = c("a", "b")) rename data.table column(s) setnames(dt, c("a", "b"), c("x", "y")) data.table Syntax DT[ i , j , by], i refers to rows. j refers to columns. by refers to adding a group data.table Syntax arguments DT[ i , j , by], with, which, allow.cartesian, roll, rollends, .SD, .SDcols, on, mult, nomatch data.table fread() function to read data, mydata = fread("https://github.com/flights_2014.csv") data.table select only 'origin' column returns a vector dat1 = mydata[ , origin] data.table select only 'origin' column returns a data.table dat1 = mydata[ , .(origin)] or dat1 = mydata[, c("origin"), with=FALSE] data.table select column dat2 =mydata[, 2, with=FALSE] data.table select column Multiple Columns dat3 = mydata[, .(origin, year, month, hour)], dat4 = mydata[, c(2:4), with=FALSE] data.table Dropping Column adding ! sign, dat5 = mydata[, !c("origin"), with=FALSE] data.table Dropping Multiple Columns dat6 = mydata[, !c("origin", "year", "month"), with=FALSE] data.table select variables that contain 'dep' use %like% operator, dat7 = mydata[,names(mydata) %like% "dep", with=FALSE] data.table Rename Variables setnames(mydata, c("dest"), c("Destination")) data.table rename multiple variables setnames(mydata, c("dest","origin"), c("Destination", "origin.of.flight")) data.table find all the flights whose origin is 'JFK' dat8 = mydata[origin == "JFK"] data.table Filter Multiple Values dat9 = mydata[origin %in% c("JFK", "LGA")] data.table selects not equal to 'JFK' and 'LGA' dat10 = mydata[!origin %in% c("JFK", "LGA")] data.table Filter Multiple variables dat11 = mydata[origin == "JFK" & carrier == "AA"] data.table Indexing Set Key tells system that data is sorted by the key column data.table setting 'origin' as a key setkey(mydata, origin), 'origin' key is turned on. data12 = mydata[c("JFK", "LGA")] data.table Indexing Multiple Columns setkey(mydata, origin, dest), key is turned on. mydata[.("JFK", "MIA")] # First key 'origin' matches “JFK” second key 'dest' matches “MIA” data.table Indexing Multiple Columns equivalent mydata[origin == "JFK" & dest == "MIA"] data.table identify the column(s) indexed by key(mydata) data.table sort data using setorder() mydata01 = setorder(mydata, origin) data.table sorting on descending order mydata02 = setorder(mydata, -origin) data.table Sorting Data based on multiple variables mydata03 = setorder(mydata, origin, -carrier) data.table Adding Columns (Calculation on rows) use := operator, mydata[, dep_sch:=dep_time - dep_delay] data.table Adding Multiple Columns mydata002 = mydata[, c("dep_sch","arr_sch"):=list(dep_time - dep_delay, arr_time - arr_delay)] data.table IF THEN ELSE Method I mydata[, flag:= 1*(min < 50)] ,set flag= 1 if min is less than 50. Otherwise, set flag =0. data.table IF THEN ELSE Method II mydata[, flag:= ifelse(min < 50, 1,0)] ,set flag= 1 if min is less than 50. Otherwise, set flag =0. data.table build a chain DT[ ] [ ] [ ], mydata[, dep_sch:=dep_time - dep_delay][,.(dep_time,dep_delay,dep_sch)] data.table Aggregate Columns mean mydata[, .(mean = mean(arr_delay, na.rm = TRUE), data.table Aggregate Columns median median = median(arr_delay, na.rm = TRUE), data.table Aggregate Columns min min = min(arr_delay, na.rm = TRUE), data.table Aggregate Columns max max = max(arr_delay, na.rm = TRUE))] data.table Summarize Multiple Columns all the summary function in a bracket, mydata[, .(mean(arr_delay), mean(dep_delay))] data.table .SD operator implies 'Subset of Data' data.table .SD and .SDcols operators calculate summary statistics for a larger list of variables data.table calculates mean of two variables mydata[, lapply(.SD, mean), .SDcols = c("arr_delay", "dep_delay")] data.table Summarize all numeric Columns mydata[, lapply(.SD, mean)] data.table Summarize with multiple statistics mydata[, sapply(.SD, function(x) c(mean=mean(x), median=median(x)))] data.table Summarize by group 'origin mydata[, .(mean_arr_delay = mean(arr_delay, na.rm = TRUE)), by = origin] data.table Summary by group useing keyby= operator mydata[, .(mean_arr_delay = mean(arr_delay, na.rm = TRUE)), keyby = origin] data.table Summarize multiple variables by group 'origin' mydata[, .(mean(arr_delay, na.rm = TRUE), mean(dep_delay, na.rm = TRUE)), by = origin], or mydata[, lapply(.SD, mean, na.rm = TRUE), .SDcols = c("arr_delay", "dep_delay"), by = origin] data.table remove non-unique / duplicate cases with unique() setkey(mydata, "carrier"), unique(mydata) data.table remove duplicated setkey(mydata, NULL), unique(mydata), Note : Setting key to NULL is not required if no key is already set. data.table Extract values within a group mydata[, .SD[1:2], by=carrier], selects first and second values from a categorical variable carrier. data.table Select LAST value from a group mydata[, .SD[.N], by=carrier] data.table window function frank() dt = mydata[, rank:=frank(-distance,ties.method = "min"), by=carrier], calculating rank of variable 'distance' by 'carrier'. data.table cumulative sum cumsum() dat = mydata[, cum:=cumsum(distance), by=carrier] data.table lag and lead with shift() shift(variable_name, number_of_lags, type=c("lag", "lead")), DT = data.table(A=1:5), DT[ , X := shift(A, 1, type="lag")], DT[ , Y := shift(A, 1, type="lead")] data.table %between% operator to define a range DT = data.table(x=6:10), DT[x %between% c(7,9)] data.table %like% to find all the values that matches a pattern DT = data.table(Name=c("dep_time","dep_delay","arrival"), ID=c(2,3,4)), DT[Name %like% "dep"] data.table Inner Join Sample Data: (dt1 = data.table(A = letters[rep(1:3, 2)], X = 1:6, key = "A")), (dt2 = data.table(A = letters[rep(2:4, 2)], Y = 6:1, key = "A")), merge(dt1, dt2, by="A") data.table Left Join merge(dt1, dt2, by="A", all.x = TRUE) data.table Right Join merge(dt1, dt2, by="A", all.y = TRUE) data.table Full Join merge(dt1, dt2, all=TRUE) Convert a data.table to data.frame setDF(mydata) convert data frame to data table setDT(), setDT(X, key = "A") data.table Reshape Data dcast.data.table() and melt.data.table() data.table Calculate total number of rows by month and then sort on descending order mydata[, .N, by = month] [order(-N)], The .N operator is used to find count. data.table Find top 3 months with high mean arrival delay mydata[, .(mean_arr_delay = mean(arr_delay, na.rm = TRUE)), by = month][order(-mean_arr_delay)][1:3] data.table Find origin of flights having average total delay is greater than 20 minutes mydata[, lapply(.SD, mean, na.rm = TRUE), .SDcols = c("arr_delay", "dep_delay"), by = origin][(arr_delay + dep_delay) > 20] data.table Extract average of arrival and departure delays for carrier == 'DL' by 'origin' and 'dest' variables mydata[carrier == "DL", lapply(.SD, mean, na.rm = TRUE), by = .(origin, dest), .SDcols = c("arr_delay", "dep_delay")] data.table Pull first value of 'air_time' by 'origin' and then sum the returned values when it is greater than 300 mydata[, .SD[1], .SDcols="air_time", by=origin][air_time > 300, sum(air_time)]extract flickr image
seek .context-thumb get background-image convert _m.jpg -> _b.jpg https://live.staticflickr.com/2941/15170815109_f81b1994d2_m.jpg https://live.staticflickr.com/2941/15170815109_f81b1994d2_b.jpgR Web Scraping
get started xpath selectors R Web Scraping Rvest webscraping-using-readlines-and-rcurl xmlTreeParse, htmlTreeParse getting web data parsing xml html error in r no applicable method for xpathapply Parse and process XML (and HTML) with xml2 ================== web_page = readLines("http://www.interestingwebsite.com") web_page = read.csv("http://www.programmingr.com/jan09rlist.html") # General-purpose data wrangling library(tidyverse) # Parsing of HTML/XML files library(rvest) # String manipulation library(stringr) # Verbose regular expressions library(rebus) # Eases DateTime manipulation library(lubridate) ================== install.packages("RCurl", dependencies = TRUE) library("RCurl") library("XML") past = getURL("http://www.iciba.com/past", ssl.verifypeer = FALSE) # getURL cannot work webpage = read_html("http://www.iciba.com/past") # getURL cannot work jan09_parsed = htmlTreeParse(jan09) ================== http://www.iciba.com/past ul class="base-list switch_part" class library('rvest') library(tidyverse) url = 'http://www.iciba.com/past' webpage = readLines(url, warn=FALSE) webpage = read_html(webpage) grappedData = html_nodes(webpage,'.base-list switch_part') parseData = htmlTreeParse(webpage) rank_data = html_text(grappedData) html_node("#mw-content-text > div > table:nth-child(18)") html_table() the function htmlParse() which is equivalent to xmlParse(file, isHTML = TRUE) output = htmlParse(webpage) class(output) To parse content into an R structure : htmlTreeParse() which is equivalent to htmlParse(file, useInternalNodes = FALSE) output = htmlTreeParse(webpage) class(output) htmlTreeParse(file) especially suited for parsing HTML content returns class "XMLDocumentContent" (R data structure) equivalent to xmlParse(file, isHTML = TRUE, useInternalNodes = FALSE) htmlParse(file, useInternalNodes = FALSE) root =xmlRoot(output) xmlChildren(output) xmlChildren(xmlRoot(output)) XMLNodeList Functions for a given node Function Description xmlName() name of the node xmlSize() number of subnodes xmlAttrs() named character vector of all attributes xmlGetAttr() value of a single attribute xmlValue() contents of a leaf node xmlParent() name of parent node xmlAncestors() name of ancestor nodes getSibling() siblings to the right or to the left xmlNamespace() the namespace (if there’s one) to parse HTML tables using R sched = readHTMLTable(html, stringsAsFactors = FALSE) The html.raw object is not immediately useful because it literally contains all of the raw HTML for the entire webpage. We can parse the raw code using the xpathApply function which parses HTML based on the path argument, which in this case specifies parsing of HTML using the paragraph tag. html.raw=htmlTreeParse('http://www.dnr.state.mn.us/lakefind/showreport.html?downum=27013300', useInternalNodes=T ) html.parse=xpathApply(html.raw, "//p", xmlValue) # evaluate input and convert to text txt = htmlToText(url) ================== url = 'http://www.iciba.com/past' webpage = readLines(url, warn=FALSE) scraping_wiki = read_html(webpage) scraping_wiki %>% html_nodes("h1") %>% html_text() url = 'testvibrate.html' webpage = readLines(url, warn=FALSE) x = read_xml(webpage) xml_name(x) =========== This cannot work in office library(rvest) Sys.setlocale(category = 'LC_ALL', 'Chinese') webpage = read_html("http://www.iciba.com/haunt") ullist = webpage %>% html_nodes("ul") content = ullist[2] %>% html_text() content = gsub("n.| |\n|adj.|adv.|prep.|vt.|vi.|&","",content) content = gsub("，|；"," ",content) %>% strsplit(split = " ") %>% unlist() %>% sort() %>% unique() paste0("past","\t",capture.output(cat(content)))R scraping html text example
scraping-html-text library(rvest) scraping_wiki = read_html("https://en.wikipedia.org/wiki/Web_scraping") scraping_wiki %>% html_nodes("h1") scraping_wiki %>% html_nodes("h2") scraping_wiki %>% html_nodes("h1") %>% html_text() scraping_wiki %>% html_nodes("h2") %>% html_text() p_nodes = scraping_wiki %>% html_nodes("p") length(p_nodes) p_text = scraping_wiki %>% html_nodes("p") %>% html_text() p_text[1] p_text[5] ul_text = scraping_wiki %>% html_nodes("ul") %>% html_text() length(ul_text) ul_text[1] substr(ul_text[2], start = 1, stop = 200) li_text = scraping_wiki %>% html_nodes("li") %>% html_text() length(li_text) li_text[1:8] li_text[104:136] all_text = scraping_wiki %>% html_nodes("div") %>% html_text() body_text = scraping_wiki %>% html_nodes("#mw-content-text") %>% html_text() # read the first 207 characters substr(body_text, start = 1, stop = 207) # read the last 73 characters substr(body_text, start = nchar(body_text)-73, stop = nchar(body_text)) # Scraping a specific heading scraping_wiki %>% html_nodes("#Techniques") %>% html_text() ## [1] "Techniques" # Scraping a specific paragraph scraping_wiki %>% html_nodes("#mw-content-text > p:nth-child(20)") %>% html_text() # Scraping a specific list scraping_wiki %>% html_nodes("#mw-content-text > div:nth-child(22)") %>% html_text() # Scraping a specific reference list item scraping_wiki %>% html_nodes("#cite_note-22") %>% html_text() # Cleaning up library(magrittr) scraping_wiki %>% html_nodes("#mw-content-text > div:nth-child(22)") %>% html_text() scraping_wiki %>% html_nodes("#mw-content-text > div:nth-child(22)") %>% html_text() %>% strsplit(split = "\n") %>% unlist() %>% .[. != ""] library(stringr) # read the last 700 characters substr(body_text, start = nchar(body_text)-700, stop = nchar(body_text)) # clean up text body_text %>% str_replace_all(pattern = "\n", replacement = " ") %>% str_replace_all(pattern = "[\\^]", replacement = " ") %>% str_replace_all(pattern = "\"", replacement = " ") %>% str_replace_all(pattern = "\\s+", replacement = " ") %>% str_trim(side = "both") %>% substr(start = nchar(body_text)-700, stop = nchar(body_text)) ################ # rvest tutorials https://blog.rstudio.com/2014/11/24/rvest-easy-web-scraping-with-r/ https://blog.gtwang.org/r/rvest-web-scraping-with-r/ https://www.rdocumentation.org/packages/rvest/versions/0.3.4 https://www.datacamp.com/community/tutorials/r-web-scraping-rvest https://stat4701.github.io/edav/2015/04/02/rvest_tutorial/ https://lmyint.github.io/post/dnd-scraping-rvest-rselenium/ ################ # parse guancha library(rvest) pageHeader="https://user.guancha.cn/main/content?id=181885" pagesource = read_html(pageHeader) ################ # parse RTHK and metroradio library(rvest) pageHeader = "http://news.rthk.hk/rthk/ch/latest-news.htm" pagesource = read_html(pageHeader) className = ".ns2-title" keywordList = html_nodes(pagesource, className) html_text(keywordList) pageHeader = "http://www.metroradio.com.hk/MetroFinance/News/NewsLive.aspx" pagesource = read_html(pageHeader) className = ".n13newslist" keywordList = html_nodes(pagesource, className) className = "a" keywordList = html_nodes(keywordList, className) html_text(keywordList) ################ # parse xhamster library(rvest) pageHeader = "https://xhamster.com/users/fredlake/photos" pagesource = read_html(pageHeader) className = ".xh-paginator-button" keywordList = html_nodes(pagesource, className) html_text(keywordList) html_name(keywordList) html_attrs(keywordList) thelist = unlist(html_attrs(keywordList)) length(keywordList) as.numeric(html_text(keywordList[length(keywordList)])) pagesource %>% html_nodes(className) %>% html_text() %>% as.numeric() for ( i in keywordList ) { qlink = html_nodes(s, ".gallery-thumb") cat("Title:", html_text(qlink), "\n") qviews = html_nodes(s, "name") cat("Views:", html_text(qviews), "\n") } ################ # parse text and href pageHeader = "http://news.rthk.hk/rthk/ch/latest-news.htm" pagesource = read_html(pageHeader) className = ".ns2-title" keywordList = html_nodes(pagesource, className) className = "a" a = html_nodes(keywordList, className) html_text(a) html_attr(a, "href") ################ # extract huanqiu.com gallery pageHeader = "https://china.huanqiu.com/gallery/9CaKrnQhXac" pagesource = read_html(pageHeader) className = "article" keywordList = html_nodes(pagesource, className) className = "img" img = html_nodes(keywordList, className) html_attr(img, "src") html_attr(img, "data-alt") ################ # html_nodes samples html_nodes(".a1.b1") html_nodes(".b1:not(.a1)") # Select class contains b1 not a1: html_nodes(".content__info__item__value") html_nodes("[class='b1']") html_nodes("center") html_nodes("font") html_nodes(ateam, "center") html_nodes(ateam, "center font") html_nodes(ateam, "center font b") html_nodes("table") %>% .[[3]] %>% html_table() html_nodes("td") html_nodes() returns all nodes html_nodes(pagesource, className) html_nodes(pg, "div > input:first-of-type"), "value") html_nodes(s, ".gallery-thumb") html_nodes(s, "name") html_nodes(xpath = '//*[@id="a"]') ateam %>% html_nodes("center") %>% html_nodes("td") ateam %>% html_nodes("center") %>% html_nodes("font") td = ateam %>% html_nodes("center") %>% html_nodes("td") td %>% html_nodes("font") if (utils::packageVersion("xml2") > "0.1.2") { td %>% html_node("font") } # To pick out an element at specified position, use magrittr::extract2 # which is an alias for [[ library(magrittr) ateam %>% html_nodes("table") %>% extract2(1) %>% html_nodes("img") ateam %>% html_nodes("table") %>% `[[`(1) %>% html_nodes("img") # Find all images contained in the first two tables ateam %>% html_nodes("table") %>% `[`(1:2) %>% html_nodes("img") ateam %>% html_nodes("table") %>% extract(1:2) %>% html_nodes("img") # XPath selectors --------------------------------------------- # If you prefer, you can use xpath selectors instead of css: html_nodes(doc, xpath = "//table//td")). # chaining with XPath is a little trickier - you may need to vary # the prefix you're using - // always selects from the root node # regardless of where you currently are in the doc ateam %>% html_nodes(xpath = "//center//font//b") %>% html_nodes(xpath = "//b") read_html() html_node() # to find the first node html_nodes(doc, "table td") # to find the all node html_nodes(doc, xpath = "//table//td")) html_name() # the name of the tag html_tag() # Extract the tag names html_text() # Extract all text inside the tag html_attr() Extract the a single attribute html_attrs() Extract all the attributes # html_attrs(keywordList) this cannot use id, just list all details # html_attr(keywordList, "id") this select the ids # html_attr(keywordList, "href") this select the hrefs html_nodes("#titleCast .itemprop span") html_nodes("#img_primary img") html_nodes("div.name > strong > a") html_attr("href") html_text(keywordList, trim = FALSE) html_name(keywordList) html_children(keywordList) html_attrs(keywordList) html_attr(keywordList, "[href]", default = NA_character_) parse with xml() then extract components using xml_node() xml_attr() xml_attrs() xml_text() and xml_name() Parse tables into data frames with html_table(). Extract, modify and submit forms with html_form() set_values() submit_form(). Detect and repair encoding problems with guess_encoding() Detect text encoding repair_encoding() repair text encoding Navigate around a website as if you’re in a browser with html_session() jump_to() follow_link() back() forward() Extract, modify and submit forms with html_form(), set_values() and submit_form() The toString() function collapse the list of strings into one. html_node(":not(#commentblock)") # exclude tags ######### demos ######### # Inspired by https://github.com/notesofdabbler library(rvest) library(tidyr) page = read_html("http://www.zillow.com/homes/for_sale/....") houses = page %>% html_nodes(".photo-cards li article") z_id = houses %>% html_attr("id") address = houses %>% html_node(".zsg-photo-card-address") %>% html_text() price = houses %>% html_node(".zsg-photo-card-price") %>% html_text() %>% readr::parse_number() params = houses %>% html_node(".zsg-photo-card-info") %>% html_text() %>% strsplit("\u00b7") beds = params %>% purrr::map_chr(1) %>% readr::parse_number() baths = params %>% purrr::map_chr(2) %>% readr::parse_number() house_area = params %>% purrr::map_chr(3) %>% readr::parse_number() ################ pagesource %>% html_nodes("table") %>% .[[3]] %>% html_table() read_html(doc) %>% html_nodes(".b1:not(.a1)") # Select class contains b1 not a1: # [1] text2 use the attribute selector: read_html(doc) %>% html_nodes("[class='b1']") # [1] text2 Select class contains both: read_html(doc) %>% html_nodes(".a1.b1") # this is 'and' operation # [1] text1 combine class and ID in CSS selector div#content.sectionA # this is 'and' operation ===================== select 2 classes in 1 tag Select class contains b1 not a1: read_html(doc) %>% html_nodes(".b1:not(.a1)") use the attribute selector: read_html(doc) %>% html_nodes("[class='b1']") Select class contains both: read_html(doc) %>% html_nodes(".a1.b1") # this is 'and' operation ===================== standard CSS selector specify either or both html_nodes(".content__info__item__value, skill") # the comma is 'or' operation {xml_nodeset (4)} [1] 5h 59m 42s [2] Beginner + Intermediate [3] September 26, 2013 [4] 82,552 # has both classes in_learning_page html_nodes(".content__info__item__value.skill") # this is 'and' operation {xml_nodeset (1)} [1] Beginner + Intermediate in_learning_page %>% html_nodes(".content__info__item__value") %>% str_subset(., "viewers") h = read_html(text) h %>% html_nodes(xpath = '//*[@id="a"]') %>% xml_attr("value") html_attr(html_nodes(pg, "div > input:first-of-type"), "value") ateam %>% html_nodes("center") %>% html_nodes("td") ateam %>% html_nodes("center") %>% html_nodes("font") td = ateam %>% html_nodes("center") %>% html_nodes("td") # When applied to a list of nodes, html_nodes() returns all nodes, # collapsing results into a new nodelist. td %>% html_nodes("font") # nodes, it returns a "missing" node if (utils::packageVersion("xml2") > "0.1.2") { td %>% html_nodes("font")sort() rank() order()
x = c(1, 8,9, 4)
sort(x)
1 4 8 9
# the original position in the sorted order
rank(x)
1 3 4 2
# the sorted position in the original position
order(x)
1 4 2 3
Bioinformatics
Bioinformatics using R bioconductor Introduction to Bioconductor:Annotation and Analysis of Genomes and Genomics Assaysa list of dataframes, 3D data arrangement
d1 = data.frame(y1=c(1,2,3),y2=c(4,5,6)) d2 = data.frame(y1=c(3,2,1),y2=c(6,5,4)) d3 = data.frame(y1=c(7,8,9),y2=c(5,2,6)) mylist = list(d1, d2, d3) names(mylist) = c("List1","List2","List3") mylist[1] # same as mylist$List1 mylist[[2]][1,2] # access an element inside a dataframe mylist[[2]][2,2] # same as mylist$List2[2,2] to concate another dataframe: d4 = data.frame(y1=c(2,5,8),y2=c(1,4,7)) mylist[[4]] = d4 to create an empty list: data = list()format time string
Sys.time() sub(".* | .*", "", Sys.time()) format(Sys.time(), '%H:%M') gsub(":", "", format(Sys.time(), '%H:%M')) format(Sys.time(), '%H%M')extract 5 digit from string
activityListCode = str_replace(activityListCode, ".*\\b(\\d{5})\\b.*", "\\1")access Components of a Data Frame
access Components of a Data Frame Components of data frame can be accessed like a list or like a matrix.Accessing like a list
We can use either[, [[ or $ operator to access columns of data frame.
> x["Name"]
Name
1 John
2 Dora
> x$Name
[1] "John" "Dora"
> x[["Name"]]
[1] "John" "Dora"
> x[[3]]
[1] "John" "Dora"
Accessing with [[ or $ is similar. However, it differs for [ in that, indexing with [ will return us a data frame but the other two will reduce it into a vector.
Accessing like a matrix
Data frames can be accessed like a matrix by providing index for row and column. To illustrate this, we use datasets already available in R. Datasets that are available can be listed with the commandlibrary(help = "datasets").
We will use the trees dataset which contains Girth, Height and Volume for Black Cherry Trees.
A data frame can be examined using functions like str() and head().
> str(trees)
'data.frame': 31 obs. of 3 variables:
$ Girth : num 8.3 8.6 8.8 10.5 10.7 10.8 11 11 11.1 11.2 ...
$ Height: num 70 65 63 72 81 83 66 75 80 75 ...
$ Volume: num 10.3 10.3 10.2 16.4 18.8 19.7 15.6 18.2 22.6 19.9 ...
> head(trees,n=3)
Girth Height Volume
1 8.3 70 10.3
2 8.6 65 10.3
3 8.8 63 10.2
We can see that trees is a data frame with 31 rows and 3 columns. We also display the first 3 rows of the data frame.
Now we proceed to access the data frame like a matrix.
> trees[2:3,] # select 2nd and 3rd row
Girth Height Volume
2 8.6 65 10.3
3 8.8 63 10.2
> trees[trees$Height > 82,] # selects rows with Height greater than 82
Girth Height Volume
6 10.8 83 19.7
17 12.9 85 33.8
18 13.3 86 27.4
31 20.6 87 77.0
> trees[10:12,2]
[1] 75 79 76
We can see in the last case that the returned type is a vector since we extracted data from a single column.
This behavior can be avoided by passing the argument drop=FALSE as follows.
> trees[10:12,2, drop = FALSE]
Height
10 75
11 79
12 76
# access first row by index, returns a data.frame
x[1,]
# access first row by "name", returns a data.frame
> x["1",]
# access first row returns a vector
use as.numeric
str(as.numeric(wAveTable["1",]))
unlist which keeps the names.
str(unlist(wAveTable["1",]))
use transpose and as.vector
str(as.vector(t(wAveTable["1",])[,1]))
use only as.vector cannot convert to vector
str(as.vector(wAveTable["1",]))
# convert dataframe to matrix
data.matrix(wAveTable)
read.csv as character
wAveTable = read.csv("wAveTable.txt", sep="\t", colClasses=c('character', 'character', 'character'))frequency manipulation
grade = c("low", "high", "medium", "high", "low", "medium", "high") # using factor to count the frequency foodfac = factor(grade) summary(foodfac) max(summary(foodfac)) min(summary(foodfac)) levels(foodfac) nlevels(foodfac) summary(levels(foodfac)) # use of table to count frequency: table(grade) sort(table(grade)) table(grade)[1] max(table(grade)) summary(table(grade)) # this locate the max item: table(grade)[which(table(grade) == max(table(grade)))] # change to dataframe and find the max item: theTable = as.data.frame(table(grade)) theTable[which(theTable$Freq == max(theTable$Freq)),] # use of the count function in plyr: library(plyr) count(grade) count(mtcars, 'gear') # use of the which function: which(letters == "g") x = c(1,5,8,4,6) which(x == 5) which(x != 5)5 must have R programming tools
1) RStudio
2) lintr
If you come from the world of Python, you’ve probably heard of linting. Essentially, linting analyzes your code for readability. It makes sure you don’t produce code that looks like this: # This is some bad R codeif ( mean(x,na.rm=T)==1) { print(“This code is bad”); } # Still bad code because this line is SO long There are many things wrong with this code. For starters, the code is too long. Nobody likes to read code with seemingly endless lines. There are also no spaces after the comma in the
mean() function, or any spaces between the
== operator.
Oftentimes data science is done hastily, but linting your code is a good reminder for creating portable and understandable code.
After all, if you can’t explain what you are doing or how you are doing it, your data science job is incomplete.
lintr is an R package, growing in popularity, that allows you to lint your code.
Once you install lintr, linting a file is as easy as
lint("filename.R") .
3) Caret
Caret, which you can find on CRAN, is central to a data scientist’s toolbox in R. Caret allows one to quickly develop models, set cross-validation methods and analyze model performance all in one. Right out of the box, Caret abstracts the various interfaces to user-made algorithms and allows you to swiftly create models from averaged neural networks to boosted trees. It can even handle parallel processing. Some of the models caret includes are: AdaBoost, Decision Trees & Random Forests, Neural Networks, Stochastic Gradient Boosting, nearest neighbors, support vector machines — among the most commonly used machine learning algorithms.4) Tidyverse
You may not have heard oftidyverse as a whole, but chances are, you’ve used one of the packages in it.
Tidyverse is a set of unified packages meant to make data science…
easyr (classic R pun).
These packages alleviate many of the problems a data scientist may run into when dealing with data, such as loading data into your workspace, manipulating data, tidying data or visualizing data.
Undoubtedly, these packages make dealing with data in R more efficient.
It’s incredibly easy to get Tidyverse, you just run
install.packages("tidyverse") and you get:
ggplot2: A popular R package for creating graphics
dplyr: A popular R package for efficiently manipulating data
tidyr: An R package for tidying up data sets
readr: An R package for reading in data
purrr: An R package which extends R’s functional programming toolkit
purrr Tutorial
tibble: An R package which introduces the
tibble (tbl_df), an enhancement of the data frame
By and large, ggplot2 and dplyr are some of the most common packages in the R sphere today, and you’ll see countless posts on StackOverflow on how to use either package.
(Fine Print: Keep in mind, you can’t just load everything with
library(tidyverse)
you must load each individually!)
5) Jupyter Notebooks or R Notebooks
Data science MUST be transparent and reproducible. For this to happen, we have to see your code! The two most common ways to do this are through Jupyter Notebooks or R Notebooks. Essentially, a notebook (of either kind) allows you to run R code block by block, and show output block my block. We can see on the left that we are summarizing the data, then checking the output. After, we plot the data, then view the plot. All of these actions take place within the notebook, and it makes analyzing both output and code a simultaneous process. This can help data scientists collaborate and ease the friction of having to open up someone’s code and understand what it does. Additionally, notebooks also make data science reproducible, which gives validity to whatever data science work you do!Honorable Mention: Git
Last but not least, I want to mention Git. Git is a version control system. So why use it? Well, it’s in the name. Git allows you to keep versions of the code you are working on. It also allows multiple people to work on the same project and allows those changes to be attributed to certain contributors. You’ve probably heard of Github, undoubtedly one of the most popular git servers. You can visit my website at www.peterxeno.com and my Github at www.github.com/peterxenoR with Javascript
R and D3 tools for working with JavaScript in R R Connecting with Javascripterror handling: tryCatch
error handling code should read up on: (expr) — evaluates expression warning(…) — generates warnings stop(…) — generates errors result =Download Image
Download Image when opened with windows image viewer it also says it is corrupt. The reason for this is that you don't have specified the mode in the download.file statement. Try this: download.file(y,'y.jpg', mode = 'wb') download.file('http://78.media.tumblr.com/83a81c41926c1da585916a5c092b4789/tumblr_or0y0vdjOP1rttk8po1_1280.jpg','y.jpg', mode = 'wb') To view the image in R library(jpeg) jj = readJPEG("y.jpg",native=TRUE) plot(0:1,0:1,type="n",ann=FALSE,axes=FALSE) rasterImage(jj,0,0,1,1)testShiny
setwd("D:/KPC/testShiny") runApp("D:/KPC/testShiny")Error in file(filename, "r", encoding = encoding)
The error indicate that either the file doesn't exist or the source() command an incorrect path.call another R program from R program
source("program_B.R")to view all the functions present in a package
To list all objects in the package use ls ls("package:Hmisc") Note that the package must be attached. To list all strings lsf.str("package:dplyr") lsf.str("package:Hmisc") To see the list of currently loaded packages use search() Alternatively calling the help would also do, even if the package is not attached: help(package = dplyr) help(package = Hmisc) Finally, use RStudio which provides an autocomplete function. So, for instance, typing Hmisc:: in the console or while editing a file will result in a popup list of all dplyr functions/objects.cut2
Function like cut but left endpoints are inclusive. install.packages("Hmisc") library(Hmisc) alist = c(-15,18,2,5,4,-7,-5,-3,-1,0,2,1,5,4,6) breaks = c(-5,-3,-1,0,1,3,5) table(cut2(alist, breaks))Reference A Data Frame Column
with the double square bracket "[[]]" operator. LastDayTable[["Vol"]] or LastDayTable$Vol or LastDayTable[,"Vol"]Writing data to a file
Problem
You want to write data to a file.Solution
Writing to a delimited text file
The easiest way to do this is to usewrite.csv(). By default, write.csv() includes row names, but these are usually unnecessary and may cause confusion.
# A sample data frame
data = read.table(header=TRUE, text='
subject sex size
1 M 7
2 F NA
3 F 9
4 M 11
')
# Write to a file, suppress row names
write.csv(data, "data.csv", row.names=FALSE)
# Same, except that instead of "NA", output blank cells
write.csv(data, "data.csv", row.names=FALSE, na="")
# Use tabs, suppress row names and column names
write.table(data, "data.csv", sep="\t", row.names=FALSE, col.names=FALSE)
Saving in R data format
write.csv() and write.table() are best for interoperability with other data analysis programs. They will not, however, preserve special attributes of the data structures, such as whether a column is a character type or factor, or the order of levels in factors. In order to do that, it should be written out in a special format for R.
Below are are three primary ways of doing this:
The first method is to output R source code which, when run, will re-create the object. This should work for most data objects, but it may not be able to faithfully re-create some more complicated data objects.
# Save in a text format that can be easily loaded in R
dump("data", "data.Rdmpd")
# Can save multiple objects:
dump(c("data", "data1"), "data.Rdmpd")
# To load the data again:
source("data.Rdmpd")
# When loaded, the original data names will automatically be used.
The next method is to write out individual data objects in RDS format. This format can be binary or ASCII. Binary is more compact, while ASCII will be more efficient with version control systems like Git.
# Save a single object in binary RDS format
saveRDS(data, "data.rds")
# Or, using ASCII format
saveRDS(data, "data.rds", ascii=TRUE)
# To load the data again:
data = readRDS("data.rds")
It’s also possible to save multiple objects into an single file, using the RData format.
# Saving multiple objects in binary RData format
save(data, file="data.RData")
# Or, using ASCII format
save(data, file="data.RData", ascii=TRUE)
# Can save multiple objects
save(data, data1, file="data.RData")
# To load the data again:
load("data.RData")
An important difference between saveRDS() and save() is that, with the former, when you readRDS() the data, you specify the name of the object, and with the latter, when you load() the data, the original object names are automatically used. Automatically using the original object names can sometimes simplify a workflow, but it can also be a drawback if the data object is meant to be distributed to others for use in a different environment.
Scatter Plot of Binary Classification Dataset
Scatter Plot of Multi-Class Classification Dataset
Scatter Plot of Imbalanced Binary Classification Dataset
How we prepared the test
In order to compare the performance of CPUs vs GPUs vs TPUs for accomplishing common data science tasks, we used the tf_flowers dataset to train a convolutional neural network, and then the exact same code was run three times using the three different backends (CPUs vs GPUs vs TPUs; GPUs were NVIDIA P100 with Intel Xeon 2GHz (2 core) CPU and 13GB RAM.
TPUs were TPUv3 (8 core) with Intel Xeon 2GHz (4 core) CPU and 16GB RAM).
The accompanying tutorial notebook demonstrates a few best practices for getting the best performance out of your TPU.
For example:
Using a dataset of sharded files (e.g., .TFRecord)
Using the tf.data API to pass the training data to the TPU
Using large batch sizes (e.g. batch_size=128)
By adding these precursory steps to your workflow, it is possible to avoid a common I/O bottleneck that otherwise prevents the TPU from operating at its full potential.
You can find additional tips for optimizing your code to run on TPUs by visiting the official Kaggle TPU documentation.
How the hardware performed
The most notable difference between the three hardware types that we tested was the speed that it took to train a model using tf.keras.
The tf.keras library is one of the most popular machine learning frameworks because tf.keras makes it easy to quickly experiment with new ideas.
If you spend less time writing code then you have more time to perform your calculations, and if you spend less time waiting for your code to run, then you have more time to evaluate new ideas (Figure 2).
tf.keras and TPUs are a powerful combination when participating in machine learning competitions!
For our first experiment, we used the same code (a modified version*** of the official tutorial notebook) for all three hardware types, which required using a very small batch size of 16 in order to avoid out-of-memory errors from the CPU and GPU.
Under these conditions, we observed that TPUs were responsible for a ~100x speedup as compared to CPUs and a ~3.5x speedup as compared to GPUs when training an Xception model (Figure 3).
Because TPUs operate more efficiently with large batch sizes, we also tried increasing the batch size to 128 and this resulted in an additional ~2x speedup for TPUs and out-of-memory errors for GPUs and CPUs.
Under these conditions, the TPU was able to train an Xception model more than 7x as fast as the GPU from the previous experiment****.
The observed speedups for model training varied according to the type of model, with Xception and Vgg16 performing better than ResNet50 (Figure 4). Model training was the only type of task where we observed the TPU to outperform the GPU by such a large margin.
For example, we observed that in our hands the TPUs were ~3x faster than CPUs and ~3x slower than GPUs for performing a small number of predictions (TPUs perform exceptionally when making predictions in some situations such as when making predictions on very large batches, which were not present in this experiment).
To supplement these results, we note that Wang et.
al have developed a rigorous benchmark called ParaDnn [1] that can be used to compare the performance of different hardware types for training machine learning models.
By using this method Wang et.
al were able to conclude that the performance benefit for parameterized models ranged from 1x to 10x, and the performance benefit for real models ranged from 3x to 6.8x when a TPU was used instead of a GPU (Figure 5).
TPUs perform best when combined with sharded datasets, large batch sizes, and large models.
Price considerations when training models
While our comparisons treated the hardware equally, there is a sizeable difference in pricing. TPUs are ~5x as expensive as GPUs ($1.46/hr for a Nvidia Tesla P100 GPU vs $8.00/hr for a Google TPU v3 vs $4.50/hr for the TPUv2 with “on-demand” access on GCP ).
If you are trying to optimize for cost then it makes sense to use a TPU if it will train your model at least 5 times as fast as if you trained the same model using a GPU.
We consistently observed model training speedups on the order of ~5x when the data was stored in a sharded format in a GCS bucket then passed to the TPU in large batch sizes, and therefore we recommend TPUs to cost-conscious consumers that are familiar with the tf.data API.
Some machine learning practitioners prioritize the reduction of model training time as opposed to prioritizing the reduction of model training costs.
For someone that just wants to train their model as fast as possible, the TPU is the best choice.
If you spend less time training your model, then you have more time to iterate upon new ideas.
But don’t take our word for it — you can evaluate the performance benefits of CPUs, GPUs, and TPUs by running your own code in a Kaggle Notebook, free-of-charge.
Kaggle users are already having a lot of fun and success experimenting with TPUs and text data: check out this forum post that describes how TPUs were used to train a BERT transformer model to win $8,000 (2nd prize) in a recent Kaggle competition.
Which hardware option should you choose?
In summary, we recommend CPUs for their versatility and for their large memory capacity.
GPUs are a great alternative to CPUs when you want to speed up a variety of data science workflows, and TPUs are best when you specifically want to train a machine learning model as fast as you possibly can.
You can get better results by optimizing your code for the specific hardware that you are using and we think it would be especially interesting to compare runtimes for code that has been optimized for a GPU to runtimes for code that has been optimized for a TPU.
For example, it would be interesting to record the time that it takes to train a gradient-boosted model using a GPU-accelerated library such as RAPIDS.ai and then to compare that to the time that it takes to train a deep learning model using a TPU-accelerated library such as tf.keras.
What is the least amount of time that one can train an accurate machine learning model? How many different ideas can you evaluate in a single day? When used in combination with tf.keras, TPUs allow machine learning practitioners to spend less time writing code and less time waiting for their code to run — leaving more time to evaluate new ideas and improve one’s performance in Kaggle Competitions.
As you can see, there are many node types in an XPath tree:
You see that there is a high correlation between the sepal length and the sepal width of the Setosa iris flowers, while the correlation is somewhat less high for the Virginica and Versicolor flowers: the data points are more spread out over the graph and don’t form a cluster like you can see in the case of the Setosa flowers.
The scatter plot that maps the petal length and the petal width tells a similar story:
You see that this graph indicates a positive correlation between the petal length and the petal width for all different species that are included into the Iris data set.
Of course, you probably need to test this hypothesis a bit further if you want to be really sure of this:
# Overall correlation `Petal.Length` and `Petal.Width`
cor(iris$Petal.Length, iris$Petal.Width)
# Return values of `iris` levels
x=levels(iris$Species)
# Print Setosa correlation matrix
print(x[1])
cor(iris[iris$Species==x[1],1:4])
# Print Versicolor correlation matrix
print(x[2])
cor(iris[iris$Species==x[2],1:4])
# Print Virginica correlation matrix
print(x[3])
cor(iris[iris$Species==x[3],1:4])
You see that when you combined all three species, the correlation was a bit stronger than it is when you look at the different species separately: the overall correlation is 0.96, while for Versicolor this is 0.79.
Setosa and Virginica, on the other hand, have correlations of petal length and width at 0.31 and 0.32 when you round up the numbers.
Tip: are you curious about ggvis, graphs or histograms in particular? Check out our histogram tutorial and/or ggvis course.
After a general visualized overview of the data, you can also view the data set by entering
# Return all `iris` data
iris
# Return first 5 lines of `iris`
head(iris)
# Return structure of `iris`
str(iris)
However, as you will see from the result of this command, this really isn’t the best way to inspect your data set thoroughly: the data set takes up a lot of space in the console, which will impede you from forming a clear idea about your data.
It is therefore a better idea to inspect the data set by executing 
This is where Machine Learning comes in.
We'll keep on feeding images of a fish to a computer with the tag "fish" until the machine learns all the features associated with a fish.
Once the machine learns all the features associated with a fish, we will feed it new data to determine how much has it learned.
In other words, Raw Data/Training Data is given to the machine, so that it learns all the features associated with the Training Data.
Once, the learning is done, it is given New Data/Test Data to determine how well the machine has learned.
Let us move ahead in this Machine Learning with R blog and understand about types of Machine Learning.
Regression and Classification are some examples of Supervised Learning.
For this example, if the first observation is given the label "Man" then it is rightly classified but if it is given the label "Woman", the classification is wrong.
Similarly for the second observation, if the label given is "Woman", it is rightly classified, else the classification is wrong.
Over here, "living_area" is the independent variable and "price" is the dependent variable i.e. we are determining how does "price" vary with respect to "living_area".
Clustering is an example of unsupervised learning.
"K-means", "Hierarchical", "Fuzzy C-Means" are some examples of clustering algorithms.
In this example, the set of observations is divided into two clusters.
Clustering is done on the basis of similarity between the observations.
There is a high intra-cluster similarity and low inter-cluster similarity i.e. there is a very high similarity between all the buses but low similarity between the buses and cars.
Description of the data-set:
Prior to building any model on the data, we are supposed to split the data into "train" and "test" sets.
The model will be built on the "train" set and it's accuracy will be checked on the "test" set.
We need to load the "caTools" package to split the data into two sets.
Let's find the error by subtracting the predicted values from the actual values and add this error as a new column to the "Final_Data":
Now, we'll go ahead and calculate "Root Mean Square Error" which gives an aggregate error for all the predictions
Going ahead, let's build another model, so that we can compare the accuracy of both these models and determine which is a better one.
We'll build a new linear regression model on the "train" set but this time, we'll be dropping the ‘x' and ‘y' columns from the independent variables i.e. the "price" of the diamonds is determined by all the columns except ‘x' and ‘y'.
The model built is stored in "mod_regress2":
We see that "rmse2" is marginally less than "rmse1" and hence the second model is marginally better than the first model.
Let's split the data into "train" and "test" sets using "sample.split()" function from "caTools" package.
Let's remove the "Species" column and create a new data-set which comprises only the first four columns from the ‘iris' data-set.
The str() function gives the structure of the kmeans which includes various parameters like withinss, betweenss, etc, analyzing which you can find out the performance of kmeans.
betweenss : Between sum of squares i.e. Intracluster similarity
withinss : Within sum of square i.e. Intercluster similarity
totwithinss : Sum of all the withinss of all the clusters i.e.Total intra-cluster similarity
A good clustering will have a lower value of "tot.withinss" and higher value of "betweenss" which depends on the number of clusters ‘k’ chosen initially.
The time is ripe to become an expert in Machine Learning to take advantage of new opportunities that come your way.
This brings us to the end of this "
Below, we discuss how inheritance is carried out for the three different class systems in R programming language.
I find that for loops can be easier to plot data, partly because there is nothing to collect (or combine) at each iteration.
But there’s still repetition there. Let’s abstract that away a bit.
Suppose we want a:
1. response variable (like Rating was)
2. grouping variable (like Season was)
3. function to apply to each level
This just writes out exactly what we had before
Pulling that into a function:
Of course, in this case, if we think in terms of vectors we can actually implement random walk using implicit vectorisation:
Which reinforces one of the advantages of thinking in terms of functions: you can change the implementation detail without the rest of the program changing.
Fig: Processing time as a function of number of iterations for a simple loop.
The processing time increases linearly with the number of iterations.
Again, processing time is not extensive for the above example.
Suppose we wanted to run the example with ten thousand iterations.
We can predict how long that would take based on the linear relationship between time and iterations.
mod=lm(times~iters) #predict times
predict(mod,newdata=data.frame(iters=1e4))/60
# 45.75964
This is all well and good if we want to wait around for 45 minutes.
Running the loop in parallel would greatly decrease this time.
I want to first illustrate the problem of running loops in sequence before I show how this can done using the foreach package.
If the above code is run with 
Fig: Resources used during sequential processing of a 
Fig: Resources used during parallel processing of a 
Figure 1.
Diagram of the base class, Agent, used for the agents in a simulation.
The methods for this class are defined in the following section.
Here we define the class and its slots, and the code to define the class is given below:
######################################################################
# Create the base Agent class
#
# This is used to represent the most basic agent in a simulation.
Agent = setClass(
# Set the name for the class
"Agent",
# Define the slots
slots = c(
location = "numeric",
velocity = "numeric",
active = "logical"
),
# Set the default values for the slots.
(optional)
prototype=list(
location = c(0.0,0.0),
active = TRUE,
velocity = c(0.0,0.0)
),
# Make a function that can test to see if the data is consistent.
# This is not called if you have an initialize function defined!
validity=function(object)
{
if(sum(object@velocity^2)>100.0) {
return("The velocity level is out of bounds.")
}
return(TRUE)
}
)
Now that the code to define the class is given we can create an object whose class is Agent.
> a = Agent()
> a
An object of class "Agent"
Slot "location":
[1] 0 0
Slot "velocity":
[1] 0 0
Slot "active":
[1] TRUE
Before we define the methods for the class a number of additional commands are explored.
The first set of functions explored are the
is.object and the isS4 commands.
The is.object command determines whether or not a variable refers to an object.
The isS4
command determines whether or not the variable is an S4 object.
The reason both are required is that the isS4 command alone cannot determine if a variable is an S3 object.
You need to determine if the variable is an object and then decide if it is S4 or not.
> is.object(a)
[1] TRUE
> isS4(a)
[1] TRUE
The next set of commands are used to get information about the data elements, or slots, within an object.
The first is the slotNames
command.
This command can take either an object or the name of a class.
It returns the names of the slots associated with the class as strings.
> slotNames(a)
[1] "location" "velocity" "active"
> slotNames("Agent")
[1] "location" "velocity" "active"
The getSlots command is similar to the slotNames command.
It takes the name of a class as a string.
It returns a vector whose entries are the types associated with the slots, and the names of the entries are the names of the slots.
> getSlots("Agent")
location velocity active
"numeric" "numeric" "logical"
> s = getSlots("Agent")
> s[1]
location
"numeric"
> s[[1]]
[1] "numeric"
> names(s)
[1] "location" "velocity" "active"
The next command examined is the getClass command.
It has two forms.
If you give it a variable that is an S4 class it returns a list of slots for the class associated with the variable.
If you give it a character string with the name of a class it gives the slots and their data types.
> getClass(a)
An object of class "Agent"
Slot "location":
[1] 0 0
Slot "velocity":
[1] 0 0
Slot "active":
[1] TRUE
> getClass("Agent")
Class "Agent" [in ".GlobalEnv"]
Slots:
Name: location velocity active
Class: numeric numeric logical
The final command examined is the slot command.
It can be used to get or set the value of a slot in an object.
It can be used in place of the “@” operator.
> slot(a,"location")
[1] 0 0
> slot(a,"location") = c(1,5)
> a
An object of class "Agent"
Slot "location":
[1] 1 5
Slot "velocity":
[1] 0 0
Slot "active":
[1] TRUE
Figure 2.
Diagram of the predator and prey classes derived from the Agent class.
The first step is to create the three new classes.
######################################################################
# Create the Prey class
#
# This is used to represent a prey animal
Prey = setClass(
# Set the name for the class
"Prey",
# Define the slots - in this case it is empty...
slots = character(0),
# Set the default values for the slots.
(optional)
prototype=list(),
# Make a function that can test to see if the data is consistent.
# This is not called if you have an initialize function defined!
validity=function(object)
{
if(sum(object@velocity^2)>70.0) {
return("The velocity level is out of bounds.")
}
return(TRUE)
},
# Set the inheritance for this class
contains = "Agent"
)
######################################################################
# Create the Bobcat class
#
# This is used to represent a smaller predator
Bobcat = setClass(
# Set the name for the class
"Bobcat",
# Define the slots - in this case it is empty...
slots = character(0),
# Set the default values for the slots.
(optional)
prototype=list(),
# Make a function that can test to see if the data is consistent.
# This is not called if you have an initialize function defined!
validity=function(object)
{
if(sum(object@velocity^2)>85.0) {
return("The velocity level is out of bounds.")
}
return(TRUE)
},
# Set the inheritance for this class
contains = "Agent"
)
######################################################################
# Create the Lynx class
#
# This is used to represent a larger predator
Lynx = setClass(
# Set the name for the class
"Lynx",
# Define the slots - in this case it is empty...
slots = character(0),
# Set the default values for the slots.
(optional)
prototype=list(),
# Make a function that can test to see if the data is consistent.
# This is not called if you have an initialize function defined!
validity=function(object)
{
if(sum(object@velocity^2)>95.0) {
return("The velocity level is out of bounds.")
}
return(TRUE)
},
# Set the inheritance for this class
contains = "Bobcat"
)
The inheritance is specified using the contains option in the
setClass command.
Note that this can be a vector allowing for multiple inheritance.
We choose not to use that to keep things simpler.
If you are feeling like you need more self-loathing in your life you should try it out and experiment.
Next we define a method, move, for the new classes.
We will include methods for the Agent, Prey, Bobcat, and Lynx classes.
The methods do not really do anything but are used to demonstrate the idea of how methods are executed.
# create a method to move the agent.
setGeneric(name="move",
def=function(theObject)
{
standardGeneric("move")
}
)
setMethod(f="move",
signature="Agent",
definition=function(theObject)
{
print("Move this Agent dude")
theObject = setVelocity(theObject,c(1,2))
validObject(theObject)
return(theObject)
}
)
setMethod(f="move",
signature="Prey",
definition=function(theObject)
{
print("Check this Prey before moving this dude")
theObject = callNextMethod(theObject)
print("Move this Prey dude")
validObject(theObject)
return(theObject)
}
)
setMethod(f="move",
signature="Bobcat",
definition=function(theObject)
{
print("Check this Bobcat before moving this dude")
theObject = setLocation(theObject,c(2,3))
theObject = callNextMethod(theObject)
print("Move this Bobcat dude")
validObject(theObject)
return(theObject)
}
)
setMethod(f="move",
signature="Lynx",
definition=function(theObject)
{
print("Check this Lynx before moving this dude")
theObject = setActive(theObject,FALSE)
theObject = callNextMethod(theObject)
print("Move this Lynx dude")
validObject(theObject)
return(theObject)
}
)
There are a number of things to note.
First each method calls the
callNextMethod command.
This command will execute the next version of the same method for the previous class in the hierarchy.
Note that
I have included the arguments (in the same order) as those called by the original function.
Also note that the function returns a copy of the object and is used to update the object passed to the original function.
Another thing to note is that the methods associated with the Lync,
Bobcat, and Agent classes arbitrarily change the values of the position, velocity, and activity for the given object.
This is done to demonstrate the changes that take place and reinforce the necessity for using the callNextMethod function the way it is used here.
Finally, it should be noted that the validObject command is called in every method.
You should try adding a print statement in the validity function.
You might find that the order is a bit odd.
You should experiment with this and play with it.
There are times you do not get the expected results so be careful!
We now give a brief example to demonstrate the order that the functions are called.
In the example we create a Bobcat object and then call the move method.
We next create a Lynx object and do the same.
We print out the slots for both agents just to demonstrate the values that are changed.
> robert = Bobcat()
> robert
An object of class "Bobcat"
Slot "location":
[1] 0 0
Slot "velocity":
[1] 0 0
Slot "active":
[1] TRUE
> robert = move(robert)
[1] "Check this Bobcat before moving this dude"
[1] "Move this Agent dude"
[1] "Move this Bobcat dude"
> robert
An object of class "Bobcat"
Slot "location":
[1] 2 3
Slot "velocity":
[1] 1 2
Slot "active":
[1] TRUE
>
>
>
> lionel = Lynx()
> lionel
An object of class "Lynx"
Slot "location":
[1] 0 0
Slot "velocity":
[1] 0 0
Slot "active":
[1] TRUE
> lionel = move(lionel)
[1] "Check this Lynx before moving this dude"
[1] "Check this Bobcat before moving this dude"
[1] "Move this Agent dude"
[1] "Move this Bobcat dude"
[1] "Move this Lynx dude"
> lionel
An object of class "Lynx"
Slot "location":
[1] 2 3
Slot "velocity":
[1] 1 2
Slot "active":
[1] FALSE
Linear Regression is mainly of two types: Simple Linear Regression and Multiple Linear Regression.
Simple Linear Regression is characterized by one independent variable.
And, Multiple Linear Regression(as the name suggests) is characterized by multiple (more than 1) independent variables.
While finding the best fit line, you can fit a polynomial or curvilinear regression.
And these are known as polynomial or curvilinear regression.
Here's a coding window to try out your hand and build your own linear regression model in Python:
source: statsexchange
In the image above, you can see that population is classified into four different groups based on multiple attributes to identify ‘if they will play or not'.
To split the population into different heterogeneous groups, it uses various techniques like Gini, Information Gain, Chi-square, entropy.
The best way to understand how decision tree works, is to play Jezzball – a classic game from Microsoft (image below).
Essentially, you have a room with moving walls and you need to create walls such that maximum area gets cleared off with out the balls.
So, every time you split the room with a wall, you are trying to create 2 different populations with in the same room.
Decision trees work in very similar fashion by dividing a population in as different groups as possible.
More: Simplified Version of Decision Tree Algorithms
Now, we will find some line that splits the data between the two differently classified groups of data.
This will be the line such that the distances from the closest point in each of the two groups will be farthest away.
In the example shown above, the line which splits the data into two differently classified groups is the black line, since the two closest points are the farthest apart from the line.
This line is our classifier.
Then, depending on where the testing data lands on either side of the line, that's what class we can classify the new data as.
More: Simplified Version of Support Vector Machine
Here,
P(c|x) is the posterior probability of class (target) given predictor (attribute).
P(c) is the prior probability of class.
P(x|c) is the likelihood which is the probability of predictor given class.
P(x) is the prior probability of predictor.
Step 3: Now, use Naive Bayesian equation to calculate the posterior probability for each class.
The class with the highest posterior probability is the outcome of prediction.
KNN can easily be mapped to our real lives.
If you want to learn about a person, of whom you have no information, you might like to find out about his close friends and the circles he moves in and gain access to his/her information!
Observing the element HTML Code, we can identify that our target text is contained in the ‘a' tag.
(highlighted in blue at the lower section of the screenshot).
Next, we need to right-click on the blue highlight.
Another box with several options opens up.
Click on “Copy” which will show us new options.
There will be an XPath option also.
Click on that.
Have a look at it in the below screenshot.
Copy this XPath in any text file and check how does it look like.
The XPath copied is /html/body/div[2]/div/div/div/div[1]/div[1]/article[1]/header/h2/a.
This will bring up a laundry list of APIs, but we only need the four pertaining to YouTube (see below) and the Freebase API.
Click on the search bar and type in YouTube and you should see four options.
Enable all of them.
After clicking on the Credentials icon, you’ll need to select the OAuth client ID option.
We’re told we’re limited to 100 sensitive scope logins until the OAuth consent screen is published.
That’s not a problem for us, so we can copy the client ID and client secret
After clicking on the copy icons, we save them into two objects in RStudio (
After signing in, you’ll be asked if the YouTube application you created can access your Google account.
If you approve, click “Allow.”
This should give you a blank page with a cryptic, 
What if we need to create dummies for multiple variables in a single shot or at once?
Well, we can then create a list of all the variables for which we need dummies using c() function and pass them as arguments through select_columns.
Example:
rm(list = ls())
#install.packages('fastDummies')
library('fastDummies')
dta = read.csv("bank-loan.csv",header=TRUE)
dim(dta)
dum = dummy_cols(dta, select_columns = c('ed','default'))
dim(dum)
Output:
Here, we have created dummies for both ‘ed' and ‘default' data columns.
> dim(dta)
[1] 850 9
> dum = dummy_cols(dta, select_columns = c('ed','default'))
> dim(dum)
[1] 850 17
Now, let’s preprocess the data by removing near zero variance columns and scaling by 
The height of the bars displays how many observations occurred with each x value.
You can compute these values manually with 
You can compute this by hand by combining 
If you wish to overlay multiple histograms in the same plot, I recommend using 
There are a few challenges with this type of plot, which we will come back to in visualising a categorical and a continuous variable.
Now that you can visualise variation, what should you look for in your plots? And what type of follow-up questions should you ask? I’ve put together a list below of the most useful types of information that you will find in your graphs, along with some follow-up questions for each type of information.
The key to asking good follow-up questions will be to rely on your curiosity (What do you want to learn more about?) as well as your skepticism (How could this be misleading?).
#
smaller <- diamonds %>% filter(carat < 3)
ggplot(data = smaller, mapping = aes(x = carat)) + geom_histogram(binwidth = 0.1)
# same result:
index = diamonds$carat < 3
smaller1 <- diamonds[index,]
barplot(table(smaller1$carat))
cut_interval(smaller1$carat, n = NULL, length = 0.01)
barplot(table(cut_interval(smaller1$carat, n = NULL, length = 0.02)))
geom_freqpoly() performs the same calculation as geom_histogram(), but uses lines instead.
ggplot(data = smaller, mapping = aes(x = carat, colour = cut)) + geom_freqpoly(binwidth = 0.1)
Clusters of similar values suggest that subgroups exist in your data.
To understand the subgroups, ask:
How are the observations within each cluster similar to each other?
How are the observations in separate clusters different from each other?
How can you explain or describe the clusters?
Why might the appearance of clusters be misleading?
The histogram below shows the length (in minutes) of 272 eruptions of the Old Faithful Geyser in Yellowstone National Park.
Eruption times appear to be clustered into two groups: there are short eruptions (of around 2 minutes) and long eruptions (4-5 minutes), but little in between.
Many of the questions above will prompt you to explore a relationship between variables, for example, to see if the values of one variable can explain the behavior of another variable.
We’ll get to that shortly.
There are so many observations in the common bins that the rare bins are so short that you can’t see them (although maybe if you stare intently at 0 you’ll spot something).
To make it easy to see the unusual values, we need to zoom to small values of the y-axis with 
(
To suppress that warning, set 
However this plot isn’t great because there are many more non-cancelled flights than cancelled flights.
In the next section we’ll explore some techniques for improving this comparison.
It’s hard to see the difference in distribution because the overall counts differ so much:
To make the comparison easier we need to swap what is displayed on the y-axis.
Instead of displaying count, we’ll display 
There’s something rather surprising about this plot - it appears that fair diamonds (the lowest quality) have the highest average price! But maybe that’s because frequency polygons are a little hard to interpret - there’s a lot going on in this plot.
Another alternative to display the distribution of a continuous variable broken down by a categorical variable is the boxplot.
A 
Let’s take a look at the distribution of price by cut using 
We see much less information about the distribution, but the boxplots are much more compact so we can more easily compare them (and fit more on one plot).
It supports the counterintuitive finding that better quality diamonds are cheaper on average! In the exercises, you’ll be challenged to figure out why.
To make the trend easier to see, we can reorder 
If you have long variable names, 
The size of each circle in the plot displays how many observations occurred at each combination of values.
Covariation will appear as a strong correlation between specific x values and specific y values.
Another approach is to compute the count with dplyr:
If the categorical variables are unordered, you might want to use the seriation package to simultaneously reorder the rows and columns in order to more clearly reveal interesting patterns.
For larger plots, you might want to try the d3heatmap or heatmaply packages, which create interactive plots.
Scatterplots become less useful as the size of your dataset grows, because points begin to overplot, and pile up into areas of uniform black (as above).
You’ve already seen one way to fix the problem: using the 
But using transparency can be challenging for very large datasets.
Another solution is to use bin.
Previously you used 
Another option is to bin one continuous variable so it acts like a categorical variable.
Then you can use one of the techniques for visualising the combination of a categorical and a continuous variable that you learned about.
For example, you could bin 
Why is a scatterplot a better display than a binned plot for this case?
Patterns provide one of the most useful tools for data scientists because they reveal covariation.
If you think of variation as a phenomenon that creates uncertainty, covariation is a phenomenon that reduces it.
If two variables covary, you can use the values of one variable to make better predictions about the values of the second.
If the covariation is due to a causal relationship (a special case), then you can use the value of one variable to control the value of the second.
Models are a tool for extracting patterns out of data.
For example, consider the diamonds data.
It’s hard to understand the relationship between cut and price, because cut and carat, and carat and price are tightly related.
It’s possible to use a model to remove the very strong relationship between price and carat so we can explore the subtleties that remain.
The following code fits a model that predicts 
Once you’ve removed the strong relationship between carat and price, you can see what you expect in the relationship between cut and price: relative to their size, better quality diamonds are more expensive.
You’ll learn how models, and the modelr package, work in the final part of the book, model.
We’re saving modelling for later because understanding what models are and how they work is easiest once you have tools of data wrangling and programming in hand.
Table 1: R packages used
The functions of these R packages will be explained when the R packages are addressed.
Before you start to build the algorithm in R, you first have to install and load the libraries of the R packages.
After installation, every R script first starts with addressing the R libraries as shown below.
library(tm)
library(SnowballC)
library(topicmodels)
You can start with retrieving the dataset (or corpus for NLP).
For this experiment, we saved three text files with bug reports from three testers in a separate directory, also being our working directory (use setwd(“directory”) to set the working directory).
#set working directory (modify path as needed)
setwd(directory)
You can load the files from this directory in the corpus:
#load files into corpus
#get listing of .txt files in directory
filenames = list.files(getwd(),pattern="*.txt") #getwd() represents working directory
Read the files into a character vector, which is a basic data structure and can be read by R.
#read files into a character vector
files = lapply(filenames,readLines)
We now have to create a corpus from the vector.
#create corpus from vector
articles.corpus = Corpus(VectorSource(files))
The demo R script and demo input text file are available on my GitHub repo (please find the link in the References section).
R has a rich set of packages for Natural Language Processing (NLP) and generating plots.
The foundational steps involve loading the text file into an R Corpus, then cleaning and stemming the data before performing analysis.
I will demonstrate these steps and analysis like Word Frequency, Word Cloud, Word Association, Sentiment Scores and Emotion Classification using various plots and charts.
Plotting the top 5 most frequent words using a bar chart is a good basic way to visualize this word frequent data.
In your R script, add the following code and run it to generate a bar chart, which will display in the Plots sections of RStudio.
# Plot the most frequent words
barplot(dtm_d[1:5,]$freq, las = 2, names.arg = dtm_d[1:5,]$word,
col ="lightgreen", main ="Top 5 most frequent words",
ylab = "Word frequencies")
The plot can be seen in Figure 3.
One could interpret the following from this bar chart:
The most frequently occurring word is "good".
Also notice that negative words like "not" don't feature in the bar chart, which indicates there are no negative prefixes to change the context or meaning of the word "good" ( In short, this indicates most responses don't mention negative phrases like "not good").
"work", "health" and "feel" are the next three most frequently occurring words, which indicate that most people feel good about their work and their team's health.
Finally, the root "improv" for words like "improve", "improvement", "improving", etc.
is also on the chart, and you need further analysis to infer if its context is positive or negative
The word cloud shows additional words that occur frequently and could be of interest for further analysis.
Words like "need", "support", "issu" (root for "issue(s)", etc.
could provide more context around the most frequently occurring words and help to gain a better understanding of the main themes.
This script shows which words are most frequently associated with the top three terms (
Diagram 1: JS-R app setup
To make R talk to JS, you need three functions:
Diagram 2: JS-R communication model
The goal is to get ourselves familiar with the JS-R communication model and the whole app development process.
The full code can be found here.
Two comments to make:
To create a slider, you need to specify these three attributes: 
# Example 1. This file explores the basic mechanism for R and JS to interact.
rm(list = ls())
library(jsReact)
library(magrittr)
my_html <- create_html() %>%
add_title("Send message") %>%
add_slider(min = "Q", max = "100", oninput = "Show_value(value)") %>%
add_title("Receive message") %>%
add_div(id = "output")
my_html %<>% add_script(
"function show_value(value) {
ws.send(value) ;
}
ws.onmessage = function(msg) {
document .getELementById("output").innerHTML = msg.data;
}")
write_html_to_file(my_html, file = "inst/sample.htm1")
r_fun <- function(msg) { print(msg) }
my_app <- create_appC"inst/sample.htmL", r_fun, insert_socket = T)
start_app(my_app)
The code is divided into four sections: Html, JS, R and others.
Thanks to the 
The code is fairly self-explanatory.
You create a empty html (
This R function prints the message it gets from the JS side (which is the slider value in our case).
For a full guide to the package, check out the package guide at CRAN
here.
For now, here are some examples of pirateplots showing off some the package updates.
Here, we can see that male pirates tend to be the tallest, but there there doesn’t seem to be a difference between those who wear headbands or not, and those who have eye patches or not.
Here, I’ll use the xmen palette to plot the distribution of the weights of chickens over time (if someone has a more suitable dataset for the xmen palette let me know!):
pirateplot(formula = weight ~ Time, data = ChickWeight, main = "Weights of chickens by Time", pal = "xmen", gl.col = "gray")
mtext(text = "Using the xmen palette!", side = 3, font = 3)
mtext(text = "*The mean and variance of chicken\nweights tend to increase over time.", side = 1, adj = 1, line = 3.5, font = 3, cex = .7)
The second palette called “pony” is inspired by the Bronys in our IT department.
# Display the pony palette
piratepal(palette = "pony", trans = .1, # Slightly transparent colors plot.result = TRUE)
Here, I’ll plot the distribution of the lengths of movies as a function of their MPAA ratings (where G is for suitable for children, and R is suitable for adults)
pirateplot(formula= time ~ rating, data = subset(movies, time > 0 & rating %in% c("G", "PG", "PG-13", "R")), pal = "pony", point.o = .05, bean.o = 1, main = "Movie times by rating", bean.lwd = 2, gl.col = "gray")
mtext(text = "Using the pony palette!", side = 3, font = 3)
mtext(text = "*Movies rated for children\n(G and PG) tend to be longer \nthan those rated for adults", side = 1, adj = 1, font = 3, line = 3.5, cex = .7)
To see all of the palettes (including those inspired by movies and a transit map of Basel), just run the function with “all” as the main argument
piratepal(palette = "all")
Of course, if you find that these color palettes give you a headache, you can always set the plot to grayscale (or any other color), by specifying a single color in the palette argument. Here, I’ll create a grayscale pirateplot showing the distribution of movie budgets by their creative type:
pirateplot(formula = budget ~ creative.type, data = subset(movies, budget > 0 & creative.type %in% c("Multiple Creative Types", "Factual") == FALSE), point.o = .02, xlab = "Movie Creative Type", main = "Movie budgets (in millions) by rating", gl.col = "gray", pal = "black")
mtext("Using a grayscale pirateplot", side = 3, font = 3)
mtext("*Superhero movies tend to have the highest budgets\n...by far!", side = 1, adj = 1, line = 3, cex = .8, font = 3)
Looks like super hero movies have the highest budgets…by far!
the function f(x)= 2x.
Histogram
rnorm(6, 0 ,1) means generating 6 random values with mean 0 and standard deviation 1.
The dev.off() command is used to close the graph.
Once chart created it will save it as the test.png file.
PHP 呼叫 R 的網頁
RStudio is an open-source tool for programming in R.
If you are interested in programming with R, it’s worth knowing about the capabilities of RStudio.
It is a flexible tool that helps you create readable analyses, and keeps your code, images, comments, and plots together in one place.
In this blog post, we’ll cover some of the best features from the free version of RStudio: RStudio Desktop.
We’ve collected some of the top RStudio tips, tricks, and shortcuts to quickly turn you into an RStudio power user!
> install.packages("Dataquest")
Start learning R today with our Introduction to R course — no credit card required!
SIGN UP
We can also access documentation in the 
You can even preview the dataset before it is loaded:
And after the dataset is loaded into RStudio, you can view it with the 
Now, each time you open RStudio, you will begin with an empty session.
None of the code generated from your previous sessions will be remembered.
The R script and datasets can be used to recreate the environment from scratch.
If you would like to use 
After you hit return to select the snippet, the 
Our other favorite is the fun snippet that provides a basic template for writing a custom function.
And you can even add snippets of your own! To learn more, check out this article on code snippets from RStudio.
After the function name is selected, the inputs and code structure needed to turn your code into a function will be added automatically.
If you have a variable that you would like to extract, highlight the variable and enter 
For full details, check out this RStudio tutorial.
Now write a SQL query to select all cars from the database with a four-cylinder engine.
This command returns a dataframe that you’ll save as 
The dataframe looks like this:
Specify 
If you want to see the SQL code that this command was converted to, you can use the 
When you’re satisfied with your query results, you use the 
There you have it! Three different approaches to querying a SQL database with similar results.
The only difference between the examples is that the 
You will then see the set of available editor commands (commands that affect the current document's contents, or the current selection), alongside RStudio commands (commands whose actions are scoped beyond just the current editor).
Each row represents a particular command binding -- the command's 
If you made a mistake, you can press 
After you've updated the bindings to your liking, click 
If two commands are bound to the same key sequence, then that conflict will be highlighted and displayed in yellow.
In this particular example, you can calculate the height of a child if you know her age:
$\text{Height} = a + \text{Age} * b$
In this case, “a” and “b” are called the intercept and the slope respectively.
With the same example, “a” or the intercept, is the value from where you start measuring.
Newborn babies with zero months are not zero centimeters necessarily; this is the function of the intercept.
The slope measures the change of height with respect to the age in months.
In general, for every month older the child is, his or her height will increase with “b”.
A linear regression can be calculated in R with the command 
As you might notice already, looking at the number of siblings is a silly way to predict the height of a child.
Another aspect to pay attention to your linear models is the p-value of the coefficients.
In the previous example, the blue rectangle indicates the p-values for the coefficients age and number of siblings.
In simple terms, a p-value indicates whether or not you can reject or accept a hypothesis.
The hypothesis, in this case, is that the predictor is not meaningful for your model.
The p-value for age is 4.34*e-10 or 0.000000000434.
A very small value means that age is probably an excellent addition to your model.
The p-value for the number of siblings is 0.85.
In other words, there’s 85% chance that this predictor is not meaningful for the regression.
A standard way to test if the predictors are not meaningful is looking if the p-values smaller than 0.05.
The idea in here is that the sum of the residuals is approximately zero or as low as possible.
In real life, most cases will not follow a perfectly straight line, so residuals are expected.
In the R summary of the lm function, you can see descriptive statistics about the residuals of the model, following the same example, the red square shows how the residuals are approximately zero.
However, it’s essential to keep in mind that sometimes a high R² is not necessarily good every single time (see below residual plots) and a low R² is not necessarily always bad.
In real life, events don’t fit in a perfectly straight line all the time.
For example, you can have in your data taller or smaller children with the same age.
In some fields, an R² of 0.5 is considered good.
With the same example as above, look at the summary of the linear model to see its R².
In the blue rectangle, notice that there’s two different R², one multiple and one adjusted.
The multiple is the R² that you saw previously.
One problem with this R² is that it cannot decrease as you add more independent variables to your model, it will continue increasing as you make the model more complex, even if these variables don’t add anything to your predictions (like the example of the number of siblings).
For this reason, the adjusted R² is probably better to look at if you are adding more than one variable to the model, since it only increases if it reduces the overall error of the predictions.
If you see the summary of your new model, you can see that it has pretty good results (look at the R²and the adjusted R²)
This can be a problem.
If you have more data, your simple linear model will not be able to generalize well.
In the previous picture, notice that there is a pattern (like a curve on the residuals).
This is not random at all.
What you can do is a transformation of the variable.
Many possible transformations can be performed on your data such as adding a quadratic term $(x^2)$, a cubic $(x^3)$ or even more complex such as ln(X), ln(X+1), sqrt(X), 1/x, Exp(X).
The choice of the correct transformation will come with some knowledge of algebraic functions, practice, trial, and error.
Let’s try with a quadratic term.
For this, add the term “I” (capital "I") before your transformation, for example, this will be the normal linear regression formula:
Now you don’t see any clear patterns on your residuals, which is good!
It is however important to understand the relation between them in order to be able to pick the right metric in a given scenario.
Gini Impurity is usually easier to compute than Entropy (involves logarithmic calculations )
Denominator of Precision is a variable : i.e the False positives (negative samples classified as positive) can vary each time in our solution.
Denominator of Recall is a constant : It represents the total Positive cases (i.e True Positive + False Negative) & hence will remain fixed throughout.
Area under this curve (called AUC) and can also be used as a performance metric.
The higher the AUC, better the model.
For a linearly separable data, we can either have a linear regression or logistic regression.
However the decision boundary cries for both.
where,
ht(x) — output of the weak classifier ‘t'
αt — weight applied to classifier ‘t'
Here, εt is the error term of the classifier ‘t'
Interpretation:
The Empirical Rule states that 99.7% of data observed following a normal distribution lies within 3 standard deviations of the mean.
Under this rule, 68% of the data falls within one standard deviation, 95% percent within two standard deviations, and 99.7% within three standard deviations from the mean.
The T distribution (also called Student's T Distribution) is a family of distributions that look almost identical to the normal distribution curve, only a bit shorter and wider/fatter.
We use the t distribution instead of the normal distribution, when we have smaller samples.
The larger the sample size, the more the t distribution looks like the normal distribution.
In fact, after 30 samples, the T- distribution is almost exactly like the normal distribution.
To view all the built-in color names which R knows about (n = 657), use the following R code :
There are 3 types of palettes : sequential, diverging, and qualitative.
Sequential palettes are suited to ordered data that progress from low to high (gradient).
The palettes names are : Blues, BuGn, BuPu, GnBu, Greens, Greys, Oranges, OrRd, PuBu, PuBuGn, PuRd, Purples, RdPu, Reds, YlGn, YlGnBu YlOrBr, YlOrRd.
Diverging palettes put equal emphasis on mid-range critical values and extremes at both ends of the data range.
The diverging palettes are : BrBG, PiYG, PRGn, PuOr, RdBu, RdGy, RdYlBu, RdYlGn, Spectral
Qualitative palettes are best suited to representing nominal or categorical data.
They not imply magnitude differences between groups.
The palettes names are : Accent, Dark2, Paired, Pastel1, Pastel2, Set1, Set2, Set3
You can also view a single RColorBrewer palette by specifying its name as follow :
Use the palettes as follow :
Missing values (“TRUE”) by feature column
We see that missingness (is.na) is FALSE for all columns, which is great.
The second line in the code chunk above tells us that there is one duplicated record, which we remove in line 3.
Histograms of all our variables in our dataset
Here you see the univariate distributions – univariate because you’re looking at one variable at a time, not at bivariate correlations at this point.
The most important variable to look at here is our Y, labeled “target”.
This is coded 1 for patients with heart disease, and 0 for healthy patients.
You can see that the data are not heavily imbalanced, there are only a few more healthy patients than patients with the disease.
Regarding the other categorical features such as sex, you can also see from this graph whether they are imbalanced (1 = male, 0 = female here).
With regard to the continuous variables such as age or maximum heart rate achieved (“thalach”), you can visually check whether they are more or less normally distributed (such as age), or whether they exhibit some distribution that might need some form of normalization or discretization (e.g., “oldpeak”).
You can also see that there are some categorical variables where certain values are represented only rarely among the patients.
For instance, “restecg” has only very few instances of values == 2.
We will deal with this in a minute.
Let’s move on to bivariate statistics.
We are plotting a correlation matrix, in order to
a) check if we have features that are highly correlated (which is problematic for some algorithms), and
b) get a first feeling about which features are correlated with the target (heart disease) and which are not:
cormat = cor(dat %>% keep(is.numeric))
cormat %>% as.data.frame %>% mutate(var2=rownames(.)) %>%
pivot_longer(!var2, values_to = "value") %>%
ggplot(aes(x=name,y=var2,fill=abs(value),label=round(value,2))) +
geom_tile() + geom_label() + xlab(") + ylab(") +
ggtitle("Correlation matrix of our predictors") +
labs(fill="Correlation\n(absolute):")
This prints the following correlation matrix:
Correlation matrix of all our predictor variables
We can see that aside from the diagonal (correlation of a variable with itself, which is 1), we have no problematically strong correlations between our predictors (strong meaning greater than 0.8 or 0.9 here).
If you have many, many features and don’t want to look at 1,000 by 1,000 correlation matrices, you can also print a list of all correlations that are greater than, say, 0.8 with the following code:
highcorr = which(cormat > .8, arr.ind = T)
paste(rownames(cormat)[row(cormat)[highcorr]],
colnames(cormat)[col(cormat)[highcorr]], sep=" vs. ") %>%
cbind(cormat[highcorr])
Now let’s look at the bivariate relations between the predictors and the outcome.
For continuous predictors and a dichotomous outcome (heart disease or no heart disease), box plots are a good way of visualizing a bivariate association:
dat %>% select(-c(sex,cp,ca,thal,restecg,slope,exang,fbs)) %>%
pivot_longer(!target, values_to = "value") %>%
ggplot(aes(x=factor(target), y=value, fill=factor(target))) +
geom_boxplot(outlier.shape = NA) + geom_jitter(size=.7, width=.1, alpha=.5) +
scale_fill_manual(values=c("steelblue", "orangered1")) +
labs(fill="Heart disease:") +
theme_minimal() +
facet_wrap(~name, scales="free")
I’ve de-selected the non-continuous variables in the first line here manually, because I haven’t transformed the categorical variables (say, “ca” or “restecg”) into factors yet, which is of course a bit lazy, but if you have hundreds of features, there are of course more flexible ways to keep only continuous variables for the following plot.
This is what the chunk above returns:
Boxplots of the associations between our continuous predictors and the outcome
You can read this graph as follows: With regard to age, the patients with heart disease (red box) are on average older compared with the patients without heart disease (blue box).
The thick horizontal line within each box denotes the median.
The box encompasses 50% of all cases (i.e.
from the 25 percentile to the 75 percentile).
The jitter points show you where all of the patients are located within each group.
So you see that, yes, heart disease patients are typically older, but you also have a couple of patients younger than 50 in the dataset who have coronary heart disease, and of course many older ones that are healthy.
But comparing the medians, you can see that age, oldpeak, and thalach are better predictors of heart disease compared with chol or trestbps, where the median values are almost equal in both groups.
For our categorical variables, we just use simple stacked barplots to show the differences between healthy and sick patients:
dat %>% select(sex,cp,ca,thal,restecg,slope,exang,fbs,target) %>%
pivot_longer(!target, values_to = "value") %>%
ggplot(aes(x=factor(value), fill=factor(target))) +
scale_fill_manual(values=c("steelblue", "orangered1")) +
geom_bar(position="fill", alpha=.7)+
theme_minimal() +
labs(fill="Heart disease:") +
facet_wrap(~name, scales="free")
Which gives us:
Associations between our categorical predictors and the outcome
Again, you can see at a glance that “fbs” is obviously not a strong predictor of heart disease, whereas “exang” definitely is.
We also see that males are overrepresented in sick patients compared with females.
So far, we have used very simple means to visualize the data.
In my experience, in many applied business use cases, you already know most of what you wanted to know at this stage! A few simple descriptive graphs and indicators most often show you what are the most important predictors, what are the important sub-groups you need to focus on in greater detail, where do you have outliers or a lot of missing data which distorts the overall picture, and so on.
Often, the complicated algorithm later on only confirms what we have seen so far.
So it’s important not to skip this step and always do visual und descriptive inspection of your data.
You might ask, if in many cases bar charts and correlation coefficients is all we need to understand what is going on, why do we need the complicated machine-learning part? That is because while 80% of the explanation is often simple and can be inferred from looking at a graph or table, the other 20% is more complicated and requires domain knowledge and/or more sophisticated statistical analysis.
Our example here perfectly illustrates this point: Older people suffer from heart disease more often than younger people; men are much more likely to get it compared with females; these findings are trivial and everyone can see that from the graphs, no PhD in statistics required.
And these simple associations can already guide clinical practice to a significant degree.
You’re an older male? You’re in a risk group.
You’re a young female? You’re probably fine.
However, there are many more complex relationships at work.
For instance, females often present themselves with different forms of chest pain compared with males.
This is an example for an interaction effect that you couldn’t easily infer from the bivariate graphs above.
Non-linearity, interaction effects, spurious correlations caused by third variables and multi-collinearity, complex data structures (e.g.
time series, nested data…) – these are examples of aspects that cause descriptive inspections to be insufficient when we not only want to find out the most obvious things (older people are more at risk than young people), but also want to get behind the more complex relations.
Decision tree on our training data
How do you read this tree? Starting from the top, the most important feature that can split the data in two most dissimilar subsets (with regard to how often heart disease occurs) is “thal2”, i.e.
wether the patient has a normal blood flow as opposed to a defect from a blood disorder called thalassemia.
If the patient has a normal blood flow (value > 0 , i.e.
1), then we continue to the right branch of the tree, if not, continue to the left.
If the blood flow is normal, then the next most important variable is “thalach”, i.e.
the maximum heart rate achieved during exercise.
You can see that if this is greater than 155 bpm, then we continue to the right where we then check for “ca1”, i.e.
whether one major vessel was colored by flouroscopy.
I’m just pretending here to understand what any of this means, but recall the bar chart above where we saw that 0 vessels colored by flouroscopy was associated with the lowest proportion of patients with coronary heart disease, whereas those with 1, 2 or 3 colored vessels were predominantly diagnosed with heart disease.
In our tree, if ca1 == 0, i.e.
not one major vessel colored, we continue to the left where we reach the end note 14 (second bar from the right).
What do the bars on the bottom of the chart mean? They show the proportion of patients in each bucket with (light grey) vs. without (dark grey) heart disease.
Meaning that end node 14 (second bar from the right) is the group of patients with the lowest risk of having coronary heart disease.
Thus, our algorithm here finds that if you:
don’t have thalassemia,
can achieve a heart rate of more than 155 bpm while exercising, and
don’t have one major vessel colored by flouroscopy,
then we predict “no heart disease” with a 98% probability (i.e.
the proportion of healthy patients in the respective bucket).
If, by contrast, you do have thalassemia, there are 1 or more colored major vessels, and your chest pain type (cp) is not “2” (2 standing for non-anginal pain), then the algorithm predicts “heart disease” with a high confidence.
You can also see that there are several end node buckets (e.g., node 5, node 15) which are quite mixed.
Patients with these combinations of features are not well understood by the algorithm and the predictions are often wrong for these groups.
Now, recall what we discussed about overfitting: Of course we could go into these groups and find more features that separated the healthy from the sick patients.
In fact, if you set the values “minsplit” (minimum number of cases separated at a split) to 1, “minbucket” (minimum number of patients in an endnote) to 0, and “mincriterion” to a small value (p-value to determine if a split is significant), you get a vastly overfitted tree.
Let’s try it out:
set.seed(2022)
tree2 = party::ctree(y_train ~ ., data=cbind(x_train, y_train),
controls = ctree_control(minsplit=1, mincriterion = 0.01,minbucket = 0))
plot(tree2)
Overfitted tree
As you can see, just like we discussed above when we were warning against the dangers of overfitting, the algorithm has come up with very specific rules that often only apply to 2 or 3 people in the dataset.
For instance, if the maximum heart rate achieved is above 109, but below 144, and the patient is male, older than 59 and does not suffer from thalassemia, the algorithm always predicts heart disease.
You can see why this type of algorithm would perform poorly with new, unseen data.
We would want to “prune” this tree of nodes that introduce classification rules that are too idiosyncratic/specific to the training data.
But of course we don’t want to prune nodes that reflect true causal relations, i.e.
the actual data-generating process (which is obviously unknown to us).
Thus, the challenge in any machine-learning model is to get an algorithm that classifies the data with as specific rules as necessary, but without getting too specific and overfit to the training data.
In your real-world application, of course, you don’t grow the second (overfitted) tree, but you can use the first one for presentation slides and as a benchmark for the models which we are about to train.
Summary of model training (random forest)
What does this mean? From the summary we can verify that we set up our dataset correctly.
There are 27 features, 211 patients, and two outcomes (1 = heart disease, 0 = no heart disease).
Then you see that five values were tried for the hyperparameter “mtry”.
With each of the values, 5-fold cross validation was performed.
If you look at the accuracy values, you can see that mtry = 10 worked best with our data.
On average (of the five runs using cross-validation), 82.4% of the validation sets (= 42 patients during each run) were classified correctly.
Although this accuracy was obtained with a train/validation split, we still have yet to judge the final evaluation score of the algorithm against the unseen test dataset, because all the patients in the training data were used to train the model at some point, so technically it’s not an “out-of-sample” accuracy.
But before the final evaluation, we want to try out a few more algorithms.
With a random forest, you can obtain a feature importance plot which tells you which of the variables were most often used as the important splits at the top of the trees.
Just run:
plot(varImp(mod), main="Feature importance of random forest model on training data")
Feature importance plot of our random forest
You can see that, unlike our single decision tree on all of the training data, where “thal2” was the most important feature before “ca0”, across an ensemble of 500 different trees, it’s actually “ca0” (= zero major vessels colored by flouroscopy, whatever that means) that ends up the most important predictor, tied with “cp0” (chest pain type 0 = asymptomatic).
Recall that a machine-learning model tuned for prediction such as a random forest cannot be interpreted as revealing causal associations between the predictors and the outcome.
Nevertheless, it can guide clinical practice knowing which features are the most useful for predicting heart disease.
This best works when enriched with domain knowledge about mechanisms and causality.
Next, let’s try out a neural network, simply because many of you will probably associate machine-learning or artificial intelligence in general with artificial neural networks, or deep learning.
In general, it is true that neural networks outperform all other machine-learning algorithms when it comes to the classification of abstract data such as images or videos.
For a more detailed tutorial about how you can build a deep learning algorithm in R, see here.
In cases with classical flat data files such as ours, on the other hand, other algorithms often work equally well or better.
Here, let’s use a simple network with as few lines of code as necessary:
set.seed(2022)
mod2 = caret::train(x_train, y_train, method="avNNet",
preProcess = c("center", "scale", "nzv"),
tuneGrid = expand.grid(size = seq(3,21,by=3), decay=c(1e-03, 0.01, 0.1,0),bag=c(T,F)),
trControl = trainControl(method="cv", number=5, verboseIter = T),
importance=T)
mod2
Here we use the pre-processing steps of centering and scaling the data because, as noted above, neural networks are optimized more easily if the features have similar numerical ranges, instead of, say, maximum heart rate being in the range of 140-200 whereas other features having values bounded by 0 and 1.
Near-zero variance (“nzv”) means that we disregard features where almost all patients have the same value.
Tree-based methods such as random forests are not as sensitive to these issues.
We have a few more tuning parameters here.
“Size” refers to the number of nodes in the hidden layer.
Our network has an input layer of 27 nodes (i.e.
the number of features) and an output layer with one node (the prediction of 1 or 0) and in between, a hidden layer where interactions between the features and non-linear transformations can be learned.
As with other hyperparameters, the optimal size of the hidden layer(s) depend on the data at hand, so we just try out different values.
Decay is a regularization parameter that causes the weights of our nodes to decrease a bit after each round of updating the values after backpropagation (i.e.
the opposite of what the learning rate does wich is used in other implementations of neural networks).
What this means is, roughly speaking, we don’t want the network to learn too ambitiously with each step of adapting its parameters to the evidence, in order to avoid overfitting.
Anyway, as you can see from the code, we have passed 7 different values for “size” to consider, 4 values for “decay”, and two for “bag” (true or false, specifying how to aggregate several networks’ predictions with various random number seeds, which is what the avNNet classifier does, bagging = bootstrap aggregating), so we have 7*4*2 = 56 combinations to try out.
The result:
Output of neural network training
Thus, our best-performing model yields 85.3% accuracy, which is a slight improvement over the random forest.
Again, we can look at a feature importance plot:
plot(varImp(mod2), main="Feature importance of neural network classifier on training data")
Feature importance plot with neural network
It’s slightly different than the plot before, but the top five features are the same, just in a different order.
Note that with “unstable” methods such as neural networks, if you run the same code 10 times, you can end up with ten (slightly) different feature importance lists, but the general pattern of which features are important and which aren’t will be the same.
Let’s try out one last algorithm.
The popular “(extreme) gradient boosted machines” (xgboost) work similar to a random forest, except they proceed sequentially: A first tree is grown, then more weight is put on the badly predicted samples before the next tree is grown.
As a result, in many cases, xgboost outperforms random forests.
Let’s see if this is the case here as well:
set.seed(2022)
mod3 = caret::train(x_train, y_train, method="xgbTree",
tuneGrid = expand.grid(nrounds=c(50,100),max_depth=c(5,7,9),
colsample_bytree=c(0.8,1),subsample=c(0.8,1),
min_child_weight=c(1,5,10),eta=c(0.1,0.3),gamma=c(0,0.5)),
trControl = trainControl(method="cv", number=5, verboseIter = T))
mod3
plot(varImp(mod3), main="Feature importance of XGBoost model on training data")
Here we have more tuning parameters compared with the random forest; I just inserted a few values that I deemed plausible into the tuning grid, but if you want to do serious hyperparameter tuning, you can of course spend a bit more time here determining which combination of parameters works best.
Xgboost is in general quite fast so even though we try out 2*3*2*2*3*2*2 = 288 parameter combinations, running this code should only take a minute at most even on a local machine.
Which means that you could tune even more.
Compare the performance of the three algorithms:
results = data.frame(Model = c(mod$method,mod2$method, mod3$method),
Accuracy = c(max(mod$results$Accuracy), max(mod2$results$Accuracy), max(mod3$results$Accuracy)))
results %>% ggplot(aes(x=Model, y=Accuracy, label=paste(round(100*Accuracy,1),"%"))) +
geom_col(fill="steelblue") + theme_minimal() + geom_label() +
ggtitle("Accuracy in the training data by algorithm")
Comparison of our used algorithms during model training
The neural network actually performed slightly better than the xgboosted tree, although the values are quite similar and if you repeat the model training a couple of times, you might get different results.
With use cases like this, I prefer to go with tree-based models such as random forests or xgboost over neural networks because with the former, I can understand better how the algorithm arrives at its predictions (see our example tree in the previous section).
You could also, of course, visualize a neural network with all the weights obtained during training displayed next to the nodes, it’s not alchemy, but it’s less easily interpreted when you want to reconstruct how the network processes a certain patient.
Anyways, let’s decide at this point that our neural network (“mod2”) was the best model and we want to move forward with it.
Confusion matrix and summary statistics of our predictions on the test set
As you can see, our out-of-sample predictive accuracy was 87.9%.
The confusion matrix tells us that 40 patients with heart disease were correctly classified, and 40 healthy patients were also correctly classified, but there were 3 patients where our model thought they had heart disease but in reality they didn’t, and, conversely, we overlooked coronary heart disease in 8 patients.
In addition to accuracy, other metrics are often used to evaluate the goodness of a machine-learning algorithm.
Keep in mind that our sample was balanced (47% have heart disease, 53% don’t), whereas in many other use cases, you often have a severe class imbalance (e.g., 99% of customers won’t buy, 1% do buy, or 99% of patients won’t die vs.
1% die), so “99% accuracy” is useless to you as an indicator in these cases.
You can resort to using sensitivity/specificity which are also given in the output (specificity = how many of the true positive cases are detected, which is a useful indicator if the positive cases are rare, and specificity = how many true negatives are correctly classified).
Which of these metrics is more important to you depends on your case, i.e.
your cost function.
In this case, I’d say it’s better to detect all true cases who have the disease, and we can live with a few false positives, so I’d look at sensitivity rather than specificity.
In other cases, you want to avoid many false positives (e.g., spam detection, it’s much more annoying if many of your important work e-mails disappear in the spam folder), so sensitivity is maybe more important.
In addition to these metrics, you also often find precision (proportion of true positive predictions relative to all “positive” predictions), and recall (proportion of true positive predictions relative to all actual positives), and F1 (harmonic mean of precision and recall).
You can get these as well with
precision(predictions, y_test)
recall(predictions, y_test)
F_meas(predictions, y_test)
Prediction of a new patient (binary and with probabilities)
The first command just predicts yes or no.
For this new patient, we predict “no heart disease”.
With the second command, we also get the probabilities to belong in each class.
We see that the new patient has a 86% probability of being healthy and a 13.9% probability of having coronary heart disease according to our algorithm.
Especially with new data I find it helpful to get the predicted probabilities to get a sense for how certain the algorithm is in assigning this prediction.
As mentioned earlier, Random forest works on the Bagging principle.
Now let’sdive in and understand bagging in detail.
Now let’s look at an example by breaking it down with the help of the following figure.
Here the bootstrap sample is taken from actual data (Bootstrap sample 01, Bootstrap sample 02, and Bootstrap sample 03) with a replacement which means there is a high possibility that each sample won’t contain unique data.
The model (Model 01, Model 02, and Model 03) obtained from this bootstrap sample is trained independently.
Each model generates results as shown.
Now the Happy emoji has a majority when compared to the Sad emoji.
Thus based on majority voting final output is obtained as Happy emoji.
For example: consider the fruit basket as the data as shown in the figure below.
Now n number of samples are taken from the fruit basket, and an individual decision tree is constructed for each sample.
Each decision tree will generate an output, as shown in the figure.
The final output is considered based on majority voting.
In the below figure, you can see that the majority decision tree gives output as an apple when compared to a banana, so the final output is taken as an apple.
# checking the oob score
classifier_rf.oob_score_
grid_search.best_score_
rf_best = grid_search.best_estimator_
rf_best
From hyperparameter tuning, we can fetch the best estimator, as shown.
The best set of parameters identified was max_depth=5, min_samples_leaf=10,n_estimators=10
from sklearn.tree import plot_tree
plt.figure(figsize=(80,40))
plot_tree(rf_best.estimators_[7], feature_names = X.columns,class_names=['Disease', "No Disease"],filled=True);
The trees created by estimators_[5] and estimators_[7] are different.
Thus we can say that each tree is independent of the other.
imp_df = pd.DataFrame({
"Varname": X_train.columns,
"Imp": rf_best.feature_importances_
})
imp_df.sort_values(by="Imp", ascending=False)
A caption
=========================
Github files; https://github.com/ghettocounselor
Useful PDF for common questions in Lectures;
https://github.com/ghettocounselor/Machine_Learning/blob/master/Machine-Learning-A-Z-Q-A.pdf
Random Forests tutorials
Decision Trees
Decision tree learning
classification and regression trees for machine learning
classification and regression trees
Random Forests
Random Forest In R
Random Forest
You start at the root node (depth 0 over 3, the top of the graph):
Each row in a confusion matrix represents an actual target, while each column represents a predicted target.
The first row of this matrix considers dead passengers (the False class): 106 were correctly classified as dead (
It is the proportion of true positive and true negative over the sum of the matrix.
With R, you can code as follow:
accuracy_Test = sum(diag(table_mat)) / sum(table_mat)
Table 1 visualizes the output of the RStudio console that got returned by the previous code and illustrates that our example data is composed of two numerical columns.
This data frame will be used to train our model.
Next, we also have to create a test data frame:
As shown in Table 2, the previous R programming syntax has created another data frame that contains only the predictor column x.
You might have noticed this output is not exact URL data that’s because it is raw data.
From the plot we can see each of the 50 states represented in a simple two-dimensional space.
The states that are close to each other on the plot have similar data patterns in regards to the variables in the original dataset.
We can also see that the certain states are more highly associated with certain crimes than others. For example, Georgia is the state closest to the variable Murder in the plot.
If we take a look at the states with the highest murder rates in the original dataset, we can see that Georgia is actually at the top of the list:
Heatmap Data Visualization
Website Heatmaps
Grid Heatmaps
Clustered Heatmaps
In the plot above, high values are in red and low values are in yellow.
It’s possible to specify a color palette using the argument col, which can be defined as follow:
Using custom colors:
Other arguments can be used including:
labRow, labCol
hclustfun: hclustfun=function(x) hclust(x, method=“ward”)
In the R code above, the bluered() function [in gplots package] is used to generate a smoothly varying set of colors.
You can also use the following color generator functions:
colorpanel(n, low, mid, high)
n: Desired number of color elements to be generated
low, mid, high: Colors to use for the Lowest, middle, and highest values.
mid may be omitted.
redgreen(n), greenred(n), bluered(n) and redblue(n)
Arguments are available for changing the default clustering metric (“euclidean”) and method (“complete”).
It’s also possible to annotate rows and columns using grouping variables.
The d3heamap() function makes it possible to:
Put the mouse on a heatmap cell of interest to view the row and the column names as well as the corresponding value.
Select an area for zooming.
After zooming, click on the heatmap again to go back to the previous display
Additional arguments:
Note that, split can be also a data frame in which different combinations of levels split the rows of the heatmap.
It’s also possible to combine km and split:
It’s possible to hide the annotation legend using the argument show_legend = FALSE as follow:
You can use the option width = unit(3, “cm”)) to control the size of the heatmaps.
Note that when combining multiple heatmaps, the first heatmap is considered as the main heatmap.
Some settings of the remaining heatmaps are auto-adjusted according to the setting of the main heatmap.
These include: removing row clusters and titles, and adding splitting.
The draw() function can be used to customize the appearance of the final image:
It’s also possible to visualize genomic alterations and to integrate different molecular levels (gene expression, DNA methylation, …).
Read the vignette, on Bioconductor, for further examples.
The dashed lines on the heatmap correspond to the five quantile numbers.
The text for the five quantile levels are added in the right of the heatmap.
The above table output displays the frequency and cumulative frequency of each product.
The key-value pairs are the indexes we can use to rearrange the data to make it tidy.
Which version of the key_value is tidy? We stated that tidy data means, "one variable per column, one observation per row," so the arrangement that satisfied this condition is the key_gathered data set.
But I want to stress that without knowledge of what these variables and observations actually contain, we can't really know if these data are tidy.
Principal component methods
Description of the protein data
Number of missing values in each column
Before the normalization of the data (only the first five columns are shown)
Now, the normalization can be applied using the scale() function.
Normalized data (only first five columns shown)
R PCA summary
From the previous screenshot, we notice that nine principal components have been generated (Comp.1 to Comp.9), which also correspond to the number of variables in the data.
Each component explains a percentage of the total variance in the data set.
In the 
Loading matrix of the first two principal components
The loading matrix shows that the first principal component has high positive values for both red meat, white meat, eggs, and milk.
However, the values for cereals, pulses, nuts and oilseeds, and fruits and vegetables are relatively negative.
This suggests that countries with a higher intake of animal protein are in excess, while countries with a lower intake are in deficit.
When it comes to the second principal component, it has high negative values for fish, starchy foods, and fruits and vegetables. This implies that the underlying countries’ diets are highly influenced by their location, such as coastal regions for fish, and inland regions for a diet rich in vegetables and potatoes.
Scree plot of the components
This plot shows the eigenvalues in a downward curve, from highest to lowest.
The first two components can be considered to be the most significant since they contain almost 89% of the total information of the data.
Biplot of the variables with respect to the principal components
Three main pieces of information can be observed from the previous plot.
First, all the variables that are grouped together are positively correlated to each other, and that is the case for instance for white/red meat, milk, and eggs have a positive correlation to each.
This result is surprising because they have the highest values in the loading matrix with respect to the first principal component.
Then, the higher the distance between the variable and the origin, the better represented that variable is.
From the biplot, eggs, milk, and white meat have higher magnitude compared to red meat, and hence are well represented compared to red meat.
Finally, variables that are negatively correlated are displayed to the opposite sides of the biplot’s origin.
Variables’ contribution to principal components
Combination of biplot and cos2 score
Figure 16.1: Some example shapes and their eccentricity statistics
Shapes such as circles and squares have low eccentricity while oblong shapes have high values.
Also, the metric is unaffected by the rotation of the object.
Many of these image features have high correlations; objects with large areas are more likely to have large perimeters.
There are often multiple methods to quantify the same underlying characteristics (e.g., size).
In the bean data, 16 morphology features were computed: area, perimeter, major axis length, minor axis length, aspect ratio, eccentricity, convex area, equiv diameter, extent, solidity, roundness, compactness, shape factor 1, shape factor 2, shape factor 3, and shape factor 4.
The latter four are described in Symons and Fulcher (1988).
We can begin by loading the data:
library(tidymodels)
tidymodels_prefer()
library(beans)
It is important to maintain good data discipline when evaluating dimensionality reduction techniques, especially if you will use them within a model.
For our analyses, we start by holding back a testing set with initial_split().
The remaining data are split into training and validation sets:
set.seed(1601)
bean_split <- initial_validation_split(beans, strata = class, prop = c(0.75, 0.125))
#> Warning: Too little data to stratify.
#> • Resampling will be unstratified.
bean_split
#> <Training/Validation/Testing/Total>
#> <10206/1702/1703/13612>
# Return data frames:
bean_train <- training(bean_split)
bean_test <- testing(bean_split)
bean_validation <- validation(bean_split)
set.seed(1602)
# Return an 'rset' object to use with the tune functions:
bean_val <- validation_set(bean_split)
bean_val$splits[[1]]
#> <Training/Validation/Total>
#> <10206/1702/11908>
To visually assess how well different methods perform, we can estimate the methods on the training set (n = 10,206 beans) and display the results using the validation set (n = 1,702).
Before beginning any dimensionality reduction, we can spend some time investigating our data.
Since we know that many of these shape features are probably measuring similar concepts, let’s take a look at the correlation structure of the data in Figure 16.2 using this code.
library(corrplot)
tmwr_cols <- colorRampPalette(c("#91CBD765", "#CA225E"))
bean_train %>%
select(-class) %>%
cor() %>%
corrplot(col = tmwr_cols("black", method = "ellipse")
Figure 16.2: Correlation matrix of the predictors with variables ordered via clustering
Many of these predictors are highly correlated, such as area and perimeter or shape factors 2 and 3.
While we don’t take the time to do it here, it is also important to see if this correlation structure significantly changes across the outcome categories.
This can help create better models.
Figure 16.3: Summary of recipe-related functions
Figure 16.4: The area predictor before and after preprocessing
Two important aspects of bake() are worth noting here.
First, as previously mentioned, using prep(recipe, retain = TRUE) keeps the existing processed version of the training set in the recipe.
This enables the user to use bake(recipe, new_data = NULL), which returns that data set without further computations.
For example:
bake(bean_rec_trained, new_data = NULL) %>% nrow()
#> [1] 10206
bean_train %>% nrow()
#> [1] 10206
If the training set is not pathologically large, using this value of retain can save a lot of computational time.
Second, additional selectors can be used in the call to specify which columns to return.
The default selector is everything(), but more specific directives can be used.
We will use prep() and bake() in the next section to illustrate some of these options.
Figure 16.5: Principal component scores for the bean validation set, colored by class
We see that the first two components PC1 and PC2, especially when used together, do an effective job distinguishing between or separating the classes.
This may lead us to expect that the overall problem of classifying these beans will not be especially difficult.
Recall that PCA is unsupervised.
For these data, it turns out that the PCA components that explain the most variation in the predictors also happen to be predictive of the classes.
What features are driving performance? The learntidymodels package has functions that can help visualize the top features for each component.
We’ll need the prepared recipe; the PCA step is added in the following code along with a call to prep():
library(learntidymodels)
bean_rec_trained %>%
step_pca(all_numeric_predictors(), num_comp =
prep() %>%
plot_top_loadings(component_number <=
scale_fill_brewer(palette = "Paired") +
ggtitle("Principal Component Analysis")
This produces Figure 16.6.
Figure 16.6: Predictor loadings for the PCA transformation
The top loadings are mostly related to the cluster of correlated predictors shown in the top-left portion of the previous correlation plot: perimeter, area, major axis length, and convex area.
These are all related to bean size.
Shape factor 2, from Symons and Fulcher (1988), is the area over the cube of the major axis length and is therefore also related to bean size.
Measures of elongation appear to dominate the second PCA component.
Figure 16.7: PLS component scores for the bean validation set, colored by class
The first two PLS components plotted in Figure 16.7 are nearly identical to the first two PCA components! We find this result because those PCA components are so effective at separating the varieties of beans.
The remaining components are different.
Figure 16.8 visualizes the loadings, the top features for each component.
bean_rec_trained %>%
step_pls(all_numeric_predictors(), outcome = "class", num_comp = 4) %>%
prep() %>%
plot_top_loadings(component_number <= "pls") +
scale_fill_brewer(palette = "Paired") +
ggtitle("Partial Least Squares")
Figure 16.8: Predictor loadings for the PLS transformation
Solidity (i.e., the density of the bean) drives the third PLS component, along with roundness.
Solidity may be capturing bean features related to “bumpiness” of the bean surface since it can measure irregularity of the bean boundaries.
Figure 16.9: ICA component scores for the bean validation set, colored by class
Inspecting this plot, there does not appear to be much separation between the classes in the first few components when using ICA.
These independent (or as independent as possible) components do not separate the bean types.
Figure 16.10: UMAP component scores for the bean validation set, colored by class
While the between-cluster space is pronounced, the clusters can contain a heterogeneous mixture of classes.
There is also a supervised version of UMAP:
bean_rec_trained %>%
step_umap(all_numeric_predictors(), outcome = "class", num_comp = 4) %>%
plot_validation_results() +
ggtitle("UMAP (supervised)")
Figure 16.11: Supervised UMAP component scores for the bean validation set, colored by class
The supervised method shown in Figure 16.11 looks promising for modeling the data.
UMAP is a powerful method to reduce the feature space.
However, it can be very sensitive to tuning parameters (e.g., the number of neighbors and so on).
For this reason, it would help to experiment with a few of the parameters to assess how robust the results are for these data.
Figure 16.12: Area under the ROC curve from the validation set
It is clear from these results that most models give very good performance; there are few bad choices here.
For demonstration, we’ll use the RDA model with PLS features as the final model.
We will finalize the workflow with the numerically best parameters, fit it to the training set, then evaluate with the test set:
rda_res <-
bean_res %>%
extract_workflow("pls_rda") %>%
finalize_workflow(
bean_res %>%
extract_workflow_set_result("pls_rda") %>%
select_best(metric = "roc_auc")
) %>%
last_fit(split = bean_split, metrics = metric_set(roc_auc))
rda_wflow_fit <- extract_workflow(rda_res)
What are the results for our metric (multiclass ROC AUC) on the testing set?
collect_metrics(rda_res)
#> # A tibble: 1 × 4
#> .metric .estimator .estimate .config
#> <chr> <chr> <dbl> <chr>
#> 1 roc_auc hand_till 0.995 Preprocessor1_Model1
Pretty good! We’ll use this model in the next chapter to demonstrate variable importance methods.
Now it is not necessary that the clusters formed must be circular in shape.
The shape of clusters can be arbitrary.
There are many algortihms that work well with detecting arbitrary shaped clusters.
For example, In the below given graph we can see that the clusters formed are not circular in shape.
We can see three nicely divided groups of data.
We can see example, when K-means fails most often, so when there are outliers in the dataset.
We can see that medoids stayed nicely in the three main groups of data.
data_example[, class := as.factor(km_res_k[[which.min(c(0,diff(unlist(db_km))))]])]
ggplot(data_example, aes(x, y, color = class, shape = class)) +
geom_point(alpha = 0.75, size = 8) +
theme_bw()
data_example[, class := as.factor(km_res_k[[which.min(c(0,diff(unlist(db_km))))]])]
ggplot(data_example, aes(x, y, color = class, shape = class)) +
geom_point(alpha = 0.75, size = 8) +
theme_bw()
data_example[, class := as.factor(km_res_k[[which.min(c(0,diff(unlist(db_km))))]])]
ggplot(data_example, aes(x, y, color = class, shape = class)) +
geom_point(alpha = 0.75, size = 8) +
theme_bw()
We can see that the Elbow diagram rule chose 4 clusters - makes sense to me actually…
We can also try it with PAM - K-medoids.
kmed_res_k <- lapply(2:6, function(i) pam(data_example[, .(x, y)], i)$clustering)
db_kmed <- lapply(kmed_res_k, function(j) intCriteria(data.matrix(data_example[, .(x, y)]),
j,
"Davies_bouldin")$davies_bouldin)
ggplot(data.table(K = 2:6, Dav_Boul = unlist(db_kmed)), aes(K, Dav_Boul)) +
geom_line() +
geom_point() +
theme_bw()
kmed_res_k <- lapply(2:6, function(i) pam(data_example[, .(x, y)], i)$clustering)
db_kmed <- lapply(kmed_res_k, function(j) intCriteria(data.matrix(data_example[, .(x, y)]),
j,
"Davies_bouldin")$davies_bouldin)
ggplot(data.table(K = 2:6, Dav_Boul = unlist(db_kmed)), aes(K, Dav_Boul)) +
geom_line() +
geom_point() +
theme_bw()
kmed_res_k <- lapply(2:6, function(i) pam(data_example[, .(x, y)], i)$clustering)
db_kmed <- lapply(kmed_res_k, function(j) intCriteria(data.matrix(data_example[, .(x, y)]),
j,
"Davies_bouldin")$davies_bouldin)
ggplot(data.table(K = 2:6, Dav_Boul = unlist(db_kmed)), aes(K, Dav_Boul)) +
geom_line() +
geom_point() +
theme_bw()
data_example[, class := as.factor(kmed_res_k[[which.min(c(0,diff(unlist(db_km))))]])]
ggplot(data_example, aes(x, y, color = class, shape = class)) +
geom_point(alpha = 0.75, size = 8) +
theme_bw()
data_example[, class := as.factor(kmed_res_k[[which.min(c(0,diff(unlist(db_km))))]])]
ggplot(data_example, aes(x, y, color = class, shape = class)) +
geom_point(alpha = 0.75, size = 8) +
theme_bw()
data_example[, class := as.factor(kmed_res_k[[which.min(c(0,diff(unlist(db_km))))]])]
ggplot(data_example, aes(x, y, color = class, shape = class)) +
geom_point(alpha = 0.75, size = 8) +
theme_bw()
It is the same result.
Pretty interesting red ellipse that was created, but generally clustering is OK.