asreml分析阈值性状笔记

asreml 分析阈值性状笔记。

理论

asreml 使用广义线性模型来分析阈值性状（通过为二分类性状，通常呈现二项分布）。对应《ASReml-R Reference Manual Version 4》书中的 39 页。

广义线性模型

广义线性模型（Generalized Linear Model, GLM）可以被人为是传统一般线性模型的延伸。

广义线性模型包括以下几个部分：

1. 因变量：因变量不同取值间相互独立，服从指定数族概率分布。即，不局限于正态分布一种情形。

2. 线性部分：该部分和传统线性部分没有什么区别：

   

3. 连接函数：用于描述因变量的期望是怎样和线性预测值 相关联的。

广义线性模型在以下方面做了扩展：

1. 放松因变量的限制，因变量分布由正态分布扩展到二项分布、Poisson 分布、负二项分布等指数分布族
2. 通过连接函数，把因变量取值变换到自变量的线性预测的取值范围。

如果观测值之间存在相关性（非独立），那么我们就需要使用广义线性混合模型（Generalized Linear Mixed Model, GLMM），模型中可以对残差项进行相关性的定义。

下图是 asreml 支持的数据类型和连接函数，这里 表示累积标准正态分布函数。（其实我感觉这个表说的不清楚，而且还少了 multinomial 分布，其可用的连接函数与 binomial 相同）

这里的写法就是补充一个 family 语句，后面 asr_* 中后* 就是因变量的分布类型的名称，之后括号中就是连接函数。

这里下面提到对于二分类性状，两个连接函数的模型公式如下（我不清楚 Logit 模型中的 exp() 函数哪里来的，和上面的公式不一样）

实操

现在我不太想自己测试一下，先把代码放上来吧。

二分类性状

数据部分：表型中为 0 和 1。

代码部分：可以选 logistic模型和probit模型 ，二者只是加一个 family 语句而已。

但是需要注意的是，logistic模型算出的方差组分，计算遗传力的时候，残差要改为 本身的值乘以 （不清楚原因，本身残差的结果均为1，应该是设定好的残差是1），举例代码如下（probit模型计算遗传力的时候不用改）

（这里直接加的 pi^2/3 ，应该写成 V2*pi^2/3 比较好）

### 二分类logistic模型
str(dat)
m2 = asreml(y1_bio ~ G + Sex + Shidai,
            random = ~ vm(ID,ainv),
            family = asr_binomial(link="logit"),
            # residual = ~ idv(units),
            na.action = na.method(x = "include",y = "include"),
            data=dat)
m2 = update(m2)
summary(m2)$varcomp
vpredict(m2, h2 ~ V1/(V1+pi^2/3)) # logit的残差为pi^2/3

提取育种值结果就和正常分析一样，举例如下

# 育种值
hblup_y1 = as.data.frame(coef(mod_y1)$random)
head(hblup_y1)

hblup_y1 = hblup_y1 %>% mutate(tid = rownames(.)) %>% separate(tid,into = c("type","ID"),sep="_") %>% select(type,ID,effect)
head(hblup_y1)

fwrite(hblup_y1,"re_hblup_y1.csv")

多分类性状

表型：性状必须转为 factor (好像012或123这两种都可以，官方貌似建议012这种标记方式)

代码：asr_binomial 改为 asr_multinomial ，举例如下

mul_logit = asreml(mulb~1,
                   random=~ Sire+Family,
                   family=asr_multinomial(link="logit"), 
                   data=data)

连续性状和阈值性状多性状模型

表型：同上

代码：和普通多性状模型类似，不同的是增加了性状的分布函数，也就是还是 family 语句部分。

举例如下，这里第一个是二分类阈值性状，第二类为服从正态分布的连续性状，因此 family 语句 中按照相应顺序去写 。这里没有指定连接函数，因此就是默认函数，二分类性状就是 logistic模型，连续性状就是 identity 模型。

mod3 = asreml(cbind(score4,norm) ~ trait + trait:(Sex + Grp + year), 

              random = ~ us(trait):Sire, 

              family = list(asr_binomial(),asr_gaussian()),

              maxit = 100,

              data=dat)

但是视频中有个问题，视频中计算遗传相关和表型相关的时候是按照正常的方式算的，没有考虑 logistic模型单性状模型中的 的系数，不清楚原因（我感觉有可能是视频中搞错了，不过不清楚残差协方差是不是和残差方差一样校正，不过好像这个广义线性模型算出来的遗传力啥的本身就可能不太可靠）。

阈值性状和阈值性状多性状模型

同上，还是改一下 family 语句部分 ，举例如下

mod2 = asreml(cbind(score4,rot) ~ trait + trait:(Sex + Grp + year), 

              random = ~ us(trait):Sire, 

              family = list(asr_binomial(),asr_binomial()),

              maxit = 100,

              data=dat)

参考文献

1.https://www.bilibili.com/video/BV1oi4y1p7uX/

2.https://www.bilibili.com/opus/965825234531778579

3.《ASReml-R Reference Manual Version 4》