第十三章 实验设计与方差分析

参考书目为安德森的《商务与经济统计》,以下为个人的学习总结,如果有错误欢迎指正。有需要本书pdf的,链接在本文末尾。(仅限个人学习使用,请勿牟利)

第十三章 实验设计与方差分析

统计研究分实验性研究和观测性研究。前者需要控制无关变量,通过实验产生我们需要的数据,后者往往通过抽样调查等方式获得。
本章介绍三种类型的实验设计:完全随机化设计、随机化区组设计和析因实验。

13.1 实验设计和方差分析简介

例子:供水过滤系统的部件组装方法有A、B和C。问题:哪种方法使每周产量最多。

在这个实验中,装备方法是独立变量因子(factor)。对应三种方法,所以这个实验有三个处理,每个处理(treatment)对应一种装配方法。并且是单因子实验(single-factor experiment),因为只涉及装配方法一个因子。也可以有多因子,因子分定性和定量的。

该实验对应三个总体:三个总体分别使用A、B和C其中一种方法。每个总体的因变量响应变量是每周装配的过滤系统的数量。

实验目的:确定三个总体的因变量是否相同。
假设我们抽取三名工人组成一个随机样本,三名工人构成实验单元,下面将使用完全随机化设计(completely randomized design),要求每种方法随机给其中一个工人,这里相当于工有A_3^3种分配方法。(随机化的概念是所有实验设计的一个重要原则

上述方法,每个装配方法只能得到一个因变量的测度,但是我们可以随机抽15个人,每种方法随机分5人。这样就得到了更多因变量的测度。这个过程叫复制。(复制的过程是实验设计的另一个重要原则。

13.1.1 数据收集

通过收集数据得到

image

我们为了了解三种装配方法的总体均值是否不同,引入\mu_1\mu_2\mu_3分别为A、B和C三种方法每周生产数量的总体均值。
假设:H_0:\mu_1=\mu_2=\mu_3 H_a:总体均值不全相等
再使用方差分析(ANOVA)统计方法来确定三个样本均值之间的差异是否大到可以拒绝H_0

13.1.2 方差分析的假定

应用方差分析需要三个假定:

  1. 对每个总体,响应变量服从正态分布。——每种装配方法每周生产的过滤系统数量服从正态分布。
  2. 响应变量的方差,记为\sigma^2,对所有总体都是相同的。——每种装配方法,每周生产的过滤系统数量的方差必须相同。
  3. 观测值必须是独立的。——对应每个工人每周生产的过滤系统数量与其它工人每周生产的过滤系统数量独立。

13.1.3 方差分析:概念性综述

样本均值彼此接近,则越支持H_0,反之支持H_a
如果原假设(H_0:\mu_1=\mu_2=\mu_3)成立,我们利用样本均值之间地变异性简历\sigma^2的一个估计。则所有样本都来自同一个总体。这些样本均值\bar x同样服从正态分布,且均值为\mu,方差为\sigma^2/n
回到过滤系统的例子中,我们假设\bar x_1=62,\bar x_2=66,\bar x_3=52都来自同一个总体(样本容量相同),\bar x抽样分布的均值的估计值为:(62+66+52)/3=60\bar x抽样分布的方差\sigma_{\bar x}^2的估计可以由三个样本均值的方差给出\sigma_{\bar x}^2=\frac{(62-60)^2+(66-60)^2+(52-60)^2}{3-1}=52
再由\sigma_{\bar x}^2=\sigma^2/n解得\sigma^2=n\sigma_{\bar x}^2=5 \times 52=260因为\sigma_{\bar x}^2是用s_{\bar x}^2作为估计量,所以这里得\sigma^2也是估计量。
所得的结果n\sigma_{\bar x}^2=260称作\sigma^2的处理间估计。

上述都是基于H_0为真的情形,如果H_0为假,且均值全不相同,则三个抽样分布来自三个总体。于是\sigma_{\bar x}^2会比较大,从而使得\sigma^2的处理间估计也变得较大。

image

当我们从每个总体抽取一个随机样本时,每个样本方差都给出了\sigma^2的一个无偏估计,我们将\sigma^2的个别估计组合或合并成一个总体估计。这种方法得到值称作\sigma^2的合并估计或处理内估计。因为这里的每个样本方差给出的\sigma^2的估计仅以每个样本内部的变异为依据。
\sigma^2的处理内估计=\frac{27.5+26.5+31.0}{3}=\frac{85}{3}=28.33

我们看到\sigma^2的处理间估计(260)远大于处理内估计(28.33),比值为9.18。
当原假设为真,处理间估计方法才是总体方差\sigma^2的一个好的估计量,
当原假设为假,处理间估计将高估总体方差\sigma^2
不过这两种情形下,处理内估计都是总方差\sigma^2的一个好的估计量。因此原假设为真,两估计量接近,比值接近1;如果原假设为假,则处理间估计将大于处理内估计,比值也会比较大。

总结
ANOVA背后的逻辑是以共同总体方差\sigma^2的两个独立的估计量为基础,即处理间估计和处理内估计。通过比较两个估计量,来确定总体均值是否相等。

13.2 方差分析和完全随机化实验设计

完全随机化实验设计中,如何用方差分析来检验k个总体均值是否相等:

  • 一般性假设:H_0:\mu_1=\mu_2=\mu_3 H_a:k个总体均值不全相等
  • \mu_j为第j个总体的均值,n_j为第j个总体的简单随机样本的容量,
  • 样本数据x_{ij}代表第j个处理(对应一个总体)第i个观测值;\bar x_j代表第j个样本的均值,s_j^2代表第j个处理的样本方差,s_j为第j个处理的样本标准差。
  • 样本数据的计算
    • \bar x_j=\frac{\sum\limits_{i=1}^{n_j} x_{ij}}{n_j}
    • s_j^2=\frac{\sum\limits_{i=1}^{n_j}(x_{ij}-\bar x_j)^2}{n_j-1}
  • 总体样本均值记作\bar{\bar x}=\frac{\sum\limits_{j=1}^{k}\sum\limits_{i=1}^{n_j}x_{ij}}{n_r}
    • 观测中总数n_r=n_1+n_2+\cdots+n_k
    • 如果每个样本容量相等,则n_r=kn,则\bar {\bar x}=\frac{\sum\limits_{j=1}^{k} \bar x_j}{k}

13.2.1 总体方差的处理间估计

我们称处理间估计的\sigma^2均方处理(mean square due to treatments, MSTR)
MSTR=\frac{\sum\limits_{j=1}^{k}n_j(\bar x_j-\bar {\bar x})^2}{k-1}
式中分子称作处理平方和(sum of squares due to treatments, SSTR)。分母k-1表示与SSTR相联系的自由度。
均方处理
MSTR=\frac{SSTR}{k-1}
SSTR=\sum\limits_{j=1}^{k}n_j(\bar x_j-\bar {\bar x})^2

H_0为真,则MSTR给出了\sigma^2的一个无偏估计。但H_0为假时,则MSTR就不是\sigma^2的无偏估计,会高估总体方差\sigma^2

回到例子:
SSTR=\sum\limits_{j=1}^{k}n_j(\bar x_j-\bar {\bar x})^2=5\times(62-60)^2+5\times(66-60)^2+5\times(52-60)^2=520
MSTR=\frac{SSTR}{k-1}=260

13.2.2 总体方差的处理内估计

\sigma^2的处理内估计称作均方误差(mean square due to error,MSE)
MSE=\frac{\sum\limits_{j=1}^{k}(n_j-1)s_j^2}{n_r-k}
分子称作误差平方和(sum of squares due to error,SSE)

均方误差
MSE=\frac{SSE}{n_r-k}
SSE=\sum\limits_{j=1}^{k}(n_j-1)s_j^2

我们注意到:MSE是以每个处理内部的变异性为依据,它不受原假设是否为真的影响。因此,MSE永远给出\sigma^2的一个无偏估计

回到例子:
SSE=\sum\limits_{j=1}^{k}(n_j-1)s_j^2=(5-1)\times27.5+(5-1)\times26.5+(5-1)\times31=340
MSE=\frac{SSE}{n_r-k}=\frac{340}{15-3}=28.33

13.2.3 方差估计量的比较:F检验

如果原假设H_0为真,则MSTR和MSE给出的\sigma^2的两个独立的无偏估计量。\sigma^2的两个独立的估计量纸币的抽样分布服从F分布。

k个总体均值相等的检验统计量:
F=\frac{MSTR}{MSE}
检验统计量服从分子自由度为k-1,分母自由度为n_r-k的F分布(ANOVA的假定要得到满足)

回到生产过滤系统的例子:在\alpha=0.05的显著水平下,进行假设实验,我们计算得到F=\frac{MSTR}{MSE}=\frac{260}{28.33}=9.18,分子自由度为2,分母自由度为12.

image

由于我们不希望F过大,所以根据EXCEL的计算,在上述自由度下,F=9.18时,上侧面积为0.0038<0.05,因此我们拒绝H_0,三个总体均值是不相等的。

当然也可以用临界值法,当\alpha=0.05时,F的临界值是3.8853<9.18。所以也拒绝H_0

总结

image

13.2.4 ANOVA表

前面的计算结果,可以使用方差分析表ANOVA表表示出来。一个完全随机化实验设计的ANOVA表的一般形式如下:

image

在上面的列表中,方差来源有SSTR和SSE,他们俩的总计被称作总平方和(SST),且SST=SSTR+SSE,并且自由度也是SSTR以及SSE的自由度的和。

总平方和SST的计算公式:SST=\sum\limits_{j=1}^{k}\sum\limits_{i=1}^{n_j}(x_{ij}-\bar {\bar x})^2
且:SST=SSTR+SSE
我们可以吧SST看作“处理平方和”与“误差平方和”的和。且自由度n_r-1也可由对应的SSTR和SSE的自由度加起来。
方差分析可以被看作将总平方和及其自由度分解成它们对应的来源(处理与误差)的一个过程。

13.2.5 方差分析的计算机输出结果

image

上图为MINITAB的计算结果,Pooled StDev是用来估计\sigma的,右下角的区间估计,三种方法的边际误差都是一样的这里的边际误差=t_{\alpha/2}\frac{s}{\sqrt{n}}=2.179\times \frac{5.323}{\sqrt{n}}=5.19。这里统一使用的误差平方和的自由度12(我也不知道为什么,我觉得好像应该用14,后面再研究)

13.2.6 k个总体均值相等的检验:一项观测性研究

例子:NCP公司对工厂员工的生产意识进行考试,共有3个工厂,每个工厂抽取6人。成绩如下:

image

假设:H_0:\mu_1=\mu_2=\mu_3 H_a:总体均值不全相等

总结

  1. 当每个样本都是n个观测值构成,则MSTR=\frac{\sum\limits_{j=1}^{k}n_j(\bar x_j-\bar {\bar x})^2}{k-1}=n\left[\frac{\sum\limits_{j=1}^{k}(\bar x_j-\bar {\bar x})^2}{k-1}\right]=ns_{\bar x}^2
  2. 当每个样本都是n个观测值构成,则MSE=\frac{\sum\limits_{j=1}^{k}(n_j-1)s_j^2}{n_r-k}=\frac{(n-1)\sum\limits_{j=1}^{k}s_j^2}{k(n-1)}=\frac{\sum\limits_{j=1}^{k}s_j^2}{k}

13.3 多重比较方法

方差分析只能告诉我们k个总体均值是否相等,但是具体哪些总体相等,哪些不相等,我们需要用多重比较方法在成对的总体均值之间进行统计比较。

13.3.1 Fisher的LSD方法

在方差分析钟拒绝了H_0,在这种情况下Fisher的最小显著性差异(least significant difference,LSD)方法可以用来确定哪些均值存在差异。

Fisher的LSD方法

H_0:\mu_i=\mu_j
H_a:\mu_i \neq \mu_j
检验统计量:
t=\frac{\bar x_i-\bar x_j}{\sqrt{MSE\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}}
拒绝法则:
p-值法:如果 p-值\leq \alpha,则拒绝H_0
临界值法:如果 t\leq -t_{\alpha/2}或者t\geq t_{\alpha/2},则拒绝H_0
其中t_{\alpha/2}是自由度为n_r-k时,t分布的上侧面积为\alpha/2的t值。
我们令\alpha=0.05,判断总体1(方法A)和总体2(方法B)的均值是否存在差异。
t=\frac{\bar x_i-\bar x_j}{\sqrt{MSE\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}}=\frac{62-66}{\sqrt{MSE\left(\frac{1}{5}+\frac{1}{5}\right)}}=-1.19

经过excel计算,t=-1.19,自由度为12时,的下侧面积为0.1285,双侧加起来即为p-值=0.2571>0.05所以,我们拒绝原假设,认为方法1和方法2的均值不相等。

基于检验统计量\bar x_i-\bar x_j的Fisher的LSD方法
H_0:\mu_i=\mu_j
H_a:\mu_i \neq \mu_j
检验统计量:
\bar x_i-\bar x_j
显著水平\alpha下的拒绝法则:如果|\bar x_i-\bar x_j|>LSD,则拒绝H_0
其中:LSD=t_{\alpha/2}\sqrt{MSE\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}

在过滤系统的例子中,通过计算得到LSD=t_{\alpha/2}\sqrt{MSE\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}=2.179\sqrt{28.33\left(\frac{1}{5}+\frac{1}{5}\right)}=7.34

计算后,我们可以把三个总体的样本均值计算出来,比如总体1和总体3的样本均值差为62-52=10>7.34,这就意味着我们拒绝认为总体1和总体3均值相等。

Fisher的LSD方法的两个总体均值之差的置信区间估计
\bar x_i-\bar x_j \pm LSD
LSD=t_{\alpha/2}\sqrt{MSE\left(\frac{1}{n_i}+\frac{1}{n_j}\right)}
其中t_{\alpha/2}是自由度为n_r-k时,t分布的上侧面积为\alpha/2的t值。

如果置信区间包含数值0,则不能拒绝两个总体均值相等的原假设。如果不包含则拒绝H_0

13.3.2 第一类错误

Fisher的LSD方法被称为保护性或限制性LSD检验,这是因为只有当我们首先找到一个用于方差分析的显著的F值时,才能使用LSD检验。
第Ⅰ类错误概率实验方式的第Ⅰ类错误概率

image

我们都是用\alpha=0.05的显著水平,对每个检验来说犯第Ⅰ类错误的概率为0.05,我们把这个概率称作比较方式的第Ⅰ类错误概率,表示单个的两两比较相联系的显著性水平。
在三次检验中至少有一次犯第Ⅰ类错误的概率为1-0.95^3=0.1426,我们称这个概率为实验方式的第Ⅰ类错误概率,记作\alpha_{EW}
当总体较多时,实验方式的第Ⅰ类错误概率就会比较大。

如何控制\alpha_{EW}呢?-使用Bonferrani修正方法
假设我们想要检验C个成对的两两比较(C=C_n^2
我们令\alpha=\frac{\alpha_{EW}}{C},例如针对5个总体,10种比较,想让实验方式的第Ⅰ类错误概率为0.05,则\alpha=\frac{\alpha_{EW}}{C}=\frac{0.05}{10}=0.005

但是一类错误和二类错误是成反比的,所以如何去权衡是个问题。也有其他方法,如Turkey方法、Duncan多重区域检验等,哪种更优有争议。

13.4 随机化区组设计

F=\frac{MSTR}{MSE}
有时外部因素(实验中没有考虑到的因素)引起MSE变大时,F将会变小。让我们误以为处理间没有差异,但是事实上是存在的。

本节将会介绍随机化区组设计(randomized block design)的实验设计。这个方法主要是通过消除MSE来自外部的变异,来达到控制变异外部来源的目的。

13.4.1 空中交通管理员工作压力测试

举例:探究不同工作系统是否产生不同的压力。现有3种设计方案,我们要探究不同方案之间有多大差异。
管理者希望管理员个人的变异性是MSE项的主要贡献者,将个人差异分离出来的一种办法是使用随机化区组设计。随机化区组需要管理员的一个单样本,分别在三个工作站接受检验。即工作站是影响因子,管理员是区组。(后面简称工作站为系统A、B和C)

每个个体都需要接受三次检验,检验顺序也需要是随机的。值是工作压力的度量。


image

截屏2020-12-27 10.19.49

13.4.2 ANOVA方法

随机化区组设计的ANOVA方法,要求我们将总平方和(SST)分解成:处理平方和(SSTR)、区组平方和(SSBL)和误差平方和(SSE)。
SST=SSTR+SSBL+SSE

image

  • k:处理的个数
  • b:区组的个数
  • n_r=kb:总样本容量

随机化区组设计,主要功能就是通过划分区组,将个人的差异从MSE中剔除。

13.4.3 计算与结论

  • x_{ij}代表区组i中对应处理j的观测值
  • \bar x_{.j}代表第j个处理的样本均值
  • \bar x_{i.}代表第i个区组的样本均值
  • \bar {\bar x}

步骤:

  1. 计算平方和SST=\sum\limits_{i=1}^{b}\sum\limits_{j=1}^{k}(x_{ij}-\bar {\bar x})^2
  2. 计算处理平方和SSTR=b\sum\limits_{j=1}^{k}(x_{.j}-\bar {\bar x})^2
  3. 计算区组平方和SSBL=k\sum\limits_{i=1}^{b}(x_{i.}-\bar {\bar x})^2
  4. 计算误差平方和SSE=SST-SSTR-SSBL

计算得到:

image

F=\frac{MSTR}{MSE}=\frac{10.5}{1.9}=5.53
分子、分母对应自由度为2和10,在F=5.53时的上侧面积即p-值为0.024<0.05则我们拒绝H_0:\mu_1=\mu_2=\mu_3

上述的例子是完全区组设计,即每个区组都要做k个处理。对应不完全区组设计,即某些(不是全部)处理被用于每个区组(如每个人都完成了系统A和B的检验,只有个别人完成了系统C的检验)

注释
由于有b个区组,使得自由度减少了b-1,所以随机化区组设计的误差自由度小雨完全随机化设计的误差自由度。如果n很小,因为误差自由度的减少,区组的潜在影响可能被掩盖;当n很大时,这种影响被最小化了。

13.5 析因实验

有时,我们需要得到一个以上变量或因子的统计结论。析因实验(factorial experiment)是一种实验设计。

举例:GMAT考试(商学院研究生考试),分数在200~800之间。现在有3种GMAT辅导课程。考生本科来自3种类型的院校。对应有9种处理组合,每个处理组合容量为2,意味着有两个复制
从种类型学校,每个学校取6人,分三组,随机分配到一个辅导课程。

image

我们希望得到的答案:

  • 主影响(因子A):辅导课程的不同是否对GMAT成绩有影响。
  • 主影响(因子B):本科院校的不同是否对GMAT成绩有影响。
  • 交互影响(因子A和B):辅导课程类型的影响是否取决于本科院校。

13.5.1 ANOVA方法

两因子析因实验的ANOVA方法要求我们将总平方和(SST)分为四个部分:因子A的平方和(SSA)、因子B的平方和(SSB)、交互作用的平方和(SSAB)、误差平方和(SSE)。
SST=SSA+SSB+SSAB+SSE

image

  • a代表因子A的水平数
  • b代表因子B的水平数
  • r代表复制的个数
  • n_r=abr观测值总数

15.5.2 计算与结论

  • x_{ijk}对应因子A处理i和因子B的处理j的第k次重复
  • \bar x_{i.}处理i(因子A)的观测值样本均值
  • \bar x_{.j}处理j(因子B)的观测值样本均值
  • \bar x_{ij}对应处理i(因子A)和处理j(因子B)的组合观测值样本均值
  • \bar {\bar x}n_r个观测值的总样本均值
  1. 计算总平方和SST=\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}\sum\limits_{k=1}^{r}(x_{ijk}-\bar {\bar x})^2
  2. 计算因子A的平方和SSA=br\sum\limits_{i=1}^{a}(\bar x_{i.}-\bar {\bar x})^2
  3. 计算因子B的平方和SSB=ar\sum\limits_{j=1}^{b}(\bar x_{.j}-\bar {\bar x})^2
  4. 计算交互作用的平方和SSAB=r\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}(\bar x_{ij}-\bar x_{i.}-\bar x_{.j}+\bar {\bar x})^2
  5. 计算误差平方和SSE=SST-SSA-SSB-SSAB

得到计算结果:


image

一般中型到大型的析因实验中涉及大量计算,需要用计算机。


image

综上,

  • 课程(因子A)对GMAT成绩影响,差异不显著。
  • 本科院校(因子B)对GMAT成绩影响,差异显著
  • 最后交互影响由于p-值>0.05,即不存在显著的交互作用的影响。

链接: https://pan.baidu.com/s/1fc0q-Q4kj3g-7Fr4MHZaqw 提取码: 333c 复制这段内容后打开百度网盘手机App,操作更方便哦

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,686评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,668评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 158,160评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,736评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,847评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,043评论 1 291
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,129评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,872评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,318评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,645评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,777评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,470评论 4 333
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,126评论 3 317
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,861评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,095评论 1 267
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,589评论 2 362
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,687评论 2 351

推荐阅读更多精彩内容