Rで重回帰分析〜実践編

アクセス分析を恒常的にしていると見えてくるものもあります。
アクセスデータの分析だけでは因子が足りず、結局机上でちょこちょこ数字を弄って終わるということ。
アクセス分析において重要なことは、以下の3つかなと考えています。

1.データのサマリー解析
2.回帰モデル構築
3.予測モデル構築
そして、これらの検証。

仮説なき分析は分析にあらず。

ということで、実戦でアクセス分析をおこなっていく際の手順についてまとめてみようと思います。

いつものように私の塩っぱいアクセスデータを使用することにします。

データは、ブログ公開から昨日までの日付付きデータ。
変数は、ページビュー、訪問数、新規訪問者の割合、平均ページ滞在時間、直帰率。
データを扱いやすくするためにこれらの変数を以下に編集します。

ページビュー ⇛PV
訪問数 ⇛Visit
新規訪問者の割合 ⇛NewerR
平均ページ滞在時間 ⇛Avetime
直帰率 ⇛BR

早速Rを使って処理をしていきます。

まずは、GoogleAnalyticsのカスタムレポートで編集した後にCSVにエクスポートしたデータをRコマンダーでインポートしていきます。

//Rコマンダーの起動
> library(Rcmdr)

Rコマンダーが起動しますので、データ > データのインポート でインポートしていきますが、GoogleAnalyticsでエクスポートしたCSVファイルはカンマ区切りのデータになるので、「フィールドの区切り記号」を「カンマ」にチェックします。

無事インポートが終了したら、データ内容をチェックします。
フィールドの区切りをしくじると一列に全データが入ってしまうので注意が必要です。

> summary(datad)
Date PV Visit NewerR
Min. :20130302 Min. : 2.00 Min. : 1.00 Min. :0.0000
1st Qu.:20130512 1st Qu.: 15.50 1st Qu.: 9.00 1st Qu.:0.5345
Median :20130723 Median : 30.00 Median : 17.00 Median :0.6667
Mean :20130739 Mean : 34.65 Mean : 18.95 Mean :0.6253
3rd Qu.:20131002 3rd Qu.: 45.00 3rd Qu.: 27.00 3rd Qu.:0.7734
Max. :20131213 Max. :381.00 Max. :176.00 Max. :1.0000
Avetime BR
Min. :0.000000 Min. :0.0000
1st Qu.:0.000515 1st Qu.:0.6793
Median :0.001065 Median :0.7826
Mean :0.001618 Mean :0.7594
3rd Qu.:0.001881 3rd Qu.:0.8571
Max. :0.018437 Max. :1.0000

こんな感じできちんと表示されれば、データが無事インポートされたことがわかります。
通常のウェブサイトでは、新規訪問者の割合や直帰率が「1」となることは極めて稀ですが、私のアクセス数ではこんなものでしょう。日付も開始日がMinとして、終了日がMaxとして間違いなく入っていることがわかります。
さらに視覚化をおこなっていきます。
視覚化は、plot関数でもできますが、私は欲張りなので、相関係数もヒストグラムもデータのバラつき具合もわかるpsychライブラリーを使って表示させています。
このライブラリーを使う前に、日付データが邪魔なので、新しいクラスを作って日付を除いたデータの塊を作ります。

> data.set <- datad[c("PV","Visit","NewerR","Avetime","BR")] > library(psych)
> pairs.panels(data.set, smooth=FALSE, density=FALSE, ellipses=FALSE, scale=TRUE)

出力結果は、このような感じ。
2013-12-14 14.29.42

見方ですが、数字の部分が相関係数、棒グラフがヒストグラム、黒い点がデータのバラつき具合になります。
まずは相関係数ですが、psychライブラリーの素晴らしいところが高い数字のフォントサイズが大きくなっているところ。
相関がもっとも高いのが、PVとVisitの0.84、続いてVisitとNewerRの0.39ですが、0.5以下はほぼ関連性はないと判断できるので(0.6以下という考えもありますが実務的に0.6も範囲に入れた方が話が早い場合も多く私の場合は0.5以下です)、これらの変数の内関係性があるのはPVとVisitくらいということがわかります。

続いて、ヒストグラムの状況ですが、PVやVisit、Avetimeが左寄り、NewerRやBRが右寄りということが見えてきます。
PVが左寄りということは、PVが少ないです、ということです(涙)。
NewerRが右寄りということは新規ユーザが非常に多いということ、BRの場合は直帰率が高いですよ、ということです(涙)。

データのバラつき度合いですが、PVとVisitの例(表の1列2行目)を例にみていくと、縦軸がPV、横軸がVisitということなので、Visit数が少ないときPVが少なくなるという傾向が見えてきます。次にNewerRとAvetimeの関係を見ていくと、縦軸がNewerRで横j区がAvetimeとすると、Avetimeが多くなるとNewerRが増えていきますが、一定のAvetimeになると徐々にNewerRが減少していくということがわかってきます。

これらでデータの傾向値を分析していき、データの状況を把握していきます。

続いて回帰モデルを作っていきます。
> dataset.lm <- lm(Visit~.,data=data.set) dataset.lmという変数に、lm関数で変数Visitを目的変数とした回帰分析を格納します。 続いて説明変数の説明度を調査していきます。 回帰分析は説明変数の組み合わせが重要となり、無駄な説明変数があるとモデルに悪影響を及ぼします。 モデルの精度を調べる方法としてAICの値による判定方法があります。 AICは値が小さいほど良い値なので、説明変数の適切な組み合わせを試行していきます。 手でAICが高い組み合わせを探していくのも良いですが、これらの組み合わせを自動でおこなってくれる関数があります。 それがStep関数です。 早速Step関数を使ってみます。 > step.lm <- step(dataset.lm) Start: AIC=1090.93 Visit ~ PV + NewerR + Avetime + BR Df Sum of Sq RSS AIC 12404 1090.9
- Avetime 1 467 12871 1099.5
- NewerR 1 805 13209 1107.0
- BR 1 2932 15336 1149.8
- PV 1 40468 52872 1505.0

出ました。通常いくつかの説明変数を試していきますが、今回用意した説明変数を全て使っていった方が適切なモデルであることがわかりました。

以上より、全ての説明変数を使ったモデル式が適切であることがわかったので、回帰モデルは以下となります。

> lm(Visit~.,data=data.set)

Call:
lm(formula = Visit ~ ., data = data.set)

Coefficients:
(Intercept) PV NewerR Avetime BR
-18.3937 0.4023 8.3269 608.0772 22.6680

これらのデータから以下の式が成り立ちます。

Visit = PV ✕ 0.4023 + NewerR ✕ 8.3269 + Avetime ✕ 608.0772 + BR ✕ 22.6680 − 18.3937

これが私のサイトの訪問数を割り出す回帰モデル式になります。
本当にこの式が適合しているかのテストをおこないます。
テストの方法は簡単で、元のデータにモデルの値に代入していくだけです。
結果、日別の状況にはバラつきが出てきますが、全期間の目的変数と実際の訪問数との差異の平均が0.06622となり適合していることがわかりました。
これらの値は当然日別に変化します。また曜日の要素が含まれていないので、曜日別に回帰モデルを設定しても良いかもしれません。日々データを追い、モデル式の調整をおこないながら進めていきましょう。

続いて予測モデルを作っていきます。

予測モデルの方法については、当ブログの「Rで将来のPVをMCMC推定する」が詳しいので詳細はこちらを参照していただくとして、今回は143日分のデータを2分割し、前半(69日分)を予測データとして使い、後半(74日分)を検証データとして使っていきます。

//訪問データをA.Visit変数へ格納
> A.Visit <- dataAB //A.Visit変数の平均値をA.meanへ格納 > A.mean <- colMeans(A.Visit) //A.Visit変数の標準偏差をA.sdへ格納 > A.sd <- sapply(A.Visit,sd) //A.meanをチェック > A.mean
V1
19.98551
//A.sdをチェック
> A.sd
V1
10.47475

//MCMC推定実施
> post.samp <- MCMCmetrop1R(llnormfun, theta.init=c(A.mean,A.sd), + x=x, + thin=1, mcmc=10000, burnin=500, + tune=0.1, + verbose=500, logfun=TRUE) //MCMC推定の結果をグラフ化 > plot(post.samp)

※グラフはこんな感じに。収束期間が長い気がする。さすがにデータが少なすぎたかも。
2013-12-14 16.21.41

//MCMC推定の結果のサマリーチェック
> summary(post.samp)

Iterations = 501:10500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean SD Naive SE Time-series SE
[1,] 12.402 2.778 0.02778 2.211
[2,] 6.048 2.736 0.02736 2.911

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%
var1 9.797 9.962 11.349 14.645 18.31
var2 3.033 3.167 5.634 8.956 10.16

//後半データ(検証用データ)の平均値を参照
> mean(dataAA)
警告: mean() is deprecated.
Use colMeans() or sapply(*, mean) instead.
V1
27.64865
//後半データ(検証用データ)の標準偏差を参照
> sd(dataAA)
警告: sd() is deprecated.
Use sapply(*, sd) instead.
V1
10.22813

う〜ん。わからん。
標準偏差は予測が10.16のところ、10.22813と近づいたけど、平均値が予測18に対して、実際は27ときた。
最近アクセス数が増がしてきていているので予測時点でのデータと差異が生じてきているのかもしれないなぁ。

ひとまずこういう感じです。

最近気になるR本はこれ。
特にRで学ぶデータサイエンスシリーズは難解だけど、相当深く理解できる逸品。
今回紹介するものは、軽量政治分析をRでおこなうもの。オバマや自民党がデータ解析を使って選挙に勝利しましたが、こういった胸熱なこともRでできるって話です。

Leave a Reply

%d人のブロガーが「いいね」をつけました。