Data yang akan digunakan berasal dari Tabel 11.2 dalam Agresti (2002) yang berkaitan dengan penelitian longitudinal. Penelitian bertujuan untuk membandingkan dua jenis obat (“new” vs. “standard”) untuk pengobatan depresi. Variabel respon adalah kondisi depresi (depression) (“normal” vs. “abnormal”). Kategori Y=1 adalah “normal”, sementara Y=0 adalah “abnormal” yang merupakan kategori referensi.

Variabel prediktor yang digunakan adalah:

  1. variabel diagnosis awal, yakni “mild” dan “severe”. Kategori “mild” merupakan kategori referensi.

  2. jenis obat untuk pengobatan depresi, yakni “new” vs. “standard”. Kategori “standard” merupakan kategori referensi.

  3. waktu pemberian obat: 0, 1 dan 2

Berikut adalah tampilan tiga baris pertama.

Dari data, kita melihat bahwa subjek 1 (id = 1) memiliki diagnosis “mild” dan menerima obat “standard” pada waktu 0, 1, dan 2.

URL <- "http://static.lib.virginia.edu/statlab/materials/data/depression.csv"
dat <- read.csv(URL, stringsAsFactors = TRUE)
dat$id <- factor(dat$id)
dat$drug <- relevel(dat$drug, ref = "standard")
head(dat, n = 3)
##   diagnose     drug id time depression
## 1     mild standard  1    0          1
## 2     mild standard  1    1          1
## 3     mild standard  1    2          1

Secara keseluruhan, pasien yang menjadi unit analisis sebanyak 340

length(unique(dat$id))
## [1] 340

Selanjutnya akan dilihat proporsi respons “Normal” (depression == 1) seiring waktu dalam setiap kategori diagnosis dan jenis obat.

Berikut adalah tabel rata-rata berdasarkan diagnosis, obat, dan waktu.

Note: Cara menghitung Proporsi “Normal” yakni misal jika di suatu kelompok ada 10 pasien, dan 7 di antaranya memiliki respons 1 (Normal), maka proporsi “Normal” adalah 7/10 = 0.7.

Data dikelompokkan berdasarkan:

Diagnosis (“mild” atau “severe”)

Obat (“standard” atau “new”)

Waktu (“0”, “1”, “2”)

Output ini mereproduksi Tabel 11.3 dalam Agresti (2002).

Dapat dilihat bahwa proporsi respons “Normal” meningkat seiring waktu untuk keempat kombinasi diagnosis dan pengobatan, dan terjadi peningkatan yang lebih drastis pada obat “new” dibandingkan obat “standard”.

library(magrittr)
with(dat, tapply(depression, list(diagnose, drug, time), mean)) %>% 
  ftable() %>% 
  round(2)
##                     0    1    2
##                                
## mild   standard  0.51 0.59 0.68
##        new       0.53 0.79 0.97
## severe standard  0.21 0.28 0.46
##        new       0.18 0.50 0.83

Pemodelan dengan Generalized Estimating Equations (GEE)

Pemodelan GEE digunakan untuk mengestimasi model marginal untuk mengetahui pengaruh diagnosis, obat, dan waktu terhadap variabel respons depresi.

Dalam penelitian ini akan dimasukkan interaksi antara obat dan waktu, untuk menunjukkan bahwa variabel respons bervariasi tergantung pada jenis obat.

Note: Argumen pada GEE dalam `gee()’ yakni:

data → Menentukan data frame yang digunakan.

family → Misalnya binomial untuk data biner.

id → Menentukan faktor pengelompokan. Dalam kasus ini, faktor pengelompokan adalah kolom id dalam data frame, yang menunjukkan setiap individu dalam penelitian.

corstr → Menentukan struktur korelasi antar pengamatan dalam satu individu.

# install.packages("gee")
library(gee) 
dep_gee <- gee(depression ~ diagnose + drug*time,
               data = dat, 
               id = id, 
               family = binomial,
               corstr = "independence", scale.fix = TRUE)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
##    (Intercept) diagnosesevere        drugnew           time   drugnew:time 
##    -0.02798843    -1.31391092    -0.05960381     0.48241209     1.01744498
summary(dep_gee)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Independent 
## 
## Call:
## gee(formula = depression ~ diagnose + drug * time, id = id, data = dat, 
##     family = binomial, corstr = "independence", scale.fix = TRUE)
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.94844242 -0.40683252  0.05155758  0.38830952  0.80242231 
## 
## 
## Coefficients:
##                   Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept)    -0.02798843  0.1639083 -0.1707566   0.1741865 -0.1606808
## diagnosesevere -1.31391092  0.1464151 -8.9738733   0.1459845 -9.0003423
## drugnew        -0.05960381  0.2222080 -0.2682343   0.2285385 -0.2608042
## time            0.48241209  0.1147626  4.2035644   0.1199350  4.0222784
## drugnew:time    1.01744498  0.1887954  5.3891398   0.1876938  5.4207709
## 
## Estimated Scale Parameter:  1
## Number of Iterations:  1
## 
## Working Correlation
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Koefisien yang ditampilkan berada dalam skala logit, sehingga interpretasi koefisien dilakukan dengan cara yang sama seperti model regresi logistik pada umumnya.

🔹 Interpretasi Koefisien Waktu time

Koefisien waktu time adalah 0.48. Interpretasi: Kecenderungan pasien yang mengkonsumsi obat standar meningkat sebesar (exp(0.48)=1.61) kali untuk mendapatkan respon “Normal” pada setiap kenaikan satu unit waktu.

🔹 Interpretasi Interaksi (drugnew:time)

Koefisien interaksi antara obat baru dan waktu (drugnew:time) adalah 1.02. Untuk menghitung efek waktu pada obat baru, kita menambahkan koefisien waktu (0.48) dengan koefisien interaksi (1.02), yakni 0.48 + 1.02 = 1.5

Jika dieksponensialkan (exp(1.5)=4.5). Artinya, kecenderungan pasien yang mengkonsumsi obat baru untuk mendapatkan respons “Normal” 4.5 kali lebih besar untuk setiap kenaikan satu unit waktu.

Pada output tersebut working correlation structure yang digunakan adalah independence, artinya diasumsikan tidak terdapat korelasi antar pengamatan di dalam subjek. Matriks berukuran 3 x 3 karena terdapat 3 observasi per subjek.

Selanjutnya akan dimodelkan jika working correlation structure yang digunakan adalah exchangeable

dep_gee2 <- gee(depression ~ diagnose + drug*time,
               data = dat, 
               id = id, 
               family = binomial,
               corstr = "exchangeable", scale.fix = TRUE)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
##    (Intercept) diagnosesevere        drugnew           time   drugnew:time 
##    -0.02798843    -1.31391092    -0.05960381     0.48241209     1.01744498
summary(dep_gee2)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Exchangeable 
## 
## Call:
## gee(formula = depression ~ diagnose + drug * time, id = id, data = dat, 
##     family = binomial, corstr = "exchangeable", scale.fix = TRUE)
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.94843397 -0.40683122  0.05156603  0.38832332  0.80238627 
## 
## 
## Coefficients:
##                   Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept)    -0.02809866  0.1637503 -0.1715945   0.1741791 -0.1613205
## diagnosesevere -1.31391033  0.1459325 -9.0035505   0.1459630 -9.0016667
## drugnew        -0.05926689  0.2221626 -0.2667725   0.2285569 -0.2593091
## time            0.48246420  0.1149581  4.1968686   0.1199383  4.0226037
## drugnew:time    1.01719312  0.1890913  5.3793750   0.1877014  5.4192084
## 
## Estimated Scale Parameter:  1
## Number of Iterations:  2
## 
## Working Correlation
##              [,1]         [,2]         [,3]
## [1,]  1.000000000 -0.003432732 -0.003432732
## [2,] -0.003432732  1.000000000 -0.003432732
## [3,] -0.003432732 -0.003432732  1.000000000

Matriks “Working Correlation” menunjukkan estimasi sebesar -0.003433 sebagai korelasi umum antara pasangan observasi dalam satu subjek. Nilai ini relatif sangat kecil, sehingga sebetulnya data ini bisa saja diperlakukan sebagai 1020 observasi independen, bukan 340 subjek dengan masing-masing 3 observasi. Selain itu, estimasi koefisien hampir identik dengan model sebelumnya yang menggunakan struktur korelasi independen.

Struktur korelasi lainnya yang dapat digunakan adalah struktur autoregresif. Struktur ini memungkinkan korelasi antara pengukuran yang diambil dalam waktu yang lebih dekat menjadi lebih tinggi dibandingkan dengan pengukuran yang diambil dalam waktu yang lebih jauh.

Struktur autoregresif yang paling umum adalah AR-1. Dalam penelitian ini berarti bahwa pengukuran pada waktu 1 dan 2 memiliki korelasi tertentu, sedangkan pengukuran pada waktu 1 dan 3 memiliki korelasi yang lebih kecil.

dep_gee3 <- gee(depression ~ diagnose + drug*time,
               data = dat, 
               id = id, 
               family = binomial,
               corstr = "AR-M", Mv = 1, scale.fix = TRUE)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
##    (Intercept) diagnosesevere        drugnew           time   drugnew:time 
##    -0.02798843    -1.31391092    -0.05960381     0.48241209     1.01744498
summary(dep_gee3)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     AR-M , M = 1 
## 
## Call:
## gee(formula = depression ~ diagnose + drug * time, id = id, data = dat, 
##     family = binomial, corstr = "AR-M", Mv = 1, scale.fix = TRUE)
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.94844464 -0.40691023  0.05155536  0.38824284  0.80236892 
## 
## 
## Coefficients:
##                   Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept)    -0.02770314  0.1643892 -0.1685216   0.1741163 -0.1591071
## diagnosesevere -1.31386512  0.1472418 -8.9231818   0.1459894 -8.9997314
## drugnew        -0.05959816  0.2226449 -0.2676826   0.2284592 -0.2608700
## time            0.48240753  0.1147578  4.2036998   0.1199244  4.0225960
## drugnew:time    1.01732678  0.1887941  5.3885527   0.1876727  5.4207510
## 
## Estimated Scale Parameter:  1
## Number of Iterations:  2
## 
## Working Correlation
##              [,1]        [,2]         [,3]
## [1,] 1.000000e+00 0.008477443 7.186704e-05
## [2,] 8.477443e-03 1.000000000 8.477443e-03
## [3,] 7.186704e-05 0.008477443 1.000000e+00

Dari hasil tersebut “Working Correlation” menunjukkan nilai yang sangat kecil. Parameter korelasi diestimasi sebesar 0.008477 (korelasi antara pengukuran 1-2 dan 2-3).

Salah satu struktur korelasi lain adalah matriks unstructured. Struktur ini memungkinkan semua korelasi bervariasi secara bebas.

dep_gee4 <- gee(depression ~ diagnose + drug*time,
               data = dat, 
               id = id, 
               family = binomial,
               corstr = "unstructured", scale.fix = TRUE)
## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate
##    (Intercept) diagnosesevere        drugnew           time   drugnew:time 
##    -0.02798843    -1.31391092    -0.05960381     0.48241209     1.01744498
summary(dep_gee4)
## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Unstructured 
## 
## Call:
## gee(formula = depression ~ diagnose + drug * time, id = id, data = dat, 
##     family = binomial, corstr = "unstructured", scale.fix = TRUE)
## 
## Summary of Residuals:
##         Min          1Q      Median          3Q         Max 
## -0.94773674 -0.40645713  0.05226326  0.38927858  0.79975454 
## 
## 
## Coefficients:
##                   Estimate Naive S.E.    Naive z Robust S.E.   Robust z
## (Intercept)    -0.02552611  0.1679741 -0.1519645   0.1726392 -0.1478581
## diagnosesevere -1.30484850  0.1461691 -8.9269772   0.1450136 -8.9981088
## drugnew        -0.05438636  0.2282121 -0.2383149   0.2271321 -0.2394481
## time            0.47587182  0.1160832  4.0994035   0.1190418  3.9975178
## drugnew:time    1.01297603  0.1887379  5.3671034   0.1865407  5.4303205
## 
## Estimated Scale Parameter:  1
## Number of Iterations:  3
## 
## Working Correlation
##             [,1]        [,2]        [,3]
## [1,]  1.00000000  0.07393977 -0.02741128
## [2,]  0.07393977  1.00000000 -0.05669559
## [3,] -0.02741128 -0.05669559  1.00000000

Bagaimana cara memilih struktur korelasi yang tepat? Menurut Agresti (2002), estimasi GEE tetap valid meskipun struktur korelasi yang digunakan tidak tepat. Agresti menyarankan untuk memulai dengan struktur exchangeable, lalu memeriksa bagaimana estimasi koefisien dan standar error berubah ketika menggunakan struktur korelasi lain. Jika perubahannya minimal, maka pilihlah struktur korelasi yang lebih sederhana.

=======================================================================================================

Sumber:

Agresti, A. (2002). Categorical data analysis (2nd ed.). Wiley

https://library.virginia.edu/data/articles/getting-started-with-generalized-estimating-equations