作者：林元震（ https://yzhlinscau2018.netlify.com）

说明：由于科学网对代码的支持不好，感兴趣的读者建议浏览我的netlify官方博客。

前言：前段时间，因疫苗问题闹得沸沸扬扬，这是非常严重的问题，关乎民生，国家应该彻查。因违法惩罚偏低，使得部分奸商为利铤而走险，不顾百姓健康。常言道，出来混的，迟早要还！我愿相信国家会在法制上想法制止这种损人利己的投机倒把。国人好斗，早已深入人心，诸如四大名著，哪部不言斗？啥时国人能跟数据好好斗斗？如果数据是原料，统计方法就是菜谱。理论上，给定一些原料，应该会得到一份合适的最佳菜谱。同理，对于给定一份数据，理论上也应有最合适的统计方法。不同统计方法，得到的结果可能会相差万里。因此，若要从事严谨的科研生产，应当关注相关行业的最新统计方法。懒散得随意分析生产数据，是一种罪过，也是对事业的不尊重，还是对人生的不尊重。今天演示的空间分析，即会看到最佳菜谱的重要性。

空间分析(spatial analysis)现已成为农林业遗传评估的常见方法之一，但其最早提出者竟是一位医生¹，这出乎我的意料。话说1840s末期，伦敦爆发了严重的霍乱疫情，死了不少人。主流理论认为霍乱通过“肮脏的空气”传染的，这错误理论导致霍乱未能得到有效控制。伦敦皇家外科医学院有位医生，名叫John Snow，虽然年纪轻轻但已事业有成，他不认同霍乱传染的主流理论，他认为水里携带的细菌才是主因。他用几年时间，走遍伦敦进行详细调研和病情记录，绘制了如下的点图：

最后，Snow将传染源锁定在一个公共抽水机，并说服政府将抽水机挪走。因为当时泰晤士河已被污染，然后通过抽水机传播霍乱。当年Snow锁定的抽水机，现已成了伦敦的地标之一：

Snow调查表明疾病的发病率跟地理位置是相关的，这就是空间分析方法的始祖。可惜，Snow英年早逝，四十多岁就累死了。

牛人已去，历史继续。Snow肯定没想到他的方法，被后人广泛运用到诸多行业，比如农林业上。现在介绍农林业上应用的空间分析原理及其示例。

1 空间分析的基本原理

遗传评估分析的通用混合线性模型如下：

$$y=xb+zb+e$$

对于空间分析²，将上式中的误差$e$进一步分解为$e=\xi +\eta$。$\xi$为空间相关的误差，$\eta$为空间不相关的随机误差。$\eta$反应了微环境差异、非加性遗传效应以及测量误差。

Authur等(1997)将$\xi$设为行、列的自回归AR1AR1，因此其方差协方差矩阵为：

$$Var\left(\xi \right)={\sigma }_{\xi }^2\left[{\Sigma }_c\right({\rho }_c\left){\Sigma }_r\right({\rho }_r\left)\right]$$

其中，为矩阵的Kronecker积(kronecker product)，${\Sigma }_c\left({\rho }_c\right)、{\Sigma }_r\left({\rho }_r\right)$为列、行的回归矩阵，列、行自回归参数分别为${\rho }_c、{\rho }_r$。

${\Sigma }_c\left({\rho }_c\right)$回归矩阵具体如下：

$$\begin{array}{ccccc}1& {\rho }_c& {\rho }_c^2& \dots & {\rho }_c^{n_{c}1}\\ {\rho }_c& 1& {\rho }_c& \dots & \\ {\rho }_c^2& {\rho }_c& 1& \dots & \\ & & & & \\ {\rho }_c^{n_{c}1}& & & \dots & 1\end{array}$$

${\Sigma }_r\left({\rho }_r\right)$回归矩阵具体如下：

$$\begin{array}{ccccc}1& {\rho }_r& {\rho }_r^2& \dots & {\rho }_r^{n_{r}1}\\ {\rho }_r& 1& {\rho }_r& \dots & \\ {\rho }_r^2& {\rho }_r& 1& \dots & \\ & & & & \\ {\rho }_r^{n_{r}1}& & & \dots & 1\end{array}$$

假定$\eta$随机误差是独立的，则误差方差协方差矩阵为：

$$Var\left(e\right)={\sigma }_{\xi }^2\left[{\Sigma }_c\right({\rho }_c\left){\Sigma }_r\right({\rho }_r\left)\right]+{\sigma }_{\eta }^{2}I$$

其中，${\sigma }_{\xi }^2$是空间相关的误差，${\sigma }_{\eta }^2$是独立的误差，$I$是单位矩阵。

数学公式不少，部分读者估计已经懵了。空间分析原理通俗点说，正如古人所言的‘近朱者赤，近墨者黑’，拿试验林为例，就是部分林地内树木都表现良好，而部分林地下却总体表现差，于是问题就来了，这些局部表现优秀或低劣的树木，到底是因为所处的微环境(比如光照、土壤养分等)还是其自身的基因型造成的？这里面事实上有两种情况：一种是基因型优秀的，但运气欠佳，落到了贫瘠的土壤里，成了‘怀才不遇型’；另一种是基因型低劣的，但运气爆棚，落到了肥沃的土壤里，成了‘滥竽充数型’。因此，空间分析就是要纠正这两种情况，通过空间自相关或者其它途径。有趣的是，基因型优秀或者低劣，是相对于试验环境而言的，即某个基因型在一种试验环境中是优秀的，但到了另一种试验环境中可能就成了低劣型。这属于基因型与环境互作的问题，本文不作这方面的分析。

现在来看看在R中如何使用breedR包进行空间分析的演示，虽然asreml包更强大，但由于它是商业软件，且费用不低，因此本文不打算介绍asreml。

2 空间模型的分析示例

首先载入所需程序包

library(breedR)
library(breedRPlus)
library(ggplot2)

(2.0) 示例数据集

示例采用agridat包的数据集stroup.nin，该数据含有小麦基因型gen、区组rep、产量yield、列号col和行号row等5个变量。

data(stroup.nin,package='agridat')
dat <- stroup.nin

数据集的变量及其特征：

names(dat)
## [1] "gen"   "rep"   "yield" "col"   "row"str(dat)
## 'data.frame':    242 obs. of  5 variables:
##  $ gen  : Factor w/ 56 levels "Arapahoe","Brule",..: 12 12 12 12 12 12 12 12 12 12 ...
##  $ rep  : Factor w/ 4 levels "R1","R2","R3",..: NA NA NA NA NA NA NA NA NA NA ...
##  $ yield: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ col  : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ row  : int  1 1 1 1 1 1 1 1 1 1 ...

目标性状yield的直方图(Fig.1)和空间趋势(Fig.2)如下：

ggplot(dat, aes(yield)) + geom_histogram(col='blue',fill='white')

Figure 1: 目标性状yield的直方图

ggplot(dat, aes(row, col, fill = yield)) +
  geom_tile(show.legend = TRUE) +
  coord_fixed() + coord_flip()+
  scale_fill_gradientn(colours = topo.colors(50))

Figure 2: 目标性状yield的空间趋势

(2.1)分析模型

(2.1.1) 模型1–RCB模型

首先采用原本的试验设计–随机完全区组设计(RCB)，进行数据分析。由于数据要针对基因型的品种表现作评估与比较，因此将其作为固定效应。具体代码如下：

RCB.bdR <- remlf90(fixed = yield ~ 1+gen,
                   random = ~ rep,
                   data = dat)

RCB模型的方差分量以及Loglook与AIC值如下：

breedRPlus::var(RCB.bdR)
##              gamma component std.error  z.ratio constraint
## rep      0.1993244    9.8829    8.7928 1.123976   Positive
## Residual 1.0000000   49.5820    5.4588 9.082949   Positive

summary(RCB.bdR)$model.fit[c(1,3)] 
##       AIC    logLik
##  1253.938 -624.9691

此外，breedR允许将区组放入空间分析方法中，但结果与RCB的一样。分析代码如下：

RCB1.bdR <- remlf90(fixed = yield ~ 1+gen,
                    spatial = list(model = 'blocks',
                                   coord = dat[, c('row','col')],
                                   id = "rep"),
                    data = dat)

RCB1模型的方差分量以及Loglik与AIC值如下：

##              gamma component std.error  z.ratio constraint
## spatial  0.1993244    9.8829    8.7928 1.123976   Positive
## Residual 1.0000000   49.5820    5.4588 9.082949   Positive


##       AIC    logLik
##  1253.938 -624.9691

(2.1.2) 模型2–SUM空间模型

现在进行空间模型分析，由于breedR包不能设定σ2ηση2=0，因此它无法进行仅含有空间相关误差的模型，这与asreml相比，多少有点遗憾，但对于初学者，这却是优点！下面展示行列均相关的空间分析模型，具体代码如下：

SUM.bdR <- remlf90(fixed = yield ~ 1+gen,
                    spatial = list(model = 'AR',
                                   coord = dat[,c('row','col')],
                                   rho =c(.66,.44)),
                    data = dat)

SUM模型的方差分量以及Loglik与AIC值如下：

##             gamma component std.error  z.ratio constraint
## spatial  6.459346   35.7570    6.8376 5.229466   Positive
## Residual 1.000000    5.5357    2.5391 2.180182   Positive


##       AIC    logLik
##  1147.881 -571.9406

SUM模型的空间相关误差趋势图如下：

plot(SUM.bdR, type = 'spatial')+ coord_flip()+
  scale_fill_gradientn(colours = topo.colors(50))

Figure 3: SUM模型的空间相关误差趋势图

(2.1.3) 模型3–ColAR空间模型

现在设定只有列存在空间相关，即行之间不存在自相关。代码如下：

ColAR.bdR <- remlf90(fixed = yield ~ 1+gen,
                     spatial = list(model = 'AR',
                                   coord = dat[,c('row','col')],
                                   rho = c(0,.44)),
                     data = dat)

ColAR模型的方差分量以及Loglik与AIC值如下：

##             gamma  component std.error   z.ratio constraint
## spatial  26138.25 43.7110000 4.7691000 9.1654610   Positive
## Residual     1.00  0.0016723 0.0082369 0.2030254   Positive


##       AIC    logLik
##  1198.539 -597.2693

ColAR模型的空间相关误差趋势图如下：

#plot(ColAR.bdR, type = 'spatial')

(2.1.4) 模型4–RowAR空间模型

现在设定只有行存在空间相关，即列之间不存在自相关。代码如下：

RowAR.bdR <- remlf90(fixed = yield ~ 1+gen,
                     spatial = list(model = 'AR',
                                   coord = dat[,c('row','col')],
                                   rho = c(0.66,0)),
                     data = dat)

RowAR模型的方差分量以及Loglik与AIC值如下：

 
 ##            gamma  component std.error    z.ratio constraint
 ## spatial  20272.6 57.0390000   6.22540 9.16230282   Positive
 ## Residual     1.0  0.0028136   0.06103 0.04610192   Positive
 

 ##       AIC    logLik
 ##  1198.989 -597.4946

RowAR模型的空间相关误差趋势图如下：

#plot(RowAR.bdR, type = 'spatial')

(2.2) 四种模型的空间相关误差比较

注意：RCB1模型是把区组放到空间分析方法里，因此所计算的空间相关误差其实就是区组方差。

compare.plots(
  list(RCB = plot(RCB1.bdR, type = 'fullspatial'),
       SUM = plot(SUM.bdR, type = 'fullspatial'),
       ColAR = plot(ColAR.bdR, type = 'fullspatial'),
       RowAR = plot(RowAR.bdR, type = 'fullspatial')
   )) + coord_flip() +
  scale_fill_gradientn(colours = topo.colors(50))

Figure 4: 四种模型的空间相关误差趋势图

从图4可以明显看出空间模型与RCB模型的图形较大，虽然3种空间模型的变化比较小，但认真查看，还是能看到微小的差异。此外，对于空间分析，还有variogram(半残差)图，但因为breedR包提供的variogram图实在太丑了，这里不做演示。

(2.3) 四种模型的拟合统计量比较

sres<-function(object){
  df<-summary(object)$model.fit[c(1,3)]
  df[,c('Spatial','Residual')]=var(object)[,2]  return(df)
}

M1<-sres(RCB1.bdR)
M2<-sres(SUM.bdR)
M3<-sres(ColAR.bdR)
M4<-sres(RowAR.bdR)

res<-rbind (M1,M2,M3,M4)
row.names(res)<-c('RCB','SUM','ColAR','RowAR')

res

从统计量结果可知，对于含有空间相关的模型，不论是行列双向自相关还是行或列单向自相关，模型拟合的效果都要优于RCB模型，因为它们的AIC值都小于RCB模型。其实，这不难理解，因为传统的RCB模型，区组往往比较大，而在空间分析模型里，通过行列自相关模型，其实是把区组细化到小区(每个行列组合)，因此模型所估计的误差方差更为精确(例如RCB模型和SUM模型，随机误差方差已从49.582下降到5.536!)，当然AIC值就会变小。根据AIC准则，SUM模型是本例中的最优模型，即行列双向自相关。

(2.4) 四种模型的gen效应值比较

mvres<-function(object){
  df<-data.frame(mv=fixef(object)$gen,mv.se=attr(fixef(object)$gen,'se'))
  df$ID<-row.names(df)
  df<-dplyr::arrange(df,desc(mv))
  df$mv.rank<-1:nrow(df)  return(df)
}

fM1<-mvres(RCB.bdR)
names(fM1)<-c('m1','m1.se','ID','m1.rank')

fM2<-mvres(SUM.bdR)
names(fM2)<-c('m2','m2.se','ID','m2.rank')

fM3<-mvres(ColAR.bdR)
names(fM3)<-c('m3','m3.se','ID','m3.rank')

fM4<-mvres(RowAR.bdR)
names(fM4)<-c('m4','m4.se','ID','m4.rank')

fM.res11<-merge(fM1,fM2,by='ID')
fM.res12<-merge(fM3,fM4,by='ID')
fM.res<-merge(fM.res11,fM.res12,by='ID')

fM.res1<-dplyr::arrange(fM.res,desc(m1))
fM.res2<-dplyr::arrange(fM.res,desc(m2))

通过上述代码获取各模型的固定gen效应值，然后按最优模型M2的效应值降序排序，并输出结果的前8名。

head(fM.res1,8)

head(fM.res2,8)

从输出结果看，OMG，结果相差甚远！尤其是RCB模型与最优模型SUM之间，基本完全不同，例如SUM中最好的基因型Buckskin在模型RCB仅排名第28名，而且产量被RCB模型严重低估了！本例中，基因型Buckskin在RCB模型中就是‘怀才不遇型’；基因型NE86503在RCB模型中就是‘滥竽充数型’，合理的空间分析模型给予了纠正。即便是3种空间模型，它们之间的差异也较大。此外，最优模型的标准误se也最小！这结果验证了空间分析的优越性，也验证了统计分析方法的重要性！

空间分析结果对人生的提示：只要是金子，迟早会发光，即便一段时间内怀才不遇；如果是狐假虎威，迟早会败露，即便一段时间内风光无限。

3 结语

3.1 从本例的分析结果来看，选择最好的软件非常有必要！此外，关注新的统计方法也非常关键！本例结果显示SUM是最优模型，那么它有没再改进的可能？答案是肯定的，感兴趣的读者可以阅读本人的一篇文章³。

3.2 空间分析的使用条件。诚如已有菜谱，但能否每人都根据菜谱做出可口的上等好菜，这需要磨练和时间。同理，空间分析模型虽好，也得要知道它的使用条件和判断标准。首先，数据集里必须有行列号；其次，对于简单数据集，模型不可过于复杂，否则可能会造成过拟合(overfit)；再者，合理的空间模型可能并唯一，但不合理的空间模型可能会导致结果高估或低估。

后话：虽然空间分析已成为农林业国际的常用方法之一，但其在国内林业遗传评估上的应用仍然有限！我最近搜了搜关于数量遗传学的著作或教材，国内外也是更新极慢，经典当属《Intrudction to quantitative genetics(4th)》(1996年)和《Genetics and Analysis of Quantitative Traits》(1998年)，历经20年了！第一本经典的作者Douglaus Falconer，可惜已于2004年逝世。第二本的作者今年出版第二部《Evolution and Selection of Quantitative Traits》，9月份出版！其实，关于统计分析史的那些大师们，有些故事还是非常耐人寻味的！多少有点遗憾，迷人的统计故事都是国外的大师造就的。

参考资料

1. https://www.zhihu.com/question/28409028

2. 林元震主编.《R与ASReml-R统计学》.中国林业出版社.2017

3. http://blog.sciencenet.cn/blog-1114360-1071888.html