做limma时要用到model.matrix函数,这是r base stats包里的建立设计矩阵的函数,找了下资料,记录备忘。
StatQuest学习笔记05——线性模型 很好的介绍线性模型的资料
Formaldehyde
carb optden
1 0.1 0.086
2 0.3 0.269
3 0.5 0.446
4 0.6 0.538
5 0.7 0.626
6 0.9 0.782
summary(fm1 <- lm(optden ~ carb, Formaldehyde))
Call:
lm(formula = optden ~ carb, data = Formaldehyde)
Residuals:
1 2 3 4 5 6
-0.006714 0.001029 0.002771 0.007143 0.007514 -0.011743
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.005086 0.007834 0.649 0.552
carb 0.876286 0.013535 64.744 3.41e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.008649 on 4 degrees of freedom
Multiple R-squared: 0.999, Adjusted R-squared: 0.9988
F-statistic: 4192 on 1 and 4 DF, p-value: 3.409e-07
model.matrix(fm1)
(Intercept) carb
1 1 0.1
2 1 0.3
3 1 0.5
4 1 0.6
5 1 0.7
6 1 0.9
attr(,"assign")
[1] 0 1
> predict(fm1)
1 2 3 4 5 6
0.09271429 0.26797143 0.44322857 0.53085714 0.61848571 0.79374286
> 0.005086+0.876286*0.1
[1] 0.0927146
设计矩阵中0,1代表着关和开,可以作为分类资料的factor因素,比如t检验,anova分析检验如果矩阵中为非0,1,代表着线性回归即直线模型,代表着计量的资料
上例formaldehyde数据的线性回归模型,通过model.matrix得出的结果就是非factor,代表每个点对应的横坐标数值。
将此横坐标代表回归方程计算出来的Y值即是通过predict函数推断出Y值。
另个例子
> a = gl(3,4)
> b = gl(4,1,12)
> a
[1] 1 1 1 1 2 2 2 2 3 3 3 3
> b
[1] 1 2 3 4 1 2 3 4 1 2 3 4
Levels: 1 2 3 4
model.matrix(~ a + b)
(Intercept) a2 a3 b2 b3 b4
1 1 0 0 0 0 0
2 1 0 0 1 0 0
3 1 0 0 0 1 0
4 1 0 0 0 0 1
5 1 1 0 0 0 0
6 1 1 0 1 0 0
7 1 1 0 0 1 0
8 1 1 0 0 0 1
9 1 0 1 0 0 0
10 1 0 1 1 0 0
11 1 0 1 0 1 0
12 1 0 1 0 0 1