一般化線形混合モデル(GLMM) -Poisson/NormalをMCMCで書き直す

GLMM(Poisson - Normal)をMCMCに！

つい、10時間ほど前に投稿したGLM（ポアソン回帰モデル）のmcmc版の作成に続けて、GLMMのポアソン/正規モデルもmcmcで書いてみました。（１０時間前の記事がこちら）

<a href="http://tomoshige-n.hatenablog.com/entry/2015/05/05/183458">ポアソン回帰モデルをMCMCに書き直してみる - Data Science by R and Python</a>tomoshige-n.hatenablog.com

こちらは、とにかく変量効果の提案分布を構成するのにとにかく手こずりました。個人的に、truncated normalが一番良さそうだと思っていますが、適当な（よく使われる）提案分布って何かあれば教えてください。（log-normal / gamma）とか的外れなのに、いろいろ試して時間を無駄にして、あげく、やはりthresholdがある提案分布が楽だなーという安易な発想のもと、、、ここに落ち着いております。

使っているアルゴリズムは、MHアルゴリズムです．相変わらず、紙に書き直して、計算をできる限り速くできるように、logを取ったり、うまいこと変数を消したりして、ある程度の効率化は図ってみました。下の方に、残骸としてはっつけております（笑。

追記（5月15日）

以下のコードには、間違いがありますので、コードについてはこちらの記事をご覧ください．<a href="http://tomoshige-n.hatenablog.com/entry/2015/05/15/195527">パッケージを使わないで、一般化線形混合モデルのMCMCアルゴリズムを1から作る． - Data Science by R and Python</a>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")