Model Regresi Poisson untuk Kasus Tabel \(2 \times 2\)

Menurut Selvin (2011), fungsi penghubung (link function) pada model Poisson adalah logaritma dari peubah dependen. Peubah dependen pada kasus ini berupa peluang kejadian. Jika faktor resiko dan kejadian yang diamati berupa peubah berkategori biner, maka model regresi dapat dituliskan sebagai berikut: \[ log(p_F)=a+bF \] dimana \(F=0\) menunjukkan bahwa tidak terdapat faktor resiko, dan \(F=1\) menunjukkan terdapat faktor resiko. Sehingga, dapat dituliskan log-peluang sebagai berikut: \[ F=0 : log(p_0)=a \]

dan

\[ F=1 : log(p_1)=a+b \]

Kedua persamaan log-peluang tersebut dapat digunakan untuk menghitung resiko relatif (\(rr\)) sebagai berikut: \[ log(p_1)-log(p_0)=(a+b)-a=b=log \bigg[ \dfrac{p_1}{p_0} \bigg] =log(rr) \]

Sehingga \[ relative \ risk=rr=e^b \]

Sebagai ilustrasi, perhatikan tabel berukuran \(2 \times 2\) berikut ini: \[\begin{matrix} \begin{array}{lrrr} \hline & & \textbf{ } & \textbf{ CHD} & \textbf{ No CHD} & \textbf{Total} \\ \hline && \text{Type-A} & 178 & 1411 & 1589 \\ && \text{Type-B} & 79 & 1486 & 1565 \\ && \textbf{Total} & 257 & 2897 & 3154 \\ \hline \end{array} \end{matrix} \]

Berdasarkan data pada tabel di atas, dapat dihitung resiko relatif sebagai rasio antara peluang kejadian CHD pada individu tipe A terhadap peluang kejadian CHD pada individu tipe B.

\[ relative \ risk=rr=\frac{P(CHD| individu \ type-A)}{P(CHD| individu \ type-B)}=\frac{p_A}{p_B} \]

Resiko Relatif

Perhitungan menggunakan program R dapat dilakukan dengan syntax berikut.

CHD.data<-matrix(c(178, 1411, 79, 1486), nrow=2, ncol=2, byrow=T)
dimnames(CHD.data)<-list(Factor=c("Type-A", "Type-B"),
                         Event=c("CHD", "No CHD"))

print(CHD.data)
##         Event
## Factor   CHD No CHD
##   Type-A 178   1411
##   Type-B  79   1486
library(DescTools)
RelRisk(CHD.data, conf.level = 0.95)
## rel. risk    lwr.ci    upr.ci 
##  2.219133  1.720682  2.864865

Pendugaan Parameter Regresi Poisson

Total<-rowSums(CHD.data)
CHD.df<-as.data.frame.table(CHD.data)
print(CHD.df)
##   Factor  Event Freq
## 1 Type-A    CHD  178
## 2 Type-B    CHD   79
## 3 Type-A No CHD 1411
## 4 Type-B No CHD 1486
CHD.df<-CHD.df[1:2,]
CHD.df$Total<-Total
CHD.df
##   Factor Event Freq Total
## 1 Type-A   CHD  178  1589
## 2 Type-B   CHD   79  1565
coba<-CHD.df
CHD.df<-within(CHD.df, Factor <- relevel(Factor, ref = 2))


regpois<-glm(Freq~Factor,offset=log(Total),family=poisson(link="log"), data=CHD.df)
summary(regpois)
## 
## Call:
## glm(formula = Freq ~ Factor, family = poisson(link = "log"), 
##     data = CHD.df, offset = log(Total))
## 
## Deviance Residuals: 
## [1]  0  0
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.9862     0.1125 -26.542  < 2e-16 ***
## FactorType-A   0.7971     0.1352   5.896 3.72e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  3.7648e+01  on 1  degrees of freedom
## Residual deviance: -1.5321e-14  on 0  degrees of freedom
## AIC: 17.23
## 
## Number of Fisher Scoring iterations: 2

Kita dapat hitung bahwa jika Type-B dianggap sebagai \(F=0\) , maka intersep pada model dapat pula dihitung sebagai \[ log(p_0)=log(p_A)=log\Bigg (\frac{79}{1565} \Bigg )=-2.986193 \]

Sementara itu dapat pula dihitung bahwa: \[ log(p_1)=a+b \\ log(p_B)=a+b \\ log\Bigg (\frac{178}{1589} \Bigg )=a+b \\ -2.189077 = -2.986193 +b \\ b=0.7971164 \]

Penduga Risiko Relatif

exp(regpois$coefficients[2])
## FactorType-A 
##     2.219133

Prediksi Peluang Kejadian

#prediksi jika diketahui objek berasal dari kelompok Type-B
#P_A=exp(a)


exp(coef(regpois)[1])
## (Intercept) 
##  0.05047923
#prediksi jika diketahui objek berasal dari kelompok Type-A
exp(coef(regpois)[1]+coef(regpois)[2])
## (Intercept) 
##   0.1120201

Kasus Data Kategorik

Mengimpor Data

CHD.smoking<-read.csv("https://raw.githubusercontent.com/raoy/data/master/chd%20smoking.csv")

CHD.smoking$Smoking<-as.factor(CHD.smoking$Smoking)
CHD.smoking$Factor<-as.factor(CHD.smoking$Factor)

CHD.smoking <- within(CHD.smoking, Smoking <- relevel(Smoking, ref = 4))
CHD.smoking <- within(CHD.smoking, Factor <- relevel(Factor, ref = 2))


levels(CHD.smoking$Factor)
## [1] "Type-B" "Type-A"
levels(CHD.smoking$Smoking)
## [1] "Nonsmoker" ">30"       "1-20"      "21-30"
CHD.smoking
##   Factor   Smoking CHD No.CHD Total
## 1 Type-A Nonsmoker  70    714   784
## 2 Type-A      1-20  45    353   398
## 3 Type-A     21-30  34    185   219
## 4 Type-A       >30  29    159   188
## 5 Type-B Nonsmoker  28    840   868
## 6 Type-B      1-20  25    382   407
## 7 Type-B     21-30  16    170   186
## 8 Type-B       >30  10     94   104

Pemodelan Regresi Poisson

regpois2<-glm(CHD~Factor+Smoking,offset=log(Total),family=poisson(link="log"), data=CHD.smoking)


summary(regpois2)
## 
## Call:
## glm(formula = CHD ~ Factor + Smoking, family = poisson(link = "log"), 
##     data = CHD.smoking, offset = log(Total))
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  0.7025  -0.3165  -0.2794  -0.3458  -1.0141   0.4445   0.4290   0.6446  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.2483     0.1347 -24.111  < 2e-16 ***
## FactorType-A   0.7472     0.1359   5.498 3.83e-08 ***
## Smoking>30     0.6955     0.1902   3.656 0.000256 ***
## Smoking1-20    0.3681     0.1565   2.352 0.018676 *  
## Smoking21-30   0.6859     0.1740   3.943 8.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 62.1571  on 7  degrees of freedom
## Residual deviance:  2.6169  on 3  degrees of freedom
## AIC: 53.951
## 
## Number of Fisher Scoring iterations: 4

Prediksi

pred<-predict(regpois2, type="response")
df.pred<-data.frame(CHD=CHD.smoking$CHD, "CHD-pred"=pred)

df.pred
##   CHD CHD.pred
## 1  70 64.28557
## 2  45 47.15689
## 3  34 35.65530
## 4  29 30.90225
## 5  28 33.71443
## 6  25 22.84311
## 7  16 14.34470
## 8  10  8.09775

Goodness of Fit Test

with(regpois2, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance df         p
## [1,]     2.616916  3 0.4545317

Berdasarkan output di atas, kita gagal menolak \(H_0\) pada taraf nyata 5%, karena p-value=0.454 > 0.05. Artinya, model dapat dianggap fit dengan data.

Latihan

Gunakan data warpbreaks pada package datasets.

# install.packages("datasets")
library(datasets) # include library datasets after installation
data <- warpbreaks
names(data)
## [1] "breaks"  "wool"    "tension"

Lakukan pemodelan regresi Poisson dengan menggunakan peubah breaks sebagai respon, serta wool dan tension sebagai peubah prediktor.

Selanjutnya, coba interpretasikan hasil yang Anda peroleh.

Referensi

Selvin, S. (2011). Statistical tools for epidemiologic research. OUP USA.