使用矩阵操作回归分析兼论学习方法

邓飞

发布于 2020-05-18 17:08:38

7820

发布于 2020-05-18 17:08:38

文章被收录于专栏：育种数据分析之放飞自我

「一朋友问我说：」

❝飞哥，你知道回归分析中利用的是最小二乘法，比如最简单的单变量回归分析，得到的有回归系数和截距，但是相关的标准误是如何计算的？？？我：……竟然讲不出来 ❞

「内心小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联系到一起的？」

❝强型解释：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 理论模型

y = X\beta + \epsilon,\ E(\epsilon) = 0 ,\ Cov(\epsilon) = \sigma^2I_n

回归系数估计：

\widehat{\beta} = (X'X)^{-1}X'y

拟合值：

\widehat{y} = X\beta

残差估计:

\widehat{\epsilon}= y - \widehat{y}

残差的平方：

\sigma^2 = \sum{\widehat{\epsilon}^2}

2. 矩阵如何操作

2.1 构建X矩阵

> 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矩阵

> # 构建Y矩阵
> y <- matrix(women$weight,,1)
> head(y)
     [,1]
[1,]  115
[2,]  117
[3,]  120
[4,]  123
[5,]  126
[6,]  129

2.3 计算

\beta

参数

第一种方法，是直接根据公式计算：

> beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y
> beta.hat
          [,1]
[1,] -87.51667
[2,]   3.45000

第二种方法，是用crossprod函数，在计算大数据时有优势

> 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

> 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

对拟合值和原始值作散点图：

plot(y.hat,y)

在这里插入图片描述

2.5 计算残差值（residual）

> 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 计算残差方差组分（残差的方差）

> sigma2 <- sum((y - y.hat)^2)/(n - p)
> sigma2
[1] 2.325641

2.7 计算方差协方差矩阵（var.beta）

> v <- solve(t(X) %*% X) * sigma2
> v
          [,1]         [,2]
[1,] 35.247305 -0.539880952
[2,] -0.539881  0.008305861

2.8 计算参数的标准误

参数的标准误是这样计算的！！！

> #standard errors of the parameter estimates
> sqrt(diag(v))
[1] 5.9369440 0.0911365

2.9 对参数进行T检验，计算T值

> # 计算T值
> t.values <- beta.hat/sqrt(diag(v))
> t.values
          [,1]
[1,] -14.74103
[2,]  37.85531

2.10 根据T值，计算p值

> 2 * (1 - pt(abs(t.values), n - p))
             [,1]
[1,] 1.711082e-09
[2,] 1.088019e-14

上面是矩阵操作计算的结果，下面我们用R语言的lm函数，对结果进行简单线性回归，得出计算结果，和矩阵的结果进行比较。

3. 用lm函数和矩阵得到的结果对比

3.1 对比参数估计

简单粗暴，两行代码：

> # 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

可以看到，和矩阵计算的截距和回归系数，一样！（肯定一样，要不然你就改代码去了，还能在这里灌水？？？如果数据不达到理想，那就严刑拷打，总能得到想要的结论------数据分析行业潜规则！！！哈哈）

> beta.hat
          [,1]
[1,] -87.51667
[2,]   3.45000

3.2 拟合值

> 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 残差值

> 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

sqrt(sigma2)

1.52500525429948

所以，如果hight是SNP分型，那么这个结果看那个指标呢？

回归系数
Pvalue

下一篇，我们模拟一个数据，比较plink的LM模型和R的LM模型的结果……结果当然是完全一样的。

当你发现，GWAS分析的知识中有线性回归分析，T检验等统计知识，是不是感觉踏实一点，起码不仅仅是各种术语和软件操作了，有一点理解的基础，是进步的不二法门。

「其它」

❝记得我刚参加工作时，要举办一个统计软件的培训（GenStat软件），我准备了很多内容，把我所知道的统统都搬上来，老板看过之后告诉我，东西太多，太深，培训把简单的内容讲透就行了，毕竟两天的培训，即使再填鸭也没有多少效果，反而让听课者畏惧，从开始到放弃。你要把简单的内容讲透，需要自己理解，让他们也建立自信，高兴而来，满意而归，这才是重点。对一件事物不畏惧，埋头下去研究，慢慢就上路了。 ❞
❝后来的工作中，我很受启发，对一件新事物，首先要消除心理的畏惧，然后像写论文综述一样，深入研究，从多个角度查阅，慢慢就会上路。后面的工作生涯或者学习生涯中，无论是对于GS，还是混合线性模型，还是Python，Julia，还是DMU，BLUPF90，都有这种规律。 ❞
❝这里，很适合引用村上春树《挪威的森林》中渡边对直子说的一句话：“我不是最聪明的，但是我不放弃，一直琢磨，肯定是理解你最深的人”（大意如此）。对待一门新知识，也应该是这种态度吧，夜半感想，与诸位共勉。 ❞

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-05-15，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

线性回归

r 语言

本文分享自育种数据分析之放飞自我微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

线性回归

r 语言

登录后参与评论

0 条评论

热度