background

เตรียมสอบปลายภาควิชา Data Analysis

สุทธิศักดิ์ เด่นดวงใจ
สุทธิศักดิ์ เด่นดวงใจ
แพทย์ประจำบ้าน สาขาอาชีวเวชศาสตร์
นิสิต ปริญญาโท สาขาการวิจัยและการจัดการด้านสุขภาพ

3. ค่า hsCRP ที่ได้มีความแตกต่างระหว่างเพศหรือไม่

พิจารณาการกระจายของข้อมูลก่อนว่า Normal Distribution หรือไม่

. histogram hsCRP, normal
(bin=36, start=.156, width=.25566667)

. graph box hsCRP, horizon

จะได้ว่าเบ้ขวา เราอาจจะต้อง Transform หรือ ใช้ Wilcoxon ranked-sum test แทน

ที่นี้หากไม่รู้ว่าจะ Transfrom ด้วย Fuction อะไร สามารถใช้ gladder เพื่อสร้างกราฟมาดูได้

. gladder hsCRP

ใช้ Independent samples T-Test

สร้างตัวแปรค่า Log ฐาน e และตรวจสอบการกระจายจะพบว่ากระจายเป็นปกติมากขึ้น

. gen log_hsCRP = log(hsCRP)

. histogram log_hsCRP , normal
(bin=36, start=-1.8578993, width=.11373179)

. graph box log_hsCRP, horizontal

จากนั้นทดสอบ Assumption ว่ามันกระจายสองกลุ่มเท่ากันไหมด้วย sdtest จากนั้นจึงใช้ ttest

. sdtest log_hsCRP,by(sex)

Variance ratio test
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
    Male |   2,038    .1035676    .0209787    .9470669    .0624256    .1447095
  Female |   2,462    .1032839    .0234997    1.166021    .0572027    .1493652
---------+--------------------------------------------------------------------
combined |   4,500    .1034124    .0159849    1.072297    .0720742    .1347506
------------------------------------------------------------------------------
    ratio = sd(Male) / sd(Female)                                 f =   0.6597
Ho: ratio = 1                                  degrees of freedom = 2037, 2461

    Ha: ratio < 1               Ha: ratio != 1                 Ha: ratio > 1
  Pr(F < f) = 0.0000         2*Pr(F < f) = 0.0000           Pr(F > f) = 1.0000

. ttest log_hsCRP,by(sex) unequal

Two-sample t test with unequal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. Err.   Std. Dev.   [95% Conf. Interval]
---------+--------------------------------------------------------------------
    Male |   2,038    .1035676    .0209787    .9470669    .0624256    .1447095
  Female |   2,462    .1032839    .0234997    1.166021    .0572027    .1493652
---------+--------------------------------------------------------------------
combined |   4,500    .1034124    .0159849    1.072297    .0720742    .1347506
---------+--------------------------------------------------------------------
    diff |            .0002837    .0315015               -.0614747     .062042
------------------------------------------------------------------------------
    diff = mean(Male) - mean(Female)                              t =   0.0090
Ho: diff = 0                     Satterthwaite's degrees of freedom =  4496.41

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.5036         Pr(|T| > |t|) = 0.9928          Pr(T > t) = 0.4964

จะพบว่าระดับ hsCRP ไม่มีความสัมพันธ์กับ เพศ

ใช้ Wilcoxon ranked-sum test

. ranksum hsCRP, by(sex)

Two-sample Wilcoxon rank-sum (Mann-Whitney) test

         sex |      obs    rank sum    expected
-------------+---------------------------------
        Male |     2038     4611317     4586519
      Female |     2462     5515933     5540731
-------------+---------------------------------
    combined |     4500    10127250    10127250

unadjusted variance   1.882e+09
adjustment for ties  -219482.36
                     ----------
adjusted variance     1.882e+09

Ho: hsCRP(sex==Male) = hsCRP(sex==Female)
             z =   0.572
    Prob > |z| =   0.5676

จะพบว่าระดับ hsCRP ไม่มีความสัมพันธ์กับ เช่นเดียวกัน

4. ค่า hsCRP ที่ได้มีความสัมพันธ์กับอายุหรือไม่

ทดสอบโดยใช้ Pearson's Correlation (ใส่ ,sig เพื่อให้ทดสอบ Significant ด้วย)

. pwcorr hsCRP age, sig

             |    hsCRP      age
-------------+------------------
       hsCRP |   1.0000 
             |
             |
         age |   0.1894   1.0000 
             |   0.0000
             |

จะเห็นว่า hsCRP มีความสัมพันธ์เชิงเส้นตรงกับ age ในระดับต่ำ อย่างมีนัยสำคัญทางสถิติ (pearson's correlation : 0.1894)

5. จำนวนคนที่มี PA_level ระดับ moderate และ high ในกลุ่มที่มี hsCRP < 1.0

โดยกำหนดให้

  • Low (0) คือ hsCRP < 1.0
  • Moderate (1) คือ hsCRP 1.0 - 3.0
  • High (2) คือ hsCRP > 3.0 แบบนี้ถ้าต้องการตอบเพียงคำถามนี้ เราใช้แค่ .. if .. ก็ได้
. tab PA_level if hsCRP <1

   PA_level |      Freq.     Percent        Cum.
------------+-----------------------------------
        Low |      1,020       45.19       45.19
   Moderate |        672       29.77       74.97
       High |        565       25.03      100.00
------------+-----------------------------------
      Total |      2,257      100.00

ก็ได้คำตอบแล้วว่า Moderate มี 672 คน High มี 565 คน

แต่ในเมื่อเราอาจจะต้องสร้างตัวแปร hsCRP_level แบ่งเป็น 3 กลุ่ม เพื่อใช้ในข้อต่อไป ก็ทำได้ ซึ่งแนะนำให้ใช้ .. cond(เงื่อนไข,สิ่งที่จะทำเมื่อ True, สิ่งที่จะทำเมื่อ False) ..

สร้างตัวแปรก่อน
. gen hsCRP_level = cond(hsCRP>3,2,cond(hsCRP<1,0,1))
ตรวจสอบผลลัพธิ์
. tab hsCRP_level

hsCRP_level |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |      2,257       50.16       50.16
          1 |      1,327       29.49       79.64
          2 |        916       20.36      100.00
------------+-----------------------------------
      Total |      4,500      100.00
ควรสร้าง Label เพื่อจะได้ไม่ต้องจำค่า
. label define HSCRP_Level 0 "Low" 1 "Moderate" 2 "High"
. label values hsCRP_level HSCRP_Level

แล้วลองใช้ tab ซึ่งก็จะได้ผลเหมือนกัน

. tab PA_level if hsCRP <1

   PA_level |      Freq.     Percent        Cum.
------------+-----------------------------------
        Low |      1,020       45.19       45.19
   Moderate |        672       29.77       74.97
       High |        565       25.03      100.00
------------+-----------------------------------
      Total |      2,257      100.00

6. จงหา Odds ration ระหว่างผู้ที่มีความเสี่ยงโรคหัวใจกับการมีกิจกรรมทางกายเพียงพอ

ในข้อนี้ ผู้ที่มีความเสี่ยงโรคหัวใจ คือผู้ที่มี hsCRP > 1.0 หรือ Moderate ขึ้นไป เราจำเป็นต้องสร้างตัวแปรใหม่ขึ้นมา โดยให้ cad_risk 1 ผู้มีความเสี่ยง 0 ผู้ไม่มีความเสี่ยง และการมีกิจกรรมทางกายเพียงพอ คือผู้ที่มี PA_level > low ดังนั้นเราจะกำหนดตามโจทย์ โดยให้ good_PA 1 ผู้มีกิจกรรมทางกายเพียงพอ 0 ผู้มีกิจกรรมทางกายไม่เพียงพอ (ที่กำหนดแบบนี้เพราะโจทย์ถามของผู้ที่มีกิจกรรมเพียงพอ ไม่ใช่ของผู้ที่ไม่เพียงพอ)

. gen cad_risk = cond(hsCRP>1,1,0)
. gen good_PA = cond(PA_level>1,1,0)

จากนั้นโจทย์ถามหา Odds ratio เราก็จะใช้ cc case_variable exposure_varible หรือ logistic dependent_variable independent_variable ก้ได้

. cc cad_risk good_PA
                                                         Proportion
                 |   Exposed   Unexposed  |      Total     Exposed
-----------------+------------------------+------------------------
           Cases |      1085        1151  |       2236       0.4852
        Controls |      1237        1027  |       2264       0.5464
-----------------+------------------------+------------------------
           Total |      2322        2178  |       4500       0.5160
                 |                        |
                 |      Point estimate    |    [95% Conf. Interval]
                 |------------------------+------------------------
      Odds ratio |         .7826276       |    .6948771    .8814537 (exact)
 Prev. frac. ex. |         .2173724       |    .1185463    .3051229 (exact)
 Prev. frac. pop |         .1187675       |
                 +-------------------------------------------------
                               chi2(1) =    16.84  Pr>chi2 = 0.0000

. logistic cad_risk good_PA

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(1)        =      16.85
                                                Prob > chi2       =     0.0000
Log likelihood =  -3110.652                     Pseudo R2         =     0.0027

------------------------------------------------------------------------------
    cad_risk | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     good_PA |   .7826276   .0467785    -4.10   0.000       .69611    .8798982
       _cons |    1.12074   .0481073     2.66   0.008     1.030309    1.219108
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

ก็จะได้ว่า Odds ratio ของความเสี่ยงโรคหัวใจในผู้ที่มีกิจกรรมทางกายเพียงพอคือ .7826276 เมื่อเทียบกับ ผู้ที่มีกิจกรรมทางกายไม่เพียงพอ

7. ศึกษาความสัมพันธ์ระหว่าง ความเสี่ยงโรคหัวใจด้วยระดับ hsCRP กับการมีกิจกรรมทางกายเพียงพอ ต้องควบคุมตัวแปร เพศ และ อายุ หรือไม่เพราะเหตุใด และเขียนสมการ regression

จะต้องควบคุมตัวแปรหรือไม่

ตัวแปร เพศ และ อายุ นั้นส่งผลต่อ Outcome ก็คือความเสี่ยงต่อโรคหัวใจ แต่ว่าแล้วจะรู้ได้อย่างไรว่า 2 ตัวแปรนี้มัน Associate กับตัวแปรต้น คือ กิจกรรมทางกาย PA_level โจทย์ข้อนี้เดาว่าน่าจะหมายถึง outcome แบบ Nominal มากกว่า Ratio (จากคำว่า ความเสี่ยงโรคหัวใจด้วยระดับ hsCRP)

วิเคราหะ์แบบ Outcome คือ เสี่ยงหรือไม่เสี่ยงโรคหัวใจ (cad_risk) : Logistic regression

เพศ

สำหรับกรณีตัวแปรเพศนั้น ก็ต้องดูการกระจายระหว่าง ตัวแปร กิจกรรมทางกาย PA_level

. tab sex good_PA, col

+-------------------+
| Key               |
|-------------------|
|     frequency     |
| column percentage |
+-------------------+

           |        good_PA
       sex |         0          1 |     Total
-----------+----------------------+----------
      Male |       900      1,138 |     2,038 
           |     41.32      49.01 |     45.29 
-----------+----------------------+----------
    Female |     1,278      1,184 |     2,462 
           |     58.68      50.99 |     54.71 
-----------+----------------------+----------
     Total |     2,178      2,322 |     4,500 
           |    100.00     100.00 |    100.00 

มันก็ต่างกันประมาณ 8% ก็แล้วแต่คุณว่าจะคิดว่ามันต่างหรือไม่ต่างหรือไม่ อาจจะลองเปรียบเทียบ 10% ดูก็ได้

. logistic cad_risk good_PA

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(1)        =      16.85
                                                Prob > chi2       =     0.0000
Log likelihood =  -3110.652                     Pseudo R2         =     0.0027

------------------------------------------------------------------------------
    cad_risk | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     good_PA |   .7826276   .0467785    -4.10   0.000       .69611    .8798982
       _cons |    1.12074   .0481073     2.66   0.008     1.030309    1.219108
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

. logistic cad_risk good_PA sex

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(2)        =      17.18
                                                Prob > chi2       =     0.0002
Log likelihood = -3110.4841                     Pseudo R2         =     0.0028

------------------------------------------------------------------------------
    cad_risk | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     good_PA |   .7805176   .0467976    -4.13   0.000     .6939802    .8778459
         sex |   .9657216   .0581315    -0.58   0.562     .8582504    1.086651
       _cons |   1.184528   .1240739     1.62   0.106     .9646854     1.45447
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

ก็จะเห็นว่า Odds ratio เปลี่ยนไป .7826276.7805176.7826276=0.269%\dfrac{.7826276-.7805176}{.7826276} = 0.269\% ซึ่งแปลว่าถึงควบคุมผลก็เปลี่ยนไปเล็กน้อย ก็อาจจะคุมหรือไม่คุมก็ได้

อายุ

สำหรับกรณีตัวแปรอายุนั้น ก็ต้องดูการกระจายระหว่าง ตัวแปร กิจกรรมทางกาย PA_level

. tabstat age,stat(mean sd) by(good_PA)

Summary for variables: age
     by categories of: good_PA 

 good_PA |      mean        sd
---------+--------------------
       0 |  45.89578  7.017314
       1 |  45.36736  6.688671
---------+--------------------
   Total |  45.62311   6.85403
------------------------------

ดูไม่ค่อยต่างกันเท่าไหร่ อาจจะลองเปรียบเทียบ 10% เล่นๆ ดูก็ได้ แต่คิดว่าไม่จำเป็น

. logistic cad_risk good_PA

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(1)        =      16.85
                                                Prob > chi2       =     0.0000
Log likelihood =  -3110.652                     Pseudo R2         =     0.0027

------------------------------------------------------------------------------
    cad_risk | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     good_PA |   .7826276   .0467785    -4.10   0.000       .69611    .8798982
       _cons |    1.12074   .0481073     2.66   0.008     1.030309    1.219108
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

. logistic cad_risk good_PA age

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(2)        =      90.63
                                                Prob > chi2       =     0.0000
Log likelihood = -3073.7625                     Pseudo R2         =     0.0145

------------------------------------------------------------------------------
    cad_risk | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     good_PA |   .7953857   .0479496    -3.80   0.000      .706746    .8951424
         age |    1.03848    .004602     8.52   0.000     1.029499    1.047539
       _cons |   .1985395   .0411825    -7.79   0.000     .1322162    .2981324
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

ก็จะเห็นว่า Odds ratio เปลี่ยนไป .7953857.7826276.7826276=1.63%\dfrac{.7953857-.7826276}{.7826276} = 1.63\% ซึ่งแปลว่าถึงควบคุมผลก็เปลี่ยนไปเล็กน้อย ก็อาจจะคุมหรือไม่คุมก็ได้

สรุป

แต่ในข้อนี้ผมเลือกที่จะคุม มันเพราะมันก็เป็น Confounder ทั้งคู่ และเรามี Event เพียงพอที่จะคุมมันได้ ที่นี้เราจะเขียนสมการเราต้องการ Coefficient เราจะใช้ logit y x หรือ logistic y x, coef ก็ได้

. logistic cad_risk good_PA sex , coef
. logit cad_risk good_PA sex

Iteration 0:   log likelihood = -3119.0752  
Iteration 1:   log likelihood = -3110.4841  
Iteration 2:   log likelihood = -3110.4841  

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(2)        =      17.18
                                                Prob > chi2       =     0.0002
Log likelihood = -3110.4841                     Pseudo R2         =     0.0028

------------------------------------------------------------------------------
    cad_risk |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     good_PA |   -.247798   .0599572    -4.13   0.000    -.3653119   -.1302842
         sex |  -.0348797   .0601949    -0.58   0.562    -.1528594    .0831001
       _cons |    .169344   .1047454     1.62   0.106    -.0359533    .3746413
------------------------------------------------------------------------------

ดังนั้นเราจะสามารถเขียนสมการได้ว่า

Log(Oddsความเสี่ยงต่อโรคหัวใจ)=.169344.0348797(เพศ).247798(การมีกิจกรรมทางกายเพียงพอ)Log (Odds ความเสี่ยงต่อโรคหัวใจ) = .169344 -.0348797(เพศ) -.247798(การมีกิจกรรมทางกายเพียงพอ)
Logit(ความเสี่ยงต่อโรคหัวใจ)=.169344.0348797(เพศ).247798(การมีกิจกรรมทางกายเพียงพอ)Logit (ความเสี่ยงต่อโรคหัวใจ) = .169344 -.0348797(เพศ) -.247798(การมีกิจกรรมทางกายเพียงพอ)

โดย เพศ : 1 = ชาย , 2 = หญิง

8. ระดับ hsCRP สัมพันธ์กับระดับการมีกิจกรรมทางกายหรือไม่ เมื่อไม่ควบคุม เพศ และ อายุ

ตัวแปรต้นเป็น Oridal ตัวแปรตามเป็น Ratio ก็คงต้องใช้ Linear Regression regress

. regress hsCRP i.PA_level

      Source |       SS           df       MS      Number of obs   =     4,500
-------------+----------------------------------   F(2, 4497)      =     47.23
       Model |  421.038286         2  210.519143   Prob > F        =    0.0000
    Residual |  20044.0209     4,497  4.45719833   R-squared       =    0.0206
-------------+----------------------------------   Adj R-squared   =    0.0201
       Total |  20465.0592     4,499  4.54880177   Root MSE        =    2.1112

------------------------------------------------------------------------------
       hsCRP |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    PA_level |
   Moderate  |   .1889763   .0724464     2.61   0.009     .0469458    .3310067
       High  |  -.6540248   .0826992    -7.91   0.000    -.8161558   -.4918938
             |
       _cons |    1.99966   .0452379    44.20   0.000     1.910972    2.088349
------------------------------------------------------------------------------
ถ้าอยากจะดู Standardized Coefficient ก็ทำได้
. regress hsCRP i.PA_level,beta

      Source |       SS           df       MS      Number of obs   =     4,500
-------------+----------------------------------   F(2, 4497)      =     47.23
       Model |  421.038286         2  210.519143   Prob > F        =    0.0000
    Residual |  20044.0209     4,497  4.45719833   R-squared       =    0.0206
-------------+----------------------------------   Adj R-squared   =    0.0201
       Total |  20465.0592     4,499  4.54880177   Root MSE        =    2.1112

------------------------------------------------------------------------------
       hsCRP |      Coef.   Std. Err.      t    P>|t|                     Beta
-------------+----------------------------------------------------------------
    PA_level |
   Moderate  |   .1889763   .0724464     2.61   0.009                 .0409595
       High  |  -.6540248   .0826992    -7.91   0.000                -.1241815
             |
       _cons |    1.99966   .0452379    44.20   0.000                        .
------------------------------------------------------------------------------

เราก็จะเห็นว่า ระดับ hsCRP มีความสัมพันธ์กับระดับการมีกิจกรรมทางกาย

9. ความสัมพันธ์ระหว่าง ความเสี่ยงโรคหัวใจด้วยระดับ hsCRP กับระดับการมีกิจกรรมทางกาย โดยควบคุม ตัวแปร เพศ และ อายุ

ข้อนี้คล้ายๆ กับข้อ 7 แต่ว่า x คือ ระดับการมีกิจกรรมทางกาย

. logistic cad_risk i.PA_level sex age

Logistic regression                             Number of obs     =      4,500
                                                LR chi2(4)        =     127.94
                                                Prob > chi2       =     0.0000
Log likelihood = -3055.1068                     Pseudo R2         =     0.0205

------------------------------------------------------------------------------
    cad_risk | Odds Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
    PA_level |
   Moderate  |    .965501   .0668805    -0.51   0.612     .8429267      1.1059
       High  |   .5754727   .0468365    -6.79   0.000     .4906222    .6749976
             |
         sex |   .8693846   .0538484    -2.26   0.024     .7699984     .981599
         age |   1.038911   .0046511     8.53   0.000     1.029835    1.048067
       _cons |    .243236   .0540969    -6.36   0.000     .1572955    .3761313
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

เฉพาะ ผู้ที่มีกิจกรรมทางกายในระดับสูง (High) จะมีแต้มต่อของความเสี่ยงโรคหัวใจ 0.575 (0.491 - 0.675) เมื่อเทียบกับผู้ที่มีกิจกรรมทางกายระดับต่ำ (Low) ส่วนผู้ที่มีกิจกรรมทางกายในระดับปานกลาง มีความเสี่ยงโรคหัวใจไม่แตกต่างกับผู้ที่มีกิจกรรมทางกายระดับต่ำ