Homework 5

#Define an ordinal or multinomial outcome variable of your choosing and define how you will recode the original variable:

brfss20sm <- readRDS("C:/Users/shahi/Dropbox/PC/Downloads/brfss20sm.rds")

names(brfss20sm) <- tolower(gsub(pattern = "_",replacement =  "",x =  names(brfss20sm)))

#My binary outcome variable is Poor Mental Health (MENT14D)

brfss20sm$poormentalhlth<-Recode(brfss20sm$ment14d, recodes="1=1; 2=2;3=3; 9=NA",
                             as.factor=T)
brfss20sm$poormentalhlth <- relevel(brfss20sm$poormentalhlth,ref = "1")

###State a research question about what factors you believe will affect your outcome variable:

**How does poor mental health affected by employment, marital status, and age category? Are divorced people suffer poorer mental health than married? Are unable or retired people suffer poorer mental health than employed? Do young adult suffer from poorer mental health than aged people?

#sex
brfss20sm$male<-as.factor(ifelse(brfss20sm$sex==1,
                                "Male",
                                "Female"))

#Age cut into intervals
brfss20sm$agec<-cut(brfss20sm$age80,
                   breaks=c(0,24,39,59,79,99))

#race/ethnicity
brfss20sm$black<-Recode(brfss20sm$racegr3,
                       recodes="2=1; 9=NA; else=0")
brfss20sm$white<-Recode(brfss20sm$racegr3,
                       recodes="1=1; 9=NA; else=0")
brfss20sm$other<-Recode(brfss20sm$racegr3,
                       recodes="3:4=1; 9=NA; else=0")
brfss20sm$hispanic<-Recode(brfss20sm$racegr3,
                          recodes="5=1; 9=NA; else=0")

brfss20sm$race_eth<-Recode(brfss20sm$racegr3,
                          recodes="1='nhwhite'; 2='nh_black'; 3='nh_other';4='nh_multirace'; 5='hispanic'; else=NA",
                          as.factor = T)

#insurance
brfss20sm$ins<-Recode(brfss20sm$hlthpln1,
                     recodes ="7:9=NA; 1=1;2=0")

#income grouping
brfss20sm$inc<-ifelse(brfss20sm$incomg==9,
                     NA,
                     brfss20sm$incomg)

#education level
brfss20sm$educ<-Recode(brfss20sm$educa,
                      recodes="1:2='0Prim'; 3='1somehs'; 4='2hsgrad'; 5='3somecol'; 6='4colgrad';9=NA",
                      as.factor=T)

#brfss20sm$educ<-relevel(brfss20sm$educ, ref='2hsgrad')

#employment
brfss20sm$employ<-Recode(brfss20sm$employ1,
                        recodes="1:2='Employed'; 2:6='nilf'; 7='retired'; 8='unable'; else=NA",
                        as.factor=T)

brfss20sm$employ<-relevel(brfss20sm$employ,
                         ref='Employed')

#marital status
brfss20sm$marst<-Recode(brfss20sm$marital,
                       recodes="1='married'; 2='divorced'; 3='widowed'; 4='separated'; 5='nm';6='cohab'; else=NA",
                       as.factor=T)

brfss20sm$marst<-relevel(brfss20sm$marst,
                        ref='married')



#BMI, in the brfss20a the bmi variable has 2 implied decimal places, so we must divide by 100 to get real bmi's

brfss20sm$bmi<-brfss20sm$bmi5/100

brfss20sm$obese<-ifelse(brfss20sm$bmi>=30,
                       1,
                       0)

brfss20sm$healthnum<-car::Recode(brfss20sm$poormentalhlth,
                                recodes="1=1;2=2;3=3; 9=NA",
                                as.factor = F)

options(survey.lonely.psu = "adjust")


library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:car':
## 
##     recode

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

sub<-brfss20sm%>%
  select(poormentalhlth,healthnum,sex,
      agec, marst, employ, mmsawt, ststr) %>%
  filter( complete.cases(.))

#First we tell R our survey design
options(survey.lonely.psu = "adjust")
des<-svydesign(ids=~1,
               strata=~ststr,
               weights=~mmsawt,
               data =sub )

Fit the ordinal or the multinomial logistic regression models to your outcome:

###present output from the model in terms of odds ratios and confidence intervals for all model parameters in a table:

To fit an ordinal logit to survey data in R, we use the svyolr function in the survey library.

library(gtsummary)

fit.solr1<-svyolr(poormentalhlth~employ+marst+agec,
                  design = des)
fit.solr1%>%
  tbl_regression()

Characteristic	log(OR)1	95% CI1
employ
Employed	—	—
nilf	0.27	0.21, 0.33
retired	-0.06	-0.14, 0.02
unable	1.0	0.94, 1.2
marst
married	—	—
cohab	0.35	0.25, 0.45
divorced	0.46	0.39, 0.53
nm	0.30	0.24, 0.37
separated	0.31	0.17, 0.45
widowed	0.35	0.24, 0.46
agec
(0,24]	—	—
(24,39]	-0.18	-0.26, -0.11
(39,59]	-0.56	-0.65, -0.48
(59,79]	-1.0	-1.1, -0.86
(79,99]	-1.4	-1.6, -1.2
1 OR = Odds Ratio, CI = Confidence Interval

##Here, we can see unable to work people are more likely to report poorer mental health than employed people. Again, divorced people have higher oods of reporting poorer mental health compared to married people. As age increases, the odds of reporting poorer mental health decreases, compared to those are less than 25 years of age.

####If you use the ordinal model, are the assumptions of the proportional odds model met? How did you evaluate the proportional odds assumption?

Now we can “examine” the proportional odds assumption by fitting logits for each change

ex1<-svyglm(I(healthnum>1)~employ+marst+agec,
            design = des,
            family="binomial")

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

ex2<-svyglm(I(healthnum>2)~employ+marst+agec,
            design = des,
            family="binomial")

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

This is NOT a statistical test, just a a rough guide. Here, I plot the coefficients of each model. If the proportional odds is correct, and the assumption is ok, then all these should be “approximately” the same values.

coefs<-data.frame(parms =c( coefficients(ex1)[-1],
                            coefficients(ex2)[-1]),
                  beta = rep(names(coefficients(ex1))[-1], 2 ), 
                  mod = c(rep("I(healthnum>1)",12 ),
                          rep("I(healthnum>2)",12 )))

coefs%>%
  ggplot()+
  geom_point(aes(x=beta, y=parms, color=mod))+
  theme(axis.text.x = element_text(angle = 45))+
  labs(title = "Graphical Summary of Proportional Odds Assumption",
       x= "Beta",
       y= "Estimate")

Based on this, the effects appear to be pretty consistent across transitions.

#Just print the odds ratios, 
round(exp(rbind(coef(ex1)[-1],
                coef(ex2)[-1])),3)

##      employnilf employretired employunable marstcohab marstdivorced marstnm
## [1,]      1.276         0.920        2.287      1.368         1.510   1.336
## [2,]      1.516         1.053        3.952      1.658         1.902   1.508
##      marstseparated marstwidowed agec(24,39] agec(39,59] agec(59,79]
## [1,]          1.216        1.360       0.781       0.516       0.349
## [2,]          1.973        1.683       0.953       0.738       0.526
##      agec(79,99]
## [1,]       0.231
## [2,]       0.357

Yes, the assumptions of the proportional odds model met.It exactly matched with the previous odd ratio.

Poorer mental health is more likely to experience by young , unable to do work and divorced people. The proportional odds assumption