Rの細かいTipsまとめ(小さいTipsの寄せ集め)


独立した記事にはならないが、それぞれ便利かつ重要な小さなRのTipsを紹介。

Contents

stratified sampling(層化抽出法)

ライブラリsamplingを使う

strata(data, stratanames=c('層化に使うカラム1', '層化に使うカラム2'), size=c(カラム1の抽出率, カラム2の抽出率))

method
'srswor': 非復元ランダムサンプリング(デフォルト)
'srswr': 復元ランダムサンプリング
'poisson': ポアソンサンプリング
'systematic': 系統抽出

Rで文字列をコマンドとして実行する(eval)

Rの繰り返し処理などでコマンドの一部のみを書き換えて同様の処理を行いたい、たとえば対象のデータフレームのみ変えて同じ処理を実行したい場合がある。

そういう場合には動的な部分(データフレームオブジェクト名)を変数として含むコマンドの文字列を生成し、その文字列をコマンドとして実行することになる。

言葉にするとややこしいが、他の言語でもあるeval処理である。

コマンドが1行の場合

最も単純なケースである。

基本形

ファイル名に動的部分を含むテキストファイル(activity_type1.txt)を読み込み、名前に動的部分を含むテーブル(t_activity_type1)に格納する例

.f <- 'activity_type1'
eval(parse(text=paste0("t_", .f, " <- fread('/path/to/datadir/", .f, ".txt', encoding='UTF-8', sep='\t', na.strings='')")))

.fが動的部分で、ファイル名である。

そして動的部分を

.f <- 'activity_type2'

などと変えても同様に実行できる。

eval(parse(text=paste0("固定文字列", 動的部分, "固定文字列")))

が基本になる。paste0()関数がパーツを結合して文字列を生成する関数で、それ以外はおまじないだと思っておけばいい。コマンドの固定文字列部分にクォーテーションマーク'を含む場合、文字列全体はダブルクォーテーションマーク"で囲む。

コマンド部分のみ文字列変数として切り出してソースを見やすくする

この記法だとソースコードの中で動的部分が見えにくく、コマンドの一部を書き換えたり繰り返し利用したりするのに向かない(特にカッコの数がわからなくなる)ので、コマンド部分を変数として取り出したのが以下になる。

.f <- 'activity_type1'
.str_cmd <- paste0("t_", .f, " <- fread('/path/to/datadir/", .f, ".csv', encoding='UTF-8', sep='\t', na.strings='')")
eval(parse(text=.str_cmd))

.str_cmdがコマンドの文字列になる。こうすると

.f <- 'activity_type1'
.str_cmd <- paste0("t_", .f, " <- fread('/path/to/datadir/", .f, ".csv', encoding='UTF-8', sep='\t', na.strings='')")
eval(parse(text=.str_cmd))
.str_cmd <- paste0("t_", .f, " %>% mutate_at(ends_with('_date'), as.Date)")
eval(parse(text=.str_cmd))

のように動的コマンドを複数生成して実行しやすくなる。

複数行のコマンドを実行する場合

上記の方法だと繰り返しeval(parse(text=.str_cmd))が必要になってしまう。動的に実行するコマンドが複数行の場合でも毎回必要になり、あまりイケてない。

そこで一時ファイルを生成して、一時ファイルに対して処理を行う。

.f <- 'activity_type1'
.str_cmd <- paste0(".tmp <- fread('/path/to/datadir/", .f, ".csv', encoding='UTF-8', sep='\t', na.strings='')")
eval(parse(text=.str_cmd)) # .tmpという一時ファイルに代入
# .tmpに対する処理
.tmp %>% 
  mutate_at(vars(ends_with('_date')), as.Date) %>% 
  mutate_at(vars(ends_with('_type')), as.factor) %>% 
  mutate_if(is.logical, as.integer) -> .tmp
# .tmpの内容を個別の永続テーブルに入れる
.str_cmd <- paste0(".tmp -> t_", .f); eval(parse(text=.str_cmd))

.tmpに対する処理の部分が動的ではない普通のコマンドとなり、ソースコードの可読性も上がり、再利用やメンテナンスもしやすくなった。

シンプルなコマンド1行であれば動的にする必要もあまりないが、複数行のコマンドを繰り返し実行したい場合にはこの方法が力を発揮する。

forループの中で実行する

以上の例を複数のファイルに対して繰り返し実行する。
動的な部分をベクトルに代入し、forループの中に先の処理を入れる。

# 動的な部分の定義
filenames <- c('activity_type1', 'activity_type7', 'activity_click', 'event_unsubscribe', 'event_old')
# 繰り返し
for (.f in filenames) {
  # いったんdata.tableオブジェクトとして取り込む
  .str_cmd <- paste0("t_", .f, " <- fread('/path/to/datadir/", .f, ".csv', encoding='UTF-8', sep='\t', na.strings='')")
  eval(parse(text=.str_cmd))
  # 作業用オブジェクトを生成
  eval(parse(text=paste(".tmp <- t_", .f, sep="")))
  .tmp %>% 
    mutate_at(vars(ends_with('_date')), as.Date) %>% 
    mutate_at(vars(ends_with('_type')), as.factor) %>% 
    mutate_if(is.logical, as.integer) -> .tmp
  .tmp[mail_type %like% '購入お礼', c('mail_type', 'mail_type2', 'mail_type3', 'channel'):=data.table('SAF', '購入お礼', NA, NA)]
  # 作業用オブジェクトから元のオブジェクトに戻す
  eval(parse(text=paste0(".tmp -> t_", .f)))
}

ソースコードがかなり洗練される。

# 条件に合致したオブジェクトに対して一括処理
for (.i in ls(pattern='t1_.*_r_')) {
  .str_cmd <- paste('class(', .i, ') <- "matrix"', sep='')
  eval(parse(text=.str_cmd))
}

Rで郵便番号マスタを作成する(データフレーム、data.table両対応)

分析用のデータの中には地域名がほしい形(都道府県名、市区町村名で区切られた形)で含まれておらず、一方で郵便番号が含まれていることもある。

その場合郵便番号から地域名を取得することになるが、その都度そのコードを書くのは面倒。
しかも郵便番号が更新されることもあるので、常に最新版を持ってくるようにしたい。

そこでコピペでそのまま郵便番号マスタを生成できるコードを用意した。
日本郵便の郵便番号データから最新版のCSVを取得してきて、郵便番号マスタを生成するものである。
これがあればRの中で常に最新版の郵便番号マスタを生成できる。
分析用データと結合して地域名を付けることも簡単である。

日本郵便の郵便番号データから最新版のCSVを取得する。
CSVファイルに含まれるカラムと形式の詳細は以下参照。

http://www.post.japanpost.jp/zipcode/dl/readme.html

必要なカラムのみ取得すればいい。
ここでは郵便番号(2列目)と市区町村コード(1列目)、都道府県名(7列目)、市区町村名(8列目)を取得する。なお郵便番号の形式はハイフンを含まない7ケタの数字列。

データフレームで読み込む場合

require(curl)
temp <- tempfile()
download.file("http://www.post.japanpost.jp/zipcode/dl/kogaki/zip/ken_all.zip", temp)
files <- unzip(temp)
df.zipcode <- read.csv(files[1], na.strings="", stringsAsFactors = FALSE, header=F, skip=0, fileEncoding = "shift-jis", colClasses = c("character", "NULL", "character", "NULL", "NULL", "NULL", "character", "character", "NULL", "NULL", "NULL", "NULL", "NULL", "NULL", "NULL"))
unlink(files)
unlink(temp)
rm(files, temp)
colnames(df.zipcode) <- c("citycode","zipcode","prefecture","city")

ちなみにdata.tableを使う場合はこんな感じ。
CSVファイルのエンコーディングがShift-JISでfread()が対応していないため、
一度read.csv()でデータフレームとして読み込んでからdata.tableに変換する。

data.tableで読み込む場合

require(curl)
require(data.table)
require(dplyr)
temp <- tempfile()
download.file("http://www.post.japanpost.jp/zipcode/dl/kogaki/zip/ken_all.zip", temp)
files <- unzip(temp)
read.csv(files[1], na.strings="", stringsAsFactors = FALSE, header=F, skip=0, fileEncoding = "shift-jis", colClasses = c("character", "NULL", "character", "NULL", "NULL", "NULL", "character", "character", "NULL", "NULL", "NULL", "NULL", "NULL", "NULL", "NULL")) %>% as.data.table -> dt.zipcode
unlink(files)
unlink(temp)
rm(files, temp)
setnames(dt.zipcode, "V1", "citycode")
setnames(dt.zipcode, "V3", "zipcode")
setnames(dt.zipcode, "V7", "prefecture")
setnames(dt.zipcode, "V8", "city")

こんなふうにすれば郵便番号を含むデータから都道府県、市区町村名も紐付けできる。
ちなみに元データの郵便番号の形式がハイフンを含む999-9999形式の場合

require(stringr)
dt.data<- dt.data %>% 
mutate(zipcode = str_replace(zipcode, "([0-9]{3})-([0-9]{4})", "\\1\\2")) %>% 
left_join(dt.zipcode)

モデルの選択

anova(単純なモデル, 複雑なモデル)

p-valueを見る

圧縮ファイルの読み込み

fread("gzip -dc df.txt")

パッケージを保持してアップグレード

https://www.r-bloggers.com/how-to-upgrade-r-without-losing-your-packages/

  1. Before you upgrade, build a temp file with all of your old packages.
tmp <- installed.packages()
installedpkgs <- as.vector(tmp[is.na(tmp[,"Priority"]), 1])
save(installedpkgs, file="installed_old.rda")
  1. Install the new version of R and let it do it’s thing.

  2. Once you’ve got the new version up and running, reload the saved packages and re-install them from CRAN.

tmp <- installed.packages()
installedpkgs.new <- as.vector(tmp[is.na(tmp[,"Priority"]), 1])
missing <- setdiff(installedpkgs, installedpkgs.new)
install.packages(missing)
update.packages()

If you had any packages from BioConductor, you will need to reload those as well.

chooseBioCmirror()
biocLite() 
load("installed_old.rda")
tmp <- installed.packages()
installedpkgs.new <- as.vector(tmp[is.na(tmp[,"Priority"]), 1])
missing <- setdiff(installedpkgs, installedpkgs.new)
for (i in 1:length(missing)) biocLite(missing[i])

多変量の可視化

表形式で―ftable()

度数で

ftable(segment ~ param1 + param2, data = data2)

              segment    1    2    3
param1 param2                       
10   10               4602  563  823
     20               5640  926  602
     30               9734  702  445
100  10               2300  256  288
     20               5435  914  375
     30               8845  545  403

比率で

round(prop.table(ftable(segment ~ param1 + param2, data = data2), margin=1) *100, 2)
  • margin=1で横方向の合計=1、margin=1で縦方向の合計=1とする。
  • useNA="no"でNAは集計対象から除外、useNA="always"でNAも集計対象に含める。

チャートで―coplot()boxplot()ggplot()library(ggplot2)

coplot()

coplot(ratio ~  click | ordered(tuningParam1) * ordered(tuningParam2), data=data1, panel = panel.smooth)

boxplot()

boxplot(ratio ~ x1 + x2, data1, outline = FALSE)
abline(h=1)

バイオリンプロット

ggplot(heightweight,aes(x=sex,y=heightIn)) + geom_violin() + geom_boxplot(width=.1,fill="black",outlier.colour=NA) + stat_summary(fun.y=median,geom="point",fill="white",shape=21,size=2.5)

高速化

  • apply()は内部でforの処理が入っている→
  • forが遅いのは、都度オブジェクトを拡張するから。for+rbind()などは最悪。
  • 最初にサイズを決めておくと速くなる。apply()を使うよりも速くなることもある。
  • ベクトル処理は高速
library(Rmpi)
library(snow)

# 1. apply (faster)
mat.matched <- mat.data==apply(mat.data, 1, max)

# 2. for (slower)
mat.matched <- matrix(nrow=nrow(mat.data), ncol=ncol(mat.data))
for(i in 1:nrow(mat.data)){ mat.matched[i,] <- mat.data[i,] == max(mat.data[i,]) }

# 3. parApply
cl <- makeCluster(4, type="SOCK")
clusterExport(cl, "mat.data")
mat.matched <- parApply(cl, mat.data, 1, max)
stopCluster(cl)

http://d.hatena.ne.jp/a_bicky/20120425/1335312593

blasの変更

yum install openblas
alternatives --install /usr/lib64/libblas.so.3 libblas.so.3 /usr/lib64/libopenblas.so.0 30
alternatives --install /usr/lib64/libblas.so.3 libblas.so.3 /usr/lib64/libblas.3.2.1 10
update-alternatives --config libblas.so.3

MPIのインストール

yum install -y openmpi-devel

事前にコマンドラインで環境変数をしてから

export LD_LIBRARY_PATH=/usr/lib64/openmpi/lib:$LD_LIBRARY_PATH
R

コンパイル時に3個の引数を渡す

install.packages("Rmpi", configure.args=c("--with-Rmpi-libpath=/usr/lib64/openmpi/lib", "--with-Rmpi-include=/usr/include/openmpi-x86_64", "--with-Rmpi-type=OPENMPI"))
library("Rmpi")
install.packages("snow")
library("snow")

参考
http://www.cybaea.net/Blogs/R-tips-Installing-Rmpi-on-Fedora-Linux.html

使い方
http://www.sfu.ca/~sblay/R/snow.html

Revolution R

  • Windowsではかなり速くなるがCentOSでは微妙
  • RevoMathはBLASにIntelのMKLを使って高速化(RROのインストール後にRevoMathをインストール)
  • 自動でマルチスレッド対応(HTは無視、純粋にコア数を見る)する。
  • 手動でやる場合は
library(RevoUtilsMath)
setMKLthreads(3)
  • getMKLthreads()でコア数を取得

RROではレポジトリが結構前の日付で固定されている。
最新版のパッケージをインストールしたい場合は

options(repos = c(CRAN = "https://mran.revolutionanalytics.com/snapshot/YYYY-MM-DD"))

を実行してからinstall.packages()を実行する。

caretが使っている

  • doMC (mac / unix)
  • doParallel (windows)

メモリ節約

sparse matix

要素のほとんどがゼロの行列→sparse matrix
普通のmatrixやデータフレーム形式で扱うとすべてのセルに対してメモリが割り当てられる
つまり行数×列数だけのメモリが消費される

ログデータを集計したもの(1000万UID×10万アイテム)などではメモリが足りないことに
行数と列数は膨大な割りに、実は値(の入るところ)が少ない→sparse matrix(疎行列)

もっと効率的にメモリを扱えないものか…値の入るセルだけメモリを使えないものか…
→これを実現したのが{Matrix}パッケージのsparse matrix形式(dgCMatrixクラスなど)

主成分分析

-(prcomp$xが主成分得点行列)
prcomp$sdevが固有値
prcomp$rotationが固有ベクトル

biplot(prcomp)でプロット

主成分得点×固有ベクトル=相関行列(scale=Tで実行した場合)

主成分得点を求める際に、主成分得点の不偏分散が固有値に一致するように計算される(prcomp()関数)か、主成分得点の標本分散が固有値に一致するように計算される(princomp()関数)かという違いがある

http://tmats.net/?p=2785

集合の便利な使い方

campaignIds.sampled <- intersect(campaignIds.sampled1, campaignIds.sampled2)

エラー処理

try{
    ...
}catch{
    ...
}

はないけど、try()関数で

try(for (i in unique(creatives.sampled$campaignId)) {
    for (j in adsize.sampled) {
        assign(paste0("result", i, "_", j, "_20161010"), myFunc(i, j, "2016-10-10"))
    }
})

正規表現にマッチしたオブジェクトのみ一括で削除

rm(list=ls(pattern="^result.*"))

ドットで始まる隠しオブジェクトを含むすべてのオブジェクトを削除

rm(list=ls(all=T))

ローカルの特定の形式に当てはまるファイル名のイメージを一括ロード

for(f in system("ls prcomp*.RData", TRUE)) load(f)

.Optionsでいろいろ

デフォルトでstringsAsFactors=FALSEにする。options()で設定しておくほうがいい

NAに対する処理

na.rm

  • NAを除外
  • 集計関数(max(), min(), range(), mean(), median(), quantile(), var())で使用。NAが入っているとデフォルトではNAを返すが、na.rm=TNAを除外してカウントする
  • cov(), cor()ではna.rmを指定できないので
cov(x[!is.na(x)], y[!is.na(y)])

とする。

na.omit

  • NAを含む行を除外
  • データフレームで使う
x <- data.frame(a=1:3, b=c(2,NA,4), c=rep(3,3))
na.omit(x)

tapply()よりaggregate()が使いやすい

aggregate(nExpected ~ category + type + audience + period, data = data1, FUN = sum)

aggregate(
    cbind(imp, click, cv) ~ segmentId + creativeId, 
    data = creatives.sampled, 
    FUN = sum
)

このようにformulaを使った集計ができるし、アウトプットもデータフレームになる。tapply()だと多次元配列になり扱いにくい。

tapply()の戻り値はarray

save()write.table()の引数の与え方の違い

似た処理なのに同じ引数の与え方で異なる挙動

paste()関数の活用法

eval(parse(text=paste(... , sep="")))
  • 代入→assign+pasteではリストの要素への代入ができない。
  • 関数の実行→do.call(fun, list(...))は実装がやや難しい。
  • eval+parse+pasteは代入も関数の実行もできるし、実装が簡単。
  • get(文字列)などという関数もあるが不要

ちょっと複雑な使い方

if (eval(parse(text=paste("is.na(", label, "[[1]][[index]])", sep=""))) == TRUE) { ... }
tbl <- eval(parse(text=paste(label, "[[1]][[", index, "]]$frequency", sep="")))
eval(parse(text=paste("result <- ", label, "[[1]]", sep="")))

ベクトルのノルム

# マンハッタンノルム
norm1 <- drop(crossprod(abs(m), rep(1, nrow(m))))
# ユークリッドノルム
norm2 <- drop(crossprod(m^2, rep(1, nrow(m))))

自作関数の検証に―オブジェクトが同一かどうかを検証する関数

  • identical(x1, x2)
  • all.equal(x1, x2)

Rのサンプルデータのわかりやすい解説

R関連記事