「一朋友问我说：」
❝飞哥，你知道回归分析中利用的是最小二乘法，比如最简单的单变量回归分析，得到的有回归系数和截距，但是相关的标准误是如何计算的？？？我：……竟然讲不出来 ❞
「内心小99」
❝作为杠精我是不服气的，就立了一个Flag，能用矩阵形式写出步骤，那么许多细节应该更加清楚了，刚好最近在学习GWAS相关理论，就继续灌水。每一步的理解，都是进步，在我最终回头总结时，希望我比现在有进步…… ❞
1.1 数据来源：来源R语言默认的数据集women
这是一个描述女性身高和体重的数据，我们以height为X变量（自变量），以weight为Y变量（因变量），进行模型的计算。
❝计算方法参考：https://stats.idre.ucla.edu/r/library/r-library-matrices-and-matrix-computations-in-r/ ❞
1.2 查看数据
数据包括两列，第一列是身高，第二列为体重。
「Excuse me? 这是怎么和GWAS联系到一起的？」
代码语言：javascript复制❝强型解释：GWAS中有一种是LM模型（一般线性模型），可以将SNP的分型重新编码为012（0为主等位基因纯合，1位杂合，2为次等位基因纯合），变为数字类型，这样可以作为X变量，Y变量为性状观测值，这就变为了LM模型跑GWAS！ ❞
> data(women)
> head(women)
height weight
1 58 115
2 59 117
3 60 120
4 61 123
5 62 126
6 63 129
1.3 理论模型
回归系数估计：
拟合值：
残差估计:
残差的平方：
2. 矩阵如何操作
2.1 构建X矩阵
代码语言：javascript复制> X <- as.matrix(cbind(1, women$height))
> n <- dim(X)[1]
> p <- dim(X)[2]
> head(X)
[,1] [,2]
[1,] 1 58
[2,] 1 59
[3,] 1 60
[4,] 1 61
[5,] 1 62
[6,] 1 63
2.2 构建y矩阵
代码语言：javascript复制> # 构建Y矩阵
> y <- matrix(women$weight,,1)
> head(y)
[,1]
[1,] 115
[2,] 117
[3,] 120
[4,] 123
[5,] 126
[6,] 129
2.3 计算
参数
第一种方法，是直接根据公式计算：
代码语言：javascript复制> beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y
> beta.hat
[,1]
[1,] -87.51667
[2,] 3.45000
第二种方法，是用crossprod函数，在计算大数据时有优势
代码语言：javascript复制> beta.hat1 <- solve(crossprod(X), crossprod(X,y)) # solve(A,B) == solve(A)%*%B
> beta.hat1
[,1]
[1,] -87.51667
[2,] 3.45000
两者结果完全一样。
2.4 计算拟合值Fitted_value
代码语言：javascript复制> y.hat <- X %*% beta.hat
> round(y.hat[1:5, 1],3) # 拟合值
[1] 112.583 116.033 119.483 122.933 126.383
> y[1:5, 1] #原始值
[1] 115 117 120 123 126
对拟合值和原始值作散点图：
代码语言：javascript复制plot(y.hat,y)
在这里插入图片描述
2.5 计算残差值（residual）
代码语言：javascript复制> residual <- y - y.hat
> head(residual)
[,1]
[1,] 2.41666667
[2,] 0.96666667
[3,] 0.51666667
[4,] 0.06666667
[5,] -0.38333333
[6,] -0.83333333
2.6 计算残差方差组分（残差的方差）
代码语言：javascript复制> sigma2 <- sum((y - y.hat)^2)/(n - p)
> sigma2
[1] 2.325641
2.7 计算方差协方差矩阵（var.beta）
代码语言：javascript复制> v <- solve(t(X) %*% X) * sigma2
> v
[,1] [,2]
[1,] 35.247305 -0.539880952
[2,] -0.539881 0.008305861
2.8 计算参数的标准误
参数的标准误是这样计算的！！！
代码语言：javascript复制> #standard errors of the parameter estimates
> sqrt(diag(v))
[1] 5.9369440 0.0911365
2.9 对参数进行T检验，计算T值
代码语言：javascript复制> # 计算T值
> t.values <- beta.hat/sqrt(diag(v))
> t.values
[,1]
[1,] -14.74103
[2,] 37.85531
2.10 根据T值，计算p值
代码语言：javascript复制> 2 * (1 - pt(abs(t.values), n - p))
[,1]
[1,] 1.711082e-09
[2,] 1.088019e-14
上面是矩阵操作计算的结果，下面我们用R语言的lm函数，对结果进行简单线性回归，得出计算结果，和矩阵的结果进行比较。
3. 用lm函数和矩阵得到的结果对比
3.1 对比参数估计
简单粗暴，两行代码：
代码语言：javascript复制> # R的lm函数
> mod <- lm(weight ~ height,data=women)
> summary(mod)$coef
Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.51667 5.9369440 -14.74103 1.711082e-09
height 3.45000 0.0911365 37.85531 1.090973e-14
可以看到，和矩阵计算的截距和回归系数，一样！（肯定一样，要不然你就改代码去了，还能在这里灌水？？？如果数据不达到理想，那就严刑拷打，总能得到想要的结论------数据分析行业潜规则！！！哈哈）
代码语言：javascript复制> beta.hat
[,1]
[1,] -87.51667
[2,] 3.45000
3.2 拟合值
代码语言：javascript复制> 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
结果也是一样的
3.3 残差值
代码语言：javascript复制> head(residuals(mod))
1 2 3 4 5 6
2.41666667 0.96666667 0.51666667 0.06666667 -0.38333333 -0.83333333
> head(residual)
[,1]
[1,] 2.41666667
[2,] 0.96666667
[3,] 0.51666667
[4,] 0.06666667
[5,] -0.38333333
[6,] -0.83333333
summary(mod)
Call:
lm(formula = weight ~ height, data = women)
Residuals:
Min 1Q Median 3Q Max
-1.7333 -1.1333 -0.3833 0.7417 3.1167
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -87.51667 5.93694 -14.74 1.71e-09 ***
height 3.45000 0.09114 37.85 1.09e-14 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared: 0.991, Adjusted R-squared: 0.9903
F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14
sigma2
2.32564102564103
代码语言：javascript复制sqrt(sigma2)
1.52500525429948
所以，如果hight是SNP分型，那么这个结果看那个指标呢？
- 回归系数
- Pvalue
下一篇，我们模拟一个数据，比较plink的LM模型和R的LM模型的结果……结果当然是完全一样的。
当你发现，GWAS分析的知识中有线性回归分析，T检验等统计知识，是不是感觉踏实一点，起码不仅仅是各种术语和软件操作了，有一点理解的基础，是进步的不二法门。
「其它」
❝记得我刚参加工作时，要举办一个统计软件的培训（GenStat软件），我准备了很多内容，把我所知道的统统都搬上来，老板看过之后告诉我，东西太多，太深，培训把简单的内容讲透就行了，毕竟两天的培训，即使再填鸭也没有多少效果，反而让听课者畏惧，从开始到放弃。你要把简单的内容讲透，需要自己理解，让他们也建立自信，高兴而来，满意而归，这才是重点。对一件事物不畏惧，埋头下去研究，慢慢就上路了。 ❞
❝后来的工作中，我很受启发，对一件新事物，首先要消除心理的畏惧，然后像写论文综述一样，深入研究，从多个角度查阅，慢慢就会上路。后面的工作生涯或者学习生涯中，无论是对于GS，还是混合线性模型，还是Python，Julia，还是DMU，BLUPF90，都有这种规律。 ❞
❝这里，很适合引用村上春树《挪威的森林》中渡边对直子说的一句话：“我不是最聪明的，但是我不放弃，一直琢磨，肯定是理解你最深的人”（大意如此）。对待一门新知识，也应该是这种态度吧，夜半感想，与诸位共勉。 ❞
