にほんごのれんしゅう

日本語として伝えるための訓練を兼ねたテクログ

みんなのRを読んだ

みんなのRを読んだ

Rを身につけなければいけないね

  • 流行りのRを身につけないと将来食っていけなくなるのではないのか、みたいな漫然とした不安感があった。
  • Rはデータを解析する上で非常に便利なソフトウェアである。Pythonのtheanoを用いたとしても、ワンラインでかけるわけでない。Rは他の言語では複数行にまたがるようなデータを解析するという視点において高い生産性を誇っている。

www.amazon.co.jp

RStudioを使うことを前提としたセットアップと基本的なシンタックスが書かれた1~14章

  • どうやって環境構築していくのか、丁寧に書かれている。バージョンの違いか、GITの設定では異なっていたが、GUIの設計などは日々進化してしまうので、特に気にしなかった。
  • 行列操作、基本的な行列の四則演算、ラベルのつけ方、文字列操作などが詳しく記されているが、MatlabOctaveの知識がほぼ同じように転用できるので、学習量的には軽かった。

Rらしいデータの見方を記してある15章から

  • 15章からがこの本の真価だと思う。
  • データの検定、回帰、クラスタリング、L1,L2正規化など実践で欲しい部分。よって15章からの理解した情報を"理系の文章"形式でメモで構成する。

理系の文章の構成の前の、個々の要素で知った部分に関して、まとめずに記していく

sample作成

  • 重複を許可して、1から10の間のサンプルを100個作る
 x <- sample(x = 1:10, size=100, replace=TRUE)  
  • xの平均値
 mean(x)
  • 1から100までのサンプルに関して、20個をNAに設定する
  • meanを取ってもNAが含まれた行列はNAを返してしまう
 y <- sample(x=1:100, size=100, replace=FALSE)
 y[sample(x=1:100, size=20, replace=FALSE)] <- NA
  • NAを取り除いて平均する
 mean(y , na.rm=TRUE)
  • 加重平均を取る
 grades <- c(10, 20, 30, 40)
 weights <- c(0.1, 0.2, 0.3, 0.4)
 weighted.mean(grades, weights)
 [1] 30
  • 分散を取る
 var(x)
 sd(x)
  • min, max, midianもある
  • いろいろ返してくれる便利関数としてsummary

    共分散

  • require(ggplot2)でeconomicsのデータを使うとする
  • 相関係数を出す関数cor $pceは個人支出、$psavertは個人貯蓄率
  cor(economics$pce, economics$psavert)
 [1] -0.9271222

 逆相関である。なんとなく納得。

  • マトリックスの中の相関を一気にとるにはこんな書き方をする
  • 全部の相関が観れるのでかなり楽
 cor(economics[,c(2,4:6)])
                 pce     psavert    uempmed    unemploy
 pce       1.0000000 -0.92712221  0.5145862  0.32441514
 psavert  -0.9271222  1.00000000 -0.3615301 -0.07641651
 uempmed   0.5145862 -0.36153012  1.0000000  0.78427918
 unemploy  0.3244151 -0.07641651  0.7842792  1.00000000
  • ロングフォーマットへ変更
 econCor <- cor(economics[,c(2,4:6)])
 econMelt <- melt(econCor, varnames=c("x","y"), value.name="Correlation")

見やすい。

 ggplot(econMelt, aes(x=x, y=y)) + geom_tile(aes(fill=Correlation)) +
 +     scale_fill_gradient2(low=muted("red"), mid="white", high="blue")

t検定

  • とってもわかりやすいt検定。
 なんのt検定をしたいかを第一引数に入れて、第三引数に検定する数値を入れれば95%信頼区間で出してくれる。
 !>t.test(tips$tip, alternative="two.sided", mu=2.5)
 One Sample t-test
 data:  tips$tip
 t = 5.6253, df = 243, p-value = 5.08e-08
 alternative hypothesis: true mean is not equal to 2.5
 95 percent confidence interval:
 2.823799 3.172758
 sample estimates:
 mean of x
 2.998279

チップの平均が2.5ドルかっていうデータであり、平均値は2.9998279ドル、95%信頼区間は3.17ドルから2.82ドルで3$あげればいいっぽい
- 男女のチップの分散はどうなの

 # アンサリブラッドリー検定で一発
 ansari.test(tip ~ sex, tips)
    Ansari-Bradley test
 data:  tip by sex
 AB = 5582.5, p-value = 0.376
 alternative hypothesis: true ratio of scales is not equal to 1
  • 2値のt検定 女性と男性でもらえるチップの分散が等しい時、有意差があるか
 t.test(tip ~ sex, data=tips, var.equal = TRUE)
 Two Sample t-test
 data:  tip by sex
 t = -1.3879, df = 242, p-value = 0.1665
 alternative hypothesis: true difference in means is not equal to 0
 95 percent confidence interval:
 -0.6197558  0.1074167
 sample estimates:
 mean in group Female   mean in group Male
            2.833448             3.089618

t検定によると無い。

分散分析

  • 他のグループとどの程度それが違うか
  • aov関数でサクッといける
 tipAnova <- aov(tip ~ day, tips)
 tipAnova$coefficients
(Intercept)      daySat      daySun     dayThur
 2.73473684  0.25836661  0.52039474  0.03671477

日曜だけ断トツで違うねってわかるのが嬉しいのかな。 より厳密に見ていこうとすると、可視化して90%信頼レベルで互いの曜日との違いを出すのがいい。

重回帰

  • 流行りの重回帰
  • これもかなり簡単にできてしまって、ええーって感じだった
  • "http://www.jaredlander.com/data/housing.csv"でデータが置いてあるので、ダウンロードして評価。ラベルが付いていないのでnames関数でつける必要がある。
 housing <- read.table("http://www.jaredlander.com/data/housing.csv", sep=",", header=TRUE, stringsAsFactors = FALSE)
 names(housing) <- c("Nh", "Class", "Units","Year","SqFt","Income","IncomSqFt", "Expence","Expence`erSqFt","NetIncome","Value","ValueperSqFt","Boro")
 #よくある、コンドミニアムの一フィートあたりの値段の関係のデータを用いる
 lm(formula = ValueperSqFt ~ Units + SqFt + Boro, data = housing)
  • これで一発重回帰
  • 正しいのかどうなっているのかよくわからないと思うので、グラフ化をお勧めする
 require(coefplot)
 要求されたパッケージ coefplot をロード中です
 coefplot(house1)

ロジスティック回帰

  • なかなかエグい世帯年収を説明するためのロジットを組む
  • "http://jaredlander.com/data/acs_ny.csv"にデータがあるので、ゲットして評価。ラベルが付いているのでname関数は必要無い
  • binary(1 or 0)にしなければならないので150000ドル以上稼いでいるか否かで判断
 acs <- read.table("http://jaredlander.com/data/acs_ny.csv",sep=",",header=TRUE, stringsAsFactors=FALSE)
 head(acs)
   Acres FamilyIncome  FamilyType NumBedrooms NumChildren NumPeople NumRooms        NumUnits
 1  1-10          150     Married           4           1         3        9 Single detached
 2  1-10          180 Female Head           3           2         4        6 Single detached
 3  1-10          280 Female Head           4           0         2        8 Single detached
 4  1-10          330 Female Head           2           1         2        4 Single detached
 5  1-10          330   Male Head           3           1         2        5 Single attached
 6  1-10          480   Male Head           0           3         4        1 Single detached
   NumVehicles NumWorkers  OwnRent   YearBuilt HouseCosts ElectricBill FoodStamp HeatingFuel
 1           1          0 Mortgage   1950-1959       1800           90        No         Gas
 2           2          0   Rented Before 1939        850           90        No         Oil
 3           3          1 Mortgage   2000-2004       2600          260        No         Oil
 4           1          0   Rented   1950-1959       1800          140        No         Oil
 5           1          0 Mortgage Before 1939        860          150        No         Gas
 6           0          0   Rented Before 1939        700          140        No         Gas
   Insurance       Language
 1      2500        English
 2         0        English
 3      6600 Other European
 4         0        English
 5       660        Spanish
 6         0        English
   acs$income <- with(acs, FamilyIncome >= 150000)
  • 以下のコマンドで一発ロジスティック回帰
 income1 <- glm(income ~ HouseCosts + NumWorkers + OwnRent + NumBedrooms + FamilyType, data=acs, family = binomial(link="logit"))
  • 結果
 summary(income1)
 Call:
 glm(formula = income ~ HouseCosts + NumWorkers + OwnRent + NumBedrooms +
    FamilyType, family = binomial(link = "logit"), data = acs)
 Deviance Residuals:
     Min       1Q   Median       3Q      Max  
 -2.8452  -0.6246  -0.4231  -0.1743   2.9503  
 Coefficients:
                       Estimate Std. Error z value Pr(>|z|)
 (Intercept)         -5.738e+00  1.185e-01 -48.421   <2e-16 ***
 HouseCosts           7.398e-04  1.724e-05  42.908   <2e-16 ***
 NumWorkers           5.611e-01  2.588e-02  21.684   <2e-16 ***
 OwnRentOutright      1.772e+00  2.075e-01   8.541   <2e-16 ***
 OwnRentRented       -8.886e-01  1.002e-01  -8.872   <2e-16 ***
 NumBedrooms          2.339e-01  1.683e-02  13.895   <2e-16 ***
 FamilyTypeMale Head  3.336e-01  1.472e-01   2.266   0.0235 *  
 FamilyTypeMarried    1.405e+00  8.704e-02  16.143   <2e-16 ***

 - できるだけみんなで働いて、結婚していて高い家に住めとしか言えない

モデル比較

  • 様々なモデルを立てて、残差二乗和が一番低いものを探すのにanova関数を用いる
 house2 <- lm(ValueperSqFt ~ Units * SqFt + Boro, data = housing)
 house3 <- lm(ValueperSqFt ~ Units + SqFt * Boro + Class, data = housing)
 house4 <- lm(ValueperSqFt ~ Units + SqFt * Boro + SqFt * Class, data = housing)
 house5 <- lm(ValueperSqFt ~ Boro + Class, data = housing)
 anova(house2, house3, house4, house5)
 Analysis of Variance Table
 Model 1: ValueperSqFt ~ Units * SqFt + Boro
 Model 2: ValueperSqFt ~ Units + SqFt * Boro + Class
 Model 3: ValueperSqFt ~ Units + SqFt * Boro + SqFt * Class
 Model 4: ValueperSqFt ~ Boro + Class
   Res.Df     RSS Df Sum of Sq       F    Pr(>F)
 1   2618 4884872
 2   2612 4619926  6    264946 25.1727 < 2.2e-16 ***
 3   2609 4576671  3     43255  8.2195  1.92e-05 ***
 4   2618 4901463 -9   -324792 20.5725 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  • house4が最も、残差二乗和が小さい
  • しかし、house4はオーバフィッティングを起こしている可能性がある
  • 赤池情報量基準AICなどで、モデルの複雑性に対してペナルティを加えてオーバフィッティングになっていないかを確認することができる
  • AICは負の値が大きいほど良い
 AIC(house2, house3, house4, house5)
        df      AIC
 house2  9 27239.94
 house3 15 27105.50
 house4 18 27086.80

クロスバリデーション

  • k個にデータを分割し、k-1を予想・学習に用いて、最後の一つを評価に用いる手法。
  • AICやAnovaよりモダンな方法でより良いとされている
  • bootパッケージに含まれている。cv.glmでクロスバリデーションを実行可能
 houseG2 <- glm(ValueperSqFt ~ Units * SqFt + Boro, data = housing)
 houseG3 <- glm(ValueperSqFt ~ Units + SqFt*Boro + Class, data=housing)
 houseG4 <- glm(ValueperSqFt ~ Units + SqFt*Boro + SqFt*Class,data=housing)
 houseG5 <- glm(ValueperSqFt ~ Boro + Class, data=housing)
 houseCV2 <- cv.glm(housing, houseG2, K=5)
 houseCV3 <- cv.glm(housing, houseG3, K=5)
 houseCV4 <- cv.glm(housing, houseG4, K=5)
 houseCV5 <- cv.glm(housing, houseG5, K=5)

bootstrap

  • bootstrapについて十分な理解ができていなかったので、再確認を行った。
  • ネット上のGoogleのランキングが高い資料では、何がいいたいのかよくわからなかったので、色々調べたところ下記の説明が一番しっくりきた。人によって理解のしやすさが異なると思うので自身の最適な資料を集めるといい。
  • https://oku.edu.mie-u.ac.jp/~okumura/stat/bootstrap.html
  • パラメトリックな分母に従うと仮定される場合、非常に強力。サンプルを取らなくても標準偏差などが求めることが可能
 require(plyr)
 baseball <- baseball[baseball$year >= 1990,]
 bat.avg <- function(data, indices=1:nrow(data), hits="h", at.bats="ab"){ sum(data[indices, hits], na.rm=TRUE)/sum(data[indices, at.bats], na.rm=TRUE)}
 avgBoot <- boot(data=baseball, statistics=bat.avg, R=1200, stype="i")
 avgBoot
 ORDINARY NONPARAMETRIC BOOTSTRAP
 Call:
 boot(data = baseball, statistic = bat.avg, R = 1200, stype = "i")
 Bootstrap Statistics :
      original        bias    std. error
 t1* 0.2729972 -1.547748e-05 0.000674947
 boot.ci(avgBoot, conf=.95, type="norm")
 BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
 Based on 1200 bootstrap replicates
 CALL :
 boot.ci(boot.out = avgBoot, conf = 0.95, type = "norm")
 Intervals :
 Level      Normal
 95%   ( 0.2717,  0.2743 )  
 Calculations and Intervals on Original Scale
 ggplot() + geom_histogram(aes(x = avgBoot$t)) + geom_vline(xintercept=avgBoot$t0 + c(-1,1)*2*sqrt(var(avgBoot$t)), linetype=2)

正規化と縮小(シュリンケージ)

ElasticNet

  • ラッソー回帰(L1正規化)とリッジ回帰(L2正規化)を動的にまぜたもの
  • ラッソー回帰で変数選択と次元削減を行う
  • リッジ回帰で係数を縮小する
  • λという変数により、ラッソー回帰とリッジ回帰のバランスを設定する
  • モデルを作るときにはmodel.matrix関数が便利。でも、多重共線性の問題があり、非常に面倒な問題とされている
  • 多重共線性をブラックボックスで解いてもらうのに、usefulパッケージのbuild.x, build.y関数が便利
acs <- read.table("http://jaredlander.com/data/acs_ny.csv", sep="," , header=TRUE, stringsAsFactors=FALSE)
acs$Income <- with(acs, FamilyIncome >= 150000)
acsX <- build.x(Income ~ NumBedrooms + NumChildren + NumPeople + NumRooms + NumUnits + NumVehicles + NumWorkers + OwnRent + YearBuilt + ElectricBill + FoodStamp + HeatingFuel + Insurance * Language -1, data=acs, contrasts=FALSE)
acsY <- build.y(Income ~ NumBedrooms + NumChildren + NumPeople + NumRooms + NumUnits + NumVehicles + NumWorkers + OwnRent +YearBuilt + ElectricBill + FoodStanmp + HeatingFuel + Insurance + Language - 1, data=acs)
require(glmnet)
acsCV1 <- cv.glmnet(x=acsX, y=acsY, family="binomial", nfold=5)
plot(acsCV1)
coef(acsCV1, s="lambda.1se")
50 x 1 sparse Matrix of class "dgCMatrix"
                                             1
(Intercept)                      -5.095309e+00
NumBedrooms                       5.620658e-02
NumChildren                       .
NumPeople                         .
NumRooms                          1.107383e-01
NumUnitsMobile home              -9.708518e-01
NumUnitsSingle attached           .
NumUnitsSingle detached           .
NumVehicles                       1.295513e-01
NumWorkers                        4.864577e-01
OwnRentMortgage                   .
OwnRentOutright                   3.184303e-01
OwnRentRented                    -2.009643e-01
YearBuilt15                       .
YearBuilt1940-1949               -4.362719e-02
YearBuilt1950-1959                .
YearBuilt1960-1969                .
YearBuilt1970-1979               -2.281293e-02
YearBuilt1980-1989                2.651037e-02
YearBuilt1990-1999                .
YearBuilt2000-2004                1.180907e-02
YearBuilt2005                     .
YearBuilt2006                     .
YearBuilt2007                     1.010983e-02
YearBuilt2008                     .
YearBuilt2009                     .
YearBuilt2010                     .
YearBuiltBefore 1939             -1.934594e-01
ElectricBill                      1.837836e-03
FoodStampNo                       7.321574e-01
FoodStampYes                      .
HeatingFuelCoal                  -3.260396e-01
HeatingFuelElectricity           -1.561999e-02
HeatingFuelGas                    .
HeatingFuelNone                   .
HeatingFuelOil                    .
HeatingFuelOther                  .
HeatingFuelSolar                  .
HeatingFuelWood                  -7.784783e-01
Insurance                         4.597842e-04
LanguageAsian Pacific             4.193924e-01
LanguageEnglish                  -1.986741e-02
LanguageOther                     .
LanguageOther European            9.601302e-02
LanguageSpanish                   .
Insurance:LanguageAsian Pacific   .
Insurance:LanguageEnglish         5.246845e-05
Insurance:LanguageOther           .
Insurance:LanguageOther European  .
Insurance:LanguageSpanish         .  

燃料に石炭や木を使うと年収が低いのは容易に想像できるが、数字として出てくるとよく判断できる。

Baysian縮小

  • 事前分布を仮定することで、標準誤差を縮小するもの。十分に設計されていないデータのモデル適合を助ける役割がある。
  • 逆に、複雑なモデルで過適合してしまう場合にはElasticNetなど次元圧縮系が有効
  • 関数名はarm::bayesglm

非線形モデル

最小二乗法

  • nls関数で説明変数を最小二乗法したいデータを用いればOK

    spline

  • 最小二乗法とN回微分の二乗をともに最小化する
  • 以下に例を乗せる
 smooth.spline(x=diamonds$carat, y=diamonds$price, df=N)

一般加法モデル(GAM)

  • 全ての変数にFx_n(X_n)のような平準化関数を噛ませて、和を取る方法
  • 書籍では、クレジットカードの支払い情報ととの因果を述べているが、フォーマットが処理しにくい

クラスタリング

  • k-means, k-medoids(pam), clusterなど
  • k-meansだけ示す
wine <- read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", header=TRUE, sep=",")
names(wine) <- c("cultivar","alcohol","Malic acid","Ash","Alcalinity of ash",  "Magnesium","Total phenols","Flavanoids","Nonflavanoid phenols","Proanthocyanins","Color intensity","Hue","OD280/OD315 of diluted wines","Proline" )
wineTrain <- wine[, which(name(wine) != "cultivar")]
wineK3 <- kmeans(x=wineTrain, centers=3)
plot(wineK3, data=wineTrain)
  • クラスタを幾つに割ったら適切なのか|下記のプロット図で0に近いほどいい
wineBest <- FitKMeans(wineTrain, max.clusters=20, nstart=25, seed=278613)
PlotHartigan(wineBest)
  • 外れ値に対して強くなるには、k-meansでなくk-medoidsなどを使う

みんなのRを読むにあたって

  • 一度、学習や復習を兼ねて通しで読んでみることをお勧めします
  • 辞書的にも使えそうだと思ったので、電子書籍でなくて実物として購入し、デスクの上に置いておくことをお勧めします
  • 業務に関係のない部分は(特に15章以降)飛ばしてもいいと思います
  • DeepLearningに使うニューラルネットワーク以外は大抵Rで完結できると感じました。労働生産性を高めるためにも積極的につかっていきたいものです。