みんなのRを読んだ
みんなのRを読んだ
Rを身につけなければいけないね
- 流行りのRを身につけないと将来食っていけなくなるのではないのか、みたいな漫然とした不安感があった。
- Rはデータを解析する上で非常に便利なソフトウェアである。Pythonのtheanoを用いたとしても、ワンラインでかけるわけでない。Rは他の言語では複数行にまたがるようなデータを解析するという視点において高い生産性を誇っている。
RStudioを使うことを前提としたセットアップと基本的なシンタックスが書かれた1~14章
- どうやって環境構築していくのか、丁寧に書かれている。バージョンの違いか、GITの設定では異なっていたが、GUIの設計などは日々進化してしまうので、特に気にしなかった。
- 行列操作、基本的な行列の四則演算、ラベルのつけ方、文字列操作などが詳しく記されているが、MatlabやOctaveの知識がほぼ同じように転用できるので、学習量的には軽かった。
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で完結できると感じました。労働生産性を高めるためにも積極的につかっていきたいものです。