环境与生态统计||线性模型

在客观世界中普遍存在着变量之间的关系，一般说来可以分为确定性的因果关系以及非确定性的相关关系。回归分析是研究相关关系的一种数学工具，它能帮助我们从一个变量取得的值去估计另一个变量所取得的值。

回归的典故

（达尔文的表兄）高尔顿最著名的发现之一是他发现了父亲的身高和儿子的身高之间存在着某种给定的关系，他通过进一步的研究发现了：事实上子辈的平均身高是其父辈平均身高以及他们所处族群平均身高的加权平均和。

他把这种趋势平均化的现象写到了自己1886年的论文中。论文的全名叫：Regression towards Mediocrity in Hereditary Stature. 这篇论文当年被发在了大不列颠以及爱尔兰人类研究学院期刊上。我们现今把论文中的这种回归现象称为：均值回归或者平庸回归（reversion to the mean/reversion to mediocrity）。背后的意义是说：哪怕单看一组父亲和孩子的身高，两个人的身高可能差异很大，但是从整个人群上来看，父亲和孩子的身高分布应该是很相近的。

什么叫线性模型

在统计学中，线性回归（Linear regression）是利用称为线性回归方程的最小二乘函数对一个或多个自变量和因变量之间关系进行建模的一种回归分析。这种函数是一个或多个称为回归系数的模型参数的线性组合。只有一个自变量的情况称为简单回归，大于一个自变量情况的叫做多元回归。

在线性回归中，数据使用线性预测函数来建模，并且未知的模型参数也是通过数据来估计。这些模型被叫做线性模型。最常用的线性回归建模是给定X值的y的条件均值是X的仿射函数。不太一般的情况，线性回归模型可以是一个中位数或一些其他的给定X的条件下y的条件分布的分位数作为X的线性函数表示。像所有形式的回归分析一样，线性回归也把焦点放在给定X值的y的条件概率分布，而不是X和y的联合概率分布（多元分析领域）。

给一个随机样本（Y_i,Y_{i1},Y_{i2},...Y_{in}）,i=1,...n，一个线性回归模型假设回归子Y_i和回归量X_i,X_{i1},X_{i2},...X_{ip}之间的关系是除了X的影响以外，还有其他的变量存在。我们加入一个误差项\epsilon_i（也是一个随机变量）来捕获除了 {\displaystyle X_{i1},\ldots ,X_{ip}}之外任何对{\displaystyle Y_{i}}的影响。所以一个多变量线性回归模型表示为以下的形式：

{\displaystyle Y_{i}=\beta _{0}+\beta _{1}X_{i1}+\beta _{2}X_{i2}+\ldots +\beta _{p}X_{ip}+\varepsilon _{i},\qquad i=1,\ldots ,n}

其他的模型可能被认定成非线性模型。一个线性回归模型不需要是自变量的线性函数。线性在这里表示{\displaystyle Y_{i}}的条件均值在参数{\displaystyle \beta }里是线性的。例如：模型{\displaystyle Y_{i}=\beta _{1}X_{i}+\beta _{2}X_{i}^{2}+\varepsilon _{i}} 在{\displaystyle \beta _{1}}和 {\displaystyle \beta _2} 里是线性的，但在{\displaystyle X_{i}^{2}}里是非线性的，它是{\displaystyle X_{i}}的非线性函数。

t检验与一元线性回归的等价性

还拿我们总磷（TP）浓度来比较两个总体的t检验。

source("D:\\Users\\Administrator\\Desktop\\RStudio\\Environmental_and_Ecological_Statistics_with_R\\DataCode\\Code\\FrontMatter.R")


### linear models ####
##################################
## Chapter 5 linear regression  ##
##################################
plotDIRch5 <- paste(plotDIR, "chapter5", "figures", sep="/")


##t.test  and linear models 
TP.reference<-read.csv("..//Data//WCA2TP.csv",header = T)
TP.reference$Month<-sapply(as.character(TP.reference$Date),function(x){strsplit(x,"/")[[1]][1]})

subI<-(TP.reference$SITE=='E4'|TP.reference$SITE=='F4')&TP.reference$Year==1994

subR<-TP.reference$Year==1994&TP.reference$SITE!='E5'&TP.reference$SITE!='F5'
x<-log(TP.reference$RESULT[subI])
y<-log(TP.reference$RESULT[subR])
t.test(x,y,alternative = "greater",var.equal =T) 



    Two Sample t-test

data:  x and y
t = -2.599, df = 157, p-value = 0.9949
alternative hypothesis: true difference in means is greater than 0
95 percent confidence interval:
 -1.102025       Inf
sample estimates:
mean of x mean of y 
-4.016832 -3.343483

同样地，我们可以用公式界面来表示t检验：

two.sample<-data.frame(TP=c(TP.reference$RESULT[subI],TP.reference$RESULT[subR]),
                       x=c(rep("I",sum(subI)),
                           rep("R",sum(subR))
                           ))
t.test(log(TP)~x,data = two.sample,,var.equal =T)

    Two Sample t-test

data:  log(TP) by x
t = -2.599, df = 157, p-value = 0.01024
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.1850781 -0.1616198
sample estimates:
mean in group I mean in group R 
      -4.016832       -3.343483

公式y \sim x在R中被约定为一个模型，符号\sim将响应变量与预测变量（解释变量）分开。在t检验中，函数t.test知道预测量是用来将数据分为两组的。但是，我们也会把公式的应用解释为把估计出的均值看做某种预测。例如我们可以利用估计出的均值预测参考站点未来的一个样本。因此双样本t检验问题可以看作是一个一元线性回归模型，其形式为：

y=\beta _{0}+\beta _{1}x+\varepsilon _{0}

模型系数\beta _{0}就是估计出的参照站点TP浓度的对数均值（-3.343483 ），也就是说，y=\beta _{0}+\beta _{1}*0+\varepsilon _{0}。假设\varepsilon \sim N(0,\sigma) ,x=0（参考站点）就等价于y \sim N(\beta _{0},\sigma)。当=1（受影响的站点），模型等价于y \sim N(\beta _{0}+\beta _{1},\sigma)。因此，\beta _{1}是受影响站点与参考站点对数均值的差值(-3.343483-（-4.0168）=0.673317)。

two.sample<-data.frame(y=c(TP.reference$RESULT[subI],TP.reference$RESULT[subR]),
                       x=c(rep(0,sum(subI)),
                           rep(1,sum(subR))
                       ))
t2lm<-lm(log(y)~x,data = two.sample)
summary(t2lm)

Call:
lm(formula = log(y) ~ x, data = two.sample)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.1780 -0.8562  0.0466  0.7262  3.5690 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -4.0168     0.2414 -16.642   <2e-16 ***
x             0.6733     0.2591   2.599   0.0102 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.106 on 157 degrees of freedom
Multiple R-squared:  0.04125,   Adjusted R-squared:  0.03514 
F-statistic: 6.755 on 1 and 157 DF,  p-value: 0.01024

正如我们预料的那样截距（Intercept，\beta _{0}）为 -4.0168，\beta _{1}为0.6733。双样本t检验问题只是一个特例，当预测变量为连续变量时，我们就看到了熟悉的线性回归问题。

作为线性模型的ANOVA

ANOVA模型中的计算工作包括各组均值\hat{\upsilon_i}的计算。假设我们对Everglandes湿地中的TP浓度的时间变化感兴趣，想看看采样6年中每年的浓度均值是否存在差异。为方便起见，我们假定每年的浓度是相互独立的。于是我们可以这样来构建模型：

anvoa.data<-data.frame(y=TP.reference$RESULT,
                       x2=ifelse(TP.reference$Year==1995,1,0),
                       x3=ifelse(TP.reference$Year==1996,1,0),
                       x4=ifelse(TP.reference$Year==1997,1,0),
                       x5=ifelse(TP.reference$Year==1998,1,0),
                       x6=ifelse(TP.reference$Year==1999,1,0))

anova.lm<-lm(log(y)~x2+x3+x4+x5+x6,data=anvoa.data)
summary(anova.lm)

Call:
lm(formula = log(y) ~ x2 + x3 + x4 + x5 + x6, data = anvoa.data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0171 -0.9545 -0.0989  0.8610  4.5633 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.50436    0.08894 -39.399  < 2e-16 ***
x2          -0.36943    0.11170  -3.307 0.000968 ***
x3          -0.16319    0.11478  -1.422 0.155345    
x4          -0.17488    0.11896  -1.470 0.141811    
x5          -0.27044    0.11323  -2.389 0.017062 *  
x6          -0.15415    0.12500  -1.233 0.217751    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.114 on 1272 degrees of freedom
Multiple R-squared:  0.01004,   Adjusted R-squared:  0.00615 
F-statistic:  2.58 on 5 and 1272 DF,  p-value: 0.0248

标记为Intercept的系数就是\beta _{0}的估计值。换个做法我们可以将变量Year换成一个因子变量，R将自动创建这些预测变量。

anova.lm<-lm(log(RESULT)~factor(Year),data=TP.reference)
summary(anova.lm)
Call:
lm(formula = log(RESULT) ~ factor(Year), data = TP.reference)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0171 -0.9545 -0.0989  0.8610  4.5633 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -3.50436    0.08894 -39.399  < 2e-16 ***
factor(Year)1995 -0.36943    0.11170  -3.307 0.000968 ***
factor(Year)1996 -0.16319    0.11478  -1.422 0.155345    
factor(Year)1997 -0.17488    0.11896  -1.470 0.141811    
factor(Year)1998 -0.27044    0.11323  -2.389 0.017062 *  
factor(Year)1999 -0.15415    0.12500  -1.233 0.217751    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.114 on 1272 degrees of freedom
Multiple R-squared:  0.01004,   Adjusted R-squared:  0.00615 
F-statistic:  2.58 on 5 and 1272 DF,  p-value: 0.0248

要查看检验的结果，可以：

anova(anova.lm)

Analysis of Variance Table

Response: log(RESULT)
               Df  Sum Sq Mean Sq F value Pr(>F)  
factor(Year)    5   16.03  3.2051  2.5805 0.0248 *
Residuals    1272 1579.88  1.2420                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

以上R输出只给出了1994年的均值，即因子预测量的第一个等级。要强制输出各个年份的均值，我们可以通过在线性模型公式右侧加一个“-1”来不让R去拟合intercept。

> anova.lm<-lm(log(RESULT)~factor(Year)-1,data=TP.reference)
> summary(anova.lm)

Call:
lm(formula = log(RESULT) ~ factor(Year) - 1, data = TP.reference)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.0171 -0.9545 -0.0989  0.8610  4.5633 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
factor(Year)1994 -3.50436    0.08894  -39.40   <2e-16 ***
factor(Year)1995 -3.87379    0.06757  -57.33   <2e-16 ***
factor(Year)1996 -3.66754    0.07255  -50.55   <2e-16 ***
factor(Year)1997 -3.67923    0.07900  -46.57   <2e-16 ***
factor(Year)1998 -3.77480    0.07007  -53.88   <2e-16 ***
factor(Year)1999 -3.65850    0.08783  -41.65   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.114 on 1272 degrees of freedom
Multiple R-squared:  0.9178,    Adjusted R-squared:  0.9174 
F-statistic:  2367 on 6 and 1272 DF,  p-value: < 2.2e-16

> anova(anova.lm)
Analysis of Variance Table

Response: log(RESULT)
               Df  Sum Sq Mean Sq F value    Pr(>F)    
factor(Year)    6 17637.9 2939.65  2366.8 < 2.2e-16 ***
Residuals    1272  1579.9    1.24                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可见，公式把第一个变量作为一个基准来对待了。

简单和多元线性回归模型

模型已定，其实主要的工作就是估计模型系数，最小平方法和最大似然法是两种最常用的方法。我们通过密歇根胡中鱼体内多氯联苯（Polychlorinated biphenyls; Polychlorodiphenyls别名：氯化联苯，PCB）浓度趋势的例子来阐述。分析鱼体内PCB的目的在于：A评估鱼体内的PCB浓度随时间的变化趋势，以确定有意义的浓度降低是否还在发生；B为鱼类消费顾问提供依据以警示公众食用被污染的鱼可能存在的风险。

为了实现这目的，采用了线性（对数线性）模型。不同年份间鲑鱼体内的PCB浓度逐年下降。

## PCB in fish

##### Lake trout ####
laketrout <- read.csv(paste(dataDIR,"laketrout2.csv", sep="/"),
                      header=T)
head(laketrout)
## dividing by 2.54 converts 60 cm to inches
## single predictor
laketrout$length <- laketrout$length*2.54
laketrout$size<- "small"
laketrout$size[laketrout$length>60] <- "large"

laketrout$large<- 0
laketrout$large[laketrout$length>60] <- 1

laketrout$len.c <- laketrout$length-mean(laketrout$length, na.rm=T)
laketrout$inches<-laketrout$length/2.54

laketrout <- laketrout[laketrout$pcb>exp(-2)&laketrout$length>0,]

head(laketrout)
  length  pcb n    lgpcb year y  size large     len.c inches
1 75.946 31.3 1 3.443618 1974 0 large     1 13.432539   29.9
2 74.930  7.9 1 2.066863 1974 0 large     1 12.416539   29.5
3 68.580 26.7 1 3.284664 1974 0 large     1  6.066539   27.0
4 66.802  8.3 1 2.116256 1974 0 large     1  4.288539   26.3
5 69.596 11.3 1 2.424803 1974 0 large     1  7.082539   27.4
6 69.088 12.6 1 2.533697 1974 0 large     1  6.574539   27.2
## 60 cm is a threshold,
##
(pcb.loess <- loess.smooth (y=log(laketrout$pcb), x=laketrout$year, degree=1, span=0.75))
par(mar=c(3,3,1,1), mgp=c(1.75,0.125,0), tck=0.01,las=1)
plot(pcb ~ year, data=laketrout, log="y", ylab="PCB Concentrations (mg/kg)",
     xlab="Year")
lines(exp(pcb.loess$y)~pcb.loess$x, lwd=2, col="gray")
dev.off()

因为尺寸大一些的鱼往往年龄也比较大，所以PCB浓度与鱼的尺寸的关系如下：

(pcb.loess <- loess.smooth(y=log(laketrout$pcb), x=laketrout$length, degree=1, span=0.75))
par(mar=c(3,3,1,1), mgp=c(1.75,0.125,0), tck=0.01,las=1)
plot(pcb ~ length, data=laketrout, log="y", ylab="PCB Concentrations (mg/kg)",
     xlab="Fish Length (cm)")
lines(exp(pcb.loess$y)~pcb.loess$x, lwd=2, col="gray")
dev.off()
用一个预测变量来回归

用于评价时间变化趋势的简单线性回归模型是对数线性模型：

\log PCB=\beta _{0}+\beta _{1}Year+\varepsilon _{0}

可以用R的函数lm()来实现：

> ## simple linear regression
> lake.lm1 <- lm(log(pcb) ~ year, data=laketrout)
> display(lake.lm1, 4)
lm(formula = log(pcb) ~ year, data = laketrout)
            coef.est coef.se 
(Intercept) 119.8467  10.9689
year         -0.0599   0.0055
---
n = 631, k = 2
residual sd = 0.8784, R-Squared = 0.16
> summary(lake.lm1)

Call:
lm(formula = log(pcb) ~ year, data = laketrout)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.60801 -0.65501 -0.01038  0.57926  2.69862 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 119.846680  10.968944   10.93   <2e-16 ***
year         -0.059903   0.005525  -10.84   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8784 on 629 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.1575,    Adjusted R-squared:  0.1561 
F-statistic: 117.6 on 1 and 629 DF,  p-value: < 2.2e-16

> (lm1.coef<-coef(lake.lm1))
 (Intercept)         year 
119.84667967  -0.05990328 
> anova(lake.lm1)
Analysis of Variance Table

Response: log(pcb)
           Df Sum Sq Mean Sq F value    Pr(>F)    
year        1  90.71  90.712  117.57 < 2.2e-16 ***
Residuals 629 485.30   0.772                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

估计出的\beta _{0}(截距，Intercept）是119.846680为预测变量为零时响应变量的期望值，\beta _{1}（斜率）是-0.059903 为每年单位变化对应的PCB对数浓度变化。拟合后的模型分为两部分，确定的部分是\beta _{1}+\beta _{1}Year，是任意给定年份的PCB对数期望或均值。随机部分是\varepsilon _{0}描述的是波动或者不确定性，描述了个体的变异。读者朋友可以尝试plot(lake.lm1)来查看拟合结果残差的响应。

给出拟合线的置信区间值

par(mar=c(3,3,0.25,0.25), mgp=c(1.25,.125,0), las=1, tck=0.01)
plot(pcb ~ year, data=laketrout, log="y",
     ylab="PCB (mg/kg)", xlab="Year", cex=0.5, col=grey(0.5))
curve(exp(lm1.coef[1] + lm1.coef[2]*(x) +0.87^2/2), add=T, lwd=2)
curve(qlnorm(0.025, lm1.coef[1] + lm1.coef[2]*(x), 0.87 ), add=T, lty=2)
curve(qlnorm(0.975, lm1.coef[1] + lm1.coef[2]*(x), 0.87 ), add=T, lty=2)
dev.off()
多元回归
> lake.lm2 <- lm(log(pcb) ~ length+I(year-1974), data=laketrout)
> display(lake.lm2, 3)
lm(formula = log(pcb) ~ length + I(year - 1974), data = laketrout)
               coef.est coef.se
(Intercept)    -1.834    0.120 
length          0.060    0.002 
I(year - 1974) -0.086    0.004 
---
n = 631, k = 3
residual sd = 0.555, R-Squared = 0.66
> summary(lake.lm2)

Call:
lm(formula = log(pcb) ~ length + I(year - 1974), data = laketrout)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.57511 -0.33308  0.01208  0.34307  1.97533 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -1.834251   0.120335  -15.24   <2e-16 ***
length          0.059608   0.001934   30.83   <2e-16 ***
I(year - 1974) -0.086193   0.003590  -24.01   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5545 on 628 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.6648,    Adjusted R-squared:  0.6637 
F-statistic: 622.7 on 2 and 628 DF,  p-value: < 2.2e-16

> (lm2.coef<-coef(lake.lm2))
   (Intercept)         length I(year - 1974) 
   -1.83425076     0.05960843    -0.08619336 
> anova(lake.lm2)
Analysis of Variance Table

Response: log(pcb)
                Df Sum Sq Mean Sq F value    Pr(>F)    
length           1 205.70 205.703  668.98 < 2.2e-16 ***
I(year - 1974)   1 177.21 177.210  576.32 < 2.2e-16 ***
Residuals      628 193.10   0.307                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> plot(lake.lm2)

为了是鱼的长度为零时公式依然是有生态学意义的，需要对长度做一下线性变化。

(mn.length <- mean(laketrout$length, na.rm=T))
(laketrout$len.c <- laketrout$length-mean(laketrout$length, na.rm=T))
lake.lm3 <- lm(log(pcb) ~ len.c+I(year-1974), data=laketrout)
display(lake.lm3, 3)
(lm3.coef<-coef(lake.lm3))
summary(lake.lm3)
anova(lake.lm3)
plot(lake.lm3)

par(mar=c(3,3,0.25,0.25), mgp=c(1.25,.125,0), las=1, tck=0.01)
plot(pcb ~ year, data=laketrout,log="y",
     ylab="PCB (mg/kg)", xlab="Year", cex=0.5, col=grey(0.5))
curve(exp( lm3.coef[1] + lm3.coef[2]*0 + lm3.coef[3]*(x-1974) +0.543^2/2),
      add=T, lwd=2)
curve(exp( lm3.coef[1] + lm3.coef[2]*(-2.6)+lm3.coef[3]*(x-1974) +0.543^2/2),
      add=T, lty=2, lwd=2)
curve(exp( lm3.coef[1] + lm3.coef[2]*(3.4) + lm3.coef[3]*(x-1974) +0.543^2/2),
      add=T, lty="13", lwd=2)
legend(x=1990, y=40, legend=c("large fish","average fish","small fish"),
       lty=c(3, 1, 2), lwd=rep(2, 3), cex=0.75, bty="n")
dev.off()
交互作用

当我们用Year和len.c作为预测变量拟合多元回归模型时，一个非常重要的假设就是年份的影响（年份的斜率）不受鱼大小的影响，并且鱼大小的影响在研究时段内是相同的。也就是说这两个预测变量是独立的。这对于多元回归模型是一个加和效应的假设。我们知道鱼的大小会随着时间的变化而变化（因为鱼会长大的嘛），也就是鱼的大小会随着时间的变化而变化，两者并不是独立的变量，这时就要用到两者的相互作用来拟合了。在模型中加入第三个变量，Year和len.c的乘积。

 lake.lm4 <- lm(log(pcb) ~ len.c*I(year-1974), data=laketrout)
> display(lake.lm4, 4)
lm(formula = log(pcb) ~ len.c * I(year - 1974), data = laketrout)
                     coef.est coef.se
(Intercept)           1.8907   0.0465
len.c                 0.0510   0.0038
I(year - 1974)       -0.0874   0.0036
len.c:I(year - 1974)  0.0008   0.0003
---
n = 631, k = 4
residual sd = 0.5520, R-Squared = 0.67
> (lm4.coef<-coef(lake.lm4))
         (Intercept)                len.c       I(year - 1974) len.c:I(year - 1974) 
        1.8907181316         0.0510365746        -0.0873925026         0.0008475978 

> anova(lake.lm4)
Analysis of Variance Table

Response: log(pcb)
                      Df  Sum Sq Mean Sq  F value  Pr(>F)    
len.c                  1 205.703 205.703 675.0042 < 2e-16 ***
I(year - 1974)         1 177.210 177.210 581.5055 < 2e-16 ***
len.c:I(year - 1974)   1   2.027   2.027   6.6523 0.01013 *  
Residuals            627 191.074   0.305                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> plot(lake.lm4)

加入交互作用项len.c:I(year - 1974) 后，模型可以描述为：

\log PCB = 1.8907181-0.0873925 Year + 0.0510366len.c + 0.0008476 Year*len.c+ \varepsilon _{0}（全模型公式）

由于乘积项的加入，模型不再是线性的。经居中调整后的长度（len.c）的斜率和年份（Year）的斜率不再是常数。尽管如此，人们还是把有交互作用的回归分析放在线性回归中去讲解。

残差和模型评估

对模型的解释是很容易的，但是就模型的形式提出问题（如，线性模型准确吗？）残差分析，即模型预测值与观测值之间差异，是回答模型合适数据这个问题的最有效的方法。在对数据和可能影响鱼体内PCB含量的因素进行初步检查后拟合出了全模型公式。使用了这个对数对数线性模型就意味着我们假定PCB浓度随时间降低是指数式的，并且认为湖中鱼尺寸增加一个单位时体内PCB浓度增加值相同。但是这些假设并不具备理论的支撑，如何在数据的基础上对这些进行假设检验呢？要回答这个问题就要明白我们建模型的目的：因果推断和预测。

好的预测模型应该是简单而准确的。用R函数lm()拟合好一个模型后，所有必须的模型总结和诊断信息都在R的结果对象中。要对模型做总体评估，我们常常用R^2和一个假设检验（F检验）来比较拟合后的模型和一个没有预测变量的模型（y=\beta_{0}+\varepsilon _{0}）。要评估每个单独的预测变量是否有存在的必要，就用t检验来评估模型变量的斜率是否与零有差异。如以上我们见到的summary()函数可以完成这项工作。

> summary(lake.lm4)

Call:
lm(formula = log(pcb) ~ len.c * I(year - 1974), data = laketrout)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.47961 -0.34105  0.01972  0.33870  1.97111 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.8907181  0.0464646  40.692   <2e-16 ***
len.c                 0.0510366  0.0038407  13.288   <2e-16 ***
I(year - 1974)       -0.0873925  0.0036045 -24.246   <2e-16 ***
len.c:I(year - 1974)  0.0008476  0.0003286   2.579   0.0101 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.552 on 627 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.6683,    Adjusted R-squared:  0.6667 
F-statistic: 421.1 on 3 and 627 DF,  p-value: < 2.2e-16

其中R_{adj}^2=1- \frac{(1-R^2)(n-1)}{n-p-1}，n为样本数，p为预测变量的个数。因为在模型中，增加多个变量，即使事实上无关的变量，也会小幅度条R平方的值，当时其是无意义，所以我们调整了下，降低R平方的值。用r square的时候，不断添加变量能让模型的效果提升，而这种提升可能是虚假的。利用adjusted r square，能对添加的非显著变量给出惩罚，也就是说随意添加一个变量不一定能让模型拟合度上升。那么被解释的变化量R^2是否具有统计显著性呢？就看下面一行的F-statistic或者 p-value。如果R_{adj}^2是一个负值，表明预测变量的解释能力还不如零模型（比如随机生成的正态分布）的解释变量。

进行模型比较时，ANOVA是一种分离总方差的重要技术。它可以用来比较具有少量预测变量的模型（简化模型）和拥有大量预测变量的模型（全模型）。不同模型比较用于推断一个预测变量是否应该包含在模型中。例如我们可以仅靠统计学来确定是否要把length加入第二个预测变量。那就是比较模型：

\log (pcb) \sim I(year-1974)

和模型：

\log (pcb) \sim I(year-1974)+ len.c

简单的方差分析为：

> summary.aov(lake.lm1)
             Df Sum Sq Mean Sq F value Pr(>F)    
year          1   90.7   90.71   117.6 <2e-16 ***
Residuals   629  485.3    0.77                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
15 observations deleted due to missingness
> summary.aov(lake.lm2)
                Df Sum Sq Mean Sq F value Pr(>F)    
length           1  205.7  205.70   669.0 <2e-16 ***
I(year - 1974)   1  177.2  177.21   576.3 <2e-16 ***
Residuals      628  193.1    0.31                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
15 observations deleted due to missingness

在模型lake.lm1中，有485.3 的残差平方和没有被预测变量year所解释，模型lake.lm2加入第二个预测变量len.c后，len.c解释了其中的292，从而获得很大的F值和很小的p值，说明加入len.c可以很大程度改善模型。

查看每个预测值的解释量（承担的方差）
library(hier.part)

mypch1.lm<-lm(formula=log(pcb) ~ len.c+year+ large, data=laketrout,na.action=na.omit)
summary(mypch1.lm)

Call:
lm(formula = log(pcb) ~ len.c + year + large, data = laketrout, 
    na.action = na.omit)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5819 -0.3416  0.0097  0.3409  1.9737 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 172.100794   7.134927  24.121   <2e-16 ***
len.c         0.060671   0.003212  18.888   <2e-16 ***
year         -0.086215   0.003593 -23.994   <2e-16 ***
large        -0.032229   0.077795  -0.414    0.679    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5549 on 627 degrees of freedom
  (15 observations deleted due to missingness)
Multiple R-squared:  0.6649,    Adjusted R-squared:  0.6633 
F-statistic: 414.6 on 3 and 627 DF,  p-value: < 2.2e-16


env1<-laketrout[, c("len.c","year","large"), drop = F]
colnames(env1)
hier.p<-hier.part(log(laketrout$pcb),env1,fam = "gaussian", gof = "Rsqu")
hier.p$IJ

              I           J     Total
len.c 0.2881845  0.06892997 0.3571145
year  0.2468852 -0.08940342 0.1574817
large 0.1297850  0.10114237 0.2309274

hier.p$IJ对应的方差之和等于mypch1.lm的r.squared

> sum(hier.p$IJ[,1])
[1] 0.6648547
> summary(mypch1.lm)$r.squared
[1] 0.6648547
残差的图形分析

对于汇总统计（尤其是假设检验的结果），残差分析应该是独立的，并且近似服从均值为0，标准差为常数的正态分布：\varepsilon_{i}\sim N(0,\sigma)。残差的图形分布不可避免地应该加入到三个阶段：正态性、独立性、常数标准差的评估中。在R中很容易实现，就是我们刚才提到的：plot(lake.lm1)。

顺时针来看。第一个图是独立性检验的模型残差和拟合值做的图，暗示预测值均偏低；第二个是正态性检验的残差Q-Q图；第三个是检测模型可能存在潜在影响的杠杆点图。数据点的Cook距离用来度量数据点的影响。Cook距离大于1就认为有影响的；第四个是检验标准差为常数的S-L图，即对残差绝对值的平方根和拟合值做散点图。

另外通过残差拟合展形图（residual - fitted - spread）我们可以看出m模型拟合值的覆盖范围与随机误差（残差）的覆盖范围重合程度。R^2测量的方差，rfs给出的是由模型所给出的相对展形。尽管R^2表明模型解释了一定的总变化，但该图形表明拟合值与残差的覆盖范围重叠情况。

rfs(lake.lm1, aspect=1)

非线性尝试：

 (lake.lm5 <- lm(log(pcb) ~ I(year-1974)*len.c + I(len.c^2), 
+                data=laketrout))

Call:
lm(formula = log(pcb) ~ I(year - 1974) * len.c + I(len.c^2), 
    data = laketrout)

Coefficients:
         (Intercept)        I(year - 1974)                 len.c            I(len.c^2)  I(year - 1974):len.c  
           1.8133187            -0.0862623             0.0590159             0.0005200             0.0004097  

> display(lake.lm5, 4)
lm(formula = log(pcb) ~ I(year - 1974) * len.c + I(len.c^2), 
    data = laketrout)
                     coef.est coef.se
(Intercept)           1.8133   0.0496
I(year - 1974)       -0.0863   0.0036
len.c                 0.0590   0.0043
I(len.c^2)            0.0005   0.0001
I(year - 1974):len.c  0.0004   0.0003
---
n = 631, k = 5
residual sd = 0.5452, R-Squared = 0.68
类型预测变量

很多情况下我们的预测变量并不是连续变量，而是离散的，比如分类变量。在我们的案例中，鱼类的大小虽然是连续的但是在拟合的过程中我们发现鱼的长度大小在60cm两侧趋势是不同的。于是我们构造一个类型预测变量size,作为因子变量，size被转化成取值为0或者1 的名义变量（Dummy）。

 (lake.lm6 <- lm(log(pcb) ~ I(year-1974)*factor(size) + 
+                  len.c * factor(size),
+                data=laketrout))

Call:
lm(formula = log(pcb) ~ I(year - 1974) * factor(size) + len.c * 
    factor(size), data = laketrout)

Coefficients:
                     (Intercept)                    I(year - 1974)                 factor(size)small                             len.c  
                       1.7394281                        -0.0846170                        -0.0646686                         0.0775689  
I(year - 1974):factor(size)small           factor(size)small:len.c  
                       0.0001259                        -0.0344983  

> display(lake.lm6, 4)
lm(formula = log(pcb) ~ I(year - 1974) * factor(size) + len.c * 
    factor(size), data = laketrout)
                                 coef.est coef.se
(Intercept)                       1.7394   0.0667
I(year - 1974)                   -0.0846   0.0044
factor(size)small                -0.0647   0.1197
len.c                             0.0776   0.0044
I(year - 1974):factor(size)small  0.0001   0.0074
factor(size)small:len.c          -0.0345   0.0063
---
n = 631, k = 6
residual sd = 0.5426, R-Squared = 0.68

简化模型

> lake.lm7 <- lm(log(pcb) ~ I(year-1974) + len.c * factor(size),
+               data=laketrout)
> display(lake.lm7, 4)
lm(formula = log(pcb) ~ I(year - 1974) + len.c * factor(size), 
    data = laketrout)
                        coef.est coef.se
(Intercept)              1.7389   0.0588
I(year - 1974)          -0.0846   0.0035
len.c                    0.0776   0.0044
factor(size)small       -0.0631   0.0779
len.c:factor(size)small -0.0345   0.0062
---
n = 631, k = 5
residual sd = 0.5422, R-Squared = 0.68

估计出的len.c的斜率和截距并未发生变化，要直接给出两类尺寸类别下的len.c的截距和斜率：

> lake.lm8 <- lm(log(pcb) ~ I(year-1974) + 
+                  len.c * factor(size)-1-len.c,
+               data=laketrout)
> display(lake.lm8, 4)
lm(formula = log(pcb) ~ I(year - 1974) + len.c * factor(size) - 
    1 - len.c, data = laketrout)
                        coef.est coef.se
I(year - 1974)          -0.0846   0.0035
factor(size)large        1.7389   0.0588
factor(size)small        1.6758   0.0795
len.c:factor(size)large  0.0776   0.0044
len.c:factor(size)small  0.0431   0.0045
---
n = 631, k = 5
residual sd = 0.5422, R-Squared = 0.83
共线性

为了量化胡泊的营养物质输入和湖内浮游植物生长的关系，线性回归模型是最常用的方法。由于磷和氮是藻类生长所需要的两种营养物质，他们也常常是藻类生长的限制因子。藻类的生长常用湖泊中叶绿素a的浓度来表示。

这关系通常表示为双对数线性回归模型：

\log Chla = \beta _{0}+\beta _{1}\log TP + \beta _{2}\log TN \varepsilon _{0}

载入数据看看变量之间的关系。

### The Finnish Lakes example ###

summer.All <- read.table(paste(dataDIR, "summerAll.csv", sep="/"), sep=",",header=T)

#names(summer.All)
# [1] "totp"  "chla"  "type"  "lake"  "year"  "totn"  "month" "depth" "surfa"
#[10] "color"

summer.All <- summer.All[log(summer.All$chla) > -20 ,]
#> names(summer.All)
# [1] "totp"  "chla"  "type"  "lake"  "year"  "totn"  "month" "depth" "surfa"
#[10] "color"

summer.All$y <- log(summer.All$chla)
summer.All$lxp<- scale(log(summer.All$totp), scale=F)
summer.All$lxn<- scale(log(summer.All$totn), scale=F)
summer.All$type.lake <- paste(summer.All$type, summer.All$lake)
summer.All$npr <- scale(log(summer.All$totn/summer.All$totp), scale=F)

### collinearity example ####
lake1 <- summer.All[summer.All$lake==19600 & summer.All$type==2,]
lake1$lxn <- scale(log(lake1$totn), scale=F)
lake1$lxp <- scale(log(lake1$totp), scale=F)
lake1$npr <- scale(lake1$npr, scale=F)
lake2 <- summer.All[summer.All$lake==1070,]
lake2$lxn <- scale(log(lake2$totn), scale=F)
lake2$lxp <- scale(log(lake2$totp), scale=F)
lake2$npr <- scale(lake2$npr, scale=F)
## lake 2 pairs
pairs(log(chla)~log(totn)+log(totp), data=lake2)

par( mgp=c(1.25,.125,0), las=1, tck=0.01)
pairs(y~lxn+lxp+npr, data=lake2, cex=0.5,
      labels=c("log Chla","log TN","log TP","log N:P"))

> Finn.lm1 <- lm(y ~ lxp + lxn, data=lake1)
> display(Finn.lm1)
lm(formula = y ~ lxp + lxn, data = lake1)
            coef.est coef.se
(Intercept)  2.35     0.02  
lxp          1.01     0.06  
lxn         -0.09     0.13  
---
n = 321, k = 3
residual sd = 0.37, R-Squared = 0.48
>

湖泊内TP和TN并不是完全独立的就要考虑他们的相互作用，换句话说要考虑两者的共线性问题。条件图是在一个变量固定的条件下考察另一个预测变量与响应变量之间的关系。

par(mar=c(3,3,3,.25), mgp=c(1.25,0.125,0), tck=0.02)
given.tn <- co.intervals(lake1$lxn, number=4, overlap=.1)
 coplot(y ~ lxp | lxn, data = lake1, given.v=given.tn, rows=1,
        xlab=c("log TP", "Given: Centered log TN"), ylab="log Chla",
        panel=panel.smooth, col="gray", col.smooth=1)

以TN为条件，对叶绿素和TP作图，在TN取值范围内，绘制叶绿素对数和TP对数之间的关系图。

给定TP条件：

par(mar=c(3,3,3,.25), mgp=c(1.25,0.125,0), tck=0.02)
given.tp <- co.intervals(lake1$lxp, number=4, overlap=.1)
coplot(y ~ lxn | lxp, data = lake1, given.v=given.tp, rows=1,
       xlab=c("log TN", "Given: Centered log TP"), ylab="log Chla",
       panel=panel.smooth, col="gray", col.smooth=1)

比较以上两个条件图可以看出，叶绿素a和TP的关系是相对稳定的，在不同TN条件下，有着相似的斜率，即暗示这是一个磷限制的湖泊。

> Finn.lm3 <- lm(y ~ lxp * lxn, data=lake1)
> display(Finn.lm3)
lm(formula = y ~ lxp * lxn, data = lake1)
            coef.est coef.se
(Intercept)  2.34     0.02  
lxp          1.01     0.06  
lxn         -0.09     0.13  
lxp:lxn      0.28     0.31  
---
n = 321, k = 4
residual sd = 0.37, R-Squared = 0.48

par(mfrow=c(1,2), mgp=c(1.25, 0.125,0), mar=c(3,3,1, 0.25), las=1, tck=0.02)
plot(y~lxp, data=lake1, axes=F, xlab="TP", ylab="Chla", cex=0.5)
axis(1, at=log(c(1, 2, 5, 10,20,50,100))-attr(lake1$lxp, which="scaled:center"),
     labels=c(1, 2, 5, 10,20,50,100))
axis(2, at=log(c(1,2,5,10,20,50)),
     labels=c(1,2,5,10,20,50))
box()
curve(coef.lm3[1]+coef.lm3[2]*x+coef.lm3[3]*quant.n[1] +
      coef.lm3[4]*x*quant.n[1], add=T, lty=1)
curve(coef.lm3[1]+coef.lm3[2]*x+coef.lm3[3]*quant.n[2] +
      coef.lm3[4]*x*quant.n[2], add=T, lty=2)
curve(coef.lm3[1]+coef.lm3[2]*x+coef.lm3[3]*quant.n[3] +
      coef.lm3[4]*x*quant.n[3], add=T, lty=3)
curve(coef.lm3[1]+coef.lm3[2]*x+coef.lm3[3]*quant.n[4] +
      coef.lm3[4]*x*quant.n[4], add=T, lty=4)
curve(coef.lm3[1]+coef.lm3[2]*x+coef.lm3[3]*quant.n[5] +
      coef.lm3[4]*x*quant.n[5], add=T, lty=5)
legend(x=log(30)-attr(lake1$lxp, which="scaled:center"),
       y=log(5),
       legend=c("2.5\\%","25\\%","50\\%","75\\%","97.5\\%"),
       cex=0.5, lty=1:5, bty="n")


plot(y~lxn, data=lake1, axes=F, xlab="TN", ylab="Chla", cex=0.5)
axis(1, at=log(c(500,1000))-attr(lake1$lxn, which="scaled:center"),
     labels=c(500,1000))
axis(2, at=log(c(1,2,5,10,20,50)),
     labels=c(1,2,5,10,20,50))
box()
curve(coef.lm3[1]+coef.lm3[3]*x+coef.lm3[2]*quant.p[1] +
      coef.lm3[4]*x*quant.p[1], add=T, lty=1)
curve(coef.lm3[1]+coef.lm3[3]*x+coef.lm3[2]*quant.p[2] +
      coef.lm3[4]*x*quant.p[2], add=T, lty=2)
curve(coef.lm3[1]+coef.lm3[3]*x+coef.lm3[2]*quant.p[3] +
      coef.lm3[4]*x*quant.p[3], add=T, lty=3)
curve(coef.lm3[1]+coef.lm3[3]*x+coef.lm3[2]*quant.p[4] +
      coef.lm3[4]*x*quant.p[4], add=T, lty=4)
curve(coef.lm3[1]+coef.lm3[3]*x+coef.lm3[2]*quant.p[5] +
      coef.lm3[4]*x*quant.p[5], add=T, lty=5)

参考：
一元线性回归的细节
机器学习之线性回归(linear regression)
恢复更新】写给大家看的机器学习书（第七篇）——线性回归（上）
一元线性回归分析
回归分析中的“回归”是什么意思？
机器学习Chapter-1（线性模型）
【统计学习3】线性回归：R方（R-squared）及调整R方（Adjusted R-Square）
加权最小二乘法与局部加权线性回归
Why You Need to Check Your Residual Plots for Regression Analysis: Or, To Err is Human, To Err Randomly is Statistically Divine

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,937评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,503评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,712评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,668评论 1赞 276
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,677评论 5赞 366
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,601评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,975评论 3赞 396
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,637评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,881评论 1赞 298
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,621评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,710评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,387评论 4赞 319
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,971评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,947评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,189评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 44,805评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,449评论 2赞 342

推荐阅读更多精彩内容