1 load package

pacman::p_load(tidyverse, MASS, HH)

2 data input

dta <- read.table("mentalSES.txt", h = T)
head(dta)
##   Id ses A B C D E mental
## 1  1   6 1 0 0 0 0      1
## 2  2   6 1 0 0 0 0      1
## 3  3   6 1 0 0 0 0      1
## 4  4   6 1 0 0 0 0      1
## 5  5   6 1 0 0 0 0      1
## 6  6   6 1 0 0 0 0      1
str(dta)
## 'data.frame':    1660 obs. of  8 variables:
##  $ Id    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ ses   : int  6 6 6 6 6 6 6 6 6 6 ...
##  $ A     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ B     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ C     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ D     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ E     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ mental: int  1 1 1 1 1 1 1 1 1 1 ...

3 data management

dta_f <- dta %>% 
  mutate(mental=dplyr::recode(mental, 'Well', 'Mild', 'Moderate', 'Impaired')) %>% 
  mutate(mental=factor(mental, levels=c('Well','Mild','Moderate','Impaired')),
        ses=factor(ses, levels=c('1', '2', '3', '4', '5', '6'), labels = c('F', 'E', 'D', 'C', 'B', 'A')))

4 approcah1

4.1 model and comparison

library(VGAM)
## Loading required package: stats4
## Loading required package: splines
## 
## Attaching package: 'VGAM'
## The following object is masked from 'package:HH':
## 
##     logit
## The following object is masked from 'package:tidyr':
## 
##     fill
m0 <- VGAM::vglm(ordered(mental) ~ ses, data=dta_f, family=cumulative(parallel=TRUE))
m1 <- VGAM::vglm(ordered(mental) ~ ses, data=dta_f, family=cumulative)
summary(m0)
## 
## Call:
## VGAM::vglm(formula = ordered(mental) ~ ses, family = cumulative(parallel = TRUE), 
##     data = dta_f)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  -2.0278     0.1346 -15.070  < 2e-16 ***
## (Intercept):2  -0.3285     0.1249  -2.630 0.008531 ** 
## (Intercept):3   0.6802     0.1260   5.398 6.72e-08 ***
## sesE            0.2571     0.1654   1.555 0.120024    
## sesD            0.5248     0.1538   3.413 0.000643 ***
## sesC            0.6157     0.1630   3.776 0.000159 ***
## sesB            0.8408     0.1696   4.958 7.13e-07 ***
## sesA            0.8238     0.1669   4.935 8.03e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<=1]), logitlink(P[Y<=2]), 
## logitlink(P[Y<=3])
## 
## Residual deviance: 4449.382 on 4972 degrees of freedom
## 
## Log-likelihood: -2224.691 on 4972 degrees of freedom
## 
## Number of Fisher scoring iterations: 4 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
##     sesE     sesD     sesC     sesB     sesA 
## 1.293161 1.690184 1.850902 2.318250 2.279240
VGAM::lrtest(m0, m1)
## Likelihood ratio test
## 
## Model 1: ordered(mental) ~ ses
## Model 2: ordered(mental) ~ ses
##    #Df  LogLik  Df  Chisq Pr(>Chisq)
## 1 4972 -2224.7                      
## 2 4962 -2220.8 -10 7.8272     0.6457

4.2 Cumulative of ses

#knitr::kable(apply(with(dta_f, table(mental, ses)), 2, cumsum))
t(with(dta_f, table(mental, ses)))
##    mental
## ses Well Mild Moderate Impaired
##   F   21   71       54       71
##   E   36   97       54       78
##   D   72  141       77       94
##   C   57  105       65       60
##   B   57   94       54       40
##   A   64   94       58       46

4.3 freq figure

HH::likert(t(with(dta_f, table(mental, ses))), main="", ylab="SES")

5 approcah2

dta_f2 <- dta_f %>% count(ses, mental, sort = F)

5.1 likert scale plotting

m <- as.data.frame(matrix(dta_f2$n, 4, 6))

names(m) <- c('F', 'E', 'D', 'C', 'B', 'A')

rownames(m) <- c("Well", "Mild", "Moderate", "Impaired")
#
print(m)
##           F  E   D   C  B  A
## Well     21 36  72  57 57 64
## Mild     71 97 141 105 94 94
## Moderate 54 54  77  65 54 58
## Impaired 71 78  94  60 40 46

5.2 freq figure

##
likert(t(m), as.percent = FALSE, main = "", ylab = "ses")

likert(t(m), as.percent = TRUE, main = "", ylab = "ses")

5.3 model

dta_f2$ses <- relevel(factor(dta_f2$ses), ref="F")
dta_f2$mental <- relevel(factor(dta_f2$mental), ref="Impaired")
str(dta_f2)
## 'data.frame':    24 obs. of  3 variables:
##  $ ses   : Factor w/ 6 levels "F","E","D","C",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ mental: Factor w/ 4 levels "Impaired","Well",..: 2 3 4 1 2 3 4 1 2 3 ...
##  $ n     : int  21 71 54 71 36 97 54 78 72 141 ...
#
summary(m0 <- polr(mental ~ ses, data = dta_f2, weight = n, 
                           method = "probit", Hess = TRUE))
## Call:
## polr(formula = mental ~ ses, data = dta_f2, weights = n, Hess = TRUE, 
##     method = "probit")
## 
## Coefficients:
##         Value Std. Error t value
## sesE -0.02871    0.09963 -0.2882
## sesD  0.01571    0.09226  0.1703
## sesC  0.10411    0.09757  1.0670
## sesB  0.15048    0.10080  1.4928
## sesA  0.12008    0.09931  1.2091
## 
## Intercepts:
##               Value   Std. Error t value
## Impaired|Well -0.6670  0.0770    -8.6663
## Well|Mild     -0.1442  0.0762    -1.8930
## Mild|Moderate  0.8383  0.0772    10.8517
## 
## Residual Deviance: 4482.826 
## AIC: 4498.826
confint(m0)
## Waiting for profiling to be done...
##            2.5 %    97.5 %
## sesE -0.22397655 0.1665702
## sesD -0.16509078 0.1965737
## sesC -0.08708925 0.2953709
## sesB -0.04705474 0.3480986
## sesA -0.07452983 0.3147750