answer ex3 - nibb-unix/gitc202402-unix GitHub Wiki

復習問題３　統計学入門

1. 以下のシミュレートしたデータを使って、線形モデルのあてはめを行ってみる。

問1-1. 生成されたデータdataMatrixの内容を確認し、barplot関数でstrain、repおよび残差で観察値observed.valueを構成しているかを確認しなさい（ヒント：行列の転置が必要です）。

下記のresidualは乱数発生し得ているので毎回異なります。よって、下記residual、observed.valueは一例です。

> dataMatrix
      genotype replicate    residual observed.value
 [1,]        3         2  0.08617822       5.086178
 [2,]        3         5 -0.36373252       7.636267
 [3,]        3         1 -0.64732017       3.352680
 [4,]        5         2 -0.13692730       6.863073
 [5,]        5         5 -0.71168865       9.288311
 [6,]        5         1  0.68113961       6.681140
 [7,]        7         2 -0.51255856       8.487441
 [8,]        7         5  0.23269009      12.232690
 [9,]        7         1 -0.13235088       7.867649

> barplot(t(dataMatrix[,1:3]), legend.text=c("residual", "replicate", "strain"))

問1-2.dataFrame、lm関数を使ってstrainを説明変数、valueを応答変数とする線形モデルのあてはめを行い、推定された係数を確認しなさい。

> fm1 <- lm(value ~ strain, data=dataFrame)
> summary(fm1)

Call:
lm(formula = value ~ strain, data = dataFrame)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0057 -1.0418 -0.7478  1.6775  2.7034 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)      5.358      1.171   4.576  0.00379
strainstrain1    2.252      1.656   1.360  0.22263
strainstrain2    4.171      1.656   2.519  0.04536

Residual standard error: 2.028 on 6 degrees of freedom
Multiple R-squared:  0.5145,	Adjusted R-squared:  0.3526 
F-statistic: 3.179 on 2 and 6 DF,  p-value: 0.1144

問1-3. 同様にstrainとrepを説明変数とする線形モデルのあてはめを行い、推定された係数を確認しなさい。

> fm2 <- lm(value ~ strain+rep, data=dataFrame)
> summary(fm2)

Call:
lm(formula = value ~ strain + rep, data = dataFrame)

Residuals:
       1        2        3        4        5        6        7        8        9 
 0.41506  0.05829 -0.47336 -0.06051 -0.54213  0.60263 -0.35456  0.48383 -0.12928 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)     4.6711     0.4467  10.456 0.000473
strainstrain1   2.2525     0.4894   4.603 0.010011
strainstrain2   4.1709     0.4894   8.523 0.001040
rep2            2.9069     0.4894   5.940 0.004029
rep3           -0.8451     0.4894  -1.727 0.159277

Residual standard error: 0.5994 on 4 degrees of freedom
Multiple R-squared:  0.9717,	Adjusted R-squared:  0.9435 
F-statistic: 34.37 on 4 and 4 DF,  p-value: 0.002353

2. t検定における帰無仮説とp値をRを使って確認してみよう。

問2-1. nullPvaluesのヒストグラムを描き、この帰無仮説のもとのp値の分布を求めなさい。

 hist(nullPvalues)

ヒストグラムのビンの幅を指定して再度ヒストグラムを描いてみる(breaks引数）。

 hist(nullPvalues, breaks=seq(0,1, by=0.05))

問2-2. nullPvaluesのうち、0.05以下の検定はいくつあったか？　length, which等を使って確認しなさい。

length(which(nullPvalues <= 0.05))

問2-3. このシミュレーションを再度行った際には同様な結果が得られるか？確認しなさい。 sapply(1:1000, function(x) generate.null.t.test(200,3))から再度行うことで確認できる。

nullPvalues <- sapply(1:1000, function(x) generate.null.t.test(200,3))
hist(nullPvalues, breaks=seq(0,1, by=0.05))
length(which(nullPvalues <= 0.05))

（もう一歩踏み込んだ場合）上記を繰り返してp値が0.05以下となる頻度を調べてみる。

NlowerThan0.05 <- sapply(1:100, function(x)
  length(which(sapply(1:1000, function(x) generate.null.t.test(200,3)) <= 0.05))
)
hist(NlowerThan0.05)

3. 以下はシロイヌナズナのトランスクリプトームデータからコサイン係数（類似度）を求め、あるコサイン係数が統計的にも有意かを調べる例です。

(コサイン係数については実践編の「多変量解析」の資料を参照）相関係数が0ではないことを検定する方法は別にありますが、この例では「このデータが取りうるランダムな相関係数とは有意に異なるか」の２つの帰無仮説について調べます。

3-1: データ中のある遺伝子の発現プロファイルと他のある遺伝子の発現プロファイルのコサイン係数が、データ内の数値から偶然得られる値とは有意に異なるかを_p_=0.05_を基準に評価する。

library(coop)　　# コサイン係数を高速に計算する関数が含まれているcoopパッケージを呼び出す

まずは関数の定義を見ていきましょう。わかりにくそうな箇所はコメントをつけました。

## function

# random sampling of values across samples:
# check whether or not observed correlation value can be observed in a random data

compute.null.cosine.gene <- function(row){
  x <- dat                             # 元の発現量比をオブジェクトxに代入（コピー）する
  #rownames(x) <- colnames(x) <- NULL  # この関数で作られる無作為化した発現量比データ(オフジェクトpermutated)は
                                       # サンプルをランダムに並び替えるので、サンプル名が残っていても意味がないので、お好みでNULLに置き換え
 
  permutated <- x                      # オブジェクトのコピー：　上記の行名・列名をNULLにしない場合はコードとして冗長　
  colnames(x) <- sample(colnames(x))　　　　 # サンプル名を無作為化

  permutated <- rbind(                 # サンプルの名前を無作為化した発現量比データの1行目を
    x[row,], permutated                # rbind関数でオブジェクトpermutatedに行として追加する
  )
  tcosine(permutated)[row,-row]        # 無作為化データからコサイン係数を計算し、計算結果を取り出す
}

# random sampling of correlations to collect correlation of random gene pairs:
# get distribution of correlation for random gene pairs to check a correlation 
# of interest is particularly high or low in the data
compute.null.cosine.genes <- function(alpha){
  rhoMatrix <- tcosine(dat)
  rhoDist   <- as.dist(rhoMatrix)

  nullCosines <- sample(rhoDist, nrow(dat), replace=TRUE)
  quantile(nullCosines, c(alpha/2, 1-alpha/2))
}

以下、計算を順に追っていきます。まず、発現量比データの読み込みです。

## main
dat <- read.delim("~/gitc/data/3_statistics/Sato_A_thaliana-P_syringae_arvRpt2_6h_expRatio_full.txt", header=TRUE, row.names=1)

次に帰無仮説となるコサイン係数の分布と実測データのコサイン係数を求めます。

nullCosineMatrixForGene  <- sapply(1:1000, function(x) compute.null.cosine.gene(1))

obs <- tcosine(dat)[1,-1]

次は、帰無仮説の相関係数（無作為化したデータから求めた相関係数）の大きい値から2.5パーセンタイルの値（ランダムに起きたと考えにくい値: p < 0.05の両側検定の場合は片側は0.025の値）を求め、閾値とします。

RhoThreshold  <- round(quantile(nullCosineMatrixForGene, 0.975), 3)

以上をまとめて可視化します。

par(mfrow=c(2,1))
par(mar=c(4,4,1,1))
histData   <- hist(
  obs, main="", 
  xlab=expression(paste("Observed correlations (", italic(R), ")"))
)
abline(v=RhoThreshold, col="magenta")
labelX      <- RhoThreshold
labelString <- as.character(RhoThreshold)
text(x=RhoThreshold+0.05, y=max(histData$count), 
     label=expression(paste(italic("R")[p=0.05])),
     adj=c(1,1))
par(mar=c(5,4,1,1))
hist(nullCosineMatrixForGene,
  main="", 
  xlab=expression(paste("Null distribution (", italic(R), ")"))
)
abline(v=RhoThreshold, col="magenta")

3-2.データ中のある遺伝子の発現プロファイルと他のある遺伝子の発現プロファイルのコサイン係数が、データ内の全ての遺伝子ペアから偶然得られる値とは有意に異なるかを_p_=0.05を基準に評価するための閾値となるコサイン係数を求める。

こちらも考え方は基本的には3-1と同じです。検定のモチベーションとして違うのは「データ内の全ての遺伝子ペアから偶然得られる値」を帰無仮説としているところです。

## function
# random sampling of correlations to collect correlation of random gene pairs:
# get distribution of correlation for random gene pairs to check a correlation 
# of interest is particularly high or low in the data
compute.null.cosine.genes <- function(alpha){
  rhoMatrix <- tcosine(dat)
  rhoDist   <- as.dist(rhoMatrix)

  nullCosines <- sample(rhoDist, nrow(dat), replace=TRUE)
  quantile(nullCosines, c(alpha/2, 1-alpha/2))
}

## main
## 上記3-1から継続
nullCosineMatrixForGenes <- sapply(1:1000, function(x) compute.null.cosine.genes(0.05))
RhoThresholdAllGenes     <- rowMeans(nullCosineMatrixForGenes)

cat("Null correlation values at p-value in two-sided:", 
    RhoThresholdAllGenes, sep="\n")

ここで得られる両側検定の値は毎回異なるが目安として-0.8587046, 0.8947758程度の値が得られる。

注目して欲しいのは異なる帰無仮説を仮定することによって、帰無仮説のコサイン係数の計算の仕方が異なり、それにより閾値として得られる値の大きさが異なる点です。相関係数を求めた際に得られるp値は「相関係数がゼロと異なるか」を検定していますが、それでは甘い（遺伝子発現プロファイルのような何らかの数字を入力として相関係数を計算してもゼロではない数値が得られる）と考える場合、3-1のようにランダムなデータから得られた値を帰無仮説とします。 3-2の場合は、ある遺伝子ペアの相関が「他の遺伝子ペアすべてに比べて」有意に大きいかを調べたい場合を想定しています。