- single-cell RNA-seq解析の手法論文を読むとき、何がハードルを上げているか、というと、「個別には理解できなくはないような、色々」を使って「総合芸術化」していることでしょうか
- いろいろな、データ処理のアイテムについて「この単語はだいたいこういうこと」という当たりがつけば、その細部はぼんやりしていても、解析戦略の全体像を誤解することはないし、誤解していなければ、手法アプリのマニュアル通りに実行できれば、大きく間違いはないのですが
- 各アイテムについて、一部だけでもまったくわからない、というものがあると、その手法を使うのは不安になるでしょう
- そんな不安解消のために、少しまとめてみました
- 英語/日本語併記です
-
- もちろん、Rmdファイルからhtml化してもＯＫ
- Rmdファイル
---
title: "SingleCellRNAseq"
author: "ryamada"
date: "Monday, January 25, 2016"
output: html_document
---
If you want to strengthen your knowledge on single-cell RNAseq data analysis, read these, for example, "Survey and Summary Single-cell RNA-seq: advances and future challenges", Nucleic Acids Research (2014) 42(12) 8845, http://www.ncbi.nlm.nih.gov/pubmed/25053837 , and "Computational and analytical challenges in single-cell transcriptomics", Nature Reviews Genetics (2015) 16(3):133 http://www.ncbi.nlm.nih.gov/pubmed/25628217 . They will help your understanding of this hands-on course.
このハンズオンコースを始める前に、１細胞RNAseqデータ解析に関して基礎固めをしておくとよいでしょう。たとえば、 "Survey and Summary Single-cell RNA-seq: advances and future challenges", Nucleic Acids Research (2014) 42(12) 8845, http://www.ncbi.nlm.nih.gov/pubmed/25053837 , と "Computational and analytical challenges in single-cell transcriptomics", Nature Reviews Genetics (2015) 16(3):133 http://www.ncbi.nlm.nih.gov/pubmed/25628217 は役に立ちます。
Scan the paper below on a method proposal.
This is one of method papers for single-cell RNAseq data analysis.
The purpose to scan this paper is to identify components of data analysis used in single-cell RNAseq studies.
Other methods use the same components.
Therefore studying the components will make you able to understand various other methodology papers.
以下の手法提案論文をざっと眺めてみましょう。１細胞RNAseq解析手法を提案する論文の一つです。
この論文をざっと眺めるのは、１細胞RNAseqスタディに用いられるデータ解析手法を構成要素がなんであるかを拾い出すためです。
他の手法も同様の要素から成り立っています。
そうして拾い出した要素について理解をすると、色々な手法論文が理解できるようになります。
Check the highlited terms in the paper "MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data", Genome Biology (2015) 16:278 http://www.genomebiology.com/2015/16/1/278 .
This comment is for a seminar and when you have to scan the paper yourself, no highlight is available.
次の論文のマーキングしたところをチェックしましょう。ただし、セミナーでハイライトしたペイパーが配布されることを意図したコメントなので、この文書を読んで勉強するだけの人は、自分で全体の流れをつかみましょう。
"MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data", Genome Biology (2015) 16:278 http://www.genomebiology.com/2015/16/1/278 .
* Survey and Summary Single-cell RNA-seq: advances and future challenges. Nucleic Acids Research (2014) 42(12) 8845 http://www.ncbi.nlm.nih.gov/pubmed/25053837
* Computational and analytical challenges in single-cell transcriptomics. Nature Reviews Genetics (2015) 16(3):133 http://www.ncbi.nlm.nih.gov/pubmed/25628217
* MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biology (2015) 16:278 http://www.genomebiology.com/2015/16/1/278 , https://github.com/RGLab/MAST/blob/master/README.md
* edgeR: differential expression analysis of digital gene expression data User's guide. Bioinformatics (2010) 26:1 R-package "edgeR" available @ Bioconductor https://bioconductor.org/packages/release/bioc/html/edgeR.html
The papers in the reference list proposed methods that are implemented as R packages, "MAST" and "edgeR". You can install them as below.
上記の引用文献にあるペイパーが手法を提案していて、それらは、Rのパッケージ "MAST"と"edgeR"として入手可能です。
インストール方法は以下の通り。
```{}
#install.packages(c("curl","Biobase"))
#install.packages('devtools')
#library(devtools)
#install_github('RGLab/MAST')
#library(MAST)
# *or* if you don't have a working latex setup
#install_github(RGLab/'MAST', build_vignettes=FALSE)
#vignette('MAST-intro')
```
```{}
#source("https://bioconductor.org/biocLite.R")
#biocLite("edgeR")
#library(edgeR)
#browseVignettes("edgeR")
```
* Tests 検定
+ t-test
+ ANOVA
+ Likelihood ratio test 尤度比検定
* Distributions 分布
+ Poisson distribution ポアソン分布
+ Negative binomial distribution 負の二項分布
+ Bimodal distribution 二峰性分布
+ Log transformation 対数変換
+ Logit transformation ロジット変換
* Regression and estimation 回帰と推定
+ Linear regression 線形回帰
+ Generalized linear regression 一般化線形回帰
+ Maximum likelihood estimation 最尤推定
+ Bayesian approach to estimate regression coefficients 回帰係数推定におけるベイジアンアプローチ
+ Shrinkage estimator with empirical Bayes 経験ベイズによるShrinkage 推定量
When you compare two groups for expression level of a gene, you may test whether the means of two groups are the same.
２群の遺伝子発現レベルを比較をするとき、発現量の平均が等しいかどうかを検定するかもしれません。
```{r}
n1 <- 10 # Gourp 1 sample size
n2 <- 12 # Gruop 2 sample size
m1 <- 10;sd1 <- 1 # mean and sd of group 1
m2 <- m1;sd2 <- sd1 # mean and sd of group 2 are the same with group 1
x1 <- rnorm(n1,m1,sd1) # random values in normal dist
x2 <- rnorm(n2,m2,sd2)
x1
x2
boxplot(x1,x2)
t.test(x1,x2,var.equal=TRUE)
```
The data can be in the shape of two vectors, a value vector and a label vector.
t.test with equal variance assumption is the same with anova-type evaluation of linear regression.
グループラベルのベクトルと値のベクトルとでデータを保持することもできます。
t-testと線形回帰をしてそれをANOVA的に評価するのとは同じ結果になります。
```{r}
x12 <- c(x1,x2)
gr12 <- c(rep("1",n1),rep("2",n2))
x12
gr12
t.test(formula=x12~gr12,var.equal=TRUE)
anova(lm(x12~gr12))
```
You may refresh your knowledge on equality of variance.
t-testにおいて、等分散を仮定するか、仮定しないか、という件が怪しい人は、要確認。
```{r}
t.test(x1,x2,var.equal=TRUE)
t.test(x1,x2,var.equal=FALSE)
```
Repeat the data simulation from null hypothesis and test them to obtain a set of p-values. Check p-value distribution.
帰無仮説データセット作成の繰り返し、それをt-testすることでp値を多数発生させ、その分布を確認すること。
```{r}
n1 <- 10
n2 <- 12
n3 <- 15
m1 <- 10;sd1 <- 1
m2 <- m1;sd2 <- sd1
m3 <- m1;sd3 <- sd1
x1 <- rnorm(n1,m1,sd1)
x2 <- rnorm(n2,m2,sd2)
x3 <- rnorm(n3,m3,sd3)
x123 <- c(x1,x2,x3)
gr123 <- c(rep("1",n1),rep("2",n2),rep("3",n3))
lm.out123 <- lm(x123~gr123)
anova(lm.out123)
```
Anova can be used for two groups. ANOVAは２群にも適用できて、それはt-testと同じになる。
```{r}
lm.out12 <- lm(x12~gr12)
anova(lm.out12)
t.test(x12~gr12,var.equal=TRUE)
```
Change number of groups from 3 to 4, 5, or 6 and generate data set and test them.
no. gr = 4,5,6に変更し、データ作成をして、テストせよ。
Check p-value distribution under null hypothesis.
帰無仮説下のp値分布を確認せよ。
$$
Y = b_1 age + b_0 + \epsilon
$$
For example, age affects on expression.
$\epsilon$ is the error term and here we assume it is in noraml distribution.
例えば、発現に年齢が影響するとする。$\epsilon$ は誤差項で、正規分布に従うとする。
```{r}
n <- 100
bage <- 1 # coefficient for age
b0 <- 50 # intercept
age <- sample(20:80,n,replace=TRUE)
age <- sort(age)
eps <- rnorm(n,0,5)
Y <- bage * age + b0 + eps
plot(age,Y,pch=20)
```
The contribution of age to gene expression in this linear model is evaluated with linear regression function, lm(), and its null hypothesis (no age effect ~ b1 = 0 ) is tested with anova() function.
発現への年齢の寄与を、lm()関数で線形回帰し、その当てはまりをanova()で評価する。
```{r}
lm.out2 <- lm(Y~age)
plot(Y ~ age,pch=20)
lines(age, lm.out2$fitted, type="l", col="red")
anova(lm.out2)
```
When our interest is not in age but groups,
年齢が発現に影響するのは当然として、それ以外に群が影響するかどうかに興味があるときには
$$
Y = b_{gr} X_{gr} + b_{age} age + b_0 + \epsilon .
$$
$b_{age}$ value can be anything but we care whether $b_{gr} = 0$ or not.
年齢の係数$b_{age}$はどんな値でも気にしない。その代りグループの寄与を表す係数$b_{gr}$は気にする。
```{r}
n1 <- 100
n2 <- 120
bage <- 1
b0 <- 50
age1 <- sample(20:80,n1,replace=TRUE)
age2 <- sample(20:80,n2,replace=TRUE)
eps1 <- rnorm(n1,0,5)
eps2 <- rnorm(n2,0,5)
bgr1 <- 0
bgr2 <- 5
Y1 <- bgr1 + bage * age1 + b0 + eps1
Y2 <- bgr2 + bage * age2 + b0 + eps2
plot(c(age1,age2),c(Y1,Y2),col=c(rep(1,n1),rep(2,n2)),pch=20)
```
```{r}
age12 <- c(age1,age2)
Y12 <- c(Y1,Y2)
gr12 <- c(rep("1",n1),rep("2",n2))
lm.out3 <- lm(Y12~gr12+age12)
anova(lm.out3)
ord.age <- order(age12)
age12.sorted <- age12[ord.age]
gr12.sorted <- gr12[ord.age]
Y12.sorted <- Y12[ord.age]
lm.out3.sorted <- lm(Y12.sorted~gr12.sorted+age12.sorted)
gr1 <- which(gr12.sorted=="1")
gr2 <- which(gr12.sorted=="2")
plot(age12.sorted,Y12.sorted,pch=20,col=gr12.sorted)
lines(age12.sorted[gr1],lm.out3.sorted$fit[gr1],type="l",col="black")
lines(age12.sorted[gr2],lm.out3.sorted$fit[gr2],type="l",col="red")
```
Now you handled one covariable, age, to evaluate effect of group on gene expression.
年齢を共変量として組み込んでグループの遺伝子発現に対する影響が評価できています。
Details may be found @
http://d.hatena.ne.jp/ryamada22/20150204/1423020380 .
詳しいことがhttp://d.hatena.ne.jp/ryamada22/20150204/1423020380 で確認できます。
The followings are the essense.
以下の内容はその要点のみです。
NGS returns you integers of read counts.
There are two discrete distribtutions to model the counts; poisson distribution and negative binomial distribution.
Poisson distribution is a discrete probability distribution that expresses the probability of a given integer counts in a fixed intervales.
NGSはリードカウントとして整数値を返します。
この整数値の分布のモデルとして、ポアソン分布と負の二項分布が知られています。
ポアソン分布は離散確率分布であり、ある固定区間に事象が生起する回数に関する分布です。
```{r}
pois.mean <- 3
x <- rpois(10^5,pois.mean)
hist(x)
table(x)
mean(x)
```
```{r}
pois.mean <- 15.3
x <- rpois(10^5,pois.mean)
hist(x)
table(x)
mean(x)
pois.mean <- 0.1
x <- rpois(10^5,pois.mean)
hist(x)
table(x)
mean(x)
```
When two groups are compared and when each group should have "expected read counts", observed counts should be fit to poisson distributions; under the null hypothesis, fitting the same poisson distribution and under the alternative distribution, fitting to different poisson distributions.
2群のデータを比較するとき、それぞれの群は「リードカウントの期待値」を持っています。観測リードカウントデータを、ポアソン分布にフィットされます。
帰無仮説では、同じポアソン分布を２群に想定し、そのようなポアソン分布にフィットさせます。
対立仮説の場合は、異なるポアソン分布が背景にあるとして、２つの異なるポアソン分布にフィットさせます。
Likelihood is the hypothetical probability that an hypothesis would yiedl an observation.
尤度は、観察を前にして、ある仮説の下でそれが起きたとしたときの「仮の確率」のこと
```{r}
obs <- c(1,0,1,4,2,0)
# hypothesis: poisson distribution with mean = m
m <- 0.5
# hypothetical probability to observe each counts
hp.each <- dpois(obs,m)
hp.each
# joint hypothetical probability
hp <- prod(hp.each)
hp
```
Because joint probability tends to be very small, logarithm of likelihood, log-likelihood, is often used instead.
同時確率(複数の観察がそれぞれ起きたときに、個々の確率を掛け合わせて算出する)はとても小さくなりがちなので、対数を取ることも多く、対数尤度と呼ぶ。
```{r}
hp.each.log <- log(dpois(obs,m))
hp.each.log
hp.log <- sum(hp.each.log)
hp.log
exp(hp.log)
hp
```
Assume the observation is random samples from a poisson disrtibution with unknown average.
There is a value of average of poisson distribution, which make the likelihood maximum.
The value is called maximum likelihood estimate of average of poisson distribution.
上記の観察がポアソン分布からのランダムサンプリングの結果だと仮定する。ただし、ポアソン分布の期待値はわからないとする。
期待値の値ごとに尤度は計算できるが、尤度を最大にするような期待値と言うものもある。そのような期待値を、ポアソン分布期待値の最尤推定量と言う。
```{r}
m.candidates <- seq(from=0.01,to=3,by=0.01)
hp.logs <- rep(NA,length(m.candidates))
for(i in 1:length(hp.logs)){
hp.logs[i] <- sum(log(dpois(obs,m.candidates[i])))
}
plot(m.candidates,hp.logs)
max.hp.log <- max(hp.logs)
max.like.est <- m.candidates[which(hp.logs==max.hp.log)]
abline(v=max.like.est,col="red")
abline(h=max.hp.log,col="red")
```
Poisson distribution is one-parameter distribution, meaning that the shape of distribution is always the same when its mean is given.
Real NGS read counts distributions do not alway fit well to poisson distributions.
Negative binomial (NB) distribution is a discrete distribution takeing 0 or more with larger flexibility in the shape and NGS read-counts distribution seems to fit to NB distribution better than Poisson distribution.
Besides the expected value, NB distribution requires one more parameter.
ポアソン分布はパラメタ１つで決まる分布である。期待値(平均値)を決めるとポアソン分布の形が一意に決まる。
NGSのリードカウントの分布は、ポアソン分布に似ているけれど、それに当てはまるかというとそうでもない。
ポアソン分布と同様に0以上の整数値を取る離散分布である負の二項分布は、ポアソン分布よりも色々な形を取ることができて、NGSリードカウントの実際の分布によくフィットすることが知られている。
NB分布には、期待値の他にもう一つパラメタを指定する。
Different k values changes the shape of NB distribution.
rの値を変えるとNB分布の形が変わる。
```{r}
Np <- 3.1 # Average number of occurence
N.max <- 20
k.max <- 0:N.max
r <- 2
P.pois <- dpois(k.max,Np)
P.negbinom <- dnbinom(k.max,r,mu=Np)
ylim = range(c(P.pois,P.negbinom))
cex <- 2/log10(N.max+1)
plot(k.max,P.pois,col=2,pch=1,cex=cex,ylim=ylim)
points(k.max,P.negbinom,col=1,cex=cex,type="h")
points(k.max,P.negbinom,col=1,pch=20,cex=cex)
legend(N.max*0.7,ylim[2]*0.9, c("pois","negbinom"), col = c(1,2),pch = c(20,1,20))
```
Poisson distribution is an extreme case of NB distribution.
ポアソン分布というのは、rの値を極端にした場合のNB分布。
```{r}
r <- 10^5 # parameter for negative binomial dist
P.pois <- dpois(k.max,Np)
P.negbinom <- dnbinom(k.max,r,mu=Np)
ylim = range(c(P.pois,P.negbinom))
cex <- 2/log10(N.max+1)
plot(k.max,P.pois,col=2,pch=1,cex=cex,ylim=ylim)
points(k.max,P.negbinom,col=1,cex=cex,type="h")
points(k.max,P.negbinom,col=1,pch=20,cex=cex)
legend(N.max*0.7,ylim[2]*0.9, c("pois","negbinom"), col = c(1,2),pch = c(20,1,20))
```
The other extreme NB distribution is "always 0".
NB分布のもう一つの極端な状態は、確率1で0であるもの。
```{r}
r <- 0 # parameter for negative binomial dist
P.pois <- dpois(k.max,Np)
P.negbinom <- dnbinom(k.max,r,mu=Np)
ylim = range(c(P.pois,P.negbinom))
cex <- 2/log10(N.max+1)
plot(k.max,P.pois,col=2,pch=1,cex=cex,ylim=ylim)
points(k.max,P.negbinom,col=1,cex=cex,type="h")
points(k.max,P.negbinom,col=1,pch=20,cex=cex)
legend(N.max*0.7,ylim[2]*0.9, c("pois","negbinom"), col = c(1,2),pch = c(20,1,20))
```
When NB distribution model is used, estimate two parameters for each NB distribution, and the null hypothesis with no difference in means among groups should be tested.
NB分布モデルの場合には、2個のパラメタを推定する必要がある。その上で、平均に差があることを帰無仮説として検定する。
ある、整数観察があったときに、それがポアソン分布から得られたのか、NB分布から得られたのかを判断したいとする。
ポアソン分布としたとき、その期待値の最尤推定値を求め、その場合の尤度を計算する。
NB分布としたときの、2つのパラメタの最尤推定値を求め、その場合の尤度を計算する。
その２つの仮説(ポアソン分布とみなす仮説とNB分布とみなす仮説)の尤度の比(Likelihood ratio)を取って、検定するのが、Likelihood ratio testである。
Assume a set of integers and you want to judge whether it comes from a poisson distribution or a NB distribution.
You can have maximum likelihood estimate for boty poisson and NB distribution.
Now you can compare maximum likelihood for both hypotheses, poisson and NB distributions.
Actually the ratio of two maximum likelihoods is concerned by likelihood ratio test.
ここでは、最尤推定値を、以下のように「探索的」に求める。
すでにやったとおり、ポアソン分布仮説での、最尤推定値とそのときの尤度は以下の通り。
We search maximum likelihood estimate.
For poisson hypothesis, we already know how to do.
```{r}
obs <- rnbinom(50,30,mu=20)
# hypothesis: poisson distribution with mean = m
m.candidates <- seq(from=15,to=25,by=0.1)
hp.logs <- rep(NA,length(m.candidates))
for(i in 1:length(hp.logs)){
hp.logs[i] <- sum(log(dpois(obs,m.candidates[i])))
}
plot(m.candidates,hp.logs)
max.hp.log <- max(hp.logs)
max.like.est <- m.candidates[which(hp.logs==max.hp.log)]
abline(v=max.like.est,col="red")
abline(h=max.hp.log,col="red")
# 最尤推定値
max.like.est
# 対数最大尤度
max.hp.log
```
For NB distribution, we should seach the highest peak of likelihood in 2-dimensional parameter space.
NB分布の場合は、２つのパラメタが作る２次元パラメタ空間を探索する。
```{r}
r.candidates <- seq(from=1,to=100,by=1)
mr.candidates <- expand.grid(m.candidates,r.candidates)
hp.logs.nb <- rep(NA,length(mr.candidates[,1]))
for(i in 1:length(hp.logs.nb)){
hp.logs.nb[i] <- sum(log(dnbinom(obs,mr.candidates[i,2],mu=mr.candidates[i,1])))
}
max.hp.log.nb <- max(hp.logs.nb)
max.like.est.nb <- mr.candidates[which(hp.logs.nb==max.hp.log.nb),]
max.like.est.nb
hp.logs.nb.matrix <- matrix(hp.logs.nb,length(m.candidates),length(r.candidates))
image(m.candidates, r.candidates,hp.logs.nb.matrix)
contour(m.candidates,r.candidates,hp.logs.nb.matrix,add=TRUE)
abline(v=max.like.est.nb[1])
abline(h=max.like.est.nb[2])
```
Ratio of two likelihood is difference of two log-likelihood.
2つの尤度の比は、2つの対数尤度の差
```{r}
log.like.ratio <- max.hp.log.nb - max.hp.log
log.like.ratio
```
This log(likelihood-ratio) should be doubled and the double value is considered a value in chi-square distribution and the value should be converted to p-value of likelihood ratio test.
To use chi-square distribution, degree of freedom is necessary.
It is difference of parameters of two hypotheses, which is 1 in this case.
尤度比検定では、この対数尤度の差(対数-尤度比)の2倍の値を、カイ二乗分布に照らして検定p値とする。
カイ二乗分布は自由度を決めないと使えない。
自由度は、2つの仮説の変数の数の差なので、ここでは、1。
```{r}
2 * log.like.ratio
pchisq(2*log.like.ratio,df=1,lower.tail=FALSE)
```
p-value is not so small, which means that the observation did not reject the hypothesis with fewer parameters, i.e., poisson model was not rejected.
p値はとても小さいわけではないので、観察データは「パラメタ数の少ない仮説」を棄却しなかった、と解釈する。パラメタ数が少なかったのは、poissonモデルなので、poisson モデルが棄却されなかった、ということになる。
Change observed data set as you like and perform likelihood ratio test several times.
観察データセットを適当に変更し、尤度比検定を実施しみよ。何度か行うこと。
Some fraction of single cells have no expression for a gene and the rest have observable expression in reality.
The fraction of zero expression is $0 \le z \le 1$.
The expression level of non-zero fraction is in some distribution.
１細胞発現解析では、細胞は発現あり・なしの２群に分かれている。発現ありの群では何がしかの分布を取る。
発現なしの割合を$0 /le z /le 1$とし、発現ありのグループに分布を定める。
```{r}
z <- 0.3
n <- 10000
zero.nonzero <- sample(0:1,n,prob=c(z,1-z),replace=TRUE)
Np <- 4
r <- 20
v.tmp <- rnbinom(n,r,mu=Np)
v <- zero.nonzero * v.tmp
hist(v)
```
We made a mixture of NB distribution and "zero-fraction".
However the mixture of zeros and somehow-close to normal distribution is good to be handled for anaylsis.
Log transformation with addition of 1 to observed values of single-cell RNAseq data does this job.
This transforms original "zero" into "zero" after the procedure.
上ではあるフラクションのゼロとNB分布との混合分布を考えた。
しかしながら、ゼロ観察と正規分布に近いような分布との混合というのは、データハンドリング上、好都合である。
現実の１細胞RNAseqデータは、観察値に1を加えて対数変換すると、それらしい分布に変わる。
この変換では、もともとの観察値が0であったものが、変換後も0であるところに特徴がある。
```{r}
hist(log(v+1))
```
Now, one fraction is zero and the other fraction is normally distributed.
It is easier to handle normal distribution than negative binomial distribution, even in bayesian estimation.
We may want to control zero.nonzero as well as value in case of nonzero with "explanatory variable(s)"
上の例では、ゼロ・ノンゼロの振り分けに$0 \le z \le 1$という変数を用いた。
この$z$も推定の対象とするとすると、$z$の取りうる値が限定されるのは、できなくはないが、面倒くさいこともある。$- \infty \le z' \le \infty$となってくれれば、取りうる範囲の制約がなくなってくれる。
そのような変換の一つが、ロジット変換である。
```{r}
z <- runif(10^5)
z. <- log(z/(1-z))
par(mfcol=c(1,2))
hist(z)
hist(z.)
par(mfcol=c(1,1))
```
Linear regression is a good method. It is simple and straight-forward and its calculation is easy at least for computers.
However the "straight" model is not always appropriate for variables.
For example, when the outcome variable is 0 or 1, "logistic regression" is famous and it is away from "linear" model.
When the dependent variables are non-negative integers, they may follow poisson distribution and it does not seem good to fit a line.
The method called "generalized linear regression" broadens the models and keep the simple strategy of linear regression.
It takes a linear formula of explanatory variables and bridge the formula with the dependent variable via a link function with an appropriate error distribution.
線形回帰はよい手法。単純でわかりやすい。計算自体も少なくともコンピュータにとっては非常に簡単。
しかしながら、そのまっすぐなモデルがいつも適切なモデルとは限らない。
例えば、帰結が0/1のときにロジスティック回帰を使うのは有名だけれど、0/1しかとならいのだから、線に乗せるのは得策ではない。
従属変数が非負整数であるような場合には、直線に回帰するよりポアソン分布的に回帰するのがよい。
一般化線形回帰と呼ばれれる手法は、対応可能なモデルの範囲を広げ、かつ、線形回帰が持つ単純な仕組みを維持している。
一般化線形回帰では、説明変数に線形式を作らせ、それを従属変数とつなぐにあたり、適当なlink関数と呼ぶ関数を使い、乱雑項も適切な分布関数を使っている。
Its implement in R, glm() function, is easy to use.
その実行は、Rではglm()関数で汎用的に実行できるようになっている。
Assume an event whose average number of occurence increases along the age.
年齢とともに平均正規回数が増えるような例を考える。
```{r}
n <- 100
age <- runif(n) * 100
age <- sort(age)
a <- 0.12
b <- 1
p <- rep(NA,n)
for(i in 1:n){
p[i] <- rpois(1,a*age[i]+b)
}
plot(age,p)
glm.out <- glm(p ~ age, family=poisson(link="log"))
glm.out
plot(age,p)
lines(age,glm.out$fit,type="l",col="red")
```
Because this is a regression method, multiple variates can be included in the model.
For example sex affects the number.
回帰なので、変数を増やすことは簡単。
例えば性別が正規回数に影響することにしてみる。
```{r}
n <- 100
age <- runif(n) * 100
age <- sort(age)
sex <- sample(0:1,n,replace=TRUE)
a <- 0.12
b <- 1
a.sex <- 3
p <- rep(NA,n)
for(i in 1:n){
p[i] <- rpois(1,a*age[i]+a.sex*sex[i]+b)
}
plot(age,p,col=sex+1)
glm.out <- glm(p ~ age + sex, family=poisson(link="log"))
glm.out
plot(age,p,col=sex+1)
sex0 <- which(sex==0)
sex1 <- which(sex==1)
lines(age[sex0],glm.out$fit[sex0],type="l",col="black")
lines(age[sex1],glm.out$fit[sex1],type="l",col="red")
```
A biomarker is related with development of a disease.
Again glm() function with a bit modification of option does the job.
ロジスティック回帰も一般化線形回帰の一つ。glm()関数のオプション指定で実行できる。
```{r}
n0 <- 60
n1 <- 40
x0 <- rnorm(n0,100,10)
x1 <- rnorm(n1,130,20)
status <- c(rep(0,n0),rep(1,n1))
x <- c(x0,x1)
glm.out <- glm(status~x,family=binomial(link="logit"))
ord.x <- order(x)
x.sorted <- x[ord.x]
status.sorted <- status[ord.x]
glm.out2 <- glm(status.sorted~x.sorted,family=binomial(link="logit"))
plot(x,status)
lines(x.sorted,glm.out2$fit,col="red")
```
# Bayesian approach to estimate regression coefficients : Bayesian glm. 回帰係数の推定におけるベイズアプローチ : ベイジアンglm
Logistic regression is good but sometimes its calculation does not go well as shown in the following example.
ロジスティック回帰はよい手法だけれど、ときに、うまく回らないことがある。次の例を見てみる。
```{r}
n <- 40
x1 <- rnorm(n)
x1 <- sort(x1)
y <- rep(0,n)
y[x1 > 0] <- 1
M1 <- glm(y~x1,family=binomial(link="logit"))
```
You can see an error message.
It means that 0/1 are perfectly separated with x1.
glm()関数によるロジスティック回帰を実行するとエラーメッセージが出る。
説明変数が0/1を完璧に分離するので、推定が収束しない、というメッセージである。
```{r}
plot(x1,y)
abline(v = 0, col="red")
```
```{r}
M1
```
Addition of penalty enables you to draw a logistic curve for observation completely separable.
One way of such penalties is a bayesian approach with a prior to make the curve not too steep, as below.
ロジスティック回帰らしいカーブに納めるためには、なんらかのペナルティを加えるとよい。
ペナルティの入れ方の一つが、ベイズ推定をするやり方で、推定前にあらかじめロジスティックカーブが屹立しないような前提(プライア)を与え、そのうえで、ベイズ推定をすると、最終出力も屹立しすぎないようになる。
```{r}
library(arm)
#n <- 40
#x1 <- c(rnorm(n/2,-1,1),rnorm(n/2,1,1))
#x1 <- sort(x1)
y <- rep(0,n)
y[x1 > 0] <- 1
M2 <- bayesglm(y~x1,family=binomial(link="logit"))
M2
plot(y ~ x1)
lines(x1, M1$fitted, type="l", col="red")
lines(x1, M2$fitted, type="l", col="blue")
```
When nothing is assumed, you have to estimate parameter values based on observation.
However when similar observations are repeated for many items, like non-zero expression level distribution among cells for many genes, we have some knowledge that cell-to-cell variation of expression level should be in some shape although the shapes should vary among genes.
For example cell-to-cell variation is measured with variance, then many variance values for many genes are available.
Now we should estimate each gene's variance based on the overall variance distribution and individual gene's observed distribution.
In the extreme case where all cells had zero expression and no cell had non-zero value, the estimate of variance of the gene if any cells expressed some, should be the expected value or most likely value of the "overall variance distribution" of all genes.
On the contrary, very many cells expressed some of another gene, it is reasonable to believe the observed distribution as the true distribution of cell-to-cell variability and its variance should be estimated based on the observation only and you may ignore the overall distribution among genes.
Therefore for genes with few cells oberved, the estimate should be closer to the value that the overall distribution tells and for genes with many cells observed, the estimate should not be so close to overall indication, as the other cases. Either case, the "estimated values" becomes close to the overall distribution-indicating value; this movement towards is called "shrinkage".
Because this approach uses the overall distribution as a prior, it is a Bayesian approch.
The following example is "success rate" estimation for many items.
何も仮定することがないときには、何かを推定するにあたって、観測データそのものを使うしかない。
しかしながら、同様の観測(たとえば複数の細胞をサンプルとして、それぞれの細胞における発現量)を、たくさんのアイテム(たとえばたくさんの遺伝子)について行ったとすると、「たくさんの細胞がどのようにばらつくか」について、何がしかの情報がつかめる。
細胞のバラツキ具合を分散で測ることにすれば、分散がどのような分布を取るかについての情報がつかめる。
いま、個々の遺伝子について、発現が細胞間でどの程度ばらつくかを分散としてとらえようとするなら、この分散の分布を使うことが出来そうである。
極端な例で言えば、もしある遺伝子では1個の細胞の観察もない(すべてが発現量0の群に落ちた)とすれば、「もし発現する細胞があったとしたら、どのくらいばらつくだろうか」ということを推定することになる。今ある情報はたくさんの遺伝子のそれであるから、たくさんの遺伝子の分散の分布に照らして、期待される値なり、最もありそうな値なりを、推定値とすることになるだろう。
逆に、非常にたくさんの細胞を観察して、cell-to-cell variabilityの分布が、この遺伝子についてしっかり観察されているなら、その情報を重点的に使って、分散がどれくらいかの推定をするのがよさそうである。
したがって、観察サンプル数が少ないなら、遺伝子全体の分布により近く推定し、観察サンプル数が多いなら、それほど遺伝子全体の分布は気に掛けなくてよい。いずれにしろ、遺伝子全体の教える値に向かって、個々の遺伝子についての情報のみが教える値から"shirink"した値を推定値ととることにしよう、という手法が、この手法である。
この方法では、全体の分布を事前分布として使うことから、ベイズ流と呼ばれる。
次の例は、たくさんのアイテムの「成功率」に関するものである。
```{r,warnng=FALSE,message=FALSE}
n.item <- 10000 # number of item
n.trial.per.item <- round(2^(seq(from=1,to=15,length=n.item))) # number of trials of each item
success.rate.per.item <- rbeta(n.item,70,30) # true success rate of each item
hist(success.rate.per.item)
obs <- matrix(0,n.item,2)
for(i in 1:n.item){
tmp.n.trial <- n.trial.per.item[i]
tmp.success.rate <- success.rate.per.item[i]
tmp <- rbinom(1,tmp.n.trial,c(tmp.success.rate,1-tmp.success.rate))
obs[i,] <- c(tmp,tmp.n.trial-tmp)
}
aves <- obs[,1]/apply(obs,1,sum)
items.many <- which(n.trial.per.item>100)
hist(aves[items.many])
m <- MASS::fitdistr(aves[items.many], dbeta,start = list(shape1 = 1, shape2 = 10))
alpha0 <- m$estimate[1]
beta0 <- m$estimate[2]
plot(dbeta(seq(from=0,to=1,length=100),alpha0,beta0))
res <- (obs[,1]+alpha0)/(apply(obs,1,sum) + alpha0 + beta0)
log.n <- log(n.trial.per.item)
cl <- log.n/max(log.n)
plot(aves,res,col=rgb(cl,1-cl,1),pch=20)
abline(h=alpha0/(alpha0+beta0))
abline(v=alpha0/(alpha0+beta0))
abline(0,1,col=2)
```
The empirical Bayes estimate is same with sample mean for the items with many trials, that is indicated by the reddish dots along the red line.
The bluish dots represent the items with relatively few trials whose shrink estimates are chrunk a lot towards the overall distribution indicating value.
経験ベイズ統計量は、トライアル数が多いアイテムの場合は、サンプル平均と同じである。それは赤味の強いドットが赤い線に沿って分布していることからわかる。
逆に、トライアル数が少ないアイテムは青っぽい点であり、それは、赤い線からのずれが大きく、アイテム全体に引きずられた値になっている。
The following article is helpful as well.
以下のウェブ記事も参考になる。
http://www.r-bloggers.com/understanding-empirical-bayes-estimation-using-baseball-statistics/
