Responsi 5 STA543 Analisis Data Kategorik
Model Log-Linear 3 Arah untuk Data Respon Poisson
setwd("D:\\Kuliah S2 IPB\\Bahan Kuliah\\Semester 2 SSD 2020\\STA543 ADK\\Responsi\\R\\UTS\\")Pendahuluan
Pada pertemuan sebelumnya kita telah mengetahui bahwa tujuan utama dari pembuatan model log linear adalah untuk menduga parameter model yang mendeskripsikan hubungan antar peubah kategorik.
Pada pertemuan kali ini, model yang kita miliki lebih kompleks karena merupakan model log linear untuk tabel kontingensi 3 arah, artinya terdapat tiga peubah pada model log linear tersebut. Dengan demikian, interaksi tertinggi yang terjadi adalah interaksi antara ketiga peubah yang berada di dalam model.
Jika pada urutan pertama hasil pengujian membuktikan ada interaksi, maka tidak perlu lanjut ke urutan kedua, dst.
Soal
Package yang digunakan
library("epitools")
library("DescTools")
library("lawstat")Input Data
##=====================##
# INPUT DATA
##=====================##
z.fund<-factor(rep(c("1fund","2mod","3lib"),each=4))
x.sex<-factor(rep(c("1M","2F"),each=2,times=3))
y.fav<-factor(rep(c("1fav","2opp"),times=6))
counts<-c(128,32,123,73,182,56,168,105,119,49,111,70)
data.frame(z.fund,x.sex,y.fav,counts) ## z.fund x.sex y.fav counts
## 1 1fund 1M 1fav 128
## 2 1fund 1M 2opp 32
## 3 1fund 2F 1fav 123
## 4 1fund 2F 2opp 73
## 5 2mod 1M 1fav 182
## 6 2mod 1M 2opp 56
## 7 2mod 2F 1fav 168
## 8 2mod 2F 2opp 105
## 9 3lib 1M 1fav 119
## 10 3lib 1M 2opp 49
## 11 3lib 2F 1fav 111
## 12 3lib 2F 2opp 70A. UJI MODEL INTERAKSI TIGA ARAH (SATURATED VS HOMOGENOUS)
Penentuan Kategori Referensi
##=============================##
# Penentuan kategori reference
##=============================##
x.sex<-relevel(x.sex,ref="2F")
y.fav<-relevel(y.fav,ref="2opp")
z.fund<-relevel(z.fund,ref="3lib")Model Saturated
#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.sex1M -0.356675 0.186263 -1.915 0.05551 .
## y.fav1fav 0.461035 0.152626 3.021 0.00252 **
## z.fund1fund 0.041964 0.167285 0.251 0.80193
## z.fund2mod 0.405465 0.154303 2.628 0.00860 **
## x.sex1M:y.fav1fav 0.426268 0.228268 1.867 0.06185 .
## x.sex1M:z.fund1fund -0.468049 0.282210 -1.659 0.09721 .
## x.sex1M:z.fund2mod -0.271934 0.249148 -1.091 0.27507
## y.fav1fav:z.fund1fund 0.060690 0.212423 0.286 0.77511
## y.fav1fav:z.fund2mod 0.008969 0.196903 0.046 0.96367
## x.sex1M:y.fav1fav:z.fund1fund 0.438301 0.336151 1.304 0.19227
## x.sex1M:y.fav1fav:z.fund2mod 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: 1.4655e-14 on 0 degrees of freedom
## AIC: 100.14
##
## Number of Fisher Scoring iterations: 3exp(model$coefficients)## (Intercept) x.sex1M
## 70.0000000 0.7000000
## y.fav1fav z.fund1fund
## 1.5857143 1.0428571
## z.fund2mod x.sex1M:y.fav1fav
## 1.5000000 1.5315315
## x.sex1M:z.fund1fund x.sex1M:z.fund2mod
## 0.6262231 0.7619048
## y.fav1fav:z.fund1fund y.fav1fav:z.fund2mod
## 1.0625694 1.0090090
## x.sex1M:y.fav1fav:z.fund1fund x.sex1M:y.fav1fav:z.fund2mod
## 1.5500717 1.3262868Model Homogenous
#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.5744 -0.3002 0.3994 0.0856 -0.1530 -0.0887 0.1129
## 9 10 11 12
## -0.4085 0.6655 0.4342 -0.5281
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.31096 0.10522 40.972 < 2e-16 ***
## x.sex1M -0.51575 0.13814 -3.733 0.000189 ***
## y.fav1fav 0.35707 0.12658 2.821 0.004788 **
## z.fund1fund -0.06762 0.14452 -0.468 0.639854
## z.fund2mod 0.33196 0.13142 2.526 0.011540 *
## x.sex1M:y.fav1fav 0.66406 0.12728 5.217 1.81e-07 ***
## x.sex1M:z.fund1fund -0.16201 0.15300 -1.059 0.289649
## x.sex1M:z.fund2mod -0.08146 0.14079 -0.579 0.562887
## y.fav1fav:z.fund1fund 0.23873 0.16402 1.455 0.145551
## y.fav1fav:z.fund2mod 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: 3Yang digunakan residual deviance.
UJI HIPOTESIS UNTUK MENGETAHUI ADA/ TIDAK HUBUNGAN TIGA ARAH (SATURATED VS HOMOGENOUS)
#pengujian hipotesis
# Deviance of Model
Deviance.model<- model2$deviance -
model$deviance
Deviance.model## [1] 1.797977# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (2 - 0)
derajat.bebas## [1] 2chi.tabel <- qchisq((1-0.05), df=derajat.bebas)
chi.tabel## [1] 5.991465Keputusan <- ifelse(Deviance.model <= chi.tabel,"Terima", "Tolak")
Keputusan## [1] "Terima"Pada taraf nyata 5%, belum cukup bukti untuk menolak 𝐻0 atau dapat dikatakan bahwa tidak ada interaksi tiga arah antara jenis kelamin, fundamentalisme, dan pendapat mengenai hukuman mati.
Ingat, dalam membuat selisih deviance, model yang menjadi pengurang adalah model yang lebih lengkap (parameter yang lebih banyak atau derajat bebasnya lebih kecil).Makin banyak parameter, makin kecil derajat bebas nya karena derajat bebas = banyaknya amatan atau perkalian dimensi tabel kontingensi (misal: 2x2x3 = 12) - banyaknya parameter atau koefisien (lihat di output R ada berapa banyak coefficients nya termasuk intercept) .
B. UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON X)
#Conditional Association 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 -1.11492 0.16142 -0.20684 0.11953 -0.21283 -0.06470 0.08220
## 9 10 11 12
## -0.74698 1.26595 -0.08912 0.11304
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.23495 0.08955 47.293 < 2e-16 ***
## x.sex1M -0.52960 0.13966 -3.792 0.000149 ***
## y.fav1fav 0.48302 0.08075 5.982 2.20e-09 ***
## z.fund1fund 0.07962 0.10309 0.772 0.439916
## z.fund2mod 0.41097 0.09585 4.288 1.81e-05 ***
## x.sex1M:y.fav1fav 0.65845 0.12708 5.181 2.20e-07 ***
## x.sex1M:z.fund1fund -0.12841 0.15109 -0.850 0.395405
## x.sex1M:z.fund2mod -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: 4yang digunakan selanjutnya adalah Residual Deviance.
#pengujian hipotesis
# Deviance of Model
Deviance.model<- model3$deviance -
model2$deviance #model3: conditional on X, model2: homogenous
Deviance.model## [1] 2.132302# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas## [1] 2chi.tabel <- qchisq((1-0.05), df=derajat.bebas)
chi.tabel## [1] 5.991465Keputusan <- ifelse(Deviance.model <= chi.tabel,"Terima", "Tolak")
Keputusan## [1] "Terima"B. UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Y)
#Conditional Association 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.89983 0.14286 0.64384 0.09764 -0.17120 -0.10112 0.12650
## 9 10 11 12
## 0.02418 0.99745 -0.02499 -0.77148
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.33931 0.09919 43.748 < 2e-16 ***
## x.sex1M -0.59345 0.10645 -5.575 2.48e-08 ***
## y.fav1fav 0.37259 0.12438 2.996 0.00274 **
## z.fund1fund -0.12516 0.13389 -0.935 0.34989
## z.fund2mod 0.30228 0.12089 2.500 0.01240 *
## x.sex1M:y.fav1fav 0.65845 0.12708 5.181 2.20e-07 ***
## y.fav1fav:z.fund1fund 0.21254 0.16205 1.312 0.18966
## y.fav1fav:z.fund2mod 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: 4yang digunakan selanjutnya adalah Residual Deviance.
#pengujian hipotesis
# Deviance of Model
Deviance.model<- model4$deviance -
model2$deviance #model4: conditional on Y, model2: homogenous
Deviance.model## [1] 1.122315# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (4 - 2)
derajat.bebas## [1] 2chi.tabel <- qchisq((1-0.05), df=derajat.bebas)
chi.tabel## [1] 5.991465Keputusan <- ifelse(Deviance.model <= chi.tabel,"Terima", "Tolak")
Keputusan## [1] "Terima"C. UJI MODEL INTERAKSI DUA ARAH (HOMOGENOUS VS CONDITIONAL ON Z)
#Conditional Association on Y
#Conditional Association 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 -2.3495 -1.3171 1.9188 1.4595 -2.2965 -1.4130 1.9781
## 9 10 11 12
## 0.7777 -1.1226 -0.7675 1.0321
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.12255 0.10518 39.195 < 2e-16 ***
## x.sex1M -0.07453 0.10713 -0.696 0.487
## y.fav1fav 0.65896 0.11292 5.836 5.36e-09 ***
## z.fund1fund -0.06540 0.15126 -0.432 0.665
## z.fund2mod 0.33196 0.13777 2.410 0.016 *
## x.sex1M:z.fund1fund -0.12841 0.15109 -0.850 0.395
## x.sex1M:z.fund2mod -0.06267 0.13908 -0.451 0.652
## y.fav1fav:z.fund1fund 0.21254 0.16205 1.312 0.190
## y.fav1fav:z.fund2mod 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: 4yang digunakan selanjutnya adalah Residual Deviance.
#pengujian hipotesis
# Deviance of Model
Deviance.model<- model5$deviance -
model2$deviance ##model5: conditional on Z, model2: homogenous
Deviance.model## [1] 27.93095# Chi Square tabel dengan alpa = 0.05
derajat.bebas <- (3 - 2)
derajat.bebas## [1] 1chi.tabel <- qchisq((1-0.05), df=derajat.bebas)
chi.tabel## [1] 3.841459Keputusan <- ifelse(Deviance.model <= chi.tabel,"Terima", "Tolak")
Keputusan## [1] "Tolak"D. PEMILIHAN MODEL TERBAIK
#Model Terbaik
#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.sex1M -0.59345 0.10645 -5.575 2.48e-08 ***
## y.fav1fav 0.48302 0.08075 5.982 2.20e-09 ***
## z.fund1fund 0.01986 0.07533 0.264 0.792
## z.fund2mod 0.38130 0.06944 5.491 4.00e-08 ***
## x.sex1M:y.fav1fav 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: 4Dari summary model diatas terlihat bahwa best model memiliki AIC yang lebih rendah dibandingkan saturated, homogeneous, dan conditional model
E. INTERPRETASI KOEFISIEN MODEL TERBAIK
#interpretasi
data.frame(koef=bestmodel$coefficients,
exp_koef=exp(bestmodel$coefficients))## koef exp_koef
## (Intercept) 4.26517861 71.1776316
## x.sex1M -0.59344782 0.5524194
## y.fav1fav 0.48302334 1.6209677
## z.fund1fund 0.01985881 1.0200573
## z.fund2mod 0.38129767 1.4641834
## x.sex1M:y.fav1fav 0.65845265 1.9318008
F. NILAI DUGAAN MODEL TERBAIK
#fitted
data.frame(Fund=z.fund,sex=x.sex,favor=
y.fav,counts=counts,
fitted=bestmodel$fitted.values)## Fund sex favor counts fitted
## 1 1fund 1M 1fav 128 125.59539
## 2 1fund 1M 2opp 32 40.10855
## 3 1fund 2F 1fav 123 117.69079
## 4 1fund 2F 2opp 73 72.60526
## 5 2mod 1M 1fav 182 180.27878
## 6 2mod 1M 2opp 56 57.57155
## 7 2mod 2F 1fav 168 168.93257
## 8 2mod 2F 2opp 105 104.21711
## 9 3lib 1M 1fav 119 123.12582
## 10 3lib 1M 2opp 49 39.31990
## 11 3lib 2F 1fav 111 115.37664
## 12 3lib 2F 2opp 70 71.17763
Keterangan: nilai miyu_ijk akan sama apapun referensi dari kategori peubahnya yang kita gunakan.
G1501202071, Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, reniamelia@apps.ipb.ac.id↩︎