「一朋友問我說：」
? 飛哥,，你知道回歸分析中利用的是最小二乘法,，比如最簡(jiǎn)單的單變量回歸分析，得到的有回歸系數(shù)和截距,，但是相關(guān)的標(biāo)準(zhǔn)誤是如何計(jì)算的,？？,？我：……竟然講不出來
? 「內(nèi)心小99」
? 作為杠精我是不服氣的,，就立了一個(gè)Flag，能用矩陣形式寫出步驟,，那么許多細(xì)節(jié)應(yīng)該更加清楚了,，剛好最近在學(xué)習(xí)GWAS相關(guān)理論，就繼續(xù)灌水,。每一步的理解,，都是進(jìn)步,，在我最終回頭總結(jié)時(shí)，希望我比現(xiàn)在有進(jìn)步……
? 這是一個(gè)描述女性身高和體重的數(shù)據(jù),，我們以height為X變量（自變量）,，以weight為Y變量（因變量），進(jìn)行模型的計(jì)算,。
? 計(jì)算方法參考：https://stats.idre./r/library/r-library-matrices-and-matrix-computations-in-r/
? 數(shù)據(jù)包括兩列,，第一列是身高，第二列為體重,。
「Excuse me? 這是怎么和GWAS聯(lián)系到一起的,？」
? 強(qiáng)型解釋：GWAS中有一種是LM模型（一般線性模型），可以將SNP的分型重新編碼為012（0為主等位基因純合,，1位雜合,，2為次等位基因純合），變?yōu)閿?shù)字類型,，這樣可以作為X變量,，Y變量為性狀觀測(cè)值，這就變?yōu)榱薒M模型跑GWAS,！
? > data(women)回歸系數(shù)估計(jì)：
擬合值：
殘差估計(jì):
殘差的平方：
> X <- as.matrix(cbind(1 , women$height)) 1 ]2 ]1 ] [,2 ]1 ,] 1 58 2 ,] 1 59 3 ,] 1 60 4 ,] 1 61 5 ,] 1 62 6 ,] 1 63 2.2 構(gòu)建y矩陣 > # 構(gòu)建Y矩陣 1 )1 ]1 ,] 115 2 ,] 117 3 ,] 120 4 ,] 123 5 ,] 126 6 ,] 129 第一種方法,，是直接根據(jù)公式計(jì)算：
> beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y1 ]1 ,] -87.51667 2 ,] 3.45000 第二種方法,，是用crossprod函數(shù),，在計(jì)算大數(shù)據(jù)時(shí)有優(yōu)勢(shì)
> beta.hat1 <- solve(crossprod(X), crossprod(X,y)) # solve(A,B) == solve(A)%*%B 1 ]1 ,] -87.51667 2 ,] 3.45000 兩者結(jié)果完全一樣,。
> y.hat <- X %*% beta.hat1 :5 , 1 ],3 ) # 擬合值 1 ] 112.583 116.033 119.483 122.933 126.383 1 :5 , 1 ] #原始值 1 ] 115 117 120 123 126 對(duì)擬合值和原始值作散點(diǎn)圖：
plot(y.hat,y)> residual <- y - y.hat1 ]1 ,] 2.41666667 2 ,] 0.96666667 3 ,] 0.51666667 4 ,] 0.06666667 5 ,] -0.38333333 6 ,] -0.83333333 > sigma2 <- sum((y - y.hat)^2 )/(n - p)1 ] 2.325641 > v <- solve(t(X) %*% X) * sigma21 ] [,2 ]1 ,] 35.247305 -0.539880952 2 ,] -0.539881 0.008305861 參數(shù)的標(biāo)準(zhǔn)誤是這樣計(jì)算的！?。?/p>> #standard errors of the parameter estimates 1 ] 5.9369440 0.0911365
> # 計(jì)算T值 1 ]1 ,] -14.74103 2 ,] 37.85531 > 2 * (1 - pt(abs(t.values), n - p))1 ]1 ,] 1.711082e-09 2 ,] 1.088019e-14 上面是矩陣操作計(jì)算的結(jié)果，下面我們用R語言的lm函數(shù),，對(duì)結(jié)果進(jìn)行簡(jiǎn)單線性回歸,，得出計(jì)算結(jié)果，和矩陣的結(jié)果進(jìn)行比較,。
3. 用lm函數(shù)和矩陣得到的結(jié)果對(duì)比 簡(jiǎn)單粗暴,，兩行代碼：
> # R的lm函數(shù) 87.51667 5.9369440 -14.74103 1.711082e-09 3.45000 0.0911365 37.85531 1.090973e-14 可以看到，和矩陣計(jì)算的截距和回歸系數(shù),，一樣?。隙ㄒ粯樱蝗荒憔透拇a去了,，還能在這里灌水,？,？？如果數(shù)據(jù)不達(dá)到理想,，那就嚴(yán)刑拷打,，總能得到想要的結(jié)論------數(shù)據(jù)分析行業(yè)潛規(guī)則！??！哈哈）
> beta.hat1 ]1 ,] -87.51667 2 ,] 3.45000 > head(fitted(mod))1 2 3 4 5 6 112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 > head(y.hat)1 ]1 ,] 112.5833 2 ,] 116.0333 3 ,] 119.4833 4 ,] 122.9333 5 ,] 126.3833 6 ,] 129.8333 結(jié)果也是一樣的
> head(residuals(mod))1 2 3 4 5 6 2.41666667 0.96666667 0.51666667 0.06666667 -0.38333333 -0.83333333 1 ]1 ,] 2.41666667 2 ,] 0.96666667 3 ,] 0.51666667 4 ,] 0.06666667 5 ,] -0.38333333 6 ,] -0.83333333 summary(mod)Call:sigma22.32564102564103
sqrt(sigma2)1.52500525429948
所以，如果hight是SNP分型,，那么這個(gè)結(jié)果看那個(gè)指標(biāo)呢,？
下一篇，我們模擬一個(gè)數(shù)據(jù),，比較plink的LM模型和R的LM模型的結(jié)果……結(jié)果當(dāng)然是完全一樣的,。
當(dāng)你發(fā)現(xiàn)，GWAS分析的知識(shí)中有線性回歸分析,，T檢驗(yàn)等統(tǒng)計(jì)知識(shí),，是不是感覺踏實(shí)一點(diǎn)，起碼不僅僅是各種術(shù)語和軟件操作了,，有一點(diǎn)理解的基礎(chǔ),，是進(jìn)步的不二法門。
「其它」
? 記得我剛參加工作時(shí),，要舉辦一個(gè)統(tǒng)計(jì)軟件的培訓(xùn)（GenStat軟件）,，我準(zhǔn)備了很多內(nèi)容，把我所知道的統(tǒng)統(tǒng)都搬上來,，老板看過之后告訴我,，東西太多，太深,，培訓(xùn)把簡(jiǎn)單的內(nèi)容講透就行了,，畢竟兩天的培訓(xùn)，即使再填鴨也沒有多少效果,，反而讓聽課者畏懼,，從開始到放棄。你要把簡(jiǎn)單的內(nèi)容講透,，需要自己理解,，讓他們也建立自信，高興而來,，滿意而歸,，這才是重點(diǎn)。對(duì)一件事物不畏懼，埋頭下去研究,，慢慢就上路了,。
? ? 后來的工作中，我很受啟發(fā),，對(duì)一件新事物,，首先要消除心理的畏懼，然后像寫論文綜述一樣,，深入研究,，從多個(gè)角度查閱，慢慢就會(huì)上路,。后面的工作生涯或者學(xué)習(xí)生涯中,，無論是對(duì)于GS，還是混合線性模型,，還是Python,，Julia，還是DMU,，BLUPF90,，都有這種規(guī)律。
? ? 這里,，很適合引用村上春樹《挪威的森林》中渡邊對(duì)直子說的一句話：“我不是最聰明的,，但是我不放棄，一直琢磨,，肯定是理解你最深的人”（大意如此）,。對(duì)待一門新知識(shí)，也應(yīng)該是這種態(tài)度吧,，夜半感想,，與諸位共勉。
?
