一般化線形混合モデル(GLMM) -Poisson/NormalをMCMCで書き直す
GLMM(Poisson - Normal)をMCMCに!
つい、10時間ほど前に投稿したGLM(ポアソン回帰モデル)のmcmc版の作成に続けて、GLMMのポアソン/正規モデルもmcmcで書いてみました。(10時間前の記事がこちら)
こちらは、とにかく変量効果の提案分布を構成するのにとにかく手こずりました。個人的に、truncated normalが一番良さそうだと思っていますが、適当な(よく使われる)提案分布って何かあれば教えてください。(log-normal / gamma)とか的外れなのに、いろいろ試して時間を無駄にして、あげく、やはりthresholdがある提案分布が楽だなーという安易な発想のもと、、、ここに落ち着いております。
使っているアルゴリズムは、MHアルゴリズムです.相変わらず、紙に書き直して、計算をできる限り速くできるように、logを取ったり、うまいこと変数を消したりして、ある程度の効率化は図ってみました。下の方に、残骸としてはっつけております(笑。
追記(5月15日)
以下のコードには、間違いがありますので、コードについてはこちらの記事をご覧ください.tomoshige-n.hatenablog.com
#poisson regression set.seed(20150505) library(dummies) library(truncnorm) #simulation x = data.frame(x1 = rnorm(100,0,1), x2 = rnorm(100,0,1)) x = as.matrix(x) z_fac = c(rep("A",30),rep("B",25),rep("C",45)) z = dummy(as.factor(z_fac)) beta = c(0.8,0.6) u = c(rnorm(30,0,0.5),rnorm(25,0,0.2),rnorm(45,0,0.8)) lambda = exp(x%*%beta + u) #z%*%u は正規分布からの乱数と考えればよい。 y = rpois(100,lambda) par(mfcol=c(1,2)) plot(x[,1],y,xlab="",ylab="",col=as.factor(z_fac)) plot(x[,2],y,xlab="",ylab="",col=as.factor(z_fac)) betaprop = function(beta_before,step_width,sigma){ p = dim(x)[2] q = dim(z)[2] n = length(beta_before) mu = c(beta_before[seq(1,p,1)],rep(0,q)) Cov = diag(c(step_width,step_width,1/sigma),p+q) sample = mvrnorm(1, mu = mu, Sigma = Cov) return(sample) } sigmaprop = function(sigma,sd_trunc){ q = dim(z)[2] mu = sigma sample = rep(0,q) for(i in 1:q){ sample[i] = rtruncnorm(1, a=0, b=Inf, mean = mu[i], sd = sd_trunc) } return(sample) } prop_lograte = function(beta_after,beta_before,sigma_after,sigma_before,p,q,sd_trunc){ beta_log = sum( (1/2)*log(sigma_after) - (sigma_after/2)*((beta_after[seq(p+1,p+q,1)])**2) ) sigma_log = sum( (1/2)*((sigma_after-sigma_before)/sd_trunc) - log(1 - pnorm(-sigma_before/sd_trunc)) ) return(beta_log + sigma_log) } #posterior prior_logsigma = function(sigma,theta){ result = sum((k-1)*log(sigma) - (1/theta)*(sigma)**2) return(result) } #mixed effectの所のlogは計算する必要がある。 prior_logmixbeta = function(beta,sigma){ #betaとsigmaは次元分あるので、最後にsumを取る result = sum(-(1/2)*((beta**2)*sigma) + (1/2)*log(sigma)) return(result) } log_post = function(y,x,z,beta,sigma,theta){ p = dim(x)[2] q = dim(z)[2] beta1 = beta[seq(1,p,1)] beta2 = beta[seq(p+1,p+q,1)] log_lambda = x%*%beta1 + z%*%beta2 likeli = sum(log_lambda*y) - sum(exp(log_lambda)) #prior calculation #beta|sigma beta_given_sigma_prior = prior_logmixbeta(beta2,sigma) #sigma prior_logsigma = prior_logsigma(sigma,theta) result = likeli + beta_given_sigma_prior + prior_logsigma return(result) } mcmc_pois = function(y,x,z,intercept = F,iter=1000,burn_in = 500,step_width=0.01,theta=0.01,sd_trunc=0.5){ p = dim(x)[2] q = dim(z)[2] initial_beta = c(rnorm(p),0,0,0) initial_sigma = rgamma(q,100,1) beta.r = matrix(0,nrow=p+q,ncol=iter) sigma.r = matrix(0,nrow=q,ncol=iter) log_likeli = rep(0,iter) accept = 0 #初期値の設定 beta.r[,1] = initial_beta sigma.r[,1] = initial_sigma log_likeli[1] = log_post(y,x,z,initial_beta,initial_sigma,theta) alpha.r = rep(0,iter) for(i in 2:iter){ #sigmaの候補を生成 candidate_sigma = sigmaprop(sigma.r[,i-1],sd_trunc) #betaの候補を生成 candidate_beta = betaprop(beta.r[,i-1],step_width,candidate_sigma) #この値に基づくlog-postを計算する log_likeli_candidate = log_post(y,x,z,candidate_beta,candidate_sigma,theta) #棄却率alphaの計算 alpha = exp( log_likeli_candidate - log_likeli[i-1] - prop_lograte(candidate_beta,beta_before=beta.r[,i-1],candidate_sigma,sigma.r[i-1],p,q,sd_trunc) + prop_lograte(beta.r[,i-1],beta_before=candidate_beta,sigma.r[i-1],candidate_sigma,p,q,sd_trunc) ) #alphaに基づく条件分岐 alpha.r[i] = alpha if(runif(1) < alpha){ beta.r[,i] = candidate_beta sigma.r[,i] = candidate_sigma log_likeli[i] = log_likeli_candidate accept = accept + 1 } else{ beta.r[,i] = beta.r[,i-1] sigma.r[,i] = sigma.r[,i-1] log_likeli[i] = log_likeli[i-1] } } print(paste("accept rate is ",accept,sep="")) return(list( sample_beta = beta.r, burn_in_beta = beta.r[,-c(1:burn_in)], sample_sigma = sigma.r, burn_in_sigma = sigma.r[,-c(1:burn_in)], alpha=alpha.r, ll = log_likeli)) } result = mcmc_pois(y,x,z,intercept = F,iter=30000,burn_in = 10000,step_width=0.3,theta=0.7,sd_trunc=0.4) #beta posterior par(mfcol=c(3,4)) hist(result$burn_in_beta[1,]) plot(result$burn_in_beta[1,],type="l") hist(result$burn_in_beta[2,]) plot(result$burn_in_beta[2,],type="l") hist(result$burn_in_beta[3,]) plot(result$burn_in_beta[3,],type="l") hist(result$burn_in_beta[4,]) plot(result$burn_in_beta[4,],type="l") hist(result$burn_in_beta[5,]) plot(result$burn_in_beta[5,],type="l") #sigma posterior par(mfcol=c(2,3)) hist(1/result$burn_in_sigma[1,]) plot(1/result$burn_in_sigma[1,],type="l") hist(1/result$burn_in_sigma[2,]) plot(1/result$burn_in_sigma[2,],type="l") hist(1/result$burn_in_sigma[3,]) plot(1/result$burn_in_sigma[3,],type="l")
追記
ハイパーパラメータの推定法は、まだ勉強できていないので、その辺りは今後の課題です。
次は、変数間に相関でも入れてみましょう。
計算の残骸
参考になりそうな文献等
相変わらず、この文献がMCMC関連はよくまとまっております。
Monte Carlo Statistical Methods (Springer Texts in Statistics)
- 作者: Christian Robert,George Casella
- 出版社/メーカー: Springer
- 発売日: 2005/08/24
- メディア: ハードカバー
- この商品を含むブログ (2件) を見る