HW4

The data set sim1.out in HW4.RData has two variables, satisfaction with health care (ordinal, 1—5, higher values indicate more satisfaction) and year (0 or 1). A simple random sample was taken from the population in year 0, and a separate simple random sample was taken from that same population in year 1, so you can assume that the two samples are independent. Use a proportional odds model to test whether satisfaction has changed from year zero to year one. Interpret the coefficient estimate associated with year. Make a histogram of satisfaction for each year, and relate it to your coefficient estimate associated with year.

load("/Users/kunaldolas/Desktop/PHC6052/HW4.rdata")
head(sim1.out)

##   satisfaction year
## 1            5    0
## 2            1    0
## 3            5    0
## 4            5    0
## 5            3    0
## 6            3    0

tail(sim1.out)

##     satisfaction year
## 195            4    1
## 196            3    1
## 197            4    1
## 198            3    1
## 199            3    1
## 200            2    1

library(nnet)
mmod<-multinom(satisfaction~year,sim1.out)

## # weights:  15 (8 variable)
## initial  value 321.887582 
## iter  10 value 299.528888
## final  value 299.303678 
## converged

mmodi<-step(mmod)

## Start:  AIC=614.61
## satisfaction ~ year
## 
## trying - year 
## # weights:  10 (4 variable)
## initial  value 321.887582 
## final  value 301.983746 
## converged
##        Df      AIC
## - year  4 611.9675
## <none>  8 614.6074
## # weights:  10 (4 variable)
## initial  value 321.887582 
## final  value 301.983746 
## converged
## 
## Step:  AIC=611.97
## satisfaction ~ 1

sim1.out1<- cbind(ID=1:nrow(sim1.out), sim1.out)
library(repolr)
repomod<-repolr(satisfaction~year, subjects="ID", data=sim1.out1, times=1, categories=5, po.test=TRUE)
summary(repomod)

## 
## repolr: 2016-02-26 version 3.4 
## 
## Call:
## repolr(formula = satisfaction ~ year, subjects = "ID", data = sim1.out1, 
##     times = 1, categories = 5, po.test = TRUE)
## 
## Coefficients: 
##          coeff    se.robust  z.robust  p.value
## cuts1|2  -2.1573   0.2662    -8.1041    0.0000
## cuts2|3  -1.3150   0.2149    -6.1191    0.0000
## cuts3|4   0.1425   0.1901     0.7496    0.4535
## cuts4|5   1.0374   0.2045     5.0729    0.0000
## year     -0.0806   0.2531    -0.3185    0.7501
## 
## PO Score Test: 5.1898 (d.f. = 3 and p.value = 0.1584)

library(MASS)
pomod<-polr(satisfaction~year, sim1.out1)
c(deviance(pomod), pomod$edf)

## [1] 603.8663   5.0000

c(deviance(mmod), mmod$edf)

## [1] 598.6074   8.0000

####### the deviance here is close to the corresponding multinomial logit model 

pomodi<-step(pomod)

## Start:  AIC=613.87
## satisfaction ~ year
## 
##        Df    AIC
## - year  1 611.97
## <none>    613.87
## 
## Step:  AIC=611.97
## satisfaction ~ 1

#######  year does not change the significance variable
deviance(pomodi)-deviance(pomod)

## [1] 0.101201

pchisq(0.101201,pomod$edf-pomodi$edf,lower=F)

## [1] 0.7503931

summary(pomod)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = satisfaction ~ year, data = sim1.out1)
## 
## Coefficients:
##        Value Std. Error t value
## year 0.08066     0.2535  0.3182
## 
## Intercepts:
##     Value   Std. Error t value
## 1|2 -2.1573  0.2667    -8.0881
## 2|3 -1.3150  0.2160    -6.0887
## 3|4  0.1425  0.1946     0.7320
## 4|5  1.0374  0.2087     4.9698
## 
## Residual Deviance: 603.8663 
## AIC: 613.8663

exp(0.08066)

## [1] 1.084002

##### The shift in significance variable from year 0 to 1 is 1.084002
###### not a significant shift in the significance variable from year 0 to 1.

inclevels<-seq(0,1, by=1)
predict(pomod, data.frame(year=inclevels,row.names = inclevels), type="probs")

##            1         2         3         4         5
## 0 0.10365219 0.1079910 0.3239140 0.2027843 0.2616586
## 1 0.09639388 0.1021035 0.3169491 0.2070146 0.2775389

###### Plotting the histogram for the satisfaction across the year 0 and 1
x <- seq(-4,4,by=0.05)
plot(x,dlogis(x),type="l")
abline(v=c(-2.1573,1.0374))
abline(v=c(-2.1573,1.0374)-1*0.08066,lty=2)

## solid lines are for year 0 and dashed lines for year1.