育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

用矩阵的方法计算回归分析参数

已有 5123 次阅读 2017-12-23 18:05 |个人分类:农学统计|系统分类:科研笔记


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 查看数据

data(women)

head(women)

图片.png

1.3 理论模型

图片.png

2.1 构建X矩阵

X <- as.matrix(cbind(1, women$height))  

n <- dim(X)[1]

p <- dim(X)[2]

head(X)

图片.png

2.2 构建y矩阵

y <- matrix(women$weight,,1)

head(y)

图片.png

2.3 计算 \beta 参数

第一种方法,是直接根据公式计算:

beta.hat <- solve(t(X) %*% X) %*% t(X) %*% y

beta.hat

图片.png

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

beta.hat1 <- solve(crossprod(X), crossprod(X,y)) # solve(A,B) == solve(A)%*%B

beta.hat1

图片.png

2.4 计算拟合值Fitted_value

y.hat <- X %*% beta.hat

round(y.hat[1:5, 1],3) # 拟合值

y[1:5, 1] #原始值

图片.png

plot(y.hat,y)

图片.png

2.5 计算残差值(residual)

residual <- y - y.hat

head(residual)

图片.png

2.6 计算残差方差组分(残差的方差)

sigma2 <- sum((y - y.hat)^2)/(n - p)

sigma2

图片.png

2.7 计算方差协方差矩阵(var.beta)

v <- solve(t(X) %*% X) * sigma2
v

图片.png

2.8 计算参数的标准误

#standard errors of the parameter estimates
sqrt(diag(v))

图片.png

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

t.values <- beta.hat/sqrt(diag(v))

t.values

图片.png

2.10 根据T值,计算p值

2 * (1 - pt(abs(t.values), n - p))

图片.png

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

3.1 对比参数估计

mod <- lm(weight ~ height,data=women)
summary(mod)$coef

图片.png

beta.hat

图片.png

3.2 拟合值

head(fitted(mod))
head(y.hat)

图片.png

3.3 残差值

head(residuals(mod))
head(residual)

图片.png

summary(mod)

图片.png

sigma2

2.32564102564103

sqrt(sigma2)

1.52500525429948

生物统计与数量遗传学公众号



https://blog.sciencenet.cn/blog-2577109-1091120.html

上一篇:如何用logistic回归曲线拟合日照和生物量积累的模型
下一篇:markdown 插入数学公式语法汇总
收藏 IP: 111.202.84.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 19:51

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部