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

重復(fù)測(cè)量診斷試驗(yàn)的ROC曲線廣義線性混合效應(yīng)模型*

2017-03-09 08:34:16馬春桃田茂再
關(guān)鍵詞:效應(yīng)測(cè)量模型

馬春桃 熊 巍 田茂再,3,4

·論著·

重復(fù)測(cè)量診斷試驗(yàn)的ROC曲線廣義線性混合效應(yīng)模型*

馬春桃1熊 巍2田茂再1,3,4

目的 針對(duì)重復(fù)測(cè)量診斷數(shù)據(jù),為同時(shí)考慮協(xié)變量對(duì)診斷試驗(yàn)準(zhǔn)確性評(píng)價(jià)的影響,度量重復(fù)測(cè)量數(shù)據(jù)間的相關(guān)性,本文探索新的ROC曲線的建模方法。方法 通過廣義線性混合效應(yīng)模型對(duì)ROC曲線進(jìn)行模擬,并采用貝葉斯參數(shù)估計(jì)方法,利用WinBUGS軟件予以實(shí)現(xiàn),進(jìn)而計(jì)算不同協(xié)變量取值下的ROC曲線下面積(AUC)以對(duì)診斷試驗(yàn)結(jié)果進(jìn)行評(píng)價(jià)。結(jié)果 實(shí)例數(shù)據(jù)分析結(jié)果表明,基于廣義線性混合效應(yīng)模型的ROC曲線建模方法可以有效地刻畫重復(fù)測(cè)量診斷試驗(yàn)數(shù)據(jù),給出更有解釋意義的回歸參數(shù),提供臨床分析的參考依據(jù)。結(jié)論 基于廣義線性混合效應(yīng)的ROC曲線模型在解決重復(fù)測(cè)量診斷試驗(yàn)的準(zhǔn)確度評(píng)價(jià)問題起著至關(guān)重要的作用。

ROC曲線 混合效應(yīng) 貝葉斯 重復(fù)測(cè)量 診斷試驗(yàn)

在醫(yī)學(xué)診斷等領(lǐng)域中,為了準(zhǔn)確評(píng)價(jià)某一診斷方法的準(zhǔn)確度,試驗(yàn)常以重復(fù)測(cè)量的方式進(jìn)行[1]。重復(fù)測(cè)量診斷試驗(yàn)是指一個(gè)患者接受多個(gè)試驗(yàn)或者某試驗(yàn)被同一患者重復(fù)接受多次的情況。同一患者的重復(fù)測(cè)量值之間是非獨(dú)立的,后一次的診斷試驗(yàn)結(jié)果可能受到前面試驗(yàn)結(jié)果的影響,診斷試驗(yàn)結(jié)果在所測(cè)量的時(shí)間范圍內(nèi)可能成趨勢(shì)性變化,不同患者之間的試驗(yàn)結(jié)果一般是獨(dú)立的[2]。為正確區(qū)分人群健康狀態(tài)、及時(shí)干預(yù)病人疾病進(jìn)程和治療,重復(fù)測(cè)量診斷試驗(yàn)的準(zhǔn)確性是研究的關(guān)鍵。為此,應(yīng)綜合考慮混合效應(yīng)模型和ROC曲線分析方法度量重復(fù)測(cè)量診斷數(shù)據(jù)間的相關(guān)關(guān)系。

ROC曲線(receiver operating characteristic curves)方法已廣泛應(yīng)用于臨床診療醫(yī)學(xué)領(lǐng)域,是描述和比較診斷試驗(yàn)的一種重要的綜合評(píng)價(jià)方法。目前,有關(guān)診斷試驗(yàn)準(zhǔn)確性最常用的指標(biāo)為靈敏度和特異度,ROC曲線通過構(gòu)圖能夠反映不同臨界點(diǎn)下靈敏度和特異度的關(guān)系,并可通過計(jì)算ROC曲線下面積(AUC)判斷診斷試驗(yàn)的效率。但是常規(guī)的ROC分析方法不能用來處理重復(fù)測(cè)量診斷數(shù)據(jù),在實(shí)際問題研究中存在些許缺陷:或是建立的隨機(jī)效應(yīng)模型不包含協(xié)變量,無法考慮多種混雜因素(協(xié)變量)對(duì)診斷試驗(yàn)準(zhǔn)確度的影響;或是建立的模型為ROC曲線的間接模型,難以給出回歸參數(shù)的合理性解釋。例如一類常用的間接模型就是先計(jì)算協(xié)變量組合的準(zhǔn)確度指標(biāo)(如AUC),再建立協(xié)變量對(duì)該指標(biāo)影響的線性回歸模型。而該類間接模型,其分析對(duì)象主要為ROC曲線的某種綜合指標(biāo)(如AUC),無法直接討論協(xié)變量的影響;并且協(xié)變量不含有“疾病狀態(tài)”這個(gè)重要的指示變量;與此同時(shí)該模型只能模擬協(xié)變量為分類變量的情形,這也要求協(xié)變量組合中有足夠的患者數(shù)據(jù)來計(jì)算準(zhǔn)確度指標(biāo)。

混合效應(yīng)模型不僅可以分析重復(fù)測(cè)量數(shù)據(jù),還能處理諸如含有缺失數(shù)據(jù)或不均衡數(shù)據(jù)的資料。混合效應(yīng)模型既考慮了固定效應(yīng)也考慮了隨機(jī)效應(yīng),能夠有效刻畫數(shù)據(jù)不獨(dú)立的情況,使得信息利用更加充分;另外該模型不僅能夠用于說明變量水平及變化趨勢(shì),也適用于分析協(xié)變量存在的情況。廣義線性混合效應(yīng)模型是廣義線性模型和混合效應(yīng)模型的推廣,對(duì)響應(yīng)變量是離散變量和分類變量的情形均適用,并通過在模型中加入隨機(jī)效應(yīng)項(xiàng)來刻畫數(shù)據(jù)內(nèi)部的相關(guān)性、異質(zhì)性及過離散等問題。因此,廣義線性混合效應(yīng)模型非常適用于處理重復(fù)測(cè)量診斷試驗(yàn)數(shù)據(jù)。

為有效解決以上難題,本文給出適合重復(fù)診斷測(cè)量試驗(yàn)的ROC模型,考慮利用廣義線性混合效應(yīng)模型對(duì)ROC曲線進(jìn)行模擬。該方法與傳統(tǒng)方法相比,具有一定優(yōu)勢(shì):一方面該方法直接對(duì)ROC曲線建模,不僅給出更有解釋意義的回歸參數(shù),還考慮了協(xié)變量和隨機(jī)效應(yīng)的多重影響;另一方面該模型適用廣泛,對(duì)離散和連續(xù)的協(xié)變量均適用,對(duì)診斷試驗(yàn)結(jié)果的分布也沒有特定要求。這為醫(yī)學(xué)診斷試驗(yàn)評(píng)價(jià)的ROC曲線分析提供了新的思路,為臨床試驗(yàn)提供一定依據(jù)。

原理與方法

1.ROC曲線廣義線性混合效應(yīng)模型

通過廣義線性混合效應(yīng)模型,已知基準(zhǔn)函數(shù)向量h(p)和協(xié)變量值x,建立隨機(jī)效應(yīng)回歸模型:

ROCx(p)=g[γ′h(p)+β′x+λ′z]

其中,g(·)是已知的連接函數(shù),且g(·)必須是嚴(yán)格單調(diào)上升或下降的函數(shù),這樣才能確保ROC曲線在(0,1)區(qū)間內(nèi)是單調(diào)遞增的;h(p)是定義ROC曲線的形狀和位置的已知的基準(zhǔn)函數(shù)向量;γ是假陽(yáng)性率FPRs影響相應(yīng)真陽(yáng)性率TPRs的未知參數(shù)向量;β是協(xié)變量對(duì)ROC曲線影響的回歸參數(shù)向量;λ為隨機(jī)效應(yīng)向量,且假定λ~N(0,D)。一般地,連接函數(shù)選取probit連接函數(shù),即g=Φ(·),基準(zhǔn)函數(shù)向量h(p)=[1,Φ-1(p)]′,Φ(·)指累積標(biāo)準(zhǔn)正態(tài)分布函數(shù),Φ-1(·)指累積標(biāo)準(zhǔn)正態(tài)分布函數(shù)的逆函數(shù)。

ROC曲線下面積(AUC)是評(píng)價(jià)診斷試驗(yàn)最常用的一個(gè)指標(biāo),它表示診斷系統(tǒng)中“患病”(陽(yáng)性)和“不患病”(陰性)診斷結(jié)果分布與“金標(biāo)準(zhǔn)”的重疊程度,體現(xiàn)了診斷試驗(yàn)的價(jià)值,面積越大診斷價(jià)值越高。理論上AUC的取值范圍為0.5~1,兩端點(diǎn)分別表示完全無價(jià)值的診斷及完善的診斷。一般ROC曲線面積在0.50~0.70,表示診斷準(zhǔn)確度低,在0.70~0.90,表示診斷準(zhǔn)確度中等,0.90以上,認(rèn)為診斷準(zhǔn)確度較高。

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

本文采用馬爾科夫蒙特卡洛(Markov chain Monte Carlo,MCMC)的貝葉斯方法來估計(jì)ROC曲線廣義線性混合效應(yīng)模型。相比于極大似然估計(jì)方法,該方法更加靈活準(zhǔn)確,且用于估計(jì)的隨機(jī)效應(yīng)變量個(gè)數(shù)可以是任意的。本文用WinBUGS軟件來進(jìn)行計(jì)算。具體的參數(shù)估計(jì)過程如下[3-5]:

第一步,給定假陽(yáng)性率集,Γ=(p)

基于Alonzo和Pepe(2001)的模擬研究,可選擇50個(gè)等間距的FPRs,即Γ=(1/51,…,50/51),不僅可以保證模型參數(shù)估計(jì)的有效性和穩(wěn)健性,還能節(jié)約模型估計(jì)的運(yùn)算時(shí)間,實(shí)現(xiàn)用較小的假陽(yáng)性率集獲得較大假陽(yáng)性率集的統(tǒng)計(jì)功效。

第二步,確定參數(shù)γ,β,D的先驗(yàn)分布

(γ|A)~U(-A,A)

(β|B)~U(-B,B)

(D-1|V,v)~Wishart(V-1,v)

第三步,MCMC過程

為獲得γ,β,D的后驗(yàn)分布,此處采用Gibbs抽樣方法,迭代更新方法如下,對(duì)于第k步,k=1,…,n,直至收斂:

(1)(γ|λ(k),p,y)

(2)(β|λ(k),y)

(3)(D|λ(k))

(4)(λ|γ(k),β(k),D(k),y)

通過數(shù)值迭代產(chǎn)生λ(k+1)。

第四步,繪制ROC曲線并計(jì)算曲線下面積AUC

根據(jù)模型參數(shù)的貝葉斯估計(jì)結(jié)果,可以繪制不同協(xié)變量組合下的ROC曲線,并計(jì)算其AUC,評(píng)價(jià)診斷試驗(yàn)的準(zhǔn)確度。

實(shí)例分析

某研究組為探討在考慮體質(zhì)指數(shù)(BMI)、每天的飲酒量等影響因素后,谷丙轉(zhuǎn)氨酶ALT(U/L)對(duì)肝膽病的診斷價(jià)值,在某醫(yī)院收集了52名疑似肝炎患者的資料。經(jīng)臨床診斷,患者31例,非患者21例,并對(duì)他們進(jìn)行了10次重復(fù)測(cè)量。數(shù)據(jù)資料的下載地址為:http://archive.ics.uci.edu/ml/datasets.html。本研究收集分析的部分?jǐn)?shù)據(jù)見表1和表2。表3 給出了幾個(gè)重要特征變量的分布。

表1 52名疑似肝炎患者的資料

表2 變量描述

表3 各指標(biāo)的統(tǒng)計(jì)描述±s)

表3表明,患者與非患者的體質(zhì)指數(shù)、酗酒量、谷丙轉(zhuǎn)氨酶的測(cè)量值是顯著不同的,經(jīng)檢驗(yàn)證實(shí),患者與非患者之間各指標(biāo)的差別均有統(tǒng)計(jì)學(xué)意義。患者的體質(zhì)指數(shù)、酗酒量、谷丙轉(zhuǎn)氨酶的測(cè)量值均高于非患者。對(duì)本實(shí)例,我們考慮幾種不同的情形深入分析。

1.不考慮協(xié)變量和數(shù)據(jù)間的相關(guān)性對(duì)肝炎診斷準(zhǔn)確性的影響。由R軟件繪制得到相應(yīng)的ROC曲線見圖1,計(jì)算AUC為0.75,可以認(rèn)為谷丙轉(zhuǎn)氨酶對(duì)肝炎的診斷的準(zhǔn)確度尚可。

圖1 未考慮協(xié)變量和隨機(jī)效應(yīng)影響的ROC曲線

表4 三種模型DIC及相關(guān)數(shù)值比較

表4列出了如上三個(gè)模型的DIC值,可以看出模型3的DIC值最小,所以選擇帶截距、時(shí)間項(xiàng)隨機(jī)效應(yīng)的ROC曲線回歸模型進(jìn)行下一步分析,即

ROCx(p)=Φ(γ1+γ2Φ-1(p)+β1x1+β2x2+λ0+λ1t)

該模型可以在WinBUGS軟件中用貝葉斯圖建模的方法表示成如圖2所示的有向關(guān)系圖。

觀察模擬結(jié)果迭代過程的蹤跡和自相關(guān)函數(shù)圖(如圖3、圖4),各參數(shù)的迭代軌跡趨于一條近似平穩(wěn)的水平線,滯后一次的自相關(guān)函數(shù)趨于0,這表明模擬結(jié)果的迭代過程基本收斂,得到的貝葉斯參數(shù)估計(jì)結(jié)果(表5)可以用于下一步的統(tǒng)計(jì)推斷。圖5的核密度函數(shù)圖給我們一些啟示,隨機(jī)效應(yīng)方差的核密度函數(shù)基本呈偏態(tài)分布,所以參數(shù)估計(jì)結(jié)果取中位數(shù)更好。故我們?nèi)≈形粩?shù)作為各參數(shù)的估計(jì)結(jié)果。

圖2 帶截距、時(shí)間項(xiàng)隨機(jī)效應(yīng)的ROC曲線 混合效應(yīng)模型有向圖

由表5可知,Φ-1(p) 和各參數(shù)的后驗(yàn)均數(shù)的 95%的置信區(qū)間都沒有包含0,表明各協(xié)變量對(duì) ALT

診斷肝膽疾病準(zhǔn)確度的影響均有統(tǒng)計(jì)學(xué)意義。不同患者間的隨機(jī)效應(yīng)方差為0.101,可以認(rèn)為不同患者的ALT對(duì)肝膽疾病診斷試驗(yàn)的準(zhǔn)確度不同,存在“個(gè)體效應(yīng)”。不同患者時(shí)間的隨機(jī)效應(yīng)方差也有統(tǒng)計(jì)學(xué)意義,方差為0.149,可以認(rèn)為不同患者的ALT對(duì)肝膽疾病的診斷準(zhǔn)確度存在著隨時(shí)間而變化的趨勢(shì)。估計(jì)方程結(jié)果顯示,在考慮了“個(gè)體效應(yīng)”后,體質(zhì)指數(shù)BMI值越小,每天的飲酒量越多,谷丙轉(zhuǎn)氨酶對(duì)肝膽疾病的診斷準(zhǔn)確性越高,也就是說,當(dāng)采用谷丙轉(zhuǎn)氨酶診斷肝膽疾病時(shí),區(qū)分非肥胖患者的準(zhǔn)確度高于肥胖患者的準(zhǔn)確度,區(qū)分酗酒患者的準(zhǔn)確度高于不酗酒患者的準(zhǔn)確度。

圖3 迭代10萬(wàn)次各參數(shù)迭代軌跡

圖4 各參數(shù)迭代10萬(wàn)次自相關(guān)函數(shù)圖

根據(jù)分層擬合的模型,可以得到不同肥胖程度不同酗酒量下,谷丙轉(zhuǎn)氨酶診斷肝膽疾病的ROC曲線。分別選取BMI值為23,26,29,酗酒量為0,5,10,20,計(jì)算相應(yīng)的AUC,見表6。

由表6可知,不同肥胖程度,不同酗酒量下的AUC不同。若不考慮協(xié)變量和隨機(jī)效應(yīng)對(duì)ROC曲線的影響,有可能導(dǎo)致ALT對(duì)于肥胖患者或者不嗜酒的患者的預(yù)測(cè)肝膽疾病的準(zhǔn)確性下降,給臨床疾病的確診和干預(yù)治療提供不良的診斷信息。所以基于廣義線性隨機(jī)效應(yīng)的ROC曲線模型可以有效地對(duì)重復(fù)測(cè)量試驗(yàn)診斷數(shù)據(jù)進(jìn)行刻畫和分析。

圖5 迭代10萬(wàn)次核密度函數(shù)圖

參數(shù)均值標(biāo)準(zhǔn)差MC誤差中位數(shù)95%置信區(qū)間gamma113660035000041360(0774,1993)gamma214360058000021436(1320,1553)beta109930049000010992(0897,1097)beta2-0181005000001-0180(-0281,-0079)sigma[1,1]01070032000010101(0065,0187)sigma[2,2]01560062000010149(0046,0299)

表6 不同肥胖程度,不同酗酒量下AUC

小 結(jié)

本文主要探討了針對(duì)重復(fù)診斷數(shù)據(jù)的ROC曲線廣義線性混合效應(yīng)模型的原理及其參數(shù)估計(jì)方法,并結(jié)合醫(yī)學(xué)實(shí)例,闡釋了該模型的軟件實(shí)現(xiàn)與應(yīng)用。重復(fù)測(cè)量診斷數(shù)據(jù)是醫(yī)藥衛(wèi)生研究中常見的數(shù)據(jù)形式,診斷試驗(yàn)的準(zhǔn)確度評(píng)價(jià)也是醫(yī)學(xué)研究不可或缺的部分。本文介紹的ROC曲線廣義線性混合效應(yīng)模型,是一種考慮協(xié)變量層次效應(yīng)對(duì)診斷試驗(yàn)準(zhǔn)確度影響的直接回歸模型,同時(shí)還度量了重復(fù)測(cè)量診斷數(shù)據(jù)間的相關(guān)性,建立了帶有隨機(jī)效應(yīng)的回歸模型。利用診斷準(zhǔn)確度的綜合指標(biāo)(ROC曲線下面積),對(duì)協(xié)變量影響和隨機(jī)效應(yīng)影響下的診斷試驗(yàn)準(zhǔn)確度進(jìn)行合理評(píng)價(jià),綜合性強(qiáng),準(zhǔn)確性高。在實(shí)際應(yīng)用中,該方法建模靈活,解釋合理,通過WinBUGS軟件編程可方便實(shí)現(xiàn),不僅可以用于診斷試驗(yàn)結(jié)果為連續(xù)性的重復(fù)測(cè)量數(shù)據(jù)和有序多分類的重復(fù)測(cè)量數(shù)據(jù)的分析,而且對(duì)于診斷試驗(yàn)結(jié)果的分布不作特定要求,根據(jù)協(xié)變量和隨機(jī)效應(yīng)對(duì)診斷試驗(yàn)準(zhǔn)確度綜合指標(biāo)的影響,繪制不同條件下的ROC曲線并計(jì)算其曲線下面積,從而形象直觀地評(píng)價(jià)診斷試驗(yàn)的準(zhǔn)確度,為臨床決策提供理論依據(jù)。

[1]李康,魏韋.醫(yī)學(xué)診斷試驗(yàn)評(píng)價(jià)的ROC分析——重復(fù)測(cè)量診斷數(shù)據(jù)的ROC曲線.中國(guó)醫(yī)院統(tǒng)計(jì),2001,8(1):3-6.

[2]Zhou XH著.宇傳華譯.診斷醫(yī)學(xué)統(tǒng)計(jì)學(xué).北京:人民衛(wèi)生出版社,2005.

[3]尉潔,宋嬌嬌,劉桂芬.基于貝葉斯估計(jì)的診斷試驗(yàn)ROC曲線回歸模型.中國(guó)衛(wèi)生統(tǒng)計(jì),2010,27(2):152-154.

[4]Zhao Y,Staudenmayer J,Coull BA,et al.General Design Bayesian Generalized Linear Mixed Models.Statistical Science,2006,21(1):35-51.

[5]O′Malley AJ,Zou KH.Bayesian multivariate hierarchical transformation models for ROC analysis.Stat Med,2006,25(3):459-479.

(責(zé)任編輯:劉 壯)

ROC Curve Based on Generalized Linear Mixed Effects Models in Repeated Diagnostic Tests

Ma Chuntao,Xiong Wei,Tian Maozai

(CenterofStatistics,theSchoolofStatistics,RenminUniversityofChina(100872),Beijing)

Objective To investigate the impact of covariates on diagnostic test and assess the correlation between repeated measurement data,this paper explores innovative modeling techniques of ROC curve.Methods We introduce the new ROC curve method based on generalized linear mixed effects model and apply Bayesian techniques to parameters estimation with Winbugs Software.Further,areas under the ROC curve(AUC) with different values of covariates could be calculated in terms of assessment.Results Cases analysis results indicate the proposed method efficiently explores the repeated measurement data and provides parameters with practical significance,serving as a golden reference.Conclusion The ROC curve based on generalized linear mixed effects models can be effectively used to solve the test accuracy evaluation problem of the repeated diagnostic trials.

ROC curve;Mixed effects;Bayesian;Repeated measurement;Diagnostic tests

*:教育部哲學(xué)社會(huì)科學(xué)研究重大課題攻關(guān)項(xiàng)目(No.15JZD015),國(guó)家自然科學(xué)基金(No.11271368),北京市社會(huì)科學(xué)基金重大項(xiàng)目(No.15ZDA17),教育部高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金(No.20130004110007),國(guó)家社會(huì)科學(xué)基金重點(diǎn)項(xiàng)目(No.13AZD064),教育部人文社會(huì)科學(xué)重點(diǎn)研究基地重大項(xiàng)目(15JJD910001),北京市社科聯(lián)項(xiàng)目決策咨詢項(xiàng)目(No.2016010021),中國(guó)人民大學(xué)科學(xué)研究基金(中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助)項(xiàng)目成果(No.15XNL008)

1.中國(guó)人民大學(xué)應(yīng)用統(tǒng)計(jì)科學(xué)研究中心,中國(guó)人民大學(xué)統(tǒng)計(jì)學(xué)院(100872)

2.對(duì)外經(jīng)濟(jì)貿(mào)易大學(xué)大數(shù)據(jù)與風(fēng)險(xiǎn)管理研究中心,統(tǒng)計(jì)學(xué)院

3.蘭州財(cái)經(jīng)大學(xué)統(tǒng)計(jì)學(xué)院

4.新疆財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與信息學(xué)院

△通信作者: 田茂再,E-mail:mztian@ruc.edu.cn

猜你喜歡
效應(yīng)測(cè)量模型
一半模型
鈾對(duì)大型溞的急性毒性效應(yīng)
懶馬效應(yīng)
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
把握四個(gè)“三” 測(cè)量變簡(jiǎn)單
滑動(dòng)摩擦力的測(cè)量和計(jì)算
滑動(dòng)摩擦力的測(cè)量與計(jì)算
應(yīng)變效應(yīng)及其應(yīng)用
3D打印中的模型分割與打包
主站蜘蛛池模板: 成人字幕网视频在线观看| 992Tv视频国产精品| 久久久久九九精品影院| 国内精品久久九九国产精品| 日韩在线播放中文字幕| 久久一级电影| 成人精品视频一区二区在线| 久久精品亚洲热综合一区二区| 在线精品亚洲一区二区古装| 亚洲中文字幕国产av| 国产精品久久久久婷婷五月| 亚洲无码视频喷水| 久久香蕉国产线看观看精品蕉| 精品一区二区三区四区五区| 狠狠亚洲五月天| 她的性爱视频| 成人在线视频一区| 国产流白浆视频| 国产精品页| 日韩成人免费网站| 视频二区国产精品职场同事| 国产一级α片| 秋霞国产在线| 91青青草视频| 国产一区二区色淫影院| 国产色伊人| 午夜不卡福利| 亚洲男人的天堂网| 亚洲中文久久精品无玛| 在线播放91| 91亚洲免费视频| 99久久精品免费观看国产| 一级香蕉人体视频| 精品综合久久久久久97超人该| 婷婷五月在线| 国产激情无码一区二区APP| 国产精品免费p区| 成年A级毛片| 日韩AV无码免费一二三区| 欧美伊人色综合久久天天| 亚洲男人天堂久久| 国产后式a一视频| 亚洲精品动漫| 亚洲男人在线天堂| 一级成人欧美一区在线观看| 内射人妻无套中出无码| 国产一级毛片网站| 亚瑟天堂久久一区二区影院| 麻豆精品视频在线原创| 日本一区二区不卡视频| 亚洲国产精品无码AV| 91探花国产综合在线精品| 国产亚洲精品97AA片在线播放| 热久久国产| 国产区福利小视频在线观看尤物| 国产精品无码久久久久AV| 久久人搡人人玩人妻精品一| 二级特黄绝大片免费视频大片| 亚洲国产欧美目韩成人综合| 亚洲综合极品香蕉久久网| 国产福利免费在线观看| 亚洲色欲色欲www网| 国产91视频免费观看| 毛片免费试看| 露脸国产精品自产在线播| 免费毛片a| 亚洲AV无码不卡无码| 亚洲中文无码av永久伊人| 四虎影视无码永久免费观看| 国产日韩久久久久无码精品 | 欧美日韩一区二区在线播放| 经典三级久久| 亚洲aaa视频| 四虎影院国产| 婷婷丁香色| 丝袜亚洲综合| 无码中文AⅤ在线观看| 四虎亚洲精品| 又爽又大又黄a级毛片在线视频| 亚洲国产成人久久77| 77777亚洲午夜久久多人| 欧美成人午夜在线全部免费|