DALS009-Matrix Algebra(矩阵代数)01-简

2019-08-16 本文已影响0人 backup备份

title: DALS009-Matrix Algebra(矩阵代数)01-简单线性回归
date: 2019-08-16 12:0:00
type: "tags"
tags:

矩阵
categories:
生物统计

前言

这篇笔记是《Data Analysis for the Life Sciences》的第4章：矩阵代数的第1部分。这一部分的主要内容涉及线性方程的简单理解。

简单案例

在这一部分中，我们将会提到3个例子，第1个来源于物理学(physics)，一个来源于遗传学(genetics)，一具是来源于小鼠实验(mouse experiment)。这3个例子明显不同，但是在分析数据中却都用了相同的统计学方法：拟合线性模型(fitting linear models)。在学习线性模型时，通常使用矩阵代数(Matrix Algebra)来进行讲解和描述。

第1案例-自由落体

先来看第1个案例，这个案例是与物理学有关。

试想，你是16世纪的伽里略(Galileo)，你想要描述一个自由落体物体的速度(velocity)。此时，一个助手爬上比萨斜塔，并在上面松开手，落下一个球，同时其他的助手记录一下球在不同时间点的位置。此时，我们就来模拟一下这个过程，在模拟的这个过程里，由于我们已经知道了自由落体运动的公式，因此我们会向这个模拟数据里添加一些检测误差，如下所示：

set.seed(1)
g <- 9.8 ##meters per second
n <- 25
tt <- seq(0,3.4,len=n) ##time in secs, note: we use tt because t is a base function
d <- 56.67 - 0.5*g*tt^2 + rnorm(n,sd=1) ##meters

现在绘制出助手们记录的不同时间点的球的位置，如下所示：

mypar()
plot(tt,d,ylab="Distance in meters",xlab="Time in seconds")

image

假如助手不知道准确的自由落体运行，但是通过观察这个图形，他大致可以推导出这个位置-时间关系应该遵循抛物线，因此可以对上述的数据进行建模，如下所示：

$Y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \varepsilon_i, i=1,\dots,n$

其中， $Y_i$ 代表了球的位置， $x_i$ 代表了时间， $\varepsilon_i$ 表示测量误差。这是一个线性模型，因为这个方程是一个已知量（也就是一系列的 $x$ ，即预测因子(predictors)或协变量(covariates)）和一个未知量（即一系列的 $\beta$ ）的线性组合。

第2案例-父子身高

再来看第2个案例，这个案例与遗传学有关。

试想，你是19世纪的Francis Galton，你收集了很多对父子的身高，你怀疑身高与遗传因素有关，你的数据如下所示：

library(UsingR)
x=father.son$fheight
y=father.son$sheight
plot(x,y,xlab="Father's height",ylab="Son's height")

image

从这个张图上我们大概能觉得，父亲身高在一定程度上影响了儿子的身高，此时我们就可以构建一个模型来描述这种情况，如下所示：
$Y_{i}=\beta_{0}+\beta_{1} x_{i}+\varepsilon, i=1, \ldots, N$

这也是一个线性模型，其中 $x_i$ 表示了第 $i$ 对父亲的身高， $Y_i$ 表示了第 $i$ 对儿子的身高。而 $\varepsilon_i$ 表示其余变量。在这个案例中，我们认为父亲的身高是预测因素，并且是固定的（不是随机的），因此我们使用小写字母来表示父亲的身高。仅仅用测量误差是无法完全解释 $\varepsilon_i$ 代表的变异程度，因为这个模型中并没有包含其它的变量，例如母亲的身高，遗传的随机因素以及环境因素也有可能影响儿子的身高。

第3案例-多个总体的抽样

第3个案例有是关小鼠实验的。还以这本书刚开始的案例说明一下，我们有两类小鼠的体重，这两类小鼠分别使用正常饲料(chow)和高脂饲料(high fat, hf)进行饲喂。现在我们从这两个种群的小鼠中各随机抽取12个样本。我们的研究重点在于，饲料的不同是否会影响小鼠的体重，如下所示：

# library(devtools)
# install_github("genomicsclass/dagdata")
dir <- system.file(package = "dagdata")
list.files(dir)
dir
filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)
# input raw data
head(dat)
str(dat)

数据如下所示：

> list.files(dir)
[1] "data"        "DESCRIPTION" "extdata"     "help"        "html"       
[6] "Meta"        "NAMESPACE"   "script"     
> dir
[1] "C:/Users/20161111/Documents/R/win-library/3.5/dagdata"
> filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
> dat <- read.csv(filename)
> # input raw data
> head(dat)
  Diet Bodyweight
1 chow      21.51
2 chow      28.14
3 chow      24.04
4 chow      23.45
5 chow      23.68
6 chow      19.79
> str(dat)
'data.frame':   24 obs. of  2 variables:
 $ Diet      : Factor w/ 2 levels "chow","hf": 1 1 1 1 1 1 1 1 1 1 ...
 $ Bodyweight: num  21.5 28.1 24 23.4 23.7 ...

查看一下小鼠的体重抽样结果，如下所示：

stripchart(Bodyweight~Diet,data=dat,vertical=TRUE,method="jitter",pch=1,main="Mi\
ce weights")

image

假如我们想要判断一下这饮食是否对小鼠的平均体重有所影响。我们可以使用t检验和置信区间来进行统计。同时我们也可以使用一个线性模型来进行精确的检验，如下所示：
$Y_{i}=\beta_{0}+\beta_{1} x_{i}+\varepsilon_{i}$
其中 $\beta_0$ 表示了正常饲料饲喂的小鼠平均体重，而 $\beta_1$ 则表示了高脂饲料(hf)饲喂的小鼠平均体重与正常饲料饲喂小鼠体重的差值，其中 $x_i = 1$ 表示第 $i$ 只小鼠吃的hf饲料， $x_i = 0$ 表示小鼠吃的是正常饲料，而 $\varepsilon_i$ 则是说明，源于相同总体的小鼠个体差异。

线性模型的一般形式

通过前面的3个案例我们知道了线性模型，包含了上述统计模型的线性模型的一般形式可以写为如下公式：
$Y_{i}=\beta_{0}+\beta_{1} x_{i, 1}+\beta_{2} x_{i, 2}+\cdots+\beta_{2} x_{i, p}+\varepsilon_{i}, i=1, \ldots, n$
即：
$Y_{i}=\beta_{0}+\sum_{j=1}^{p} \beta_{j} x_{i, j}+\varepsilon_{i}, i=1, \ldots, n$
在这个公式里，我们使用了一个预测量 $p$ ，矩阵代数提供了一种精练的语言和数学框架来计算以及推导那些满足上述框架的任何线性模型。

估计参数(Estimating parameters)

当我们用上述的模型来估计未知量 $\beta$ 时非常有用。例如，在第1个案例中，我们通过线性模型来计算其中的未知的参数。在第2个案例中，通过估计其中的未知参数，我们可以更好地理解遗传因素(inheritance)在父子身高方面的影响。在第3个案例中，我们想研究饲料是否对小鼠的体重有影响，也就是说我们要计算出 $\beta_{1} \neq 0$ 。

在科学研究中，标准的做法就是找出能拟合这些数据的线性模型中的那些数值（也就是上面的 $\beta_{j}$ ），找到的这个线性方向与这些点的距离最小。这个过程通常是通过最小二乘(Least Squares, LS)的思想来实现的，下面的这个公式被称为最小平方(Least squares)方程：
$\sum_{i=1}^{n}\left\{Y_{i}-\left(\beta_{0}+\sum_{j=1}^{p} \beta_{j} x_{i, j}\right)\right\}^{2}$

一旦我们找到了这些值，我们就称这些值为最小二乘估计(Least Squares Estimates, LSE)，并将其命名为 $\hat{\beta}$ 。当我们在估计值处计算出了最小二乘方程后，上面的那个公式就称为残差平方和(Residual Sum of Squares, RSS)。由于RSS所有的量都是取决于Y，因此这个值是随机变量。

自由落体案例回顾

在高中物理课本中我们就知道，自由落体运行的方程为：
$d=h_{0}+v_{0} t-0.5 \times 9.8 t^{2}$
其中， $h_{0}$ 和 $v_{0}$ 分别是起始的高度与速度。在我们模拟的自由落体运行方程中， $v_{0}=0$ ， $h_{0}=56.67$ ，R代码如下所示：

g <- 9.8 ##meters per second
n <- 25
tt <- seq(0,3.4,len=n) ##time in secs, t is a base function
f <- 56.67 - 0.5*g*tt^2
y <- f + rnorm(n,sd=1)

png(file="../Figs/falling.png", width=1200, height=1200, res=300,
    pointsize=10)
plot(tt,y,ylab="Distance in meters",xlab="Time in seconds")
lines(tt,f,col=2)
dev.off()

我们模拟的自由落体运行图形如下所示：

image

但是，当我们假设自己是伽利略时，这个方程还没有被推导出来，因此我们并不知道其中相应的参数，这些数据只是表现得像一条抛物线，因此我们可以对这些数据进行建模，如下所示：
$Y_{i}=\beta_{0}+\beta_{1} x_{i}+\beta_{2} x_{i}^{2}+\varepsilon, i=1, \ldots, n$
此时，我们如何找到LSE？

lm()函数

R中的lm()函数可以很好地计算线性方程，现在我们先来简单计算一下这个方程，如下所示：

g <- 9.8 ##meters per second
tt <- seq(0,3.4,len=n) ##time in secs, t is a base function
y <- f + rnorm(n,sd=1)
tt2 <-tt^2
fit <- lm(y~tt+tt2)
summary(fit)$coef

计算结果如下所示：

> summary(fit)$coef
              Estimate Std. Error     t value     Pr(>|t|)
(Intercept) 57.1399311  0.5484054 104.1928697 3.897090e-31
tt          -0.1997275  0.7470438  -0.2673572 7.916846e-01
tt2         -4.9193723  0.2122243 -23.1800643 5.973048e-17

这样就计算出了LSE，以及标准误，p值。

最小二乘估计(The Least squares estimate, LSE)

现在我们写一个函数，这个函数的功能在于计算任任 $\beta$ 的RSS，如下所示：

rss <- function(Beta0,Beta1,Beta2){
    r <- y - (Beta0+Beta1*tt+Beta2*tt^2)
    return(sum(r^2))
}

所以对于任何三维向量，我们得到一个RSS。下图是当我们保持另外两个 $\beta$ 固定时，关于 $\beta_{2}$ 的函数的RSS的曲线，如下所示：

rss <- function(Beta0,Beta1,Beta2){
    r <- y - (Beta0+Beta1*tt+Beta2*tt^2)
    return(sum(r^2))
}

Beta2s<- seq(-10,0,len=100)
plot(Beta2s,sapply(Beta2s,rss,Beta0=55,Beta1=0),
     ylab="RSS",xlab="Beta2",type="l")
##Let's add another curve fixing another pair:
Beta2s<- seq(-10,0,len=100)
lines(Beta2s,sapply(Beta2s,rss,Beta0=65,Beta1=0),col=2)

image

反复计算在这里行不通，相反，我们会使用微积分(calculus)的方法来进行计算，通过取偏导数(partial derivatives)，将它们设为0，然后求解。当然，如果我们有很多参数，这些方程会变得相当复杂。线性代数为这个问题提供了精练且通用的方法。

Galton拓展

Galton在研究父子身高问题时，他通过探索性数据分析有了一个新颖的发现，如下所示：

image

Galton发现，如果他将父子的身高数据制成表格，然后将那些x和y值的和都相同的点连接起来，就构成了一个椭圆。通过计算还会进一步发现，这些数据服从二元正态分布（具体的高尔顿图可以参考《于忠义, YUZhong-yi. 高尔顿发现相关与回归的历史回顾与反思[J]. 统计与信息论坛, 2009, 24(9):17-25.》）如下所示：
$\operatorname{Pr}(X<a, Y<b)=\\ \int_{-\infty}^{a} \int_{-\infty}^{b} \frac{1}{2 \pi \sigma_{x} \sigma_{y} \sqrt{1-\rho^{2}}} \exp \left\{\frac{1}{2\left(1-\rho^{2}\right)}\left[\left(\frac{x-\mu_{x}}{\sigma_{x}}\right)^{2}-2 \rho\left(\frac{x-\mu_{x}}{\sigma_{x}}\right)\left(\frac{y-\mu_{y}}{\sigma_{y}}\right)+\left(\frac{y-\mu_{y}}{\sigma_{y}}\right)^{2}\right]\right\}$
从上面的公式我们知道，当我们保持 $X$ 不变时（也就是设定好 $x$ ）， $Y$ 的分布是服从正态分布的，即服从均值为 $\mu_{x}+\sigma_{y} \rho\left(\frac{x-\mu_{x}}{\sigma_{x}}\right)$ ，标准差为 $\sigma_{y} \sqrt{1-\rho^{2}}$ 的正态分布。 $\rho$ 是 $X$ 和 $Y$ 的相关系数，这就说明，如果我们设定 $X=x$ ， $Y$ 实际上是一个线性模型。在我们的这个简单线性模型中， $\beta_{0}$ 和 $\beta_{1}$ 参数表示为 $\mu_{x}, \mu_{y}, \sigma_{x}, \sigma_{y}$ 和 $\rho$ 。

练习

P152

参考资料

Rafael A Irizarry and Michael I Love.Data Analysis for the Life Sciences.This book is for sale at http://leanpub.com/dataanalysisforthelifesciences
于忠义, YUZhong-yi. 高尔顿发现相关与回归的历史回顾与反思[J]. 统计与信息论坛, 2009, 24(9):17-25.

DALS009-Matrix Algebra(矩阵代数)01-简

前言

简单案例

第1案例-自由落体

第2案例-父子身高

第3案例-多个总体的抽样

线性模型的一般形式

估计参数(Estimating parameters)

自由落体案例回顾

lm()函数

最小二乘估计(The Least squares estimate, LSE)

Galton拓展

练习

参考资料

猜你喜欢

热点阅读