Responsi 6 STA543 Analisis Data Kategorik
Model Regresi Logit Binomial-I (Peubah Bebasnya Kontinu)
setwd("D:\\Kuliah S2 IPB\\Bahan Kuliah\\Semester 2 SSD 2020\\STA543 ADK\\Responsi\\R\\UTS\\")Pendahuluan
Regresi logistik biner digunakan untuk memodelkan hubungan antara peubah respon yang terdiri dari dua kategori dengan satu atau lebih peubah penjelas. Peubah penjelasanya bisa berupa data kontinu atau kategorik.
Pendugaan parameter pada model regresi logistik biner menggunakan metode kemungkinan maksimum dan diselesaikan dengan iterasi Newton Raphson. Dalam pengujian parameter model dilakukan dengan likelihood ratio test atau uji G (uji simultan) dan statistik uji Wald (uji parsial).
Soal 1
Input Data
##=====================##
# INPUT DATA
##=====================##
dataku <- read.csv("datacrab.csv",sep=";")
head(dataku)## C S W Wt Sa
## 1 2 3 28.3 3.05 8
## 2 3 3 26.0 2.60 4
## 3 3 3 25.6 2.15 0
## 4 4 2 21.0 1.85 0
## 5 2 3 29.0 3.00 1
## 6 1 2 25.0 2.30 3dim(dataku)## [1] 173 5Keterangan Data:
This data from a study of nesting horseshoe crabs.Each female horseshoe crab had a male crab resident in her nest. The study invertigated factors affecting whether the female crab had any other males, called satellites, residing nearby. Explanatory variables are the female crabs’s color, spine sondition, weight, and carapace width.
Data have 173 obs. of 5 variables:
C: Colour (1. Light medium, 2. Medium, 3. Dark medium, 4. Dark)
S: Spine Condition (1. both good, 2. One worn or broken, 3. Both worn or broken)
W: Carapace Width (cm)
Wt: Weight (kg)
Sa: Number of satellites
Source: Data courtesy of Jane Brockmann, Zoology Departmen, University of Floridak Study described in Ethology 102: 1-21 (1996).
Pendefinisian Peubah
str(dataku) #struktur data## 'data.frame': 173 obs. of 5 variables:
## $ C : int 2 3 3 4 2 1 4 2 2 2 ...
## $ S : int 3 3 3 2 3 2 3 3 1 3 ...
## $ W : num 28.3 26 25.6 21 29 25 26.2 24.9 25.7 27.5 ...
## $ Wt: num 3.05 2.6 2.15 1.85 3 2.3 1.3 2.1 2 3.15 ...
## $ Sa: int 8 4 0 0 1 3 0 0 8 6 ...summary(dataku)## C S W Wt
## Min. :1.000 Min. :1.000 Min. :21.0 Min. :1.200
## 1st Qu.:2.000 1st Qu.:2.000 1st Qu.:24.9 1st Qu.:2.000
## Median :2.000 Median :3.000 Median :26.1 Median :2.350
## Mean :2.439 Mean :2.486 Mean :26.3 Mean :2.437
## 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:27.7 3rd Qu.:2.850
## Max. :4.000 Max. :3.000 Max. :33.5 Max. :5.200
## Sa
## Min. : 0.000
## 1st Qu.: 0.000
## Median : 2.000
## Mean : 2.919
## 3rd Qu.: 5.000
## Max. :15.000##=====================##
# PENDEFINISIAN PEUBAH
##=====================##
dataku$c<- factor(dataku[,1]) #datanya adalah data kategorik, diubah as factor
dataku$s<- factor(dataku[,2]) #datanya adalah data kategorik, diubah as factor
dataku$w<- dataku[,3]
dataku$wt<- dataku[,4]
dataku$sa<- dataku[,5]
#pengkategorian peubah respon
# Y=1= number of satellites (sa)>0 Y=0=sa=0
dataku$y<- c(1:173)
for (i in 1:length(dataku$sa))
{
if(dataku$sa[i]>0)(dataku$y[i]=1)else(dataku$y[i]=0)
}
head(dataku)## C S W Wt Sa c s w wt sa y
## 1 2 3 28.3 3.05 8 2 3 28.3 3.05 8 1
## 2 3 3 26.0 2.60 4 3 3 26.0 2.60 4 1
## 3 3 3 25.6 2.15 0 3 3 25.6 2.15 0 0
## 4 4 2 21.0 1.85 0 4 2 21.0 1.85 0 0
## 5 2 3 29.0 3.00 1 2 3 29.0 3.00 1 1
## 6 1 2 25.0 2.30 3 1 2 25.0 2.30 3 1A. MODEL REGRESI LOGISTIK DENGAN PEUBAH BEBAS WIDTH (w)
##=====================##
# MODEL REGRESI LOGISTIK BINER DENGAN W SEBAGAI PEUBAH BEBAS
##=====================##
model <- glm(y~w, data=dataku, family=binomial("link"=logit))
summary(model)##
## Call:
## glm(formula = y ~ w, family = binomial(link = logit), data = dataku)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.0281 -1.0458 0.5480 0.9066 1.6942
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -12.3508 2.6287 -4.698 2.62e-06 ***
## w 0.4972 0.1017 4.887 1.02e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 225.76 on 172 degrees of freedom
## Residual deviance: 194.45 on 171 degrees of freedom
## AIC: 198.45
##
## Number of Fisher Scoring iterations: 4exp(model$coefficients)## (Intercept) w
## 4.326214e-06 1.644162e+00A. MEMBANDINGKAN OUTPUT R DENGAN OUTPUT SAS
Output SAS:
Hasil dugaan parameter pada kedua output program baik R maupun SAS adalah sama. Namun demikian, pada output R tidak menampilkan hasil uji Wald.
Untuk mendapatkan hasil uji Wald dari hasil output program R kita dapat menghitungnya secara manual menggunakan informasi yang terdapat pada hasil output atau dengan cara menambahkan sintaks sebagai berikut:
##=====================##
# UJI WALD
##=====================##
#install.packages('car') ,jika belum install package car, di install dulu
library('car')## Loading required package: carDataAnova (model,type='II',test='Wald')## Analysis of Deviance Table (Type II tests)
##
## Response: y
## Df Chisq Pr(>Chisq)
## w 1 23.887 1.021e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Hasilnya sama dengan SAS.
A. INTERPRETASI KOEFISIEN REGRESI
##=====================##
# INTERPRETASI KOEFISIEN MODEL
#REGRESI LOGISTIK BINER DENGAN W
# SEBAGAI PEUBAH BEBAS
##=====================##
exp(model$coefficients)## (Intercept) w
## 4.326214e-06 1.644162e+00Interpretasi dugaan parameter model:
-Interpretasi b0 : perbandingan peluang
tanpa memperhatikan lebar cangkang, peluang kepiting betina untuk menarik kepiting jantan adalah 4.32*10^-6 dibandingkan peluang untuk tidak menarik kepiting jantan
-Interpretasi b1 (slope) : perbandingan odds
odds pada kepiting betina yang mampu menarik kepiting jantan akan meningkat sebesar 1.64 kali jika lebar cangkang naik sebesar satu satuan
B. MODEL REGRESI LOGISTIK DENGAN PEUBAH BEBAS WIDTH DAN AKAR WIDTH
##=====================##
# MODEL REGRESI LOGISTIK BINER DENGAN W DAN
#AKAR W SEBAGAI PEUBAH BEBAS
##=====================##
dataku$w1<- sqrt(dataku$w)
model2 <- glm(y~w+w1, data=dataku, family=binomial("link"=logit))
summary(model2)##
## Call:
## glm(formula = y ~ w + w1, family = binomial(link = logit), data = dataku)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.1128 -1.0417 0.5088 0.9451 1.5447
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 89.073 118.713 0.750 0.453
## w 4.435 4.630 0.958 0.338
## w1 -39.991 46.901 -0.853 0.394
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 225.76 on 172 degrees of freedom
## Residual deviance: 193.70 on 170 degrees of freedom
## AIC: 199.7
##
## Number of Fisher Scoring iterations: 5P_Value untuk intercept, w, dan w1 atau akar w > alpha=0.05, berarti tidak tolak Ho. Semuanya tidak nyata atau tidak signifikan pada taraf alpha 0.05.
B. UJI WALD
##=====================##
# UJI WALD
##=====================##
Anova(model2,type='II',test='Wald')## Analysis of Deviance Table (Type II tests)
##
## Response: y
## Df Chisq Pr(>Chisq)
## w 1 0.9176 0.3381
## w1 1 0.7270 0.3938Nilai-p yang diperoleh dari output program R untuk kedua peubah baik width maupun akar width adalah lebih dari 0.05, sehingga tidak tolak H0, dengan kata lain peubah width dan akar width tidak berpengaruh nyata. Hasil dari output program R sama dengan hasil perhitungan manual,
C. PERBANDINGAN MODEL DI POIN A DAN MODEL POIN B
Nilai AIC untuk model (a) 198.45, sedangkan AIC untuk model (b) 199.7.
Dugaan parameter pada model (a) berpengaruh nyata, sedangkan parameter pada model (b) tidak berpengaruh nyata pada α = 5%.
Dengan demikian model terbaik adalah model pada (a).
Soal 2
Input Data
##=====================##
# INPUT DATA
##=====================##
gss<-read.csv("Data gss 2006.csv", sep=";")
dim(gss)## [1] 1834 5head(gss)## No educ sexeduc age sei
## 1 1 14 1 27 50.7
## 2 2 9 2 67 35.7
## 3 3 12 1 50 38.4
## 4 4 16 1 35 54.7
## 5 5 18 1 32 80.3
## 6 6 16 1 60 54.2Cek Missing Data
##=====================##
#cek missing data
##=====================##
apply(gss,2,function(x) sum(is.na(x)))## No educ sexeduc age sei
## 0 0 0 0 0##=====================##Supaya perbandingan model valid maka perlu membuang observasi yang mengandung missing value.Di sini tidak ada missing value.
Pendefinisian Peubah
Variabel respon (sexeduc) yang berkode 1 dan 2 akan di konvert menjadi 1 dan 0 untuk memudahkan proses selanjutnya. Nilai ini tidak perlu dibuat factor karena nanti tidak sesuai di hosmernya (yang dibutuhkan numerik).
##=====================##
#convert y 0,1
##=====================##
gss$sexeduc<-ifelse(gss$sexeduc==1,1,0)
head(gss)## No educ sexeduc age sei
## 1 1 14 1 27 50.7
## 2 2 9 0 67 35.7
## 3 3 12 1 50 38.4
## 4 4 16 1 35 54.7
## 5 5 18 1 32 80.3
## 6 6 16 1 60 54.2A. Lakukan uji rasio likelihood untuk membandingkan model dengan tiga prediktor terhadap model dengan hanya prediktor AGE dan EDUC. Laporkan dan interpretasi (secara substansive) hasil dari pengujian ini
Model 3 Prediktor
#Model 3 predictor
model1<-glm(sexeduc~age+educ+sei, data=gss,
family=binomial("link"=logit))
summary(model1)##
## Call:
## glm(formula = sexeduc ~ age + educ + sei, family = binomial(link = logit),
## data = gss)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4638 0.3757 0.4300 0.4904 0.8376
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.248968 0.402101 5.593 2.23e-08 ***
## age -0.018709 0.004631 -4.040 5.34e-05 ***
## educ 0.042218 0.027922 1.512 0.131
## sei 0.005889 0.004961 1.187 0.235
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1208.1 on 1833 degrees of freedom
## Residual deviance: 1182.3 on 1830 degrees of freedom
## AIC: 1190.3
##
## Number of Fisher Scoring iterations: 5Persamaan yang terbentuk adalah
Yang selanjutnya digunakan adalah Residual Deviance.
Model 2 Prediktor
#Model 2 predictor: Age dan EDUC
model2<-glm(sexeduc~age+educ, data=gss,
family=binomial("link"=logit))
summary(model2)##
## Call:
## glm(formula = sexeduc ~ age + educ, family = binomial(link = logit),
## data = gss)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4250 0.3769 0.4319 0.4899 0.8609
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.255340 0.396930 5.682 1.33e-08 ***
## age -0.018170 0.004607 -3.944 8.00e-05 ***
## educ 0.061286 0.022508 2.723 0.00647 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1208.1 on 1833 degrees of freedom
## Residual deviance: 1183.8 on 1831 degrees of freedom
## AIC: 1189.8
##
## Number of Fisher Scoring iterations: 5Persamaan yang terbentuk
Uji Rasio Likelihood Model 3 Prediktor dan 2 Prediktor
db Deviance model 2 pred = 1834-3 = 1831 (angka 3 diperoleh dari banyaknya parameter model 2 termasuk intercept)
db Deviance model 3 pred = 1834 - 4 = 1830 (angka 3 diperoleh dari banyaknya parameter model 3 termasuk intercept)
model2$deviance-model1$deviance## [1] 1.423566(-2*logLik(model2))-(-2*logLik(model1))## 'log Lik.' 1.423566 (df=3)library("lmtest")## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlrtest(model2,model1)## Likelihood ratio test
##
## Model 1: sexeduc ~ age + educ
## Model 2: sexeduc ~ age + educ + sei
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -591.88
## 2 4 -591.17 1 1.4236 0.2328B. Uji Rasio Likelihood Model 2 Prediktor dan 1 Prediktor
Model 1 Prediktor
#Model 1 predictor: Age
model3<-glm(sexeduc~age, data=gss,
family=binomial("link"=logit))
summary(model3)##
## Call:
## glm(formula = sexeduc ~ age, family = binomial(link = logit),
## data = gss)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3721 0.3859 0.4348 0.4982 0.6549
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.11517 0.24832 12.545 < 2e-16 ***
## age -0.01914 0.00461 -4.153 3.29e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1208.1 on 1833 degrees of freedom
## Residual deviance: 1190.9 on 1832 degrees of freedom
## AIC: 1194.9
##
## Number of Fisher Scoring iterations: 5Persamaan Yang Terbentuk
Uji likelihood rasio test
model3$deviance-model2$deviance## [1] 7.178502(-2*logLik(model3))-(-2*logLik(model2))## 'log Lik.' 7.178502 (df=2)lrtest(model3,model2)## Likelihood ratio test
##
## Model 1: sexeduc ~ age
## Model 2: sexeduc ~ age + educ
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 2 -595.47
## 2 3 -591.88 1 7.1785 0.007378 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1C. Model yang paling sederhana untuk data gss
Berdasarkan hasil dari poin (a) dan (b), model yang paling sederhana untuk data gss adalah model dengan prediktor umur (AGE) dan tahun tertinggi pendidikan yang ditamatkan (EDUC). Hal ini karena indeks sosial ekonomi (SEI) tidak berpengaruh terhadap pendapat mengenai pendidikan sex di sekolah (poin a) sehingga SEI tidak perlu dimasukkan dalam model. Selain itu, tahun tertinggi pendidikan yang ditamatkan berpengaruh nyata terhadap pendapat mengenai pendidikan sex di sekolah (poin b), sehingga perlu ditambahkan dalam model. Jika dilihat dari nilai AIC model, maka model dengan prediktor AGE dan EDUC memiliki nilai AIC terkecil dibandingkan dua model lainnya. Semakin kecil nilai AIC suatu model maka semakin baik model tersebut.
model1$aic #3 peubah## [1] 1190.335model2$aic #2 peubah## [1] 1189.758model3$aic #1 peubah## [1] 1194.937D Interpretasikan secara lengkap nilai, tanda, dan uji signifikansi dari penduga parameter pada model yang paling sederhana.
#interpretasi model yang paling sederhana
variabel<-function(x){
koef<-x$coefficients
se<-summary(x)$coefficients[,2]
wald.stat<-(koef/se)^2
pval<-pchisq(wald.stat,1,lower.tail = F)
exp.koef<-exp(koef)
print(round(data.frame(koef,se,wald.stat,pval,exp.koef),4
))
}
variabel(model2)## koef se wald.stat pval exp.koef
## (Intercept) 2.2553 0.3969 32.2847 0.0000 9.5385
## age -0.0182 0.0046 15.5584 0.0001 0.9820
## educ 0.0613 0.0225 7.4140 0.0065 1.0632Intercept
Interpretasi
Tanpa memperhatikan umur (AGE) dan tahun tertinggi pendidikan yang ditamatkan, peluang responden untuk mendukung pendidikan sex di sekolah adalah 9,534 kali dibandingkan peluang menolak pendidikan sex di sekolah. Dengan kata lain, responden akan lebih cenderung untuk mendukung pendidikan sex di sekolah.
Uji Signifikansi
Age
Interpretasi
Odd mendukung pendidikan sex di sekolah akan meningkat 0.9820 kali jika umur responden bertambah 1 tahun dan tahun tertinggi pendidikan yang ditamatkan tetap.
Tanda negatif pada penduga parameter age menunjukkan adanya hubungan yang berlawanan. Semakin tua umur responden maka akan cenderung untuk menolak pendidikan seks di sekolah (kecenderungan mendukung pendidikan sex di sekolah semakin menurun).
Uji Signifikansi
Educ
Interpretasi
Odd mendukung pendidikan sex di sekolah akan meningkat 1.0623 kali jika tahun pendidikan tertinggi yang ditamatkan responden bertambah 1 tahun dan umur tetap.
Tanda positif pada penduga parameter EDUC menunjukkan adanya hubungan yang searah. Semakin tingi pendidikan tertinggi yang ditamatkan responden maka akan cenderung untuk mendukung pendidikan seks di sekolah (kecenderungan mendukung pendidikan sex di sekolah semakin meningkat).
Uji Signifikansi
E. Ketika sebuah model adalah yang paling parsimoni dibandingkan dengan beberapa model yang diberikan, model tersebut mungkin atau mungkin tidak mem-fitkan nilai pengamatan secara baik. Lakukan uji kecocokan (yang membandingkan nilai pengamatan dan nilai prediksi) untuk model paling parsimoni dan interpretasikan hasilnya.
#hosmer lemeshow test
#install.packages("ResourceSelection")
library(ResourceSelection)## ResourceSelection 0.3-5 2019-07-22hoslem.test(gss$sexeduc,model2$fitted.values)##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: gss$sexeduc, model2$fitted.values
## X-squared = 11.985, df = 8, p-value = 0.1519
F. Berdasarkan Model dengan prediktor hanya AGE dan EDUC, apa dugaan peluang bahwa seorang respondent akan mendukung pendidikan sex di sekolah jika responden tersebut berumur 40 tahun dan telah menyelesaikan 12 tahun pendidikan.
#f. predict
new=data.frame(age=40, educ=12)
predict(model2,newdata=new,type="response")## 1
## 0.9058473Jika seseorang berumur 40 tahun dan telah menempuh 12 tahun pendidikan, maka peluang orang tersebut akan mendukung pendidikan sex di sekolah adalah 0,9058 atau dengan kata lain orang tersebut akan diprediksi untuk mendukung adanya pendidikan seks di sekolah (tergolong kategori 1 dengan cutting point = 0,5)
G1501202071, Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University, reniamelia@apps.ipb.ac.id↩︎