线性回归 | Typhon

一般线性回归模型（LM）

我们在小学就学过了二元一次方程，其中 y 为因变量，x 为自变量。当 x 发生变化时，y 的值也会按照一定的规律发生改变。
一般的我们认为其表达式为:

$y = β_{0} + β_{1}x$

在现实生活中我们有成对的（x，y），想构建一条直线能够让这些点均匀的分布在这条直线的周围。所以我们需要对 β 的值（参数）进行估计。统计世界中的公式会有一个小尾巴，用来代表误差，即：

$y = β_0+β_{1}x+μ$

损失函数

有了成对的（x，y），我们可以任意地在散点间画一条直线，那么哪条直线才是最准确的呢？我们先从残差说起。残差说白了就是真实值和预测值间的差值（也可以理解为差距、距离），用公式表示是:

$e = y-\hat{y}$

我们把所有的残差先平方（消除负号），然后再都加起来得到一个值，用这个值来评判这条直线是否是最好的。残差平方和越接近0，那就说明真实值和预测值越相近。
残差平方和公式:

$Q = \sum(y_{i}-\hat{y}_{i})^{2} = \sum(y_{i}-(\hat{β}_{0}+\hat{β}_{1}x_{i}))^{2}$

还记得微积分知识的话，就知道导数为0时，Q 取最小值，因此我们对两个 β 求偏导并令其为0:

$\frac{\partial Q}{\partial{β}_{0}} = 2\sum(y_{i}-\hat{β}_{0}-\hat{β}_{1}x_{i})=0$ $\frac{\partial Q}{\partial{β}_{1}} = 2\sum(y_{i}-\hat{β}_{0}-\hat{β}_{1}x_{i})x_{i}=0$

将所有的（x，y）代入上述式子就得到了一个二元一次方程，求解就可以得到两个参数的值。这个方法叫做最小二乘法。

限制

一般线性模型要求因变量 y 要符合正态分布
误差要符合正态分布
方差齐
（x，y）线性相关
如果是多元线性回归，多个自变量之间需不相关

广义线性模型（GLM）

由于一般线性模型的限制较多，无法适应现实中的所有情况，如分类型因变量（动植物是否抗病、药物是否有效），由此提出了广义线性模型。广义线性模型主要包含两种模型：逻辑回归模型和泊松分布模型，前者假设因变量为分类数据（二分类），后者假设因变量为计数类数据。
广义线性模型对于一般线性模型来说，主要是多了一个连接函数，连接函数将因变量和线性预测器。

Lasso回归和岭回归

我们在用实际数据进行分析时，可能会遇到以下两种问题

过拟合, overfitting
欠拟合, underfitting

而 Lasso 回归和 Ridge 回归是在损失函数中分别加入一个约束项，增加模型的泛化能力。
岭回归损失函数如下（第二项为所有参数平方和，即L2范数）:

$L(x,y) = \sum(y_{i}-βx_{i})^{2}+λ\sumβ^{2}$

Lasso 回归损失函数如下（第二项为所有参数绝对值之和，即L1范数）:

$L(x,y) = \sum(y_{i}-βx_{i})^{2}+λ\sum|β|$

上述两个公式中第二项为添加的约束项(正则项)。

则求解β 参数可以用下述表示：

$\hatβ = argmin||y−Xβ||^{2}_{2}+λ||β||^{2}_{2}$ $\hatβ = argmin||y−Xβ||^{2}_{2}+λ||β||_{1}$

范数

L1范数：表示向量x中非零元素的绝对值之和。

$||X||_{1} = \sum|x_{i}|$

L2范数：表示向量元素的平方和再开方。

$||X||_{2} = \sqrt{\sum x_{i}^{2}}$

线性混合模型（LMM）

混合线性模型也被称作为随机效应模型（random effects model）,它假定要分析的数据来自不同总体的层次，这些总体的差异与该层次有关。其假设正好适用于动植物遗传进化研究，不同的个体很有可能具有亲缘关系。模型通式如下：

$y = Xβ+Zμ+ε$

可以看出，其比一般线性模型多了 Zμ 项。其中 β 称作固定效应，μ 称作随机效应，ε 为随机误差。若将上述表达式中的 Zμ+ε 作为一个整体的随机误差项来看，混合模型标准形式就变成传统的线性模型。在求解混合线性模型参数时，不再使用最小二乘法，而是使用最大似然估计（ML）或者约束最大似然估计（REML）。

混合模型方程（MME) 如下:

$\begin{bmatrix} X^{T}R^{-1}X & X^{T}R^{-1}Z \\ Z^{T}R^{-1}X & Z^{T}R^{-1}Z+G^{-1} \end{bmatrix} \begin{bmatrix}β\\μ\end{bmatrix} = \begin{bmatrix}X^{T}R^{-1}Y\\Z^{T}R^{-1}Y\end{bmatrix}$

其中：

$R = Iσ^{2}_{e},G = Iσ^{2}_{μ}$

混合线性模型依然假设因变量为正态分布，所以无法拟合分类型变量。

模型GWAS应用

以下代码均在R语言环境下运行。

library(data.table)
geno = fread("c.raw")
phe = fread("phe.txt")
pca = fread("pca.txt")
plink = fread("pca_cov.txt",header=F,sep=" ")
dd = data.frame(phe = phe$V3,pca1 = pca$V1,pca2 = pca$V1,geno)
mymodel = lm(phe ~ pca1+pca2 + X(单个位点名字),data=dd) #pca1+pca2作为模型协变量
summary(mymodel)

GLM(Logistic回归)

library(data.table)
geno = fread("c.raw",)
phe = fread("phe.txt") #表型为 0-1分类
dd = data.frame(phe$V3,geno)
mymodel = glm(phe$V3 ~ X(单个位点名字) ,family = "binomial",data=dd )

MLM

M <- matrix(rep(0,200*1000),1000,200)
for (i in 1:200) {
M[,i] <- ifelse(runif(1000)<0.5,-1,1)
}
colnames(M) <- 1:200
geno <- data.frame(marker=1:1000,chrom=rep(1,1000),pos=1:1000,M,check.names=FALSE)
QTL <- 100*(1:5) #pick 5 QTL
u <- rep(0,1000) #marker effects
u[QTL] <- 1
g <- as.vector(crossprod(M,u))
h2 <- 0.5
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))
pheno <- data.frame(line=1:200,y=y) #表型后几列可以增加固定效应（性别，种群结构）
kk <-A.mat(geno)
scores <- GWAS(pheno,geno,plot=TRUE,K=kk) #结果第四列为log p值