R – ページ 4 – かものはしの分析ブログ

RStanで学部時代の研究を振り返ってみる

研究概要

大学時代に実験経済学で行った実験結果のデータがUSBに入っていたので、振り返って分析などをしてみたいと思います。

研究目的
　ピア効果に関して、競争相手が自分よりも秀でた人がいいのか劣った人がいいのかを確かめる。

実験方法
・1分間で100マス計算を2セット解いてもらう。（めちゃ速い人には3枚目も渡した）
・実験開始後、実験対象のクラスによって、途中で「平均的なクラスは○○マスまで進んでいます！」とアナウンスします。アナウンスすることで、競争相手のレベルを知り、焦るなり余裕を感じるなりしてもらおうという計画です。
なお、対照群はアナウンスをしていません。アナウンス内容は「平均告知（１８秒）、上告知（１５秒）、超上告知（１２秒）、下告知（２０秒）」と4パターンとなります。
・計算が間違っているものは加点しません。

実験対象
　某国立大学の経済学部生の1~2年の必修科目履修者217名（先生に交渉して授業の開始5分を頂いて実験を行いました。）
　内部進学やスポ専などがない分、計算能力的にある程度近い集団ではないかと思われます。

検証方法
　アナウンスごとに100マス計算の点数の水準が変わりうるのかを回帰分析などで判断。

データ可視化

以下、実験カテゴリごとの略記です。
下告知（２０秒）・・・slow_20
上告知（１５秒）・・・fast_15
超上告知（１２秒）・・・fastest_12
平均告知（１８秒）・・・average_18
対照群・・・baseline

データ構造の確認です。

> str(dataset)
'data.frame':	217 obs. of  4 variables:
 $ categories    : chr  "baseline" "baseline" "baseline" "baseline" ...
 $ points        : int  72 79 81 81 98 99 100 101 102 104 ...
 $ errors        : int  0 0 0 0 2 1 0 0 1 0 ...
 $ genuine_points: int  72 79 81 81 96 98 100 101 101 104 ...

> str(dataset)

'data.frame': 217 obs. of 4 variables:

$ categories : chr "baseline" "baseline" "baseline" "baseline" ...

$ points : int 72 79 81 81 98 99 100 101 102 104 ...

$ errors : int 0 0 0 0 2 1 0 0 1 0 ...

$ genuine_points: int 72 79 81 81 96 98 100 101 101 104 ...

平均値、中央値、標準偏差、サンプルサイズを出してみます。

> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points),
+                                                median=median(genuine_points),
+                                                stdev=sd(genuine_points),
+                                                sample=n())
# A tibble: 5 × 5
  categories  average median    stdev sample
       <chr>    <dbl>  <dbl>    <dbl>  <int>
1 average_18 120.8605    117 28.08455     43
2   baseline 121.0000    121 22.19109     46
3    fast_15 123.0357    126 24.32355     56
4 fastest_12 123.6774    125 24.09756     31
5    slow_20 126.3902    126 23.20547     41

> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points),

+ median=median(genuine_points),

+ stdev=sd(genuine_points),

+ sample=n())

# A tibble: 5 × 5

categories average median stdev sample

1 average_18 120.8605 117 28.08455 43

2 baseline 121.0000 121 22.19109 46

3 fast_15 123.0357 126 24.32355 56

4 fastest_12 123.6774 125 24.09756 31

5 slow_20 126.3902 126 23.20547 41

中央値で見てみると、baselineに対してわずかですが点数に違いがありそうに見えます。

実験種別で点数に関するヒストグラムと確率密度関数を確認してみます。

library(ggplot2)

g <- ggplot(data = dataset,
            aes(x = genuine_points,
                y = ..density..)) +
            geom_histogram(alpha = 0.5,position = "identity") +
            geom_density(alpha = 0)
g + facet_wrap(~categories,nrow=5)

library(ggplot2)

g <- ggplot(data = dataset,

aes(x = genuine_points,

y = ..density..)) +

geom_histogram(alpha = 0.5,position = "identity") +

geom_density(alpha = 0)

g + facet_wrap(~categories,nrow=5)

baselineが多峰性がありそうなのが気になります。average_18は低そうに見えますね。

RStanで重回帰

『StanとRでベイズ統計モデリング』にあるコードを参考にしています。正規分布を事前分布にした線形回帰モデルです。
被説明変数が点数、説明変数が実験種別のダミー変数だけからなります。

library(rstan)
library(dummies)

dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories"))

analytical_dataset <- dummies %>% select(categories_average_18,
                                         categories_fast_15,
                                         categories_fastest_12,
                                         categories_slow_20,
                                         genuine_points)


data <- list(N=nrow(analytical_dataset),
             average_18=analytical_dataset$categories_average_18,
             fast_15=analytical_dataset$categories_fast_15,
             fastest_12=analytical_dataset$categories_fastest_12,
             slow_20=analytical_dataset$categories_slow_20,
             genuine_points=analytical_dataset$genuine_points)


stan_code <- "
data{
int N; //the number of student
int<lower=0> genuine_points[N];
real<lower=0, upper=1> average_18[N];
real<lower=0, upper=1> fast_15[N];
real<lower=0, upper=1> fastest_12[N];
real<lower=0, upper=1> slow_20[N];
}

parameters{
real b1;
real b2;
real b3;
real b4;
real b5;
real<lower=0> sigma;
}

transformed parameters{
real mu[N];
for(n in 1:N)
mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n];
}

model{
for(n in 1:N)
genuine_points[n] ~ normal(mu[n], sigma);
}
"

fit <- stan(model_code =stan_code, data=data, seed=1234)
fit.summary <-data.frame(summary(fit)$summary)
head(fit.summary,6)

library(rstan)

library(dummies)

dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories"))

analytical_dataset <- dummies %>% select(categories_average_18,

categories_fast_15,

categories_fastest_12,

categories_slow_20,

genuine_points)

data <- list(N=nrow(analytical_dataset),

average_18=analytical_dataset$categories_average_18,

fast_15=analytical_dataset$categories_fast_15,

fastest_12=analytical_dataset$categories_fastest_12,

slow_20=analytical_dataset$categories_slow_20,

genuine_points=analytical_dataset$genuine_points)

stan_code <- "

data{

int N; //the number of student

int<lower=0> genuine_points[N];

real<lower=0, upper=1> average_18[N];

real<lower=0, upper=1> fast_15[N];

real<lower=0, upper=1> fastest_12[N];

real<lower=0, upper=1> slow_20[N];

}

parameters{

real b1;

real b2;

real b3;

real b4;

real b5;

real<lower=0> sigma;

}

transformed parameters{

real mu[N];

for(n in 1:N)

mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n];

}

model{

for(n in 1:N)

genuine_points[n] ~ normal(mu[n], sigma);

}

fit <- stan(model_code =stan_code, data=data, seed=1234)

fit.summary <-data.frame(summary(fit)$summary)

head(fit.summary,6)

結果

traceplot(fit)でMCMCのサンプリング結果を確認する。

収束しているように見えます。

以下推定結果ですが、残念ながらベイズ予測区間において符号の逆転が起きていないものはなかったので、アナウンスによる効果があるとは言えないようです。ただ、slow_20の係数がおしいですね。少なくとも他の実験種別よりも、アナウンス効果があるかもしれないという考察に止まりそうです。

> head(fit.summary,6)
             mean   se_mean       sd      X2.5.       X25.        X50.       X75.     X97.5.    n_eff      Rhat
b1    121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152
b2     -0.1195318 0.1326887 4.957153 -10.009827  -3.372795  -0.1304577   3.117195   9.694184 1395.715 1.0004034
b3      1.9440496 0.1248433 4.661571  -7.494857  -1.043060   1.9756929   5.007503  10.827640 1394.228 1.0000830
b4      2.5460694 0.1436941 5.447930  -8.026616  -1.128587   2.4997953   6.220891  13.044467 1437.426 1.0008688
b5      5.3029965 0.1351662 4.988323  -4.479865   2.004280   5.3108476   8.543292  15.458379 1361.988 1.0004238
sigma  24.6158717 0.0211252 1.150481  22.501612  23.811858  24.5900106  25.368057  26.970197 2965.905 1.0017414

> head(fit.summary,6)

mean se_mean sd X2.5. X25. X50. X75. X97.5. n_eff Rhat

b1 121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152

b2 -0.1195318 0.1326887 4.957153 -10.009827 -3.372795 -0.1304577 3.117195 9.694184 1395.715 1.0004034

b3 1.9440496 0.1248433 4.661571 -7.494857 -1.043060 1.9756929 5.007503 10.827640 1394.228 1.0000830

b4 2.5460694 0.1436941 5.447930 -8.026616 -1.128587 2.4997953 6.220891 13.044467 1437.426 1.0008688

b5 5.3029965 0.1351662 4.988323 -4.479865 2.004280 5.3108476 8.543292 15.458379 1361.988 1.0004238

sigma 24.6158717 0.0211252 1.150481 22.501612 23.811858 24.5900106 25.368057 26.970197 2965.905 1.0017414

分布でも見てみます。アナウンス効果が0を確実に超えているとは言えないですね。

library (reshape)
library(dplyr)
library(ggplot2)
post   <- extract (fit, permuted = F)
m.post <- melt (post)
m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5"))
graph  <- ggplot (m.post, aes(x = value))
graph  <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw() 
plot (graph)

library (reshape)

library(dplyr)

library(ggplot2)

post <- extract (fit, permuted = F)

m.post <- melt (post)

m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5"))

graph <- ggplot (m.post, aes(x = value))

graph <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw()

plot (graph)

結局、学部時代のレポートと結論は変わらないのですが、係数が0よりも大きい確率という観点で結果に向き合えたのは良かったと思います。

参考文献

StanとRでベイズ統計モデリング (Wonderful R)

PythonやRを用いたアルゴリズム取引・株価分析まとめ

はじめに

まわりでシステムトレードや株価の機械学習による予測などに関心が高まってきたので、私も少し調べてみようと思いPythonやRで行われた分析・実装の事例を集めてみました。自分の資産を突っ込む気にはなれないですが、事例を知っておくだけ知っておきたいですね。

調査法

Google検索において以下のクエリで上位に表示されたサイトを中心にまとめました。
「python 機械学習株価」
「python 機械学習為替」
「python アルゴリズム取引」
「python machine learning stock price」
（同様にRも調べました。）

アルゴリズムトレードの理論から学ぶ

システムトレードで億万長者になるぞ! coursera で Computational Investing Part I を受けた
こちらはCourseraで開かれていたシステムトレード向けの講座を受けた方のブログです。Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Riskはもっとも重要なアルゴリズムトレードの本とされているようです。出版年が1999年と結構古いですが。
Active Portfolio Management: A Quantitative Approach for Producing Superior Returns and Selecting Superior Returns and Controlling Risk (McGraw-Hill Library of Investment and Finance)

Rを用いた株価データのモニタリング

Stock Analysis using R
こちらはアルゴリズム取引とかではないですが、Rで株価を簡単に取得できるquantmodライブラリについて紹介されています。

library(quantmod)
getSymbols("MSFT")
chartSeries(MSFT, subset="last 3 months")
addBBands()

library(quantmod)

getSymbols("MSFT")

chartSeries(MSFT, subset="last 3 months")

addBBands()

以上のコードを実行しただけで時系列データがボックスプロットで描写されました。非常に便利そうです。複数のグラフを並べてディスプレーで見てみたいですね。図ではテクニカル分析で使われるボリンジャーバンドを一発で出してくれています。傾向を掴むための分析に役立つのではないでしょうか。

RSIで株価の連動性を見る

過去のデータからビッグデータ分析で株価を予測する
ここでは、主にpandasを用いて株価の分析を行っています。RSI（Relative Strength Index）という「一定の期間変動幅の中でどれ位株価が上昇しているのか、下落しているのかをはかるもの」を計算して日経平均との相関を見ています。データさえあれば、個別株の市場連動性を見るぶんにはpandasで十分に分析できそうです。

アルゴリズムトレードのシステム構成

自動トレードボット
こちらでは、仮想通貨取引のための自動トレードボットの作成のための手順などが書かれています。
工程としては
・仮想通貨価格データ取得
・バックテストの実施
・明日の価格予想
・学習パラメータの最適化
・結果のメール送信
・APIを使った成行トレード
・ジョブのスケジューリング
などとなっているようです。
完全自動化しようと思うと、作るのは大変そうですね。

WEB屋の自分が機械学習株価予想プログラムを開発した結果
こちらはよりシンプルで
「チャートを見る限り、誰がどうみても今日値上がりする銘柄」を検索して、毎朝Slackに通知してくれるシステム
をPythonで作成されています。最終的に資産を溶かす形になっているようです。

アルゴリズムトレード向けのPythonライブラリ

pythonのアルゴリズムトレードライブラリ
こちらでは
zipline（ http://www.zipline.io/ ）
PyAlgoTrade( http://gbeced.github.io/pyalgotrade/ )
pybacktest( https://github.com/ematvey/pybacktest )
backtrader( https://github.com/mementum/backtrader )
というライブラリが紹介されています。

遺伝的アルゴリズムで為替の自動売買

pythonと遺伝的アルゴリズムで作るFX自動売買システムその1
遺伝的アルゴリズムでFX自動売買その2 進化する売買AIの実装
 遺伝的アルゴリズムでFX自動売買その3 OandaAPIで実際に取引
こちらは、外国為替の取引をPythonで自動化させた試みです。最終的に負けてしまってはいますがシステム構成についても詳しく書かれているので、勉強になりました。データはOANDAというサービスが提供しているAPIを用いて取得し、GA（遺伝的アルゴリズム）を使って為替の売買タイミングを決めているようです。実際の売買にもOandaAPIというものを利用して完全に自動化させています。

決定木による株価リターンの予測

機械学習で未来を予測する – scikit-learn の決定木で未来の株価を予測
こちらはPythonでscikit-learnを用いて、決定木による株価の予測をされています。目的変数としてはリターンインデックスをおいています。前処理においてpandasが使われています。

ランダムフォレストによる株価の予測

Pythonで機械学習を使った株価予測のコードを書こう
こちらでは、ランダムフォレストを用いた機械学習で、ETFなどのデータを特徴量にして個別の株価予測を行っています。

アカデミックサイドでの研究事例

Stock Price Forecasting Using Information from Yahoo Finance and Google Trend
UC　Berkeleyの経済学部生の研究です。こちらはR言語で、ヤフーファイナンスとGoogleトレンドの情報を用いて株価を予測している研究です。これまでの時系列手法よりもパフォーマンスが良いそうです。

最後に

事例を集めてみて、どこまでのレベルのものに手を出すべきか悩ましいと思いました。完全自動化は維持するコストがかかりそうですし、ロスが増大しないか心配です。自分としては情報を自動で取得し、リターンが発生する確率の高そうな銘柄をサジェストしてくれるレベルで十分な気がします。

追記

人工知能学会の全国大会2017でAIを用いた株式運用に関する研究がなされていました。深層学習を用いたものが多そうです。
株価変動パターンの類似性を用いた株価予測
 深層学習と高頻度データを用いた株式注文状況の推定
 LSTM を用いた株価変動予測
 深層学習を用いた株価予測の分析

Rで株価の時系列データを簡単に集計する

はじめに

かねてからYahoo!ファイナンスなどの情報に、株価について期間を指定したリターンや標準偏差や回帰分析のあてはまりなどが無いので、そのようなデータを簡単に見れるようにしたいという願望がありました。加えて、できるだけ多くの銘柄の情報をみたいという願望もありました。Rを使えばその全てが簡単に実現できるのでチャレンジしてみました。

・データ収集
・データテーブルの整形
・成長率の計算
・標準偏差の計算
・自己回帰の計算
と非常にシンプルな工程ですが、3600銘柄分のこれらの情報を一気に手に入れれるというのは魅力的だなぁと思います。

データの収集

RFinanceYJというパッケージをインストールして、以下のrfinance_module.Rのスクリプトを置いておきます。

library(RFinanceYJ)

#API
quoteStockTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.stock &lt;- function(quote.table.item){
    if( xmlSize(quote.table.item) &lt; 5) return(NULL) 
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o &lt;- as.number(xmlValue(quote.table.item[[2]]))
    h &lt;- as.number(xmlValue(quote.table.item[[3]]))
    l &lt;- as.number(xmlValue(quote.table.item[[4]]))
    c &lt;- as.number(xmlValue(quote.table.item[[5]]))
    v &lt;- ifelse(xmlSize(quote.table.item) &gt;= 6,as.number(xmlValue(quote.table.item[[6]])),0)
    a &lt;- ifelse(xmlSize(quote.table.item) &gt;= 7,as.number(xmlValue(quote.table.item[[7]])),0)
    return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))
  }
  return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))
}
quoteFundTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.fund &lt;- function(quote.table.item){
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    if(time.interval=='monthly'){
      d &lt;- endOfMonth(d)
    }
    c &lt;- as.number(xmlValue(quote.table.item[[2]]))
    v &lt;- as.number(xmlValue(quote.table.item[[3]]))
    return(data.frame(date=d,constant.value=c,NAV=v))
  }
  return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))
}
quoteFXTsData &lt;- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval &lt;- substr(time.interval,1,1)
  function.fx &lt;- function(quote.table.item){
    d &lt;- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o &lt;- as.number(xmlValue(quote.table.item[[2]]))
    h &lt;- as.number(xmlValue(quote.table.item[[3]]))
    l &lt;- as.number(xmlValue(quote.table.item[[4]]))
    c &lt;- as.number(xmlValue(quote.table.item[[5]]))
    return(data.frame(date=d,open=o,high=h,low=l,close=c))
  }
  return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))
}
######  private functions  #####
#get time series data from Yahoo! Finance.
quoteTsData &lt;- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){
  r &lt;- NULL
  result.num &lt;- 51
  financial.data &lt;- data.frame(NULL)
  start &lt;- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&amp;sy=\\1&amp;sm=\\2&amp;sd=\\3",since))
  end   &lt;- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&amp;ey=\\1&amp;em=\\2&amp;ed=\\3",date.end))
  
  if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")
  
  extractQuoteTable &lt;- function(r,type){
    if(type %in% c("fund","fx")){
      tbl &lt;- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]
    }
    else{
      tbl &lt;- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]
    }
    return(tbl)
  }
  
  #while( result.num &gt;= 51 ){
  while(1){
    start.num &lt;- start.num + 1
    quote.table &lt;- NULL
    quote.url &lt;- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&amp;p=',start.num,'&amp;tm=',substr(time.interval,1,1),sep="")
    try(r&lt;-htmlParse(quote.url))
    if( is.null(r) ) stop(paste("Can not access :", quote.url))
    
    try( quote.table &lt;- xpathApply(r,"//table")[[2]], TRUE )
    
    quote.size&lt;-xmlSize(quote.table)
    #cat(paste("size:",quote.size))
    if(xmlSize(quote.table)&lt;=1){
      return (financial.data)
    }
    if( is.null(quote.table) ){
      if( is.null(financial.data) ){
        stop(paste("Can not quote :", x))
      }else{
        financial.data &lt;- financial.data[order(financial.data$date),]
        return(financial.data)
      }
    }
    
    size &lt;- xmlSize(quote.table)
    for(i in 2:size){
      financial.data &lt;- rbind(financial.data,function.financialproduct(quote.table[[i]]))
    }
    
    #result.num &lt;- xmlSize(quote.table)
    Sys.sleep(1)
  }
  financial.data &lt;- financial.data[order(financial.data$date),]
  return(financial.data)  
}
#convert string formart date to POSIXct object
convertToDate &lt;- function(date.string,time.interval)
{
  #data format is different between monthly and dialy or weekly
  if(any(time.interval==c('d','w'))){
    result &lt;- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)
  }else if(time.interval=='m'){
    result &lt;- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)
  }
  return(as.POSIXct(result))
}
#convert string to number.
as.number &lt;- function(string)
{
  return(as.double(as.character(gsub("[^0-9.]", "",string))))
}
#return end of month date.
endOfMonth &lt;- function(date.obj)
{
  startOfMonth     &lt;- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")
  startOfNextMonth &lt;- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")
  return(startOfNextMonth-1)
}</pre>

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

library(RFinanceYJ)

#API

quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.stock <- function(quote.table.item){

if( xmlSize(quote.table.item) < 5) return(NULL)

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)

a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)

return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))

}

return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))

}

quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fund <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

if(time.interval=='monthly'){

d <- endOfMonth(d)

}

c <- as.number(xmlValue(quote.table.item[[2]]))

v <- as.number(xmlValue(quote.table.item[[3]]))

return(data.frame(date=d,constant.value=c,NAV=v))

}

return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))

}

quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fx <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

return(data.frame(date=d,open=o,high=h,low=l,close=c))

}

return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))

}

###### private functions #####

#get time series data from Yahoo! Finance.

quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){

r <- NULL

result.num <- 51

financial.data <- data.frame(NULL)

start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))

end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))

if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")

extractQuoteTable <- function(r,type){

if(type %in% c("fund","fx")){

tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]

}

else{

tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]

}

return(tbl)

}

#while( result.num >= 51 ){

while(1){

start.num <- start.num + 1

quote.table <- NULL

quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")

try(r<-htmlParse(quote.url))

if( is.null(r) ) stop(paste("Can not access :", quote.url))

try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )

quote.size<-xmlSize(quote.table)

#cat(paste("size:",quote.size))

if(xmlSize(quote.table)<=1){

return (financial.data)

}

if( is.null(quote.table) ){

if( is.null(financial.data) ){

stop(paste("Can not quote :", x))

}else{

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

size <- xmlSize(quote.table)

for(i in 2:size){

financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))

}

#result.num <- xmlSize(quote.table)

Sys.sleep(1)

}

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

#convert string formart date to POSIXct object

convertToDate <- function(date.string,time.interval)

{

#data format is different between monthly and dialy or weekly

if(any(time.interval==c('d','w'))){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)

}else if(time.interval=='m'){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)

}

return(as.POSIXct(result))

}

#convert string to number.

as.number <- function(string)

{

return(as.double(as.character(gsub("[^0-9.]", "",string))))

}

#return end of month date.

endOfMonth <- function(date.obj)

{

startOfMonth <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")

startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")

return(startOfNextMonth-1)

}</pre>

以下のコードで2016-01-01から現在に至るまでの、指定した株価コードの株価情報を集めてきます。3600銘柄の株価コードを用意する必要があるので、集めてcodelist.csvで保存しておきます。銘柄コードの情報はここでは掲載しません。（適当なサイトにころがっています。）

source("rfinance_module.R")

codelist <- read.csv("codelist.csv")
since    <- "2016-01-01"

stock_price_database <- 0
for(i in 1:nrow(codelist)){
  if(i ==1){
    dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)
    code    <- matrix(codelist[i,1], nrow(dataset), 1)
    dataset <- cbind(code, dataset) 
    stock_price_database <- dataset
  }
  else{
    dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)
    code    <- matrix(codelist[i,1], nrow(dataset), 1)
    dataset <- cbind(code, dataset) 
    stock_price_database <- rbind(stock_price_database,dataset)    
  }
}

source("rfinance_module.R")

codelist <- read.csv("codelist.csv")

since <- "2016-01-01"

stock_price_database <- 0

for(i in 1:nrow(codelist)){

if(i ==1){

dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)

code <- matrix(codelist[i,1], nrow(dataset), 1)

dataset <- cbind(code, dataset)

stock_price_database <- dataset

}

else{

dataset <- quoteStockTsData(paste(codelist[i,1],".t",sep = ""),since=since)

code <- matrix(codelist[i,1], nrow(dataset), 1)

dataset <- cbind(code, dataset)

stock_price_database <- rbind(stock_price_database,dataset)

}

以下のコードを実行して、元のデータを調整済み終値について、縦が日付で横が銘柄のテーブルに整形します。

#make stock price data table
library(dplyr)
code_list <- data.frame(code=distinct(stock_price_database,code)$code)
date_list <- data.frame(date=distinct(stock_price_database,date)$date)

stock_dataset <- date_list
stock_price_database_adj <- stock_price_database[,which(colnames(stock_price_database) %in% c("code","date","adj_close"))]

for( i in 1:nrow(code_list)){ #nrow(code_list)
  each_stock <- filter(stock_price_database_adj, code == code_list$code[i])
  each_stock <- each_stock[,-1]
  each_stock <- distinct(each_stock, date,.keep_all = TRUE)
  colnames(each_stock)[2] <- paste("sdj",code_list$code[i],sep = "") 
  stock_dataset <- left_join(stock_dataset, each_stock, by = "date")
}

#make stock price data table

library(dplyr)

code_list <- data.frame(code=distinct(stock_price_database,code)$code)

date_list <- data.frame(date=distinct(stock_price_database,date)$date)

stock_dataset <- date_list

stock_price_database_adj <- stock_price_database[,which(colnames(stock_price_database) %in% c("code","date","adj_close"))]

for( i in 1:nrow(code_list)){ #nrow(code_list)

each_stock <- filter(stock_price_database_adj, code == code_list$code[i])

each_stock <- each_stock[,-1]

each_stock <- distinct(each_stock, date,.keep_all = TRUE)

colnames(each_stock)[2] <- paste("sdj",code_list$code[i],sep = "")

stock_dataset <- left_join(stock_dataset, each_stock, by = "date")

}

成長率

指定した期間での成長率を計算します。

#calculate growth rate

#start_date
start_date <- "2016-01-04"
#end_date
end_date <- "2017-02-17"

#return rate datatable
return_data <- filter(stock_dataset, date==start_date |date==end_date)
return_data <- return_data[,-1]
return_data <- rbind(return_data, (as.numeric(return_data[1,]) / as.numeric(return_data[2,])-1)*100)
return_data <- data.frame(code=colnames(return_data),t(return_data))
colnames(return_data) <- c("code",end_date,start_date,"return_rate")
return_data <- return_data[order(return_data$return_rate,decreasing = TRUE),]

#calculate growth rate

#start_date

start_date <- "2016-01-04"

#end_date

end_date <- "2017-02-17"

#return rate datatable

return_data <- filter(stock_dataset, date==start_date |date==end_date)

return_data <- return_data[,-1]

return_data <- rbind(return_data, (as.numeric(return_data[1,]) / as.numeric(return_data[2,])-1)*100)

return_data <- data.frame(code=colnames(return_data),t(return_data))

colnames(return_data) <- c("code",end_date,start_date,"return_rate")

return_data <- return_data[order(return_data$return_rate,decreasing = TRUE),]

標準偏差

ここではシンプルに標準偏差の計算を行います。

#volatility
volatility_data <- stock_dataset[,-1]
volatility_data <- log(volatility_data[-nrow(volatility_data),]) - log(volatility_data[-1,])

volatility <- data.frame(matrix(0,1,2))
colnames(volatility) <- c("code", "stdev")

for(i in 1:length(volatility_data)){
  volatility[i,1] <- colnames(volatility_data)[i]
  volatility[i,2] <- sd(volatility_data[,i],na.rm = TRUE)
}

#volatility

volatility_data <- stock_dataset[,-1]

volatility_data <- log(volatility_data[-nrow(volatility_data),]) - log(volatility_data[-1,])

volatility <- data.frame(matrix(0,1,2))

colnames(volatility) <- c("code", "stdev")

for(i in 1:length(volatility_data)){

volatility[i,1] <- colnames(volatility_data)[i]

volatility[i,2] <- sd(volatility_data[,i],na.rm = TRUE)

}

自己回帰による決定係数

自己回帰でフィッティングの良い銘柄はどれなのかを単純に知りたいので、3600銘柄分計算させます。

#AR(1) estimation

ar1_estimation <- stock_dataset[,-1]
ar1_rsuares <- data.frame(matrix(0,1,2))
colnames(ar1_rsuares) <- c("code", "adj_rsquare")

for(i in 1:length(ar1_estimation)){
  #objective variable
  ar1_y <- ar1_estimation[-nrow(ar1_estimation),i]
  #input variable
  ar1_x <- ar1_estimation[-1,i]
  #combine data
  ar1_data <- data.frame(cbind(ar1_y, ar1_x))
  
  #estimation
  if(class(try(lm(data = ar1_data, formula = ar1_y ~. )))=="lm"){
    model <- lm(data = ar1_data, formula = ar1_y ~. )
    model_sum <- summary(model)
    
    ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]
    ar1_rsuares[i,2] <- model_sum$adj.r.squared    
  }
  else{
    ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]
    ar1_rsuares[i,2] <- "NULL"
  }
}

#AR(1) estimation

ar1_estimation <- stock_dataset[,-1]

ar1_rsuares <- data.frame(matrix(0,1,2))

colnames(ar1_rsuares) <- c("code", "adj_rsquare")

for(i in 1:length(ar1_estimation)){

#objective variable

ar1_y <- ar1_estimation[-nrow(ar1_estimation),i]

#input variable

ar1_x <- ar1_estimation[-1,i]

#combine data

ar1_data <- data.frame(cbind(ar1_y, ar1_x))

#estimation

if(class(try(lm(data = ar1_data, formula = ar1_y ~. )))=="lm"){

model <- lm(data = ar1_data, formula = ar1_y ~. )

model_sum <- summary(model)

ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]

ar1_rsuares[i,2] <- model_sum$adj.r.squared

}

else{

ar1_rsuares[i,1] <- colnames(ar1_estimation)[i]

ar1_rsuares[i,2] <- "NULL"

}

以上より、見てみたかった指標の3600銘柄分を一気に計算することができました。

可視化

試しに一番成長率の高かった株価をプロットしてみます。

#visualization
library(ggplot2)
stock_code <- "sdj7612"

df <- data.frame(stock_dataset[,which(colnames(stock_dataset) %in% c("date", stock_code))])

g <- ggplot(df,aes(x= date, y= df[,2]))+
     geom_line(colour = "magenta",linetype = 2,size = 0.5)+
     geom_smooth (method = "lm")+
     xlab("Date")+
     ylab("Stock Price")+
     ggtitle("Stock Price")
plot(g)

#visualization

library(ggplot2)

stock_code <- "sdj7612"

df <- data.frame(stock_dataset[,which(colnames(stock_dataset) %in% c("date", stock_code))])

g <- ggplot(df,aes(x= date, y= df[,2]))+

geom_line(colour = "magenta",linetype = 2,size = 0.5)+

geom_smooth (method = "lm")+

xlab("Date")+

ylab("Stock Price")+

ggtitle("Stock Price")

plot(g)

おわりに

株価データにカジュアルにアクセスできる環境を手に入れることができたので、今後は様々な指標や時系列解析・業界ごとの分析を試してみたいと思います。

顧客生涯価値(CLV)の計算 with R

顧客生涯価値(CLV：Customer Lifetime Value)を計算してくれるRのコード（Calculating Customer Lifetime Value with Recency, Frequency, and Monetary (RFM)）があったので、今更感がありますが取り上げたいと思います。

目次

・顧客生涯価値の数式
・データセット
・関数
・データセットの読み込みと加工
・再購買率とRFMとの関係
・再購買率の推定
・顧客生涯価値の計算
・参考情報

顧客生涯価値の数式

まず、顧客生涯価値の数式は以下の通りです。

t：年や月などの期間
n：顧客が解約するまでの期間合計
r：保持率（1-解約率）
P(t)：t期に顧客から得られる収益
d：割引率

rは数式上では固定ですが、実際にはデモグラ属性（年齢、地理情報、職種など）や行動（RFMなど）や在職中かどうかなどの要因により変わりうるものだと考えられます。参考文献のブログでは、このrのロジスティック回帰による推定がなされています。

データセット

データ名：CDNow
概要：1997年の第一四半期をスタート時点とした顧客の購買行動データ
期間：1997年1月〜1998年6月
顧客数：23570
取引レコード数：69659
変数：顧客ID、購入日、購入金額
入手方法：DatasetsでCDNOW dataset (full dataset)をダウンロード

関数

参考文献にはgetDataFrame関数、getPercentages関数、getCLV関数の三つの関数が出てきますが、CLVの計算に必要なのはgetDataFrame関数、getCLV関数の二つです。getPercentages関数はRecencyなどに応じて細かく分析する際に用います。

getDataFrame関数・・・生のデータセットから、指定した期間に応じたRecencyのデータを作成する関数です。

getDataFrame <- function(df,startDate,endDate,tIDColName="ID",tDateColName="Date",tAmountColName="Amount"){
  
  #日付について降順でデータフレームを並び替える
  df <- df[order(df[,tDateColName],decreasing = TRUE),]
  
  #開始期間の前と終了期間の後のデータを取り除く
  df <- df[df[,tDateColName]>= startDate,]
  df <- df[df[,tDateColName]<= endDate,]
  
  #重複したIDの行を除外して新しいデータフレームを生成する
  newdf <- df[!duplicated(df[,tIDColName]),]
  
  # 終了期間に対してのRecencyを計算する（最も小さい値は最も最近であることを示す）
  Recency <- as.numeric(difftime(endDate,newdf[,tDateColName],units="days"))
  
  # newdfデータフレームに日付の列を加える
  newdf <- cbind(newdf,Recency)
  
  # table関数やtapply関数の返す順番と揃うようにデータフレームをIDごとに並べ替える
  newdf <- newdf[order(newdf[,tIDColName]),]
  
  # 頻度を計算する
  fre <- as.data.frame(table(df[,tIDColName]))
  Frequency <- fre[,2]
  newdf <- cbind(newdf,Frequency)
  
  # 取引ごとの金額を計算する
  m <- as.data.frame(tapply(df[,tAmountColName],df[,tIDColName],sum))
  Monetary <- m[,1]/Frequency
  newdf <- cbind(newdf,Monetary)
  
  return(newdf)
  
}

getDataFrame <- function(df,startDate,endDate,tIDColName="ID",tDateColName="Date",tAmountColName="Amount"){

#日付について降順でデータフレームを並び替える

df <- df[order(df[,tDateColName],decreasing = TRUE),]

#開始期間の前と終了期間の後のデータを取り除く

df <- df[df[,tDateColName]>= startDate,]

df <- df[df[,tDateColName]<= endDate,]

#重複したIDの行を除外して新しいデータフレームを生成する

newdf <- df[!duplicated(df[,tIDColName]),]

# 終了期間に対してのRecencyを計算する（最も小さい値は最も最近であることを示す）

Recency <- as.numeric(difftime(endDate,newdf[,tDateColName],units="days"))

# newdfデータフレームに日付の列を加える

newdf <- cbind(newdf,Recency)

# table関数やtapply関数の返す順番と揃うようにデータフレームをIDごとに並べ替える

newdf <- newdf[order(newdf[,tIDColName]),]

# 頻度を計算する

fre <- as.data.frame(table(df[,tIDColName]))

Frequency <- fre[,2]

newdf <- cbind(newdf,Frequency)

# 取引ごとの金額を計算する

m <- as.data.frame(tapply(df[,tAmountColName],df[,tIDColName],sum))

Monetary <- m[,1]/Frequency

newdf <- cbind(newdf,Monetary)

return(newdf)

}

getPercentages関数・・・Recencyなどの回数に応じて、購入した顧客の割合などを計算するための関数です。

require(plyr)
getPercentages <- function(df,colNames){
  #指定した列の名前とBuyからなるベクトルの作成
  Var <- c(colNames,"Buy")
  #指定した変数とBuyが合致する列だけを抽出
  df <- df[,names(df) %in% Var,drop=F]
  
  #購買した回数に応じたデータ数を返す
  a <- ddply(df,Var,summarize,Number=length(Buy))
  #回数に応じたデータ数に対して、全体の割合を返す
  b <- ddply(a,
             .(),
             .fun=function(x){
               transform(x, Percentage=with(x,round(ave(Number,a[,names(a) %in% Var,drop=F],FUN=sum)/ave(Number,a[,names(a) %in% colNames,drop=F],FUN=sum),2)))
             })
  #1列目を除外する
  b <- b[b$Buy==1,-1]
  
  return(b)
  
}

require(plyr)

getPercentages <- function(df,colNames){

#指定した列の名前とBuyからなるベクトルの作成

Var <- c(colNames,"Buy")

#指定した変数とBuyが合致する列だけを抽出

df <- df[,names(df) %in% Var,drop=F]

#購買した回数に応じたデータ数を返す

a <- ddply(df,Var,summarize,Number=length(Buy))

#回数に応じたデータ数に対して、全体の割合を返す

b <- ddply(a,

.(),

.fun=function(x){

transform(x, Percentage=with(x,round(ave(Number,a[,names(a) %in% Var,drop=F],FUN=sum)/ave(Number,a[,names(a) %in% colNames,drop=F],FUN=sum),2)))

})

#1列目を除外する

b <- b[b$Buy==1,-1]

return(b)

}

getCLV関数・・・Recency、Frequency、Monetary、購入者の数（1人と置く）、コスト（0としている）、期間、割引率、推定したモデルをもとにCLVを計算する関数です。

getCLV<-function(r,f,m,n,cost,periods,dr,pModel){
  
  df<-data.frame(period=c(0),r=c(r),f=c(f),n=c(n),value=c(0))
  
  for(i in 1:periods){
    backstep<-df[df$period==i-1,]
    nrow<-nrow(backstep)
    for(j in 1:nrow){
      r<-backstep[j,]$r
      f<-backstep[j,]$f
      n<-backstep[j,]$n
      p<-predict(pModel,data.frame(Recency=r,Frequency=f),type='response')[1]
      buyers<-n*p
      df<-rbind(df,c(i,0,f+1,buyers,buyers*(m-cost) / (1+dr)^i))
      df<-rbind(df,c(i,r+1,f,n-buyers,(n-buyers)*(-cost)  / (1+dr)^i ))
    }
  }
  
  return(sum(df$value))
  
}

getCLV<-function(r,f,m,n,cost,periods,dr,pModel){

df<-data.frame(period=c(0),r=c(r),f=c(f),n=c(n),value=c(0))

for(i in 1:periods){

backstep<-df[df$period==i-1,]

nrow<-nrow(backstep)

for(j in 1:nrow){

r<-backstep[j,]$r

f<-backstep[j,]$f

n<-backstep[j,]$n

p<-predict(pModel,data.frame(Recency=r,Frequency=f),type='response')[1]

buyers<-n*p

df<-rbind(df,c(i,0,f+1,buyers,buyers*(m-cost) / (1+dr)^i))

df<-rbind(df,c(i,r+1,f,n-buyers,(n-buyers)*(-cost) / (1+dr)^i ))

}

return(sum(df$value))

}

データセットの読み込みと加工

#　CDNOW_SAMPLE.txtを読み込む
df <- read.table(file = "CDNOW_master.txt",header=F)

# 取引ごと且つ顧客ごとに、ID、取引日、購入金額の列からなるデータフレームを生成する
df <- as.data.frame(cbind(df[,1],df[,2],df[,4]))

# 列名の変更
names <- c("ID","Date","Amount")
names(df) <- names

#日付の型の変更
df[,2] <- as.Date(as.character(df[,2]),"%Y%m%d")

#　CDNOW_SAMPLE.txtを読み込む

df <- read.table(file = "CDNOW_master.txt",header=F)

# 取引ごと且つ顧客ごとに、ID、取引日、購入金額の列からなるデータフレームを生成する

df <- as.data.frame(cbind(df[,1],df[,2],df[,4]))

# 列名の変更

names <- c("ID","Date","Amount")

names(df) <- names

#日付の型の変更

df[,2] <- as.Date(as.character(df[,2]),"%Y%m%d")

> head(df)
  ID       Date Amount
1  1 1997-01-01  11.77
2  2 1997-01-12  12.00
3  2 1997-01-12  77.00
4  3 1997-01-02  20.76
5  3 1997-03-30  20.76
6  3 1997-04-02  19.54

> head(df)

ID Date Amount

1 1 1997-01-01 11.77

2 2 1997-01-12 12.00

3 2 1997-01-12 77.00

4 3 1997-01-02 20.76

5 3 1997-03-30 20.76

6 3 1997-04-02 19.54

再購買率とRFMとの関係

まず初めにデータセットを加工します。ロジットの推定における説明変数用のデータとして19970101〜19980228のデータを用い、被説明変数にあたる購入したかどうかのデータを19980301〜19980430のデータを用いて作ります。

# 取引期間の範囲を指定する
startDate_history <- as.Date("19970101","%Y%m%d")
endDate_history <- as.Date("19980228","%Y%m%d")

# 2ヶ月間の購入サイクルにあたる予測期間を指定する
startDate_forcast <- as.Date("19980301","%Y%m%d")
endDate_forcast <- as.Date("19980430","%Y%m%d")

# getDataFrame関数でRFMデータに変換する
history <- getDataFrame(df,startDate_history,endDate_history)
forcast <- getDataFrame(df,startDate_forcast,endDate_forcast)

# 購入サイクルを60日とし、Recencyを離散化する
history$Recency <- history$Recency %/% 60 

# Monetaryを10ドル単位で離散化する
breaks <- seq(0,round(max(history$Monetary)+9),by=10)
history$Monetary <- as.numeric(cut(history$Monetary,breaks,labels=FALSE))

# RFMのデータフレームに購入・非購入データを加える
Buy <- rep(0,nrow(history))
history <- cbind(history,Buy)

# Buyに関して、予測期間（1998/3/1~1998/4/30）の期間で再度購入しているならば1を返すようにする
history[history$ID %in% forcast$ID, ]$Buy<-1
train　<-　history

# 取引期間の範囲を指定する

startDate_history <- as.Date("19970101","%Y%m%d")

endDate_history <- as.Date("19980228","%Y%m%d")

# 2ヶ月間の購入サイクルにあたる予測期間を指定する

startDate_forcast <- as.Date("19980301","%Y%m%d")

endDate_forcast <- as.Date("19980430","%Y%m%d")

# getDataFrame関数でRFMデータに変換する

history <- getDataFrame(df,startDate_history,endDate_history)

forcast <- getDataFrame(df,startDate_forcast,endDate_forcast)

# 購入サイクルを60日とし、Recencyを離散化する

history$Recency <- history$Recency %/% 60

# Monetaryを10ドル単位で離散化する

breaks <- seq(0,round(max(history$Monetary)+9),by=10)

history$Monetary <- as.numeric(cut(history$Monetary,breaks,labels=FALSE))

# RFMのデータフレームに購入・非購入データを加える

Buy <- rep(0,nrow(history))

history <- cbind(history,Buy)

# Buyに関して、予測期間（1998/3/1~1998/4/30）の期間で再度購入しているならば1を返すようにする

history[history$ID %in% forcast$ID, ]$Buy<-1

train　<-　history

データを確認します。

> head(train)
   ID       Date Amount Recency Frequency Monetary Buy
1   1 1997-01-01  11.77       7         1        2   0
2   2 1997-01-12  12.00       6         2        5   0
8   3 1997-11-25  20.96       1         5        3   0
13  4 1997-12-12  26.48       1         4        3   0
24  5 1998-01-03  37.47       0        11        4   0
25  6 1997-01-01  20.99       7         1        3   0

> head(train)

ID Date Amount Recency Frequency Monetary Buy

1 1 1997-01-01 11.77 7 1 2 0

2 2 1997-01-12 12.00 6 2 5 0

8 3 1997-11-25 20.96 1 5 3 0

13 4 1997-12-12 26.48 1 4 3 0

24 5 1998-01-03 37.47 0 11 4 0

25 6 1997-01-01 20.99 7 1 3 0

# getPercentages関数を用いて、Recency変数に基づいて購入者の割合を計算する
colNames <- c("Recency")
p <- getPercentages(train,colNames)

# getPercentages関数を用いて、Recency変数に基づいて購入者の割合を計算する

colNames <- c("Recency")

p <- getPercentages(train,colNames)

> p
   Recency Buy Number Percentage
2        0   1   1180       0.45
4        1   1    581       0.28
6        2   1    279       0.22
8        3   1    198       0.17
10       4   1    163       0.14
12       5   1    249       0.05
14       6   1    316       0.03
16       7   1     13       0.03

> p

Recency Buy Number Percentage

2 0 1 1180 0.45

4 1 1 581 0.28

6 2 1 279 0.22

8 3 1 198 0.17

10 4 1 163 0.14

12 5 1 249 0.05

14 6 1 316 0.03

16 7 1 13 0.03

60日以内に購入した顧客（Recency=0）のうち、45%が再び購入しているようです。

# getPercentages関数を用いて、Monetary変数に基づいて購入者の割合を計算する
colNames <- c("Monetary")
p <- getPercentages(train,colNames)

# getPercentages関数を用いて、Monetary変数に基づいて購入者の割合を計算する

colNames <- c("Monetary")

p <- getPercentages(train,colNames)

> head(p,10)
   Monetary Buy Number Percentage
2         1   1     43       0.05
4         2   1    727       0.09
6         3   1    760       0.15
8         4   1    525       0.17
10        5   1    333       0.16
12        6   1    218       0.17
14        7   1    117       0.16
16        8   1     87       0.16
18        9   1     45       0.13
20       10   1     44       0.19

> head(p,10)

Monetary Buy Number Percentage

2 1 1 43 0.05

4 2 1 727 0.09

6 3 1 760 0.15

8 4 1 525 0.17

10 5 1 333 0.16

12 6 1 218 0.17

14 7 1 117 0.16

16 8 1 87 0.16

18 9 1 45 0.13

20 10 1 44 0.19

100ドル購入した顧客（Monetary=10）のうち、19%が再び購入しているようです。

# getPercentages関数を用いて、Frequency変数に基づいて購入者の割合を計算する
colNames<-c("Frequency")
p<-getPercentages(train,colNames)

# getPercentages関数を用いて、Frequency変数に基づいて購入者の割合を計算する

colNames<-c("Frequency")

p<-getPercentages(train,colNames)

> head(p,10)
   Frequency Buy Number Percentage
2          1   1    399       0.03
4          2   1    455       0.11
6          3   1    411       0.18
8          4   1    299       0.23
10         5   1    293       0.32
12         6   1    212       0.38
14         7   1    164       0.42
16         8   1    123       0.41
18         9   1    108       0.50
20        10   1     70       0.49

> head(p,10)

Frequency Buy Number Percentage

2 1 1 399 0.03

4 2 1 455 0.11

6 3 1 411 0.18

8 4 1 299 0.23

10 5 1 293 0.32

12 6 1 212 0.38

14 7 1 164 0.42

16 8 1 123 0.41

18 9 1 108 0.50

20 10 1 70 0.49

10回購入したことのある顧客（Frequency=10）は49%が再び購入しているようです。

再購買率の推定

RFM（Recency、Frequency、Monetary）のデータに基づいて、再購買率をロジスティック回帰によって推定し、予測確率を用いて顧客生涯価値を計算します。

# 再購買の割合をRecencyでロジスティック回帰
r.glm = glm(Percentage~Recency,family=quasibinomial(link='logit'),data=p)
p_r <- p

# 再購買の割合をRecencyでロジスティック回帰

r.glm = glm(Percentage~Recency,family=quasibinomial(link='logit'),data=p)

p_r <- p

> summary(r.glm)

Call:
glm(formula = Percentage ~ Recency, family = quasibinomial(link = "logit"), 
    data = p)

Deviance Residuals: 
      Min         1Q     Median         3Q        Max  
-0.085865  -0.067551   0.000333   0.041373   0.116927  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.28538    0.11395  -2.504   0.0462 *  
Recency     -0.47037    0.04035 -11.657  2.4e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.005805806)

    Null deviance: 1.052105  on 7  degrees of freedom
Residual deviance: 0.034439  on 6  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 6

> summary(r.glm)

Call:

glm(formula = Percentage ~ Recency, family = quasibinomial(link = "logit"),

data = p)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.085865 -0.067551 0.000333 0.041373 0.116927

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -0.28538 0.11395 -2.504 0.0462 *

Recency -0.47037 0.04035 -11.657 2.4e-05 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.005805806)

Null deviance: 1.052105 on 7 degrees of freedom

Residual deviance: 0.034439 on 6 degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 6

# 再購買の割合を頻度でロジスティック回帰
f.glm = glm(Percentage~Frequency,family=quasibinomial(link='logit'),data=p)
p_f <- p

# 再購買の割合を頻度でロジスティック回帰

f.glm = glm(Percentage~Frequency,family=quasibinomial(link='logit'),data=p)

p_f <- p

> summary(f.glm)

Call:
glm(formula = Percentage ~ Frequency, family = quasibinomial(link = "logit"), 
    data = p)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.79432  -0.03744   0.10447   0.19982   0.52309  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.48217    0.21490  -6.897 6.05e-09 ***
Frequency    0.12601    0.01139  11.067 1.67e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.08078336)

    Null deviance: 25.652  on 55  degrees of freedom
Residual deviance:  4.053  on 54  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 7

> summary(f.glm)

Call:

glm(formula = Percentage ~ Frequency, family = quasibinomial(link = "logit"),

data = p)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.79432 -0.03744 0.10447 0.19982 0.52309

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -1.48217 0.21490 -6.897 6.05e-09 ***

Frequency 0.12601 0.01139 11.067 1.67e-15 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.08078336)

Null deviance: 25.652 on 55 degrees of freedom

Residual deviance: 4.053 on 54 degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 7

# 再購買の割合をMonetaryでロジスティック回帰
m.glm = glm(Percentage~Monetary,family=quasibinomial(link='logit'),data=p)
p_m <- p

# 再購買の割合をMonetaryでロジスティック回帰

m.glm = glm(Percentage~Monetary,family=quasibinomial(link='logit'),data=p)

p_m <- p

> summary(m.glm)

Call:
glm(formula = Percentage ~ Monetary, family = quasibinomial(link = "logit"), 
    data = p)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.40487  -0.16209   0.00298   0.15201   0.96271  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.61381    0.27229  -9.599 2.53e-09 ***
Monetary     0.07482    0.01266   5.908 6.04e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.06566381)

    Null deviance: 4.2966  on 23  degrees of freedom
Residual deviance: 1.8027  on 22  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 5

> summary(m.glm)

Call:

glm(formula = Percentage ~ Monetary, family = quasibinomial(link = "logit"),

data = p)

Deviance Residuals:

Min 1Q Median 3Q Max

-0.40487 -0.16209 0.00298 0.15201 0.96271

Coefficients:

Estimate Std. Error t value Pr(>|t|)

(Intercept) -2.61381 0.27229 -9.599 2.53e-09 ***

Monetary 0.07482 0.01266 5.908 6.04e-06 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasibinomial family taken to be 0.06566381)

Null deviance: 4.2966 on 23 degrees of freedom

Residual deviance: 1.8027 on 22 degrees of freedom

AIC: NA

Number of Fisher Scoring iterations: 5

RecencyとFrequencyによるロジスティック回帰

model <- glm(Buy~Recency+Frequency,family=quasibinomial(link='logit'),data=train)
pred <- predict(model,data.frame(Recency=c(0),Frequency=c(1)),type='response')

1 2	model <- glm(Buy~Recency+Frequency,family=quasibinomial(link='logit'),data=train) pred <- predict(model,data.frame(Recency=c(0),Frequency=c(1)),type='response')

> pred
        1 
0.2579282

> pred

0.2579282

顧客生涯価値の計算

推定したロジットを用いて、生涯価値を計算します。

# 3期間に関する顧客生涯価値の計算（Recency=0、Frequency=1、平均収益=100ドル、割引率=2%）
v<-getCLV(0,1,100,1,0,3,0.02,model)

1 2	# 3期間に関する顧客生涯価値の計算（Recency=0、Frequency=1、平均収益=100ドル、割引率=2%） v<-getCLV(0,1,100,1,0,3,0.02,model)

> v
[1] 63.91906

1 2	> v [1] 63.91906

それではさっそく、1998年5月〜6月のデータを用いて、今回推定した顧客生涯価値が妥当なのかどうかを確かめたいと思います。

#トレーニングデータをもとにCLVを計算する
train_pred <- cbind(train,lifetimevalue=matrix(,nrow(train),1))

pb <- txtProgressBar(min = 1, max = nrow(train_pred), style = 3)

for(i in 1:nrow(train_pred)){
  train_pred$lifetimevalue[i] <- getCLV(train_pred$Recency[i],train_pred$Frequency[i],100,1,0,3,0.02,model)
  setTxtProgressBar(pb, i) 
}

#1998年5月1日〜6月30日の間の取引データからIDと購入金額を集計する。
startDate_test <- as.Date("19980501","%Y%m%d")
endDate_test <- as.Date("19980630","%Y%m%d")
test <- getDataFrame(df,startDate_test,endDate_test)
test_monetary <- dplyr::select(test,ID,Monetary)

#予測したCLVと実際の取引金額のデータをマージする。
train_pred_2 <- data.frame(train_pred[,colnames(train_pred) %in% c("ID","lifetimevalue")])
train_pred_2 <- merge(train_pred,test_monetary,by = "ID",all=T)
train_pred_2$Monetary.y[is.na(train_pred_2$Monetary.y)] <- 0

#トレーニングデータをもとにCLVを計算する

train_pred <- cbind(train,lifetimevalue=matrix(,nrow(train),1))

pb <- txtProgressBar(min = 1, max = nrow(train_pred), style = 3)

for(i in 1:nrow(train_pred)){

train_pred$lifetimevalue[i] <- getCLV(train_pred$Recency[i],train_pred$Frequency[i],100,1,0,3,0.02,model)

setTxtProgressBar(pb, i)

}

#1998年5月1日〜6月30日の間の取引データからIDと購入金額を集計する。

startDate_test <- as.Date("19980501","%Y%m%d")

endDate_test <- as.Date("19980630","%Y%m%d")

test <- getDataFrame(df,startDate_test,endDate_test)

test_monetary <- dplyr::select(test,ID,Monetary)

#予測したCLVと実際の取引金額のデータをマージする。

train_pred_2 <- data.frame(train_pred[,colnames(train_pred) %in% c("ID","lifetimevalue")])

train_pred_2 <- merge(train_pred,test_monetary,by = "ID",all=T)

train_pred_2$Monetary.y[is.na(train_pred_2$Monetary.y)] <- 0

予測したCLVと実際の取引金額データを散布図で描き、回帰線を引く。

library(ggplot2)
g <- ggplot(train_pred_2,aes (x = lifetimevalue,y = Monetary.y))
g <- g +  geom_point(shape = 20,size = 0.8,na.rm = TRUE)
g <- g + geom_smooth(method = glm) + xlab("lifetimevalue") + ylab("Monetary")
plot(g)

library(ggplot2)

g <- ggplot(train_pred_2,aes (x = lifetimevalue,y = Monetary.y))

g <- g + geom_point(shape = 20,size = 0.8,na.rm = TRUE)

g <- g + geom_smooth(method = glm) + xlab("lifetimevalue") + ylab("Monetary")

plot(g)

CLVが上がれば、1998年5月1日〜6月30日の間（未来）の取引金額が増すような傾向が出ています。

参考情報

RFM Customer Analysis with R Language
Calculating Customer Lifetime Value with Recency, Frequency, and Monetary (RFM)

Japan.R 2016のスライドまとめ

まだ手に入れていないスライドもあるので随時更新しますが、Japan.R 2016(connpass)のスライドをまとめています。後日、登場したパッケージなどのサンプルコードも載せていく予定です。

目次

・石田基広さんのキーノート
・ホクソエムとは何だったのか（ホクソエムさん）
・Rと探索的データ分析で、国連での日本の立ち位置を可視化する（安田洋介さん）
・マウス操作でかんたん予測分析（鈴木了太さん）
・高速・省メモリにlibsvm形式でダンプする方法を研究してみた（@hskksk）
・Rでてんしょくかつどう（@Med_KU）
・RStudio vs Emacs（@y__mattu）
・randomforestで高次元の変数重要度見る（@siero5335）
・Rで本を作りたい（前田和寛さん）
・28歳でプログラミングを始めた話（市川太祐さん）
・LDA-Visパッケージのご紹介（@doradora09）
・【e2d3R】E2D3からDot-Bar-Chartのご紹介（楠本一哲さん）
・このIRのグラフがすごい！上場企業2016（@ito_yan）
・Rでカルマンフィルタをしたい（@tetsuroito）
・PPAP（仮）（@yutannihilation）
・スライド未公開、ユーザーの状態遷移に関する分析のお話（@sanoche16）
・私とR（高栁慎一さん）
・めくってもめくってもサンプル画像（服部恵美さん）
・木と電話と選挙(causalTree)（安井翔太さん）
・スライド未公開、dplyrの話（@tomomoto）
・てかLINEやってる？(仮)（@wonder_zone）
・心理学における「再現性」の問題とBayes Factor（@NSushi）

・石田基広さんのキーノート

スライド未公開です。

・Linux使い
・ヘブライ語の意味構造を代数学でやっていた
・S/R言語の生みの親はJohn Chambers
　以下の二つは最近書かれた本だそうです。
　Software for Data Analysis: Programming with R (Statistics and Computing)
　Extending R (Chapman & Hall/CRC The R Series)
・S→S-plus→Rの順番で発展
・purrrを最近使い始めたそうです。
・XLConnectパッケージを使って、大学教員の採点活動を効率化しているそうです。

・ホクソエムとは何だったのか（ホクソエムさん）

匿名技術者集団ホクソエムの2016年の成果
・densratio（ densratio: Density Ratio Estimation ）
・githubinstall
（ githubinstall: A Helpful Way to Install R Packages Hosted on GitHub ）
・healthplanet（ Wrapper package for healthplanet api ）
・RODBCDBI
（ RODBCDBI: Provides Access to Databases Through the ODBC Interface ）
・jpmesh（ jpmesh: Utilities for Japanese Mesh Code ）

起業されたとのことです。懸命に頑張って下さい！
株式会社ホクソエム

awesomeな人材が必要とのことで、awesomeな方はアプライしてみてはいかがでしょうか。

・Rと探索的データ分析で、国連での日本の立ち位置を可視化する（安田洋介さん）

スライド未公開です。
国連のデータを使って、Exploratoryを用いた探索的データ分析の実演をされていました。

・マウス操作でかんたん予測分析（鈴木了太さん）

R AnalyticFlow
Rで実践！データサイエンス～初めの一歩から高度な応用まで～

・高速・省メモリにlibsvm形式でダンプする方法を研究してみた（@hskksk）

高速・省メモリにlibsvm形式でダンプする方法を研究してみた from Keisuke Hosaka

・Rでてんしょくかつどう（@Med_KU）

Rmd でreveal.js のhtml スライドプレゼンテーション

・RStudio vs Emacs（@y__mattu）

RStudio vs Emacs Japan.R 2016

・randomforestで高次元の変数重要度見る（@siero5335）

Randomforestで高次元の変数重要度を見る #japanr LT from Akifumi Eguchi

・Rで本を作りたい（前田和寛さん）

Rで本を作りたい

・28歳でプログラミングを始めた話（市川太祐さん）

・医療関連のアプリ開発でデータサイエンスを駆使しようとしているそうです。

スライド未公開です。
スライドがシェアされ次第載せます。

・LDA-Visパッケージのご紹介（@doradora09）

20161127 doradora09 japanr2016_lt from Nobuaki Oshiro

・【e2d3R】E2D3からDot-Bar-Chartのご紹介（楠本一哲さん）

スライドは未公開です。
E2D3をRで表示する試みのようです。
Experiments with e2d3 in R

・このIRのグラフがすごい！上場企業2016（@ito_yan）

スライド未公開です。後日シェアしていただけるようです。

今年のプレゼンはネタを集めるにしても、これが限界でした…楽しんで頂けたでしょうか？スライドは後日アップします。 #JapanR

— いとーちゃん (@ito_yan) 2016年11月27日

・Rでカルマンフィルタをしたい（@tetsuroito）

・PPAP（仮）（@yutannihilation）

Tidyverseとは from yutannihilation

・スライド未公開、ユーザーの状態遷移に関する分析のお話（@sanoche16）

スライドがシェアされ次第載せます。

・私とR（高栁慎一さん）

RjpWiki
統計・データ解析
 統計解析フリーソフト R の備忘録頁 ver.3.1
seekR(R限定の検索エンジン)
からだにいいもの
アブラタニブログってなんでしょう。油谷さんのブログ？

・めくってもめくってもサンプル画像（服部恵美さん）

Rのサンプルコードはあるけれども、どんな図ができるのかはわからない。そこで、サンプルコードとグラフを大量にまとめているサイトを作ったそうです。検索性は未知数ですが、暇なときに眺めておきたいですね。
R Graphical Manual

・木と電話と選挙(causalTree)（安井翔太さん）

木と電話と選挙(causalTree) from Shota Yasui

・スライド未公開、dplyrの話（@tomomoto）

スライドがシェアされ次第載せます。

・てかLINEやってる？(仮)（@wonder_zone）

てかLINEやってる？ (Japan.R 2016 LT) #JapanR from cancolle

・心理学における「再現性」の問題とBayes Factor（@NSushi）

スライドは後日公開とのことです。

今日はjapan.Rお疲れ様でしたー！楽しかったですありがとうございましたー！
ネット環境がないので心理学の再現性およびベイズファクターおよびにこにースライドは明日アップしますー。

— numba shushi (@NSushi) 2016年11月27日

『マーケティング・サイエンス入門』に出てくる手法をRで実行してみる

友人に『マーケティング・サイエンス入門』がおすすめと言われて読んだんですが、やっぱり実行できないとモヤモヤしてしまいますよね。そこで、登場する手法に関連したRのコードやらを集めてみました。

・BASSモデル
・多次元尺度法
・因子分析
・ロジット&プロビット
・分散分析
・クラスター分析
・判別分析
・決定木
・コンジョイント分析
・RFM分析
・共分散構造分析

BASSモデル

市場全体の規模が動的にどのように変化するかを予測するために使われるモデル。
R を使ってバスモデルを当てはめてみた – 廿TT
こちらにRのコードや適用例がいくつか載っています。

早速、私も携帯電話の加入契約数の時系列データを用いて、コードを実行してみました。データは平成25年版の総務省の情報通信白書の表から得ました。（第2部　情報通信の現況・政策の動向）

> fit_logis$value
[1] 69833.02
> fit_bass$value
[1] 68588.92

> fit_logis$value

[1] 69833.02

> fit_bass$value

[1] 68588.92

当てはまりはわずかながら、BASSモデルの方が良いようです。

多次元尺度法

多次元尺度法で遊んでみる（オレ流 R入門）
こちらのブログで山手線の駅間の距離データの可視化がなされています。
各駅ごとの距離からなる行列さえ用意すれば、cmdscale()関数を実行することで可能なようです。

今回はContaminatedMixtパッケージに含まれているワインのデータセットを使って多次元尺度法を適用してみようと思います。

データはこんな感じです。

以下のコードで実行しました。

library(ContaminatedMixt)

data(wine)
wine.dist <- dist(wine[,-1])
wine.cmd  <- cmdscale(wine.dist)

plot(wine.cmd,type = "n")
wine.lab <- factor(wine$Type)
text(wine.cmd,labels = wine.lab,col = unclass(wine.lab))

library(ContaminatedMixt)

data(wine)

wine.dist <- dist(wine[,-1])

wine.cmd <- cmdscale(wine.dist)

plot(wine.cmd,type = "n")

wine.lab <- factor(wine$Type)

text(wine.cmd,labels = wine.lab,col = unclass(wine.lab))

Barbera（バルベーラ）・・・基本的にはタンニンをあまり含まず、酸味の強い色の濃い赤ワインで庶民的。
Barolo（バローロ）・・・アルコール度数が高く、非常に重厚な味わいのワインでワインの王様と呼ばれる。
Grignolino（グリニョリーノ）・・・僅かにタンニンを感じるサッパリとした辛口の赤ワインで庶民的。

庶民と王様のワインは成分においても違いがありそうですね。

因子分析

psychパッケージというものがあるようです。こちらのサイトを参考にして進めます。（スナック菓子の食感についてRで因子分析してみた）
今回は大好きなwiskyのデータセットを使ってみます。（ Classification of whiskies ）

library(psych)
library(GPArotation)

whiskies <- read.csv("whiskies.txt", row.names = 1, stringsAsFactors = FALSE)
rownames(whiskies) <- whiskies[,1]
whiskies <- whiskies[,-1]
whiskies <- whiskies[,-15]
whiskies <- whiskies[,-14]
whiskies <- whiskies[,-13]

MAPminres <- vss(whiskies, fm="minres")
print(MAPminres)

MAPml <- vss(whiskies, fm="ml")
print(MAPml)

parallel <- fa.parallel(whiskies)
print(parallel)

res <- fa(whiskies, nfactors=2, fm="minres", rotate="oblimin")
print(res, digits=3)

biplot(res, labels=rownames(whiskies))

library(psych)

library(GPArotation)

whiskies <- read.csv("whiskies.txt", row.names = 1, stringsAsFactors = FALSE)

rownames(whiskies) <- whiskies[,1]

whiskies <- whiskies[,-1]

whiskies <- whiskies[,-15]

whiskies <- whiskies[,-14]

whiskies <- whiskies[,-13]

MAPminres <- vss(whiskies, fm="minres")

print(MAPminres)

MAPml <- vss(whiskies, fm="ml")

print(MAPml)

parallel <- fa.parallel(whiskies)

print(parallel)

res <- fa(whiskies, nfactors=2, fm="minres", rotate="oblimin")

print(res, digits=3)

biplot(res, labels=rownames(whiskies))

グレンフィディックやカリラやタリスカーがイメージ通りにプロットされています。ラガブーリンやアードベッグがはみ出しているのが残念ですが。

ロジット・プロビット

これらの手法はビルトインの関数でできてしまいますが、せっかくウイスキーのデータがあるので、薬っぽさに繋がりそうな変数を見つけてみます。

library(dplyr)

whiskies2 <- whiskies %>% mutate(target=ifelse(whiskies$Medicinal>0,1,0))
whiskies2 <- whiskies2[,-which(colnames(whiskies2)=="Medicinal")]

#ロジット
model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))
summary(model_logit)

#プロビット
model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))
summary(model_probit)

library(dplyr)

whiskies2 <- whiskies %>% mutate(target=ifelse(whiskies$Medicinal>0,1,0))

whiskies2 <- whiskies2[,-which(colnames(whiskies2)=="Medicinal")]

#ロジット

model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))

summary(model_logit)

#プロビット

model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))

summary(model_probit)

推定結果はこちらです。スモーキーさが関係しているのは納得です。

> model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))
> summary(model_logit)

Call:
glm(formula = target ~ ., family = binomial(link = "logit"), 
    data = whiskies2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0173  -0.6313  -0.3509   0.2206   2.4214  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  4.39412    2.47716   1.774   0.0761 .
Body        -0.06226    0.51985  -0.120   0.9047  
Sweetness   -0.59932    0.46458  -1.290   0.1970  
Smoky        1.35983    0.58327   2.331   0.0197 *
Tobacco      1.75440    1.10595   1.586   0.1127  
Honey       -0.39101    0.46604  -0.839   0.4015  
Spicy       -0.22916    0.44331  -0.517   0.6052  
Winey       -0.86015    0.49829  -1.726   0.0843 .
Nutty       -0.56436    0.42843  -1.317   0.1878  
Malty       -1.17076    0.66243  -1.767   0.0772 .
Fruity      -0.33978    0.43754  -0.777   0.4374  
Floral      -0.64322    0.47291  -1.360   0.1738  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.023  on 85  degrees of freedom
Residual deviance:  64.882  on 74  degrees of freedom
AIC: 88.882

Number of Fisher Scoring iterations: 6

> model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))
> summary(model_probit)

Call:
glm(formula = target ~ ., family = binomial(link = "probit"), 
    data = whiskies2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9365  -0.6399  -0.3491   0.2064   2.4188  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  2.53606    1.41370   1.794   0.0728 .
Body        -0.05535    0.28834  -0.192   0.8478  
Sweetness   -0.36200    0.26740  -1.354   0.1758  
Smoky        0.76890    0.32649   2.355   0.0185 *
Tobacco      1.01419    0.64201   1.580   0.1142  
Honey       -0.19843    0.26056  -0.762   0.4463  
Spicy       -0.13537    0.25214  -0.537   0.5913  
Winey       -0.50543    0.27073  -1.867   0.0619 .
Nutty       -0.30335    0.24152  -1.256   0.2091  
Malty       -0.65838    0.36494  -1.804   0.0712 .
Fruity      -0.20000    0.25160  -0.795   0.4267  
Floral      -0.35983    0.26557  -1.355   0.1754  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.02  on 85  degrees of freedom
Residual deviance:  65.06  on 74  degrees of freedom
AIC: 89.06

Number of Fisher Scoring iterations: 7

> model_logit <- glm(target ~ ., data=whiskies2, family=binomial(link="logit"))

> summary(model_logit)

Call:

glm(formula = target ~ ., family = binomial(link = "logit"),

data = whiskies2)

Deviance Residuals:

Min 1Q Median 3Q Max

-2.0173 -0.6313 -0.3509 0.2206 2.4214

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 4.39412 2.47716 1.774 0.0761 .

Body -0.06226 0.51985 -0.120 0.9047

Sweetness -0.59932 0.46458 -1.290 0.1970

Smoky 1.35983 0.58327 2.331 0.0197 *

Tobacco 1.75440 1.10595 1.586 0.1127

Honey -0.39101 0.46604 -0.839 0.4015

Spicy -0.22916 0.44331 -0.517 0.6052

Winey -0.86015 0.49829 -1.726 0.0843 .

Nutty -0.56436 0.42843 -1.317 0.1878

Malty -1.17076 0.66243 -1.767 0.0772 .

Fruity -0.33978 0.43754 -0.777 0.4374

Floral -0.64322 0.47291 -1.360 0.1738

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 107.023 on 85 degrees of freedom

Residual deviance: 64.882 on 74 degrees of freedom

AIC: 88.882

Number of Fisher Scoring iterations: 6

> model_probit <-glm(target ~ ., data=whiskies2, family=binomial(link="probit"))

> summary(model_probit)

Call:

glm(formula = target ~ ., family = binomial(link = "probit"),

data = whiskies2)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.9365 -0.6399 -0.3491 0.2064 2.4188

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) 2.53606 1.41370 1.794 0.0728 .

Body -0.05535 0.28834 -0.192 0.8478

Sweetness -0.36200 0.26740 -1.354 0.1758

Smoky 0.76890 0.32649 2.355 0.0185 *

Tobacco 1.01419 0.64201 1.580 0.1142

Honey -0.19843 0.26056 -0.762 0.4463

Spicy -0.13537 0.25214 -0.537 0.5913

Winey -0.50543 0.27073 -1.867 0.0619 .

Nutty -0.30335 0.24152 -1.256 0.2091

Malty -0.65838 0.36494 -1.804 0.0712 .

Fruity -0.20000 0.25160 -0.795 0.4267

Floral -0.35983 0.26557 -1.355 0.1754

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 107.02 on 85 degrees of freedom

Residual deviance: 65.06 on 74 degrees of freedom

AIC: 89.06

Number of Fisher Scoring iterations: 7

ちなみに、多項ロジットに関しては、mlogitパッケージを使えばできるようです。（多項ロジット（Multinomial Logit）, R – mlogit 使用メモ）大学院時代に多項ロジットはSTATAでよく使っていましたが、Rだとこのパッケージなんですかね。推定した係数の値の解釈が若干複雑だったりします。

分散分析

分散分析もビルトインの関数で実行することができます。今回はワインのデータを用いて、銘柄から30個ランダムサンプリングをした上で、アルコールに関して群間の母平均値が同じかどうかを確かめてみます。コードはこちらを参考にしました。（ R による分散分析（一元配置））

A <- wine %>% filter(Type=="Barbera")
sampleNum_A <- sample(nrow(A),30)
A <- A[sampleNum_A,]

B <- wine %>% filter(Type=="Barolo")
sampleNum_B <- sample(nrow(B),30)
B <- B[sampleNum_B,]

C <- wine %>% filter(Type=="Grignolino")
sampleNum_C <- sample(nrow(C),30)
C <- C[sampleNum_C,]

d <- data.frame(A = A$Alcohol, B = B$Alcohol, C = C$Alcohol )
library(reshape)
x <- melt(d, variable_name = "group")
res <- aov(value ~ factor(group), data = x)
summary(res)

A <- wine %>% filter(Type=="Barbera")

sampleNum_A <- sample(nrow(A),30)

A <- A[sampleNum_A,]

B <- wine %>% filter(Type=="Barolo")

sampleNum_B <- sample(nrow(B),30)

B <- B[sampleNum_B,]

C <- wine %>% filter(Type=="Grignolino")

sampleNum_C <- sample(nrow(C),30)

C <- C[sampleNum_C,]

d <- data.frame(A = A$Alcohol, B = B$Alcohol, C = C$Alcohol )

library(reshape)

x <- melt(d, variable_name = "group")

res <- aov(value ~ factor(group), data = x)

summary(res)

推定結果はこちらです。アルコールに関しては、3群間において差があるようです。

> summary(res)
              Df Sum Sq Mean Sq F value Pr(>F)    
factor(group)  2  38.12  19.059   91.22 <2e-16 ***
Residuals     87  18.18   0.209                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(res)

Df Sum Sq Mean Sq F value Pr(>F)

factor(group) 2 38.12 19.059 91.22 <2e-16 ***

Residuals 87 18.18 0.209

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

クラスター分析

クラスター分析もビルトインの関数で実行可能です。ここでは参考文献（ K-means Clustering 86 Single Malt Scotch Whiskies ）のウイスキーのサンプルで取り上げられたK-mean法をそのまま紹介します。

whiskies_k <- scale(whiskies)
ssPlot <- function(data, maxCluster = 9) {
  # Initialize within sum of squares
  SSw <- (nrow(data) - 1) * sum(apply(data, 2, var))
  SSw <- vector()
  for (i in 2:maxCluster) {
    SSw[i] <- sum(kmeans(data, centers = i)$withinss)
  }
  plot(1:maxCluster, SSw, type = "b", xlab = "Number of Clusters", ylab = "Within groups sum of squares")
}
ssPlot(whiskies_k)

fit <- kmeans(whiskies_k, 4) 
whiskies <- data.frame(whiskies, fit$cluster)
whiskies$fit.cluster <- as.factor(whiskies$fit.cluster)

subset(whiskies, fit.cluster == 2)

whiskies_k <- scale(whiskies)

ssPlot <- function(data, maxCluster = 9) {

# Initialize within sum of squares

SSw <- (nrow(data) - 1) * sum(apply(data, 2, var))

SSw <- vector()

for (i in 2:maxCluster) {

SSw[i] <- sum(kmeans(data, centers = i)$withinss)

}

plot(1:maxCluster, SSw, type = "b", xlab = "Number of Clusters", ylab = "Within groups sum of squares")

}

ssPlot(whiskies_k)

fit <- kmeans(whiskies_k, 4)

whiskies <- data.frame(whiskies, fit$cluster)

whiskies$fit.cluster <- as.factor(whiskies$fit.cluster)

subset(whiskies, fit.cluster == 2)

気になるクラスターの結果ですが、どうやらアイラ島系のウイスキーのクラスターを作れたようです。

> subset(whiskies, fit.cluster == 2)
          Body Sweetness Smoky Medicinal Tobacco Honey Spicy Winey Nutty Malty Fruity
Ardbeg       4         1     4         4       0     0     2     0     1     2      1
Caol Ila     3         1     4         2       1     0     2     0     2     1      1
Clynelish    3         2     3         3       1     0     2     0     1     1      2
Lagavulin    4         1     4         4       1     0     1     2     1     1      1
Laphroig     4         2     4         4       1     0     0     1     1     1      0
Talisker     4         2     3         3       0     1     3     0     1     2      2
          Floral fit.cluster
Ardbeg         0           2
Caol Ila       1           2
Clynelish      0           2
Lagavulin      0           2
Laphroig       0           2
Talisker       0           2

> subset(whiskies, fit.cluster == 2)

Body Sweetness Smoky Medicinal Tobacco Honey Spicy Winey Nutty Malty Fruity

Ardbeg 4 1 4 4 0 0 2 0 1 2 1

Caol Ila 3 1 4 2 1 0 2 0 2 1 1

Clynelish 3 2 3 3 1 0 2 0 1 1 2

Lagavulin 4 1 4 4 1 0 1 2 1 1 1

Laphroig 4 2 4 4 1 0 0 1 1 1 0

Talisker 4 2 3 3 0 1 3 0 1 2 2

Floral fit.cluster

Ardbeg 0 2

Caol Ila 1 2

Clynelish 0 2

Lagavulin 0 2

Laphroig 0 2

Talisker 0 2

判別分析

MASSパッケージで実行可能です。線形識別関数の実行例がこちらの参考文献に載っていたので、ワインのデータで試してみます。（【Rによるデータサイエンス】線形判別分析）

library(MASS)
wine2 <- data.frame(wine[,-1],type=wine.lab)
even.n <- 2*(1:(nrow(wine2)/2))-1
wine2.train <- wine2[even.n,]
wine2.test  <- wine2[-even.n,]

#train
Z.lda<-lda(type~.,data=wine2.train)
Z.lda$means
Z.lda$scaling
Z.lda$means%*%Z.lda$scaling
apply(Z.lda$means%*%Z.lda$scaling,2,mean)
table(wine2.train[,ncol(wine2)],predict(Z.lda)$class)
plot(Z.lda,dimen=1)
plot(Z.lda,dimen=2)

#test
testResult <- predict(Z.lda,wine2.test[,-ncol(wine2)])

table(wine2.test[,ncol(wine2)],testResult$class)
testResult$class

plot(testResult$x,type="n")
text(testResult$x,labels=wine2.test$type)

library(MASS)

wine2 <- data.frame(wine[,-1],type=wine.lab)

even.n <- 2*(1:(nrow(wine2)/2))-1

wine2.train <- wine2[even.n,]

wine2.test <- wine2[-even.n,]

#train

Z.lda<-lda(type~.,data=wine2.train)

Z.lda$means

Z.lda$scaling

Z.lda$means%*%Z.lda$scaling

apply(Z.lda$means%*%Z.lda$scaling,2,mean)

table(wine2.train[,ncol(wine2)],predict(Z.lda)$class)

plot(Z.lda,dimen=1)

plot(Z.lda,dimen=2)

#test

testResult <- predict(Z.lda,wine2.test[,-ncol(wine2)])

table(wine2.test[,ncol(wine2)],testResult$class)

testResult$class

plot(testResult$x,type="n")

text(testResult$x,labels=wine2.test$type)

判別関数得点

さすがパッケージ用のデータセットだけあって、綺麗に分類できたようです。誤分類は2件だけです。

> table(wine2.test[,ncol(wine2)],testResult$class)
            
             Barbera Barolo Grignolino
  Barbera         24      0          0
  Barolo           0     29          0
  Grignolino       0      2         34

> table(wine2.test[,ncol(wine2)],testResult$class)

Barbera Barolo Grignolino

Barbera 24 0 0

Barolo 0 29 0

Grignolino 0 2 34

決定木

決定木はrpartパッケージで実行します。ウイスキーのデータを使って、薬っぽさを決める条件を探してみます。コードはこちらを参考にしました。（ R言語で決定木分析）

library(rpart)
library(partykit)
library(epitools)

result <- rpart(target~., data = whiskies2)
plot(as.party(result))

library(rpart)

library(partykit)

library(epitools)

result <- rpart(target~., data = whiskies2)

plot(as.party(result))

コンジョイント分析

conjointパッケージなるものがあるようです。こちらの参考文献を元に紹介します。（ Rでコンジョイント分析）

まずは直交表を作ってみます。

library(conjoint)

experiment <- expand.grid(
  麺の量   =c("マシ","普通","少なめ"),
  ニンニク =c("マシ","普通","少なめ"),
  野菜     =c("マシ","普通","少なめ"),
  あぶら   =c("マシ","普通","少なめ")
)

design.ort <- caFactorialDesign(data = experiment,type = "orthogonal")
design.ort

caEncodedDesign(design.ort)
cor(caEncodedDesign(design.ort))

library(conjoint)

experiment <- expand.grid(

麺の量 =c("マシ","普通","少なめ"),

ニンニク =c("マシ","普通","少なめ"),

野菜 =c("マシ","普通","少なめ"),

あぶら =c("マシ","普通","少なめ")

)

design.ort <- caFactorialDesign(data = experiment,type = "orthogonal")

design.ort

caEncodedDesign(design.ort)

cor(caEncodedDesign(design.ort))

> design.ort
   麺の量 ニンニク   野菜 あぶら
5    普通     普通   マシ   マシ
10   マシ     マシ   普通   マシ
27 少なめ   少なめ 少なめ   マシ
34   マシ   少なめ   マシ   普通
42 少なめ     普通   普通   普通
47   普通     マシ 少なめ   普通
57 少なめ     マシ   マシ 少なめ
71   普通   少なめ   普通 少なめ
76   マシ     普通 少なめ 少なめ

> caEncodedDesign(design.ort)
   麺の量 ニンニク 野菜 あぶら
5       2        2    1      1
10      1        1    2      1
27      3        3    3      1
34      1        3    1      2
42      3        2    2      2
47      2        1    3      2
57      3        1    1      3
71      2        3    2      3
76      1        2    3      3
> cor(caEncodedDesign(design.ort))
         麺の量 ニンニク 野菜 あぶら
麺の量        1        0    0      0
ニンニク      0        1    0      0
野菜          0        0    1      0
あぶら        0        0    0      1

> design.ort

麺の量ニンニク野菜あぶら

5 普通普通マシマシ

10 マシマシ普通マシ

27 少なめ少なめ少なめマシ

34 マシ少なめマシ普通

42 少なめ普通普通普通

47 普通マシ少なめ普通

57 少なめマシマシ少なめ

71 普通少なめ普通少なめ

76 マシ普通少なめ少なめ

> caEncodedDesign(design.ort)

麺の量ニンニク野菜あぶら

5 2 2 1 1

10 1 1 2 1

27 3 3 3 1

34 1 3 1 2

42 3 2 2 2

47 2 1 3 2

57 3 1 1 3

71 2 3 2 3

76 1 2 3 3

> cor(caEncodedDesign(design.ort))

麺の量ニンニク野菜あぶら

麺の量 1 0 0 0

ニンニク 0 1 0 0

野菜 0 0 1 0

あぶら 0 0 0 1

残念ながら、面白そうなデータがないので、サンプルについているお茶のデータを使ってみます。

library(conjoint)

data(tea)

print(tprof)
print(tlevn)
print(tprefm)

Conjoint(tprefm,tprof,tlevn)

library(conjoint)

data(tea)

print(tprof)

print(tlevn)

print(tprefm)

Conjoint(tprefm,tprof,tlevn)

価格が効用に与える影響の可視化です。

RFM分析

ほくそ笑むの親分がeasyRFMパッケージを作っていたようです。（ RFM 分析を簡単に実行できる R パッケージ easyRFM を作った）都合良く取引データがなかったので、kaggleの掲示板で落ちていたデータを使いました。（ Sample of transaction data ）

library(easyRFM)

transaction <- read.csv("transactions-sample.csv")
transaction <- transaction %>% mutate(payment=purchasequantity*purchaseamount)
transaction <- transaction[,which(colnames(transaction) %in% c("id","payment","date"))]
transaction$date <- as.character(transaction$date)

result <- rfm_auto(transaction, breaks=3)
head(result$rfm)
result$classes

library(easyRFM)

transaction <- read.csv("transactions-sample.csv")

transaction <- transaction %>% mutate(payment=purchasequantity*purchaseamount)

transaction <- transaction[,which(colnames(transaction) %in% c("id","payment","date"))]

transaction$date <- as.character(transaction$date)

result <- rfm_auto(transaction, breaks=3)

head(result$rfm)

result$classes

結果は以下のとおりです。

> head(result$rfm)
        id    Recency Frequency  Monetary RecencyClass FrequencyClass MonetaryClass
1    86246 2013-04-23     12609  89333.08            1              3             3
2    86252 2013-03-26     12087 104436.35            1              3             3
3 12262064 2013-06-25      1101   5116.08            3              1             1
4 12277270 2013-06-19      1361  10820.49            3              2             2
5 12332190 2013-06-13       684   8709.90            2              1             1
> result$classes
$recency_class
[1] "2013-03-26 00:00:00 to 2013-05-11" "2013-05-11 00:00:01 to 2013-06-18"
[3] "2013-06-18 00:00:01 to 2013-06-26"

$recency_class_days
[1] "92 to 46" "45 to 8"  "7 to 0"  

$frequency_class
[1] "681 to 1200"   "1201 to 8600"  "8601 to 13000"

$monetary_class
[1] "5100 to 9500"    "9501 to 64000"   "64001 to 110000"

> head(result$rfm)

id Recency Frequency Monetary RecencyClass FrequencyClass MonetaryClass

1 86246 2013-04-23 12609 89333.08 1 3 3

2 86252 2013-03-26 12087 104436.35 1 3 3

3 12262064 2013-06-25 1101 5116.08 3 1 1

4 12277270 2013-06-19 1361 10820.49 3 2 2

5 12332190 2013-06-13 684 8709.90 2 1 1

> result$classes

$recency_class

[1] "2013-03-26 00:00:00 to 2013-05-11" "2013-05-11 00:00:01 to 2013-06-18"

[3] "2013-06-18 00:00:01 to 2013-06-26"

$recency_class_days

[1] "92 to 46" "45 to 8" "7 to 0"

$frequency_class

[1] "681 to 1200" "1201 to 8600" "8601 to 13000"

$monetary_class

[1] "5100 to 9500" "9501 to 64000" "64001 to 110000"

共分散構造分析

semパッケージで実行可能です。こちらの参考文献のデータを用います。（ Rによるパス解析）データはこちらにあります。（練習用データ）

library(sem)

data01	<-	read.csv(file="regpath.csv",header=T,
                   fileEncoding="Shift-JIS")

model01	<-	specifyModel()
意欲	->	学力,	b1,	NA
意欲	->	適応,	b2,	NA
適応	->	学力,	b3, NA
学力	->	自尊,	b4,	NA
学力	<->	学力,	v1,	NA
自尊	<->	自尊,	v2,	NA
適応	<->	適応,	v3, NA
意欲	<->	意欲,	v4,	NA

cov01	<-	cov(data01[,2:5])
cov01

fit01	<-	sem(model=model01,	S=cov01,	N=nrow(data01))

summary(fit01,rsquare=T,fit.indices=c("GFI","AGFI","SRMR","RMSEA","AIC","BIC"))
stdCoef(fit01)

library(sem)

data01 <- read.csv(file="regpath.csv",header=T,

fileEncoding="Shift-JIS")

model01 <- specifyModel()

意欲 -> 学力, b1, NA

意欲 -> 適応, b2, NA

適応 -> 学力, b3, NA

学力 -> 自尊, b4, NA

学力 <-> 学力, v1, NA

自尊 <-> 自尊, v2, NA

適応 <-> 適応, v3, NA

意欲 <-> 意欲, v4, NA

cov01 <- cov(data01[,2:5])

cov01

fit01 <- sem(model=model01, S=cov01, N=nrow(data01))

summary(fit01,rsquare=T,fit.indices=c("GFI","AGFI","SRMR","RMSEA","AIC","BIC"))

stdCoef(fit01)

Rを使った分析（SEM）
こちらの方がパスの図も出力できるので、良いかもしれません。

参考文献

マーケティング・サイエンス入門 — 市場対応の科学的マネジメント新版 (有斐閣アルマ)

Cloud Vision APIをRで動かしてみる

Google Cloud Vision APIを使ってAKBのデータを把握する
こちらの所沢義男さん（偽名）のブログにカジュアルにVison APIを使う方法が載せられていたので、早速私も触ってみようと思います。（Vision APIのアカウント作成には住所や氏名や電話番号、クレカ情報が求められます。）

ラーメン二郎のメニュー表

紹介されていたコードを実行してみます。

CROWD_VISION_KEY <- "YOUR KEY"

library("httr")

f <- "http://livedoor.blogimg.jp/musashinonikki/imgs/9/b/9b62aa78.jpg"
tmp <- tempfile()
download.file(f, tmp)
img <- readBin(tmp, "raw", file.info(tmp)[1, "size"])
u <- paste0("https://vision.googleapis.com/v1/images:annotate?key=", CROWD_VISION_KEY)
body <- list(requests = list(image=list(content=img), 
                             features=list(type="TEXT_DETECTION")
                             )
             )
res <- POST(url=u,
            encode="json",
            body=body,
            content_type_json()
)
res <- content(res)
cat(res$responses[[1]]$textAnnotations[[1]]$description[1])

CROWD_VISION_KEY <- "YOUR KEY"

library("httr")

f <- "http://livedoor.blogimg.jp/musashinonikki/imgs/9/b/9b62aa78.jpg"

tmp <- tempfile()

download.file(f, tmp)

img <- readBin(tmp, "raw", file.info(tmp)[1, "size"])

u <- paste0("https://vision.googleapis.com/v1/images:annotate?key=", CROWD_VISION_KEY)

body <- list(requests = list(image=list(content=img),

features=list(type="TEXT_DETECTION")

)

res <- POST(url=u,

encode="json",

body=body,

content_type_json()

)

res <- content(res)

cat(res$responses[[1]]$textAnnotations[[1]]$description[1])

メニュー
券売機番号
券売機番号
11小ラーメン
730円14大ラーメン
830円
12小ラーメン豚入り880円15大ラーメン誘入り 980円
17小つけ麺
880円18大つけ麺
980円
1 19小つけ麺豚入り 1030円20大つけ麺豚入り 1130円
13ミニラーメン 680円
21ビール
500円24黒ウーロン茶
200円
22うずら卵(5個) 100円23かつお魚粉
100円
02刻み生姜
50円
初めてのお客様は
かミニラーメンがお勧めです。
0000
3883
! 1
11ナナ
1つ1
売4580
4-3
券1112
22
00000000
38838005
7880651
睡! | |
| |
券11
1191220

メニュー

券売機番号

11小ラーメン

730円14大ラーメン

830円

12小ラーメン豚入り880円15大ラーメン誘入り 980円

17小つけ麺

880円18大つけ麺

980円

1 19小つけ麺豚入り 1030円20大つけ麺豚入り 1130円

13ミニラーメン 680円

21ビール

500円24黒ウーロン茶

200円

22うずら卵(5個) 100円23かつお魚粉

100円

02刻み生姜

50円

初めてのお客様は

かミニラーメンがお勧めです。

0000

3883

! 1

11ナナ

1つ1

売4580

4-3

券1112

00000000

38838005

7880651

睡! | |

| |

券11

1191220

以下の表は実際の表記との比較ですが、光の反射で見にくくなっているところ以外は正しい表現を抽出できているようです。

ExploratoryをさわってRでGoogleSearchConsoleのデータを集計・可視化

先日知り合った方から、Exploratoryの紹介をしていただき、Rをカジュアルに利用できる可視化ツールとのことで関心を持っていたのですが、さっそくこのブログのアクセス関連のデータを集計・可視化してみようと思います。

【目次】
・Exploratoryについて
・ダウンロードとインストール
・データの読み込み（Built in）と可視化
・データの読み込み（With R）と可視化
・dplyrの活用
・使ってみて思うところ

Exploratoryについて

ExploratoryはRを用いた、前処理・可視化・分析のためのデスクトップアプリケーションです。Amazon Redshift、Google BigQuery、PostgreSQL、MySQL、MongoDB、Presto、Google Analytics、Google Spreadsheet、Twitter、Web Page Scraping、CSV、Excel、JSON、Statistical files、R scriptなどの豊富なデータソースを利用できるようです。特に、Rのscriptでデータを取得できるのは魅力的だなぁと思いました。加えて、前処理に適したRの関数（dplyrなど）をカジュアルに扱えるようなコマンドラインもありますので、前処理しながら可視化するという作業がしやすいと思われます。分析に関しては、クラスター分析、分類、回帰、アソシエーションルール、相関、類似度、自然言語処理、文書感情分析などができるようです。まぁ、やりたい分析に関するパッケージは新しくRにインストールしてしまえばいいので、R使いにとって自由度の高い環境と言えます。

ダウンロードとインストール

こちらのExploratory公式サイト（ https://exploratory.io/download ）でemailアドレスなどのアカウント情報を入力してダウンロードできます。あとはインストールしてしまえば良いです。Rに関してもExploratoryに対応したバージョンをインストールしてくれます。後は、アプリケーションを開いて、Create Newのボタンを押して分析スタートです。

データの読み込み（Built in）と可視化

私はGoogle Analytics(GA)のデータをよく使うので、さっそくExploratoryで読み込んでみました。

自分のGAのプロパティやビューを選んで、任意の期間を選び、好きなDimensionやMetricsを選択して抽出されたデータを確認してOKすれば読み込みは完了です。非常に簡単でしたし、データを確認しながらDimensionやMetricsを選べるのは試行錯誤しやすく便利です。

ピボットテーブルが非常に使いやすく、デバイス別・ページ別の直帰率の集計も楽にできました。

データの読み込み（With R）と可視化

私が魅力に感じたのはデータ読み込みの際に、Rのscriptをそのまま使えることです。RMeCabを使えばテキストデータを形態素解析して可視化することもできますし、RFinanceなどで取得した株価のデータとかも使えるわけです。これまでに使ってたコードをそのまま使えるのは嬉しい限りです。今回は恥ずかしながら私のブログのGoogle検索結果のクエリをGoogle Search Console APIを使って取得し、可視化にチャレンジしてみます。

Google Search Consoleは以下のコマンドでAPIからデータ取得が可能です。

library("googleAuthR")
library("searchConsoleR")

gar_auth()
sc_websites <- list_websites()

own_site <- "http://kamonohashiperry.com/"
time_period_start <- "2016-07-01" ###任意の開始日を記入
time_period_end   <- "2016-09-20" ###任意の終了日を記入
filter_url <- c("page~~/")

dataset <- search_analytics(siteURL = sc_websites[match(own_site,sc_websites[,1]),1],
                 startDate = as.character(as.Date(time_period_start) ),
                 endDate = as.character(as.Date(time_period_end) ),
                 dimensions = c("query","page","device","date"),
                 dimensionFilterExp = filter_url,
                 searchType = "web", rowLimit = 5000)

library("googleAuthR")

library("searchConsoleR")

gar_auth()

sc_websites <- list_websites()

own_site <- "http://kamonohashiperry.com/"

time_period_start <- "2016-07-01" ###任意の開始日を記入

time_period_end <- "2016-09-20" ###任意の終了日を記入

filter_url <- c("page~~/")

dataset <- search_analytics(siteURL = sc_websites[match(own_site,sc_websites[,1]),1],

startDate = as.character(as.Date(time_period_start) ),

endDate = as.character(as.Date(time_period_end) ),

dimensions = c("query","page","device","date"),

dimensionFilterExp = filter_url,

searchType = "web", rowLimit = 5000)

Google検索クエリの検索結果上での私のブログのデバイス別、順位・CTR・インプレッション数・クリック数を3D散布図に描いてみました。

こんな感じで結果のグラフなどをシェアできるようです。触れば動きます。

時系列データに関してもサクサク作れたので、私のブログのデバイス別の平均順位を時系列で描写してみました。

カーソルを合わせたらデータの詳細が表示されます。

dplyrの活用

細かい前処理がしたい場合は、画面上にあるcommandのウインドウでdplyrのコマンドを叩けば痒いところに手が届かせることも可能です。

使ってみて思うところ

RやPythonを使っている我々にとって、Tableauを使うモチベーションが非常に低かったんですが、Rとつないで好き放題できる環境というのは非常に魅力的だと思いました。私は個人的にWebスクレイピングを結構しているので、テキストデータの可視化などにもチャレンジしてみたいです。数をこなして可視化のセンスを磨きたい。

統計的因果推論に関するスライドとRのサンプルコード

導入・入門から実践までのスライドと実践するためのコードを載せています。因果推論を実務でカジュアルに使えるまで上達させたいものです。

以下、
・統計的因果推論に関するスライド
・資料
・用語整理
・統計的因果推論に関するRの実行例
・書籍情報
からなります。

統計的因果推論に関するスライド

統計的因果推論の学習を進める上でのヒントが記されていると思われるスライドです。カジュアルに回帰分析を行うことに関する注意も書かれています。

統計的因果推論勉強会　第1回 from Hikaru GOTO

因果推論の歴史や利用することのモチベーションが非常にわかりやすく書かれています。
多重共線性を引き合いに出していたりしているのも理解が捗ります。

相関と因果について考える：統計的因果推論、その(不)可能性の中心 from takehikoihayashi

以下の2つは星野先生の『調査観察データの統計科学』通称赤本の1~3章に関するスライドです。

星野「調査観察データの統計科学」第1＆2章 from Shuyo Nakatani

星野「調査観察データの統計科学」第3章 from Shuyo Nakatani

Rによる実践に関してのスライドです。コードが書かれているので実行してみると良いでしょう。

傾向スコア：その概念とRによる実装 from takehikoihayashi

傾向スコアの計算をRで実践しているスライドです。

Rで学ぶ傾向スコア解析入門 – 無作為割り当てが出来ない時の因果効果推定 – from Yohei Sato

資料

こちらの資料は定義が書かれていて助かります。
2016/01/23 (Sat) 第 3 回因果推論を学ぶ会

こちらは論文や教科書の紹介もしています。
プロペンシティスコア（Propensity score; PS）（１）－PSの正しい使い方

カーネギーメロン大学の機械学習における因果推論に関しての講義資料です。
Lecture Notes 17 Causal Inference

Googleのハル・バリアンがまとめたペーパーも良いです。
Causal Inference in Economics and Marketing

用語整理

ATE（Average Treatment Effect：平均処置効果（因果効果））

例で述べるとするならば、母集団全てのユーザーにバナーを見せた場合のブランド名検索による訪問数の平均と、全てのユーザーにバナーを見せなかった場合のブランド名検索による訪問数の平均の差として表されます。バナーを見てしまったユーザーにとって、バナーを見なかったら、という反実仮想なデータは当然存在しないので、そのままでは計算できません。ただし、RCT（Randomized Control Trial）、無作為化比較対照実験、の状況ではバイアスなく推定できるとされています。RCTなケースは限られていると思いますが。

ATT（Average Treatment Effect on the Treated：処置群における平均処置効果）

バナーを見せたユーザーにおける、バナーを見せた場合と見せなかった場合の差の期待値。マーケティングにおける施策のROIを計算する際に使うことが望ましいとされています。ATEと同じくRCTにおいてバイアスなく推定できるとされています。

ATU（Average Treatment Effect on the Untreated：対照群における平均処置効果）

バナーを見せていないユーザーにおける、バナーを見せた場合と見せなかった場合の差の期待値。マーケティング施策を拡大させるか否かを判断する際に使うことができます。ATEと同じくRCTにおいてバイアスなく推定できるとされています。

強い意味での無視可能性

共変量に対し求める強い仮定のことで、「バナーを見たか見てないかのバイナリーな変数」や「ブランド名検索での訪問数」などに影響を与えるような共変量に対し、共変量自体で条件をつけて期待値をとると、「バナーを見たか見てないかのバイナリーな変数」と「潜在的なブランド名検索での訪問数」が独立するような特徴が求められています。「バナーを見たか見てないかのバイナリーな変数」が「過去のサイト訪問数（共変量）」や「特定ページへの接触（共変量）」で、配信対象を割り振られている場合は、そのバイナリーな変数は「潜在的なブランド名検索での訪問数」に影響を与えないとされています。

マッチング

バナーを見せられたユーザーの持つ、共変量（サイトへの訪問数や、見たページのカテゴリなど）の値と同じ（完全マッチング）、あるいは近い（距離を使ったマッチング）共変量を持っているが、バナーを見せられていない他のユーザーを「同じ人」と見なして、「バナーを見た・見てない」の与える「ブランド名検索での訪問数」への因果効果を推定します。

傾向スコア（Propensity score）

処置への割り当ての確率。つまり、上述の例でいうところの、バナーを見せられる確率。確率なので、当然0〜1の間の値をとります。推定には2項ロジットモデルが使われているようです。真の傾向スコアを推定できれば、ATE・ATT・ATUを計算することが可能になるそうです。この理屈はベイズの定理より導くことができるようです。詳しくは資料の”第 3 回因果推論を学ぶ会”を見てみてください。

統計的因果推論に関するRの実行例

“Rで学ぶ傾向スコア解析入門 – 無作為割り当てが出来ない時の因果効果推定”で紹介されていたコードを以下に掲載します。

#Matchingパッケージを呼び出す
library(Matching)

#サンプルデータを呼び出す
data(lalonde)

#データの確認
head(lalonde)

# age educ black hisp married nodegr re74 re75     re78 u74 u75 treat
# 1  37   11     1    0       1      1    0    0  9930.05   1   1     1
# 2  22    9     0    1       0      1    0    0  3595.89   1   1     1
# 3  30   12     1    0       0      0    0    0 24909.50   1   1     1
# 4  27   11     1    0       0      1    0    0  7506.15   1   1     1
# 5  33    8     1    0       0      1    0    0   289.79   1   1     1
# 6  22    9     1    0       0      1    0    0  4056.49   1   1     1

#ロジスティック回帰モデルを用いて傾向スコアを推定する。
#re78に対する効果を調べたいので、説明変数から除外しておく。
#treatはCMなどを見せたかどうかのバイナリーデータ。
logi <- glm(treat~., data=lalonde[, -9], family = binomial)

#因果効果を計算するために、マッチングを推定する。
#今回の効果を見るためのre78をYに指定して、CMなどを見せたかどうかのtreatをTrに指定する。
#共変量としてロジスティック回帰の予測値logi$fittedを用いている。
nsw1 <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =logi$fitted)

summary(nsw1)

# Estimate...  2138.6 
# AI SE......  797.76 
# T-stat.....  2.6807 
# p.val......  0.0073468 
# 
# Original number of observations..............  445 
# Original number of treated obs...............  185 
# Matched number of observations...............  185 
# Matched number of observations  (unweighted).  322

#マッチングのペアを確認する。
lalonde2 <- lalonde
lalonde2$id <- 1:nrow(lalonde2)
lalonde2$score <- logi$fitted

pair.df <- cbind(lalonde2[nsw1$index.treated, c("id","score")],
                 lalonde2[nsw1$index.control, c("id", "score")])

names(pair.df) <- c("t.id", "t.score", "c.id", "c.score")

head(pair.df)

# t.id   t.score c.id   c.score
# 1      1 0.3927536  357 0.3935865
# 2      2 0.2271642  231 0.2242716
# 3      3 0.5313484  261 0.5313484
# 4      4 0.3285956  254 0.3285956
# 4.1    4 0.3285956  328 0.3286097
# 4.2    4 0.3285956  333 0.3286097
#確かに近そうだ。

#キャリパーマッチング（ペアが特定の距離以上になる時はマッチングしないマッチング）
nsw2 <- Match(Y =lalonde2$re78, Tr=lalonde2$treat, X =logi$fitted, caliper = T)
summary(nsw2)

# Estimate...  2138.6 
# AI SE......  797.76 
# T-stat.....  2.6807 
# p.val......  0.0073468 
# 
# Original number of observations..............  445 
# Original number of treated obs...............  185 
# Matched number of observations...............  185 
# Matched number of observations  (unweighted).  322 
# 
# Caliper (SDs)........................................   TRUE 
# Number of obs dropped by 'exact' or 'caliper'  0

#IPW(Inverse Probability Weight)の計算
ivec1 <- lalonde$treat
ivec2 <- rep(1, nrow(lalonde)) - ivec1

ivec <- cbind(ivec1, ivec2)
head(ivec)

# ivec1 ivec2
# [1,]     1     0
# [2,]     1     0
# [3,]     1     0
# [4,]     1     0
# [5,]     1     0
# [6,]     1     0

iestp1 <- (ivec1/logi$fitted)*(length(ivec1)/sum(ivec1))
iestp2 <- (ivec2/logi$fitted)*(length(ivec2)/sum(ivec2))

iestp <- iestp1 + iestp2
head(iestp)

# 1         2         3         4         5         6 
# 6.124465 10.588840  4.526984  7.320260  6.053291  6.745951 

ipwe <- lm(re78~ivec - 1, weights=iestp, data = lalonde)
summary(ipwe)

# Call:
#   lm(formula = re78 ~ ivec - 1, data = lalonde, weights = iestp)
# 
# Weighted Residuals:
#   Min     1Q Median     3Q    Max 
# -17241  -9985  -4048   6201 129799 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(&gt;|t|)    
# ivecivec1   6213.0      462.2   13.44   &lt;2e-16 ***
#   ivecivec2   4589.4      436.4   10.52   &lt;2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 15090 on 443 degrees of freedom
# Multiple R-squared:  0.3967,  Adjusted R-squared:  0.3939 
# F-statistic: 145.6 on 2 and 443 DF,  p-value: &lt; 2.2e-16

#因果効果
causal_effect <- ipwe$coefficients[1]- ipwe$coefficients[2]
causal_effect

# ivecivec1 
# 1623.581

#標準誤差
std_causal_effect <- sqrt(summary(ipwe)$coefficients[3]^2 + summary(ipwe)$coefficients[4]^2)
std_causal_effect

# [1] 635.7069

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

#Matchingパッケージを呼び出す

library(Matching)

#サンプルデータを呼び出す

data(lalonde)

#データの確認

head(lalonde)

# age educ black hisp married nodegr re74 re75 re78 u74 u75 treat

# 1 37 11 1 0 1 1 0 0 9930.05 1 1 1

# 2 22 9 0 1 0 1 0 0 3595.89 1 1 1

# 3 30 12 1 0 0 0 0 0 24909.50 1 1 1

# 4 27 11 1 0 0 1 0 0 7506.15 1 1 1

# 5 33 8 1 0 0 1 0 0 289.79 1 1 1

# 6 22 9 1 0 0 1 0 0 4056.49 1 1 1

#ロジスティック回帰モデルを用いて傾向スコアを推定する。

#re78に対する効果を調べたいので、説明変数から除外しておく。

#treatはCMなどを見せたかどうかのバイナリーデータ。

logi <- glm(treat~., data=lalonde[, -9], family = binomial)

#因果効果を計算するために、マッチングを推定する。

#今回の効果を見るためのre78をYに指定して、CMなどを見せたかどうかのtreatをTrに指定する。

#共変量としてロジスティック回帰の予測値logi$fittedを用いている。

nsw1 <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =logi$fitted)

summary(nsw1)

# Estimate... 2138.6

# AI SE...... 797.76

# T-stat..... 2.6807

# p.val...... 0.0073468

# Original number of observations.............. 445

# Original number of treated obs............... 185

# Matched number of observations............... 185

# Matched number of observations (unweighted). 322

#マッチングのペアを確認する。

lalonde2 <- lalonde

lalonde2$id <- 1:nrow(lalonde2)

lalonde2$score <- logi$fitted

pair.df <- cbind(lalonde2[nsw1$index.treated, c("id","score")],

lalonde2[nsw1$index.control, c("id", "score")])

names(pair.df) <- c("t.id", "t.score", "c.id", "c.score")

head(pair.df)

# t.id t.score c.id c.score

# 1 1 0.3927536 357 0.3935865

# 2 2 0.2271642 231 0.2242716

# 3 3 0.5313484 261 0.5313484

# 4 4 0.3285956 254 0.3285956

# 4.1 4 0.3285956 328 0.3286097

# 4.2 4 0.3285956 333 0.3286097

#確かに近そうだ。

#キャリパーマッチング（ペアが特定の距離以上になる時はマッチングしないマッチング）

nsw2 <- Match(Y =lalonde2$re78, Tr=lalonde2$treat, X =logi$fitted, caliper = T)

summary(nsw2)

# Estimate... 2138.6

# AI SE...... 797.76

# T-stat..... 2.6807

# p.val...... 0.0073468

# Original number of observations.............. 445

# Original number of treated obs............... 185

# Matched number of observations............... 185

# Matched number of observations (unweighted). 322

# Caliper (SDs)........................................ TRUE

# Number of obs dropped by 'exact' or 'caliper' 0

#IPW(Inverse Probability Weight)の計算

ivec1 <- lalonde$treat

ivec2 <- rep(1, nrow(lalonde)) - ivec1

ivec <- cbind(ivec1, ivec2)

head(ivec)

# ivec1 ivec2

# [1,] 1 0

# [2,] 1 0

# [3,] 1 0

# [4,] 1 0

# [5,] 1 0

# [6,] 1 0

iestp1 <- (ivec1/logi$fitted)*(length(ivec1)/sum(ivec1))

iestp2 <- (ivec2/logi$fitted)*(length(ivec2)/sum(ivec2))

iestp <- iestp1 + iestp2

head(iestp)

# 1 2 3 4 5 6

# 6.124465 10.588840 4.526984 7.320260 6.053291 6.745951

ipwe <- lm(re78~ivec - 1, weights=iestp, data = lalonde)

summary(ipwe)

# Call:

# lm(formula = re78 ~ ivec - 1, data = lalonde, weights = iestp)

# Weighted Residuals:

# Min 1Q Median 3Q Max

# -17241 -9985 -4048 6201 129799

# Coefficients:

# Estimate Std. Error t value Pr(>|t|)

# ivecivec1 6213.0 462.2 13.44 <2e-16 ***

# ivecivec2 4589.4 436.4 10.52 <2e-16 ***

# ---

# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Residual standard error: 15090 on 443 degrees of freedom

# Multiple R-squared: 0.3967, Adjusted R-squared: 0.3939

# F-statistic: 145.6 on 2 and 443 DF, p-value: < 2.2e-16

#因果効果

causal_effect <- ipwe$coefficients[1]- ipwe$coefficients[2]

causal_effect

# ivecivec1

# 1623.581

#標準誤差

std_causal_effect <- sqrt(summary(ipwe)$coefficients[3]^2 + summary(ipwe)$coefficients[4]^2)

std_causal_effect

# [1] 635.7069

書籍情報

データ分析の力　因果関係に迫る思考法 (光文社新書)
「原因と結果」の経済学―――データから真実を見抜く思考法
 岩波データサイエンス Vol.3
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

Tokyo.R#53で得たパッケージ情報とその実践

第53回のTokyo.Rで気になったパッケージの情報と実行例をいくつかあげました。スライドなどもろもろの発表はこちらの方のブログ「第53回R勉強会@東京で発表してきた」が非常に詳しく書かれています。

【目次】
・ggradarパッケージ
・proxyパッケージ
・因果推論（CBPSパッケージ）
・MXNetパッケージ
・missForestパッケージ
・RFinanceパッケージ

ggradarパッケージ

簡単にレーダーチャートを作れるパッケージです。こちらのブログを参考にしています。

install.packages("devtools")
devtools::install_github("ricardo-bion/ggradar")

1 2	install.packages("devtools") devtools::install_github("ricardo-bion/ggradar")

企業の職場環境に関してまとめられた某口コミサイトから4個ほどデータを拝借してきました。

> CompanyVoiceData
  company growth stability salary rewarding idea difficulty welfare education
1  google    5.0       5.0    4.9       5.0  4.3        5.0     5.0       4.6
2   yahoo    3.9       5.0    3.2       3.8  3.7        3.9     3.1       3.3
3 recruit    4.4       4.8    5.0       5.0  5.0        5.0     4.0       5.0
4  amazon    5.0       5.0    4.2       4.0  4.2        5.0     3.6       3.3

> CompanyVoiceData

company growth stability salary rewarding idea difficulty welfare education

1 google 5.0 5.0 4.9 5.0 4.3 5.0 5.0 4.6

2 yahoo 3.9 5.0 3.2 3.8 3.7 3.9 3.1 3.3

3 recruit 4.4 4.8 5.0 5.0 5.0 5.0 4.0 5.0

4 amazon 5.0 5.0 4.2 4.0 4.2 5.0 3.6 3.3

ggradarをそのまま使おうとすると、Circular Air Lightというフォントが必要だと怒られるので、参考のブログにある通り、OSXの場合はこちらをダブルクリックでインストールして再起動します。

先ほどのデータに対して、以下のコードを実行すれば非常に簡単にレーダーチャートが作れました。

library("ggradar")
CompanyVoiceData <- data.frame(read.csv(file ="company_voice.csv",header = TRUE))

ggradar(CompanyVoiceData, 
        grid.max = max(CompanyVoiceData[, 2:ncol(CompanyVoiceData)]),
        background.circle.colour = "#ffdd99", #背景色の指定
        background.circle.transparency = 1, #背景色の透明度を指定
        group.line.width = 2, #線の太さの指定
        group.point.size = 6, #シンボルの大きさの指定
        axis.label.size = 5, #軸ラベルサイズの指定
        gridline.min.colour = "#4b61ba", #最小円の線色の指定
        gridline.mid.colour = "#a87963", #中円の線色の指定
        gridline.max.colour = "#e1e6ea", #最大円の線色の指定
        grid.line.width = 1.5, #各円の線の太さの指定
        gridline.min.linetype = "longdash", #線種の指定
        gridline.mid.linetype = "longdash", #線種の指定
        gridline.max.linetype = "longdash") #線種の指定

library("ggradar")

CompanyVoiceData <- data.frame(read.csv(file ="company_voice.csv",header = TRUE))

ggradar(CompanyVoiceData,

grid.max = max(CompanyVoiceData[, 2:ncol(CompanyVoiceData)]),

background.circle.colour = "#ffdd99", #背景色の指定

background.circle.transparency = 1, #背景色の透明度を指定

group.line.width = 2, #線の太さの指定

group.point.size = 6, #シンボルの大きさの指定

axis.label.size = 5, #軸ラベルサイズの指定

gridline.min.colour = "#4b61ba", #最小円の線色の指定

gridline.mid.colour = "#a87963", #中円の線色の指定

gridline.max.colour = "#e1e6ea", #最大円の線色の指定

grid.line.width = 1.5, #各円の線の太さの指定

gridline.min.linetype = "longdash", #線種の指定

gridline.mid.linetype = "longdash", #線種の指定

gridline.max.linetype = "longdash") #線種の指定

proxyパッケージ

距離や類似度を計算するパッケージです。
先ほどのデータに対して類似度と距離を計算してみます。

library(proxy)
> simil(CompanyVoiceData[,-1])
          1         2         3
2 0.2286639                    
3 0.6373648 0.1339713          
4 0.6499133 0.5787506 0.4188571
> dist(CompanyVoiceData[,-1])
         1        2        3
2 3.522783                  
3 1.435270 3.401470         
4 2.269361 1.989975 2.393742

library(proxy)

> simil(CompanyVoiceData[,-1])

1 2 3

2 0.2286639

3 0.6373648 0.1339713

4 0.6499133 0.5787506 0.4188571

> dist(CompanyVoiceData[,-1])

1 2 3

2 3.522783

3 1.435270 3.401470

4 2.269361 1.989975 2.393742

こんな感じで、類似度や距離の計算ができます。

因果推論

こちらはパッケージとかそういうものではなく、既存の関数などで計算できるようです。
こちらのブログ、「調査観察データにおける因果推論(3) – Rによる傾向スコア，IPW推定量，二重にロバストな推定量の算出」に詳しく書かれています。
・glm関数での傾向スコアの算出
・傾向スコアを共変量としてlm関数で回帰分析
・コードを愚直に書いてIPW推定量の算出
・期待値の標準誤差を出すための関数を作成
・DR推定量の算出をするための関数を作成
などで、推定自体は実現できるようです。

ただし、CBPS(Covariate Balancing Propensity Score)というパッケージがあるらしく、このパッケージを用いれば因果推論の計算を行えるようです。

Package ‘CBPS’
以下のようなExampleコードが載っていたので、実行してみましたが、なかなか結果が返ってこなかったので不安になりました。計算が終わるまで10分以上はかかったと思います。

library(CBPS)
data(Blackwell)

form1<-"d.gone.neg ~ d.gone.neg.l1 + d.gone.neg.l2 + d.neg.frac.l3 + camp.length + camp.length +
deminc + base.poll + year.2002 + year.2004 + year.2006 + base.und + office"

##Fitting the models in Imai and Ratkovic (2014)
##Warning: may take a few mintues; setting time.vary to FALSE
##Results in a quicker fit but with poorer balance
fit1 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",
            iterations = NULL, twostep = TRUE, msm.variance = "full", time.vary = TRUE)
fit2 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",
            iterations = NULL, twostep = TRUE, msm.variance = "approx", time.vary = TRUE)

##Assessing balance
bal1 <- balance(fit1)
bal2 <- balance(fit2)

##Effect estimation: Replicating Effect Estimates in
##Table 3 of Imai and Ratkovic (2014)
lm1 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$glm.weights)
lm2 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$weights)
lm3 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit2$weights)
lm4 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$glm.weights)
lm5 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$weights)
lm6 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit2$weights)

library(CBPS)

data(Blackwell)

form1<-"d.gone.neg ~ d.gone.neg.l1 + d.gone.neg.l2 + d.neg.frac.l3 + camp.length + camp.length +

deminc + base.poll + year.2002 + year.2004 + year.2006 + base.und + office"

##Fitting the models in Imai and Ratkovic (2014)

##Warning: may take a few mintues; setting time.vary to FALSE

##Results in a quicker fit but with poorer balance

fit1 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",

iterations = NULL, twostep = TRUE, msm.variance = "full", time.vary = TRUE)

fit2 <- CBMSM(formula = form1, time=Blackwell$time,id=Blackwell$demName,data=Blackwell, type="MSM",

iterations = NULL, twostep = TRUE, msm.variance = "approx", time.vary = TRUE)

##Assessing balance

bal1 <- balance(fit1)

bal2 <- balance(fit2)

##Effect estimation: Replicating Effect Estimates in

##Table 3 of Imai and Ratkovic (2014)

lm1 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$glm.weights)

lm2 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit1$weights)

lm3 <- lm(demprcnt[time==1]~fit1$treat.hist,data=Blackwell,weights=fit2$weights)

lm4 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$glm.weights)

lm5 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit1$weights)

lm6 <- lm(demprcnt[time==1]~fit1$treat.cum,data=Blackwell,weights=fit2$weights)

MXNet

XGBoostのパッケージを作ったチームが手がけているパッケージで、深層学習を実行できます。

インストール方法はここに書かれています。
Deep Learning for R

install.packages("drat", repos="https://cran.rstudio.com")
drat:::addRepo("dmlc")
install.packages("mxnet")

install.packages("drat", repos="https://cran.rstudio.com")

drat:::addRepo("dmlc")

install.packages("mxnet")

あれ、OSXではエラーが返ってきてライブラリが読み込めないですね。どうやら私のためにあるようなブログ「Installing mxnet for R on Yosemite」があったので、時間を見つけてチャレンジしてみようと思います。

ディープラーニングを用いた回帰分析については、Neural Network with MXNet in Five Minutesにコードがもろもろ載っていますので、チャレンジしてみると良いと思います。

リンク先に載っているのですが、一応コードを以下に記しておきます。

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)
train.x = data.matrix(BostonHousing[train.ind, -14])
train.y = BostonHousing[train.ind, 14]
test.x = data.matrix(BostonHousing[-train.ind, -14])
test.y = BostonHousing[-train.ind, 14]

# Define the input data
data <- mx.symbol.Variable("data")
# A fully connected hidden layer
# data: input source
# num_hidden: number of neurons in this hidden layer
fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# Use linear regression for the output layer
lro <- mx.symbol.LinearRegressionOutput(fc1)

preds = predict(model, test.x)

## Auto detect layout of input matrix, use rowmajor..
sqrt(mean((preds-test.y)^2))

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {
  res <- mean(abs(label-pred))
  return(res)
})

mx.set.seed(0)
model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y,
                                     ctx=mx.cpu(), num.round=50, array.batch.size=20,
                                     learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)

data(BostonHousing, package="mlbench")

train.ind = seq(1, 506, 3)

train.x = data.matrix(BostonHousing[train.ind, -14])

train.y = BostonHousing[train.ind, 14]

test.x = data.matrix(BostonHousing[-train.ind, -14])

test.y = BostonHousing[-train.ind, 14]

# Define the input data

data <- mx.symbol.Variable("data")

# A fully connected hidden layer

# data: input source

# num_hidden: number of neurons in this hidden layer

fc1 <- mx.symbol.FullyConnected(data, num_hidden=1)

# Use linear regression for the output layer

lro <- mx.symbol.LinearRegressionOutput(fc1)

preds = predict(model, test.x)

## Auto detect layout of input matrix, use rowmajor..

sqrt(mean((preds-test.y)^2))

demo.metric.mae <- mx.metric.custom("mae", function(label, pred) {

res <- mean(abs(label-pred))

return(res)

})

mx.set.seed(0)

model <- mx.model.FeedForward.create(lro, X=train.x, y=train.y,

ctx=mx.cpu(), num.round=50, array.batch.size=20,

learning.rate=2e-6, momentum=0.9, eval.metric=demo.metric.mae)

missForest

ランダムフォレストを用いて、欠損値補完を行うためのパッケージです。目的変数が欠損していても適用できるようです。
詳しくは、スライドを見ていただいた方がいいですが、以下のプログラムで実行できました。ちなみにスライドはこちら、「Imputation of Missing Values using Random Forest」

library(missForest)
library(dplyr)

#ggplot2のデータセットを読み込む
data(diamonds, package = "ggplot2")
dia.sample <- sample_n(diamonds, size=2000)
dia.sample <- as.data.frame(dia.sample)

#既存データセットに5%の欠損を与える
dia.mis <- prodNA(dia.sample, noNA=0.05)

#補完の実行
dia.imp <- missForest(dia.mis, verbose=TRUE)
dia.imp %>% str(max.level=1)

#補完精度の推定
dia.imp$OOBerror
dia.imp <- missForest(dia.mis, verbose=TRUE, variablewise=TRUE)

#補完精度の検証
mixError(ximp = dia.imp$ximp,
         xmis = dia.mis,
         xtrue = dia.sample)

library(missForest)

library(dplyr)

#ggplot2のデータセットを読み込む

data(diamonds, package = "ggplot2")

dia.sample <- sample_n(diamonds, size=2000)

dia.sample <- as.data.frame(dia.sample)

#既存データセットに5%の欠損を与える

dia.mis <- prodNA(dia.sample, noNA=0.05)

#補完の実行

dia.imp <- missForest(dia.mis, verbose=TRUE)

dia.imp %>% str(max.level=1)

#補完精度の推定

dia.imp$OOBerror

dia.imp <- missForest(dia.mis, verbose=TRUE, variablewise=TRUE)

#補完精度の検証

mixError(ximp = dia.imp$ximp,

xmis = dia.mis,

xtrue = dia.sample)

RFinanceYJ

Yohei Sato, Nobuaki Oshiro, Shinichi Takayanagiさんたちが作った、Yahoo!ファイナンスの株価データを取得できるパッケージです。だいぶ前からあったようですが、使って分析している人は初めて見ました。どうやらYahoo!ファイナンスの仕様によって書き換えていかないといけないようです。「2015-01-20 Rでチャートを書いてみる(9)」のブログに実行可能なプログラムがあります。以下、実行可能なコードを転載いたします。

library(RFinanceYJ)

#API
quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.stock <- function(quote.table.item){
    if( xmlSize(quote.table.item) < 5) return(NULL) 
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o <- as.number(xmlValue(quote.table.item[[2]]))
    h <- as.number(xmlValue(quote.table.item[[3]]))
    l <- as.number(xmlValue(quote.table.item[[4]]))
    c <- as.number(xmlValue(quote.table.item[[5]]))
    v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)
    a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)
    return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))
  }
  return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))
}
quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.fund <- function(quote.table.item){
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    if(time.interval=='monthly'){
      d <- endOfMonth(d)
    }
    c <- as.number(xmlValue(quote.table.item[[2]]))
    v <- as.number(xmlValue(quote.table.item[[3]]))
    return(data.frame(date=d,constant.value=c,NAV=v))
  }
  return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))
}
quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')
{
  time.interval <- substr(time.interval,1,1)
  function.fx <- function(quote.table.item){
    d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)
    o <- as.number(xmlValue(quote.table.item[[2]]))
    h <- as.number(xmlValue(quote.table.item[[3]]))
    l <- as.number(xmlValue(quote.table.item[[4]]))
    c <- as.number(xmlValue(quote.table.item[[5]]))
    return(data.frame(date=d,open=o,high=h,low=l,close=c))
  }
  return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))
}
######  private functions  #####
#get time series data from Yahoo! Finance.
quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){
  r <- NULL
  result.num <- 51
  financial.data <- data.frame(NULL)
  #start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&c=\\1&a=\\2&b=\\3",since))
  #end   <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&f=\\1&d=\\2&e=\\3",date.end))
  start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))
  end   <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))
  
  if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")
  
  extractQuoteTable <- function(r,type){
    if(type %in% c("fund","fx")){
      tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]
    }
    else{
      tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]
    }
    return(tbl)
  }
  
  #while( result.num >= 51 ){
  while(1){
    start.num <- start.num + 1
    quote.table <- NULL
    quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")
    #cat(quote.url)
    #try( r <- xmlRoot(htmlTreeParse(quote.url,error=xmlErrorCumulator(immediate=F))), TRUE)  # これだと取得時にエラーが出た。。
    try(r<-htmlParse(quote.url))
    if( is.null(r) ) stop(paste("Can not access :", quote.url))
    
    #try( quote.table <- r[[2]][[1]][[1]][[16]][[1]][[1]][[1]][[4]][[1]][[1]][[1]], TRUE )
    #try( quote.table <- extractQuoteTable(r,type), TRUE )
    try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )
    
    quote.size<-xmlSize(quote.table)
    #cat(paste("size:",quote.size))
    if(xmlSize(quote.table)<=1){
      return (financial.data)
    }
    if( is.null(quote.table) ){
      if( is.null(financial.data) ){
        stop(paste("Can not quote :", x))
      }else{
        financial.data <- financial.data[order(financial.data$date),]
        return(financial.data)
      }
    }
    
    size <- xmlSize(quote.table)
    for(i in 2:size){
      financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))
    }
    
    #result.num <- xmlSize(quote.table)
    Sys.sleep(1)
  }
  financial.data <- financial.data[order(financial.data$date),]
  return(financial.data)  
}
#convert string formart date to POSIXct object
convertToDate <- function(date.string,time.interval)
{
  #data format is different between monthly and dialy or weekly
  if(any(time.interval==c('d','w'))){
    result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)
  }else if(time.interval=='m'){
    result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)
  }
  return(as.POSIXct(result))
}
#convert string to number.
as.number <- function(string)
{
  return(as.double(as.character(gsub("[^0-9.]", "",string))))
}
#return end of month date.
endOfMonth <- function(date.obj)
{
  startOfMonth     <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")
  startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")
  return(startOfNextMonth-1)
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

library(RFinanceYJ)

#API

quoteStockTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.stock <- function(quote.table.item){

if( xmlSize(quote.table.item) < 5) return(NULL)

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

v <- ifelse(xmlSize(quote.table.item) >= 6,as.number(xmlValue(quote.table.item[[6]])),0)

a <- ifelse(xmlSize(quote.table.item) >= 7,as.number(xmlValue(quote.table.item[[7]])),0)

return(data.frame(date=d,open=o,high=h,low=l,close=c,volume=v, adj_close=a))

}

return(quoteTsData(x,function.stock,since,start.num,date.end,time.interval,type="stock"))

}

quoteFundTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fund <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

if(time.interval=='monthly'){

d <- endOfMonth(d)

}

c <- as.number(xmlValue(quote.table.item[[2]]))

v <- as.number(xmlValue(quote.table.item[[3]]))

return(data.frame(date=d,constant.value=c,NAV=v))

}

return(quoteTsData(x,function.fund,since,start.num,date.end,time.interval,type="fund"))

}

quoteFXTsData <- function(x, since=NULL,start.num=0,date.end=NULL,time.interval='daily')

{

time.interval <- substr(time.interval,1,1)

function.fx <- function(quote.table.item){

d <- convertToDate(xmlValue(quote.table.item[[1]]),time.interval)

o <- as.number(xmlValue(quote.table.item[[2]]))

h <- as.number(xmlValue(quote.table.item[[3]]))

l <- as.number(xmlValue(quote.table.item[[4]]))

c <- as.number(xmlValue(quote.table.item[[5]]))

return(data.frame(date=d,open=o,high=h,low=l,close=c))

}

return(quoteTsData(x,function.fx,since,start.num,date.end,time.interval,type="fx"))

}

###### private functions #####

#get time series data from Yahoo! Finance.

quoteTsData <- function(x,function.financialproduct,since,start.num,date.end,time.interval,type="stock"){

r <- NULL

result.num <- 51

financial.data <- data.frame(NULL)

#start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&c=\\1&a=\\2&b=\\3",since))

#end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&f=\\1&d=\\2&e=\\3",date.end))

start <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&sy=\\1&sm=\\2&sd=\\3",since))

end <- (gsub("([0-9]{4,4})-([0-9]{2,2})-([0-9]{2,2})","&ey=\\1&em=\\2&ed=\\3",date.end))

if(!any(time.interval==c('d','w','m'))) stop("Invalid time.interval value")

extractQuoteTable <- function(r,type){

if(type %in% c("fund","fx")){

tbl <- r[[2]][[2]][[7]][[3]][[3]][[9]][[2]]

}

else{

tbl <- r[[2]][[2]][[7]][[3]][[3]][[10]][[2]]

}

return(tbl)

}

#while( result.num >= 51 ){

while(1){

start.num <- start.num + 1

quote.table <- NULL

quote.url <- paste('http://info.finance.yahoo.co.jp/history/?code=',x,start,end,'&p=',start.num,'&tm=',substr(time.interval,1,1),sep="")

#cat(quote.url)

#try( r <- xmlRoot(htmlTreeParse(quote.url,error=xmlErrorCumulator(immediate=F))), TRUE) # これだと取得時にエラーが出た。。

try(r<-htmlParse(quote.url))

if( is.null(r) ) stop(paste("Can not access :", quote.url))

#try( quote.table <- r[[2]][[1]][[1]][[16]][[1]][[1]][[1]][[4]][[1]][[1]][[1]], TRUE )

#try( quote.table <- extractQuoteTable(r,type), TRUE )

try( quote.table <- xpathApply(r,"//table")[[2]], TRUE )

quote.size<-xmlSize(quote.table)

#cat(paste("size:",quote.size))

if(xmlSize(quote.table)<=1){

return (financial.data)

}

if( is.null(quote.table) ){

if( is.null(financial.data) ){

stop(paste("Can not quote :", x))

}else{

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

size <- xmlSize(quote.table)

for(i in 2:size){

financial.data <- rbind(financial.data,function.financialproduct(quote.table[[i]]))

}

#result.num <- xmlSize(quote.table)

Sys.sleep(1)

}

financial.data <- financial.data[order(financial.data$date),]

return(financial.data)

}

#convert string formart date to POSIXct object

convertToDate <- function(date.string,time.interval)

{

#data format is different between monthly and dialy or weekly

if(any(time.interval==c('d','w'))){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-\\5",date.string)

}else if(time.interval=='m'){

result <- gsub("^([0-9]{4})([^0-9]+)([0-9]{1,2})([^0-9]+)","\\1-\\3-01",date.string)

}

return(as.POSIXct(result))

}

#convert string to number.

as.number <- function(string)

{

return(as.double(as.character(gsub("[^0-9.]", "",string))))

}

#return end of month date.

endOfMonth <- function(date.obj)

{

startOfMonth <- as.Date(format(date.obj,"%Y%m01"),"%Y%m%d")

startOfNextMonth <- as.Date(format(startOfMonth+31,"%Y%m01"),"%Y%m%d")

return(startOfNextMonth-1)

}

このコードでYahoo!ジャパンの株価を見てみましょう。ちなみに番号は4689です。どうやら上手く取れているようです。

> quoteStockTsData("4689.t",since="2016-01-01")
         date open high low close   volume adj_close
1  2016-05-02  476  483 475   478 18498100       478
2  2016-04-28  504  508 493   496 11966300       496
3  2016-04-27  505  511 495   497 12973800       497
4  2016-04-26  507  508 495   500  7712600       500
5  2016-04-25  513  515 506   509  7350600       509
6  2016-04-22  515  517 509   514  8908900       514
7  2016-04-21  512  517 506   514 13249900       514
8  2016-04-20  511  515 493   506 14455700       506
9  2016-04-19  516  523 511   516 13345800       516
10 2016-04-18  503  509 499   503 10275900       503
11 2016-04-15  504  519 504   513 16962900       513

> quoteStockTsData("4689.t",since="2016-01-01")

date open high low close volume adj_close

1 2016-05-02 476 483 475 478 18498100 478

2 2016-04-28 504 508 493 496 11966300 496

3 2016-04-27 505 511 495 497 12973800 497

4 2016-04-26 507 508 495 500 7712600 500

5 2016-04-25 513 515 506 509 7350600 509

6 2016-04-22 515 517 509 514 8908900 514

7 2016-04-21 512 517 506 514 13249900 514

8 2016-04-20 511 515 493 506 14455700 506

9 2016-04-19 516 523 511 516 13345800 516

10 2016-04-18 503 509 499 503 10275900 503

11 2016-04-15 504 519 504 513 16962900 513