Model Log-Linear untuk Data Respon Poisson

Melakukan Set Direktori

setwd("D:/KULIAH IPB/SEMESTER 3/STK 654 ANALISIS DATA KATEGORIK/UAS/PRAKTIKUM")

Penyelesaian untuk soal dari Buku “Categorical Data Analysis for the Behavioral and Social Sciences” (Razia Azen dan Cindy M.Walker ) Problem 7.1

Tabel 7.20 menggambarkan data yang diambil dari data administrasi GSS tahun 1994 dengan tiga peubah yaitu: jenis kelamin, tingkat fundamentalism, dan setuju atau tidak tentang hukuman mati untuk pembunuh. Duga dan interpretasikan model log-linear yang paling parsimoni. Deskripsikan proses yang Anda lakukan untuk memperoleh model dan dugaan terbaik. Apakah ada hubungan? Uraikan bagaimana memperoleh nilai dugaan dari model yang menggambarkan asosiasi antarpeubah

Memasukkan Data Tabel Kontingensi ke R

x.sex <- factor(rep(c("M","F"),6))
y.fav <- factor(rep(c("Fav","Opp"), rep(6,2)))
z.fund <- rep((factor(rep(c("Fund", "Mod","Lib"),rep(2,3)))),2)
counts <- c(128,123,182,168,119, 111,32,73,56,105,49,70)
data.frame(z.fund,x.sex,y.fav,counts)
##    z.fund x.sex y.fav counts
## 1    Fund     M   Fav    128
## 2    Fund     F   Fav    123
## 3     Mod     M   Fav    182
## 4     Mod     F   Fav    168
## 5     Lib     M   Fav    119
## 6     Lib     F   Fav    111
## 7    Fund     M   Opp     32
## 8    Fund     F   Opp     73
## 9     Mod     M   Opp     56
## 10    Mod     F   Opp    105
## 11    Lib     M   Opp     49
## 12    Lib     F   Opp     70
  • data.frame dalam hal ini digunakan untuk mengecek apakah data yang diinput dalam R dan pada data asli sama. Dalam hal ini adalah sama, dimana jumlah data adalah 12

Menentukan Kategori Referensi

x.sex<-relevel(x.sex,ref="F")
y.fav<-relevel(y.fav,ref="Opp")
z.fund<-relevel(z.fund,ref="Lib")

Membentuk Model Saturated

model <- glm(counts~x.sex+y.fav+z.fund+x.sex*y.fav+x.sex*z.fund+ 
y.fav*z.fund+x.sex*y.fav*z.fund,family=poisson("link"=log))
summary(model)
## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund + x.sex * y.fav * z.fund, 
##     family = poisson(link = log))
## 
## Deviance Residuals: 
##  [1]  0  0  0  0  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 4.248495   0.119523  35.545  < 2e-16 ***
## x.sexM                     -0.356675   0.186263  -1.915  0.05551 .  
## y.favFav                    0.461035   0.152626   3.021  0.00252 ** 
## z.fundFund                  0.041964   0.167285   0.251  0.80193    
## z.fundMod                   0.405465   0.154303   2.628  0.00860 ** 
## x.sexM:y.favFav             0.426268   0.228268   1.867  0.06185 .  
## x.sexM:z.fundFund          -0.468049   0.282210  -1.659  0.09721 .  
## x.sexM:z.fundMod           -0.271934   0.249148  -1.091  0.27507    
## y.favFav:z.fundFund         0.060690   0.212423   0.286  0.77511    
## y.favFav:z.fundMod          0.008969   0.196903   0.046  0.96367    
## x.sexM:y.favFav:z.fundFund  0.438301   0.336151   1.304  0.19227    
## x.sexM:y.favFav:z.fundMod   0.282383   0.301553   0.936  0.34905    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.4536e+02  on 11  degrees of freedom
## Residual deviance: 2.6645e-14  on  0  degrees of freedom
## AIC: 100.14
## 
## Number of Fisher Scoring iterations: 3

Membentuk Homogenous Model

model2 <- glm(counts~x.sex+y.fav+z.fund+x.sex*y.fav+ 
x.sex*z.fund+y.fav*z.fund,family=poisson("link"=log))
summary(model2)
## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund + y.fav * z.fund, family = poisson(link = log))
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  0.2996  -0.3002   0.0856  -0.0887  -0.4085   0.4342  -0.5744   0.3994  
##       9       10       11       12  
## -0.1530   0.1129   0.6655  -0.5281  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.31096    0.10522  40.972  < 2e-16 ***
## x.sexM              -0.51575    0.13814  -3.733 0.000189 ***
## y.favFav             0.35707    0.12658   2.821 0.004788 ** 
## z.fundFund          -0.06762    0.14452  -0.468 0.639854    
## z.fundMod            0.33196    0.13142   2.526 0.011540 *  
## x.sexM:y.favFav      0.66406    0.12728   5.217 1.81e-07 ***
## x.sexM:z.fundFund   -0.16201    0.15300  -1.059 0.289649    
## x.sexM:z.fundMod    -0.08146    0.14079  -0.579 0.562887    
## y.favFav:z.fundFund  0.23873    0.16402   1.455 0.145551    
## y.favFav:z.fundMod   0.13081    0.14951   0.875 0.381614    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:   1.798  on  2  degrees of freedom
## AIC: 97.934
## 
## Number of Fisher Scoring iterations: 3

Pengujian Hipotesis Saturated vs Homogenous

Deviance.model_1 <- model2$deviance - model$deviance
Deviance.model_1
## [1] 1.797977
# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (2 - 0)
derajat.bebas
## [1] 2
chi.tabel <- qchisq((1-0.05), df=derajat.bebas) 
chi.tabel
## [1] 5.991465
Keputusan_1 <- ifelse(Deviance.model_1 <= chi.tabel, "Terima", "Tolak")
Keputusan_1
## [1] "Terima"

Pada taraf nyata 5%, belum cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa tidak ada interaksi tiga arahantara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati

Membentuk Conditional Association based on X

model3<-glm(counts~x.sex+y.fav+z.fund+x.sex*y.fav+ 
x.sex*z.fund,family=poisson("link"=log))
summary(model3)
## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     x.sex * z.fund, family = poisson(link = log))
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
##  0.60542   0.16142   0.11953  -0.06470  -0.74698  -0.08912  -1.11492  -0.20684  
##        9        10        11        12  
## -0.21283   0.08220   1.26595   0.11304  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        4.23495    0.08955  47.293  < 2e-16 ***
## x.sexM            -0.52960    0.13966  -3.792 0.000149 ***
## y.favFav           0.48302    0.08075   5.982 2.20e-09 ***
## z.fundFund         0.07962    0.10309   0.772 0.439916    
## z.fundMod          0.41097    0.09585   4.288 1.81e-05 ***
## x.sexM:y.favFav    0.65845    0.12708   5.181 2.20e-07 ***
## x.sexM:z.fundFund -0.12841    0.15109  -0.850 0.395405    
## x.sexM:z.fundMod  -0.06267    0.13908  -0.451 0.652274    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   3.9303  on  4  degrees of freedom
## AIC: 96.067
## 
## Number of Fisher Scoring iterations: 4

Pengujian Hipotesis Conditional on X vs Homogenous

Deviance.model_2 <- model3$deviance - model2$deviance
Deviance.model_2
## [1] 2.132302
# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas
## [1] 2
chi.tabel <- qchisq((1-0.05), df=derajat.bebas) 
chi.tabel
## [1] 5.991465
Keputusan_2 <- ifelse(Deviance.model_2 <= chi.tabel, "Terima", "Tolak")
Keputusan_2
## [1] "Terima"

Pada taraf nyata 5%, belum cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa tidak ada interaksi antara pendapat mengenai hukuman mati dan fundamentalisme

Membentuk Conditional Association based on Y

model4<-glm(counts~x.sex+y.fav+z.fund+x.sex*y.fav+y.fav*z.fund, 
family=poisson("link"=log))
summary(model4)
## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav + 
##     y.fav * z.fund, family = poisson(link = log))
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7         8  
## -0.13887   0.14286   0.09764  -0.10112   0.02418  -0.02499  -0.89983   0.64384  
##        9        10        11        12  
## -0.17120   0.12650   0.99745  -0.77148  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.33931    0.09919  43.748  < 2e-16 ***
## x.sexM              -0.59345    0.10645  -5.575 2.48e-08 ***
## y.favFav             0.37259    0.12438   2.996  0.00274 ** 
## z.fundFund          -0.12516    0.13389  -0.935  0.34989    
## z.fundMod            0.30228    0.12089   2.500  0.01240 *  
## x.sexM:y.favFav      0.65845    0.12708   5.181 2.20e-07 ***
## y.favFav:z.fundFund  0.21254    0.16205   1.312  0.18966    
## y.favFav:z.fundMod   0.11757    0.14771   0.796  0.42606    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   2.9203  on  4  degrees of freedom
## AIC: 95.057
## 
## Number of Fisher Scoring iterations: 4

Pengujian Hipotesis Conditional on Y vs Homogenous

Deviance.model_3 <- model4$deviance - model2$deviance
Deviance.model_3
## [1] 1.122315
# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (4 -2)
derajat.bebas
## [1] 2
chi.tabel <- qchisq((1-0.05), df=derajat.bebas) 
chi.tabel
## [1] 5.991465
Keputusan_3 <- ifelse(Deviance.model_3 <= chi.tabel, "Terima", "Tolak")
Keputusan_3
## [1] "Terima"

Pada taraf nyata 5%, belum cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa tidak ada interaksi antara jenis kelamin dan fundamentalisme

Membentuk Conditional Association based on Z

model5<-glm(counts~x.sex+y.fav+z.fund+x.sex*z.fund+y.fav*z.fund, 
family=poisson("link"=log))
summary(model5)
## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * z.fund + 
##     y.fav * z.fund, family = poisson(link = log))
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  1.3998  -1.3171   1.4595  -1.4130   0.7777  -0.7675  -2.3495   1.9188  
##       9       10       11       12  
## -2.2965   1.9781  -1.1226   1.0321  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          4.12255    0.10518  39.195  < 2e-16 ***
## x.sexM              -0.07453    0.10713  -0.696    0.487    
## y.favFav             0.65896    0.11292   5.836 5.36e-09 ***
## z.fundFund          -0.06540    0.15126  -0.432    0.665    
## z.fundMod            0.33196    0.13777   2.410    0.016 *  
## x.sexM:z.fundFund   -0.12841    0.15109  -0.850    0.395    
## x.sexM:z.fundMod    -0.06267    0.13908  -0.451    0.652    
## y.favFav:z.fundFund  0.21254    0.16205   1.312    0.190    
## y.favFav:z.fundMod   0.11757    0.14771   0.796    0.426    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.361  on 11  degrees of freedom
## Residual deviance:  29.729  on  3  degrees of freedom
## AIC: 123.87
## 
## Number of Fisher Scoring iterations: 4

Pengujian Hipotesis Conditional on Z vs Homogenous

Deviance.model_4 <- model5$deviance - model2$deviance
Deviance.model_4
## [1] 27.93095
# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (3 - 2)
derajat.bebas
## [1] 1
chi.tabel <- qchisq((1-0.05), df=derajat.bebas) 
chi.tabel
## [1] 3.841459
Keputusan_4 <- ifelse(Deviance.model_4 <= chi.tabel, "Terima", "Tolak")
Keputusan_4
## [1] "Tolak"

Pada taraf nyata 5%, cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa ada interaksi antara jenis kelamin dan pendapat mengenai hukuman mati

Penentuan Best Model

bestmodel<-glm(counts~x.sex+y.fav+z.fund+x.sex*y.fav, 
family=poisson("link"=log))
summary(bestmodel)
## 
## Call:
## glm(formula = counts ~ x.sex + y.fav + z.fund + x.sex * y.fav, 
##     family = poisson(link = log))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.32758  -0.24954  -0.01277   0.14946   1.48613  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      4.26518    0.07794  54.721  < 2e-16 ***
## x.sexM          -0.59345    0.10645  -5.575 2.48e-08 ***
## y.favFav         0.48302    0.08075   5.982 2.20e-09 ***
## z.fundFund       0.01986    0.07533   0.264    0.792    
## z.fundMod        0.38130    0.06944   5.491 4.00e-08 ***
## x.sexM:y.favFav  0.65845    0.12708   5.181 2.20e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 245.3612  on 11  degrees of freedom
## Residual deviance:   4.6532  on  6  degrees of freedom
## AIC: 92.79
## 
## Number of Fisher Scoring iterations: 4

Dari summary model diatas terlihat bahwa best model memiliki AIC yang lebih rendahdibandingkan saturated, homogeneous, dan conditional model

Interpretasi

data.frame(koef=bestmodel$coefficients,exp_koef=exp(bestmodel$coefficients))
##                        koef   exp_koef
## (Intercept)      4.26517861 71.1776316
## x.sexM          -0.59344782  0.5524194
## y.favFav         0.48302334  1.6209677
## z.fundFund       0.01985881  1.0200573
## z.fundMod        0.38129767  1.4641834
## x.sexM:y.favFav  0.65845265  1.9318008

Interpretasi Model Terbaik a. exp −0.59345 = 0.552. Interpretasi: dengan mengabaikan status fundamentalism dan opini terhadap hukuman mati bagi pembunuh, peluang seseorang berjenis kelamin laki-laki adalah 0.552 kali dibanding perempuan. Interpretasi akan menjadi lebih menarik apabila perbandingan dibalik menjadi seperti ini: dengan mengabaikan status fundamentalism dan opini terhadap hukuman mati bagi pembunuh, peluang seseorang berjenis kelamin perempuan adalah 1.8116 kali dibanding laki-laki.Angka 1.8116 diperoleh dari 1 dibagi 0.552. b. exp 0.48302 = 1.62. Interpretasi: dengan mengabaikan status fundamentalism dan jenis kelamin, peluang seseorang mendukung hukuman mati bagi pembunuh adalah 1,62 kali dibanding yang menolak. c.exp 0.01986 = 1.02. Interpretasi: dengan mengabaikan jenis kelamin dan opini tentang hukuman mati, peluang seseorag fundamentalist adalah 1,02 dibanding liberal. d. exp 0.38130 = 1.46. Interpretasi: dengan mengabaikan jenis kelamin dan opini tentang hukuman mati, peluang seseorang moderat adalah 1,02 dibanding liberal. e. exp 0.65845 = 1.93. Interpretasi: tanpa melihat fundamentalism, odds yang mendukung hukuman mati dibandingkan dengan yang menolak jika dia laki-laki adalah 1,93 dibanding odd sama jika dia perempuan

Nilai Dugaan

data.frame(Fund=z.fund,sex=x.sex,favor=y.fav,counts=counts, fitted=bestmodel$fitted.values)
##    Fund sex favor counts    fitted
## 1  Fund   M   Fav    128 125.59539
## 2  Fund   F   Fav    123 117.69079
## 3   Mod   M   Fav    182 180.27878
## 4   Mod   F   Fav    168 168.93257
## 5   Lib   M   Fav    119 123.12582
## 6   Lib   F   Fav    111 115.37664
## 7  Fund   M   Opp     32  40.10855
## 8  Fund   F   Opp     73  72.60526
## 9   Mod   M   Opp     56  57.57155
## 10  Mod   F   Opp    105 104.21711
## 11  Lib   M   Opp     49  39.31990
## 12  Lib   F   Opp     70  71.17763