翻译 | 反应时数据的多种方式 (附R代码)

[ray306]

(原文地址在最后)

原文链接：https://lindeloev.github.io/shiny-rt/

校译：COSN 翻译小分队成员↓↓

游小波 广西大学

Mel 京都大学

Blitz? 波士顿学院

杨姝 香港中文大学

郭治辰 江西师范大学

刘庆芳 俄亥俄州立大

二次校对：Ted

研究者们经常使用正态分布来拟合反应时数据。但其实我们有更好的选择，希望这篇文章可以减轻使用其他分布拟合反应时数据时面临的负担，因为正态分布对于反应时数据的拟合非常糟糕（如你将会看到的），转向其他分布的回报将会很大！

下面的速查表提供了一个概述（PNG 和 PDF 格式）。你也可以直接跳到交互部分（或在浏览器中打开：https://lindeloev.shinyapps.io/shiny-rt/）或者应用代码范例 (https://lindeloev.shinyapps.io/shiny-rt/#section-code)。

1 观点：三种重要的参数

对于反应时，笔者认为可以用三种参数来描述它们的分布：难度（Difficulty），偏移（Onset）和尺度（Scale），在下图中，蓝色曲线显示的是提高这些参数后的影响。读者可以在原网页 “3 Drag those sliders” 部分进行操作。

1.1?难度 (Difficulty)

反应时的一个主要特征就是其标准差随着平均值的增长呈近似的线性增长（Wagenmakers et al., 2007）。当反应时较长时（比如 stroop task 中的不一致试次 / incongruent trials），其分布会更分散；当反应时较短时，其分布会更集中（比如 stroop task 中的一致试次 / congruent trials）。从方便解释的角度来讲，一个能同时捕捉平均值和标准差变化的参数要比分别使用平均值和标准差这两个参数更方便，因为使用均值和标准差时，我们通常需要分别对它们进行线性回归并解释它们的交互效应。理想情况下，难度参数的值应该在分布的中心附近（最高密度区域 / highest density region，例如反应时间的中位数或众数）。

1.2 起始位置（Location onset）和中心 (Center)

这两个参数决定了反应时分布的左移（即更长反应时）和右移（即更短反应时）而不改变其分布形状。

（1）起始 (Onset)：这可能是反应时的开始阶段（在描述性模型中一般称为偏移 / shift）或者是决策过程的最早开始时间（在机制模型中一般称为非决策时间 / non-decision time）。对于人类而言，这个时间一般在 200 毫秒左右。为什么要设置这样一个参数？因为所有动物的反应时都得大于 0.0 毫秒。没有起始参数的模型会赋予更高的可能性(positive probability) 给短的出奇的反应时(impossibly short RTs)。

（2）中心 (Center)：靠近反应时分布中间的某处（比如正态分布中的均值μ）。笔者认为可以采用类似 difficulty 那样的参数(a difficulty-like parameter) 来描述反应时分布的中间。

1.3 范围 (Scale)

尺度是指作用于难度之外的对分布的拉伸或压缩，且不会（严重地）改变分布的位置。这可以用来调整分布以适应特定人群或特定实验的特征。例如，你可以采用形状更平缓的反应时分布来拟合患者样本的数据和双任务范式的数据。

1.4 哪种分布更好?

理想情况下，一个参数应该只需要描述一种分布特征。如果一个参数可以同时描述好几种分布特征，或者多个参数同时对应一种分布特征，这种情况我把它们称之为混乱 (messy)，因为你需要调整多个参数来实现上述效果。这种混乱会让人难以理解。不过，如果你有一个潜在的理论来驱动这种行为 (messy)，并且在你的理论中每一个参数都有意义，那么在这种情况下，这些参数可以称之为机制参数 (mechanism parameters)

有了这个术语，这里是我们的三条准则:

• 不好的反应时分布: 带有混乱参数 (messy parameters) 的分布。

• 好的反应时分布: 带有起始 (Onset)、难度和尺度参数的分布。

• 好的反应时分布: 带有机制参数的分布。

根据这些标准，DDM（drift diffusion model）、LBA(Linear Ballistic Accumulator)、(Shifted)Log-normal distribution 和 Inverse Gaussian distribution 都可以算是好的反应时分布。如果我们还考虑解释参数取值的难易程度，我认为 DDM 和偏移对数正态分布是最优的，当然这有待讨论。

2 选择数据

以下所有数据集均来自单一条件下的单一被试。

选择一个数据集。

辨别单词 / 非单词——准确（Wagenmakers et al., (2008)）

... 或粘贴你自己的 RTs（分隔符 = 空格 /,/;/tab）

最小 RT 最大 RT 组距

0.18

2

70

3 选择数据

这里是一些有趣的挑战：

• **单参数冠军：**用速度 / 准确度和 MRT 数据，尝试（1）将参数拟合到一个数据集，（2）改变数据集，（3）现在看看只改变一个参数你可以达到多接近数据。

• **拟合冠军：**如果你可以很好地拟合数据使得均方误差（MSE）低于 40，你将获得一个特别奖励。用键盘的方向键来微调参数。为了提供易读的数字，图中的均方误差被乘了 10.000。

• **最优冠军：**如文中代码段展示的那样使用 brms 来找到最优拟合。实证数据集在 rtdists 包里。使用 library(rtdists’)，data(‘rr98’)，和 data(‘speed_acc’) 来载入数据。在这里下载 MRT 数据。你也尝试用（频率学派的 / frequentist approach，即计算 p-values）gamlss 包，尽管这样结果解读起来会更难。

3.1 高斯分布（正态分布）(normal distribution)

高斯分布的特征：单个平均值加上向两边对称分散的分布，对应更快或者更慢的反应时。**多数情况下，正态分布是用来拟合反应时数据的最糟糕分布。**比如说，对于典型的反应时数据，在基于正态分布的拟合模型中会有 5% 的反应时比真实反应时数据中最短的反应时还要快！

但是，高斯模型是到目前为止在分析反应时数据中中用得最多的，因为它是大多数统计模型 (比如 t 检验，方差分析，线性回归)，代码包和分析软件(比如 lm, lme4::lmer, brms::brm) 的默认选择分布。

另外，如果你的预处理中包括给每个被试或者每种实验条件 (condition) 计算平均值，并且你把那些平均值当作数据来使用，那么你就已经使用了高斯模型（并且扔掉了很多有用的数据）。

参数：

(1) μ（中心）：平均值。这个受分布两侧尾端的影响大。

(2) σ（尺度）：标准差。这个描述了反应时数据的 (对称) 离散程度。

通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest) 中获得对这些参数值的估计。在这里的示例模型中，我们把实验条件(condition) 作为预测变量并以此估计μ (mu) 值，同时保持σ(sigma) 的恒定。详细内容见建模笔记和应用范例。

library(brms)

fit = brm(rt ~ condition, rt_data, family=gaussian())

3.2 指数高斯模型 (Ex-gaussian)

如果你认为反应是由两个独立的过程混合造成的，即（1）一个高斯模型 和（2）一个指数衰减模型 (decaying exponential)，那么这就是所谓的指数高斯模型。大多数人不这么认为，但他们却往往将指数高斯模型当成一种纯描述性模型来使用，因为它通常对于反应时数据有很好的拟合。见维基百科。

参数:

(1)?μ（中心）：高斯分布的平均值。描述了更短或更长的平均反应时。

(2)?σ（尺度）：高斯分布的标准差。这是对μ左右的对称离散性建模。

(3)?λ（混乱） ：指数模型的衰减率。这是对长反应时间的尾端分布建模。更高的λ意味高斯分布中尾端分布占有更多主导。

(4)?通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest) 中获得对这些参数值的估计。在这里的示例模型中，我们把实验条件 (condition) 作为预测变量并以此估计μ (mu)值，同时保持σ(sigma)的恒定。详细内容见建模笔记和应用范例。

library(brms)

fit = brm(rt ~ condition, rt_data, family=exgaussian())

3.3 偏态分布 (Skew normal)

偏态分布就是给高斯分布新加一个描述偏态程度的参数。根据笔者的经验，偏态分布对反应时数据的拟合通常较差，虽然它好于高斯分布。

参数：

（1）?μ（中心）：平均值。

（2）?σ（尺度）：标准差。改变σ不会影响原有的平均值μ，但会改变中位数。注意当α≠0 的时候，分布不是对称的，所以不要将它解读为正态分布的标准差。

（3） α（混乱）：偏度。当α=0，这是一个高斯分布。

（4）通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest) 中获得对这些参数值的估计。在这里的示例模型中，我们把实验条件 (condition) 作为预测变量并以此估计μ (mu)值，同时保持σ(sigma)的恒定。。。详细内容见建模笔记和应用范例。

library(brms)

fit = brm(rt ~ condition, rt_data, family=skew_normal())

3.4 (偏移) 对数正态分布 (Shifted Log-normal)

使用对数时，你是对增加或减少的百分比而非绝对值差异 (absolute difference) 进行建模。也就是说，此时你的参数含义不是 “RTs 在不一致试次慢 40 毫秒”，而是“RTs 在不一致试次中慢 7%”。这比较符合人们的直觉: 当反应时本身较短时，40 毫秒的差异是巨大的，但这种差异对于 10 秒左右的 RTs 来说却微不足道。从数学上讲，它是一种乘积关系(multiplicative relationships) 而不是加法关系(additive relationships)，而对于讨论 RTs，这比较说得通。所以笔者认为对数正态分布是介于描述性模型和机制性模型之间的。一般来说，当因变量只能取正值时(positive-only values)，比如反应时数据，对数模型往往是非常适用的.

为什么叫 “对数正态分布”？是因为其参数是经对数变换的 RTs 的平均值和标准差，而我们假设这些经过对数转换的参数都服从正太分布 (即高斯分布)。例如，对于简单的 RTs，您会经常看到接近 exp(-1.5) = 0.22 秒左右的截距。

参数:

（1）μ(难度) : 对数正态分布的平均值。μ的均值代表反应时数据的中位数。换句话讲，反应时数据的中位数 = shift+exp(μ)（译者注：偏移量加上难度系数的指数）。

（2）σ(尺度) : 对数正态分布的标准差。会增加μ（译者注：这里μ是指难度参数，即对数正态分布的均值，而非 RT 的均值）的均值但不会增加μ的中位数。

（3）shift(偏移) : 最短开始反应的时间。当偏移 = 0 时，则是一个只带有两个参数的直上的（译者注：即从 t=0 就立即开始反应）对数正态分布。

通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest)中获得对这些参数值的估计。在这里的示例模型中，我们把实验条件 (condition) 作为预测变量并以此估计μ (mu)值，同时保持σ(sigma)的恒定。。。详细内容见建模笔记和应用范例。??

**library**(brms)

fit = brm(rt ~ condition, rt_data, family=

shifted_lognormal())

3.5 (偏移) 逆高斯分布（Shifted Inverse Gaussian）

也被称为 (偏移的) 瓦尔德（Wald）。它并不是倒高斯分布。它与维纳 / 决策扩散模型 (Weiener/Decision model) 和 LBA 模型密切相关(引自维基百科）。

“这个名字可能会引起误解: 它有 “倒” 称谓，只是因为高斯分布描述了一个固定时间的布朗运动的水平，而倒高斯分布描述了一个带有正向漂移的布朗运动达到一个固定正值水平所需时间的分布”

换句话说，如果你用 X 轴表示时间，并且在 t0=shift 处开始了很多 “随机漫步” 高斯函数描述了它们在一个固定时间 (t_i) 的垂直位置分布，而倒高斯函数描述了它们越过 y 轴上某个阈值时的时间分布。像维纳和 LBA 一样，倒高斯函数可以被看作是对一个决策过程的建模，而不需要整合决策反应的正确性。

参数:

（1）μ?(难度) : 均值。注意，此参数受到长 RTs 的影响很大。

（2）λ?(尺度) : 减小方差，但很难解读，因为σ=(μ3/λ)。

（3）shift (起始时间) : 可能的最早反应时。当 shift=0 时，它是纯倒高斯分布或瓦尔德分布，即只有两个参数

我们可以很容易地推断出 (无偏移的) 逆高斯 /  瓦尔德函数的参数。通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest) 中获得对这些参数值的估计。在这里的示例模型中，我们把实验条件 (condition) 作为预测变量并以此估计μ (mu)值，同时保持σ(sigma)的恒定。详细内容见建模笔记和应用范例。

**library**(brms)

fit = brm(rt ~ condition, rt_data,

family=inverse.gaussian())

3.6 维纳扩散 (Wiener diffusion)

也称为扩散决策模型 (DDM)，它可能是最常见的用来描述机制的模型。它适用于结果变量为二选一迫选的数据，并且要求数据中既包括反应时也包括正确率(accuracy, scored as correct/wrong)。它将 RTs 建模为当随机漫步者跨越两个阈值(边界) 之一时所做的决策。很明显，机制性解释性是一种优势，但理解这个模型需要花一些时间和精力。阅读 Ratcliff 和 McKoon(2008)以获得关于这个模型的细致介绍。笔者本人喜欢这种视觉化呈现（译者注：读者可进入原文链接进行操作）。

值得一提的是，扩散的概念也是构成倒高斯分布和 LBA 模型的基础。

参数:

所有这些都是机制性的：

（1）α(难度): 边界分离。阈值之间的距离。距离越远 - > 路径变异越大 - >RTs 越分散。

（2）τ(起始): 非决策时间 (ndt)。随机漫步者(random walker) 开始的时间。

（3）β(难度): 偏量。漫步者的起始点到上阈值间的距离 (β=0.5 即为中间值)。参数值越高 = 起始点越接近上阈值 = RTs 越快越不离散。

（4）δ(难度): 漂移速率。漫步者速度。参数值越高，RTs 越快越不离散。

（1）如果你想在 brms::brm 中使用该模型来拟合反应时数据，请参考这一篇指南 (http://singmann.org/wiener-model-analysis-with-brms-part-i/)。不过可能需要花上几天甚至数周的时间。

（2）除了 brms::brm，你还可以使用 Wagenmakers(2007) 等人开发的 EZ diffusion calculator (https://www.ejwagenmakers.com/EZ.html)。

（3）如果你想在 Python 当中使用该模型拟合反应时数据，那你可以使用 HDDM python package (http://ski.clps.brown.edu/hddm_docs/)。

3.7?Linear Ballistic Accumulator

这个模型与维纳扩散模型有一些相似之处，因为它模拟了向决策阈值逼近的潜在 "运动"，而且非常注重机制。它同样适用于強制性二选一（two-alternative-force-choice）的情況，同样既需要反应时数据也需要正确率的数据 (scored as correct/incorrect response)）。虽然它涉及到很多参数，但如 Brown & Heathcote（2008）解释的，模型本身（至少应该是）比扩散模型更简单。该论文免费阅读见此处

参数：

所有参数都是机制性的。

（1）A: 起点的上界。起点值为在这个区间内随机抽取的值。

（2）b（难度）：反应阈值（通常，b>A）。

（3）t0（偏移）："运动" 启动前的时间，即非决策时间。

（4）mean_v（难度）：漂移速率的均值以及此均值的离散性（SD）。

（5）sd_v：漂移速率的平均标准差及此平均标准差的离散性（SD）。

使用 glba::lba() 对其进行拟合。在 rtdists 包中有一些有用的分布。

3.8 伽马分布（Gamma）

Gamma 分布在统计学中非常重要，也延展到了 RT 研究中。虽然它可以拟合数据，但很难解释涉及到 RT 的参数，笔者这里暂不详述。理论上讲，人们可以很容易地做出一个偏移的 Gamma (shifted Gamma)，但笔者尚未见过它被用于反应时数据的分析。

参数：

这两个参数单独来看都是类似难度的，所以它们的相互作用决定了尺度和位置。笔者认为这使得它们变得混乱。

（1）α（混乱的）：一个描述分布形状参数。

（2）β（混乱的）：一个描述分布频率参数。

通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest) 中获得对这些参数值的估计。在这里的示例模型中，我们把实验条件 (condition) 作为预测变量并以此估计 μ (mu)值，同时保持σ(sigma)的恒定。。。详细内容见建模笔记和应用范例。

详见建模说明和应用实例。

**library**(brms)

fit = brm(rt ~ condition, rt_data,

family=weibu11())

3.9 (偏移) 威布尔 (Shifted Weibull)

...... 当其存在一个偏移参数时，偏移 Weibull 也被称作三参数 Weibull。RT 可以被理解为存在一个潜在的速率，随着时间的推移，潜在的非反应过程会以此数率 “转变成” 反应。如果λ=0，这个速率是不变的。如果λ>1，速率会随着时间的推移而降低。不过参数值并不对应任何直观的东西。

参数：

（1）λ (难度): 因为某些原因，它被称为尺度参数， 但它影响平均值，中位数，和众数。

（2）k（混乱的）。一个描述分布形状参数。

（3）偏移 (起始时间)：当 shift=0 时，就是一个带有上述两个参数的直上的 Weibull 分布。

从数据中很容易推断出（非偏移的）Weibull 的参数。这里，我们通过条件对λ进行建模。它是对数，所以使用 exp() 进行反向变换。k 在各个条件下是固定的。详见建模说明和应用实例。

**library**(brms)

fit = brm(rt ~ condition, rt_data,

family=gamma())

4 应用代码范例

我们载入 Wagenmakers & Brown (2007) 的数据然后预处理它：

# Load data_（载入数据）

_library(rtdists)

data(speed_acc)

# Remove bad trials_（除去不好的试次）

_library(tidyverse)

data = speed_acc %>%?

filter(rt> 0.18, rt < 3) %>%? # Between 180 and 3000 ms_（在 180-3000 毫秒之间）_?

mutate_at(c('stim_cat', 'response'), as.character) %>%? # from factor to char_（从 factor 到 character）

filter(response != 'error', stim_cat == response) # Only correct responses（只保留正确反应试次，去掉所有错误反应试次）_

现在我们用偏移对数正态分布来问这样一个问题：“反应时间的中位数在各个实验条件中差异显著吗？” 换句话说，通过使用贝叶斯模型并结合数据，我们可以从后验分布 (posterior distribution for parameters of interest) 中获得对μ此外让尺度σ和起始偏移保持不变。我们将使用贝叶斯混合模型(Bayesian mixed models)，包括一个被试水平上的随机截距(a random intercept per participant)。

见笔者的反应时间模型公式笔记 或者更详细的信息。在这个例子中，我们有大约 28000 个数据点所以运行代码需要花点时间。

library(brms)

fit = brm(formula = rt ~ condition + (1|id),

dat = data,? ? ? ? ?

family = shifted_lognormal(),?

file = 'fit_slog')? # Save all that hard work.__（保存所有工作。）

完成！

4.1 检验模型是否拟合

和上面的互动小程序差不多，我们可以看到参数是如何预测真实数据的：

_#plot(fit)(_给拟合作图) _# Check convergence and see posteriors(__检查收敛和后验分布)

_pp_check(fit)? # Posterior predictive check （后验预测检查）

这比高斯分布要好得多，但仍然有改进空间。这基本上是因为偏移（shift）接近整个数据集中最短的反应时间，这个最短反应时间是某个被试在某一条件中的反应时，0.18 秒（被试 id #8）——刚刚高于数据预处理后的下限。平均的最短反应时间是 0.28 秒，因此这名被试是没有代表性的。我们可以在混合模型中给每个被试设定一个偏移（shift）更合理：

formula = bf(rt ~ condition + (1|id),

ndt ~ id)

另外，我们还可以去掉每个被试最短反应时试次，如果你可以说服其他人这些最短反应时是意外反应的话：

data %>%

group_by(id) %>%

arrange(rt) %>%?

slice(4:n())? # Remove t__hree shortest RTs_（对每个被试都除去三个最短的反应时试次）_

4.2 了解参数估计

先把这个放在一边，我们来看看我们对偏移对数正态模型的参数有了怎样的了解。距和条件速度一起定义了μ（难度）你也可以看到 sigma(尺度) 和 ndt(开始)，所以我们的三个参数都有后验分布 (posterior distribution)!

summary(fit)??#?See?p__arameters

Population-level?Effects:

Estimate?Est.Error?l-95%?CI?u-95%?CI?Rhat?Bulk_ESS?Tail_ESS

Intercept?????????-0.69??????0.04????-0.76????-0.61?1.01??????508??????672

conditionspeed????-0.40??????0.00????-0.41????-0.39?1.00?????2184?????2143

Family?Specific?Parameters:

Estimate?Est.Error?l-95%?CI?u-95%?CI?Rhat?Bulk_ESS?Tail_ESS

sigma????0.37?????0.00????0.36????0.37?1.00?????1839?????2088

ndt??????0.17?????0.00????0.17????0.18?1.00?????1576?????1724

因为μ 和σ定义了对数正态分布，即他们是对数级的单位。这意味着，对 “精度”(截距) 条件，我们预期 ndt+exp(μ)=0.17+exp(?0.69)=0.50 秒的中位数反应时，对 “速度” 条件，预期 exp(- 69 -0.40) = 0.34 秒的中位数反应时。如果您对这些转换感到困惑，不用担心。brms 可以帮你恢复，不用考虑其分布:

marginal_effects(fit)?# Back-transformed parameter estimates

marginal_effects(fit, method='predict')? # Same, but for responses

我们非常肯定被试在 “速度” 条件下的反应更快，因为置信间的上界远小于零。我们可以更直接地测试:

hypothesis(fit, '**conditionspee****d** **< 0**')?

_# Numerically_

4.3 先验值和描述性统计

为了简洁起见，我在这里没有使用默认的先验值，而是指定了一些合理的先验值，然后运行 brm(..., prior=prior)。主要是作为正则化。数据太多，会导致大多数先验值被淹没。

# 看看你可以为这个模型设置哪些先验值。

get_prior(rt ~ condition + (1|id), speed_acc, family = shifted_lognormal())

#定义它们

prior = c(

set_prior('normal(-1, 0.5)',), # 大约 exp(-1) = 0.36 秒。

set_prior('normal(0.4, 0.3)',), #单个 RTS 的 SD，以对数单位表示。

set_prior('normal(0, 0.3)',), # “大约快或慢 exp(-1)-exp(-1 + 0.3)=150 毫秒的反应時应? set_prior('normal(0.3, 0.1)',) # 参与者之间存在一定的差异。)

我可能还会通过绘制几个参与者的反应时分布图来检查预处理是否有效。

data_plot = filter(data, id%in % c(1, 8, 15))

ggplot(data_plot, aes(x=rt)) +。? ?geom_histogram(aes(fill=condition), alpha=0.5, bins=60) +。

facet_grid(~id) + # 每个 id 有一个面板。

coord_cartesian(xlim=c(0, 1.6)) # 放大

5 关于好的反应时间模型的最后说明

贯穿整篇文章，我演示代码用了一个简单的模型：rt ~ condition。对于真实世界的数据，你大概会对数据产生的过程了解得多很多。下面是一个更好的典型方法：

rt ~ condition + covariates + (condition | id) + (1|nuisance)

下面是一个使用此方法的例子：

rt ~ congruent + trial_no + bad_lag1 + (congruent|id) +

(1|counterbalancing)

我们来解读一下：

condition（条件）是我们感兴趣的对照。这可以被看作一个一致或不一致的变量，视觉搜索中项目的数量，或是其他的变量。它也可以是一个交互效应。我刚做了一个分析其中包含 congruent:condition 这一项来对每种条件的一致性效应建模。

covariates（协变量）可以是，举例来说，trial_no + bad_lag1，换个说法就是时间效应（人们一般随练习而变得更快）和每组测试中的第一次测试所花的较长的时间（你需要给这个变量编码）。通常对协变量进行缩放（scale）和居中（center，对于连续和二元变量都可能）会比较好，前者有助于模型拟合，后者使得 condition 产生的效应代表其他条件的平均值，而不是一个各种条件的某种组合。你可以用 scale 来进行缩放和居中。

(1 | nuisance) 是一种让模型 "理解" 数据中的结构，而不影响固定参数估计的方法。如果你进行的是因子设计，那些非模型设计变量一般会放在这里。

注：

1."难度" 是借用 Wagenmakers & Brown(2007) 的术语，但它不是标准术语。

2. 有些分布是用速率作为参数的，速率是位置中心的倒数。

3. 许多分布都包含形状参数，但由于所有 RT 分布的形状都大致相同，所以我更愿意选择一个本来就具有 "正确" 形状的分布，而不是使用一个形状参数来建模。

排版：李会玲、王薇薇

审核：胡传鹏

原文地址