bayesmを使った階層ベイズモデル（線形回帰）

通常の重回帰分析を階層ベイズモデルに拡張するには、Rのbayesmを使うのが簡単です。
bayesmは、いろいろな階層モデルが用意されていますが、重回帰分析のモデルは
rhierLinearModel
となります。

MCMCの計算自体はとても簡単です。

階層モデルのパラメータDelta.meanと観測モデルのパラメータbeta.meanの行列計算が間違えるポイントかもしれません。
パラメータが2次元の行列なのと、さらに、MCMCごとのデータがあるので、3次元の構造となっています。
どこの軸を集計すると、良いのか、慣れるまで時間がかかるかもしれません。

[観測モデル]
y=bo+b1*x1+b2*x2+b3*x3+ε

[階層モデル]
b0=a1*z1+a2*z2+ε
b1=a1*z1+a2*z2+ε
b2=a1*z1+a2*z2+ε

簡単な使い方の例

library(bayesm)
R=20000

# 個人数
nreg <- length(member_id) # num of id

# Zは、階層モデル用のデータ
nz <- ncol(Z) # カラム数

# 観測モデル用のデータ
regdata=NULL
for (reg in 1:nreg) {
regdata[[reg]] <- list(y=y, X=X)
}

# MCMC
Data1 = list(regdata=regdata, Z=Z)
Mcmc1 = list(R=R, keep=1)
Mcmc_out = rhierLinearModel(Data=Data1, Mcmc=Mcmc1)

# Burn-In の設定
mcmc_start <- R/2
mcmc_end <- R

# パラメータの事後平均のデビアンス
for(j in 1:nvar_para){
for(i in 1:nreg){beta.mean[i,j] <- mean(Mcmc_out$betadraw[i, j, (mcmc_start:mcmc_end)])}
}

## 階層モデルのパラメータ
Delta.mean <- matrix(0, nz, nvar_para)

for(j in 1:nvar_para) {
for(i in 1:nz) {
Delta.mean[i,j] <- mean(Mcmc_out$Deltadraw[(mcmc_start:mcmc_end), i+(j-1)*nz])
}
}
Delta.mean

階層ベイズモデルの切片項の扱い

今、下記の様なモデルを考えたとします。

[観測モデル]
y=bo+b1*x1+b2*x2+b3*x3+b4*x4+ε

[階層モデル]
b0=a0+a1*z1+a2*z2+a3*z3+ε
b1=a0+a1*z1+a2*z2+a3*z3+ε
b2=a0+a1*z1+a2*z2+a3*z3+ε
b3=a0+a1*z1+a2*z2+a3*z3+ε
b4=a0+a1*z1+a2*z2+a3*z3+ε

ここで、zは、仮に血液型として、
z1：O型
z2：A型
z3：B型
z4：AB型
とします。

当初、観測モデルにも観測モデルにも切片項を入れていましたが、微妙にしっくりこない結果となります。
どうやら、両方に切片項が入っているとうまく収束しない問題があるらしい。

そこで、切片項は片方だけに入れるモデルに修正して、代わりに、z4を入れるモデルに変更しました。
通常は、観測部分に切片項を入れて、階層部分は入れない。

[観測モデル]
y=bo+b1*x1+b2*x2+b3*x3+b4*x4+ε

[階層モデル]
階層モデル
b0=a1*z1+a2*z2+a3*z3+a4*z4+ε
b1=a1*z1+a2*z2+a3*z3+a4*z4+ε
b2=a1*z1+a2*z2+a3*z3+a4*z4+ε
b3=a1*z1+a2*z2+a3*z3+a4*z4+ε
b4=a1*z1+a2*z2+a3*z3+a4*z4+ε

↑
こちらの修正したモデルの方が、感覚的にしっくりくるパラメータ推定となっていました。

jupyter notebookを使って機械学習のモデルを作成したい

ようやく、jupyter notebookで機械学習をするための環境が整いました。

初めてのPython入門みたいな書物は多いですが、入門すぎて読んでも使えなかった。
いくつか、つまずいたポイントを備忘録としてメモ。

1. データの準備
練習用のお試しだと、ローカルのcsvファイルを読み込んで、それを分析してという流れになりますが、企業のデータは、csvファイルを使っているわけではない。
おそらく、多くのデータは、データベース（google big queryなど）に入っています。
そこで、データベースから、ローカルマシンにデータを持ってくる必要があるのですが、セキュリティ認証の設定とか、そもそものPythonコードの書き方をどうするの？問題がありました。

2. ディレクトリの構成
sqlをPythonの中にベタ書きするのではなく、sqlなどのフォルダを作ってそこに書いておく。
そして、Pythonからクエリを読み込み、実行してデータを取得する、という自然な流れなのですが、ここも、少々面倒くさい。

MacにBigQueryAPIクライアントライブラリをインストール その4 pandas-gbqをインストールする
https://skellington.blog.ss-blog.jp/2021-01-28
↑
query_1 = read_query('/sqls/test.sql')
read_queryでsqlを変数に取り込んで、

data_frame = pd.read_gbq(query_1, PROJECT_ID)
read_gbqで実行する。

3. クラス定数の設定とか、初期化スクリプトの実行とか
このあたりの設定は、慣れると簡単ですが、最初は、何をどう設定すれば良いのか、？？？状態でした。ｗ

新型コロナウィルスの抗体保有率 東京都は0．91％

東京の抗体保有率が発表されて、0.91%でした。
イギリス・ロンドンだと20%くらいで、インド・デリーは56%を超えているとか。
日本は、まだまだ少ないんですね。

新型コロナウィルスは、かれこれ1年に経ち、いろいろなことがわかってきたと思います。
いくつか整理すると
・東南アジアの人は、感染する確率が欧米に比べて1/10程度と低い。
・感染すると、そのあとの重篤化率や致死率はどこも同じ。
・ワクチンは重症化するのを抑えたり、感染を抑えたりする効果はある程度あるが0にはならない。
・人が密になると、感染する率はあがる。
・暴飲暴食を避けて、しっかり寝ることが大切。
などなど。

また、抗体は数ヶ月で消えるけど、一度感染していると、B細胞やT細胞が攻撃パターンを覚えていて（その数、無量大数らしい）重症化を避けることができる。

ここから考えられることとして、今、日本は抗体保有率が低いので、経済を再開すると、また、感染者数がわ〜っと増えると思います。
ただ、東南アジアの人は重症化しにくいので、そこまで死者は増えない（基礎疾患がある人を除く）。
そして、これらを繰り返し、数年後には、ほとんどの人が一度は感染するウイルスになると考えられます。
RSウイルスとかヘルペスウイルスとか、ほとんどの人が感染している状態。

おそらく、重症化しにくい子供のうちに感染しておくと、B細胞に記憶され、ただの風邪になる。

東南アジアの人が感染しにくい理由としては、過去によく似たウィルスが既に入ってきており抗体ができているとか遺伝子的に強いとか、いろいろな理由が考えられますが、実際どうなのかは、これからの研究を待たないとわからないです。

首相 高所得者に理解求める 児童手当廃止問題

1200万円以上が廃止になるとのこと。
正直、呆れてしまった。。。

高所得者は何かと税金取られます。
取りやすいところから取りましょう、ってことでしょうかね。
選挙でも高所得者の割合は少ないので、影響は少なそうだし。
一方、公平にというのであれば、現役の高所得者世代だければなく、高齢世代も含めて広く負担を分かち合った方が良い。
でも、選挙のことを考えたら、高齢者の負担は言えないんだろうな。。。

高所得者といっても、何もせずに高所得者になったわけでない。
多くの人は人よりも努力をして、人よりも多く働いて、人よりも多くリスクと取り、人よりも多く責任（ストレス）があることをして、年収を上げてきたんだと思います。

今回だけでなく、これまで安直に高所得者から税金を取る政策ばかりでした。
頑張った人が働く意欲をなくす政策より、頑張ろうとする政策ってできないんだろうか？

また、少子化というのは、みんなが考えている以上に深刻だと思う。
新型コロナウィルスの影響で少子化がさらに加速しています。
よくわからない対策や、経済を止めたことによる損失を考えると、今回の高所得者から得られる税金はわずかなもの。

どこにお金を使うべきか、今一度考えて欲しいものです。

【モンスト】ワールドトリガーコラボ、太刀川慶【超究極】

ワールドトリガーコラボですが、今回の超究極は2つ。
・ヴィザ
・太刀川慶

ヴィザは、いったん後回しで、太刀川慶を周回していました。
その他3体でクリアというミッションがあるのですが、アルセーヌ（アナーキー/獣神化）が強すぎるので、後は、なんでも良い気がします。w

挑んだPTは、こんな感じです。
1.光ゼレフ（神化）
2.志々雄真実（進化）
3.オールフォーワン（進化）
4. アルセーヌ（アナーキー/獣神化）

後で、攻略をみたら、
カレン（進化）
も良さそうです。
オールフォーワン（進化）の代わりにカレン（進化）が良いかもしれません。

【映画】銀魂 THE FINAL

「銀魂 THE FINAL」の上映時間を観てきました。

銀魂 THE FINAL
https://wwws.warnerbros.co.jp/gintamamovie/

相変わらず、パロディー全開で、お笑いが中心の内容ですが、師との迫力ある戦闘シーンもなかなかの見所です。
ネタバレしないように、感想を書くと、良い終わり方だったんじゃないかなって思います。

最後には、どこかで撮った舞台挨拶も流れていました。
せっかくなので、そこまで観ておきたいかなと。

本当にfinalかどうかはわかりませんが・・・w
でも、これで終わりってのも寂しいですね。

統計の「超」勉強法

嫁が、daigoの超勉強法を動画を見て、子供の受験対策で何か活用できないか、と話題になった。
「超」勉強法の話を聞いていると、自分が統計を身につけた体験と結構似ていて、知らず知らずに自分も統計を「超」勉強法っぽい感じで学習していたようです。

せっかくなので、統計の「超」勉強法っぽいことを書きたいと思います。

まずは、基本原則から。
1. 「面白いことを勉強する」
統計を勉強する際に、推定と検定が割と教科書の最初に出てくるのですが、一番ややこしい部分です。
なので、さらっと流して、モデリング（回帰）から始めるのが良いかと思います。
中でも、自分が必要としているモデリングからやるのが良い。
面白いというよりかは、必要に迫られているからものすごく集中して身につけようとするから。

2.「全体から理解する」
統計の細かい数式を理解しようとすると、挫折します。
ただ、細かい数式を理解できると幅が出てくるので、最初はあまり細部を見なくて、何度か繰り返し学習していく中で少しずつ細部を理解していくのが良いかと思います。

全体を理解するってことですが、まずは感覚的に理解すること。
数式がわからない人にどうやったら理解してもらえるか？を考えると、「要は、こういうことなんだよね。」という例え話で理解すると良いかもしれません。

3.「8割までをやる」
よくある単元を完璧に理解するまで、先に進めない人がいます。
なかなか100%ってのはないので、7割から8割くらいで次に進むのが良いかと思います。


########################################
ここから先は、経験的に、自分が思う学習法として載せておきます。
########################################

4. 人に教える
人に教える、つまり、口に出して話すことで、意外とうろ覚えだったり、話している中で「どういう意味だっけ？」となります。
その部分について、徹底的に調べて理解する。
2,3の逆ですが、一つ一つ潰していくことで、7割〜8割を9割、10割りに近づけていく。
これをやらないと、いつまで経っても、8割のまま。

また、人に教えている中で、必ず質問がきます。
この質問が大切で、自分が気がつかなかった視点を教えている人から教えてもらうことができます。
当然、すぐにわからない部分もたくさんあるので、ここに関しても、再学習です。

5. 忘れたら覚える。
忘れる前に覚え直すことはあまり良くないらしい。
忘れてしまった！という後悔が記憶にいいらしいので、忘れたら覚え直すを繰り返すと頭に定着します。

6. 別の角度から見つめ直す。
同じ教科書からだと、同じ側からしか見えない。
いろいろな人の講義を聞いたり、いろいろな書物を読むことで、平面だったものが立体に見えてくる。


########################################
さらに上級者になるために
########################################

7. 真理の扉の先に真理の扉がある
もう、この部分に関しては、相当知識がついて、ほぼ完璧にわかっている！と思っていても、実は、まだまだ、自分はわかっていなかったって気がつく部分があります。

これが最後の真理の扉かなと思って開いたら、さらに、その奥に広い部屋があって、まだ、次の扉を開いていく。
そんな感覚です。

1.〜6.を繰り返していくと、世の中的にも相当できる人になっているかと思います。
そこで、驕れることなく、まだ、この先に、自分が知らない世界があるはずだという気持ちを忘れなければ、きっと、扉は向こうからやってくることでしょう。

「これを知る者はこれを好む者に如かず。これを好む者はこれを楽しむ者に如かず。」ってやつですかね。

最短の時間で最大の成果を手に入れる　超効率勉強法

• 作者: メンタリストＤａｉＧｏ
• 出版社/メーカー: 学研プラス
• 発売日: 2019/03/05
• メディア: Kindle版

appleの学割キャンペーン2021、開始！

初売りセール2021が終わってすぐですが、appleの学割キャンペーン2021が始まりました。
開催期間：2021年 2月4日(木) 〜 4月7日(水)

自分は、おもいっきり学生なので、学割が使えるわけですが、さて、何を買おうかと。
初売りと違って、じっくりと考える余裕はあります。w

初売りは、iPhoneとかありましたが、学割はMacbookとかipadが中心。
ipadに至っては、家に3台くらいあるので、さすがに不要。
PCも、MacPCだけで、5台くらい転がっています。
何かを捨てないと、買えなさそうですが・・・。

個人的には、M1チップ (Apple Silicon)搭載マシンが欲しいので、選択はかなり限られている気がします。。。w

対応分析（コレスポンデンス分析）の結果をバブルチャートで作図したい

統計の講師をしていて、なるほど！と思う質問を受けることがあります。
せっかくなので、その中からピックアップして紹介できればと思います。

【質問】
対応分析（コレスポンデンス分析）の結果をバブルチャートで作図したいけど、どうすれば良いでしょうか？

【回答】
まずは、対応分析の結果を出します。

library(MASS)
csv.movie <- "./data/movie_count.csv"
movie <- read.csv(file=csv.movie, header=TRUE, sep=",")

ca.res <- corresp(movie.x, nf=2)
ca.res

biplot(ca.res, cex=0.6)

対応分析の結果をca.resに入れておきます。

# バブルの大きさの定義
c <- colSums(movie.x)
r <- rowSums(movie.x)

# バブルプロットの描画
symbols(ca.res$cscore[, 1], ca.res$cscore[, 2], circle = c, inches = 0.4, fg = "white", bg = "lightblue", xlab = "Dim 1", ylab = "Dim 2")
text(ca.res$cscore[, 1], ca.res$cscore[, 2], colnames(movie.x))
par(new=T)
symbols(ca.res$rscore[, 1], ca.res$rscore[, 2], circle = r, inches = 0.4, fg = "white", bg = "pink", xlab = "Dim 1", ylab = "Dim 2")
text(ca.res$rscore[, 1], ca.res$rscore[, 2], rownames(movie.x))

R stanのエラー： コネクションを開くことができません

統計の講師をしていて、なるほど！と思う質問を受けることがあります。
せっかくなので、その中からピックアップして紹介できればと思います。

【質問】
R stanでエラーが出ますが、どうすれば良いですか？

Error in file(fname, "rt") : コネクションを開くことができません
追加情報: 警告メッセージ:
1: normalizePath(file) で:
path[1]="./stan/Binom.stan": No such file or directory
2: file(fname, "rt") で:
ファイル './stan/Binom.stan' を開くことができません: No such file or directory
get_model_strcode(file, model_code) でエラー:
cannot open model file "./stan/Binom.stan"

【回答】
stanを動かすには、なかなかプログラムを書くという難しがあります。
データを入れて、コマンドを書けば、それっぽい結果が出てくるというその辺りのツールとは違うので、取り扱うにはそれなりのスキルと、理論的な背景の理解が必要だったりします。

Error in file(fname, "rt") : コネクションを開くことができません
とか
No such file or directory
というメッセージをみると、多くは、ファイルの場所の問題です。

作業ディレクトリは、getwd()で確認できます。
上記の場合、作業ディレクトリの配下に、stanというフォルダがあり、その中に、Binom.stanというstanファイルがある前提です。

エラーの多くは、全然別の場所を指定している場合がほとんどなので、
setwd(" 作業ディレクトリの場所 ") で変更してあげると良いでしょう。

個人的な感覚ですが、エラーで動きませんと言われて、けっこうな確率で、このファイルの場所問題である場合が多いです。

傾向スコアを使ったアップリフトモデルについて考える

ABテストみたいにランダム実験ができる環境だと良いのですが、例えば、薬を投与する/投与しないといったことをランダムでやるには人道的に反しますし、ビジネスにおいてもあえて効果がない群にキャンペーンを実施することはNGとなる場合が多いです。

そこで傾向スコアみたいな考え方が出てきて、あえてランダム実験をしていなくてもアップリフトが計算したいというニーズはありそうです。

傾向スコアでもでりんぐしてみるとこんな感じです。

知りたいのは、
目的変数: re78
介入効果: treat
の効果が知りたい。

# Propensity Score
install.packages("Matching")
library(Matching)
data(lalonde)
head(lalonde)

# どんな人に介入をしたのか（treat）をロジスティック回帰分析で傾向をつかんでおきます。
model.glm <- glm(treat ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74+ u75,
family=binomial, data=lalonde)
summary(model.glm)

# 傾向スコアを使わない場合
model2 <- lm(re78~treat, data=lalonde)
summary(model2)

## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4554.8 408.0 11.162 < 2e-16 ***
## treat 1794.3 632.9 2.835 0.00479 **
##
## 介入効果は、1794.3
## ただし、他の変数の影響を受けている

# 傾向スコアを使う場合
model3 <- Match(Y=lalonde$re78, Tr=lalonde$treat, X = model.glm$fitted)
summary(model3)

## 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

2138.6というのがtreatの効果だと言っています。

ここで、ふと気がついたのですが・・・
マーケティングにおけるキャンペーンって、人によって、また、タイミングによって施策の効果が変わってくるはずです。
知りたいのは、平均的な消費者像に対するキャンペーン効果ではないはず。

となると、傾向スコアを使ったアップリフトモデルというのは、ちょっと厳しそう。
実は、ここの異質性を考慮した傾向スコアというのも出せるのかもしれませんが、ちょっと研究不足です。m(_ _)m

R stanを実行する方法 モデルを読み込む2つの方法

stanを実行する際に必要なものが2つあります。
・データ
・stanのモデル

stanのモデルを読み込む方法として、
- 外部にテキストで保存してそれを読み込む方法
- Rでモデルを書いてそれを読み込む方法
があります。

?stan
をRで実行するとヘルプが出てくるのですが、それぞれ引数が書かれています。

Usage
stan(file, model_name = "anon_model", model_code = "", fit = NA, 
  data = list(), pars = NA,
  chains = 4, iter = 2000, warmup = floor(iter/2), thin = 1, 
  init = "random", seed = sample.int(.Machine$integer.max, 1), 
  algorithm = c("NUTS", "HMC", "Fixed_param"), 
  control = NULL, sample_file = NULL, diagnostic_file = NULL, 
  save_dso = TRUE, verbose = FALSE, include = TRUE,
  cores = getOption("mc.cores", 1L),
  open_progress = interactive() && !isatty(stdout()) &&
                  !identical(Sys.getenv("RSTUDIO"), "1"),
  ...,
  boost_lib = NULL, eigen_lib = NULL
  )

file
The path to the Stan program to use. file should be a character string file name or a connection that R supports containing the text of a model specification in the Stan modeling language.
A model may also be specified directly as a character string using the model_code argument, but we recommend always putting Stan programs in separate files with a .stan extension.
The stan function can also use the Stan program from an existing stanfit object via the fit argument. When fit is specified, the file argument is ignored.

model_code
A character string either containing the model definition or the name of a character string object in the workspace. This argument is used only if arguments file and fit are not specified.

- 外部にテキストで保存してそれを読み込む方法
file= 引数を使う
res_mcmc_bin <- stan(file="./stan/Binom.stan", data=data_stan_bin)

- Rでモデルを書いてそれを読み込む方法
model_code=　引数を使う
stan_code <- '
data{
int n;
int x;
}

parameters{
real p; } model{ x ~ binomial(n, p); } ' res_mcmc_bin <- stan(model_code=stan_code, data=data_stan_bin) あらかじめRでモデルを書く方法ですが、複数行のコードなので、シングルクォートでくくります。 stan_code <- ' ***** ***** ***** '

二項ロジスティック回帰分析の結果が全て同じ結果になる！？

統計の講師をしていて、なるほど！と思う質問を受けることがあります。
せっかくなので、その中からピックアップして紹介できればと思います。

【質問】
二項ロジスティック回帰分析を行ったところ、全ての説明変数で同じ数値（係数やp値）になってしまいました。

【回答】
この場合は、多重共線性が疑われます。

### vifで多重共線性を確認
library(car)
model <- glm(y~., family = binomial, data=dat)
vif(model)

### AICによる変数選択
library(MASS)
## 増減法"both"
modelAIC <- stepAIC(model, direction="both")

傾向スコアを使ったアップリフトモデルについて考える その2

傾向スコアを使ったアップリフトモデルについて考える
https://skellington.blog.ss-blog.jp/2021-02-12
↑
こちらの続き。

まずは、前回のおさらいから。

傾向スコアの手順
1. データを3つに分けて考えます。
目的変数:y
説明変数その1:CTL or TG
説明変数その2:他の変数X（変数Xの中にyや"説明変数その1"は入れない）

2. CTLとTGを他の変数Xを使ってロジスティック回帰分析で変数XがCTL/TGにあたえる影響をモデリングする。

3. 目的変数yとCTL/TGの傾向スコアをみて、yとCTL/TGの関係をモデル化する。

アップリフトモデルの思想としては、
a. CTL と TGの効果を正しく見積もりたい
b. 人によって、施策効果（TGの効果）は違うはずなので、効率の良いセグメントを見つけたい
です。

しかし、この方法だと、a.の効果は分かっても、b.の効果はわからない。
CTL/TGの効果は、平均的な消費者像を扱っていることになります。

b.に関して、人ごとの効果のモデリングはは難しいとして、ある程度のセグメントごとにその傾向をみたいと考えるわけです。

簡単なセグメントの場合（一般顧客、優良顧客くらい）だと、データを一般と優良の2つに分けることで、それぞれのCTL/TGの効果をみることはできそうです。
しかし、セグメントをものすごく細かくしていった場合どうなるか？

完全にランダムデータの場合（CTLとTGがランダムに選ばれている場合）は、わざわざ傾向スコアを持ち出さなくても、簡単に分析できます。
マーケティングの施策（CTLとTG）は、何かの戦略によって、TGになったりCTLになったりしているので、セグメントを細かく指定いくと、ほとんどCTLだったりTGだったりとなる可能性出てきます。
このような状況の場合は、各セグメント後の施策効果を正しく判断するのはできなくなります。

いっけん遠回りに思えるかもしれないが、最初にランダムな実験をしておいて、質の良いトレーニングデータを作っておくことで、その後のモデル化が効率化されたり、各セグメントの効果を正しく見積もれるようになると考えています。

【モンスト】ヴィザ【超究極】〜サイドエフェクトなしでクリア

サイドエフェクトありだと簡単にクリアできるのですが、なしだとなかなかクリアできない。
実の厳選とか紋章適当とかだったので、仕方ないですかね。

最終日ということで、一度、盛り盛りに持った状態でやってみたところ、2回目で勝つことができました。

基本的な戦略はこんな感じ。

1.ブリュンヒルデ（獣神化）
2.ブリュンヒルデ（獣神化）
3.爆豪勝己（獣神化改）
4.空閑遊真（獣神化）

ボスの2戦目で、爆豪勝己のssを打ってから、空閑遊真のssで大幅に削れます。
ゲージ飛ばすことできなかったけど、3戦目はかなり楽になります。

ボス1戦目もssを使いたいところだけど、ここは我慢。
空閑遊真が2体いたらもっと楽だったかもしれません。

クリアしてもらえるのは、称号くらいなんですがね・・・w