10個(gè)令人相見恨晚的R語言包
大約3年前我開始使用R,起初進(jìn)展很慢,與我習(xí)慣的語言相比,語法更加直觀也比較簡(jiǎn)單,而且需要一段時(shí)間才能習(xí)慣于細(xì)微的差別。我還不清楚語言的力量與社區(qū)和各種包的密切關(guān)系。
和其他語言(比如Python和Java)相比,R可以更模糊和麻煩。好消息是,有大量的包可以在R基礎(chǔ)庫上提供簡(jiǎn)單和熟悉的界面。這篇文章是我喜歡和每天使用的10個(gè)包,并且我希望自己能早些知道他們。
1. sqldf
R語言學(xué)習(xí)曲線中最陡峭的一部分就是語法,我花了一段時(shí)間才習(xí)慣使用<-代替=。我聽到很多人問如何實(shí)現(xiàn)VLOOKUP?!?R 對(duì)于一般的數(shù)據(jù)粗加工任務(wù)非常有用,但需要一段時(shí)間才能掌握??梢哉J(rèn)為sqldf是我的R”輔助輪子”。
sqldf讓你在R數(shù)據(jù)框上執(zhí)行SQL查詢。來自SAS的人會(huì)發(fā)現(xiàn)它非常熟悉,任何具有基本SQL技能的人都可以輕松的使用它—sqldf使用SQLite語法。
- install.packages("sqldf")
- library(sqldf)
- sqldf("SELECT
- day
- , avg(temp) as avg_temp
- FROM beaver2
- GROUP BY
- day;")
- # day avg_temp
- #1 307 37.57931
- #2 308 37.71308
- #beavers1 和 beavers2 是R base 自帶的兩個(gè)數(shù)據(jù)集,記錄了兩種海貍的體溫序列
- beavers <- sqldf("select * from beaver1
- union all
- select * from beaver2;")
- #head(beavers)
- # day time temp activ
- #1 346 840 36.33 0
- #2 346 850 36.34 0
- #3 346 900 36.35 0
- #4 346 910 36.42 0
- #5 346 920 36.55 0
- #6 346 930 36.69 0
- movies <- data.frame(
- title=c("The Great Outdoors", "Caddyshack", "Fletch", "Days of Thunder", "Crazy Heart"),
- year=c(1988, 1980, 1985, 1990, 2009)
- )
- boxoffice <- data.frame(
- title=c("The Great Outdoors", "Caddyshack", "Fletch", "Days of Thunder","Top Gun"),
- revenue=c(43455230, 39846344, 59600000, 157920733, 353816701)
- )
- sqldf("SELECT
- m.*
- , b.revenue
- FROM
- movies m
- INNER JOIN
- boxoffice b
- ON m.title = b.title;")
- # title year revenue
- #1 The Great Outdoors 1988 43455230
- #2 Caddyshack 1980 39846344
- #3 Fletch 1985 59600000
- #4 Days of Thunder 1990 157920733
如果你喜歡sqldf,可以使用pandasql包來查詢pandas中的DataFrame,通過SQL。
2. forecast
我不經(jīng)常做時(shí)間序列分析,但是當(dāng)我做的時(shí)候forecast包是我的選擇。forecast對(duì)ARIMA,ARMA,AR,指數(shù)平滑等時(shí)間序列模型的預(yù)測(cè)簡(jiǎn)單的令人難以置信。
- install.packages("forecast")
- library(forecast)
- # mdeaths: 英國(guó)每月死于肺病的人數(shù)
- fit <- auto.arima(mdeaths)
- #定制你的置信區(qū)間
- forecast(fit, level=c(80, 95, 99), h=3)
- # Point Forecast Lo 80 Hi 80 Lo 95 Hi 95 Lo 99 Hi 99
- #Jan 1980 1822.863 1564.192 2081.534 1427.259 2218.467 1302.952 2342.774
- #Feb 1980 1923.190 1635.530 2210.851 1483.251 2363.130 1345.012 2501.368
- #Mar 1980 1789.153 1495.048 2083.258 1339.359 2238.947 1198.023 2380.283
- plot(forecast(fit), shadecols="oldstyle")
我最喜歡的特性是產(chǎn)生預(yù)測(cè)的時(shí)序圖。
3. plyr
當(dāng)我***次使用R時(shí),我用基本的控制運(yùn)算來操縱數(shù)據(jù)(for, if, while, etc.)。我很快知道這是一個(gè)業(yè)余的做法,并且有更好的方法去實(shí)現(xiàn)。
在R中,apply函數(shù)族是在對(duì)列表或者向量每個(gè)元素調(diào)用函數(shù)的***方法。雖然R基礎(chǔ)庫中有這些函數(shù),但它們的使用可能難以掌握。我發(fā)現(xiàn)plyr包 是一個(gè)對(duì)R基礎(chǔ)庫中諸如split,apply, combine的泛函的更好用的替代。
plyr 給予你一些函數(shù) (ddply, daply, dlply, adply, ldply)按照常見的藍(lán)圖:將數(shù)據(jù)結(jié)構(gòu)分組拆分,對(duì)每個(gè)組應(yīng)用一個(gè)函數(shù),將結(jié)果返回到數(shù)據(jù)結(jié)構(gòu)中。
ddply 拆分一個(gè)數(shù)據(jù)框(data frame)并且返回一個(gè)數(shù)據(jù)框 (所以是 dd)。 daply 拆分一個(gè)數(shù)據(jù)框并且返回一個(gè)數(shù)組(array) (所以是 da)。希望你明白這個(gè)想法。
譯者注:plyr包包含了12個(gè)命名與功能相關(guān)的函數(shù),均以..ply命名,***個(gè).表示輸入的數(shù)據(jù)類型(a數(shù)組 d數(shù)據(jù)框 l列表),第二個(gè).表示輸出的數(shù)據(jù)類型(_表示不輸出)
- install.packages("plyr")
- library(plyr)
- # 按照 Species 拆分?jǐn)?shù)據(jù)庫,匯總一下,然后轉(zhuǎn)換結(jié)果
- # 到數(shù)據(jù)框
- ddply(iris, .(Species), summarise,
- mean_petal_length=mean(Petal.Length)
- )
- # Species mean_petal_length
- #1 setosa 1.462
- #2 versicolor 4.260
- #3 virginica 5.552
- # 按照 Species 拆分?jǐn)?shù)據(jù)庫,匯總一下,然后轉(zhuǎn)換結(jié)果
- # 到數(shù)組
- unlist(daply(iris[,4:5], .(Species), colwise(mean)))
- # setosa.Petal.Width versicolor.Petal.Width virginica.Petal.Width
- # 0.246 1.326 2.026
4. stringr
我發(fā)現(xiàn)R基礎(chǔ)庫的字符串功能使用起來非常困難和麻煩。Hadley Wickham編寫的另一個(gè)包, stringr,提供了一些非常需要的字符串運(yùn)算符。很多函數(shù)使用那些做基礎(chǔ)分析時(shí)不常用的數(shù)據(jù)結(jié)構(gòu)。
stringr 非常易于使用。幾乎所有的(和所有的重要功能)都以”str”為前綴,所以很容易記住。
- install.packages("stringr")
- library(stringr)
- names(iris)
- #[1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species"
- names(iris) <- str_replace_all(names(iris), "[.]", "_")
- names(iris)
- #[1] "Sepal_Length" "Sepal_Width" "Petal_Length" "Petal_Width" "Species"
- s <- c("Go to Heaven for the climate, Hell for the company.")
- str_extract_all(s, "[H][a-z]+ ")
- #[[1]]
- #[1] "Heaven " "Hell "
5. 數(shù)據(jù)庫驅(qū)動(dòng)的包
- install.packages("RPostgreSQL")
- install.packages("RMySQL")
- install.packages("RMongo")
- install.packages("RODBC")
- install.packages("RSQLite")
每個(gè)人(包括我自己)開始的時(shí)候都會(huì)這樣做。你剛在***的SQL編輯器中寫了一個(gè)很棒的查詢。一切都是***的—列名都是snake case(譯者注:表示單詞之間用下劃線連接。單詞要么全部大寫,要么全部小寫。),日期有正確的數(shù)據(jù)類型,***調(diào)試出了"must appear in the GROUP BY clause or be used in an aggregate function"的問題。你現(xiàn)在準(zhǔn)備在R中進(jìn)行一些分析,因此你可以在SQL編輯器中運(yùn)行查詢,將結(jié)果復(fù)制到csv(或者……xlsx)并讀入R,你并不需要這樣做!
R對(duì)于幾乎每一個(gè)可以想到的數(shù)據(jù)庫都有好的驅(qū)動(dòng)。當(dāng)你在偶爾使用不具有獨(dú)立驅(qū)動(dòng)程序的數(shù)據(jù)庫(SQL Server)時(shí),你可以隨時(shí)使用RODBC。
- library(RPostgreSQL)
- drv <- dbDriver("PostgreSQL")
- db <- dbConnect(drv, dbname="ncaa",
- user="YOUR USER NAME", password="YOUR PASSWORD")
- q <- "SELECT
- *
- FROM
- game_scores;"
- data <- dbGetQuery(db, q)
- head(data)
- #id school game_date spread school_score opponent opp_score was_home
- #1 45111 Boston College 1985-11-16 6.0 21 Syracuse 41 False
- #2 45112 Boston College 1985-11-02 13.5 12 Penn State 16 False
- #3 45113 Boston College 1985-10-26 -11.0 17 Cincinnati 24 False
- #4 45114 Boston College 1985-10-12 -2.0 14 Army 45 False
- #5 45115 Boston College 1985-09-28 5.0 10 Miami 45 True
- #6 45116 Boston College 1985-09-21 6.5 29 Pittsburgh 22 False
- nrow(data)
- #[1] 30932
- ncol(data)
- #[1] 8
下次你完成了***的查詢后,只需要粘貼到R里面,即可使用 RPostgreSQL, RMySQL,RMongo, SQLite, 或 RODBC執(zhí)行。不僅可以避免生成數(shù)以百計(jì)的CSV文件,在R中運(yùn)行查詢還可以節(jié)省I/O和轉(zhuǎn)換數(shù)據(jù)類型的時(shí)間。日期,時(shí)間等會(huì)自動(dòng)設(shè)置為R中的等價(jià)表示。它還使你的R腳本可重復(fù),因此你或你團(tuán)隊(duì)中的其他人可以輕松獲得相同的結(jié)果。
6. lubridate
在R中處理日期我從來沒有幸運(yùn)過。我從來沒有完全掌握用POSIXs和R內(nèi)建日期類型合作的方法。請(qǐng)用 lubridate。
lubridate 是那些似乎完全按照你期望的神包之一。這些函數(shù)都有易懂的名字如 year,month, ymd, 和 ymd_hms。對(duì)于熟悉javascript的人來說,它類似于Moment.js 。
- install.packages("lubridate")
- library(lubridate)
- year("2012-12-12")
- #[1] 2012
- day("2012-12-12")
- #[1] 12
- ymd("2012-12-12")
- #1 按照 %Y-%m-%d 的方式解析
- #[1] "2012-12-12 UTC"
這是我在一篇文章發(fā)現(xiàn)的非常方便的參考卡片。它涵蓋了你處理日期時(shí)可能想要做的一切事情。 我還發(fā)現(xiàn)了這個(gè)日期速查表也可以作為一個(gè)方便的參考。
7. ggplot2
另一個(gè)Hadley Wickham的包,也許是他最知名的一個(gè)。 ggplot2在每個(gè)人“喜愛的R包”的列表中排名很高。 它很容易使用,并且產(chǎn)生一些很好看的圖像。 這是介紹你的工作的好方法,有很多資源可以幫助你開始使用。
- ggplot2:數(shù)據(jù)分析與圖形藝術(shù) by Hadley Wickham (Amazon)
- 從Excel到ggplot的羅塞塔石碑 (Yaksis Blog)
- Hadley Wickham在Google的ggplot2演講 (youtube)
- R數(shù)據(jù)可視化手冊(cè) by Winston Chang (Amazon)
8. qcc
- install.packages("qcc")
- library(qcc)
- # 均值為10的序列,加上白噪聲
- x <- rep(10, 100) + rnorm(100)
- # 測(cè)試序列,均值為11
- new.x <- rep(11, 15) + rnorm(15)
- # qcc 會(huì)標(biāo)記出新的點(diǎn)
- qcc(x, newdata=new.x, type="xbar.one")
qcc 是用于 統(tǒng)計(jì)質(zhì)量控制的庫。 早在上世紀(jì)五十年代,現(xiàn)已不復(fù)存在的西方電氣公司正在尋找一種更好檢測(cè)電話線和電線問題的方法。他們提出了一系列 規(guī)則 來幫助識(shí)別有問題的線。規(guī)則觀察一系列數(shù)據(jù)點(diǎn)的歷史平均值,并且基于標(biāo)準(zhǔn)差的偏差,該規(guī)則有助于判斷一組新的點(diǎn)是否經(jīng)歷均值漂移。
典型的例子是監(jiān)控生產(chǎn) 螺母的機(jī)器。假設(shè)機(jī)器應(yīng)該生產(chǎn)2.5英寸長(zhǎng)的螺母。我們測(cè)量一系列的螺母: 2.48, 2.47, 2.51, 2.52, 2.54, 2.42, 2.52, 2.58, 2.51。機(jī)器出故障了嗎?很難說,但上述規(guī)則可以幫助描述。
雖然你可能不會(huì)監(jiān)控電話線,qcc可幫你監(jiān)控你網(wǎng)站的交易量,數(shù)據(jù)庫的訪問者或者登錄名,以及其他許多流程。
9. reshape2
我經(jīng)常發(fā)現(xiàn),任何分析中最難的部分是把數(shù)據(jù)轉(zhuǎn)化成正確的格式。 reshape2 正是Hadley Wickham的另一個(gè)軟件包,專門用于 “寬”數(shù)據(jù)表 和“窄”數(shù)據(jù)表 的轉(zhuǎn)換。我一般會(huì)和ggplot2 及 plyr一起使用它。
- install.packages("reshape2")
- library(reshape2)
- # 為每一行生成唯一的ID; 這樣我們可以稍后轉(zhuǎn)回到寬格式
- iris$id <- 1:nrow(iris)
- iris.lng <- melt(iris, id=c("id", "Species"))
- head(iris.lng)
- # id Species variable value
- #1 1 setosa Sepal.Length 5.1
- #2 2 setosa Sepal.Length 4.9
- #3 3 setosa Sepal.Length 4.7
- #4 4 setosa Sepal.Length 4.6
- #5 5 setosa Sepal.Length 5.0
- #6 6 setosa Sepal.Length 5.4
- iris.wide <- dcast(iris.lng, id + Species ~ variable)
- head(iris.wide)
- # id Species Sepal.Length Sepal.Width Petal.Length Petal.Width
- #1 1 setosa 5.1 3.5 1.4 0.2
- #2 2 setosa 4.9 3.0 1.4 0.2
- #3 3 setosa 4.7 3.2 1.3 0.2
- #4 4 setosa 4.6 3.1 1.5 0.2
- #5 5 setosa 5.0 3.6 1.4 0.2
- #6 6 setosa 5.4 3.9 1.7 0.4
- library(ggplot2)
- # 為數(shù)據(jù)集中每個(gè)數(shù)值列繪制直方圖
- p <- ggplot(aes(x=value, fill=Species), data=iris.lng)
- p + geom_histogram() +
- facet_wrap(~variable, scales="free")
這是一個(gè)快速查看數(shù)據(jù)集并且獲得轉(zhuǎn)接的方法。你可以使用 melt 函數(shù)將寬數(shù)據(jù)轉(zhuǎn)換為窄數(shù)據(jù), 使用 dcast 將窄數(shù)據(jù)轉(zhuǎn)換為寬數(shù)據(jù)。
10. randomForest
如果這個(gè)列表不包括至少一個(gè)能你的朋友震驚的機(jī)器學(xué)習(xí)包就不會(huì)完整。隨機(jī)森林 是一個(gè)很好的算法。它很容易使用,可以進(jìn)行監(jiān)督學(xué)習(xí)或者無監(jiān)督學(xué)習(xí),它可以與許多不同類型的數(shù)據(jù)集一起使用,但最重要的是它的高效率!這是它在R中的使用方法。
- install.packages("randomForest")
- library(randomForest)
- # 下載泰坦尼克號(hào)幸存者數(shù)據(jù)集
- data <- read.table("http://math.ucdenver.edu/RTutorial/titanic.txt", h=T, sep="\t")
- # 將Survived列轉(zhuǎn)為yes/no因子
- data$Survived <- as.factor(ifelse(data$Survived==1, "yes", "no"))
- # 拆分為訓(xùn)練集和測(cè)試集
- idx <- runif(nrow(data)) <= .75
- data.train <- data[idx,]
- data.test <- data[-idx,]
- # 訓(xùn)練一個(gè)隨機(jī)森林
- rf <- randomForest(Survived ~ PClass + Age + Sex,
- data=data.train, importance=TRUE, na.action=na.omit)
- # 模型中每個(gè)變量的重要程度
- imp <- importance(rf)
- o <- order(imp[,3], decreasing=T)
- imp[o,]
- # no yes MeanDecreaseAccuracy MeanDecreaseGini
- #Sex 51.49855 53.30255 55.13458 63.46861
- #PClass 25.48715 24.12522 28.43298 22.31789
- #Age 20.08571 14.07954 24.64607 19.57423
- # 混淆矩陣 [[真反例, 假正例], [假反例, 正正例]]
- table(data.test$Survived, predict(rf, data.test), dnn=list("actual", "predicted"))
- # predicted
- #actual no yes
- # no 427 16
- # yes 117 195