1 Data Management

Load data

pacman::p_load(tidyverse, MASS, HH)

dta3 <- read.table("C:/Users/ASUS/Desktop/data/mentalSES.txt", h=T)
head(dta3)
##   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(dta3)
## '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 ...
dta3_1<- dta3 %>% 
  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')))

2 Analysis

2.1 m1

library(VGAM)
## Warning: package 'VGAM' was built under R version 4.0.3
## 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
m1 <- vglm(ordered(mental) ~ ses, data=dta3_1,family=cumulative(parallel=TRUE))
summary(m1)
## 
## Call:
## vglm(formula = ordered(mental) ~ ses, family = cumulative(parallel = TRUE), 
##     data = dta3_1)
## 
## 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

2.2 m2

m2 <- vglm(ordered(mental) ~ ses, data=dta3_1, family = cumulative)

2.3 Model comparison

VGAM::lrtest(m1, m2)
## 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

2.4 Visdualization

with(dta3_1, ftable(ses, mental))
##     mental Well Mild Moderate Impaired
## ses                                   
## 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
HH::likert(t(with(dta3_1, table(mental, ses))), 
           main="", ylab="SES")