Praktikum 10 - Model Log Linear untuk Data Respon Poisson
Nofrida Elly
5/14/2021
Mengatur Direktori Kerja
setwd ("D:/PASCASARJANA/Semester 3/1. ADK/02. Praktikum/10")Problem 8.1 (Agresti hlm.347)
Input data ke Tabel Kontingensi
x.io <- factor(rep(c("1sup","2opp"),each=2,times=2))
y.ho <- factor(rep(c("1sup","2opp"),times=4))
z.gender <- factor(rep(c("1M", "2F"), each=4))
counts <- c(76,160,6,25,114,181,11,48)
data.frame(x.io, y.ho, z.gender, counts)## x.io y.ho z.gender counts
## 1 1sup 1sup 1M 76
## 2 1sup 2opp 1M 160
## 3 2opp 1sup 1M 6
## 4 2opp 2opp 1M 25
## 5 1sup 1sup 2F 114
## 6 1sup 2opp 2F 181
## 7 2opp 1sup 2F 11
## 8 2opp 2opp 2F 48Menentukan Kategori Referensi
x.io <- relevel(x.io, ref = "1sup")
y.ho <- relevel(y.ho, ref = "1sup")
z.gender <- relevel(z.gender, ref = "1M")- Fit loglinear models (GH, GI), (GH, HI), (GI, HI), and (GH, GI, HI). Show that models that lack the HI term fit poorly.
Untuk melihat model yang fit poorly, dapat dilihat dari nilai AIC yang terkecil dari masing-masing model. Berikut adalah modelnya :
Model (GH, GI) : Conditional on G(z)
model_a <- glm(counts~x.io+y.ho+z.gender+x.io*z.gender+y.ho*z.gender, family = poisson("link"=log))
summary(model_a)##
## Call:
## glm(formula = counts ~ x.io + y.ho + z.gender + x.io * z.gender +
## y.ho * z.gender, family = poisson(link = log))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.4103 -0.2763 -1.2251 0.7402 0.9489 -0.7181 -2.3699 1.5298
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2833 0.1126 38.027 < 2e-16 ***
## x.io2opp -2.0298 0.1910 -10.625 < 2e-16 ***
## y.ho2opp 0.8136 0.1327 6.133 8.63e-10 ***
## z.gender2F 0.3627 0.1458 2.488 0.0128 *
## x.io2opp:z.gender2F 0.4204 0.2384 1.763 0.0778 .
## y.ho2opp:z.gender2F -0.2082 0.1731 -1.203 0.2290
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 445.823 on 7 degrees of freedom
## Residual deviance: 11.666 on 2 degrees of freedom
## AIC: 69.048
##
## Number of Fisher Scoring iterations: 4data.frame(info=x.io, health=y.ho, gender=z.gender, counts=counts, fitted=model_a$fitted.values)## info health gender counts fitted
## 1 1sup 1sup 1M 76 72.479401
## 2 1sup 2opp 1M 160 163.520599
## 3 2opp 1sup 1M 6 9.520599
## 4 2opp 2opp 1M 25 21.479401
## 5 1sup 1sup 2F 114 104.166667
## 6 1sup 2opp 2F 181 190.833333
## 7 2opp 1sup 2F 11 20.833333
## 8 2opp 2opp 2F 48 38.166667Model (GH, HI) : Conditional on H(y)
model_b <- glm(counts~x.io+y.ho+z.gender+x.io*y.ho+y.ho*z.gender, family = poisson("link"=log))
summary(model_b)##
## Call:
## glm(formula = counts ~ x.io + y.ho + z.gender + x.io * y.ho +
## y.ho * z.gender, family = poisson(link = log))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## 0.08450 0.61232 -0.28835 -1.39207 -0.06863 -0.55869 0.22653 1.16424
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.3210 0.1124 38.453 < 2e-16 ***
## x.io2opp -2.4138 0.2532 -9.535 < 2e-16 ***
## y.ho2opp 0.7053 0.1362 5.179 2.23e-07 ***
## z.gender2F 0.4216 0.1421 2.967 0.00301 **
## x.io2opp:y.ho2opp 0.8724 0.2841 3.071 0.00214 **
## y.ho2opp:z.gender2F -0.2082 0.1731 -1.203 0.22903
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 445.8233 on 7 degrees of freedom
## Residual deviance: 4.1267 on 2 degrees of freedom
## AIC: 61.509
##
## Number of Fisher Scoring iterations: 4data.frame(info=x.io, health=y.ho, gender=z.gender, counts=counts, fitted=model_b$fitted.values)## info health gender counts fitted
## 1 1sup 1sup 1M 76 75.26570
## 2 1sup 2opp 1M 160 152.37923
## 3 2opp 1sup 1M 6 6.73430
## 4 2opp 2opp 1M 25 32.62077
## 5 1sup 1sup 2F 114 114.73430
## 6 1sup 2opp 2F 181 188.62077
## 7 2opp 1sup 2F 11 10.26570
## 8 2opp 2opp 2F 48 40.37923Model (GI, HI) -> conditional on I(x)
model_c <- glm(counts~x.io+y.ho+z.gender+x.io*y.ho+x.io*z.gender, family = poisson("link"=log))
summary(model_c)##
## Call:
## glm(formula = counts ~ x.io + y.ho + z.gender + x.io * y.ho +
## x.io * z.gender, family = poisson(link = log))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.93493 0.67971 0.05945 -0.02883 0.81131 -0.61817 -0.04336 0.02087
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.43609 0.08728 50.828 < 2e-16 ***
## x.io2opp -2.66870 0.29595 -9.017 < 2e-16 ***
## y.ho2opp 0.58486 0.09053 6.460 1.04e-10 ***
## z.gender2F 0.22314 0.08733 2.555 0.01062 *
## x.io2opp:y.ho2opp 0.87239 0.28411 3.071 0.00214 **
## x.io2opp:z.gender2F 0.42041 0.23840 1.763 0.07782 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 445.8233 on 7 degrees of freedom
## Residual deviance: 2.3831 on 2 degrees of freedom
## AIC: 59.765
##
## Number of Fisher Scoring iterations: 3data.frame(info=x.io, health=y.ho, gender=z.gender, counts=counts, fitted=model_c$fitted.values)## info health gender counts fitted
## 1 1sup 1sup 1M 76 84.444444
## 2 1sup 2opp 1M 160 151.555556
## 3 2opp 1sup 1M 6 5.855556
## 4 2opp 2opp 1M 25 25.144444
## 5 1sup 1sup 2F 114 105.555556
## 6 1sup 2opp 2F 181 189.444444
## 7 2opp 1sup 2F 11 11.144444
## 8 2opp 2opp 2F 48 47.855556Model (GH, GI, HI)
model_d <- glm(counts~x.io+y.ho+z.gender+z.gender*y.ho+z.gender*x.io+y.ho*x.io, family = poisson("link"=log))
summary(model_d)##
## Call:
## glm(formula = counts ~ x.io + y.ho + z.gender + z.gender * y.ho +
## z.gender * x.io + y.ho * x.io, family = poisson(link = log))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.10362 0.07183 0.39073 -0.17923 0.08516 -0.06730 -0.26626 0.13173
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.3426 0.1120 38.763 < 2e-16 ***
## x.io2opp -2.7147 0.3035 -8.945 < 2e-16 ***
## y.ho2opp 0.7269 0.1353 5.374 7.68e-08 ***
## z.gender2F 0.3856 0.1434 2.689 0.00717 **
## y.ho2opp:z.gender2F -0.2516 0.1749 -1.438 0.15035
## x.io2opp:z.gender2F 0.4636 0.2406 1.927 0.05401 .
## x.io2opp:y.ho2opp 0.8997 0.2852 3.155 0.00160 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 445.82335 on 7 degrees of freedom
## Residual deviance: 0.30072 on 1 degrees of freedom
## AIC: 59.683
##
## Number of Fisher Scoring iterations: 4data.frame(info=x.io, health=y.ho, gender=z.gender, counts=counts, fitted=model_d$fitted.values)## info health gender counts fitted
## 1 1sup 1sup 1M 76 76.90689
## 2 1sup 2opp 1M 160 159.09311
## 3 2opp 1sup 1M 6 5.09311
## 4 2opp 2opp 1M 25 25.90689
## 5 1sup 1sup 2F 114 113.09311
## 6 1sup 2opp 2F 181 181.90689
## 7 2opp 1sup 2F 11 11.90689
## 8 2opp 2opp 2F 48 47.09311Dengan kategori referensi diatas, diperoleh keempat model diatas dengan nilai AIC yang terkecil adalah sebesar 59.683 yaitu model_d (GH, GI, HI), jika dibandingkan dengan ketiga model lainnya. Sehingga dapat dikatakan bahwa model tersebut merupakan model yang paling sesuai diantara model lainnya. Sedangkan model dengan nilai AIC yang terbesar adalah 69.048 yaitu model_a (GH, GI) atau model yang tidak mengandung HI. Sehingga model ini merupakan model yang paling tidak sesuai (fit poorly).
Pemeriksaan Model Terbaik
Dengan cara berikut, juga dapat diperoleh model terbaiknya adalah :
best_model <- glm(counts~x.io+y.ho+z.gender+x.io*y.ho+ x.io*z.gender+y.ho*z.gender,family=poisson("link"=log))
summary(best_model)##
## Call:
## glm(formula = counts ~ x.io + y.ho + z.gender + x.io * y.ho +
## x.io * z.gender + y.ho * z.gender, family = poisson(link = log))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.10362 0.07183 0.39073 -0.17923 0.08516 -0.06730 -0.26626 0.13173
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.3426 0.1120 38.763 < 2e-16 ***
## x.io2opp -2.7147 0.3035 -8.945 < 2e-16 ***
## y.ho2opp 0.7269 0.1353 5.374 7.68e-08 ***
## z.gender2F 0.3856 0.1434 2.689 0.00717 **
## x.io2opp:y.ho2opp 0.8997 0.2852 3.155 0.00160 **
## x.io2opp:z.gender2F 0.4636 0.2406 1.927 0.05401 .
## y.ho2opp:z.gender2F -0.2516 0.1749 -1.438 0.15035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 445.82335 on 7 degrees of freedom
## Residual deviance: 0.30072 on 1 degrees of freedom
## AIC: 59.683
##
## Number of Fisher Scoring iterations: 4- For model (GH, GI, HI), show that 95% wald confidence intervals equal (0.55, 1.10) for the GH conditional odds ratio and (0.99, 2.55) for the GI conditional odds ratio. Is it plausible that gender has no effect on oponion for these issues?
summary(model_d)##
## Call:
## glm(formula = counts ~ x.io + y.ho + z.gender + z.gender * y.ho +
## z.gender * x.io + y.ho * x.io, family = poisson(link = log))
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.10362 0.07183 0.39073 -0.17923 0.08516 -0.06730 -0.26626 0.13173
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.3426 0.1120 38.763 < 2e-16 ***
## x.io2opp -2.7147 0.3035 -8.945 < 2e-16 ***
## y.ho2opp 0.7269 0.1353 5.374 7.68e-08 ***
## z.gender2F 0.3856 0.1434 2.689 0.00717 **
## y.ho2opp:z.gender2F -0.2516 0.1749 -1.438 0.15035
## x.io2opp:z.gender2F 0.4636 0.2406 1.927 0.05401 .
## x.io2opp:y.ho2opp 0.8997 0.2852 3.155 0.00160 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 445.82335 on 7 degrees of freedom
## Residual deviance: 0.30072 on 1 degrees of freedom
## AIC: 59.683
##
## Number of Fisher Scoring iterations: 4Untuk Model GH, GI, HI (dari model_d)
OR=exp(cbind(OR = coef(model_d), confint(model_d)))## Waiting for profiling to be done...OR## OR 2.5 % 97.5 %
## (Intercept) 76.90688954 61.24698047 95.0692783
## x.io2opp 0.06622437 0.03528729 0.1164387
## y.ho2opp 2.06864575 1.59271144 2.7085301
## z.gender2F 1.47051989 1.11251765 1.9537031
## y.ho2opp:z.gender2F 0.77754740 0.55071660 1.0939520
## x.io2opp:z.gender2F 1.58980672 0.99908275 2.5733079
## x.io2opp:y.ho2opp 2.45892936 1.43846301 4.4275459exp(model_d$coefficients[-1]) ## x.io2opp y.ho2opp z.gender2F y.ho2opp:z.gender2F
## 0.06622437 2.06864575 1.47051989 0.77754740
## x.io2opp:z.gender2F x.io2opp:y.ho2opp
## 1.58980672 2.45892936Interpretasi : Dari hasil diatas dapat terlihat bahwa untuk Wald Confidence Intervals untuk Odd Ratio dari :
GH Conditional
(y.ho2opp:z.gender2F) yaitu (0.55071660, 1.0939520) sehingga dapat dikatakan bahwa : dengan tingkat kepercayaan sebesar 95% kita yakin bahwa jika seorang responden berjenis kelamin perempuan, odd mendukung health care costs berada diantara 0.55071660 sampai dengan 1.0939520 kali dibandingkan jika responden berjenis kelamin laki-laki.
GI Conditional
(x.io2opp:z.gender2F) yaitu (0.99908275, 2.5733079) sehingga dapat dikatakan bahwa : dengan tingkat kepercayaan 95%, kita yakin bahwa jika seorang responden berjenis kelamin perempuan, odd mendukung information opinion berada diantara 0.99908275 sampai dengan 2.5733079 kali dibandingkan odd yang sama jika seorang responden berjenis kelamin laki-laki. Kesimpulan : Dari kedua selang kepercayaan yang dihasilkan baik dari model GI dan GH Conditional, mencakup nilai 1 sehingga dapat dikatakan bahwa jenis kelamin seorang responden tidak berpengaruh terhadap opini dalam menangani AIDS melalui information opinion dan health care cost.