Linear models

最后发布时间 : 2022-08-18 20:49:29 浏览量 :

本文是将线性模型用于显著性测试的直观理解，主要包括以下两个方面：

显示线性模型如何等效于t-test
线性模型的灵活性

从t-test开始

Question

如下所示，从RWSBIM2122软件包加载GEXP数据。它提供了两组（这里是1和2）中感兴趣的一个基因的log2表达量，暂时忽略ID列。

install.packages("remotes")
remotes::install_git("https://gitee.com/bioinfoFungi/rWSBIM2122.git")

library(rWSBIM2122)
data(gexp)
gexp

##    expression group ID
## 1         2.7     1  1
## 2         0.4     1  2
## 3         1.8     1  3
## 4         0.8     1  4
## 5         1.9     1  5
## 6         5.4     1  6
## 7         5.7     1  7
## 8         2.8     1  8
## 9         2.0     1  9
## 10        4.0     1 10
## 11        3.9     2  1
## 12        2.8     2  2
## 13        3.1     2  3
## 14        2.1     2  4
## 15        1.9     2  5
## 16        6.4     2  6
## 17        7.5     2  7
## 18        3.6     2  8
## 19        6.6     2  9
## 20        5.4     2 10

使用t检验，比较第1组和第2组的基因表达。特别是，报告各组基因的平均表达、 log2 fold-change、相关p-value，并可视化数据。

线性模型的等价性

现在，我们使用线性回归和lm()函数来对第1组和2组中的基因表达进行建模：
$y=\beta_0+\beta_1+x+\epsilon$

y是我们要建模的基因表达量
x代表两组，这里编码为0和1

本例中的基因表达可以这样描述：

当x=0时，第一组中的基因表达为 $\beta_0$ （第一组的平均基因表达）
当x=1时，第二组的基因表达为 $\beta_0+\beta_1$ （\beta_1表示第一组和第二组之间的差异）

使用R语言建模如下：

mod1 <- lm(expression ~ group, data = gexp)
mod1

## 
## Call:
## lm(formula = expression ~ group, data = gexp)
## 
## Coefficients:
## (Intercept)       group2  
##        2.75         1.58

这里可以看到：

截距对应于第1组中的平均值
模型的斜率（group2 的系数）对应于log2 fold-change ，正如期待的那样，第二组的平均值等于斜率+截距

我们可以通过以下方式将该线性模型可视化：

ggplot(data=gexp,aes(group,expression))+
	geom_boxplot()+
	geom_jitter()+
	geom_point(data = data.frame(x = c(1, 2), y = c(2.75, 4.33)),
               aes(x = x, y = y),
               colour = "red", size = 5) +
    geom_abline(intercept = coefficients(mod1)[1] - coefficients(mod1)[2],
                slope = coefficients(mod1)[2])

图片alt

请注意，上面，我们需要在绘制线性模型时从截距中减去斜率，因为组1和2在模型中编码为0和1，但在图上绘制为1和2。

我们还可以提取与系数相关的p值：

summary(mod1)

## 
## Call:
## lm(formula = expression ~ group, data = gexp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.430 -1.305 -0.580  1.455  3.170 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.7500     0.6004   4.580 0.000232 ***
## group2        1.5800     0.8491   1.861 0.079187 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.899 on 18 degrees of freedom
## Multiple R-squared:  0.1613, Adjusted R-squared:  0.1147 
## F-statistic: 3.463 on 1 and 18 DF,  p-value: 0.07919

that are computer as

$t_{score} = \frac{\beta_i - 0}{SE_{\beta_i}}$
如果零假设为真，则具有n-2个自由度的t-distribution。我们使用t-test来测试斜率是否不同于0：

$H_0:\beta_i = 0$
$H_0:\beta_i \neq 0$

斜率 $\beta_1$ 的p值略有不同，因为t-test中未设置的均匀性假设(homoscedasticity assumption)，我们可以设置var.equal = TRUE以确保与上面结果相同

t.test(expression ~ group, data = gexp, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  expression by group
## t = -1.8608, df = 18, p-value = 0.07919
## alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
## 95 percent confidence interval:
##  -3.363874  0.203874
## sample estimates:
## mean in group 1 mean in group 2 
##            2.75            4.33

配对t-test

dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design= ~ batch + condition)

通过在设计中添加变量，可以控制计数的额外变化。例如，如果条件样本在实验批次之间是平衡的，则通过batch在设计中包含该因子，可以提高发现差异的敏感性condition。

DESeq2

# create random data of 12 samples and 4800 genes
x <- round(matrix(rexp(480 * 10, rate=.1), ncol=12), 0)
rownames(x) <- paste("gene", 1:nrow(x))
colnames(x) <- paste("sample", 1:ncol(x))

coldata <- data.frame(
  condition = factor(c(
    rep("ctl", 3),
    rep("A", 3),
    rep("B", 3),
    rep("C", 3))))
	
dds <- DESeqDataSetFromMatrix(
  countData = x,
  colData = coldata,
  design= ~ condition)
  
  dds <- DESeq(dds, betaPrior = FALSE)

Linear models

从t-test开始

Question

线性模型的等价性

配对t-test

参考