Statistics for Genomics Naoki Kimoto (@_kimoton) 【第３回】ゼロから始めるゲノム解析

統計の背後の考え ● 実験の技術的偏り ○ 技術者の腕 ○ シーケンサー ● 実験環境 ○ 温度 ○ 湿度 ● サンプルの自然変動 ○ サンプル劣化 FIGURE 3.1: PAX6遺伝子から 20回遺伝子発現量を取得した結果 Why？ 生物学や他の多くの分野では、データは実験によって収集される。 様々な要因により、毎回同じ測定値を取得することはできない。 観測値がばらつく原因 2

平均と中央値 平均 median()関数 numeric型ベクトルの中央値を計算 中央値 （n: 奇数） （n: 偶数） mean()関数 numeric型ベクトルの平均を計算 > runif(10) > mean(x) ## [1] 0.3738963 > runif(10) > median(x) ## [1] 0.3277896 生物学や他の多くの分野では、データは実験によって収集される。 中央値は外れ値の影響を受けづらいという特徴がある。 runif(n) は一様分布からn 個乱数生成する関数 runif(n) は一様分布からn 個乱数生成する関数 3

分散と標準偏差 分散 sd()関数 numeric型ベクトルの標準偏差を計算 標準偏差 var()関数 numeric型ベクトルの分散を計算 > x=rnorm(20,mean=6,sd=0.7) > var(x) ## [1] 0.2531495 > x=rnorm(20,mean=6,sd=0.7) > sd(x) ## [1] 0.5031397 データを要約する他の指標として、値の変動性を数値化することが考えられる。 「分散」と「標準偏差」はこの指標として用いられる。 rnorm(n) は正規分布からn 個乱数生成する関数 rnorm(n) は正規分布からn 個乱数生成する関数 4

四分位範囲と分位数 値の大部分がどこに分布しているのかを調べることも役に立つ。 四分位範囲は分散と異なり、外れ値の影響を受けにくいという特徴を持つ。 分位数 quantile()関数 numeric型ベクトルの 分位数を計算 四分位範囲 IQR()関数 numeric型ベクトルの 四分位範囲を計算 > x=rnorm(20,mean=6,sd=0.7) > quantile(x) ## 0% 25% ## 5.437119 5.742895 ## 50% 75% 100% ## 5.860302 6.243991 6.558112 > x=rnorm(20,mean=6,sd=0.7) > IQR(x) ## [1] 0.5010954 > boxplot(x,horizontal = T) 四分位範囲 5

正規分布 分布は平均や分散といったパラメータを持ちます。 最も頻繁に現れる分布として、正規分布が挙げられます。正規分布はベルのような形状をしており、 平均に近い値が発生しやすい分布となっています。 正規分布のパラメータごとの分布 確率変数Xが 平均 、分散 の正規分布に従うことを示す。 6

Rで正規分布 * norm関数のファミリー（rnorm、dnorm、qnorm、およびpnorm）を用いることで正規分布から得 られる様々な値を取得することができます。 # 平均0、標準偏差2の正規分布から、 X=-2 における確率密度関数を取得 > dnorm(-2, mean=0, sd=2) ## [1] 0.1209854 # 平均0、標準偏差2の正規分布から、 P(X =< -2)を取得 > pnorm(-2, mean=0, sd=2) ## [1] 0.1586553 # 平均0、標準偏差2の正規分布から 5つの乱数を取得 > rnorm(5, mean=0, sd=2) ## [1] -1.8109030 -1.9220710 -0.5146717 0.8216728 -0.7900804 # 平均0、標準偏差2の正規分布から P(X =< t)=0.15 を満たすtを取得 > qnorm(0.15, mean=0, sd=2) ## [1] -2.072867 7

その他Rで使える代表的な分布関数 # 成功確率が 0.3 である試行を 10 回行ったときに、 4 回成功する確率を計算。 > dbinom(4, n=10, p=0.3) ## [1] 0.2001209 # 期待値λ=1のポアソン分布でを取得 > dpois(1, lambda=1) ## [1] 0.3678794 # 平均0、標準偏差2の負の二項分布から 5つの乱数を取得 > dnbinom(1, size=1, prob=0.1) ## [1] 0.09 # 自由度(10,20)のF分布から、X=1 における確率密度関数の値を取得 > df(1, df1=10, df2=20) ## [1] 0.7143568 # 自由度1のχ二乗分布から、 X=1 における確率密度関数の値を取得 > dchisq(1, df=1) ## [1] 0.2419707 負の二項分布は、シーケンス読み取りカウ ントなどのモデル化に登場！ 二項分布は、二値（0, 1）がある確率で生じ る場合に登場！（ex: メチル化） F分布やχ二乗分布は線形モデルや一般 化線形モデルを使用する際に登場！ 8

標本 母集団の推定 標本統計量が母数をどれだけ精度よく表せているか、 通常母集団が得られていることはないため、標本から母集団を推定することになる。 母集団 推定！ 9

標本 ブートストラップリサンプリング 1つの標本から復元抽出を繰り返して大量の標本を生成し、それらの標本から推定値 を計算し、母 集団の性質やモデルの推測の誤差などを分析する方法をブートストラップサンプリングという。 母集団 推定！ n回観測 統計量を記録 10

信頼区間とは サンプリングを行った際に定義される区間。 「繰り返し無作為抽出を実行し、そのたびに信頼区間を構成すると、そのうちの約 95%の信頼区間が 母平均を含む」 ※母平均は定数であって確率変数ではない。確率的に分布するのは標本である。 母平均 11

# 必要なパッケージの読み込み library(mosaic) # シード値の設定 (毎回同じ結果にするための設定 ) set.seed(21) sample1= rnorm(50,20,5) # サンプルを取得 # 重複を許容し、ブートストラップサンプリングを実行 boot.means= do(1000) * mean(resample(sample1)) # 上記で得られたサンプルからパーセンタイルを取得 q=quantile(boot.means[,1],p=c(0.025,0.975)) # ヒストグラムをプロット hist(boot.means[,1], col="cornflowerblue", border="white", xlab="sample means") # パーセンタイル値に赤線を引く abline(v=c(q[1], q[2] ),col="red") text(x=q[1],y=200,round(q[1],3),adj=c(1,0)) text(x=q[2],y=200,round(q[2],3),adj=c(0,0)) 信頼区間の導出（リサンプリング使用） 1000個のブートストラップサンプルを使用した サンプル平均の精度推定 12

中心極限定理（CLT） 平均 、分散　 の母集団 から無作為抽出した 個の標本の平均　 は 平均 、分散　　　　の正規分布に近づく。 　　　は標準正規分布に従うため、 Zスコアが算出可能 Zスコアに変換した値が特定の範囲内に存在する確率か ら、信頼区間を導出することが可能！ 標本平均は、母集団分布に関係なく正規分布に従う 平均が0、標準偏差 (SD) が1に なるように変換した得点 13

信頼区間の導出（数値的導出） 中心極限定理から、標本平均は標準正規分布に従うため、 Zスコアの信頼区間から標本平均の信頼区間を数理的に導出できる。 95％信頼区間の導出 α ％信頼区間 14

Rで標準正規分布から信頼区間の導出 # 中心極限定理から数理的に 95%信頼区間を推定 sample1= rnorm(50,20,5) # サンプルを取得 alpha=0.05 sd=5 n=50 # 標準正規分布から信頼区間を推定 mean(sample1) + qnorm(c(alpha/2, 1-alpha/2)) * sd/sqrt(n) ## [1] 19.23638 22.00819 ブートストラップで得られる信頼区 間 [19.21, 21.989] 中心極限定理から数理的に 得られる信頼区間 [19.23638, 22.00819] qnorm() はZスコアを算出する関数 15

母集団の標準偏差分からない問題 一般に母集団の統計量（母数）は分からない … ↓ でも標本の標準偏差（　）を使用して標準偏差を推定したい！ ↓ 信頼区間の計算に標準正規分布の代わりに t分布（→） を使用すれば、　　　　　　　　から信頼区間が導出可能！ 標本から母集団の信頼区間を求める際は t分布を使用する　 標本 母集団 推定！ 16

Rでt分布から信頼区間の導出 ブートストラップで得られる信頼区 間 [19.21, 21.989] # 中心極限定理から数理的に 95%信頼区間を推定 sample1= rnorm(50,20,5) # サンプルを取得 alpha=0.05 sd=5 n=50 # t分布から信頼区間を推定 mean(sample1) + qt(c(alpha/2, 1-alpha/2), df=n-1) * sd/sqrt(n) ## [1] 19.20130 22.04327 qt() はt分布の分位点を算出する関数 t分布から数理的に得られる信頼区 間 [19.20130, 22.04327] 17

3.2 How to test for differences between samples サンプルセット間の違いを検証する方法 18

サンプルセット間の比較を行いたい事例 ①野生型サンプルセットの発現が変異体のサンプルセットと異なる場合 ②健康なサンプルセットと疾患があるサンプルセットの測定可能な特徴（血球数、遺伝 子発現、特定の遺伝子座のメチル化）が異なる場合 A群とB群に差はある？ C群とD群に差はある？ 知りたいこと 19

どのようにサンプルセット間比較を行うか 最も簡単な方法：2つのサンプルの平均値の差によって判断する。 ⇒ 平均値に差があればサンプルセット同士には差があると判断し、 　 差がなければサンプルセット同士には差がないと判断する。 平均値に差があれば サンプルセット同士に差がある？？ A群の 平均値 B群の 平均値 20

どのようにサンプルセット間比較を行うか 最も簡単な方法：2つのサンプルの平均値の差によって判断する。 ⇒ 平均値に差があればサンプルセット同士には差があると判断し、 　 差がなければサンプルセット同士には差がないと判断する。 知っておかなければならない点：測定値にはばらつき(誤差)がある。 (測定器に依存した誤差、測定者に依存した誤差、環境条件に依存した誤差) そのため、サンプル比較する際にはその誤差を考慮する必要がある。 ⇒平均値の差には誤差が含まれることを理解する。 　(ある程度の差が生まれてしまうことは仕方がない。) 21

誤差ではないと言える平均値の差とは 2つのサンプルの平均値の差が誤差によるものかを判断するためには、 誤差と言える閾値を設定し、平均値の差がそれ以上の値であるかを判断すればよい。 ⇒「仮説検定」と呼ばれる手法の概念 知りたいこと この差は誤差によるものなのか、そう ではないのか？ A群の 平均値 B群の 平均値 22

仮説検定の手順 ①帰無仮説と対立仮説を設定する。 「帰無仮説(H 0 )」：検証したい仮説 ⇒今回の場合、帰無仮説は、「サンプル間に差がない」こと。 「対立仮説(H 1 )」：帰無仮説に対する仮説 (本来証明したい仮説) ⇒今回の場合、対立仮説は、「サンプル間に差がある」こと。 ②帰無仮説の真偽を検証するための統計量 (今回の場合、差があると言える閾値 )を決定する。 ③統計値を計算し、それを参照値 (今回の場合、平均値の差)と比較して、意味のある差(有意差)がある のかをP値を用いて判断する。 ④結果、帰無仮説H 0 を棄却するか棄却しないかを判断する。 23

3.2.1 Randomization-based testing for difference of the means ランダム化を用いた平均値の差の検証方法 24

ランダム化を用いた平均値の差の検証方法 概念 サンプルセット間に差がないと思われる場合は、 サンプルセットのラベル(例: “健康” と “病気”)には意味がないはず。 2つのサンプルセットをごちゃ混ぜにして、 再度ランダムに2つのサンプルセットに分けた時の平均値の差と 実際のサンプルセットの平均値の差を比較して、 サンプルセットのラベルには意味がなかったのかを検証する。 ⇒ 実際にサンプルセット間には違いがなかったのか（帰無仮説）を検証する。 25

ランダム化を用いた平均値の差の検証方法 手順 1) サンプルセットにラベルをランダムに割り当てて平均値 の差を計算することを繰り返す（可能であれば全ての 組み合わせを試す）。 ⇒平均値の差の分布（帰無分布）が計算できる。 2) 帰無仮説が成り立つ場合に、実際の平均値の差の値 を得る確率（頻度）を計算する。 3) 得た確率が設定した有意水準以下であれば、有意な 差があると判断する。（帰無仮説を棄却する。） ランダムセットの 平均値の差の分布(帰無分布) 有意 水準 実際の 平均値の差 26

# シード値を設定 (毎回同じ結果にするための設定 ) set.seed(100) # 2つの異なる分布から 2つのサンプルセットを取得 gene1=rnorm(30,mean=4,sd=2) gene2=rnorm(30,mean=2,sd=2) # 平均値の差を計算 org.diff=mean(gene1)- mean(gene2) gene.df=data.frame(exp=c(gene1,gene2), group=c(rep("test",30),rep("control",30))) # 帰無分布を取得 exp.null <- do(1000) * diff(mosaic::mean(exp ~ shuffle(group), data=gene.df)) # ヒストグラムを描画 hist(exp.null[, 1],xlab="null distribution | no difference in samples" , main=expression(paste(H[0]," :no difference in means" )), xlim=c(-2,2),col="cornflowerblue" ,border="white") abline(v=quantile(exp.null[, 1],0.95),col="red" ) abline(v=org.diff,col="blue" ) text(x=quantile(exp.null[, 1],0.95), y=200, "0.05",adj=c(1,0), col="red") text(x=org.diff,y=200,"org. diff." ,adj=c(1,0),col="blue") Rでの実装 27

Rでの実装 # P値（実際の平均値の差より高くなる確率）の算出 p.val=sum(exp.null[,1]>org.diff)/length(exp.null[,1]) p.val ## [1] 0.001 28

3.2.2 Using t-test for difference of the means between two samples t検定による2つのサンプルセット間の平均値の比較 29

t検定 ランダム化による平均値比較の欠点 ①サンプルが少なすぎて十分なランダム化テストを実行できない場合がある。 ②ランダム化に時間がかかる。 ⇒これらの欠点を解決するのがt検定 t検定とは：t分布に依存する検定 中心極限定理から平均値の差がt分布に従うことを示すことができる。 t検定にはいくつかの種類がある。 30

スチューデントのt検定 (母分散が等しいと仮定する場合) 母分散が等しいと仮定すると、次の式が成り立つ。 上式は、平均値の差が自由度n 1 +n 2 -2のt分布に従うことを意味する。 スチューデントのt検定はRで実装されており、自由度n 1 +n 2 -2のt分布における、計算し た統計量tの値のP値が取得できる。 , 31

ウェルチのt検定 (母分散が等しいと仮定できない場合) 母分散が等しいと仮定できない場合、次の統計量を用いる。 d.f.は “degrees of freedom”の略で、自由度を意味する。 ウェルチのt検定はRで実装されており、容易に利用可能である。 , , 32

Rで実装 (ウェルチのt検定) # Welch’s t-test stats::t.test(gene1, gene2) # ## ## Welch Two Sample t-test ## ## data: gene1 and gene2 ## t = 3.7653, df = 47.552, p-value = 0.0004575 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.872397 2.872761 ## sample estimates: ## mean of x mean of y ## 4.057728 2.185149 33

Rで実装 (母分散が等しいと仮定したt検定) # t-test with equal variance assumption stats::t.test(gene1, gene2, var.equal=TRUE ) # 出力 ## ## Two Sample t-test ## ## data: gene1 and gene2 ## t = 3.7653, df = 58, p-value = 0.0003905 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.8770753 2.8680832 ## sample estimates: ## mean of x mean of y ## 4.057728 2.185149 34

t検定 (補足) t検定は、サンプルの母集団が正規分布に従うことを想定しているが、 2つのサンプルの分布が同じ方向に適度に歪んでいる場合、 t検定は正規分布からの逸脱を許容出来ることが示されている。 これは、サンプルサイズが大きい場合、母集団の分布に関係なく、 サンプルの平均値が正規分布に従うという中心極限定理によるものである。 標本サイズが多ければ多いほど、標準検定値である は Z値に近似することになる。 35

3.2.3 Multiple testing correction 多重検定の修正 36

検定の評価方法 TN/(FP+TN): 特異度 TP/(TP+FN): 感度(検定の検出力) 帰無仮説を 採択した 帰無仮説を 棄却した 帰無仮説が 真であった 帰無仮説が 偽であった (帰無仮説が真である数 ) (対立仮説が真である数 ) (第2種の過誤) (第1種の過誤) 高くすることが目標 37

多重検定の修正 検定の数が増えると… 第1種の過誤(False Positives)が増えると予想される。 例えば、5％の有意水準で検定を実行すると、 帰無仮説が真である場合、帰無仮説を誤って棄却する可能性が5%ある。 帰無仮説が全て当てはまる1000回の検定を行うと、 誤った棄却の平均の数は50になる。 また、少なくとも誤った棄却が1回はほぼ100%の確率で発生する。 38 10,000遺伝子についてマイクロ アレイの発現量を比較する場合 など

多重検定の修正法① 多重検定における第1種の過誤を防ぐための統計手法は複数ある。 ①ボンフェローニ補正 全検定中でひとつでも偽陽性が得られる確率がある目標値となるように、 　複数の検定から得られたP値をより高い値に補正する。 　具体的には、個々のP値が十分に低い場合に、 　個々のP値(p i )に検定数(m)を掛ける(m×p i )。 問題点： 陽性判定の条件が厳しくなりすぎるため、帰無仮説を棄却しづらい 39

多重検定の修正法② ②Benjamini-Hochberg法 ある程度第1種の過誤を許容して、第2種の過誤を起こす可能性を小さくする方 法の一つとして、False Discovery Rateを調整するという方法が考案された。 P値をランク付けし、検定数(m) × 個々のP値(p i ) をランクiで割る (m×p i / i) 　上記の手順により、「FDR調整済みP値（=q-value）」を取得し、 　特定のテストでの誤検出の割合を見積もることができる。 FDR調整済みP値が0.05であることは、 棄却された帰無仮説の5％が誤りであることを意味する。 問題点： p-valueの分布を一様分布と仮定している。 False Discovery Rate (FDR) : 真の帰無仮説のうち誤って棄却した割合 FP / (TP + FP) 40

# ライブラリ読み込み library(qvalue) data(hedenfalk) qvalues <- qvalue(hedenfalk$p)$q bonf.pval=p.adjust(hedenfalk$p,method ="bonferroni") fdr.adj.pval=p.adjust(hedenfalk$p,method ="fdr") plot(hedenfalk$p,qvalues,pch=19,ylim=c(0,1), xlab="raw P-values" ,ylab="adjusted P-values" ) points(hedenfalk$p,bonf.pval,pch=19,col="red") points(hedenfalk$p,fdr.adj.pval, pch=19,col="blue") legend("bottomright",legend=c("q-value","FDR (BH)","Bonferroni"), fill=c("black","blue","red")) Rでの実装 出力 41

3.3 Relationship between variables: Linear models and correlation 変数間の関係：線形モデルと相関 42

変数間の関係を測定/モデル化したい in Genomics 薬を X 投与！ 肝臓である遺伝子 発現量が Y 上昇！ X 箇所がメチル化！ メチル化された遺伝 子発現量が Y 減少！ X と Y の関係をモデル化したい！ 43

線形回帰（1変数） 目的（応答）変数 説明（従属）変数 誤差項 H3K4me3ヒストン修飾と遺伝子発現の関係 H3K4me3ヒストン修飾の増加は遺伝子発現の増加に 関連しているように見える。 線形回帰では、説明変数、目的変数及び 誤差項を用いてモデリングされる。 44

線形回帰（2変数） 目的（応答）変数 説明（従属）変数 誤差項 遺伝子発現とH3K4me3、H3K27me3ヒストン修飾との関係 H3K4me3ヒストン修飾の増加は遺伝子発現の増加に、 H3K27me3ヒストン修飾の増加は遺伝子発現の減少に、 関連しているように見える。 線形回帰では、説明変数、目的変数及び 誤差項を用いてモデリングされる。 45

線形回帰（多変数） より多くの変数を説明変数に持つモデルの場合、行列式を用いることで簡潔に記述でき る。 46

3.3.1 How to fit a line 線形フィッティングの方法 47

βの値によるコスト関数の大きさ コスト関数/損失関数を用いたアプローチ 「コスト関数」もしくは「損失関数」と呼ばれる関数を最適化する。 一般的にこの関数には、予測値 と観測値 の二乗誤差が使用される。 最適化手順（最急降下法） 1. ランダムな開始点（ , ）を選択する。 2. コスト関数の偏導関数から、コスト関数をどの方向に最適化するかを決定。 3. コスト関数を最小化する方向に向かって最適化を進める。ステップの大きさはパラ メータとして与える。 4. 収束するまで手順2、3を繰り返す。 48

平均 、分散 に従う正規分布を仮定した場合、応答変数を観測する 確率は下記となる。 最尤関数を用いたアプローチ① 最尤推定法では、与えられたデータセット内の全ての応答変数を観測する確率を最大 化するような を求める。 尤度関数 は、全てのデータポイントについて観測される確率の乗算となる。 対数変換 MAX y=ln(x) は連続 49

最尤関数を用いたアプローチ② についての導関数を取得することにより、 を最小化する は、 下記を満たす となる。 について、 コスト関数を用いたアプローチと全く同じ結果が得られる。 最尤推定を用いたアプローチでは、前提の分布に正規分布以外の分布を仮定すること で、様々なフィッティングが可能となる。 50

Rで線形回帰 # 乱数シードを設定 set.seed(32) # 1から100までの50個のX値を取得 x = runif(50,1,100) # b0とb1、分散（sigma）を設定 b0 = 10 b1 = 2 sigma = 20 # 正規分布から誤差項をシミュレート eps = rnorm(50,0,sigma) ＃線形方程式と誤差項の追加からy値を取得 y = b0 + b1*x+ eps mod1=lm(y~x) # データをプロット plot(x,y,pch=20, ylab="Gene Expression", xlab="Histone modification score") # 回帰直線をプロット abline(mod1,col="blue") 線形回帰によってモデル化された遺伝 子発現とヒストン修飾スコア lmはlinear modelの略 51

3.3.2 How to estimate the error of the coefficients 回帰係数の誤差を推定する方法 52

回帰係数の推定標準誤差 標本を使用して係数を推定しているため、 ランダムサンプリングごとに異なる推定値が得られる。 誤差項の分散が一定で平均がゼロであると仮定すると、標本 から母集団の特性を推測するための特性として、 回帰係数の不確実性を下記の通りモデル化することができる （証明は略）。 回帰係数は、ランダムサンプルごとに異なる 回帰係数の変動をヒストグラムで表示 は RSE、及びXの 分散に依存 53

Rで線形回帰（結果の確認） # 線形回帰結果を取得 > mod1=lm(y~x) > summary(mod1) ## Call: ## lm(formula = y ~ x) ## Residuals: ## Min 1Q Median 3Q Max ## -38.50 -12.62 1.66 14.26 37.20 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.40951 6.43208 0.064 0.95 ## x 2.08742 0.09775 21.355 <2e-16 *** ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## Residual standard error: 16.96 on 48 degrees of freedom ## Multiple R-squared: 0.9048, Adjusted R-squared: 0.9028 ## F-statistic: 456 on 1 and 48 DF, p-value: < 2.2e-16 # βの信頼区間を取得 > confint(mod1) ## 2.5 % 97.5 % ## (Intercept) -12.523065 13.342076 ## x 1.890882 2.283952 # βの推定値を抽出 > coef(mod1) ## (Intercept) x ## 0.4095054 2.0874167 推定値、標準誤差、t値、p値 を表示 54

3.3.3 Accuracy of the model モデルの精度 55

RSE（Residual standard error） 残差標準誤差。誤差項の二乗和を自由度で割った平方根。 pは自由度であり、1変数の線形回帰の場合p=2（直線の傾きと y 切片）となる。モデル の適合性を評価する指標となる。大きいほど適合度が悪い。 56

決定係数（ ） 決定係数は０から１までの値をとり、総変動のうち推定されたモデルでどの程度説明さ れているか。残差変動が全変動に対してどれだけ少ないかを表す。 総変動 残差変動 相関係数との関係 相関係数の二乗！ 57

F統計量（分散分析） TSS - RSSは「説明された変動性」、RSSは「説明されていない変動性」 帰無仮説が成り立つとき、この比率は説明された分散の自由度と説明されていない分 散の自由度の2つのパラメーターを使用したF分布に従う。 全ての係数を考慮した検定を行える。 → 多重検定問題を避けられる！ 帰無仮説は全ての係数が0 58

3.3.4 Regression with categorical variables カテゴリ変数による回帰 59

カテゴリ変数を用いた回帰 # 線形回帰結果を取得 set.seed(100) gene1=rnorm(30,mean=4,sd=2) gene2=rnorm(30,mean=2,sd=2) gene.df=data.frame(exp=c(gene1,gene2), group=c(rep(1,30),rep(0,30))) mod2=lm(exp~group,data=gene.df) summary(mod2) ## ## Call: ## lm(formula = exp ~ group, data = gene.df) ## ## Residuals: ## Min 1Q Median 3Q Max ## -4.7290 -1.0664 0.0122 1.3840 4.5629 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.1851 0.3517 6.214 6.04e-08 *** ## group 1.8726 0.4973 3.765 0.000391 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.926 on 58 degrees of freedom ## Multiple R-squared: 0.1964, Adjusted R-squared: 0.1826 ## F-statistic: 14.18 on 1 and 58 DF, p-value: 0.0003905 require(mosaic) plotModel(mod2) 60

3.3.5 Regression pitfalls 回帰の落とし穴 61

線形回帰における落とし穴 非線形性 真の関係性が線形からかけ離れている場合、予測精度が低下する 場合によって , , 等の変換を検討する必要がある。 説明変数の相関 （多重共線性） 説明変数間の相関が高い場合、係数を適切に推定できない。 係数の推定標準偏差が大きくなる傾向が見られる。 誤差項の相関 誤差項が互いに無相関であることを前提としている。 相関がある場合、係数の信頼区間が狭くなる傾向がある。 誤差項の非一定分散 予測変数の値に関係なく、異なる応答変数の誤差に同じ分散がある。 エラーが一定でない場合（例： X値が増加するとエラーが増加します）、 標準誤差の推定値が信頼できなくなる。 外れ値と高レバレッジ 点 外れ値：Y に関する極値、高レバレッジ点：異常な X値を指す。 これらは、近似直線と標準誤差に影響を与える可能性がある。 場合によってはデータから削除して適合性を高めるべき。 62