应用统计学与R语言实现学习笔记(九)——线性回归

Chapter 9 Linear Regression

本篇是第九章,内容是回归分析(主要以线性回归为主)。回归分析是数理统计、数理分析中最基础(也可以说是最重要)的一个分析,所以这一章内容相对来说也较多。

1 变量间的关系

  • 确定型关系vs不确定型关系
    函数关系——一一对应的确定型关系设有两个变量x和y,变量y随变量x一起变化, 并完全依赖于x,当变量x取某个数值时,y依确定的关系取相应的值,则称y是x的函数,记为y=f(x),其中x称为自变量,y称为因变量各观测点落在一条线上。
    相关关系(correlation)——变量间关系不能用函数关系精确表达。一个变量的取值不能由另一个变量唯一确定。当变量x取某个值时, 变量y的取值可能有几个。各观测点分布在直线周围。

相关关系包括了线性相关(正相关、负相关)、非线性相关、完全相关(正相关、负相关)、不相关。

除了如上的图,可以看下面的链接——关于相同统计量不同数据的一篇外文。

https://www.autodeskresearch.com/publications/samestats

相关系数(correlation coefficient)

  • 对变量之间关系密切程度的度量(只关心密切程度,无关因果关系);
  • 对两个变量之间线性相关程度的度量称为简单相关系数;
  • 若相关系数是根据总体全部数据计算的,称为总体相关系数,记为ρ;
  • 若是根据样本数据计算的,则称为样本相关系数,记为 r。

总体相关系数的计算公式

相关系数特点

  • 无量纲(Unitfree);
  • ρ的取值范围是 [-1,1];
  • |ρ|=1,为完全相关(ρ=1为完全正相关;ρ=-1为完全负相关);
  • ρ=0,不存在线性相关关系;
  • -1≤ρ<0,为负相关,0<ρ≤1,为正相关;
  • |ρ|越趋于1表示线性关系越密切;|ρ|越趋于0表示线性关系越不密切;
  • 若X与Y相互独立,则ρ=0,但ρ=0,X与Y不一定相互独立;
  • 若ρ= 0,且X与Y服从正态分布,则X与Y相互独立。

样本相关系数计算公式

样本相关系数特点

  • 无量纲(Unitfree);
  • r的取值范围是 [-1,1];
  • |r|=1,为完全相关(r=1为完全正相关;r=-1为完全负相关);
  • r=0,不存在线性相关关系;
  • -1≤r<0为负相关,0<r≤1为正相关;
  • |r|越趋于1表示线性关系越密切;|r|越趋于0表示线性关系越不密切;

对变量之间关系密切程度的度量,只关心密切程度,无关因果关系。
比如撑伞的人数和降雨量的相关系数非常高。但是我们不能说因为撑伞的人多了,所以降雨量大。

r的抽样分布
r的抽样分布随总体相关系数和样本容量的大小而变化。当样本数据来自服从正态分布的总体时,随着n的增大,r的抽样分布趋于正态分布,尤其是在总体相关系数ρ很小或接近0时,趋于正态分布的趋势非常明显。而当ρ远离0时,除非n非常大,否则r的抽样分布呈现一定的偏态。当ρ为较大的正值时, r呈现左偏分布;当ρ为较小的负值时, r 呈现右偏分布。只有当ρ接近于0,而样本容量n很大时,才能认为r是接近于正态分布的随机变量。
相关系数的显著性检验步骤
检验两个变量之间是否存在线性相关关系,等价于对回归系数β1的检验。采用R. A. Fisher提出的t检验。
检验的步骤为:

2 回归分析和简单线性回归分析

2.1 回归分析

什么是回归分析(Regression)?

从一组样本数据出发,确定变量之间的数学关系式。对这些关系式的可信程度进行各种统计检验,并从影响某一特定变量的诸多变量中找出哪些变量的影响显著, 哪些不显著。利用所求的关系式,根据一个或几个变量的取值来预测或控制另一个特定变量的取值, 并给出这种预测或控制的精确程度。

回归分析与相关分析的区别

相关分析中,变量x变量y处于平等的地位;回归分析中,变量y称为因变量,处在被解释的地位,x称为自变量,用于预测因变量的变化;
相关分析中所涉及的变量x和y都是随机变量;回归分析中,因变量y是随机变量,自变量x可以是随机变量,也可以是非随机的确定变量;
相关分析主要是描述两个变量之间线性关系的密切程度;回归分析不仅可以揭示变量x对变量y的影响大小,还可以由回归方程进行预测和控制。

回归模型(regression model)——回答“变量之间是什么样的关系?”方程中运用1个数值型因变量(响应变量)作为被预测的变量;1个或多个数值型或分类型自变量 (解释变量)作为用于预测的变量。主要用于预测和估计。回归模型的类型包括一元回归模型(线性和非线性)和多元回归模型(线性和非线性)。
接下来先从简单线性回归分析讲起。

2.2 简单线性回归分析

简单线性回归(Simple Linear Regression)——涉及一个自变量的回归,因变量y与自变量x之间为线性关系。被预测或被解释的变量称为因变量(dependent variable),用y表示;用来预测或用来解释因变量的一个或多个变量称为自变量(independent variable),用x表示。因变量与自变量之间的关系用一个线性方程来表示。
描述因变量y如何依赖于自变量x和误差项ε的方程称为回归模型(Regression Model,定义如前)。
(1)简单线性回归模型的表示形式

y是x的线性函数(部分)加上误差项(residual/random error term)。线性部分反映了由于x的变化而引起的y的变化。误差项ε是随机变量。反映了除x和y之间的线性关系之外的随机因素对y的影响,是不能由x和y之间的线性关系所解释的变异性。β0和β1称为模型的参数(interception, slope)。
(2)简单线性回归模型的基本假定
误差项ε是一个期望值为0的随机变量,即E(ε)=0。对于一个给定的x值,y的期望值为

(3)简单线性回归方程(regression equation)
描述y的平均值或期望值如何依赖于x的方程称为回归方程;简单线性回归方程的形式如下

方程的图示是一条直线,也称为直线回归方程。β0是回归直线在y轴上的截距(interception),是当x=0时y的期望值。β1是直线的斜率(slope),称为回归系数,表示当x每变动一个单位时,y的平均变动值。
(4)估计的回归方程(estimated regression equation)

(5)最小二乘估计

在r语言中,简单线性回归的代码如下:

modele<-lm(e~a)

(7)回归直线的拟合优度
变差
因变量 y 的取值是不同的, y 取值的这种波动称为变差。 变差来源于两个方面:

离差平方和的分解(三个平方和的关系与意义)

从左至右分别为SST,SSR,SSE。
所以就有SST=SSR+SSE。
总平方和(SST)——反映因变量的 n 个观察值与其均值的总离差;
回归平方和(SSR)——反映自变量 x 的变化对因变量 y 取值变化的影响,或者说,是由于x与y之间的线性关系引起的y的取值变化,也称为可解释的平方和;
残差平方和(SSE)——反映除x以外的其他因素对y取值的影响,也称为不可解释的平方和或剩余平方和。

判定系数R²(coefficient of determination)
回归平方和占总离差平方和的比例。

估计标准误差(standard error of estimate)

显著性检验

  • 线性关系的显著性检验:检验自变量与因变量之间的线性关系是否显著,即检验x与y之间是否具有线性关系,或者说,检验自变量x对因变量y的影响是否显著;
  • 回归系数的显著性检验:检验回归系数是否不等于0;
  • 在简单线性回归中,线性关系的显著性检验等价于回归系数的显著性检验。
    线性关系的检验
    将回归均方(MSR)同残差均方(MSE)加以比较, 应用F检验来分析二者之间的差别是否显著。
    回归均方:回归平方和SSR除以相应的自由度(自变量的个数p);
    残差均方:残差平方和SSE除以相应的自由度(n-p-1)。

回归系数的检验(检验步骤)

显著性检验的几点注意
显著性关系的结论不意味着因果关系。显著性关系的结论也不能推出线性关系的结论,仅能说在x的样本观测之范围内,x和y是相关的,而且一个线性关系只揭示了y的变异的主要部分。当样本容量很大时,对于小的b1值也能得到统计上是显著的结果。

3 利用回归方程进行估计和预测

根据自变量x的取值估计或预测因变量y的取值。
估计或预测的类型

(1)点估计:y的平均值的点估计,y的个别值的点估计;
(2)区间估计:y的平均值的置信区间估计,y的个别值的预测区间估计。

(1)点估计
对于自变量x的一个给定值x0,根据回归方程得到因变量y的一个估计值^y0。
点估计值有y的平均值的点估计y的个别值的点估计。在点估计条件下,平均值的点估计和个别值的的点估计是一样的,但在区间估计中则不同。

(2)区间估计
点估计不能给出估计的精度, 点估计值与实际值之间是有误差的, 因此需要进行区间估计。对于自变量x的一个给定值$x_0$,根据回归方程得到因变量y的一个估计区间。区间估计有两种类型:置信区间估计(confidence interval estimate)预测区间估计(prediction interval estimate)

影响区间宽度的因素

其实在R语言里主要用predict.lm函数来进行区间估计。代码样例如下:

con<-predict.lm(modele,h,interval="confidence",level=0.95)

其中interval控制是置信区间(参数填confidence)、预测区间(参数填prediction)或者是不做区间估计,level是置信水平,接着用R绘制一个简单的回归和置信区间的图,这里先给出如何绘制置信区间band的代码,完整代码还是老规矩,在这一部分笔记写完后给出。

polygon(c(h[,1], rev(h[,1])), c(con[,3], rev(con[,2])),border="red",lwd=1,lty = c("dashed", "solid"))

4 残差分析

残差(residual)——因变量的观测值与根据估计的回归方程求出的预测值之差,用e表示。

反映了用估计的回归方程去预测而引起的误差。
残差检验的目的

  • 检验线性的假设是否成立;
  • 确定有关误差项ε的假定是否成立(正态分布;方差为常数;独立性)。
  • 检测有影响的观测值。

残差图(residual plot)

  • 表示残差的图形(关于x的残差图,关于y的残差图,标准化残差图)。
  • 用直方图或正态概率图检验正态性。

标准化残差(standardized residual)

标准化残差图
用以直观地判断误差项服从正态分布这一假定是否成立。

  • 若假定成立, 标准化残差的分布也应服从正态分布。
  • 在标准化残差图中, 大约有95%的标准化残差在-2到+2之间。

变换
数据变换的问题在前面第七章拟合优度检验提过,那么什么时候做变换?
如果从散点图观察发现残差是自变量的函数,通过变换可能可以解决问题。
做什么变换?观察残差与因变量观测值的均值的关系:

  • 如果残差的标准差与因变量观测值的均值有线性关系,用log变换;
  • 如果残差的方差与因变量观测值的均值有线性关系,用square root变换;
  • 如果残差的标准差与因变量观测值的均值的平方有线性关系,用inverse变换;
  • 如果残差的标准差与因变量观测值的均值的幂有线性关系,用power变换。

序列相关(自相关)
当数据是按时间顺序采集的,有可能引起误差项之间的相关(Serial correlation,autocorrelation)。
这里介绍一个相关的杜宾-瓦特森(Durbin-Watson)检验统计量:

是否遗漏了重要的对因变量有时序影响的自变量,有时可通过引入度量观测次数的自变量解决该问题。这部分属于时间序列分析的范畴,这里就不进一步阐述了。

在R语言中,线性回归方程残差图绘制非常简单。模型拟合过程会自动给出四个残差可视化相关的图。绘制方法如下:

layout(matrix(c(1,2,3,4),nrow=2,byrow=T))
plot(modele)

结果如图

异常值(outlier)与识别
如果某一个点与其他点所呈现的趋势不相吻合,这个点就有可能是异常点。

  • 如果异常值是一个错误的数据, 比如记录错误造成的, 应该修正该数据, 以便改善回归的效果;
  • 如果是由于模型的假定不合理, 使得标准化残差偏大, 应该考虑采用其他形式的模型,比如非线性模型;
  • 如果完全是由于随机因素而造成的异常值, 则应该保留该数据。

在处理异常值时, 若一个异常值是一个有效的观测值, 不应轻易地将其从数据集中予以剔除。

  • 异常值也可以通过标准化残差来识别;
  • 如果某一个观测值所对应的标准化残差较大, 就可以识别为异常值;
  • 一般情况下,当一个观测值所对应的标准化残差小于-2或大于+2时,就可以将其视为异常值。

有影响的观测值
如果某一个或某一些观测值对回归的结果有强烈的影响,那么该观测值或这些观测值就是有影响的观测值。
一个有影响的观测值可能是:一个异常值, 即有一个值远远偏离了散点图中的趋势线;对应一个远离自变量平均值的观测值;或者是这二者组合而形成的观测值。
如果有影响的观测值是一个错误的数据,比如记录错误造成的, 应该修正该数据,以便改善回归的效果。
如果有影响的观测值是一个有效的数据则应该保留它, 可以帮助我们分析模型的假定是否合理。
杠杆率点(leverage point)
如果自变量存在一个极端值, 该观测值则称为高杠杆率点(high leverage point),在简单回归中,第i个观测值的杠杆率用$h_i$表示,其计算公式为:

如果一个观测值的杠杆率hi>n/6,就可以将该观测值识别为有高杠杆率的点;
一个有高杠杆率的观测值未必是一个有影响的观测值, 它可能对回归直线的斜率没有什么影响。

5 多元线性回归(multiple regression model)

多元线性回归(multiple regression model)

多元回归模型的基本假定

多元回归方程(multiple regression equation)

二元回归方程的几何表达——回归面。

估计的多元回归的方程(estimated multiple regression equation)

参数的最小二乘法

多重判定系数(multiple coefficient of determination)
回归平方和占总平方和的比例,计算公式为

因变量取值的变差中, 能被估计的多元回归方程所解释的比例。
修正多重判定系数(adjusted multiple coefficient of determination)

估计标准误差s
对误差项ε的标准差σ的一个估计值。衡量多元回归方程的拟合优度。计算公式为

这里写图片描述

线性关系检验
检验因变量与所有自变量之间的线性关系是否显著,也被称为总体的显著性检验。检验方法是将回归均方和(MSR)同离差均方和(MSE)加以比较,应用F检验来分析二者之间的差别是否显著。

  • 如果是显著的, 因变量与自变量之间存在线性关系;
  • 如果不显著, 因变量与自变量之间不存在线性关系。

回归系数的检验(检验步骤)

  • 线性关系检验通过后,对各个回归系数进行检验。
  • 对每一个自变量单独应用 t 检验统计量进行检验。

回归系数的推断(置信区间)
回归系数在(1-α)%置信水平下的置信区间为

回归系数的抽样标准差

6 多重共线性(multicollinearity)

回归模型中两个或两个以上的自变量彼此相关。多重共线性带来的问题有:可能会使回归的结果造成混乱, 甚至会把分析引入歧途;可能对参数估计值的正负号产生影响, 特别是各回归系数的正负号有可能同我们预期的正负号相反。
多重共线性的识别

  • 检测多重共线性的最简单的一种办法是计算模型中各对自变量之间的相关系数, 并对各相关系数进行显著性检验;
    若有一个或多个相关系数显著, 就表示模型中所用的自变量之间相关,存在着多重共线性。
  • 如果出现下列情况,暗示存在多重共线性:
    模型中各对自变量之间显著相关。
    当模型的线性关系检验(F检验)显著时,几乎所有回归系数的t检验却不显著。
    回归系数的正负号与预期的相反。

检测多重共线性(Variance Inflationary Factor)

多重共线性(问题的处理)
将一个或多个相关的自变量从模型中剔除,使保留的自变量尽可能不相关。
如果要在模型中保留所有的自变量,则应避免根据t统计量对单个参数进行检验,对因变量值的推断(估计或预测)的限定在自变量样本值的范围内。

7 定性自变量的回归

虚拟变量(dummy variable)
定性自变量————只有两个水平的定性自变量或有两个以上水平的定性自变量。
虚拟变量——用数字代码表示的定性自变量。
虚拟变量的取值为0,1。
虚拟变量的个数

8 非线性回归

(1)二阶回归模型(Quadratic Regression Model)——当散点图如下所示,可考虑二次回归模型。

二阶回归模型的显著性检验

(2)交互作用
交互作用——两个自变量共同作用对因变量产生的潜在影响。

交互作用显著性检验

(3)其他非线性回归
因变量y与x之间不是线性关系,可通过变量代换转换成线性关系,用最小二乘法求出参数的估计值。但是并非所有的非线性模型都可以化为线性模型。

  • 双曲线
  • 幂函数曲线
  • 对数曲线
  • 指数曲线
  • S型曲线

9 建立回归模型

得到描述因变量与一个或一个以上自变量之间关系的估计的回归方程。目的是建立一个基于最好自变量集合的模型。找到一个适合的描述变量关系之间关系的函数。选择模型应包含的变量。

  • 俭约的模型–用尽可能少的变量来提供足够精度的预测。
  • 将不重要的变量除去更容易对模型进行解释。
  • 发生多重共线性的可能变小。

变量选择Variable Selection

有些变量的作用不是很大,SSE 不会随着变量个数的增加而增加,但MSE=SSE/(n-k-1) 有可能会随着变量
个数的增加而增加。最小的MSE可作为最优变量选择的一个准则,但需考虑所有子集 (2^p个)。

检验增加变量是否适宜的F统计

F越大,说明增加变量减少预测误差的效果越显著。
变量选择过程

  • 向前选择(Forward Selection)
  1. 从没有自变量的模型开始。
  2. 如果所有的F统计量的p-值大于预先设定的终止值,说明增加任一变量效果不显著,停止。
  3. 否则,加入具有最大F统计量值的变量。
  4. 重新回归, Go to Step 2。
  • 后向消元(Backward Elimination)
  1. 从包含所有自变量的模型开始。
  2. 如果所有的F统计量的p-值小于预先设定的终止值,说明减少任一变量效果显著,停止。
  3. 否则,删除具有最小F统计量值的变量。
  4. 重新回归, Go to Step 2。
  • 逐步回归(Stepwise regression procedure)
    向前选择和后向消元的结合。
    1.先检查是否有变量需从模型中删除。
    2.再检查增加一个变量是否能改善模型。
    3.重复以上过程。
    注意: α进≤α出,否则F进<F<F出,会导致无限循环。
  • 最佳子集回归(Best-subset approach)
    对所有可能的自变量组合进行估计。找出具有最大的修正判定系数$adj.R^2$和最小的估计误差标准差$s_ε$。

10 回归中的常见错误

(1)没有检验线性关系假设

画散点图。
如果不是线性的,检验其它非线性。
用线性关系描述非线性关系会引起误导。

(2)只看结果不看图表

要将画散点图作为回归分析的一部分。
检验回归直线与实际观测值间的关系。
对自动回归来说这一步更为重要。

(3)用回归系数判定变量的重要性

回归系数依赖于自变量的量纲,因此系数的大小与变量的重要性无关。
例如,将秒变为微秒没有改变任何事实,但是变量的系数却有所改变。

(4)没有确定置信区间

观察值是随机样本,所以回归结果有一定随机性。
不确定置信区间,不可能理解参数的真正含义。

(5)没有计算判定系数

没有$R^2$,很难确定多少变异是由回归解释的。
即使$R^2$看起来很好,安全起见还应做F-test。

(6)错误解释相关系数

判定系数是R²。
相关系数是R。
R²给出变异由回归解释的百分比,不是R。
如:R =0.5,R²=0.25——回归解释了25%的变异,不是50%。

(7)使用强相关的自变量

模型同时包括两强相关的自变量会降低回归模型的显著性。
要尽可能的了解自变量间的关系。

(8)用回归模型预测观测值范围之外的区域

回归是基于某一特定观测样本的。
在样本观测值范围内能提供较为精确的估计。

(9)观测值取值范围太小

回归只有在观测值取值范围附近预测的结果比较好。
如果不在常用的范围内取值,回归模型用处不大。

(10)包括太多的自变量

变量越多的模型不一定越好。
有可能出现多重共线性。

(11)认为好的预测变量是好的控制变量
相关关系不一定因果关系:A与B相关,并不意味着可以通过改变A来控制B。

(12)线性回归结果会给人以误导

为了提供一个简练的总结,回归过程中舍弃了一些信息。
有时一些重要的特征也舍弃了——看图形表示可以告诉我们是否有问题。

11 Logistic 回归

Logistic回归提出的目的是为了解决二值化数据的回归问题。那么为什么简单线性回归模型不适合二值化数据的回归呢?详细原因可见如下图。

二值化变量是“yes”或者"no"的数据。可以被编码为1和0,也就是说不会有其他的变异数值。所以对于这种情况模型的要求是:模型的边界为0和1,模型可以输出的是一个在这类或者另一类的概率。我们想要的是一个实际值落入这类或者另一类的概率大小。而理想的模型是很好的估计0和1,或者换句话说,结果是0或1。所以解决方案就是Logistic回归。

Logistic的基本形式为

典型案例:
城市增长问题,城市化预测模拟,

常见的问题

  • 都有一个二值化(或分类)变量:
  • 都涉及到预测的思想机会,概率,比例或百分比。
  • 不像其他的预测情况,y值是有界的。

Logistic 回归与简单线性回归

logistic回归是一种统计技术,可以用二值化变量问题中。回归虽有相似之处,但它不同于普通最小二乘法。识别重要和相似之处是两种技术的区别。

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

推荐阅读更多精彩内容