第1章 时空数据的例子

#翻译书籍(仅供学术交流)

#目录


#1.1 引文

      本章介绍了几个已被实际用于执行和说明时空建模的各个方面的激励数据集。这里的主要目的是让读者熟悉数据集,以便在引入建模概念时能够掌握这些概念。一个典型的读者,主要对应用建模感兴趣,而不是对理论发展感兴趣,能够快速浏览这里呈现的所有数据集,以便选择一个更接近他们的研究方向的数据集,用于后续针对特定数据集的建模发展。

     本章由三个主要部分组成。第1.2节讨论了三种广泛的时空数据类型。第1.3节提供了6个点引用数据的例子,第1.4节提供了5个面积单元数据集的例子。本书的不同章节使用了两个特定的数据集作为运行示例,一个是纽约州的空气污染监测数据集(第1.3节),另一个是英国2020年全球大流行期间Covid-19死亡人数数据集(第1.4节)。

      时空建模背后的思想可以根据:(a)它们的动机,(b)它们的潜在目标,和(c)数据的规模广泛地交叉分类。(a)模型的动力可分为四类:(i)时间序列方法向空间的扩展,(ii)随机场和成像技术向时间的扩展,(iii)时间和空间方法的相互作用,(iv)物理模型。(b)项下的主要目标可以看作是数据缩减或预测。最后,在(c)项下,可用数据可能在时间上是稀疏的,也可能在空间上是密集的,建模方法通常会考虑到这种数据规模。此外,数据可以在空间和/或时间上连续索引或离散索引。根据这些考虑,特别是第(一)-(三)项,建立统计模型并加以执行。

#1.2 时空数据类型

       只要每个数据都带有一个位置和一个时间戳,就称为时空数据。即使我们排除了两种极端退化的可能性,即数据要么在一个时间点,要么在一个空间位置被观察到,这也留下了大量时空数据类型的可能性。当分析师对空间或时间的变化不感兴趣时,他们通常会选择两种退化可能性中的一种。这可以简单地通过适当地聚合被忽略类别中的变异来实现。例如,当可以获得每日数据时,可以在一个网络中报告不同监测点的年均空气污染水平。在时间或空间上进行聚合会减少要分析的数据的可变性,也会限制可以做出的推论的范围。例如,仅仅通过分析年度总量来检测或分摊月度趋势是不可能的。本教材将假定数据中存在时空变化,尽管它将根据需要讨论研究空间或时间数据的重要概念。

      分析时空数据的首要任务之一是选择用于数据建模的时空分辨率。要考虑的主要问题是研究的推论目标。在这里,我们需要确定必须进行推断的最高空间和时间分辨率。例如,这些考虑将在很大程度上决定我们是否需要使用每日、月度或年度数据。这里还有其他的问题。例如,在较细的分辨率下,变量之间的关系可能比在较粗的分辨率下更容易理解,例如,小时空气污染水平与小时风速的关系可能比年平均空气污染水平与年平均风速的关系更强。因此,聚集可能导致相关性退化。然而,过于精细的时间和/或空间分辨率可能会给数据处理、建模和分析带来巨大的挑战,而不会添加太多额外的信息。因此,在决定要分析和建模的数据的时空分辨率时,需要保持一种平衡。这些决策必须在正式建模和分析开始之前的数据预处理阶段做出。

      假设数据的时间分辨率已经确定,我们用符号t表示每个时间点,我们假设总共有T个有规则间隔的时间点。这是一个限制性的假设,因为经常有不规则间隔的时间数据的情况。例如,一个病人可能在不同的时间间隔被跟踪;在定期检查期间,可能会有错过的约会。以空气污染监测为例,在空气污染严重时,可能需要进行更多的监测。例如,可能有必要决定是按小时还是按日的时间尺度模拟空气污染。在这种情况下,根据研究的主要目的,有两种主要的建模策略。

     第一种方法是在可能的最高规则时间分辨率下进行建模,在这种情况下,对于所有未观测到的时间点,将会有许多遗漏的观测结果。建模将有助于估计那些缺失的观测数据及其不确定性。例如,以每小时的分辨率模拟空气污染将导致所有未监测小时的观测数据缺失。第二种策略避免了第一种策略中许多缺失的观测,即在一个聚集的时间分辨率下建模,其中在一个规则的时间窗口内的所有可用观测都适当地平均,以得到将对应于该聚集时间的观测。例如,将所有每小时观测到的空气污染记录取平均值,得出每日空气污染平均值。尽管该策略成功地避免了许多遗漏的观测值,但在建模中仍可能存在挑战,因为聚合响应值可能显示不相等的方差,因为这些值是在聚合时间窗口期间从不同数量的观测值获得的。总之,在选择建模的时间分辨率时需要仔细考虑。

     为了对时空数据进行建模,显然需要跟踪在区域D中的空间位置s和观测的时间点t(即随机变量Y(s, t)的观测值y(s, t))。可能还有额外的协变量信息可用,我们将其表示为x(s, t)。不同的数据类型是由在D中观察点s的方式产生的。

     现在的注意力转向探索时间固定时数据集的空间方面。在空间分析中,有三种不同的数据类型:

    (i)点参考数据,

   (ii)区域单位数据,

   (iii)点模式数据。

      典型的点参考数据出现在每个数据点由地图上的一个非随机点描述的情况下,例如一个纬度-经度对。举个例子,在一个网络内的特定监测站点观察到了一个响应,例如空气污染水平值。对于点参考数据,它在固定的研究区域D内连续变化。举个例子,我们可能观察到了在一组n个位置si(i=1,...,n)和T个时间点t(t=1,...,T)处的响应y(s, t)和协变量。这组空间位置可以是固定的监测站,就像在空气污染的例子中一样,也可以随时间变化,例如从一艘研究船测量海洋特征时获得的数据,当它在海洋中移动时。

     区域单位数据通常在空间上聚合的领域中观察到,例如邮政编码、县或区等行政地理。对于这种数据,空间参考是地图中的一个区域。公共健康机构在疾病地图中广泛使用这种数据类型,例如Covid-19死亡率。第1.4节提供了几个例子。

      第三种数据类型是点模式,当感兴趣的事件(例如疾病的爆发)在随机位置发生时产生。因此,观察位置的随机性将点模式数据与点参考数据区分开来,尽管点模式中的位置由点来引用,每个点都可以由一个纬度-经度对来描述。通常,响应变量的概念在点模式数据集中通常是不存在的,尽管可能存在几个解释点的图案的解释变量。正如之前所指出的,本教科书不会进一步讨论或分析点模式数据。有兴趣的读者可以参考许多优秀的教科书,如Baddeley等人(2015年)和Diggle(2014年)。

#1.3 本书使用的点参考数据

##1.3.1 纽约市空气污染数据集

      我们使用了一个真实数据集,该数据集之前由Sahu和Bakar(2012a)进行了分析。这个数据集记录了2006年7月和8月每天最高8小时平均地面臭氧浓度,在纽约州的28个监测站点进行观测。我们考虑了三个重要的协变量:最高温度(xmaxtemp,以摄氏度为单位)、风速(xwdsp,以海里为单位)和相对湿度的百分比平均值(xrh),用于建立臭氧浓度的时空模型。关于协变量值及其用于预测的空间插值的更多详细信息可以在Bakar(2012)中找到。该数据集作为数据框nysptime包含在bmstdr包中。该数据集有12列和1736行。R帮助命令“?nysptime”提供了各列的详细描述。

     在这本书中,我们还将使用这个数据集的时间聚合空间版本。该数据集作为bmstdr包中的nyspatial数据集提供,它包含28个站点的数据,共有28行和9列。R语言中的帮助命令“?nyspatial”提供了各列的详细描述。该数据集中响应变量和协变量的值是nysptime中相应列的简单平均值(在删除缺失观测值之后)。图1.1展示了纽约州的地图以及28个监测位置。

```

library(bmstdr)

library(ggsn)

nymap <− map_data(database="state",regions="new york")

p <− ggplot() +

geom_polygon(data=nymap, aes(x=long, y=lat, group=group),color="black

", size = 0.6, fill=NA) +

geom_point(data =nyspatial, aes(x=Longitude,y=Latitude)) +

labs( title= "28 air pollution monitoring sites in New York", x="

Longitude", y = "Latitude") +

scalebar(data =nymap, dist = 100, location = "bottomleft", transform

=T, dist_unit = "km",

st.dist = .05, st.size = 5, height = .06, st.bottom=T,

model="WGS84") +

north(data=nymap, location="topleft", symbol=12)

p

```


图1.1    纽约的28个空气污染监测点

在第六章和第七章中,nyspatial和nysptime这两个数据集将作为所有点参考数据模型的示例。时空数据集nysptime还用于第九章中的预测说明。第三章对这些两个数据集进行了一些初步的探索性分析。

##1.3.2英格兰和威尔士的空气污染数据

        这个例子展示了从英格兰和威尔士的144个地点获得的氮氧化物(NO2)、臭氧(O3)和小于10微米(PM10)以及2.5微米(PM2.5)的颗粒物的日平均浓度的建模。这144个地点被分为三种站点类型:农村(16个)、城市(81个)和道路及路缘侧(RKS)(47个)。这些144个地点的位置如图1.2所示。

        Mukhopadhyay和Sahu(2018)对这一数据集进行了全面分析,并获得了英格兰和威尔士每个地方当局和单一当局的聚合年度预测。要重现Mukhopadhyay和Sahu(2018)所做的工作需要大量的计算和数据处理,因此我们在这里不尝试这样做。相反,我们展示了2011年一年中NO2的时空建模。我们还基于第8.2节中的每日数据建模获得了一个年度平均预测图。

图1.2  英格兰及威尔士地图,显示144个AURN空气污染监测点的位置。

##1.3.3 美国东部空气污染

     这个例子来自Sahu和Bakar(2012b),其中我们考虑对美国东部691个监测站点的每日最大8小时平均臭氧浓度数据进行建模,如图1.3所示。这些污染监测站点由646个城市和郊区监测站点组成,被称为国家空气监测站/州和地方空气监测站(NAMS/SLAMS),以及由清洁空气状况和趋势网络(CASTNET)监测的45个农村站点。

     我们从每年5月至9月的T=153天中分析每日数据,因为这是美国臭氧浓度最高的季节。我们考虑了从1997年到2006年的10年期间的数据,以便研究臭氧浓度水平的趋势。因此,我们总共有1,057,230个观测值,其中有大约10.44%是缺失的,我们假设这些缺失是随机的,尽管在这个缺失率中有一些年度变化。

      这里进行建模的主要目的是评估对主要臭氧标准的遵守情况,该标准规定年度第4高的最大每日8小时平均臭氧浓度水平的3年滚动平均值不应超过85 ppb,参见例如Sahu等人(2007)。图1.4绘制了第4高的最大值及其3年滚动平均值,并在85处叠加了一条水平线。如预期的那样,滚动平均值的绘图比年度第4高的最大值的绘图更平滑。这些绘图显示许多站点符合标准,但还有许多其他站点不符合。此外,滚动平均值的绘图显示出一个非常缓慢的下降趋势。这两个绘图都显示存在一些异常站点,这可能是由于空气污染的特定问题,例如由于森林火灾等自然灾害。这个数据集将在第8.3节进行分析。

##1.3.4哈伯德布鲁克降水资料

综合空间和时间的降水总量测量对于许多环境和生态原因,如空气和水的质量、洪水和干旱风险的时空趋势、林业管理和城市规划决策具有重要意义。


图1.3  美国东部691个臭氧监测点的地图。


图1.4   来自691个站点的臭氧浓度摘要的时间序列图。左面板:年度第4高的最大值,右面板:年度第4高的最大值的3年滚动平均值。y轴值为85的实黑色线是当时政府规定的臭氧标准。

       哈伯德布鲁克生态系统研究(HBES)位于美国新罕布什尔州,成立于1955年,持续观察许多环境结果变量,如温度、降水量、水流中的养分体积。HBES基于8000英亩的哈伯德布鲁克实验森林(见https://hubbardbrook.org/),是政策制定者、公众成员、学生和科学家的重要科学信息来源。这里感兴趣的是一个关于每周降水量的时空数据集,从1997年到2015年收集了来自22个雨量计的数据。

图1.5显示了23个雨量计(不包括RG22)的位置,它们被划分为南北面向集水区。RG1RG11量规位于地图北部面向南的流域,其余量规位于地图南部面向北的流域的一部分。南向的集水区由六个流域组成,标记为W1-W6,而北向的集水区由其余三个流域W7-W9组成。这里的主要模拟目标是研究流域降水量的时空趋势。我们没有复制该地图所需的数据。该数据集将在第8.4节中进行分析。


图1.5  雨量计和流域图。

##1.3.5 海洋叶绿素数据


图1.6  显示世界上23个重要海洋区域的地图

    这个例子摘自Hammond等人(2017),研究了海洋中叶绿素(chl)水平的长期趋势,chl是浮游植物(海藻)的代理指标。浮游植物处于食物链的最底层,是所有海洋生态系统的基础。浮游植物的丰度影响着营养物质的供应和光照。全球变暖可能会影响浮游植物的分布和丰度,因此研究影响浮游植物丰度的chl的长期变化趋势具有重要的科学意义。图1.6显示了我们通过卫星观测到的23个值得关注的海洋区域的地图。这里的主要模拟目标是研究这23个大洋区chl水平的长期趋势。第8.5节评估了这些趋势值。

##1.3.6  大西洋水温盐度数据集

   这个例子来自Sahu和Challenor(2008)关于从漫游的Argo浮标中建模深海温度数据的研究。Argo浮标项目,例如,可以查看http://www.argo.ucsd.edu,旨在全球范围内测量海洋上层2公里的温度和盐度。这些浮标记录实际测量值,与卫星数据形成对比,如第1.3.5节中的海洋叶绿素示例中使用的卫星数据,它们提供了许多缺失观测值的不太准确的观测结果。每个Argo浮标都被编程为下沉到一公里的深度,在该深度上漂浮约10天。在此之后,浮标进一步下沉一公里,达到两公里的深度,并通过调整其浮力上升到表面,在上升过程中测量温度和电导率(从中推导出盐度测量值)。一旦到达表面,数据和浮标的位置通过卫星传输。这使科学家能够获得接近实时的数据。传输数据后,浮标返回其“静止”深度一公里,并在另一个位置测量另一个温度和盐度剖面之前再漂浮10天。Argo数据可以通过国际Argo项目办公室免费获取,请参阅上述网站。

我们考虑北大西洋洋域在纬度20°至60°北和经度10°至50°西之间观察到的数据。图1.7显示了深海中Argo浮标的位置。该图显示了Argo浮标在每个月的运动性质。这里的主要建模目标是构建一个深海温度的年度地图以及其不确定性。观测数据的时间点不是等间隔的,我们在建模过程中不假设这一点。在第8.6节中进行了产生北大西洋海洋年度温度图所需的建模工作。


图1.7  2003年在深海中移动的Argo浮标位置。

#1.4  本书中使用的区域单元数据集

##1.4.1  英国Covid-19死亡率数据

      这个数据集展示了2020年3月13日至7月31日期间,英格兰的313个地方当局地区、县和统一管理局(LADCUA)中因Covid-19导致的死亡人数。请参考图1.8。在这段时间的20周内,每周记录的死亡人数为49,292人。图1.9显示了每个LADCUA的死亡人数和每10万人的死亡率的地图。对比这两个图表,可以明显看出右面板中的每10万人死亡率存在较大的空间变异。图1.10所示的每周死亡率箱线图显示了在第15周和第16周(4月10日至4月23日)出现的第一个峰值,以及峰值之后死亡人数非常缓慢地下降。这里的主要目的是对死亡率进行时空变异建模。该数据集将作为第10章中所有区域单元数据模型的示例使用。第3章提供了对该数据集的一些进一步的初步探索性分析。


图1.8   英国地方当局和九个行政区域的地图。空气污染监测点在地图上显示为蓝点。


图1.9  Covid死亡的原始人数(左)和每10万人中的Covid死亡人数(右)。

##1.4.2  肯尼亚儿童疫苗接种覆盖率

     人口和健康调查(DHS)计划定期收集多个数据集,用于全球范围内的健康监测。这个例子基于2014年东非国家肯尼亚的疫苗接种覆盖率数据集。该数据集包含在2014年调查之前任何时间接受过含麻疹疫苗(MCV1)的第一剂的12-23个月龄儿童的数量。图1.11绘制了2014年的观察到的疫苗接种比例。Utazi等人(2021)对这一数据和几个相关数据集进行了大量分析。在第11.2节中建模的这个例子旨在评估肯尼亚不同县的疫苗接种覆盖率。


图1.11   2014年观察到的疫苗接种比例。

##1.4.3  美国的癌症发病率

    美国疾病控制和预防中心提供不同地理级别(如50个州)的可下载癌症发病率数据。这样的数据集可以与各种信息一起下载,例如性别、种族和癌症类型。然而,由于数据可识别性和数据保护的原因,一些较小的费率计数(这是由于按因素进行更细的分类而产生的)费率没有公开。因此,为了说明本书的目的,我们的目标是在州一级建立汇总的年度数据模型。完整的数据集提供了2003年至2017年各州因各种原因导致的癌症年发病率。图1.12显示了从2003年到2017年,48个相邻州每10万人中各种原因的癌症总发病率。这是等值线图的一个例子,它使用颜色深浅或灰度将值分为几个大类,就像直方图一样。该数字表明,东北地区的总发病率高于西南地区。佛罗里达州也显示出较高的比率,这可能是由于该州退休老人的比例较大。完整的时空数据集将在11.3节中进行分析。


图1.12    2003-2017年因各种原因导致的每10万人癌症发病率的等值线图

      图1.13显示了10个选定州的观察到的标准化死亡率(见关于如何获得这些数据的第11.3节的讨论)。这些状态是手工挑选的,以表示SMR值的全部范围。这里有趣的研究问题是,“在考虑了时空相关性和任何其他重要的固定效应协变量后,这些比率是否有上升趋势?”第11.3节对此进行了研究。


图1.13  美国10个选定州的标准化癌症死亡率。

##1.4.4  英国住院数据

      根据Lee等人(2017)的研究,在2007年至2011年的5年期间,英格兰的323个地方和单位当局(LUA)每月报告的因呼吸系统疾病住院的人数是可获得的。这些计数取决于受风险人群的规模和人口统计学特征,通过计算预期的入院人数Eit来进行调整,使用所谓的间接标准化方法,参见第2.12节,从英格兰全国年龄和性别特定的住院率中得出。

在本例中,研究区域是英国的英格兰,分为i = 1,...,n = 323个地方和单位当局(LUA),数据可用于t = 1,...,T = 60个月,时间范围为2007年至2011年。对于LUA i和月份t的呼吸系统住院人数的计数由Yit表示,其中i = 1,...,323和t = 1,...,60,其中位数值为111,范围从6到2485。每月的时间尺度与Greven等人(2011)的研究相匹配,而大多数研究如Lee等人(2009)则使用年度数据。月度尺度的一个优点是它需要较少的数据聚合远离个体水平,但这也意味着Yit可能包括由慢性和急性污染暴露驱动的入院。

图1.14显示了标准化发病率比SMRit = Yit/Eit的空间(左面板)和时间(底部面板)模式,其中值1.2对应于相对于Eit增加20%的风险。该图显示了英格兰中部和北部城市的最高风险,如伯明翰、利兹和曼彻斯特,而时间模式具有强烈的季节性,冬季由于流感流行和低温等因素导致入院风险较高。这个数据集在本书籍的第11.4节中作为示例使用。



图1.14   顶部面板:平均SMR的地图。底部面板:5年内每月SMR值的箱线图

##1.4.5  伦敦的儿童贫困状况

      该数据集关注的是英国32个行政区和伦敦市每年儿童贫困水平的监测。我们所考虑的儿童贫困衡量标准是2010年英国儿童贫困法案中规定的相对低收入儿童贫困的广泛代表。这里的目的是在本地级别执行分析。该指标显示了家庭中领取失业救济金(经济情况调查)或税收抵免的儿童所占的比例,这些儿童报告的收入低于英国收入中值的60%。这个例子将分析从2006年到2015年10年的儿童贫困数据。图1.15显示了平均水平,图1.16显示了贫困水平的下降趋势。将对章节11.5中的数据集进行时空建模。


图1.15   伦敦市和伦敦32个行政区的地图。1 =伦敦市,2 =陶尔哈姆莱茨,3 =伊斯灵顿,4 =卡姆登,5 =威斯敏斯特,6 =肯辛顿和切尔西,7 =汉默史密斯和富勒姆。


图1.16  伦敦贫困儿童比例

#1.5 结论

      本章概述了几个数据示例,用于说明点参考和区域单元数据的时空模型拟合和验证。每个示例的目标是向读者介绍这些数据集和相关的建模问题,以激发他们研究后续章节中介绍的方法。同时,鼓励有实际需求的读者仅阅读与他们自己的实际建模问题相似的示例。除了纽约和肯尼亚的示例外,本章只提供关于如何绘制地图的评论。为了简洁起见,它没有给出实际的代码来重现图形。这些代码行可以在线从github(https://github.com/sujit-sahu/bookbmstdr)获取,以便读者可以直接在R中进行代码实验,而不是在这里阅读那些繁琐的文本块。然而,所提供的评论涉及到重现图形所需的主要步骤和技巧。这些技巧在其他地方也会很有用。

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

推荐阅读更多精彩内容