チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

研究について,最近思うこと

 研究を行う上で非常に重要なのは,やはりどのような問題に決着をつけるべきなのかその問題を見極めることだと思う。問題を設定するために,はじめは指導教員や先輩の研究者などに教えを請うのが望ましい。それによって,何が研究としてあるべき形なのかを身につけることができると考えられる。
 もちろん,自分で地道に考えていくことで研究関心を見つけていくのもなしではないが,効率が悪いのみならず外した後のリスクが大きくなるように感じられる。また,自分一人で考えたことは自分にこだわりがあまりにも強く出過ぎてしまう可能性もある。
 自分の問題ではなく,その分野の発展に資する問題関心を理解することが必要だろう。ただし,あまりに王道の研究関心は競争相手も非常に多く,マンパワー金のパワーの勝負になりやすい。そういう研究室に所属していて,指導教員の研究の一端を担うならばそれでもいいかもしれないが,個人の研究としてことを進めるのであれば,それはおすすめできない。
 狙うべきはニッチな領域である。ニッチな領域を押さえてそこでてっぺんを取ることが必要。人間はやはり成果が出ないことをやり続けることは難しいが,ニッチなことであれば相対的に成果を出しやすくモチベーションも維持しやすい。モチベーションを維持し続けられるような研究関心にしつつ,成果にも繋がる部分というものを見定めるためには,他人の力も借りて自分の考えを相対化する必要があろう。
 また,ここまで未定義に使ってきた問題という言葉であるが,何が「問題」なのかは注意して考える必要がある。個人的には「今解決されないと多くの人が困ること」を問題とするのがわかりやすいように感じている。(もちろん,2つ以上の集団で決着がついていない根本問題などを扱ってもいいかもしれないが,これを扱うのは必ずしも容易ではない。)もう少しラフに,「今こういうことがわかるとみんながハッピーになること」というような点に注目すると良いかと思う。
ちなみに,「今」というのもポイントだ。ある問題について今解決する必要がないというのは緊急性が低く,人に重要性を伝えられない。ある問題が解決されればよいけど今は別の問題の方が大事,ということもある。問題は一過性のものかもしれないが,「今」だからこそ解決する価値があり,時が立てば意味がなくなる問題も往々にして存在するということには注意したい。
 ちなみに自分の専門の統計学に近づけていえば,統計を使う人の役に立つこと,というのが研究としては大切なのではないかと思っている。新しい手法を提案して,新しい知見を生み出す道具を提供し続けること,今までの道具の不備を解消してより使いやすい方法を提案することなどなどやるべきことがいくらでもある。統計の特定の分野に関していえば,こうした発想の勝負として研究テーマは数多くでてくるように思われる。
 さて,問題を「今解決するとうれしいこと」とみなすと問題設定と同程度に重要な「ストーリー」を組み立てるのが容易になる。論文の導入ではなぜその研究を行う必要があるのか,特定の問題関心についての前提知識を読者と共有し,問題が解決されないとどれほど良くないことがあるのかを示す。それが伝わればあとは妥当は方法論を使って研究を遂行すれば良い。このストーリーが伝わらないと厄介である。
 なぜその研究を行うのか,どう嬉しいのかという説得を行うことができなければ話を理解することができない。その意味で,論文を書くことは一種の説得活動であるともいえる。通常の説得と違うのは単なる感情論での説得は意味をなさず,先行研究を使って論理的に説得するということである。
 論理的な説得を行うためには問題関心に関連した大きなメッセージ,その論文を読んだ人がどのような行動をしてほしいのかといった事がいえなければならない。この大きなメッセージのために,先行研究を積上げ,必然的に疑問に至らざるを得ないように組み立てる。
 どのような「問題」であれば重要とみなされるのか,どのようなストーリーであれば説得的か,論理的飛躍が小さいか,といった点は少なからず感覚的なものにならざるをえない。(論理的な飛躍は,数学などの分野では当然非常に厳密だが。。。)なので,必要なのはそういうことが上手い人の仕事を間近で見ることが最大の学びになる。もちろん,自分で試行錯誤し査読のコメントから学ぶということも多くある。しかし,王道の展開やよくある論理のパタンなどをしっかり身につけるのも感覚がない人には難しい(自分がかなりそう)。
 いつでもどこでも「問題」,「ストーリー」,「説得」ということは頭に入れて研究を行いたいものだ。このあたりについて,もう少し考えが深まって体系的な学習方法や教育方法がわかってきたら,まとまった形にして残しておきたい。

(殴り書き終了)

JAGSで予測分布の作成その2

IRTモデルで予測分布や対数尤度の計算などをおこなってみる。
対数密度の記述は以下の頁の質問に少し書かれている。
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/c728e87b/

The general rule is that if y ~ dfoo(alpha, beta, gamma, ...) then the contribution to the log likelihood of y is given by lik[i] <- logdensity.foo(y, alpha, beta, gamma, ...).

とのことなので,これに従って書いてみよう。(binomial分布はlogdensity.binなので注意)
また,bugsモジュールも読み込んでおく。
(どうでもいいが,T(0,)という関数はサンプリングする範囲を制限するために使うっぽいので,活用すること。)

#
# IRTモデルでの予測分布の作成など。 
#
library("irtoys")
library("jagsUI")
# help(package="jagsUI")
# help(package="irtoys")
#
# 適当な真値
#
Scored2pl

# 適当なデータ
data <- list(x = Scored,
             I = nrow(Scored),
             J = ncol(Scored)
)


# JAGS code 
cat("model{
    # likelihood
    for(i in 1:I){
      for(j in 1:J){
        x[i, j] ~ dbern(prob[i, j])
        logit(prob[i, j]) <- alpha[j] *(theta[i] - delta[j])
        # 予測分布
        x_pred[i, j] ~ dbern(prob[i, j])
        
        # 個々人の個々の項目の尤度
        logLik[i,j]  <- logdensity.bern(x[i,j],prob[i,j])
        #logLik_temp[i,j]  <- logdensity.bern(x[i,j],prob[i,j])
      }
      #logLik[i] <- sum(logLik_temp[i,])

      # prior
      theta[i] ~ dnorm(0.0, 1.0)
    }
    
  # prior
    for(j in 1:J){
      delta[j] ~ dnorm(0 , pr.delta)
      alpha[j] ~ dnorm(0 , pr.alpha)T(0,)
  }
    pr.delta <- pow(100 , -2)
    pr.alpha <- pow(100 , -2)
    }", file="irt_2pl.jag")

fit_jag <- jags(data=data,
                parameters.to.save=c("alpha", "delta","x_pred","logLik","theta"), 
                model.file="irt_2pl.jag",
                modules=c("glm","bugs","basemod"),
                n.chains = 1,
                n.adapt = 0,
                n.iter = 2000,
                n.burnin = 1000,
                n.thin = 2)
fit_jag
Scored2pl$est

JAGSでの予測分布の作成

やり方は簡単で,予測分布を未知の変数としてモデルに書けばいい。
Stanではgenerated quantitiyブロックでrng系の関数を使えばいいが,JAGSではそのまま"x_pred ~ dnorm(mu, sigma)"などとするとRでいうところのx_pred <- rnom(1,mu, sigma)のように乱数を生成してくれる。
以下では母相関係数の推定を行ったあとに,2変量のデータの予測分布を作成してみる。

#
# 予測分布の作成
#

library("jagsUI")

#
# 相関係数については,以下の頁のものを参考にした。
# https://kyougokumakoto.blogspot.jp/2016/09/jags_91.html
#
library(MASS)
set.seed(1234) 
mu <- c(0,0) 
Sig <- matrix(c(1,0.75,0.75,1),2,2) 
x <- mvrnorm(100,mu,Sig) 
cor.test(x[,1],x[,2])

data = list(x = x, N = nrow(x)) 

# JAGS code 
cat("model{
    # likelihood
    for(i in 1:N){
    x[i,1:2] ~ dmnorm(mu[], prec[ , ]) 
    
    # 対数尤度
    logLik[i] <- logdensity.mnorm(x[i,1:2], mu, prec) 
    } 
    
    # 予測分布
    x_pred[1:2] ~ dmnorm(mu[], prec[ , ]) 

    # corresponding the precision matrix
    prec[1:2, 1:2] <- inverse(cov[,]) 
    
    # constructing the covariance matrix
    cov[1,1] <- sigma[1]^2
    cov[1,2] <- sigma[1]*sigma[2]*rho
    cov[2,1] <- sigma[1]*sigma[2]*rho
    cov[2,2] <- sigma[2]^2 
    
    #priors
    rho ~ dunif(-1,1) 
    mu[1] ~ dnorm(0, 1.0e-4) 
    mu[2] ~ dnorm(0, 1.0e-4)
    sigma[1] ~ dt(0, 1, 1)T(0,)
    sigma[2] ~ dt(0, 1, 1)T(0,) 
    
    }", file="cor_pred.jag")

fit_jag <- jags(data=data,
                parameters.to.save=c("mu", "sigma","rho","x_pred","logLik"), 
                model.file="cor_pred.jag",
                n.chains = 1,
                n.adapt = 100,
                n.iter = 2000,
                n.burnin = 1000,
                n.thin = 2)

# 予測分布の相関
cor(fit_jag$samples[[1]][,c("x_pred[1]","x_pred[2]")])
# 予測分布のヒストグラムの重ね合わせ
hist(fit_jag$samples[[1]][,c("x_pred[1]")],freq = F,col=rgb(1,0,1,alpha=0.5))
hist(x[,1],freq = F,add=T,col=rgb(1,0,0,alpha=0.5))

JAGSで対数尤度を得る方法(Stanも)

WAICなどを計算するには,対数尤度のMCMCサンプルが必要。
個々のMCMCサンプルからがんばって計算する方法もあるが,より簡単にはMCMCサンプラーのなかで対数尤度を得ると楽。
Stanではgenerated quantityが容易されているが,JAGSでは少し工夫が必要。

参考
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/a3bef062/
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/87902622/

その他参考
ito-hi.blog.so-net.ne.jp
hamada.hatenablog.jp



JAGSのモデルの中でのようなものでも対数尤度が得られる
一変量の分布に関しては,密度関数などがRのように用意されている。

model{
for(i in 1:I) logLik[i] <- log(dnorm(x[i], mu, 1/sigma))
}

そのほか, logdensity.**系の関数が用意されているが,これはマニュアルに書かれていないので必要に応じて調べる。
多分,すべての分布に対してこれらの関数が利用できる?
(未確認)

#
# 相関係数については,以下の頁のものを参考にした。
# https://kyougokumakoto.blogspot.jp/2016/09/jags_91.html
# 
# logdensity.mnormで多変量正規分布の対数密度を得られる。
# https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/a3bef062/
#
library(MASS)
set.seed(1234) 
mu <- c(0,0) 
Sig <- matrix(c(1,0.75,0.75,1),2,2) 
x <- mvrnorm(100,mu,Sig) 
cor.test(x[,1],x[,2])

data = list(x = x, N = nrow(x)) 

# JAGS code 
cat("model{
# likelihood
for(i in 1:N){
x[i,1:2] ~ dmnorm(mu[], prec[ , ]) 

logLik[i] <- logdensity.mnorm(x[i,1:2], mu, prec) 
} 

# corresponding the precision matrix
prec[1:2, 1:2] <- inverse(cov[,]) 

# constructing the covariance matrix
cov[1,1] <- sigma[1]^2
cov[1,2] <- sigma[1]*sigma[2]*rho
cov[2,1] <- sigma[1]*sigma[2]*rho
cov[2,2] <- sigma[2]^2 

#priors
rho ~ dunif(-1,1) 
mu[1] ~ dnorm(0, 1.0e-4) 
mu[2] ~ dnorm(0, 1.0e-4)
sigma[1] ~ dt(0, 1, 1)T(0,)
sigma[2] ~ dt(0, 1, 1)T(0,) 

}", file="cor.jag")

fit_jag <- jags(data=data,
            parameters.to.save=c("mu", "sigma","rho","logLik"), 
            model.file="cor.jag",
            n.chains = 3,
            n.adapt = 100,
            n.iter = 5000,
            n.burnin = 2500,
            n.thin = 2)

library("loo")

# チェインがたくさんある場合には,長い行列をとしてみる?
log_lik1 <- sapply(1:nrow(x), function(i) do.call("rbind",fit_jag$samples)[,paste0("logLik[",i,"]")])
waic(log_lik1)

log_lik2 <- sapply(1:nrow(x), function(i) fit_jag$samples[[1]][,paste0("logLik[",i,"]")])
waic(log_lik2)


library(rstan)

#
# Stanのコード参考:https://qiita.com/hoxo_m/items/0f1b05681f5d6c4b560a
# ベクトルの配列について注意。y[i]で得られるのはi番目の長さが2のベクトル。
# y[i][1]でi番目の配列にあるベクトルの1番目の要素を取り出すということになる。 
#

cat("data{
 int<lower=0> N;
    vector[2] y[N];
    }
    parameters{
    vector[2] mu;
    real<lower=0> var_x;
    real<lower=0> var_y;
    real cov;
    }
    
    transformed parameters{
    matrix[2,2] sigma;
    sigma[1,1] = var_x;
    sigma[2,1] = cov;
    sigma[1,2] = cov;
    sigma[2,2] = var_y;
    }
    
    model{
    for(i in 1:N){
    target += multi_normal_lpdf(y[i]| mu, sigma);
    //target += normal_lpdf(y[i][2]|0, 1);
    //target += normal_lpdf(y[i][1]|0, 1);
    //y[i] ~ multi_normal(mu, sigma);
    }
    }
    
    
    generated quantities {
    real cor;
    real log_lik[N];
    cor = cov/sqrt(var_x * var_y);
    for(i in 1:N) log_lik[i] = multi_normal_lpdf(y[i] | mu, sigma);
    }

", file="cor.stan")
datastan <- list(N=nrow(x), y = x)
fit <- stan(file = "cor.stan", data=datastan, iter = 5000, chains = 2)

log_lik_stan <- extract_log_lik(fit)
str(log_lik_stan)
waic(log_lik_stan)
waic(log_lik1)
waic(log_lik2)

# 値は大体一致している。

MplusのMCMCサンプルをGetする方法(+Mplusの結果をRでプロットする方法)その1

気が向いたらStanやJAGSなどとの対応関係をきちんと取りならが検算を行っていきたい。


MCMCサンプルは「.gh5」ファイルに入っているので,これを読み込む。

以下の資料のp.65あたりの「mplus.get.bayesian.」を利用すると,ギブスサンプラーMCMCサンプルが得られるはず。
http://www.statmodel.com/mplus-R/Mplus%20R%20tutorial.pdf


以下のinpファイルを作成して,Mplusで推定を行う。
データは
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat
にあるものを,ダウンロードしてくる。

TITLE:	this is an example of an exploratory 
		factor analysis with continuous 
		factor indicators
DATA:		FILE IS ex4.1a.dat;
VARIABLE:	NAMES ARE y1-y12;
ANALYSIS:

TYPE = general;
ESTIMATOR = BAYES;
PROCESSORS = 2;
BITERATIONS = (2000);

model:

y1-y12 with y1-y12;
y1-y12;
[y1-y12];


OUTPUT: CINTERVAL;
plot: type = plot1 plot2 plot3;
library("MplusAutomation")
# MacではRを使ってプロットをする必要がある。
# 以下のsourceを読み込むと必要な関数がいろいろと読み込まれる。
source("https://www.statmodel.com/mplus-R/mplus.R")

# アウトプットの読み込み
hoge1 <- readModels("ex4.1bayes.out")
# 結果の表示
hoge1$parameters

# Rstudioのコンソールでmplusと入力して,タブを押すと色々な関数が候補として出てくる。

# プロットできるもの一覧を表示する。
mplus.view.plots("ex4.1bayes.gh5")
# 変数名を見る
(v_list <- mplus.list.variables("ex4.1bayes.gh5"))
# ヒストグラムの作成
mplus.plot.histogram("ex4.1bayes.gh5","Y1")
# 密度曲線を追加
mplus.plot.densityplot("ex4.1bayes.gh5","Y1")
# qqplotの追加
mplus.plot.qqnorm("ex4.1bayes.gh5","Y1")
# 散布図の作成
mplus.plot.scatterplot("ex4.1bayes.gh5","Y1","Y2")
# データの取得:1変数だけ。
mplus.get.data("ex4.1bayes.gh5","Y1")

ベイズ推定の結果をプロットするためのもの。
まだ利用しきれていない関数もあるので,今後追記する。

# 
# ベイズ推定の結果をいろいろと得るためのもの
# mplus.get.bayesian.●●
#
# パラメタ名を取得する
# ついでに,パラメタ名が必要な関数も表示してくれる。
(par_list <- mplus.list.bayesian.parameters("ex4.1bayes.gh5"))
#
# 特定のパラメタのMCMCサンプルを取得する。
# 列がチェインに対応している。
#
par_data <- mplus.get.bayesian.parameter.data("ex4.1bayes.gh5", par_list[1])
# トレースプロット
matplot(par_data, type="l")
mplus.plot.bayesian.traceplot("ex4.1bayes.gh5",paramstr = par_list[1])
# 自己相関
mplus.plot.bayesian.autocorrelation("ex4.1bayes.gh5",paramstr = par_list[1])
# 事後分布
mplus.plot.bayesian.distribution("ex4.1bayes.gh5",paramstr = par_list[1])


Mplusは便利なので,もう少し使いこなしていきたいところ。

よく必要になるけど,よく忘れるR(とExcel)の小技メモ

文字(数字)の桁数を揃えたい。
参考URL:
http://d.hatena.ne.jp/teramonagi/20100902/1283429021

formatC("hoge", widh = 5, flag="0")

flagは '0 + - #' だけが指定可能とのこと。


不要な空白を消したい。

# いらない空白を全て削除する
# 全角半角の指定に注意。
gsub(" ", "", x, fixed = TRUE)
# 文字列の前後の空白を削除する
stringr::str_trim(x, side="both")

全角半角を変換したい。
Nipponパッケージを利用するといいらしいが,上手く動作しない場合がある。
そんなときは,Excelを使う。
参考URL:全角を半角に、半角を全角に変換する エクセルの関数技

文字を半角にするにはASC関数を使い,半角文字を全角にするにはJIS関数を使う。

不要な全角半角スペースを削除するには,=SUBSTITUTE(SUBSTITUTE(A1," ","")," ","")という関数を利用するとよい。
(trimは文字列前後の空白を削除するための関数。)


Rで図を出力する場合のあれこれ。
参考URL
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html
R のグラフで日本語を使う - 知に至る病
https://oku.edu.mie-u.ac.jp/~okumura/stat/first.html
Rのグラフで軸の目盛りの刻み幅を変更する方法 - Rプログラミングの小ネタ

# 文字のフォントを指定する。
# これを指定しておけば,とりあえず文字化けしない。
pdf(file="hoge.pdf",  family="Japan1GothicBBB")

# これはMac向き
par(family="HiraKakuProN-W3")  # ヒラギノ角ゴシックProN W3

# 空のプロットを作成する。
 plot.new()

# tclで軸を内側に向ける。
# 正の値が内側で,枠の大きさに対する割合で指定で指定する。0.3ぐらいがけっこうちょうどいい。
plot(0,1,tcl = 0.3)

# メモリを独自に振り直したいときに使う。
#メモリを消す
curve(dnorm, xlim=c(-3, 3), xaxt="n") 
# atに書き込みたい目盛りのベクトルを指定する
axis(side=1, at=-3:3) 
# sideは下から時計回りに1,2,3,4になっている。

# ラベルを各軸に並行して描く
las = 0
# ラベルをすべて水平に描く
las = 1
# ラベルを軸に対して垂直に描く
las = 2
# ラベルをすべて垂直に描く
las = 3

Mplusでターゲット行列を指定した回転を行う方法

Analysisコマンドで,rotation = targetを指定し,
modelコマンド内で目標とする値を~で指定する。
また,各因子の最後の(*t)などとつけて,各因子負荷が1つの因子負荷行列で表すことができることを示さないといけない。
Mplusのマニュアルは,versin 8のp681などに記載されている。正直,いろいろと回転がありすぎて何が何だかわからない。
本文をtargetで検索するのがはやい。
これを使えば,任意の値を指定したターゲット行列に向けた回転を行ってくれる。
(つまり,プロクラステス(procrustes)回転ができる。)
確認的因子分析よりもゆるく,探索的因子分析よりも仮説を反映した分析ができるので,もう少し利用されてもいい方法のように感じる。


使用したデータ
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat

TITLE:
DATA:
FILE IS ex4.1a.dat;
VARIABLE:
NAMES ARE y1-y12;
ANALYSIS: 
ESTIMATOR IS ML;
 ROTATION=TARGET;
 MODEL:
f1 BY y1-y12 y1~0 y2~0 y3~0 y4~0 y5~1 y6~1(*t);
f2 BY y1-y12 y7~0 y8~0 y9~0 y10~0 y11~1 y12~1(*t);
  OUTPUT:
MODINDICES CINTERVAL;
  

参考PDF
http://www.restore.ac.uk/appliedpsychometrics/Tutorial.materials/Introduction.to.Mplus/Introduction%20to%20Mplus%20session%202.pdf