基于纵向数据半参数混合效应模型的网络营销影响因素研究

Research on the Influencing Factors of Online Marketing Based on Longitudinal Data Semi Parametric Mixed Effects Model

摘要: 本文建立了一个纵向数据的半参数混合效应模型，采用B样条方法建立纵向指标与时间的关系，为模型中的所有未知参数和随机效应指定先验分布，在共轭先验条件下，对该模型进行贝叶斯分析，通过Gibbs抽样方法获取未知参数、随机效应和非参数函数的贝叶斯估计。最后将该模型应用于饮食类非遗产品网络营销的影响因素进行分析，实证结果显示本文建立的指标对淘宝店铺的网络销售量产生了不同程度的影响，证明了模型的有效性。

Abstract: This article establishes a semi parametric mixed effects model for longitudinal data, using the B-spline method to establish the relationship between longitudinal indicators and time. Prior distributions are specified for all unknown parameters and random effects in the model, and Bayesian analysis of the model is considered under conjugate priors. Bayesian estimates of unknown parameters, random effects, and non parametric functions are extracted based on Gibbs sampling, and the model is applied to the study of influencing factors in online marketing. The empirical results show that the indicators established in this article have varying degrees of impact on the online sales volume of Taobao stores, proving the effectiveness of the model.

1. 引言

网络营销是指利用互联网及其相关技术，通过各种在线渠道和工具进行产品或服务的推广、销售和品牌建设的过程。它包括但不限于搜索引擎优化(SEO)、搜索引擎营销(SEM)、社交媒体营销、电子邮件营销、内容营销和在线广告等策略。网络营销的目标是吸引潜在客户、提高品牌知名度、促进销售转化，并与消费者建立长期关系。通过数据分析和用户反馈，企业能够实时调整营销策略，以适应市场变化和消费者需求。在大数据网络时代，网络营销面临着前所未有的变革与机遇。目前，主要的网络营销方式包括搜索引擎营销、微信营销、数据库营销和病毒式营销等[1]-[4]。

纵向数据是指在不同时间点上对同一组个体或观察单位进行重复测量的数据。这种数据类型常用于研究时间效应、个体变化和动态过程，能够揭示个体随时间的变化趋势。在大数据背景下，纵向数据因其特点在网络营销过程中占据了重要地位。例如，在收集某搜索引擎平台上黄金的销售数据时，这些数据不仅受到黄金公司自身在产品质量、样式和促销等方面的影响，还受到国家对黄金价格的调控政策的影响。如果仅依赖截面数据，虽然可以分析各黄金公司受到多种变量的不同影响，但无法揭示国家调控政策对黄金销售的具体影响。而如果仅使用时间序列数据，则只能分析黄金价格和销售情况随时间的变化，而无法有效体现不同公司在黄金销售中受到的具体变量影响[5]。

对于纵向数据通常采用混合效应模型建模，混合效应模型中包含固定效应和随机效应。而线性混合效应模型在一些情况中拟合效果较差。因此，许多统计学家开始采用半参数混合效应模型对纵向数据进行建模。半参数混合效应模型兼具参数模型和非参数模型的优点，在有效利用数据信息的同时考虑了信息不完整的变量，增强了模型的灵活性。Angelo等[6]通过引入B样条对光滑样条模型进行了优化，提出了一种新的方法以拟合非线性混合效应模型。阙烨等[7]采用广义最小二乘法对未知参数、随机效应及方差成分进行了估计，并证明了所提出估计量的渐进性质。为提高估计的效率，Lindley等[8]将贝叶斯方法应用于线性混合模型中，而齐培艳等[9]则研究了具有变点的半参数非线性混合效应模型的多重估计方法。

半参数混合效应模型中的参数估计通常有两种方法，一种是频率学派的极大似然方法，一种是贝叶斯学派的MCMC (Markov Chain Monte Carlo)方法。与频率学派方法相比，贝叶斯方法建模时最重要的优点就是不需要渐进逼近，由于所有参数均为随机变量，因此可以采用一些基于抽样的方法(例如蒙特卡洛方法)来构建复杂模型。同时，贝叶斯方法可以有效减轻计算负担，使得模型求解更容易[10]。

本文使用半参数混合效应模型对纵向数据进行建模，模型中纵向指标与时间的关系采用B样条逼近，并对模型中的所有未知参数、随机效应给出先验分布，基于贝叶斯MCMC方法对模型中的未知参数进行估计。最后将所提出的方法运用于网络营销影响因素的研究中，从而验证模型的有效性。

2. 模型

2.1. 半参数混合效应模型

假设有n个个体，每个个体之间相互独立，对每个个体分别进行$m_i$ 次观测，本文采用以下半参数混合效应模型：

$Y_{ij}=X_{ij}^{\text{T}}\beta +Z_{ij}^{\text{T}}{b}_i+f\left(t_{ij}\right)+{\epsilon }_{ij}$ (1)

其中，$i=1,\cdots ,n$ ；$j=1,\cdots ,m_i$ ；${\displaystyle \sum _{i=1}^n{m}_i=M}$ ；$Y_i\left(t_{ij}\right)$ 为第i个个体在第j个观测时间$t_{ij}$ 观察到的响应变量；$X_i^{\text{T}}\left(t_{ij}\right)$ 和$Z_i^{\text{T}}\left(t_{ij}\right)$ 分别是固定效应$\beta$ 和随机效应$b_i$ 的设计向量，固定效应$\beta ={\left({\beta }_1,\cdots ,{\beta }_p\right)}^{\text{T}}$ ，随机效应$b_i$ 服从多元正态分布$N_q\left(0,D\right)$ 且相互独立；$f\left(t_{ij}\right)$ 是一个未知的光滑函数；${\epsilon }_i\left(t_{ij}\right)$ 是第i个个体在第j个观察点$t_{ij}$ 的测量误差，假定${\epsilon }_i\left(t_{ij}\right)$ 服从正态分布$N\left(0,{\sigma }^2\right)$ 且与$b_i$ 相互独立。

2.2. B样条

参考Rupper等[11]，我们在贝叶斯框架下使用B样条来逼近模型(1)中的未知光滑函数。具体而言，我们考虑以下形式的光滑函数：

$f\left(t\right)={\lambda }_0^{\left(1\right)}+{\lambda }_1^{\left(1\right)}t+\cdots +{\lambda }_s^{\left(1\right)}{t}^s+{\displaystyle \sum _{l=1}^L{\lambda }_l^{\left(2\right)}{\left(t-{\tau }_l\right)}_{+}^s}$ (2)

其中，L为光滑函数的节点数，s为样条的自由度，将光滑函数划分为$L+1$ 个区间，我们考虑按照以下规则来选取节点和样条：假设${\tau }_l$ 为第l个节点，则选取样本的第$l+1/L+2$ 分位数，并满足${\tau }_0=t_{\min }<{\tau }_1<\cdots <{\tau }_l=t_{\max }$ ，${\left(t-{\tau }_l\right)}_{+}^s={\left(\max \left\{0,t-{\tau }_l\right\}\right)}^s$ ，令$\lambda ={\left({\lambda }^{\left(1\right)\text{T}},{\lambda }^{\left(2\right)\text{T}}\right)}^{\text{T}}={\left({\lambda }_0^{\left(1\right)},{\lambda }_1^{\left(1\right)},\cdots ,{\lambda }_k^{\left(1\right)},{\lambda }_1^{\left(2\right)},\cdots ,{\lambda }_L^{\left(2\right)}\right)}^{\text{T}}$ ，截断的幂基函数$\phi \left(t\right)={\left({\phi }_1{\left(t\right)}^{\text{T}},{\phi }_2{\left(t\right)}^{\text{T}}\right)}^{\text{T}}={\left(1,t,\cdots ,t^s,{\left(t-{\tau }_1\right)}_{+}^s,\cdots ,{\left(t-{\tau }_L\right)}_{+}^s\right)}^{\text{T}}$ ，则式(2)可写为：

$f\left(t\right)={\phi }_1{\left(t\right)}^{\text{T}}{\lambda }^{\left(1\right)}+{\phi }_2{\left(t\right)}^{\text{T}}{\lambda }^{\left(2\right)}=\phi {\left(t\right)}^{\text{T}}\lambda$

记$Y_i={\left(Y_{i1},\cdots ,Y_{i{m}_i}\right)}^{\text{T}}$ ，$Y={\left(Y_1^{\text{T}},\cdots ,Y_n^{\text{T}}\right)}^{\text{T}}$ ，以同样方式定义$X,Z,x,\epsilon$ ，设$Z=\text{diag}\left(Z_1,\cdots ,Z_n\right)$ 为对角元素为$Z_1,\cdots ,Z_n$ 的对角矩阵，则式(1)可以表示为：

$Y=X_{*}{\beta }_{*}+Zb+Wa+\epsilon$ (3)

其中，$X_{*}=\left(X,{\phi }_1^{\text{T}}\right)$ ，$W={\phi }_2^{\text{T}}{\left(t_{ij}\right)}_{M\times L}$ ，$a={\lambda }^{\left(2\right)}$ ，${\beta }_{*}={\left({\beta }^{\text{T}},{\lambda }^{\left(1\right)\text{T}}\right)}^{\text{T}}$ ，$\lambda ={\left({\lambda }^{\left(1\right)\text{T}},{\lambda }^{\left(1\right)\text{T}}\right)}^{\text{T}}~N_{s+1+L}\left(0,\xi \right)$ ，$\xi =\text{diag}\left({\xi }_1,{\xi }_2\right)$ ，${\xi }_1=\text{diag}\left({\sigma }_{{\xi }_1}^2,\cdots ,{\sigma }_{{\xi }_1}^2\right)$ ，${\xi }_2=\text{diag}\left({\sigma }_{{\xi }_2}^2,\cdots ,{\sigma }_{{\xi }_2}^2\right)$ 。

3. 贝叶斯分析

令$\theta =\left\{{\beta }_{*},D,{\sigma }^2,\xi \right\}$ 为未知参数的集合，为了对所有未知参数进行贝叶斯估计，需要从所有未知参数的联合似然函数中获得样本观测，由(3)式可得似然表达式如下：

$L\left(Y,b,a|\theta \right)={\left(2\pi {\sigma }^2\right)}^{-n/2}\times \exp \left\{-{\left(2{\sigma }^2\right)}^{-2}\times {\left(Y-X_{*}{\beta }_{*}-Zb-Wa\right)}^{\text{T}}\times \left(Y-X_{*}{\beta }_{*}-Zb-Wa\right)\right\}$ (4)

3.1. 先验的设置

联合模型的估计可以通过两种方式进行，一种是频率学派的极大似然方法，另一种是贝叶斯学派的MCMC方法。与频率学派方法相比，贝叶斯方法建模时最重要的优点就是不需要渐进逼近，且因为所有的参数都是随机变量，因此可以使用基于抽样的方法(如蒙特卡洛方法)来构建复杂模型。本文使用贝叶斯推理方法，基于指定先验分布的联合似然，通过马尔可夫链蒙特卡罗过程在完全贝叶斯框架下同时估计联合模型的参数。在贝叶斯框架下，为模型中的所有未知参数和随机效应指定先验分布如下：

$\begin{array}{l}{\beta }_{*}~N_{p+2}\left({\varphi }_1,{\varphi }_2\right),b_i~N_q\left(0,D\right),a~N_{r-2}\left(0,\tau I\right),\\ D~I{W}_q\left(o_o,R_o\right),\tau ~IW\left(o_1,R_1\right),{\sigma }^2~1/{\sigma }^2\end{array}$

其中，$IW$ 表示逆Wishart分布，${\sigma }^2$ 采用扩散先验分布，${\varphi }_1,{\varphi }_2,o_o,R_o,o_1,R_1$ 是超参数。

3.2. 条件后验的推导

(1) 为固定效应${\beta }_{*}$ 指定先验分布为${\beta }_{*}~N_{p+2}\left({\varphi }_1,{\varphi }_2\right)$ ，根据似然表达式(4)，通过多元配方法[12]，推导出${\beta }_{*}$ 的条件后验分布：

$\pi \left({\beta }_{*}|Y,b,a,{\sigma }^2\right)\propto \exp \left\{-\frac{1}{2}{\left({\beta }_{*}-\left({\varphi }_1-{\beta }_{*}\right)\right)}^{\text{T}}\left(\frac{1}{{\sigma }^2}{X}_{*}^{\text{T}}{X}_{*}+{\varphi }_2^{-1}\right)\left({\beta }_{*}-\left({\varphi }_1-{\beta }_{*}\right)\right)\right\}$

其中${\beta }_{*}={\left(X_{*}^{\text{T}}{X}_{*}+{\sigma }^2{\varphi }_2^{-1}\right)}^{-1}{X}_{*}^{\text{T}}{X}_{*}\left({\varphi }_1-\beta \right)$ ，可推导出${\beta }_{*}$ 的后验分布：

${\beta }_{*}|Y,b,a,{\sigma }^2~N_{p+2}\left(\left({\varphi }_1-{\beta }_{*}\right),{\left({\sigma }^{-2}{X}_{*}^{\text{T}}{X}_{*}+{\varphi }_2^{-1}\right)}^{-1}\right)$ (5)

(2) 为随机效应$b_i$ 指定先验分布为$b_i~N_q\left(0,D\right)$ ，令$w_i=Y_i-X_i{\beta }_{*}-Z^{\left(2\right)}a$ ，根据似然表达式(4)，通过多元配方法，推导出第i个个体的随机效应$b_i$ 的后验分布：

$\begin{array}{c}\pi \left(b_i|Y_i,\beta ,D,{\sigma }^2\right)\propto \exp \left\{-\frac{1}{2{\sigma }^2}{\left(b_i-{\widehat{b}}_i\right)}^{\text{T}}{\left(Z_1\right)}^{\text{T}}{Z}_i^{\left(1\right)}\left(b_i-{\widehat{b}}_i\right)\right\}\exp \left\{-\frac{1}{2}{b}_i^{\text{T}}{D}^{-1}{b}_i\right\}\\ \propto \exp \left\{-\frac{1}{2}{\left(b_i-b_i^{*}\right)}^{\text{T}}\left(\frac{1}{{\sigma }^2}{\left(Z_i^{\left(1\right)}\right)}^T{Z}_i^{\left(1\right)}+D^{-1}\right)\left(b_i-b_i^{*}\right)\right\}\end{array}$

其中，${\widehat{b}}_i={\left({\left(Z_i^{\left(1\right)}\right)}^{\text{T}}{Z}_i^{\left(1\right)}\right)}^{-1}{\left(Z_i^{\left(1\right)}\right)}^{\text{T}}{w}_i$，$b_i^{*}={\left({\left(Z_i^{\left(1\right)}\right)}^{\text{T}}{Z}_i^{\left(1\right)}+{\sigma }^2{D}^{-1}\right)}^{-1}{\left(Z_i^{\left(1\right)}\right)}^{\text{T}}\left(Y_i-X_i\beta -Z^{\left(2\right)}a\right)$ ，可推导出随机效应$b_i$ 的后验分布：

$b_i|Y_i,\beta ,D,{\sigma }^2~N_q\left(b_i^{*},{\left({\sigma }^{-2}{\left(Z_i^{\left(1\right)}\right)}^{\text{T}}{Z}_i^{\left(1\right)}+D^{-1}\right)}^{-1}\right)$ (6)

(3) 为随机效应a指定先验分布为$a~N_{r-2}\left(0,\tau I\right)$ ，根据似然表达式(4)，通过多元配方法，推导出a的条件后验分布：

$\pi \left(a|Y,{\beta }_{*},b,D,\tau ,{\sigma }^2\right)\propto \exp \left\{-\frac{1}{2}{\left(a-a^{*}\right)}^{\text{T}}\left(\frac{1}{{\sigma }^2}{\left(Z^{\left(2\right)}\right)}^{\text{T}}{Z}^{\left(2\right)}+{\left(\tau I\right)}^{-1}\right)\left(a-a^{*}\right)\right\}$

其中$a^{*}={\left({\left(Z^{\left(2\right)}\right)}^{\text{T}}{Z}^{\left(2\right)}+{\sigma }^2{\left(\tau I\right)}^{-1}\right)}^{-1}{\left(Z^{\left(2\right)}\right)}^{\text{T}}\left(Y-X_{*}{\beta }_{*}-Z^{\left(1\right)}b\right)$ ，可推导出随机效应a的条件后验分布：

$a|Y,{\beta }_{*},b,D,\tau ,{\sigma }^2~N_{r-2}\left(a^{*},{\left({\sigma }^{-2}{\left(Z^{\left(2\right)}\right)}^{\text{T}}{Z}^{\left(2\right)}+{\tau }^{-1}I\right)}^{-1}\right)$ (7)

(4) 为随机效应方差D指定先验分布为$D~I{W}_q\left(o_o,R_o\right)$ ，根据似然表达式(4)，通过多元配方法，推导出D的条件后验分布：

$\begin{array}{c}\pi \left(D|b\right)\propto {\left|D\right|}^{-\frac{m}{2}}\exp \left\{-\frac{1}{2}{\displaystyle \sum _{i=1}^m{b}_i^{\text{T}}{D}^{-1}{b}_i}\right\}{\left|D\right|}^{-\frac{o_0+q+1}{2}}\exp \left\{-\frac{1}{2}tr\left(D^{-1}{R}_0\right)\right\}\\ \propto {\left|D\right|}^{-\frac{m+o_0+q+1}{2}}\exp \left\{-\frac{1}{2}tr\left(D^{-1}\left(R_0+{\displaystyle \sum _{i=1}^m{b}_i{b}_i^{\text{T}}}\right)\right)\right\}\end{array}$

则随机效应方差D的条件后验分布：

$D|b_i~I{W}_q\left(m+o_0,R_o+{\displaystyle \sum _{i=1}^m{b}_i{b}_i^{\text{T}}}\right)$ (8)

(5) 为$\tau$ 指定先验分布为$\tau ~IW\left(o_1,R_1\right)$ ，根据似然表达式(4)，通过多元配方法，推导出$\tau$ 的条件后验分布：

$\pi \left(\tau |a\right)\propto {\left|\tau \right|}^{-\frac{r-2+o_1+1+1}{2}}\exp \left\{-\frac{1}{2}tr\left({\tau }^{-1}\left(R_1+a^{\text{T}}a\right)\right)\right\}$

可得到$\tau$ 的条件后验分布：

$\tau |a~IW\left(r-2+o_1,R_1+a^{\text{T}}a\right)$ (9)

(6) 为${\sigma }^2$ 指定先验分布${\sigma }^2~1/{\sigma }^2$ ，根据似然表达式(4)，推导出${\sigma }^2$ 的条件后验分布：

$\pi \left({\sigma }^2|{\beta }_{*},b,a,Y\right)\propto {\left({\sigma }^2\right)}^{-\left(\frac{n}{2}+1\right)}\exp \left\{-\frac{1}{{\sigma }^2}\left(\frac{1}{2}{u}^{\text{T}}u\right)\right\}$

其中$u=Y-X_{*}{\beta }_{*}-Z^{\left(1\right)}b-Z^{\left(2\right)}a$ 。

上述结果为逆Gamma分布，则推导出${\sigma }^2$ 的条件后验分布：

${\sigma }^2|Y,{\beta }_{*},b,a~IG\left(\frac{n}{2},\frac{1}{2}{u}^{\text{T}}u\right)$ (10)

3.3. 贝叶斯估计

式(5)~(10)中参数的后验分布为标准形式，可采用Gibbs抽样，Gibbs抽样是一种马尔可夫链蒙特卡洛(MCMC)方法，用于从多维概率分布中生成样本，具体步骤如下：

(1) 初始化：选择初始值${\theta }^{\left(0\right)}$ 作为参数的起始点。

(2) 迭代抽样：对于每个参数${\theta }_i$ (在$\theta =\left({\theta }_1,{\theta }_2,\cdots ,{\theta }_k\right)$ 中)，执行以下步骤：

a) 从条件分布$p\left({\theta }_i|{\theta }_{-i}\right)$ 中抽样，其中${\theta }_{-i}$ 表示除${\theta }_i$ 以外的所有参数；

b) 更新${\theta }_i$ 的值为新抽样的结果。

(3) 重复迭代：重复步骤(2)，直至达到预定的迭代次数或收敛标准。

设$\left\{{\beta }^{\left(k\right)},D^{\left(k\right)},{\tau }^{\left(k\right)},{\sigma }^{2\left(k\right)}:k=1,\cdots ,K\right\}$ 是由Gibbs抽样抽取的随机样本，则它们的贝叶斯估计可以表示如下：

$\widehat{\beta }=\frac{1}{K}{\displaystyle \sum _{k=1}^K{\beta }^{\left(k\right)}},\widehat{{\sigma }^2}=\frac{1}{K}{\displaystyle \sum _{k=1}^K{\sigma }^2{}^{\left(k\right)}},\widehat{D}=\frac{1}{K}{\displaystyle \sum _{k=1}^K{D}^{\left(k\right)}},\widehat{\tau }=\frac{1}{K}{\displaystyle \sum _{k=1}^K{\tau }^{\left(k\right)}}$

4. 实证分析

为了验证所提出方法的有效性，本节将基于所提出的半参数混合效应模型对饮食类非遗产品网络营销的影响因素进行分析。在166个国家级饮食类非物质文化遗产项目中，选择了在淘宝设有官方旗舰店的样本，总计80个，共进行了12次统计(2023年1月~12月) [13]研究，定义网络销量为响应变量Y，并选取淘宝官方旗舰店近半年综合销量排名第一的产品的月销量。以网络搜索量($X_1$ )、短视频传播效果($X_2$ )、店铺经营年限($X_3$ )、店铺粉丝数量($X_4$ )和店铺好评率($X_5$ )作为解释变量建立如下半参数混合效应模型：

$Y_{ij}={\beta }_0+{\beta }_1{X}_{1ij}+{\beta }_2{X}_{2ij}+{\beta }_3{X}_{3ij}+{\beta }_4{X}_{4ij}+{\beta }_5{X}_{5ij}+f\left(t_{ij}\right)+b_i+{\epsilon }_{ij}$ (11)

其中，$i=1,\cdots ,80$ ，$j=1,\cdots ,12$ ，$f\left(t\right)$ 是一个以t为自变量的函数，$b_i$ 是随机截距，服从正态分布$N\left(0,D\right)$ 且相互独立，${\epsilon }_{ij}$ 为误差项，假定其服从正态分布$N\left(0,{\sigma }^2\right)$ 。

将数据代入模型(11)，通过贝叶斯方法估计出模型中的参数，估计结果见表1。

Table 1. Parameter estimation results

表1. 参数估计结果

参数	估计值	标准差
${\beta }_0$	24.833	0.004
${\beta }_1$	0.652	0.142
${\beta }_2$	0.877	0.084
${\beta }_3$	0.538	0.008
${\beta }_4$	1.083	0.001
${\beta }_5$	1.314	0.001
$D$	0.217	0.023
${\sigma }^2$	0.209	0.019

从表1可以看出，截距项为24.833，网络搜索量系数为0.652，代表网络搜索量对于店铺的网络营销量有一个正向影响，即网络搜索量越高，店铺营销量越高。短视频传播效果系数为0.877，代表短视频传播效果对于网络营销量有一个较强的正向影响，即短视频传播效果越好，营销量越高。店铺经营年限系数为0.538，代表店铺经营年限对于店铺的网络营销量有一个正向影响，即店铺经营年限越长，网络营销量越高。店铺粉丝数量系数为1.083，代表店铺粉丝量对店铺的网络营销量有一个较强的正向影响，即店铺粉丝量越多，店铺的网络营销量越高。店铺好评率系数为1.314，说明店铺好评率对店铺的网络营销量有一个很强的正向影响，即店铺好评率越高，店铺的网络营销量越高。以上情况数据合理地解释了网络营销中的部分营销因素，与实际情况大致相符。

5. 结束语

本文构建了一个纵向数据的半参数混合效应模型，并运用贝叶斯方法进行了统计推断。研究发现，模型中的固定效应、随机效应及其方差，以及在先验分布下的随机误差的后验分布均符合标准分布，从而有效简化了模型估计的统计计算。同时，所提出的模型被应用于分析网络营销中的营销因素，实证结果显示这些因素对淘宝店铺的网络销售量产生了不同程度的影响。实证分析结论也与网络营销的实际情况相符，说明了模型的可行性和有效性。但是在实际情况中，店铺的网络营销量不仅受到文中所述的解释变量影响，还受到诸多因素的共同影响，如：支付方式、物流与配送、品牌信誉、市场营销策略等。如果忽视这些因素，可能会导致模型的准确性下降。因此，若能对这些因素进行量化并将其纳入模型进行计算，结果将会更加精准。这也是本文未来研究的一个重要方向。