R语言使用矩阵操作回归分析 |
您所在的位置：网站首页 › r语言生成矩阵的两种方法是什么 › R语言使用矩阵操作回归分析 |
R语言使用矩阵操作回归分析
|
一朋友问我说： 飞哥，你知道回归分析中利用的是最小二乘法，比如最简单的单变量回归分析，得到的有回归系数和截距，但是相关的标准误是如何计算的？？？ 我：……竟然讲不出来 内心小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 β + ϵ , E ( ϵ ) = 0 , C o v ( ϵ ) = σ 2 I n y = X\beta + \epsilon,\ E(\epsilon) = 0 ,\ Cov(\epsilon) = \sigma^2I_n y=Xβ+ϵ, E(ϵ)=0, Cov(ϵ)=σ2In 回归系数估计： β ^ = ( X ′ X ) − 1 X ′ y \widehat{\beta} = (X'X)^{-1}X'y β =(X′X)−1X′y 拟合值： y ^ = X β \widehat{y} = X\beta y =Xβ 残差估计: ϵ ^ = y − y ^ \widehat{\epsilon}= y - \widehat{y} ϵ =y−y 残差的平方： σ 2 = ∑ ϵ ^ 2 \sigma^2 = \sum{\widehat{\epsilon}^2} σ2=∑ϵ 2 2.1 构建X矩阵 > X n p 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 head(y) [,1] [1,] 115 [2,] 117 [3,] 120 [4,] 123 [5,] 126 [6,] 129 2.3 计算 β \beta β 参数第一种方法，是直接根据公式计算： > beta.hat beta.hat [,1] [1,] -87.51667 [2,] 3.45000第二种方法，是用crossprod函数，在计算大数据时有优势 > beta.hat1 beta.hat1 [,1] [1,] -87.51667 [2,] 3.45000两者结果完全一样。 2.4 计算拟合值Fitted_value > y.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)
参数的标准误是这样计算的！！！ > #standard errors of the parameter estimates > sqrt(diag(v)) [1] 5.9369440 0.0911365 2.9 对参数进行T检验，计算T值 > # 计算T值 > t.values 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 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 sigma22.32564102564103 sqrt(sigma2)1.52500525429948 所以，如果hight是SNP分型，那么这个结果看那个指标呢？ 回归系数Pvalue下一篇，我们模拟一个数据，比较plink的LM模型和R的LM模型的结果……结果当然是完全一样的。 当你发现，GWAS分析的知识中有线性回归分析，T检验等统计知识，是不是感觉踏实一点，起码不仅仅是各种术语和软件操作了，有一点理解的基础，是进步的不二法门。 其它 记得我刚参加工作时，要举办一个统计软件的培训（GenStat软件），我准备了很多内容，把我所知道的统统都搬上来，老板看过之后告诉我，东西太多，太深，培训把简单的内容讲透就行了，毕竟两天的培训，即使再填鸭也没有多少效果，反而让听课者畏惧，从开始到放弃。你要把简单的内容讲透，需要自己理解，让他们也建立自信，高兴而来，满意而归，这才是重点。对一件事物不畏惧，埋头下去研究，慢慢就上路了。 后来的工作中，我很受启发，对一件新事物，首先要消除心理的畏惧，然后像写论文综述一样，深入研究，从多个角度查阅，慢慢就会上路。后面的工作生涯或者学习生涯中，无论是对于GS，还是混合线性模型，还是Python，Julia，还是DMU，BLUPF90，都有这种规律。 这里，很适合引用村上春树《挪威的森林》中渡边对直子说的一句话：“我不是最聪明的，但是我不放弃，一直琢磨，肯定是理解你最深的人”（大意如此）。对待一门新知识，也应该是这种态度吧，夜半感想，与诸位共勉。 |
【本文地址】
今日新闻 |
推荐新闻 |