2016-09-11

勝手にStan

Stan rstan

Stanを使えば事後分布が乱択的に得られる
けれど、単純なモデルではつまらないし、普通に回帰等をする方がよっぽどよい
どれくらいモデルが複雑だとStanがよいかを言い切るほどよくわかっていないけれど、階層ベイズはそうなのだろう(こちら)
SNPの遺伝因子評価に使うとどうなるか、いじってみる
ひとつは、性・年齢が単純でなく関わっているとき。ただしStanでのモデル化は単純なので、結局、glm()と大差ないはず
もうひとつは、SNP-phenotype関連と、phenotype 発現量関連と、SNP-eqtl関連を統合するようなもので、phenotypeが入っている関連では個人が共通で、eqtlは全く別の個人セットとする(とおもったけれど、もう少し複雑にしないとやりたいことをやっていないようだった…要修正)
ひとつめ

functions{}

data{
	int<lower=0> nsample;
	int<lower=0,upper=1> phenotype[nsample];
	int<lower=0,upper=2> homo[nsample];
	int<lower=0,upper=1> hetero[nsample];
	int<lower=0,upper=150> age[nsample];
	int<lower=0,upper=1> sex[nsample];
}

transformed data{
}

parameters{
	real rhomo;
	real<lower=0,upper=1> rhetero;
	real rsex;
	real rage;
}

transformed parameters{
	vector[nsample] risk;
	for(i in 1:nsample){
		risk[i] <- rhomo * homo[i] + rhetero * rhomo * hetero[i] + rsex * sex[i] + rage * age[i];
	}
}

model{
	rhomo ~ normal(0,1000);
	rhetero ~ uniform(0,1);
	rsex ~ normal(0,1000);
	rage ~ normal(0,1000);
	
	for(i in 1:nsample){
		phenotype[i] ~ bernoulli(inv_logit(risk[i]));
	}
}

generated quantities{}

gp <- c(0.49,0.42,0.09)
n <- 10^3
g <- sample(0:2,n,replace=TRUE,prob=gp)
homo <- hetero <- rep(0,n)
homo[which(g==2)] <- 1
hetero[which(g==1)] <- 1
sex <- sample(0:1,n,replace=TRUE,prob=c(0.3,0.7))
library(MCMCpack)
age <- sample(20:70,n,replace=TRUE,prob=rdirichlet(1,rep(1,length(20:70))))
age. <- age-40
age.[which(age.<0)] <- 0

ph <- (g+1)^2 + sex * 2 + (age./max(age.))^1.8 * 5 + ((age-30)/max(age-30))^2*(-3.4)
ph. <- ph + rnorm(n,0,mean(ph))

ph.obs <- (ph.-min(ph.))/(max(ph.-min(ph.)))
ph.obs. <- ph.obs > mean(ph.obs)
glm(ph.obs. ~ g,family="binomial")

PH <- as.numeric(ph.obs.)

dat <- list(phenotype=PH,nsample=n,age=age,sex=sex,homo=homo,hetero=hetero)
stan.out <- stan(file="SNPcovariate.stan",data=dat)

stan_hist(stan.out,pars=c("rhomo","rhetero","rage","rsex"))
glm.out <- glm(PH~g+hetero+age+sex,family="binomial")

ふたつめ

functions{}

data{
	int<lower=0> nsample;
	int<lower=0,upper=1> phenotype[nsample];
	int<lower=0,upper=2> homo[nsample];
	int<lower=0,upper=1> hetero[nsample];
	int<lower=0,upper=150> age[nsample];
	int<lower=0,upper=1> sex[nsample];
}

transformed data{
}

parameters{
	real rhomo;
	real<lower=0,upper=1> rhetero;
	real rsex;
	real rage;
}

transformed parameters{
	vector[nsample] risk;
	for(i in 1:nsample){
		risk[i] <- rhomo * homo[i] + rhetero * rhomo * hetero[i] + rsex * sex[i] + rage * age[i];
	}
}

model{
	rhomo ~ normal(0,1000);
	rhetero ~ uniform(0,1);
	rsex ~ normal(0,1000);
	rage ~ normal(0,1000);
	
	for(i in 1:nsample){
		phenotype[i] ~ bernoulli(inv_logit(risk[i]));
	}
}

generated quantities{}

ngenotype <- 1000
g <- sample(0:2,ngenotype,replace=TRUE,prob=c(0.49,42,0.09))
pre.ph <- g * 3 + 1 + rnorm(ngenotype,1)
ph <- as.numeric(pre.ph > mean(pre.ph))

nexpression <- 500
gexpression <- sample(0:2,nexpression,replace=TRUE)
expression <- rep(0,nexpression)
for(i in 1:nexpression){
	expression[i] <- rnorm(1,gexpression[i],1)
}

blood <- rep(0,ngenotype)
for(i in 1:ngenotype){
	blood[i] <- rnorm(1,g[i],1)
}

dat <- list(ngenotype=ngenotype,neqtl=nexpression,genotype=g,phenotype=ph,expression=expression,gexpression=gexpression,blool=blood)

stan.out <- stan(file="GEBP.stan",data=dat)