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

        ?

        樣本量估計及其在nQuery和SAS軟件上的實現
        ——均數比較(九)*

        2019-05-24 01:38:52南方醫(yī)科大學公共衛(wèi)生學院生物統計學系510515錢若遠陳平雁
        中國衛(wèi)生統計 2019年2期
        關鍵詞:參數設置樣本量協方差

        南方醫(yī)科大學公共衛(wèi)生學院生物統計學系(510515) 王 巖 錢若遠 陳平雁

        1.3 多樣本的均數比較

        1.3.1 差異性檢驗

        1.3.1.4 協方差分析(ANCOVA)

        (1-43)

        計算樣本量時,先設定樣本量初始值,然后迭代樣本量直到所得的檢驗效能滿足設定值為止,此時的樣本量即研究所需的樣本量[2]。

        【例1-27】某研究欲比較三種不同教學干預方法(有聲思維教學法、指導性閱讀思維教學法、指導性閱讀教學法)對小學四年級學生閱讀理解能力的影響。采用平衡設計,在干預前后分別評估受試者的錯誤檢測任務(error detection task,EDT)得分。以干預前EDT得分和干預后理解監(jiān)控問卷(comprehension monitoring questionnaire)得分為2個協變量,干預后EDT得分為因變量進行協方差分析。根據既往研究數據,三組均值分別為8.2220,9.8148,6.1904,標準差為2.3788,ρ2為0.4434,設定檢驗水準為0.05,試以檢驗效能分別為0.8和0.9進行樣本量估計。

        nQuery Advanced實現:設定檢驗水準α=0.05;檢驗效能分別取80%和90%。依據上述基礎數據可知,G=3,c=2,ρ2=0.4344,σ=2.3788,r1=r2=r3=1,μ={μ1,μ2,μ3}={8.2220,9.8148,6.1904}計算可得V=2.20,在nQueryAdvanced 主菜單選擇:

        Design:⊙Fixed Term

        Goal:⊙Means

        No.of Groups:⊙>Two

        Analysis Method:⊙Test

        方法框中選擇Analysis of Covariance(ANCOVA)。

        在彈出的樣本量計算窗口將各參數鍵入,如圖1-65所示,檢驗效能取80%的結果N=18,即總樣本量為18例,每組6例。檢驗效能取90%的結果N=22,即總樣本量為22例,三組例數分別為8例、7例和7例(因總例數不是組數的整倍數,故各組例數不絕對相等)。

        圖1-65 nQuery Advanced 關于例1-27樣本量估計的參數設置與計算結果

        SAS 9.4軟件實現:

        %let u={8.2220 9.8148 6.1904}; /*各組均值*/

        %let r={1 1 1}; /*各組樣本量與第一組(參照組)樣本量的比值*/

        proc IML;

        start MGT4_1(a,G,sd,c,R2,power);

        error=0;

        if ( a>=1 | a<=0 ) then do; error=1; print "Error" "Test significance Level must be in 0-1"; end;

        if ( sd<=0 ) then do; error=1; print "Error" "Common standard deviation must be >=0"; end;

        if (G<=0 | ceil(G)^=G) then do; error=1; print "Error" "The Number of groups must be positive integer "; end;

        if (c<=0 | ceil(c)^=c) then do; error=1; print "Error" "The Number of coveriates must be positive integer "; end;

        if(power>=100 | power<1) then do; error=1; print "error" "Power(%) must be in 1-100"; end;

        if(R2>1 | R2<0) then do; error=1; print "error" "R_squared value must be in 0-1";end;

        if(error=1) then stop;

        if(error=0) then do;

        n1=ceil(c/G+1); n=n1#&r.;sum_N=n[1,+]-0.01;

        do until(pw>=power/100);

        sum_n=ceil(sum_n+0.01); n1=sum_N/&r.[1,+]; n=n1#&r.;

        u_mean=(n#&u.)[,+]/sum_N;

        v=(&r.#(&u.-u_mean)##2)[,+]/&r.[,+];

        n_mean=sum_N/G;

        sigma_e2=(1-R2)*sd##2;

        lamda=n_mean*G*v/sigma_e2;

        df1=G-1;df2=sum_N-G-c;

        f=finv(1-a,df1,df2);

        pw=1-probf(f,df1,df2,lamda); pw1=100*pw;

        end;

        end;

        print a[label="Test Significance level"]

        G[label="Number of groups"]

        V[label="Variance of means"]

        sd[label="Common standard deviation"]

        c[label="Number of coveriates"]

        R2[label="R_squared value between the reponse and the covariates"]

        Pw1[label="Power(%)"]

        sum_N[label="Total sample size"];

        finish MGT4_1;

        run MGT4_1(0.05,3,2.3788,2,0.4434,80);

        run MGT4_1(0.05,3,2.3788,2,0.4434,90);

        quit;

        SAS 9.4運行結果:

        圖1-66a SAS 9.4關于例1-27樣本量估計的參數設置與計算結果(power=80%)

        圖1-66b SAS 9.4關于例1-27樣本量估計的參數設置與計算結果(power=90%)

        1.3.1.5 多變量方差分析(MANOVA)

        方法:多變量方差分析針對分析中有多個反應變量的情況,目前有多種常用統計量可用于分析,nQuery給出了其中三種供用戶選擇,分別是Wilks’似然比統計量、Pillai-Bartlett Trace統計量、Hotelling-Lawley Trace統計量。Muller and Barton (1989)[3],以及 Muller,LaVange,Ramey(1992)[4]等給出了多變量方差分析的樣本量和檢驗效能的估計方法,各主效應及其交互效應的檢驗效能估計是建立在各自的自由度及非中心參數的F分布上。其檢驗效能的計算公式為:

        1-β=1-ProbF(F1-α,df1,df2,df1,df2,λ)

        (1-44)

        式中,df1,df2代表F分布的自由度,λ為非中心參數,根據不同的統計量,其各參數計算方式如下:

        首先,我們用p代表因變量個數,q代表所研究的分組變量的水平數,X代表設計矩陣,r代表設計矩陣的秩,M代表均值矩陣,Σ代表協方差矩陣,C代表對比矩陣,不同效應的檢驗可以構建不同對比矩陣C,N代表總樣本量。Wilks’似然比統計量、Pillai-Bartlett Trace統計量、Hotelling-Lawley Trace統計量可基于如下矩陣構建 :

        H=(CM)′[C(X′X)-1C′]-1(CM)

        (1-45)

        E=Σ(N-r)

        (1-46)

        T=H+E

        (1-47)

        (1)Wilks’ Lambda

        利用矩陣(1-46), (1-47)可得Wilks’似然比統計量:

        W=|ET-1|

        (1-48)

        該檢驗統計量轉化為近似的F統計量為:

        (1-49)

        其中,

        df1=ap,

        df2=g[(N-r)-(p-a+1)/2]-(ap-2)/2,

        (2) Pillai-Bartlett Trace

        利用矩陣(1-45),(1-47)可得Pillai-Bartlett Trace統計量:

        PBT=tr(HT-1)

        (1-50)

        該檢驗統計量轉化為近似的F統計量為:

        (1-51)

        其中,

        df1=ap,

        df2=s[(N-r)-p+s],

        (3)Hotelling-Lawley Trace

        利用矩陣(1-45), (1-46)可得Hotelling-Lawley Trace統計量:

        HLT=tr(HE-1)

        (1-52)

        該檢驗統計量轉化為近似的F統計量為:

        (1-53)

        其中,

        df1=ap,

        df2=s[(N-r)-p-1]+2,

        對于上述各統計量,其非中心參數λ為:

        λ=F·df1

        (1-54)

        計算樣本量時,先設定樣本量初始值,然后選擇相應統計方法并迭代樣本量直到所得的檢驗效能滿足設定值為止,此時的樣本量即研究所需的樣本量[2]。

        【例1-28】某研究欲比較不同教學方法在培養(yǎng)中學生良好學習行為中的作用,教學方法分為傳統教學組、視聽教學組和對照組三個組,因變量包括4個衡量學習行為的指標:學習環(huán)境、學習習慣、筆記能力和總結能力。根據既往研究,均值矩陣如矩陣M所示,協方差矩陣如矩陣Σ所示,設定檢驗水準為0.05,檢驗效能為80%,試采用Wilks’ Lambda方法據此參數進行樣本量估計。

        nQuery Advanced實現:設定檢驗水準為α= 0.05;檢驗效能取80%。依據上述基礎數據可知,p=4,試驗包含1個影響因素,其水平數為3,均值矩陣如矩陣M所示,協方差矩陣Σ如矩陣所示。在nQueryAdvanced 主菜單選擇:

        Design:⊙Fixed Term

        Goal:⊙Means

        No.of Groups:⊙>Two

        Analysis Method:⊙Test

        方法框中選擇Multivariate Analysis of Variance(MANOVA)。

        在彈出的樣本量計算窗口將各參數鍵入,如圖1-67a所示,均值矩陣和協方差矩陣如圖1-67b和圖1-67c所示,結果n和N分別為37和111,即每組樣本量為37例,總樣本量為111例。

        圖1-67a nQuery Advanced 關于例1-28樣本量估計的參數設置與計算結果

        圖1-67b nQuery Advanced 關于例1-28樣本量估計的參數設置與計算結果(均值矩陣)

        圖1-67c nQueryAdvanced 關于例1-28樣本量估計的參數設置與計算結果(協方差矩陣)

        SAS 9.4軟件實現:

        %let f_level={3}; /*各分組變量的水平數*/

        %let m= { 8.75 7.79 8.14,8.10 7.71 7.19,17.83 17.32 16.67,18.90 19.18 17.69};

        /*因變量在所有分組變量各個水平組合下的均值*/

        prociml;

        start MGT3(alpha,power,p,f,rho,sd,sigma_def,method);

        f_level=&f_level.;

        m=t(&m.);

        /*parameter check*/

        error=0;

        if ( alpha>=1 | alpha<=0 ) then do; error=1; print "Error" "Test significance Level must be in 0-1"; end;

        /*參數f代表分組變量個數,不超過3個*/

        if (f^=1 & f^=2 & f^=3) then do; error=1; print "Error""The number of factors must be in 1,2,3";end;

        /*參數power代表各個效應的檢驗效能向量;當研究包含A和B兩個分組變量時,可依次定義A、B及其交互效應AB的檢驗效能;當研究包含A、B、C三個分組變量時,可依次定義A、B、C、AB、AC、BC、ABC的檢驗效能*/

        if f=1 & nrow(power)^=1 & ncol(power)^=1 then do;error=1;print "Error""Please input power of factor A";end;

        if f=2 & nrow(power)^=3 & ncol(power)^=3 then do;error=1;print "Error""Please input power of factor A,factor B,factor AB (input "." if missing)";end;

        if f=3 & nrow(power)^=7 & ncol(power)^=7 then do; error=1; print "Error""Please input power of factor A,factor B,factor C,factor AB,factor AC,factor BC,factor ABC (input "." if missing)"; end;

        if (p<=0 | ceil(p)^=p) then do; error=1; print "Error" "The number of response variables must be a positive integer "; end;

        if error=0 then do; q=1;do i=1 to f; q= q * &f_level.[i]; end;end;else stop;

        if nrow(f_level)^=f then do;error=1; print "Error""The row number of factor levels should be equal to the number of factors";end;

        if ncol(m)^=p then do; error=1; print "Error" "The row number of mean matrix should be equal to the number of response variables"; end;

        if nrow(m)^=q then do; error=1; print "Error" "The column number of mean matrix should be equal to the product of factor levels"; end;

        /*協方差矩陣可以由參數rho和sd生成,也可以由參數sigma_def直接定義*/

        if (rho^=.& (rho>=1 | rho<=0)) then do; error=1;print "error" "Correlation must be in 0-1";end;

        if (sd^=.& ( sd<=0 )) then do;error=1; print "Error" "Standard deviation at each level must be >=0"; end;

        if (sigma_def^=.& (ncol(sigma_def)^=p | nrow(sigma_def)^=p)) then do;

        error=1; print "Error" "The row and column numbers of covariance matrix should be equal to the number of response variables"; end;

        /*參數method代表所選用的檢驗統計量,其中1=Wilks' Lambda統計量,2=Pillai-Bartlett Trace統計量,3=Hotelling-Lawley Trace統計量。*/

        if (method^=1 & method^=2 & method^=3) then do; error=1; print "Error""method must be in 1,2,3";end;

        if (rho^=.& sd^=.) then do;

        covariance=rho*sd**2;sigma=j(p,p,covariance);

        do i=1 to p; sigma[i,i]=sd**2; end;

        end;

        else sigma=sigma_def;

        if(error=1) then stop;

        if(error=0)then do;

        /*C matrix*/

        max_level=&f_level.[<>,];

        origin_c=j(max_level-1,max_level,0);

        origin_j=j(max_level-1,max_level,.);

        do i=2 to max_level;

        e1=-(i-1)/sqrt(i*(i-1)); e2=1/sqrt(i*(i-1)); e3=j(1,i-1,e2); e4=1/sqrt(i);

        e5=j(1,i,e4);

        index=max_level+1-i;

        origin_c[index,index]=e1;

        origin_c[index,index+1:max_level]=e3;

        origin_j[index,index:max_level]=e5;

        end;

        f1=&f_level.[1]; index1=max_level-f1+1;

        C1=origin_c[index1:max_level-1,index1:max_level];

        J1=origin_j[index1,index1:max_level];C_A=C1; power_a=power[1];

        if f>1 then do;

        f2=&f_level.[2]; index2=max_level-f2+1; C2=origin_c[index2:max_level-1,index2:max_level];

        J2=origin_j[index2,index2:max_level];

        C_A=C1@J2;C_B=J1@C2;C_AB=C1@C2;

        power_b=power[2];power_ab=power[3];

        end;

        if f>2 then do;

        f3=&f_level.[3];

        index3=max_level-f3+1;

        C3=origin_c[index3:max_level-1,index3:max_level];

        J3=origin_j[index3,index3:max_level];

        C_A=C1@J2@J3;C_B=J1@C2@J3;C_C=J1@J2@C3;

        C_AB=C1@C2@J3;C_AC=C1@J2@C3;C_BC=J1@C2@C3;

        C_ABC=C1@C2@C3;

        power_c=power[3];

        power_ab=power[4];power_ac=power[5];power_bc=power[6];

        power_abc=power[7];

        end;

        /*power*/

        %macro pw(c,pw_exp);

        n_orig=j(1,q,1);n=j(1,q,1);

        if &pw_exp.^=.then do;

        if method=1 then do;

        test="Wilks' Lambda";

        do until (pw>=&pw_exp.);

        n=n+n_orig;xx=t(n)#I(q);E=sigma * (n[<>,+]-q);

        theta=&C.* m; H=t(theta)*inv((&C.*inv(XX)*t(&C.)))*theta;

        T=H+E;df1=a*p;W=det(E*inv(T));fmm=a**2+p**2-5;

        if fmm>0 then g=sqrt((a**2*p**2-4)/(a**2+p**2-5)); else g=1;

        eta=1-W**inv(g);

        df2=g*((n[<>,+]-q)-(p-a+1)/2)-(a*p-2)/2;

        if df2>0 then do;

        F_statistics=(eta/df1)/((1-eta)/df2);

        f_c=finv(1-alpha,df1,df2);lambda=df1*F_statistics;

        pw=(1-probf(f_c,df1,df2,lambda))*100; sum_n=n[<>,+];

        end;

        else pw=0;

        end;

        end;

        if method=2 then do;

        test="Pillai-Bartlett Trace";

        do until (pw>=&pw_exp.);

        n=n+n_orig;xx=t(n)#I(q);E=sigma * (n[<>,+]-q);

        theta=&C.* m; H=t(theta)*inv((&C.*inv(XX)*t(&C.)))*theta;

        T=H+E;df1=a*p;

        ht=h*inv(t);PBT=trace(ht);s=min(a,p);

        eta=pbt/s;df2=s*((n[<>,+]-q)-p+s);

        if df2>0 then do;

        F_statistics=(eta/df1)/((1-eta)/df2);

        f_c=finv(1-alpha,df1,df2);lambda=df1*F_statistics;

        pw=(1-probf(f_c,df1,df2,lambda))*100; sum_n=n[<>,+];

        end;

        else pw=0;

        end;

        end;

        if method=3 then do;

        test="Hotelling-Lawley Trace";

        do until (pw>=&pw_exp.);

        n=n+n_orig;xx=t(n)#I(q);E=sigma * (n[<>,+]-q);

        theta=&C.* m; H=t(theta)*inv((&C.*inv(XX)*t(&C.)))*theta;

        T=H+E;df1=a*p;

        he=h*inv(E);hlt=trace(he);

        s=min(a,p);eta=(hlt/s)/(1+hlt/s);

        df2=s*((n[<>,+]-q)-p-1)+2;

        if df2>0 then do;

        F_statistics=(eta/df1)/((1-eta)/df2);

        f_c=finv(1-alpha,df1,df2);lambda=df1*F_statistics;

        pw=(1-probf(f_c,df1,df2,lambda))*100; sum_n=n[<>,+];

        end;

        else pw=0;

        end;

        end;

        end;

        else do; pw=.; sum_n=0; end;

        %mend;

        /*compute total sample size*/

        a=&f_level.[1]-1; %pw(C_A,power_a);sum_n_a=sum_n;Smp_size=sum_n_a;

        if f>1 then do;

        a=&f_level.[2]-1;%pw(C_B,power_b);sum_n_b=sum_n;

        a=(&f_level.[1]-1)*(&f_level.[2]-1);%pw(C_AB,power_ab);

        sum_n_ab=sum_n;

        Smp_size=j(3,1,.);

        Smp_size[1]=sum_n_a;

        Smp_size[2]=sum_n_b;

        Smp_size[3]=sum_n_ab;

        end;

        if f>2 then do;

        a=&f_level.[3]-1;%pw(C_C,power_c);pw_c=pw;sum_n_c=sum_n;

        a=(&f_level.[1]-1)*(&f_level.[3]-1);

        %pw(C_AC,power_ac);sum_n_ac=sum_n;

        a=(&f_level.[2]-1)*(&f_level.[3]-1);

        %pw(C_BC,power_bc);sum_n_bc=sum_n;

        a=(&f_level.[1]-1)*(&f_level.[2]-1)*(&f_level.[3]-1);

        %pw(C_ABC,power_abc);sum_n_abc=sum_n;

        Smp_size=j(7,1,.);

        Smp_size[1]=sum_n_a;

        Smp_size[2]=sum_n_b;

        Smp_size[3]=sum_n;

        Smp_size[4]=sum_n_ab;

        Smp_size[5]=sum_n_ac;

        Smp_size[6]=sum_n_bc;

        Smp_size[7]=sum_n_abc;

        end;

        Total_n=Smp_size[<>,];grp_size=Total_n/q;

        /*compute power for total sample size*/

        %macro pw2(c,pw_exp);

        if &pw_exp.^=.then do;

        n=j(1,q,grp_size);xx=t(n)#I(q);E=sigma * (n[<>,+]-q);

        theta=&C.* m; H=t(theta)*inv((&C.*inv(XX)*t(&C.)))*theta;

        T=H+E;df1=a*p;

        if method=1 then do;

        W=det(E*inv(T));fmm=a**2+p**2-5;

        if fmm>0 then g=sqrt((a**2*p**2-4)/(a**2+p**2-5)); else g=1;

        eta=1-W**inv(g);

        df2=g*((n[<>,+]-q)-(p-a+1)/2)-(a*p-2)/2;

        F_statistics=(eta/df1)/((1-eta)/df2);

        end;

        if method=2 then do;

        ht=h*inv(t);PBT=trace(ht);s=min(a,p);

        eta=pbt/s;df2=s*((n[<>,+]-q)-p+s);

        F_statistics=(eta/df1)/((1-eta)/df2);

        end;

        if method=3 then do;

        he=h*inv(E);hlt=trace(he);s=min(a,p);

        eta=(hlt/s)/(1+hlt/s);df2=s*((n[<>,+]-q)-p-1)+2;

        F_statistics=(eta/df1)/((1-eta)/df2);

        end;

        f_c=finv(1-alpha,df1,df2);lambda=df1*F_statistics;

        pw=(1-probf(f_c,df1,df2,lambda))*100;

        end;

        else pw=.;

        %mend;

        Level= &f_Level.;

        factor="Factor A";

        a=&f_level.[1]-1; %pw2(C_A,power_a); pw_a=pw;power=pw_a;Alpha1=alpha;

        if f>1 then do;

        factor={"Factor A","FactorB","Factor AB"};

        a=&f_level.[2]-1; %pw2(C_B,power_b); pw_b=pw;

        a=(&f_level.[1]-1)*(&f_level.[2]-1); %pw2(C_AB,power_ab); pw_ab=pw;

        power=j(3,1,.);power[1]=pw_a; power[2]=pw_b;power[3]=pw_ab;

        Alpha1=j(3,1,alpha);

        end;

        if f>2 then do;

        factor={"Factor A","FactorB","FactorC","FactorAB","FactorAC","FactorBC","Factor ABC"};

        a=&f_level.[3]-1; %pw2(C_C,power_c); pw_c=pw;

        a=(&f_level.[1]-1)*(&f_level.[3]-1); %pw2(C_AC,power_ac); pw_ac=pw;

        a=(&f_level.[2]-1)*(&f_level.[3]-1); %pw2(C_BC,power_bc); pw_bc=pw;

        a=(&f_level.[1]-1)*(&f_level.[2]-1)*(&f_level.[3]-1);

        %pw2(C_ABC,power_abc);pw_abc=pw;

        power=j(7,1,.);power[1]=pw_a; power[2]=pw_b;power[3]=pw_c;

        power[4]=pw_ab;power[5]=pw_ac;power[6]=pw_bc;power[7]=pw_abc;

        Alpha1=j(7,1,alpha);

        end;

        Mean_Matrix=&m.;

        end;

        print test[Label="Test"]

        p[Label="Number of Response Variables"]

        sd[Label="Common Standard Deviation"]

        rho[Label="Correlation"]

        grp_size[Label="Group Size"]

        Total_n[Label="Total Sample Size"];

        Print Factor [label="Factor"]

        Level [Label="Level"]

        Alpha1 [Label="Alpha"]

        Power [Label="Power(%)"];

        Print Mean_Matrix;

        Print sigma[label="Covariance Matrix"];

        finish MGT3;

        %let sigma1={3.641 1.274 2.641 4.555,1.274 2.623 1.947 2.722,2.641 1.947 9.548 7.001,4.555 2.722 7.001 15.914};

        %let power={80};

        run MGT3(0.05,&power.,4,1,.,.,&sigma1.,1);

        quit;

        SAS 9.4運行結果:

        圖1-68 SAS 9.4關于例1-28樣本量估計的參數設置與計算結果

        猜你喜歡
        參數設置樣本量協方差
        醫(yī)學研究中樣本量的選擇
        內蒙古統計(2021年4期)2021-12-06 02:49:20
        航空裝備測試性試驗樣本量確定方法
        測控技術(2018年4期)2018-11-25 09:46:52
        Sample Size Calculations for Comparing Groups with Binary Outcomes
        蟻群算法求解TSP中的參數設置
        不確定系統改進的魯棒協方差交叉融合穩(wěn)態(tài)Kalman預報器
        自動化學報(2016年8期)2016-04-16 03:38:55
        一種基于廣義協方差矩陣的欠定盲辨識方法
        動車環(huán)境下U900異頻切換參數設置探討
        基于MATLAB仿真的井下變壓器參數設置研究
        縱向數據分析中使用滑動平均Cholesky分解對回歸均值和協方差矩陣進行同時半參數建模
        關于協方差的U統計量檢驗法
        日本一区二区三区一级免费| 欧美gv在线观看| 久久麻豆精品国产99国产精| 蜜臀aⅴ永久无码一区二区| 久久精见国产亚洲av高清热| 国产高清在线观看av片| 国产精品亚洲二区在线观看| 中国精品视频一区二区三区| 久久久国产精品三级av| 无套熟女av呻吟在线观看| 国产精品对白刺激久久久| 最新精品国偷自产在线婷婷| 日本女优爱爱中文字幕| 激情综合婷婷色五月蜜桃| 成人免费看吃奶视频网站| 国产精品视频免费的| 超碰青青草手机在线免费观看| 999精品无码a片在线1级| 999国内精品永久免费视频| 高清国产美女av一区二区| 美女被内射很爽的视频网站| 天堂在线资源中文在线8| 中文字幕熟妇人妻在线视频| 精品一区二区三区影片| 男女激情视频网站免费在线| 东京热无码av一区二区| 久久中文字幕av一区二区不卡| 日本国主产一区二区三区在线观看| 中文字幕亚洲乱码熟女1区| 久久99久久99精品中文字幕| 五月婷婷激情六月| 国产av一啪一区二区| 久久人人爽爽爽人久久久| 乱子真实露脸刺激对白| 大香蕉久久精品一区二区字幕| 日本免费大片一区二区| 男同gay毛片免费可播放| 99久久精品一区二区三区蜜臀| 青青草手机视频免费在线播放| 中文字幕乱码无码人妻系列蜜桃| 亚洲专区欧美|