在"Infectivity, susceptibility, and risk factors associated with SARS-CoV-2 transmission under intensive contact tracing in Hunan, China"[1]一文中，作者使用了广义线性混合模型来研究SARS-Cov2易感性和传染性的潜在驱动因素。因此，针对该模型进行深入学习。
文章目录
线性回归模型
线性模型(linear model)试图学得一个通过属性的线性组合来进行预测的函数，既可以用于分类,也可以用于回归。在机器学习的术语中，当预测值为连续值时，称为“回归问题”，离散值时为“分类问题”。
对于线性回归模型，从统计学角度来看，因变量Y与设计矩阵X间存在统计关系：𝑌=𝑋𝛽+𝜀，其中Y为因变量的测量值向量，X为固定效应自变量的设计矩阵， β是与X对应的固定效应参数向量，𝜀为剩余误差向量。Xβ为在X条件下的Y的平均值向量，即𝐸[𝑌|𝑋]=𝑋𝛽.𝜀假定为独立、等方差及均值为0的正态分布，即𝜀∼𝑁(0,𝜎^2 )。
因此，使用一般线性模型时，是需要满足以下3点假设的：
1.正态性，因变量y符合正态分布。即 Y i N ( μ i , σ 2 ) , E [ Y i ] = μ i = X i β Y_i~N(μ_i,σ^2 ),E[Y_i ]=μ_i=X_i β Yi N(μi,σ2),E[Yi]=μi=Xiβ
2.独立性，不同类别y的观察值之间相互独立，相关系数为零
3.方差齐性，不同类别y的方差相等
广义线性模型
当Y不再服从正态分布，而是其他已知分布时，需要建立Y和𝜂=𝑋𝛽之间的联系，从而转化为线性模型的形式。在广义线性模型的理论框架中，则假设目标变量Y服从指数分布族。
指数分布族
指数分布族的概率密度函数为：
其中，𝜃𝑖是指数族的自然参数,𝜑𝑖称为尺度参数或讨厌参数；b(𝜃𝑖)和c(y𝑖,𝜑𝑖)是依据不同指数族而确定的函数。指数分布族的均值
E
(
Y
i
)
=
b
′
(
θ
i
)
;
方差
V
a
r
(
Y
i
)
=
E
[
Y
i
−
b
′
(
θ
i
)
]
2
=
φ
i
b
′
′
(
θ
i
)
E(Y_i )=b'(θ_i );方差Var(Y_i )=E[Y_i-b'(θ_i )]^2=φ_i b''(θ_i )
E(Yi)=b′(θi);方差Var(Yi)=E[Yi−b′(θi)]2=φib′′(θi)
下面以佰努利分布为例，证明其是否属于指数分布族：
已知佰努利分布B(1,P)的概率密度函数为：
f
(
y
∣
p
)
=
p
y
(
1
−
p
)
(
1
−
y
)
,
y
∈
0
,
1
f(y|p)=p^y (1-p)^{(1-y)},y∈{0,1}
f(y∣p)=py(1−p)(1−y),y∈0,1
我们将其转换为指数分布族的形式：
f
(
y
∣
θ
)
=
e
x
p
(
y
θ
−
l
o
g
(
1
+
e
θ
)
)
f(y|θ)=exp(yθ-log(1+e^θ ))
f(y∣θ)=exp(yθ−log(1+eθ))
其中，
θ
=
l
o
g
p
1
−
p
,
φ
=
1
,
b
(
θ
)
=
l
o
g
(
1
+
e
θ
)
,
c
(
y
,
φ
)
=
0
θ=log \frac{p}{1-p},φ=1,b(θ)=log(1+e^θ ),c(y,φ)=0
θ=log1−pp,φ=1,b(θ)=log(1+eθ),c(y,φ)=0
因此佰努利分布属于指数分布族。
正态分布同样属于指数分布族，因此线性回归可以看作是广义线性模型的特例。
广义线性模型结构
随机部分：随机部分的Y由N个观测值𝑌={
𝑌
1
𝑌_1
Y1,…,
𝑌
𝑁
𝑌_𝑁
YN }组成，每个
𝑌
𝑖
𝑌_𝑖
Yi相互独立且由相同的一维未知参数𝜃决定的指数分布族。结合指数分布族的均值公式可得：
E
[
Y
i
]
=
μ
i
=
∂
b
i
(
θ
)
∂
θ
E[Y_i ]=μ_i=\frac{∂b_i (θ)}{∂θ}
E[Yi]=μi=∂θ∂bi(θ)
系统部分：系统部分是一组解释变量构成的线性模型：
θ
i
=
X
i
β
θ_i=X_i β
θi=Xiβ。其中𝜃𝑖为自然参数，即指数族分布中的自然参数𝜃和X呈线性关系；
𝑋
𝑖
𝑋_𝑖
Xi 为p维矩阵(𝑖=1,…,𝑁),是对解释变量的一组观测值；𝛽是一组p维未知常向量。
连接函数：若
E
[
Y
i
]
=
μ
i
E[Y_i ]=μ_i
E[Yi]=μi,则连接函数可写为
g
(
μ
i
)
=
θ
i
=
X
i
β
,
i
=
1
,
…
,
N
g(μ_i )=θ_i=X_i β,i=1,…,N
g(μi)=θi=Xiβ,i=1,…,N
逻辑回归的广义线性模型推导
logistic回归是一种广义的线性回归分析模型，其因变量是二分类的。下面针对逻辑回归进行广义线性模型的推导：
因变量服从指数分布族？
逻辑回归是在对一个二分类问题进行建模，随机变量y取值为 0 或 1，即y服从伯努利分布。满足“在给定X的条件下，假设Y服从以η为自然参数的指数族分布，即𝑦|𝑥∼𝐸𝑥𝑝𝑜𝑛𝑒𝑛𝑡𝑖𝑎𝑙𝐹𝑎𝑚𝑖𝑙𝑦(𝜂)”
模型能否预测Y的期望值？
①由于y服从伯努利分布，则其均值为𝐸[𝑌]=0⋅𝑝(𝑦=0|𝑥)+1⋅𝑝(𝑦=1|𝑥)=𝑝。
②在前面证明伯努利分布属于指数分布族时，得出
𝑏
(
θ
)
=
𝑙
𝑜
𝑔
(
1
+
𝑒
θ
)
𝑏(θ)=𝑙𝑜𝑔(1+𝑒^θ )
b(θ)=log(1+eθ),𝜃=
𝑙
𝑜
𝑔
p
1
−
p
𝑙𝑜𝑔\frac{p}{1-p}
log1−pp，则𝐸[𝑌]=𝑏′ (𝜃)=
𝑒
θ
1
+
𝑒
θ
\frac{𝑒^θ}{1+𝑒^θ}
1+eθeθ=𝑝。
③满足“在给定X的条件下，我们的目标是得到一个模型h(x)能预测出Y的的期望值，即ℎ(𝑥)=𝐸[𝑦|𝑥]=𝜇𝑖”
自然参数与X呈线性关系？
①指数分布的自然参数𝜃与X呈线性关系，即
θ
=
𝑤
𝑇
𝑥
θ=𝑤^𝑇 𝑥
θ=wTx。
②在逻辑回归中
h
(
𝑥
)
=
𝐸
[
𝑦
∣
𝑥
]
=
𝑝
=
𝑒
θ
1
+
𝑒
θ
ℎ(𝑥)=𝐸[𝑦|𝑥]=𝑝=\frac{𝑒^θ}{1+𝑒^θ}
h(x)=E[y∣x]=p=1+eθeθ,所以存在连接函数g，使得
𝑔
(
𝑒
θ
1
+
𝑒
θ
)
=
θ
=
𝑤
𝑇
𝑥
𝑔(\frac{𝑒^θ}{1+𝑒^θ})=θ=𝑤^𝑇 𝑥
g(1+eθeθ)=θ=wTx.
③连接函数为
𝑔
(
𝑡
)
=
𝑙
𝑛
(
𝑡
1
−
𝑡
)
𝑔(𝑡)=𝑙𝑛(\frac{𝑡}{1−𝑡})
g(t)=ln(1−tt),即“sigmoid”函数，它将t值转化为一个接近0或1的y值。
线性混合模型
线性混合模型的提出
在上述线性模型中，我们知道因变量y是满足正态性、独立性、方差齐性的。当Y具有群体特性时，如在抽样调查中，被调查者会来自不同的城市、不同的学校，这就形成一个层次结构，同一城市或同一学校的学生各方面的特征应当更加相似。也就是基本的观察单位聚集在更高层次的不同单位中，如同一城市的学生数据具有相关性，不能满足观察值间的独立性。
如果我们对不同的群体分别建立各自的回归模型，当群体数较少，群体内样本容量较大，传统的分析方法可能是有效的。但是如果我们把这些群体看成是从总体中抽样来的一个样本(例如多阶段抽样和重复测度数据)，并想分析不同群体之间的总体差异，那么简单地使用传统的统计方法是不够的。
因此回归模型的提出既保留了传统线性模型中的正态性假定条件，又对独立性和方差齐性不作要求，从而扩大了传统线性模型的适用范围
线性混合模型结构
线性模型结构为𝑌=𝑋𝛽+𝜀，而混合线性模型将一般线性模型扩展为： Y = X β + Z Γ + ε Y=Xβ+ZΓ+ε Y=Xβ+ZΓ+ε
- 其中Z为随机效应变量构造的设计矩阵，其构成方式与X相同。
- Γ为随机效应参数向量，服从均值为0，方差为G的正态分布，即Γ∼𝑁(0，𝐺)。
- 随机误差𝜀服从均值为0，方差为R的正态分布，即𝜀~𝑁(0,𝑅)。放宽了对𝜀的限制条件，其元素不必为独立同分布，即没有𝐶𝑜𝑣(𝜀𝑖,𝜀𝑗 )=0的假定。
- G和R无相关关系，即𝐶𝑜𝑣(𝐺,𝑅)=0,此时Y的方差为𝑉𝑎𝑟(𝑌)=𝑍𝐺𝑍+𝑅,期望值为𝐸(𝑌)=𝑋𝛽。当Z=0，𝑅=𝜎^2 时，混合线形模型转变为一般线形模型。
相比简单线性模型，多出了Z这一项，这一项称之为随机效应。在一般线性模型中，其自变量全部为固定效应自变量，而线性混合模型中，除了固定效应自变量外，还包含了随机效应自变量。
线性混合模型关键之处在于判定自变量的类别，如果抽样数据集中的自变量可以包含该自变量的所有情况，则作为固定效应，如果只能代表总体的一部分，则作为随机效应。
广义线性混合模型
广义线性混合模型GLMM，可以看做是线性混合模型LMM的扩展形式，使得可以同时包含固定效应和随机效应；也可以看作是GLM的扩展形式，使得因变量不再要求满足正态分布。
在"Infectivity, susceptibility, and risk factors associated with SARS-CoV-2 transmission under intensive contact tracing in Hunan, China"[1]一文中，作者使用了广义线性混合模型来研究SARS-Cov2易感性和传染性的潜在驱动因素。模型定义如下：
- 𝑔(∙)为logit连接函数，即 𝑔 ( 𝑡 ) = 𝑙 𝑛 ( 𝑡 1 − 𝑡 ) 𝑔(𝑡)=𝑙𝑛(\frac{𝑡}{1−𝑡}) g(t)=ln(1−tt)
- 𝛼为随机误差
- 固定效应分别为感染者年龄组、接触者年龄组、接触类型、传播代数、感染者的密接数、感染者性别、接触者性别、感染者有无症状、观察时期。 β 𝑖 ( 𝑖 = 1 , … , 9 ) β_𝑖 (𝑖=1,…,9) βi(i=1,…,9)为固定效应的参数向量。
- 𝑢 0 𝑢_0 u0和 𝑢 1 𝑢_1 u1 是感染者个体和集群的随机效应。 𝑢 𝑖 = 𝐸 [ 𝑌 ∣ 𝑢 0 , 𝑢 1 ] 𝑢_𝑖=𝐸[𝑌|𝑢_0,𝑢_1 ] ui=E[Y∣u0,u1]是给定随机效应值的𝑌𝑖 的平均值(因为样本中感染者个体id与所在集群id不能覆盖总体所有取值，所以作为随机效应)
模型实现-glmer函数
广义线性混合模型的实现利用了R语言中lme4包中的glmer函数（也可利用MASS包建立回归模型）。函数基本表达式如下：
glmer(formula = DV ~ Fixed_Factor + (Random_intercept + Random_Slope | Random_Factor), data = , family = ，control = glmerControl())
- Formula参数中设定模型的固定效应和随机效应：
·DV为因变量
·Fixed_Factor固定因子，即考察的自变量
·Random_intercept为随机截距，即认为不同群体的因变量的分布不同
·Random_Slope为随机斜率，即认为不同群体受固定因子的影响是不同的
·Random_Factor为随机因子 - data参数用来选择数据集
- family参数用来设定模型链接函数的分布族
- control参数用来设定模型优化器
数据集描述
①数据集中因变量为infectee，表示个体的感染状态，有yes和no两种取值状态。
②自变量中agegroup_index为患者所在年龄组、agegroup为个体所在年龄组、contact_typr1为接触类型、generation为传播代数、gender_index为患者性别、gender为个体性别、clinic_index为患者是否出现症状。它们都是分类变量。
③自变量中no.persons为患者的密接数、logage_index为患者年龄的对数、logage为个体年龄对数。它们都是连续变量。
④id_case为患者id、cluster_id为集群id，为随机效应。
回归模型中的分类变量
分类变量是取值是有限的类别值，如性别：男、女。分类变量是不能直接用到回归模型中的，即使用 1 表示男，用 0 表示女，这个 1 和 0 仍然只能是起类别区分的作用，如果不加处理让它们当数值 1 和 0 使用了，那么整个模型的逻辑和结果都是不正确的。R中分类变量只要是因子型或字符型，当加入回归模型时，不需要做任何额外操作将自动处理成虚拟变量用进模型。
以contact_type为例，该分类变量包含四个类别：”Household contacts”,” Social contacts”,” Relative contacts”,” Other contacts”。虚拟变量是一种二值变量（0-1），只表示是否。二分类或多分类变量，可以这样转化为虚拟变量，即contact_type变量转化为多个二值变量：
contact_type是否为Household contacts，是否为Social contacts，是否为Relative contacts，是否为Other contacts
若contact_type=Household contacts，要用上述4个二值变量表示的话，就是分别为1, 0, 0，0
由于4个虚拟变量存在冗余，需要将一个类别作为参照组。比如去掉contact_type是否为Household contacts，就相当于“Household contacts”组是参照组，另外3组与参照组做比较。模型将变为y=β1Social+β2Relative+β3*Other。当接触类型为social时，即Social=1，与Household相比，Social的传播风险增加了β1.
代码实现
引入数据和包
library(lme4)
library(optimx)
library(tidyverse)
library(broom.mixed)
setwd(".../1. input/")
data<-read.csv("data_model.csv")
str(data)
设置分类因子
在R语言中需要将字符型变量转化为分类因子并设置参照才能放入模型
data$Infectee = factor(data$Infectee)
data$agegroup_index = relevel(factor(data$agegroup_index),ref="15-64y")
data$agegroup = relevel(factor(data$agegroup),ref="15-64y")
data$contact_type1 = relevel(factor(data$contact_type1),ref="Household contacts")
data$gender = relevel(factor(data$gender),ref="Male")
data$clinic_index2 = relevel(factor(data$clinic_index2),ref="Symptomatic infectors")
建立回归模型
model= glmer(Infectee ~ agegroup_index + agegroup + no.persons + gender + contact_type1 +
clinic_index2 + (1|id_case) + (1|cluster_id)
,family = binomial(link = "logit"), data = data,
control=glmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
result<-tidy(model,conf.int=TRUE,exponentiate=TRUE,effects="fixed")
result1<-summary(model)
- (1|id_case)和(1|cluster_id)表示感染者个体和所在群体的随机效应
- family表示连接函数选用二分类logit函数
- 将optimx作为优化器
- tidy函数可得出回归系数对应的OR值
- summary函数可将模型结果输出。
回归结果
- Estimate为相应的自变量改变一个单位，应变量的改变量，即自变量前面的系数（回归系数）。
- Std. Error为标准误差。在估计解释变量对于被解释变量的影响时，实际上是估计的解释变量对于被解释变量的数学期望的影响，E(y)=a+b*x，是一个均值估计。在估计的时候都是用样本估计的，抽取一个样本就可以得到一个估计系数，所以估计系数本身就是随机变量。而通过抽样获得的随机变量的标准差就叫做标准误差。
- Z value是回归系数除以其标准误差。它有时也称为 z 统计量。如果 z 值的幅度太大（即太正或太负），对应的自变量很重要。
- Pr(>|z|)是p值用概率记法的表示。p值估计了系数不显著的可能性，越小越好。如果p值很大，说明不显著的可能性很高。（由于数据集不是论文中真实数据，而是实例。带入回归模型中拟合效果不是很好）
衡量模型拟合效果
文章建立了两个回归模型，因此需要衡量两个模型的好坏：
modelA：患者和密接者年龄为分类变量，即所在年龄组
modelB：患者和密接者年龄为连续变量，即真实年龄的对数。
两者比较如下：
- AIC信息准则是衡量统计模型拟合优良性的一种标准。其表达式为AIC=-2ln(L) + 2k。其中k为参数数量；L是似然函数。AIC越小，模型拟合效果越好。
- BIC是贝叶斯信息准则。它相比于AIC考虑了样本的数量，可有效防止模型精度过高造成的模型复杂度过高。表达式为BIC=ln(n)k-2ln(L)。
- Loglikelihood为对数似然函数值，用来反映模型的拟合程度，越小越好。即上述AIC、BIC表达式中的ln(L)。
OR值
由于文章结果分析中利用变量间OR值进行对比，因此对OR值进行简要介绍。
odds: 称为几率、比值、比数，是指某事件发生的可能性(概率)与不发生的可能性（概率）之比。用p表示事件发生的概率，则：odds = p/(1-p)。
OR:比值比，为实验组的事件发生几率(odds1)/对照组的事件发生几率(odds2)
- OR值大于1，表示该因素对事件发生起促进作用；
- OR值小于1，表示该因素对事件发生起抑制作用；
- OR值等于1，表示该因素对事件发生无影响。
以家庭接触对患病影响为例, 则 𝑂 𝑅 = 𝑎 / 𝑐 𝑏 / 𝑑 = 𝑎 𝑑 / 𝑏 𝑐 𝑂𝑅=\frac{𝑎/𝑐}{𝑏/𝑑}=𝑎𝑑/𝑏𝑐 OR=b/da/c=ad/bc:
回归系数指数化就得到对应的自变量改变一个单位的优势比OR,即 𝑂 𝑅 = 𝑒 β 𝑂𝑅=𝑒^β OR=eβ。β为对应自变量的回归系数。
文章结果分析
在文章回归结果表格中，以年龄因素为例，对于年龄为0-14岁的儿童，其接触者患病概率为a/193；不患病的概率为(193-a)/193，所以其接触者患病比率odd1=a/(193-a)。相应的，对于年龄在15-64岁的成人，其接触者患病概率为b/6833；不患病的概率为(6833-b)/6833，所以其接触者患病比率odd2=b/(6833-b)。0-14岁年龄区间相对于15-64岁年龄区间的OR值为odd1/odd2=a*(6833-b)/(b*(193-a))=0.17。可以认为，0-14岁儿童相比15-64岁成人，接触者患病概率低83%，即成人感染性更强。
从回归结果来看，0-14岁年龄组变量的回归系数为ln（0.17）=-1.77。
其中a为0-14岁的儿童中接触者患病人数，b为15-64岁的成人中接触者患病人数。
上图则为回归结果的可视化。
a.15 岁以下的接触者和 65 岁以上的接触者的相对易感性（参考组 [红线] 是 15-64 岁的接触者）
b.接触过不同类型感染者的接触者感染 SARS-CoV-2 的相对风险（参考组 [红线] 是在家庭中接触过感染者的接触者）
c.接触过不同代 SARS-CoV-2 传播的感染者的接触者感染 SARS-CoV-2 的相对风险（参考组 [红线] 是接触过指示病例的接触者）
d.在特定环境中，特定年龄的接触者感染 SARS-CoV-2 的概率。
(a 到 c 中的箱线图显示了与参考组相比，SARS-CoV-2 感染的相对风险的点估计值和 95% 置信区间。d 中的线 和阴影区域分别代表 SARS-CoV-2 感染概率的点估计值和 95% 置信区间。)
总结
回归模型在探索SARS-Cov2传播潜在因素中广泛应用。Ge[2]等人使用具有稳健 SE 方差的修正混合效应泊松回归来所有分析,得出与接触无症状指示患者相比，接触有症状患者感染风险更高等结论。Hall[3]等人通过逻辑回归分析估计了不同组在指示病例发生后 2 周内成为继发病例的几率，得出人数较多的家庭更有可能出现继发病例等结论。Li[4]等人使用链二项式传播模型来估计二次攻击率。该模型还用于评估传染性和感染易感性的决定因素,得出60岁或以上的人感染 SARS-CoV-2 的风险更高等结论。
建立回归模型的步骤可分为：
参考文献
[1] Hu S, Wang W, Wang Y, et al. Infectivity, susceptibility, and risk factors associated with SARS-CoV-2 transmission under intensive contact tracing in Hunan, China[J]. Nature communications, 2021, 12(1): 1-11.
[2] Ge Y, Martinez L, Sun S, et al. COVID-19 transmission dynamics among close contacts of index patients with COVID-19: a population-based cohort study in Zhejiang province, China[J]. JAMA Internal Medicine, 2021, 181(10): 1343-1350.
[3] Hall J A, Harris R J, Zaidi A, et al. HOSTED—England’s Household Transmission Evaluation Dataset: preliminary findings from a novel passive surveillance system of COVID-19[J]. International journal of epidemiology, 2021, 50(3): 743-752.
[4] Li F, Li Y Y, Liu M J, et al. Household transmission of SARS-CoV-2 and risk factors for susceptibility and infectivity in Wuhan: a retrospective observational study[J]. The Lancet Infectious Diseases, 2021, 21(5): 617-628.
6895
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
