999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于傾向得分匹配與加權調整的非概率抽樣統計推斷方法研究

2016-12-20 03:30:54金勇進
統計與決策 2016年21期
關鍵詞:方法

劉 展,金勇進

(中國人民大學a.統計學院;b.應用統計科學研究中心,北京100872)

基于傾向得分匹配與加權調整的非概率抽樣統計推斷方法研究

劉 展a,金勇進b

(中國人民大學a.統計學院;b.應用統計科學研究中心,北京100872)

文章針對非概率抽樣統計推斷問題,提出了一種解決方法:首先采用傾向得分匹配選擇樣本,然后采用傾向得分逆加權、加權組調整和事后分層調整三種方法對匹配樣本進行加權調整來估計目標總體,并比較不同方法估計的效果。蒙特卡羅模擬與實證研究表明:當網絡訪問固定樣本大小與目標樣本大小的比率小于3時,三種加權方法估計的效果均比未加權時匹配樣本的估計效果好;當網絡訪問固定樣本大小與目標樣本大小的比率不小于3時,傾向得分事后分層調整與未加權的匹配樣本估計效果較好。

傾向得分;匹配;加權調整;非概率

0 引言

在當今的大數據時代,大數據還不能代替總體,但大數據下的抽樣仍然是必要的。只是由于數據的大體量、非結構,且數據來源復雜化,難以構造抽樣框,使得一些抽取的樣本屬于非概率樣本,難以將傳統的抽樣推斷理論應用到非概率樣本中。與此同時,隨著網絡調查特別是網絡訪問固定樣本的發展,非概率抽樣重新引起了人們廣泛的關注與重視,如何解決非概率抽樣的統計推斷問題,不僅是大數據背景下抽樣調查面臨的嚴重挑戰,也是網絡調查發展的迫切需求。

解決非概率抽樣統計推斷問題的一種方法就是樣本匹配。樣本匹配多年來一直被用于觀察性研究中,主要目的是根據一個或多個協變量找到與處理組相匹配(近似)的對照組,從而減少處理效應估計的偏差,直到最近才被提倡用于網絡訪問固定樣本的相關調查中[1]。所謂網絡訪問固定樣本[2]就是愿意完成網絡調查的網絡訪問(上網)人群,這就意味著存在一個潛在的受訪者的樣本數據庫,在未來的數據收集中,如果他們被選擇為調查對象,他們將愿意配合完成調查。目前已有研究者對網絡訪問固定樣本調查中的樣本匹配問題進行了一些研究。Rivers[3]在2006年就提出可使用樣本匹配從網絡訪問固定樣本中抽取代表性樣本。Vavreck與Rivers[4]在2008年采用了這種方法,他們首先從公布的美國社區調查文檔中抽取了一個38000人的隨機樣本,并對每個隨機樣本單元從網絡訪問固定樣本中找到了最近的匹配單元,并利用匹配單元的調查數據來估計總體。Rivers與Bailey(2009)[5]使用了2008年美國總統選舉的數據評價了利用從網絡訪問固定樣本中選擇的匹配樣本進行推斷的效果,提出由于不完美的匹配,有必要在匹配后對匹配樣本進一步加權調整。Terhanian與Bremer(2012)[6]基于平行調查(一個隨機數字撥號電話調查和一個網絡訪問固定樣本調查)使用傾向得分來選擇匹配樣本,相對于直接基于協變量的匹配來說,傾向得分匹配具有將匹配的維度降為一維的優點,這極大的簡化了匹配的過程,因而受到廣泛的使用。事實上,傾向得分已經被用于非概率抽樣的加權調整階段(Terhanian et al.,2001;Lee&Valliant,2009),直到最近才被引入非概率調查的匹配抽樣之中(Rivers,2007;Terhanian&Bremer, 2012)。基于此,本文考慮在非概率樣本選擇階段采用傾向得分來選擇匹配樣本,同時為了減少不完美匹配造成的偏差,進一步采用傾向得分對匹配樣本進行加權調整,從而降低估計的偏差,提高估計的精度。

本文針對非概率抽樣的統計推斷問題,提出在樣本選擇階段采用傾向得分匹配選擇匹配樣本,進一步對匹配后的樣本采用傾向得分逆加權、加權組調整和事后分層調整三種方法進行加權調整,估計目標總體,并比較各種方法估計的效果,以期豐富非概率抽樣的統計推斷問題研究,促進非概率抽樣的廣泛應用,也給大數據時代下海量數據的處理與分析提供一定的啟示。

1 非概率抽樣的統計推斷方法

對于非概率抽樣的統計推斷,可分為兩個階段來考慮,第一個階段是非概率樣本的選擇階段,可考慮采用傾向得分匹配來選擇樣本;第二個階段是加權調整階段,主要是對第一階段選擇的匹配樣本進行傾向得分加權調整,最終實現對目標總體的估計。

1.1 基于傾向得分匹配的樣本選擇

傾向得分為在給定協變量Xi的條件下,個體i接受處理的條件概率[7]。假設是否接受處理為Di(接受處理,Di=1;否則Di=0),則第i個單元的傾向得分定義為:p(Xi)=P(Di=1|Xi)。根據傾向得分來進行匹配的方法就是傾向得分匹配(Propensity Score Matching,PSM)。這里,所謂的匹配就是假設個體i屬于處理組,找到屬于對照組的某個體 j,使得個體 j與個體i的協變量或傾向得分取值盡可能相似,即Xi?Xj或p(Xi)?p(Xj)。為了刻畫個體j與個體i之間的相似(相近)程度,產生不同的距離函數定義方法,也由此產生不同的匹配方法,常用的匹配方法有:最近鄰匹配、卡鉗與半徑匹配、分層或區間匹配、核與局部線性匹配等。

基于傾向得分匹配樣本選擇的基本思想是首先從抽樣框(包含一系列的協變量)中抽取一個概率樣本,其中目標總體單元個數為N,樣本量為n,每個單元i都有一個目標變量Yi和一些協變量組成的向量 Xi=(Xi1,Xi2,...,Xip),i=1…n,p為協變量的個數,稱這個樣本為目標樣本,因為這是一個概率樣本,因此具有總體的代表性。為了便于討論,假設從抽樣框中抽取的概率樣本為簡單隨機樣本。本來將調查目標樣本中的對象,但是當發現目標樣本的調查比較困難時,并不要求用目標樣本中的對象完成調查,相反地尋找網絡訪問固定樣本(包含了與抽樣框相同的一系列協變量)中與目標樣本對象相似(近)的成員,稱為匹配單元,匹配單元的集合稱為匹配樣本,邀請匹配樣本單元完成調查。假設匹配樣本的回答率為100%,可以看到匹配樣本本質上屬于非概率的樣本。這里,在尋找網絡訪問固定樣本中與目標樣本對象相似(近)的成員,即選擇匹配樣本時,采用傾向得分匹配方法。引入二值示性變量Di,如果單元i在目標樣本中,則Di=1(相當于接受處理),否則Di=0,在本文中Di=0表示單元i在網絡訪問固定樣本中。傾向得分匹配方法具體步驟如下:

(1)估計傾向得分

傾向得分常常需要估計,在估計 p(Xi)=P(Di=1|Xi)時,可使用參數估計(probit或logit)或非參數估計,最流行的方法為logit,即常常將示性變量(D)作為因變量,單元的協向量 X作為解釋變量建立Logistic回歸模型來得到。具體地,假設Xi都經過中心化變換,則有:

由式(1)可得:

(2)匹配樣本的選擇

當選擇匹配樣本時,對目標樣本的每個成員,有必要在網絡訪問固定樣本中找到最近的匹配單元,可以采用一些匹配方法,方法之一就是最近鄰匹配。最近鄰匹配(Nearest Neighbor Matching,NNM)[8]包括單一無放回最近鄰匹配、單一有放回最近鄰匹配和多重最近鄰匹配,本文主要采取單一無放回最近鄰匹配。設集合

為目標樣本(D=1)中每個單元i的一個鄰域,|| ||是一個范數,如1-范數、2-范數、¥-范數等。單一無放回最近鄰匹配就是將與Xi最近的Xj對應的一個網絡訪問固定樣本單元(Dj=0)選擇為匹配單元,且該匹配單元僅能匹配一個目標樣本單元。當然匹配方法不一定在所有的情形下都是非常有效的,當它滿足以下兩個條件[9]時最可能有效:

假設1:可忽略性。假設網絡訪問固定樣本關于用于匹配的變量(簡稱匹配變量)是可忽略的。這意味著如果本文檢驗匹配變量值完全相同的網絡訪問固定樣本成員和非網絡訪問固定樣本成員,則通常這兩者在調查中的回答將沒有區別。

假設2:共同的支撐。對于網絡訪問固定樣本成員和非網絡訪問固定樣本成員,匹配變量的分布應該有重疊。

1.2 傾向得分加權調整

在基于傾向得分匹配的樣本選擇中,已將目標樣本與網絡訪問固定樣本結合,建立了Logistic回歸模型,估計出匹配樣本第k個單元在目標樣本中的傾向得分為(Xk),則相應的第k個單元在匹配樣本中的傾向得分估計為假設采用傾向得分匹配方法最終選擇的匹配樣本為sM,可以利用匹配樣本單元的傾向得分估計進一步對匹配樣本進行加權調整,從而實現對目標總體的估計??紤]采取以下三種傾向得分加權調整方法:

(1)傾向得分逆加權

其中Y0k為匹配樣本單元k的觀察值,dMk為匹配樣本單元k的基礎權數,由于匹配樣本為非概率樣本,沒有基礎權數,但因其與目標樣本相匹配,所以匹配樣本單元的基礎權數可采用與之相匹配的目標樣本單元的基礎權數,即為目標樣本單元入樣概率的倒數。因本文主要討論的目標樣本為簡單隨機樣本,故每個目標樣本單元的基礎權數均為

(2)傾向得分加權組調整

其中sMg是在第g組中匹配樣本單元的集合。同理,若目標樣本為簡單隨機樣本,則

(3)傾向得分事后分層調整

其中YˉMg為在第g層中匹配單元目標變量均值。

2 模擬研究

為了檢驗采用傾向得分匹配來選擇樣本,并對匹配樣本進行傾向得分加權調整后的估計效果,采用蒙特卡羅模擬進行分析。

2.1 方法

數據生成過程如下:

目標總體協變量:X1~N(0,1),X1?[-1,1];X2~N(0,1),X2?[-1,1];ρ(X1,X2)=-0.6

網絡訪問固定樣本協變量:X1~N(0.8,0.4),X1? [-1,1];X2~N(0.7,0.35),X2?[-1,1];Cov(X1,X2)=0.3

由上可知,網絡訪問固定樣本關于協變量(匹配變量)是可忽略的,且目標總體與網絡訪問固定樣本的協變量分布有重疊,但協變量的分布不同。在模擬中,目標樣本是從目標總體中簡單隨機抽取的樣本量n=1000的樣本,因為E(X1)=E(X2)=0,θ0=E(Y)=0為總體參數真值,從目標總體中抽取的簡單隨機樣本將產生θ0的一個無偏估計;同時考慮網絡訪問固定樣本的規模分別為1500、2000、3000、5000、10000。因為目標變量依賴于協變量X1、X2,因此若X1、X2來自目標總體協變量的分布,則示性變量D=1,若X1、X2來自網絡訪問固定樣本協變量的分布,則示性變量D=0,并進一步由協變量與示性變量估計傾向得分,利用傾向得分匹配選擇匹配樣本,匹配方法為單一無放回最近鄰匹配。最后,采用傾向得分逆加權、加權組調整和事后分層調整三種方法對匹配后的樣本進行加權調整,估計目標總體均值θ0。重復進行蒙特卡羅模擬1000次,并計算1000次模擬數據集上目標總體均值估計的均值(Mean)、標準差(SD)和均方誤差(RMSE)。此外,數據的生成、傾向得分匹配與加權調整均在R軟件中進行。特別地,單一無放回最近鄰匹配采用的是Package Matching中的Match函數進行編程。

2.2 結果

利用未加權匹配樣本以及對匹配樣本采用傾向得分逆加權、加權組調整、事后分層調整后,計算目標總體均值估計的均值(Mean)、標準差(SD)和均方誤差(RMSE),見表1。由表1可見,利用未加權匹配樣本估計目標總體均值,隨著網絡訪問固定樣本規模的增大,其均值絕對值從0.285逐漸減少至0.003;標準差先減小后增大,基本穩定在0.030~0.045之間,估計比較穩定。對匹配樣本采用傾向得分逆加權估計目標總體均值,隨著網絡訪問固定樣本規模的增大,其均值絕對值從0.012先增大后減少至0.071;標準差在0.030~0.055之間波動。對匹配樣本采用傾向得分加權組調整估計目標總體均值,隨著網絡訪問固定樣本規模的增大,其均值絕對值從0.076先減小再增大再減小至0.061;標準差在0.030~0.050之間波動。對匹配樣本采用傾向得分事后分層調整估計目標總體均值,隨著網絡訪問固定樣本規模的增大,其均值絕對值從0.001開始波動最終升至0.003;標準差先減小后逐漸增大,基本穩定在0.025~0.045之間??偟膩碚f,隨著網絡訪問固定樣本規模的增大,各種方法下估計的均值絕對值均在0.001~ 0.285之間,并且各種方法下估計的波動較小,比較穩定。

表1 匹配樣本未加權與三種方法加權調整后的模擬結果

進一步結合表1與圖1來分析各種方法下的RMSE。首先從不同方法的網絡訪問固定樣本規模變化來看RMSE,使用未加權的匹配樣本估計目標總體均值,RMSE隨著網絡訪問固定樣本規模的增大,先急劇減小后趨于穩定(穩定于0.040左右),估計效果變好并逐漸穩定。采用傾向得分逆加權估計目標總體均值,RMSE隨著網絡訪問固定樣本規模的增大,先增大至最高(N=3000,RMSE= 0.160)后逐漸減小,估計效果先變差然后逐漸變好。采用傾向得分加權組調整估計目標總體均值,RMSE隨著網絡訪問固定樣本規模的增大,先減小至最低(N=2000,RMSE=0.038),后增大至最高(N=3000,RMSE=0.124),此后逐漸下降,估計效果呈波動狀態。采用傾向得分事后分層調整估計目標總體均值,RMSE隨著網絡訪問固定樣本規模的增大,變化較為平穩,基本穩定在0.025~0.045之間,估計效果非常好并且非常穩定。

另一方面,從不同網絡訪問固定樣本規模來比較四種方法估計的RMSE,當N=1500和2000時,未加權匹配樣本估計的RMSE均最大,傾向得分事后分層調整估計的RMSE均最小,傾向得分加權組調整與逆加權估計的RMSE均處于兩者之間,可見傾向得分事后分層調整估計的效果最好,傾向得分加權組調整與事后分層調整估計效果次之,未加權匹配樣本估計效果最差。當N=3000和5000時,傾向得分逆加權估計的RMSE均最大,估計效果均最差,其次是傾向得分加權組調整估計,此時未加權匹配樣本估計的RMSE均較小,估計效果較好,但估計效果最好的仍然是傾向得分事后分層調整估計。當N=10000時,傾向得分逆加權與加權組調整估計的RMSE仍然較高,未加權匹配樣本估計與傾向得分事后分層調整估計的RMSE基本相同,估計效果差不多。總的來說,當網絡訪問固定樣本規模較小時,三種加權方法下的估計效果均好于未加權的匹配樣本估計;當網絡訪問固定樣本規模較大時,未加權的匹配樣本估計效果較好;無論網絡訪問固定樣本規模如何變化,傾向得分事后分層調整估計的效果始終都是最好的。這與Rivers(2007)[10]所提到的“當網絡訪問固定樣本規模較小時,對匹配樣本的事后分層調整是有用的,有助于減少不完美匹配所造成的偏差。”想法是一致的。同樣的思想Rivers和Bailey(2009)[5]在“2008年美國全國選舉的匹配樣本推斷”一文中也提到,并且他們進一步指出當網絡訪問固定樣本大小與目標樣本大小的比率比5還小時,有必要對匹配后的樣本加權。本文的模擬結果與之不同的是當網絡訪問固定樣本大小與目標樣本大小的比率比3(非5)小時,有必要對匹配后的樣本加權,可采取傾向得分逆加權、加權組調整和事后分層調整三種方法加權,估計效果均比未加權的匹配樣本估計好。存在不同的原因可能是本文在樣本匹配時采用的匹配方法(傾向得分匹配)與Rivers和Bailey不同,說明若將傾向得分用于樣本選擇階段,可能會出現當網絡訪問固定樣本大小與目標樣本大小的比率(大于等于3)比5還小時,仍可直接利用匹配樣本估計目標總體均值,而無需進一步加權調整。

圖1 網絡訪問固定樣本規模與RMSE之間的關系

3 實證分析

本文采用2014年美國社區調查中阿拉斯加州的數據http://factfinder.census.gov/faces/nav/jsf/pages/searchresults. xhtml?refresh=t,進一步驗證本文所提出的方法。選取阿拉斯加州數據中的年齡、公民身份、婚姻狀況、教育、性別、就業狀況記錄、種族、過去12個月的工資或工資收入共8個變量,6787個個案。其中年齡在1歲以下的記為0,1到99歲分別記為1,…,99;公民身份的取值為1、2、3、4、5,1表示在美國出生,2表示出生于波多黎各、關島、美屬維爾京群島或北馬里亞納群島,3表示出生在國外的美國父母,4表示加入美國國籍的美國公民,5表示非美國公民;婚姻狀況的取值為1、2、3、4、5,1表示已婚,2表示喪偶,3表示離異,4表示分居,5表示未婚或者15歲以下;教育的取值為0到24,0表示小于3歲,1表示沒有上學,2表示幼兒園,3表示學前班,4到9分別表示小學一年級到六年級,10到12分別表示初中一年級到三年級,13、14分別表示高中一年級和二年級,15表示高中三年級但無文憑,16表示普通高中文憑,17表示GED證書或者同等學歷證明,18表示上了大學但不到1年,19表示上了1年或以上的大學但沒有學位,20表示準學士學位或專科畢業證書,21表示學士學位,22表示碩士學位,23表示除了學士學位還有專業學位,24表示博士學位;性別的取值為1、2,1表示男性,2表示女性;就業狀況記錄取值為0到6,0表示小于16歲,1表示聘用的文職人員且正在工作,2表示聘用的文職人員有一個職位,但并沒有工作,3表示失業,4表示武裝部隊且正在工作,5表示武裝部隊有一個職位,但并沒有工作,6表示無勞動能力;種族取值為1到9,1表示白人,2表示黑人或非裔美國人,3表示美國印第安人,4表示阿拉斯加本地人,5表示既是美國印第安人又是阿拉斯加本地人或者不確定是美國印第安人還是阿拉斯加本地人,6表示亞洲人,7表示夏威夷原住民和其他太平洋島民,8表示一些其他種族,9表示兩個或多個種族。過去12個月的工資或工資收入取值為0到999999美元,該變量取值的最大值為366000,最小值為0,均值為20516.589,標準差為33308.300。

假設該6787個單元構成目標總體,從目標總體中簡單隨機抽取一個樣本量為1000的樣本(目標樣本),并以過去12個月的工資或工資收入為目標變量,其他變量為協變量。假設網絡訪問固定樣本為從目標總體中除去目標樣本后的剩余單元中選取(通過網絡招募),其規模分別為1500、2000和3000。在不同的網絡訪問固定樣本規模下,利用傾向得分匹配選擇樣本,匹配方法為單一無放回最近鄰匹配;并對匹配后的樣本采用傾向得分逆加權、加權組調整和事后分層調整三種方法進行加權調整,估計目標總體均值,每種方法重復進行10次,最終計算的目標總體均值估計的均方誤差(RMSE)見表2。

表2 匹配樣本未加權與三種方法加權調整后的RMSE

由表2可見,當網絡訪問固定樣本規模為1500和2000時,采用傾向得分逆加權、加權組調整和事后分層調整三種方法估計目標總體均值的效果均優于未加權時匹配樣本的估計效果,并且傾向得分事后分層調整的估計效果最好,傾向得分加權組調整的估計效果次之。當網絡訪問固定樣本規模為3000時,采用傾向得分事后分層調整的估計效果最好,其次是未加權時匹配樣本的估計,兩者估計的RMSE相差不大,并且采用傾向得分逆加權和加權組調整的估計效果均不如未加權時匹配樣本的估計效果??梢钥吹剑瑢嵶C結果與模擬結果基本一致。

4 結論

本文提出采用傾向得分匹配選擇樣本,并對匹配后的樣本利用傾向得分逆加權、加權組調整和事后分層調整三種方法進行加權調整,從而提高估計的精度,實現對目標總體的統計推斷,并進一步采用蒙特卡羅模擬和實際數據比較各種方法的估計效果。蒙特卡羅模擬與實證研究表明:當網絡訪問固定樣本大小與目標樣本大小的比率比3小時,有必要對匹配后的樣本加權,可采取傾向得分逆加權、加權組調整和事后分層調整三種方法加權,估計效果均比未加權的匹配樣本估計的效果好;當網絡訪問固定樣本大小與目標樣本大小的比率不小于3時,傾向得分事后分層調整與未加權的匹配樣本估計效果較好,可以考慮利用未加權的匹配樣本或者傾向得分事后分層調整進行估計。

本文所提出的方法將傾向得分同時運用于非概率的樣本選擇階段與加權調整階段,提高了估計的精度,豐富了非概率抽樣的統計推斷方法。同時,由于網絡訪問固定樣本具有成本較低,且能得到較快回答的特點,使得本文所提出的方法具有較強的可操作性。

[1]Baker R,Brick J M,Bates N A,et al.Summary Report of the AAPOR Task Force on Nonprobability Sampling[J].Journal of Survey Statis?tics and Methodology,2013,1(2).

[2]Svensson J.Web Panel Surveys--Can They Be Designed and Used in a Scientifically Sound Way?[C].HongKong:59th World Statistics Congress,2013.

[3]Rivers D.Sample Matching--representative Sampling From Internet Panels[J].Polimetrix White Paper Series,2006.

[4]Vavreck L,Rivers D.The 2006 Cooperative Congressional Election Study[J].Journal of Elections,Public Opinion&Parties,2008,18(4). [5]Rivers D,Bailey D.Inference From Matched Samples in the 2008 U. S.NationalElections[J].American Association of Public Re?search--JSM,2009.

[6]Terhanian G,Bremer J.A Smarter Way to Select Respondents for Sur?veys?[J].International Journal of Market Research,2012,54(6).

[7]Rosenbaum P R,Rubin D B.The Central Role of the Propensity Score in Observational Studies for Causal Effects[J].Biometrika,1983,70 (1).

[8]Smith J A,Todd P E.Does Matching Overcome LaLonde’s Critique of Nonexperimental Estimators?[J].Journal of Econometrics,2005,125 (2).

[9]Caliendo M,Kopeinig S.Some Practical Guidance for the Implementa?tion of Propensity Score Matching[J].Journal of Economic Surveys, 2008,22(1).

[10]Rivers D.Sampling for Web Surveys[C].The 2007 Joint Statistical Meetings,2007.

(責任編輯/易永生)

C811

A

1002-6487(2016)21-0004-05

國家社會科學基金資助項目(15BTJ014);中國人民大學2015年度拔尖創新人才培育資助計劃項目

劉 展(1981—),女,湖北宜昌人,博士研究生,研究方向:抽樣調查技術與數據分析。

金勇進(1953—),男,北京人,教授,博士生導師,研究方向:抽樣調查技術與數據分析。

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲性网站| 国产精品无码制服丝袜| 精品成人一区二区三区电影| 欧美、日韩、国产综合一区| 在线免费a视频| 亚洲女同欧美在线| 亚洲精品黄| 99精品免费欧美成人小视频| 国产一级二级三级毛片| 国产成人综合欧美精品久久| 精品国产www| 亚洲免费福利视频| 亚洲视频四区| 538国产视频| 成人免费午夜视频| 思思热在线视频精品| 国产精品女主播| 日韩国产黄色网站| 国产乱人视频免费观看| 国产美女一级毛片| 欧美国产在线看| 久久黄色一级视频| 国产无码精品在线播放| 亚洲精品手机在线| 中文字幕日韩欧美| 国产成人在线无码免费视频| 久久综合干| 四虎免费视频网站| 精品一区二区三区自慰喷水| 国产精品观看视频免费完整版| 亚洲欧美日韩成人高清在线一区| 亚洲无码在线午夜电影| 亚洲人成在线精品| 精品国产乱码久久久久久一区二区| 日韩在线视频网| 自拍偷拍欧美日韩| 亚洲第一黄色网| 51国产偷自视频区视频手机观看| 欧美三级自拍| 久久中文字幕2021精品| 午夜高清国产拍精品| 日韩在线永久免费播放| 亚洲最大情网站在线观看| 国产另类视频| 中文字幕在线视频免费| 久久精品国产免费观看频道| 亚洲国产精品VA在线看黑人| 免费A∨中文乱码专区| 亚洲无码高清一区| 5555国产在线观看| 老司机aⅴ在线精品导航| 免费久久一级欧美特大黄| 欧美黄网在线| 99精品视频九九精品| 最新国产午夜精品视频成人| 国内丰满少妇猛烈精品播| 无码精品国产VA在线观看DVD| 亚洲一区二区三区国产精品 | 国产成人精品一区二区免费看京| 成人在线观看不卡| 在线观看精品自拍视频| 精品免费在线视频| 新SSS无码手机在线观看| 日韩色图在线观看| 中文字幕无线码一区| 青青青亚洲精品国产| 久久久久人妻一区精品| 综合久久久久久久综合网| 怡春院欧美一区二区三区免费| 亚洲三级网站| 欧美亚洲一区二区三区在线| 最新国产高清在线| 美女免费黄网站| 午夜在线不卡| 无码电影在线观看| 国产视频一区二区在线观看| 亚洲日韩高清无码| 国产视频资源在线观看| 免费毛片网站在线观看| 天堂av综合网| 欧美区在线播放| 免费在线成人网|