Rでデータマイニング

Rで疑問に思ったことなどを

forcats

最近見つけた私にとっての便利なパッケージforcats

もともとfactorの扱いが苦手。

今回困ったのはggplotで棒グラフを描いたときに件数順にソートしたくてさっとはできなかった。

使った関数はfct_infreq

levelsを件数順にしてくれます。

もともとの変数

library(knitr)
table(chickwts$feed)
## 
##    casein horsebean   linseed  meatmeal   soybean sunflower 
##        12        10        12        11        14        12
levels(chickwts$feed)
## [1] "casein"    "horsebean" "linseed"   "meatmeal"  "soybean"   "sunflower"

件数順にする

library(forcats)
rep <- fct_infreq(chickwts$feed)
table(rep)
## rep
##   soybean    casein   linseed sunflower  meatmeal horsebean 
##        14        12        12        12        11        10
levels(rep)
## [1] "soybean"   "casein"    "linseed"   "sunflower" "meatmeal"  "horsebean"

ggplot2 windows環境でフォントを指定したい

グラフを書くときにフォントに悩みのはもう終わりにしたい。

メイリオフォントを指定する。

でもメイリオはあまり好みじゃないな。

library(ggplot2)
windowsFonts("MEI"=windowsFont("Meiryo"))
qplot(
  data=iris,
  Sepal.Length,
  Sepal.Width,
  main="メイリオフォントだよ"
) + theme_bw(base_size = 16, base_family = "MEI")

私の環境でインストールしているTakaoフォントを試してみる。

windowsFonts("TP"=windowsFont("Takaoゴシック"))
windowsFonts("TPG"=windowsFont("Takao Pゴシック"))
windowsFonts("TM"=windowsFont("Takao明朝"))
windowsFonts("TPM"=windowsFont("Takao P明朝"))

qplot(
  data=iris,
  Sepal.Length,
  Sepal.Width,
  main="Takao P明朝だよ"
) + theme_bw(base_size = 16, base_family = "TPM")

開発環境を少し変更

開発環境を再考して少しだけ自分の好みに近づいた。

  1. jupyter R
  2. vim key binding

1. jupyter R

いままではRstudioを使用していたが満足していなかった。一つ目の理由は余計な機能が多すぎること。2つ目はpythonと異なる2つの開発環境を併用するのが面倒。

これで一つの環境でpythonとRの開発ができる。

2. vim key binding

jupyterの拡張機能vimのキーに変更できることを知ったので導入。

グラフについてとても参考になった。 ダメ出し:ちゃんとしたグラフを描こう

偶然見つけた下記のページ。

ダメ出し:ちゃんとしたグラフを描こう

グラフを見ると明らかに分かりやすくなっている。 ggplot2はだめなのか? ggplot2で下のグラフはつくれないのか?

話は変わって愚痴っぽくなるがRで新しいことを試そうとすると躓くことが多くてなかなか捗らない。 一方pythonは不慣れなので時間を要するかと思いきやすんなり終わることが多々ある。

もしかしてRよりpythonのほうが生産性が高い?

まあ使い分けということにしておこう。

tm:言語処理とちょっとだけ嵌る

嵌ったこと

言言処理の前処理が必要でtmパッケージを使用しました。
参考情報は下記

Basic Text Mining in R

stemDocument関数で英単語が集計しやすい形で出力されるはずなのにされない。

参考情報に従って問題を確認

データはacqを使用。

前処理内容
1. ピリオド除去
2. 数字除去
3. 小文字化
4. 前置詞、代名詞など集計に必要のない単語の除去
5. 複数形など除去
6. 余分な空白や改行など除去

実行結果を確認すると元々PlainTextDocumentクラスがcharacterクラスになり、
複数形などの除去が失敗しています。

library(tm)  
## Loading required package: NLP
# 使用するデータacqの確認  
data(acq)  
class(acq[[1]])  
## [1] "PlainTextDocument" "TextDocument"
acq[[1]]$content  
## [1] "Computer Terminal Systems Inc said\nit has completed the sale of 200,000 shares of its common\nstock, and warrants to acquire an additional one mln shares, to\n<Sedio N.V.> of Lugano, Switzerland for 50,000 dlrs.\n    The company said the warrants are exercisable for five\nyears at a purchase price of .125 dlrs per share.\n    Computer Terminal said Sedio also has the right to buy\nadditional shares and increase its total holdings up to 40 pct\nof the Computer Terminal's outstanding common stock under\ncertain circumstances involving change of control at the\ncompany.\n    The company said if the conditions occur the warrants would\nbe exercisable at a price equal to 75 pct of its common stock's\nmarket price at the time, not to exceed 1.50 dlrs per share.\n    Computer Terminal also said it sold the technolgy rights to\nits Dot Matrix impact technology, including any future\nimprovements, to <Woodco Inc> of Houston, Tex. for 200,000\ndlrs. But, it said it would continue to be the exclusive\nworldwide licensee of the technology for Woodco.\n    The company said the moves were part of its reorganization\nplan and would help pay current operation costs and ensure\nproduct delivery.\n    Computer Terminal makes computer generated labels, forms,\ntags and ticket printers and terminals.\n Reuter"
# ピリオド除去  
acqr <- tm_map(acq, removePunctuation)    
  
# 数字除去  
acqr <- tm_map(acqr, removeNumbers)    
  
# 小文字化  
acqr <- tm_map(acqr, tolower)  
  
# 前置詞、代名詞など集計に必要のない単語の除去  
acqr <- tm_map(acqr, removeWords, stopwords("english"))  
  
# 複数形など除去  
library(SnowballC)   
acqr <- tm_map(acqr, stemDocument)  
  
# 余分な空白や改行など除去  
acqr <- tm_map(acqr, stripWhitespace)  
  
# 前処理結果確認  
class(acqr[[1]])  
## [1] "character"
acqr[[1]]  
## [1] "computer terminal systems inc said completed sale shares common stock warrants acquire additional one mln shares sedio nv lugano switzerland dlrs company said warrants exercisable five years purchase price dlrs per share computer terminal said sedio also right buy additional shares increase total holdings pct computer terminals outstanding common stock certain circumstances involving change control company company said conditions occur warrants exercisable price equal pct common stocks market price time exceed dlrs per share computer terminal also said sold technolgy rights dot matrix impact technology including future improvements woodco inc houston tex dlrs said continue exclusive worldwide licensee technology woodco company said moves part reorganization plan help pay current operation costs ensure product delivery computer terminal makes computer generated labels forms tags ticket printers terminals reut"

原因

原因はstemDocument関数ではなく、tolower関数。
tmパッケージ内の関数はPlainTextDocumentクラスに適用する関数で問題ないが
tolower関数はbaseパッケージ内の関数。返り値が“character”になってしまう。

具体例

library(tm)  
  
# テスト用の文字列を作成  
testtext <- PlainTextDocument("aAa 12 aFoiewjfa")  
class(testtext)  
## [1] "PlainTextDocument" "TextDocument"
# 数字を除去してみる  
r1 <- removeNumbers(testtext)  
class(r1)  
## [1] "PlainTextDocument" "TextDocument"
r1$content  
## [1] "aAa  aFoiewjfa"
# 小文字化してみる  
  
r2 <- tolower(testtext)  
class(r2)  
## [1] "character"
r2  
## [1] "aaa 12 afoiewjfa"

解決策

content_transformerをかますことで解決できます。

data(acq)  
  
# ピリオド除去  
acqr <- tm_map(acq, removePunctuation)    
  
# 数字除去  
acqr <- tm_map(acqr, removeNumbers)    
  
# 小文字化  
acqr <- tm_map(acqr, content_transformer(tolower))  
  
# 前置詞、代名詞など集計に必要のない単語の除去  
acqr <- tm_map(acqr, removeWords, stopwords("english"))  
  
# 複数形など除去  
library(SnowballC)   
acqr <- tm_map(acqr, stemDocument)  
  
# 余分な空白や改行など除去  
acqr <- tm_map(acqr, stripWhitespace)  
  
# 前処理結果確認  
class(acqr[[1]])  
## [1] "PlainTextDocument" "TextDocument"
acqr[[1]]$content  
## [1] "comput termin system inc said complet sale share common stock warrant acquir addit one mln share sedio nv lugano switzerland dlrs compani said warrant exercis five year purchas price dlrs per share comput termin said sedio also right buy addit share increas total hold pct comput termin outstand common stock certain circumst involv chang control company compani said condit occur warrant exercis price equal pct common stocks market price time exceed dlrs per share comput termin also said sold technolgi right dot matrix impact technolog includ future improv woodco inc houston tex dlrs said continu exclusive worldwid license technolog woodco compani said move part reorganization plan help pay current oper cost ensure product delivery comput termin make comput generat label forms tag ticket printer terminals reuter"

今更欠損値対応:mlr

欠損値対応

TREEベースのアルゴリズムでモデルを構築することが多いので欠損値対応は必要ないと思っていた。 今回はモデルを構築する前にクラスタリングすることにしてkmeansのため欠損値対応が必要になった。

平均値に置換

難しいことをしたいわけではなく、単純に平均値に置換したい。 平均値に置換するのは良い方法ではないみたいだけどわかりやすい。

欠損値対応のパッケージ

毎回、コードを書くのは大変なのでパッケージで対応したい。 miやmiceはもっと小難しいことをしているみたいで平均値置換にむいてない?

mlrが良さそう.

mlrが良さそう。 mlrはcaretと似たようなパッケージでモデル構築のためのパッケージ。 前処理で欠損値対応が可能。

欠損値対応

サンプルコード

# https://mlr-org.github.io/mlr-tutorial/release/html/create_imputation/index.html
library(mlr)
## Loading required package: BBmisc
## Loading required package: ggplot2
## Loading required package: ParamHelpers
data(airquality)
summary(airquality)
##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
##  3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
##  NA's   :37       NA's   :7                                       
##      Month            Day      
##  Min.   :5.000   Min.   : 1.0  
##  1st Qu.:6.000   1st Qu.: 8.0  
##  Median :7.000   Median :16.0  
##  Mean   :6.993   Mean   :15.8  
##  3rd Qu.:8.000   3rd Qu.:23.0  
##  Max.   :9.000   Max.   :31.0  
## 
imp <- impute(
    airquality, 
    cols = list(
        Ozone = imputeMean(), 
        Solar.R = imputeMean()
    )
)
head(imp$data, 10)
##       Ozone  Solar.R Wind Temp Month Day
## 1  41.00000 190.0000  7.4   67     5   1
## 2  36.00000 118.0000  8.0   72     5   2
## 3  12.00000 149.0000 12.6   74     5   3
## 4  18.00000 313.0000 11.5   62     5   4
## 5  42.12931 185.9315 14.3   56     5   5
## 6  28.00000 185.9315 14.9   66     5   6
## 7  23.00000 299.0000  8.6   65     5   7
## 8  19.00000  99.0000 13.8   59     5   8
## 9   8.00000  19.0000 20.1   61     5   9
## 10 42.12931 194.0000  8.6   69     5  10
str(imp)
## List of 2
##  $ data:'data.frame':    153 obs. of  6 variables:
##   ..$ Ozone  : num [1:153] 41 36 12 18 42.1 ...
##   ..$ Solar.R: num [1:153] 190 118 149 313 186 ...
##   ..$ Wind   : num [1:153] 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##   ..$ Temp   : int [1:153] 67 72 74 62 56 66 65 59 61 69 ...
##   ..$ Month  : int [1:153] 5 5 5 5 5 5 5 5 5 5 ...
##   ..$ Day    : int [1:153] 1 2 3 4 5 6 7 8 9 10 ...
##  $ desc:List of 9
##   ..$ target              : chr(0) 
##   ..$ features            : chr [1:6] "Ozone" "Solar.R" "Wind" "Temp" ...
##   ..$ lvls                : Named list()
##   ..$ impute              :List of 2
##   .. ..$ Ozone  :List of 2
##   .. .. ..$ impute:function (data, target, col, const)  
##   .. .. ..$ args  : num 42.1
##   .. ..$ Solar.R:List of 2
##   .. .. ..$ impute:function (data, target, col, const)  
##   .. .. ..$ args  : num 186
##   ..$ dummies             : chr(0) 
##   ..$ dummy.type          : chr "factor"
##   ..$ force.dummies       : logi FALSE
##   ..$ impute.new.levels   : logi TRUE
##   ..$ recode.factor.levels: logi TRUE
##   ..- attr(*, "class")= chr "ImputationDesc"
summary(imp$data)
##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
##  1st Qu.: 21.00   1st Qu.:120.0   1st Qu.: 7.400   1st Qu.:72.00  
##  Median : 42.13   Median :194.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
##  3rd Qu.: 46.00   3rd Qu.:256.0   3rd Qu.:11.500   3rd Qu.:85.00  
##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
##      Month            Day      
##  Min.   :5.000   Min.   : 1.0  
##  1st Qu.:6.000   1st Qu.: 8.0  
##  Median :7.000   Median :16.0  
##  Mean   :6.993   Mean   :15.8  
##  3rd Qu.:8.000   3rd Qu.:23.0  
##  Max.   :9.000   Max.   :31.0

そろそろ最適化問題を勉強したい:optim

最適化問題をおろそかにしてきた

実務ではパラメータ推計やモデル構築で事足りて最適化問題をおろそかにしてきた。

最適化が必要な状況

順序予測を回帰でモデル構築した後に予測値を離散化して順序にしたい。
離散化のカットオフ値を決めるときに最適化が使えるらしい。
モデルの予測値が得られた後に現実に合わせるのに最適化は使えそう。

難しいことは置いといて簡単な考え方

難しいことはわからないので簡単イメージ。 高校で勉強した極値問題と考えると自分のなかではイメージしやすい。

簡単な問題:2次関数

簡単な問題。

x^2+3x+1
3x3

極値を求める。

# 関数定義
f <- function(x) {
  x ^ 2 + 3 * x + 1
  }
f(3)
## [1] 19
library(ggplot2)
p <- ggplot(data=data.frame(X=c(-3,3)), aes(x=X))
p <- p + stat_function(fun=f)
p

# 最小化 グラフを見ると下に凸でx=-1.5ぐらいでy=-1.5ぐらいが最小値。
実際にoptim関数で計算する。

optim(par = -2, fn = f, method = "L-BFGS-B", lower = -3, upper = 3)
## $par
## [1] -1.5
## 
## $value
## [1] -1.25
## 
## $counts
## function gradient 
##        4        4 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

結果はx=-1.5, y=-1.25。

最大化

グラフを見ると下に凸でx=3でy=19が最小値。
実際にoptim関数で計算する。
最大化する場合は control=list(fnscale=-1) を指定する。

optim(par = -2, fn = f, method = "L-BFGS-B", lower = -3, upper = 3, control=list(fnscale=-1))
## $par
## [1] -3
## 
## $value
## [1] 1
## 
## $counts
## function gradient 
##        3        3 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

結果はx=-3, y=1。結果が間違い。
これは初期値 par = -2 を指定したことが原因。
極値複数ある場合は局所最適に陥り全体最適にならない場合がある。
このような場合は初期値を複数試してみて全体最適を探る必要がある。
初期値 x = 0 で試してみる。

optim(par = 0, fn = f, method = "L-BFGS-B", lower = -3, upper = 3, control=list(fnscale=-1))
## $par
## [1] 3
## 
## $value
## [1] 19
## 
## $counts
## function gradient 
##        2        2 
## 
## $convergence
## [1] 0
## 
## $message
## [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

疑問:離散の最適化

ふと、離散の最適化について疑問に思った。
調べてみるとやはり難しい問題のようで、おいおい必要に応じて勉強したい。

今、私が求めている最適化問題

  1. 目的変数が順序尺度
  2. 回帰でモデルを構築
  3. モデルで計算した予測値を順序尺度にしたい
  4. 最も適したカットオフ値を求めたい
  5. モデルの評価指標は重み付きkappa

続きはまた後で。今回はここまで