在"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.