Input Data

##=====================##
# INPUT DATA
##=====================##
dataku<-read.csv("data_mental.csv", sep=";", header=TRUE)
dim(dataku)

## [1] 40  5

head(dataku)

##   Mental_Imp degree degree_2 x2 x1
## 1       well      1        4  1  1
## 2       well      1        4  1  9
## 3       well      1        4  1  4
## 4       well      1        4  1  3
## 5       well      1        4  0  2
## 6       well      1        4  1  0

View(dataku)
mental.degree<-factor(dataku$degree)
ses<-factor(dataku$x2)
life<-dataku$x1
data.frame(mental.degree,ses,life)

##    mental.degree ses life
## 1              1   1    1
## 2              1   1    9
## 3              1   1    4
## 4              1   1    3
## 5              1   0    2
## 6              1   1    0
## 7              1   0    1
## 8              1   1    3
## 9              1   1    3
## 10             1   1    7
## 11             1   0    1
## 12             1   0    2
## 13             2   1    5
## 14             2   0    6
## 15             2   1    3
## 16             2   0    1
## 17             2   1    8
## 18             2   1    2
## 19             2   0    5
## 20             2   1    5
## 21             2   1    9
## 22             2   0    3
## 23             2   1    3
## 24             2   1    1
## 25             3   0    0
## 26             3   1    4
## 27             3   0    3
## 28             3   0    9
## 29             3   1    6
## 30             3   0    4
## 31             3   0    3
## 32             4   1    8
## 33             4   1    2
## 34             4   1    7
## 35             4   0    5
## 36             4   0    4
## 37             4   0    4
## 38             4   1    8
## 39             4   0    8
## 40             4   0    9

Fitting Model

##=====================## 
# FITTING MODEL 
##=====================##
library('foreign')
library('MASS')
library('nnet')
model<-polr(mental.degree~ses+life, method="logistic")
summary(model)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = mental.degree ~ ses + life, method = "logistic")
## 
## Coefficients:
##        Value Std. Error t value
## ses1 -1.1112     0.6109  -1.819
## life  0.3189     0.1210   2.635
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -0.2819  0.6423    -0.4389
## 2|3  1.2128  0.6607     1.8357
## 3|4  2.2094  0.7210     3.0644
## 
## Residual Deviance: 99.0979 
## AIC: 109.0979

A. Bandingkan hasilnya dengan output SAS pada buku Agresti tersebut serta berikan interpretasi pada tiap nilai dugaan parameter model.

##=====================## 
# FITTING MODEL 
##=====================##

model<-polr(mental.degree~ses+life, method="logistic")
summary(model)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = mental.degree ~ ses + life, method = "logistic")
## 
## Coefficients:
##        Value Std. Error t value
## ses1 -1.1112     0.6109  -1.819
## life  0.3189     0.1210   2.635
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -0.2819  0.6423    -0.4389
## 2|3  1.2128  0.6607     1.8357
## 3|4  2.2094  0.7210     3.0644
## 
## Residual Deviance: 99.0979 
## AIC: 109.0979

Output SAS:

##Interpretasi prameter#
coefmodel<-c(-model$coefficients,model$zeta)
data.frame(coefmodel,expcoef=exp(coefmodel)) 

##       coefmodel   expcoef
## ses1  1.1112310 3.0380960
## life -0.3188613 0.7269764
## 1|2  -0.2819031 0.7543468
## 2|3   1.2127926 3.3628626
## 3|4   2.2093721 9.1099947

B. Berdasarkan hasil pada poin (a) diatas, tentukan nilai dugaan 𝑃(𝑌)=1,𝑃(𝑌)=3 , dan 𝑃(𝑌)>2

##=====================## 
# predict model
##=====================## 
new<-data.frame(ses=factor(c(1)),life=4) 
predik<-predict(model, newdata=new, type="probs") #peluang point
peluang<-data.frame(cumulatif=cumsum(predik),point=predik) 
peluang

##   cumulatif     point
## 1 0.3902843 0.3902843
## 2 0.7405018 0.3502175
## 3 0.8854574 0.1449556
## 4 1.0000000 0.1145426

#peluang p(>2)
1-peluang[1,1] 

## [1] 0.6097157

C. Tentukan model terbaik

#model interaksi life*ses
model1<-polr(mental.degree~life*ses,method="logistic") 
summary(model1)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = mental.degree ~ life * ses, method = "logistic")
## 
## Coefficients:
##             Value Std. Error t value
## life       0.4204     0.1864  2.2551
## ses1      -0.3709     1.1361 -0.3265
## life:ses1 -0.1813     0.2383 -0.7608
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2  0.0981  0.8195     0.1197
## 2|3  1.5925  0.8396     1.8968
## 3|4  2.6066  0.9035     2.8849
## 
## Residual Deviance: 98.50444 
## AIC: 110.5044

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

lrtest(model1,model)#cek interaksi

## Likelihood ratio test
## 
## Model 1: mental.degree ~ life * ses
## Model 2: mental.degree ~ ses + life
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   6 -49.252                     
## 2   5 -49.549 -1 0.5935     0.4411

#model hanya peubah life 
model2<-polr(mental.degree~life,method="logistic") 
summary(model2)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = mental.degree ~ life, method = "logistic")
## 
## Coefficients:
##       Value Std. Error t value
## life 0.2879     0.1175   2.451
## 
## Intercepts:
##     Value  Std. Error t value
## 1|2 0.2614 0.5639     0.4636 
## 2|3 1.6563 0.6171     2.6838 
## 3|4 2.5876 0.6938     3.7296 
## 
## Residual Deviance: 102.5271 
## AIC: 110.5271

lrtest (model2,model)#cek ses

## Likelihood ratio test
## 
## Model 1: mental.degree ~ life
## Model 2: mental.degree ~ ses + life
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1   4 -51.264                       
## 2   5 -49.549  1 3.4292    0.06405 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#model hanya peubah ses
model3<-polr(mental.degree~ses,method="logistic") 
summary(model3)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = mental.degree ~ ses, method = "logistic")
## 
## Coefficients:
##        Value Std. Error t value
## ses1 -0.8553     0.5865  -1.458
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -1.3638  0.5038    -2.7070
## 2|3 -0.0417  0.4448    -0.0938
## 3|4  0.8306  0.4659     1.7827
## 
## Residual Deviance: 106.8744 
## AIC: 114.8744

lrtest (model3,model)#cek ses

## Likelihood ratio test
## 
## Model 1: mental.degree ~ ses
## Model 2: mental.degree ~ ses + life
##   #Df  LogLik Df  Chisq Pr(>Chisq)   
## 1   4 -53.437                        
## 2   5 -49.549  1 7.7765   0.005293 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

D. Misalkan seorang individu diketahui bahwa Life Events (x1=8) dan SES (x2=1), berdasarkan model pada poin (c) tentukan dugaan “Mental Impairment”

#prediksi data baru
new=data.frame(life=8,ses=as.factor("1"))
predict(model3,newdata=new,"probs")

##         1         2         3         4 
## 0.3755634 0.3173147 0.1508126 0.1563094

Berdasarkan hasil output R di atas maka dugaan mental impairment-nya adalah well karena memiliki peluang paling besar

A. Persamaan regresi logistic untuk memprediksi odd bahwa seorang anak mencapai second proficiency level (level 1) or less

B. prediksi odd seorang anak laki-laki akan mencapai second proficiency level (level 1) or less dan prediksi odd yang sama untuk perempuan

C. jelaskan bagaimana jawaban anda pada part (a) dan (b) berkaitan dengan interpretasi pendugaan parameter untuk gender

D. ulangi bagian (a), (b) dan (c) untuk odd bahwa seorang anak mencapai first proficiency level (level 0) atau dibawahnya

A.Berapa prediksi peluang kumulatif untuk seorang anak laki-laki mencapai fourth proficiency level (level 3) or less? Berapa prediksi peluang kumulatif untuk seorang anak laki-laki mencapai third proficiency level (level 2) or less?

B. Berapa prediksi peluang seorang anak laki-laki mencapai fourth proficiency level (level 3)? Berapa dugaan peluang yang sama untuk anak perempuan?

C. Berapa prediksi peluang seorang anak laki-laki mencapai the highest proficiency level (level 5)? Berapa dugaan peluang yang sama untuk anak perempuan?