亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        處理零計(jì)數(shù)過(guò)多數(shù)據(jù)的兩部模型方法介紹*

        2016-12-27 08:49:33哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室150081
        關(guān)鍵詞:泊松參數(shù)估計(jì)對(duì)數(shù)

        哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)

        姜 博 門(mén)志紅 劉匆提 劉 艷△

        ·綜述·

        處理零計(jì)數(shù)過(guò)多數(shù)據(jù)的兩部模型方法介紹*

        哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計(jì)學(xué)教研室(150081)

        姜 博 門(mén)志紅 劉匆提 劉 艷△

        在公共衛(wèi)生、社會(huì)科學(xué)等領(lǐng)域的調(diào)查研究中,定量數(shù)據(jù)可能會(huì)包含過(guò)多零值,定性數(shù)據(jù)也有可能存在著零膨脹[1],如果仍按照一般模型(如泊松回歸或線(xiàn)性回歸)的數(shù)據(jù)分布假設(shè)進(jìn)行分析,會(huì)導(dǎo)致參數(shù)估計(jì)產(chǎn)生偏倚[1-2]。針對(duì)此類(lèi)問(wèn)題,兩部模型法是一種比較常用的擬合模型方式。目前,兩部模型法包括應(yīng)用于定性數(shù)據(jù)的 Hurdle模型、零膨脹模型(zero-inflated model)等,以及應(yīng)用于定量數(shù)據(jù)的半連續(xù)數(shù)據(jù)(semicontiunuous data)兩部模型等[3]。

        定性數(shù)據(jù)兩部模型

        對(duì)于定性數(shù)據(jù)的零膨脹問(wèn)題,可以采用混合泊松回歸模型、birth process模型、兩部模型等[4]。兩部模型的基本思想是將數(shù)據(jù)分成兩個(gè)部分,第一部分考慮事件是否發(fā)生,第二部分考慮事件發(fā)生的次數(shù),主要包括 hurdle模型[5]與零膨脹模型[6]。如 Neelon等調(diào)查某一地區(qū)內(nèi)發(fā)生急診科就醫(yī)的情況,因其中大部分被調(diào)查者急診科就醫(yī)為零計(jì)數(shù),故使用hurdle模型進(jìn)行分析[7];Rose等使用零膨脹模型分析疫苗不良事件發(fā)生情況[8]。兩種模型主要的區(qū)別在于零值來(lái)源的假設(shè)處理不同[9],Hurdle模型更適合處理數(shù)據(jù)的零計(jì)數(shù)僅來(lái)自于抽樣的零,而零膨脹模型更適合處理數(shù)據(jù)的零計(jì)數(shù)來(lái)自結(jié)構(gòu)的零與抽樣的零,實(shí)際應(yīng)用中還要根據(jù)研究目的、專(zhuān)業(yè)解釋及數(shù)據(jù)的具體分布來(lái)選擇合適的兩部模型[10]。

        1.Hurdle模型

        Mullahy對(duì)Hurdle模型進(jìn)行了深入討論[5],第一部分考慮零計(jì)數(shù)是否發(fā)生,服從二項(xiàng)分布假設(shè),連接函數(shù)一般為 logit、clog、probit等[9];第二部分再對(duì)非零計(jì)數(shù)部分進(jìn)行分析,一般假設(shè)其服從泊松分布,擬合泊松回歸模型,即泊松hurdle模型。此外,第二部分還有其他的分布假設(shè)類(lèi)型,如負(fù)二項(xiàng)分布[10]、廣義泊松分布[11]等。Hurdle模型參數(shù)估計(jì)是將兩個(gè)部分作為相互獨(dú)立的模型進(jìn)行估計(jì),其具體方法包括多種,較常用的是極大似然法,也有使用廣義估計(jì)方程(generalized estimating equations,GEE)方法分別對(duì)兩部分進(jìn)行參數(shù)估計(jì)[12]。

        以泊松hurdle模型為例,其第一部分為p(y=0|x)=π,0≤π≤1,第二部分為p(y>0|x)=(1-中π為發(fā)生零計(jì)數(shù)的概率,μ為事件發(fā)生的平均次數(shù),y為事件發(fā)生次數(shù)[5];其參數(shù)估計(jì)的對(duì)數(shù)似然函數(shù)為:

        2.零膨脹模型(zero-inflated model)

        Lambert[6,10]首次提出零膨脹泊松模型,模型的假設(shè)是數(shù)據(jù)服從一種混合分布[3],數(shù)據(jù)分為兩個(gè)部分,同時(shí)數(shù)據(jù)中的“0”也被分成兩種組成部分,第一部分仍然考慮是否發(fā)生零計(jì)數(shù)情況,假設(shè)服從二項(xiàng)分布,連接函數(shù)一般為logit,該部分的零值稱(chēng)為不可能發(fā)生的零或結(jié)構(gòu)的零,即不存在發(fā)生可能所導(dǎo)致的零值;第二部分考慮可能發(fā)生事件的情況,即抽樣的零,也就是有可能發(fā)生但并未發(fā)生事件,常用的模型假設(shè)為泊松分布,負(fù)二項(xiàng)分布[10]、廣義泊松分布[13]、廣義冪級(jí)數(shù)分布(generalized power series distribution)[14]等。零膨脹模型的參數(shù)估計(jì)依然較常采用極大似然法,可使用最大期望算法(expectation maximization algorithm,EM)[6]、牛頓迭代法(Newton-raphson algorithm)[15]進(jìn)行極大似然運(yùn)算,同樣也可以使用廣義估計(jì)方程方法估計(jì)相應(yīng)參數(shù)[16],或應(yīng)用貝葉斯方法進(jìn)行參數(shù)估計(jì)[17]。

        以零膨脹泊松模型[6,10]為例,其第一部分為 p意義與 hurdle模型中相同;其對(duì)數(shù)似然函數(shù)表示為[6,10]exp(x1β2)]。

        3.模型發(fā)展

        伴隨著定性數(shù)據(jù)兩部模型方法的不斷發(fā)展,兩部模型的應(yīng)用也不斷擴(kuò)展。已有一些研究在零膨脹模型的基礎(chǔ)上提出了半?yún)?shù)零膨脹模型[18]、修正零膨脹泊松模型[1]等模型,還有一些研究也提出適合于其他類(lèi)型數(shù)據(jù)應(yīng)用的模型,如縱向定性數(shù)據(jù)的零膨脹增長(zhǎng)曲線(xiàn)模型(zero-inflated growth curve model)[19]、層次結(jié)構(gòu)數(shù)據(jù)的多水平零膨脹泊松模型(multi-level zero-inflated poisson regression modelling)[20]、空間泊松 hurdle模型(spatial poisson hurdle model)[7]、時(shí)空 hurdle模型(spatiotemporal hurdle models)[21]等等。

        定量數(shù)據(jù)兩部模型

        在大型隊(duì)列研究或橫斷面調(diào)查中,醫(yī)療費(fèi)用經(jīng)常呈現(xiàn)為一種半連續(xù)數(shù)據(jù),即在調(diào)查人群中存在大量的調(diào)查對(duì)象未發(fā)生醫(yī)療行為即不產(chǎn)生費(fèi)用,發(fā)生醫(yī)療行為調(diào)查對(duì)象的費(fèi)用呈正偏態(tài)分布。對(duì)于這類(lèi)半連續(xù)數(shù)據(jù),過(guò)多的零值導(dǎo)致費(fèi)用分布右偏嚴(yán)重,一般傳統(tǒng)模型假設(shè)分布類(lèi)型難以滿(mǎn)足這類(lèi)半連續(xù)數(shù)據(jù)[22],可以選用Tobit模型[23]、樣本選擇(sample selection)模型[24]、Cox比例風(fēng)險(xiǎn)回歸模型[25]、兩部模型[26]進(jìn)行數(shù)據(jù)分析。如Finkelstein等使用美國(guó)MEPS(medical expenditure panel surveys)數(shù)據(jù),應(yīng)用兩部模型分析肥胖與醫(yī)療費(fèi)用的關(guān)系[27];Bock等使用隊(duì)列人群的橫斷面調(diào)查數(shù)據(jù),應(yīng)用兩部模型分析德國(guó)老年人口自費(fèi)健康醫(yī)療服務(wù)的不公平性問(wèn)題[28]。

        1.模型組成

        與定性數(shù)據(jù)的hurdle模型類(lèi)似,定量數(shù)據(jù)的兩部模型依然是將半連續(xù)數(shù)據(jù)分成兩個(gè)部分,第一部分將應(yīng)變量作為二項(xiàng)分布處理,如是否發(fā)生醫(yī)療行為,構(gòu)建二項(xiàng)分布概率模型,常用logistic、probit回歸模型分析是否發(fā)生醫(yī)療行為的影響因素;第二部分對(duì)發(fā)生醫(yī)療行為即大于零的數(shù)據(jù)部分,一般常假設(shè)費(fèi)用滿(mǎn)足對(duì)數(shù)正態(tài)分布[29],進(jìn)行數(shù)據(jù)轉(zhuǎn)換后擬合廣義線(xiàn)性回歸模型,分析醫(yī)療費(fèi)用的影響因素,gama分布也可以作為第二部分模型的分布假設(shè)[30]。當(dāng)假設(shè)費(fèi)用數(shù)據(jù)服從對(duì)數(shù)正態(tài)分布時(shí),如果存在個(gè)體真實(shí)費(fèi)用為零的情況時(shí),由于零值不能進(jìn)行對(duì)數(shù)轉(zhuǎn)化,若將其歸為第二部分,則軟件運(yùn)算時(shí)會(huì)將其當(dāng)作缺失值而忽略。

        2.參數(shù)估計(jì)

        半連續(xù)數(shù)據(jù)兩部模型的參數(shù)估計(jì)方法較多,選擇也較為靈活,主要是根據(jù)調(diào)查目的及各部分所選用的模型而定,較常用的依然是極大似然法[3]。但隨著兩部模型的應(yīng)用拓展,不同種類(lèi)兩部模型不斷被開(kāi)發(fā)出來(lái),其參數(shù)估計(jì)方法也各不相同,如MK Olsen與JL Schafer應(yīng)用兩部模型分析縱向數(shù)據(jù),擬合兩部隨機(jī)效應(yīng)模型,使用極大似然法估計(jì)固定效應(yīng),基于高階拉普拉斯法的類(lèi)Fisher評(píng)分法進(jìn)行方差估計(jì)[31];Neelon等使用貝葉斯兩部潛類(lèi)別模型(two-part latent class model),應(yīng)用馬氏鏈蒙特卡羅(MCMC)法擬合兩部模型[32]。

        以第一部分使用probit回歸模型、第二部分使用對(duì)數(shù)變換后線(xiàn)性模型為例,擬合兩步法模型,其各部函數(shù)及似然函數(shù)表示為[29]:第一部分probit回歸模型為yi=β1xi+e1i,e1i~N(0,1),第二部分對(duì)數(shù)轉(zhuǎn)化線(xiàn)性模型為 log(yi|yi>0)=β2xi+e2i,e2i~N(0,σ2),表示第i個(gè)觀(guān)測(cè)。

        3.模型預(yù)測(cè)

        半連續(xù)數(shù)據(jù)兩部模型可以用來(lái)預(yù)測(cè)個(gè)體的醫(yī)療費(fèi)用情況,具體方法是通過(guò)第一部分的模型預(yù)測(cè)個(gè)體可能發(fā)生醫(yī)療行為的概率,再通過(guò)第二部分模型預(yù)測(cè)所產(chǎn)生費(fèi)用的期望,兩個(gè)部分模型的估計(jì)值相乘即可得出個(gè)體醫(yī)療費(fèi)用的估計(jì)值,其函數(shù)可以表示為E(yi|xi)=prob(yi>0|xi)·E(yi|xi,yi>0)[33]。通常,第一部分產(chǎn)生費(fèi)用的概率估計(jì)值采用probit回歸或logistic回歸模型進(jìn)行計(jì)算,而第二部分具體費(fèi)用的估計(jì)值則可根據(jù)數(shù)據(jù)的分布選擇合適的計(jì)算方式。若假設(shè)第二部分?jǐn)?shù)據(jù)服從對(duì)數(shù)正態(tài)分布,可依據(jù)最小二乘法(OLS)估計(jì)相應(yīng)統(tǒng)計(jì)量,再對(duì)估計(jì)值取指數(shù),表示為E(yi|χi,yi>0)=exp(β2xi+σ2/2);若應(yīng)用對(duì)數(shù)轉(zhuǎn)化后,誤差項(xiàng)仍無(wú)法滿(mǎn)足正態(tài)分布,使用最小二乘法進(jìn)行估計(jì)將出現(xiàn)偏倚,故Duan提出一種非參數(shù)的估計(jì)方法,即 Smearing估計(jì)法[34]。Smearing估計(jì)法對(duì)數(shù)據(jù)分布沒(méi)有特定假設(shè),僅要求誤差項(xiàng)獨(dú)立同分布,模型的估計(jì)值可表示為[29]

        最小二乘法、Smearing估計(jì)法、GLMs模型三種估計(jì)方式各有相應(yīng)的適用條件,需根據(jù)數(shù)據(jù)分布的具體情況選擇合適的方法,若數(shù)據(jù)轉(zhuǎn)換后符合相應(yīng)分布(如對(duì)數(shù)正態(tài)分布)或峰度較高,則OLS法估計(jì)更為準(zhǔn)確;若不符合相應(yīng)分布,可選用Smearing法;轉(zhuǎn)換后數(shù)據(jù)峰度較低時(shí),可使用 GLMs[3]。

        4.模型發(fā)展

        兩部模型的簡(jiǎn)單與靈活性決定其被廣泛應(yīng)用于半連續(xù)數(shù)據(jù)分析[26,36]。1981年 Manning等將兩部模型應(yīng)用于衛(wèi)生經(jīng)濟(jì)領(lǐng)域,研究醫(yī)療保險(xiǎn)費(fèi)用[37]。為滿(mǎn)足多種數(shù)據(jù)分布類(lèi)型的發(fā)展,兩部模型不斷拓展,第二部分模型已提出廣義 gamma分布、box-cox變換等[38];處理縱向數(shù)據(jù)、層次結(jié)構(gòu)數(shù)據(jù)的兩部模型已提出兩部隨機(jī)效應(yīng)模型[26]、多水平兩部模型[39]、基于貝葉斯方法的雙變量?jī)刹磕P?(bivariate two-part model)[40]等;甚至在兩部模型理論的基礎(chǔ)上,提出三部模型[2]、四部模型[29]等等。

        相關(guān)軟件

        在統(tǒng)計(jì)學(xué)軟件方面,R軟件可以使用pscl等程序包[41],SAS軟件可以通過(guò) GENMOD、NLMIXED等過(guò)程,Stata軟件可以通過(guò) ZIP、HPLOGIT等命令、Win-BUGS軟件[42]通過(guò)貝葉斯算法擬合構(gòu)建定性數(shù)據(jù)兩部模型;SAS軟件通過(guò)GLIMMIX和MCMC(貝葉斯算法)等過(guò)程、Stata通過(guò)GLM等命令、WinBUGS軟件[43](貝葉斯算法)構(gòu)建定量數(shù)據(jù)兩部模型。

        結(jié) 論

        兩部模型可有效地解決公共衛(wèi)生、社會(huì)科學(xué)等領(lǐng)域調(diào)查研究中出現(xiàn)的零計(jì)數(shù)過(guò)多問(wèn)題,其靈活的數(shù)據(jù)分布類(lèi)型假設(shè)與參數(shù)估計(jì)方式,使其擁有較為廣泛的適用范圍。定性數(shù)據(jù)兩部模型的使用選擇,要充分考慮數(shù)據(jù)中零計(jì)數(shù)的來(lái)源方式、數(shù)據(jù)分布類(lèi)型、參數(shù)估計(jì)方式、研究目的及專(zhuān)業(yè)知識(shí)等方面從而有效減小偏倚;定量數(shù)據(jù)兩部模型使用時(shí),要選擇合適的數(shù)據(jù)分布類(lèi)型、參數(shù)估計(jì)方式及模型預(yù)測(cè)的計(jì)算方式。目前,仍有較多研究和項(xiàng)目關(guān)注于兩部模型法的計(jì)算與使用,使其在零計(jì)數(shù)過(guò)多問(wèn)題的處理上繼續(xù)保持較大的應(yīng)用價(jià)值及推廣意義。

        [1]郭念國(guó).零膨脹泊松模型的改進(jìn)在零次索賠建模中的應(yīng)用.統(tǒng)計(jì)與信息論壇,2010,25(7):22-25.

        [2]葉玲瓏.葉玲瓏.基于兩部模型的家庭醫(yī)療需求與消費(fèi)結(jié)構(gòu)研究.廈門(mén):廈門(mén)大學(xué),2014.

        [3]Neelon B.Two-Part Models for Zero-Modified Count and Semicontinuous Data.Duke University,2013.

        [4]Ridout M,Demétrio CG,Hinde J.Models for count data with many zeros//Proceedings of the XIXth International Biometric Conference,1998,19:179-192.

        [5]Mullahy J.Specification and testing of some modified count data models.J Econom,1986,33(3):341-365.

        [6]Lambert D.Zero-inflated Poisson regression,with an application to defects in manufacturing.Technometrics,1992,34(1):1-14.

        [7]B Neelon PG,Loebs PF.A spatial Poisson hurdle model for exploring geographic variation in emergency department visits,2013,176(2):389-413.

        [8]CE Rose SWM,Wannemuehler KA.On the use of zero-inflated and hurdle models for modeling vaccine adverse event count data,2006,16(4):463-481.

        [9]曾平,趙晉芳,劉桂芬.居民就診次數(shù)的Hurdle模型分析.中國(guó)衛(wèi)生統(tǒng)計(jì),2010,27(6):603-605.

        [10]原靜,劉桂芬,薛玉強(qiáng).零膨脹計(jì)數(shù)資料模型選擇與比較.中國(guó)衛(wèi)生統(tǒng)計(jì),2011,28(4):354-356,360.

        [11]SE Saffari RA,Greene W.Investigating the impact of excess zeros on hurdle-generalized Poisson regression model with right censored count data.Statistica Neerlandica,2013,67(1):67-80.

        [12]MJDobbie AHW.Theory&Methods:Modelling Correlated Zero-inflated Count Data.2001,43(4):431-444.

        [13]H Joe RZ.Generalized Poisson distribution:the property of mixture of Poisson and comparison with negative binomial distribution.Biometrical Journal,2005,47(2):219-229.

        [14]Statistics YZF-Ci.Model Selection of Zero-Inflated Generalized Power Series Distribution with Missing Responses.Theory and Methods,2012,41(6):1013-1028.

        [15]閆鳴宇.零膨脹模型的若干問(wèn)題研究.蘇州:蘇州大學(xué),2011.

        [16]DB Hall ZZ.Marginal models for zero inflated clustered data.Statistical Modelling,2004,4(3):161-180.

        [17]Statistics JR-Ci.Bayesian analysis of zero-inflated distributions.Theory and Methods,2003,32(2):281-289.

        [18]KF Lam HX,Cheung YB.Semiparametric Analysis of Zero-Inflated Count Data.Biometrics,2006,62(4):996-1003.

        [19]Liu H.Growth curve models for zero-inflated count data:An application to smoking behavior,2007,14(2):247-279.

        [20]AH Lee KW,JA Scott KY.Multi-level zero-inflated Poisson regression modelling of correlated count data with excess zeros,2006,15(1):47-61.

        [21]B Neelon HHC,Ling Q.Spatiotemporal hurdle models for zero-inflated count data:Exploring trends in emergency department visits.Stat Methods Med Res,2014:0962280214527079.

        [22]Liu L,Strawderman RL,Cowen ME,et al.A flexible two-part random effects model for correlated medical costs.JHealth Econ,2010,29(1):110-23.

        [23]Tobin J.Estimation of relationships for limited dependent variables.Econometrica:journal of the Econometric Society,1958,26(1):24-36.

        [24]Madden D.Sample selection versus two-part models revisited:The case of female smoking and drinking.Journal of Health Economics,2008,27(2):300-307.

        [25]A Basu WGM,Mullahy J.Comparing alternative models:log vs Cox proportional hazard.Health Econ,2004,13(8):749-765.

        [26]Tian L,Huang J.A two-part model for censored medical cost data.Stat Med,2007,26(23):4273-92.

        [27]EA Finkelstein JGT,JW Cohen WD.Annual medical spending attributable to obesity:payer-and service-specific estimates.Health Aff(Millwood),2009,28(5):w822-w831.

        [28]JO Bock HM,Brenner H.Inequalities in out-of-pocket payments for health care services among elderly Germans-results of a populationbased cross-sectional study.Int JEquity Health,2014,13(1):3.

        [29]N Duan JPN,CN Morris WGM.A comparison of alternative models for the demand for medical care:Rand Corporation,1982:15-32.

        [30]周虹,余松林,向蕙云.兩部模型的基本原理與SAS GLIMMIX過(guò)程的配合.中國(guó)衛(wèi)生經(jīng)濟(jì),2008,27(10):19-22.

        [31]MK Olsen JLS.A two-part random-effects model for semicontinuous longitudinal data.2001,96(454):730-745.

        [32]Neelon B,O'Malley AJ,Normand SL.A bayesian two-part latent class model for longitudinal medical expenditure data:assessing the impact of mental health and substance abuse parity.Biometrics,2011,67(1):280-289.

        [33] Buntin MB,Zaslavsky AM.Too much ado about two-part models and transformation?Comparing methods of modeling Medicare expenditures.JHealth Econ,2004,23(3):525-542.

        [34]Duan N.Smearing estimate:a nonparametric retransformation method.J Am Stat Assoc,1983,78(383):605-610.

        [35] WG Manning JM.Estimating log models:to transform or not to transform.J Health Econ,2001,20(4):461-494.

        [36]Mullahy J.Much ado about two:reconsidering retransformation and the two-part model in health econometrics.J Health Econ,1998,17(3):247-281.

        [37]WG Manning CNM.A two-part model of the demand for medical care:preliminary results from the health insurance study.Health,economics,and health economics,1981:103-123.

        [38]L Liu RLS,Johnson BA.Analyzing repeated measures semi-continuous data,with application to an alcohol dependence study.Stat Methods Med Res,2012:0962280212443324.

        [39]L Liu JZM,Johnson BA.A multi-level two-part random effects model,with application to an alcohol-dependence study.Stat Med,2008,27(18):3528-3539.

        [40]Su L,Brown S,Ghosh P,et al.Modelling Household Debt and Financial Assets:A Bayesian Approach to a Bivariate Two-Part Model.The University of Sheffield,Department of Economics,2012.

        [41]A Zeileis CK,Jackman S.Regression models for count data in R.Journal of Statistical Software,2008,27(8).[42]BH Neelon AJO,Normand S.A Bayesian model for repeated measures zero-inflated count data with application to outpatient psychiatric service use.Stat Modelling,2010,10(4):421-439.

        [43]P Ghosh PSA.A Bayesian analysis for longitudinal semicontinuous data with an application to an acupuncture clinical trial.Computational statistics&data analysis,2009,53(3):699-706.

        *:國(guó)家自然科學(xué)基金(81172741;30972537)

        △通信作者:劉艷,E-mail:liuyan@ems.hrbmu.edu.cn

        鄧 妍)

        猜你喜歡
        泊松參數(shù)估計(jì)對(duì)數(shù)
        基于泊松對(duì)相關(guān)的偽隨機(jī)數(shù)發(fā)生器的統(tǒng)計(jì)測(cè)試方法
        含有對(duì)數(shù)非線(xiàn)性項(xiàng)Kirchhoff方程多解的存在性
        基于新型DFrFT的LFM信號(hào)參數(shù)估計(jì)算法
        指數(shù)與對(duì)數(shù)
        帶有雙臨界項(xiàng)的薛定諤-泊松系統(tǒng)非平凡解的存在性
        指數(shù)與對(duì)數(shù)
        對(duì)數(shù)簡(jiǎn)史
        Logistic回歸模型的幾乎無(wú)偏兩參數(shù)估計(jì)
        基于向前方程的平穩(wěn)分布參數(shù)估計(jì)
        基于競(jìng)爭(zhēng)失效數(shù)據(jù)的Lindley分布參數(shù)估計(jì)
        女人被男人爽到呻吟的视频| 人妻av一区二区三区高| 精品国产三级国产av| 国产毛片视频一区二区| 美女把尿囗扒开让男人添| 免费一区二区三区久久| 久久无码高潮喷水抽搐| 日本视频一区二区三区在线| 中文字幕亚洲精品无码| 亚洲欧美精品伊人久久| 国产福利小视频91| 尤物精品国产亚洲亚洲av麻豆 | 国产精品一区二区三区卡| 特黄特色的大片观看免费视频 | 精品无码日韩一区二区三区不卡 | 人人妻人人澡人人爽人人精品浪潮 | 国产免费一区二区三区在线视频| 伊人情人色综合网站| 天天做天天爱天天爽综合网 | 在线观看亚洲你懂得| 亚洲国产av高清一区二区三区| 欧美大胆性生话| 人妻少妇av无码一区二区| 精品日韩欧美一区二区三区在线播放| 免费人妻精品一区二区三区| 人妻洗澡被强公日日澡电影| 亚洲欧美精品伊人久久| 亚洲毛片av一区二区三区| 国产精品专区第一页天堂2019| 精品日产卡一卡二卡国色天香| 国内精品人妻无码久久久影院94| 国产精品av免费网站| 亚洲成av人片一区二区密柚| 无码粉嫩虎白一线天在线观看 | 精品中文字幕制服中文| 免费观看在线视频播放| 国产精品白浆在线观看免费| 1000部精品久久久久久久久| 91国产自拍视频在线| 中文字幕亚洲乱码熟女1区| 无码骚夜夜精品|