Sex Differences Associated with Opioid Use Disorder, Mental Health Treatment, and Mental Distress: The R Code
Antonia Young
2025-06-05
Read In Data/Packages
load("/Users/antoniayoung/Downloads/NSDUH_2023.Rdata")
library("psych")
library("survey")## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchartlibrary("svydiags")## Loading required package: MASSlibrary("dplyr")##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("tidyr")##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpacklibrary("ggplot2")##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alphalibrary("reghelper")##
## Attaching package: 'reghelper'## The following object is masked from 'package:psych':
##
## ICC## The following object is masked from 'package:base':
##
## betalibrary("ppcor")
library("interactions")
library("gtsummary")##
## Attaching package: 'gtsummary'## The following object is masked from 'package:MASS':
##
## selectlibrary("car")## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:psych':
##
## logitlibrary("lmtest")## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary("ggstatsplot")## You can cite this package as:
## Patil, I. (2021). Visualizations with statistical details: The 'ggstatsplot' approach.
## Journal of Open Source Software, 6(61), 3167, doi:10.21105/joss.03167Public Use Dataset
nsduh <- puf2023_102124
rm(puf2023_102124)Re-Coding Predictors
nsduh$IRSEX <- replace(nsduh$IRSEX, nsduh$IRSEX==1, 0)
nsduh$IRSEX <- replace(nsduh$IRSEX, nsduh$IRSEX==2, 1)
nsduh$IRHERFM[nsduh$IRHERFM==93] <- 0 #Represents lifetime use
nsduh$IRHERFM[nsduh$IRHERFM==1.5] <- 1
nsduh$IRHERFM[nsduh$IRHERFM==91] <- NA
nsduh$SVYROPIANY[nsduh$SVYROPIANY==4] <- 0
nsduh$IRPNRNM30FQ[nsduh$IRPNRNM30FQ==91] <- NA
nsduh$IRPNRNM30FQ[nsduh$IRPNRNM30FQ==93] <- 0
nsduh <- nsduh %>%
mutate(opioid30 = rowSums(cbind(IRHERFM, IRPNRNM30FQ), na.rm=TRUE))
nsduh <- nsduh %>%
mutate(mhtused = IRMHTOUTDOC + IRMHTOUTHOSP + IRMHTOUTMHCR + IRMHTOUTRHAB + IRMHTOUTSCHL + IRMHTOUTTHRP + MHTOUTOTPY)Correlation Matrix
nsduhcor <- subset(nsduh,select=c(SVYROPIANY,opioid30,mhtused,IRSEX,KSSLR6MONED), !is.na(KSSLR6MONED) & SVYROPIANY!=0)
nsduhcor$SVYROPIANY <- as.numeric(nsduhcor$SVYROPIANY)
psych::corr.test(nsduhcor)## Call:psych::corr.test(x = nsduhcor)
## Correlation matrix
## SVYROPIANY opioid30 mhtused IRSEX KSSLR6MONED
## SVYROPIANY 1.00 0.39 0.14 -0.02 0.25
## opioid30 0.39 1.00 -0.03 -0.05 0.02
## mhtused 0.14 -0.03 1.00 0.12 0.26
## IRSEX -0.02 -0.05 0.12 1.00 0.11
## KSSLR6MONED 0.25 0.02 0.26 0.11 1.00
## Sample Size
## [1] 813
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## SVYROPIANY opioid30 mhtused IRSEX KSSLR6MONED
## SVYROPIANY 0.0 0.00 0 1.00 0.00
## opioid30 0.0 0.00 1 0.59 1.00
## mhtused 0.0 0.35 0 0.01 0.00
## IRSEX 0.6 0.15 0 0.00 0.01
## KSSLR6MONED 0.0 0.53 0 0.00 0.00
##
## To see confidence intervals of the correlations, print with the short=FALSE optiontable(nsduhcor$SVYROPIANY)##
## 1 2 3
## 469 140 204Assigning Weights and Centering
nsduh$MildOUD <- ifelse(nsduh$SVYROPIANY=="1",1,
ifelse(nsduh$SVYROPIANY=="2",-140/469, -204/469))
nsduh$ModerateOUD <- ifelse(nsduh$SVYROPIANY=="2",1,0)
nsduh$SevereOUD <- ifelse(nsduh$SVYROPIANY=="3",1,0)
nsduh$opioid30C <- scale(nsduh$opioid30,center=TRUE,scale=FALSE)
nsduh$mhtC <- scale(nsduh$mhtused,center=TRUE,scale=FALSE)
nsduh$IRSEXC <- scale(nsduh$IRSEX,center=TRUE,scale=FALSE)Setting NSDUH Weights
options(survey.lonely.psu="adjust")
options(survey.adjust.domain.lonely=TRUE)
design.nsduh2 <- svydesign(strata=~VESTR_C,id=~VEREP,weights=~ANALWT2_C,data=nsduh,nest=TRUE)
design.nsduh <- subset(design.nsduh2, !is.na(KSSLR6MONED) & SVYROPIANY!=0)K6 Full Dataset and Subsample Reliability and CFA
nsduhK6 <- subset(design.nsduh2$variables,select=c(DSTCHR30, DSTEFF30, DSTHOP30, DSTNGD30, DSTNRV30, DSTRST30))
psych::alpha(nsduhK6)## Number of categories should be increased in order to count frequencies.##
## Reliability analysis
## Call: psych::alpha(x = nsduhK6)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 1 1 1 0.99 960 6.9e-06 25 39 0.99
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 1 1 1
## Duhachek 1 1 1
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## DSTCHR30 1 1 1 0.99 725 9.3e-06 4.0e-06 0.99
## DSTEFF30 1 1 1 0.99 976 6.9e-06 9.0e-07 0.99
## DSTHOP30 1 1 1 0.99 792 8.6e-06 3.8e-06 0.99
## DSTNGD30 1 1 1 0.99 776 8.8e-06 5.2e-06 0.99
## DSTNRV30 1 1 1 0.99 817 8.3e-06 3.9e-06 0.99
## DSTRST30 1 1 1 0.99 755 9.0e-06 4.6e-06 0.99
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## DSTCHR30 56705 1 1 1.00 1.00 25 39
## DSTEFF30 56705 1 1 0.99 0.99 25 40
## DSTHOP30 56705 1 1 1.00 1.00 25 39
## DSTNGD30 56705 1 1 1.00 1.00 25 39
## DSTNRV30 56705 1 1 1.00 1.00 25 40
## DSTRST30 56705 1 1 1.00 1.00 25 39psych::omega(nsduhK6, title = "K6 Scale CFA")## Loading required namespace: GPArotation
## K6 Scale CFA
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 1
## G.6: 1
## Omega Hierarchical: 0.16
## Omega H asymptotic: 0.16
## Omega Total 1
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 h2 u2 p2 com
## DSTCHR30 0.40 0.92 1.00 1.00 0.00 0.16 1.36
## DSTEFF30 0.40 0.91 0.99 0.99 0.01 0.16 1.38
## DSTHOP30 0.38 0.92 0.99 0.99 0.01 0.15 1.34
## DSTNGD30 0.40 0.91 1.00 1.00 0.00 0.16 1.38
## DSTNRV30 0.38 0.92 0.99 0.99 0.01 0.15 1.34
## DSTRST30 0.39 0.92 1.00 1.00 0.00 0.15 1.35
##
## With Sums of squares of:
## g F1* F2* F3* h2
## 0.93 5.03 0.00 0.00 5.93
##
## general/max 0.16 max/min = Inf
## mean percent general = 0.16 with sd = 0.01 and cv of 0.04
## Explained Common Variance of the general factor = 0.16
##
## The degrees of freedom are 0 and the fit is 0.05
## The number of observations was 56705 with Chi Square = 2582.56 with prob < NA
## The root mean square of the residuals is 0
## The df corrected root mean square of the residuals is NA
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9 and the fit is 21.29
## The number of observations was 56705 with Chi Square = 1206944 with prob < 0
## The root mean square of the residuals is 0.84
## The df corrected root mean square of the residuals is 1.08
##
## RMSEA index = 1.538 and the 10 % confidence intervals are 1.536 1.54
## BIC = 1206846
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.44 0.92 0.52 0
## Multiple R square of scores with factors 0.19 0.85 0.27 0
## Minimum correlation of factor score estimates -0.62 0.70 -0.46 -1
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 1.00 1.00 NA NA
## Omega general for total scores and subscales 0.16 0.16 NA NA
## Omega group for total scores and subscales 0.84 0.84 NA NAnsduhK62 <- subset(design.nsduh$variables,select=c(DSTCHR30, DSTEFF30, DSTHOP30, DSTNGD30, DSTNRV30, DSTRST30))
psych::alpha(nsduhK62)##
## Reliability analysis
## Call: psych::alpha(x = nsduhK62)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.9 0.91 0.9 0.61 9.5 0.0052 3.3 1.1 0.61
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.89 0.9 0.91
## Duhachek 0.89 0.9 0.92
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## DSTCHR30 0.88 0.88 0.86 0.59 7.1 0.0069 0.0042 0.61
## DSTEFF30 0.90 0.90 0.89 0.64 8.7 0.0057 0.0063 0.61
## DSTHOP30 0.88 0.88 0.87 0.60 7.4 0.0066 0.0068 0.61
## DSTNGD30 0.88 0.88 0.87 0.60 7.4 0.0066 0.0064 0.61
## DSTNRV30 0.90 0.90 0.88 0.65 9.1 0.0055 0.0052 0.64
## DSTRST30 0.89 0.89 0.88 0.62 8.1 0.0061 0.0098 0.63
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## DSTCHR30 813 0.88 0.88 0.86 0.82 3.5 1.3
## DSTEFF30 813 0.78 0.78 0.71 0.68 3.0 1.3
## DSTHOP30 813 0.85 0.86 0.83 0.78 3.5 1.3
## DSTNGD30 813 0.86 0.86 0.83 0.78 3.4 1.3
## DSTNRV30 813 0.76 0.76 0.69 0.65 3.1 1.3
## DSTRST30 813 0.81 0.82 0.76 0.73 3.2 1.3
##
## Non missing response frequency for each item
## 1 2 3 4 5 miss
## DSTCHR30 0.10 0.13 0.22 0.25 0.30 0
## DSTEFF30 0.17 0.20 0.25 0.20 0.18 0
## DSTHOP30 0.09 0.12 0.25 0.24 0.29 0
## DSTNGD30 0.12 0.15 0.23 0.21 0.29 0
## DSTNRV30 0.14 0.16 0.28 0.25 0.17 0
## DSTRST30 0.13 0.14 0.29 0.23 0.20 0psych::omega(nsduhK62, title = "K6 Scale Subsample CFA")
## K6 Scale Subsample CFA
## Call: omegah(m = m, nfactors = nfactors, fm = fm, key = key, flip = flip,
## digits = digits, title = title, sl = sl, labels = labels,
## plot = plot, n.obs = n.obs, rotate = rotate, Phi = Phi, option = option,
## covar = covar)
## Alpha: 0.91
## G.6: 0.9
## Omega Hierarchical: 0.87
## Omega H asymptotic: 0.93
## Omega Total 0.94
##
## Schmid Leiman Factor loadings greater than 0.2
## g F1* F2* F3* h2 h2 u2 p2 com
## DSTCHR30 0.91 0.83 0.83 0.17 1.00 1.01
## DSTEFF30 0.68 0.45 0.68 0.68 0.32 0.69 1.75
## DSTHOP30 0.84 0.72 0.72 0.28 0.98 1.04
## DSTNGD30 0.83 0.71 0.71 0.29 0.96 1.07
## DSTNRV30 0.63 0.77 1.00 1.00 0.00 0.40 1.92
## DSTRST30 0.71 0.21 0.57 0.57 0.43 0.90 1.22
##
## With Sums of squares of:
## g F1* F2* F3* h2
## 3.60 0.00 0.66 0.25 3.49
##
## general/max 1.03 max/min = Inf
## mean percent general = 0.82 with sd = 0.24 and cv of 0.29
## Explained Common Variance of the general factor = 0.8
##
## The degrees of freedom are 0 and the fit is 0
## The number of observations was 813 with Chi Square = 1.16 with prob < NA
## The root mean square of the residuals is 0
## The df corrected root mean square of the residuals is NA
##
## Compare this with the adequacy of just a general factor and no group factors
## The degrees of freedom for just the general factor are 9 and the fit is 0.18
## The number of observations was 813 with Chi Square = 142.55 with prob < 3.1e-26
## The root mean square of the residuals is 0.06
## The df corrected root mean square of the residuals is 0.07
##
## RMSEA index = 0.135 and the 10 % confidence intervals are 0.116 0.155
## BIC = 82.25
##
## Measures of factor score adequacy
## g F1* F2* F3*
## Correlation of scores with factors 0.96 0 0.97 0.61
## Multiple R square of scores with factors 0.92 0 0.95 0.38
## Minimum correlation of factor score estimates 0.84 -1 0.90 -0.25
##
## Total, General and Subset omega for each subset
## g F1* F2* F3*
## Omega total for total scores and subscales 0.94 NA 0.92 0.81
## Omega general for total scores and subscales 0.87 NA 0.83 0.70
## Omega group for total scores and subscales 0.06 NA 0.09 0.11Weighted Data Distributions
layout(matrix(c(1,2,3,4,5,6),nrow=2,ncol=3))
svyhist(~opioid30,design.nsduh,ylim=c(0,0.2),main="Opioid misuse frequency",xlab="Times used")
svyboxplot(~opioid30~1,design.nsduh,all.outliers=T)
svyhist(~mhtused,design.nsduh,ylim=c(0,1.5),main="Number of MHTs used",xlab="Quantity")
svyboxplot(~mhtused~1,design.nsduh,all.outliers=T)
svyhist(~KSSLR6MONED,design.nsduh,ylim=c(0,0.1),xlim=c(0,25),main="K6 scores (past month)",xlab="Total")
svyboxplot(~KSSLR6MONED~1,design.nsduh,all.outliers=T,ylim=c(0,25))
svyciprop(~opioid30==0,design.nsduh)## 2.5% 97.5%
## opioid30 == 0 0.812 0.771 0.847Testing the Moderated Multiple Linear Regression Model
mod <- svyglm(KSSLR6MONED ~ SevereOUD + ModerateOUD + MildOUD + opioid30C + mhtC + IRSEXC + SevereOUD*IRSEXC + ModerateOUD*IRSEXC,family=gaussian(),design=design.nsduh)
summary(mod)##
## Call:
## svyglm(formula = KSSLR6MONED ~ SevereOUD + ModerateOUD + MildOUD +
## opioid30C + mhtC + IRSEXC + SevereOUD * IRSEXC + ModerateOUD *
## IRSEXC, design = design.nsduh, family = gaussian())
##
## Survey design:
## subset(design.nsduh2, !is.na(KSSLR6MONED) & SVYROPIANY != 0)
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.14101 0.37945 18.819 < 2e-16 ***
## SevereOUD 3.08163 0.99713 3.091 0.003499 **
## ModerateOUD 0.50282 0.82179 0.612 0.543851
## opioid30C -0.02770 0.04653 -0.595 0.554751
## mhtC 1.07917 0.28183 3.829 0.000413 ***
## IRSEXC -0.59284 0.60345 -0.982 0.331388
## SevereOUD:IRSEXC 2.93149 1.74489 1.680 0.100202
## ModerateOUD:IRSEXC 6.57379 1.70547 3.855 0.000382 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 31.31447)
##
## Number of Fisher Scoring iterations: 2confint(mod)## 2.5 % 97.5 %
## (Intercept) 6.3757754 7.90624934
## SevereOUD 1.0707269 5.09253350
## ModerateOUD -1.1544813 2.16012837
## opioid30C -0.1215250 0.06613002
## mhtC 0.5107979 1.64754444
## IRSEXC -1.8098184 0.62413132
## SevereOUD:IRSEXC -0.5874242 6.45040270
## ModerateOUD:IRSEXC 3.1343716 10.01320232The svyglm() function does not provide R-squared.
GG Plots with Continuous Covariates
ggplot(design.nsduh$variables, aes(x = opioid30, y = mod$fitted.values)) + geom_point() + geom_smooth(method="lm")## `geom_smooth()` using formula = 'y ~ x'
ggplot(design.nsduh$variables, aes(x = mhtused, y = KSSLR6MONED)) + geom_point(alpha=0.7) + geom_smooth(method="lm") + geom_point(aes(x=mhtused,y=mod$fitted.values), color='blue') + geom_segment(aes(xend=mhtused,yend=mod$fitted.values), color='red', linetype='dashed') + ggtitle("K6 Scores by MHTs Used")## `geom_smooth()` using formula = 'y ~ x'
ggplot(design.nsduh$variables, aes(x = opioid30, y = KSSLR6MONED)) + geom_point(alpha=0.7) + geom_smooth(method="lm") + geom_point(aes(x=opioid30,y=mod$fitted.values), color='blue') + geom_segment(aes(xend=opioid30,yend=mod$fitted.values), color='red', linetype='dashed') + ggtitle("K6 Scores by Opioid Misuse Frequency")## `geom_smooth()` using formula = 'y ~ x'
Testing Significant Interaction
mod_2 <- svyglm(KSSLR6MONED ~ SevereOUD + ModerateOUD + MildOUD + opioid30 + mhtused + IRSEX + SevereOUD*IRSEX + ModerateOUD*IRSEX,family=gaussian(),design=design.nsduh)
sim_slopes(mod_2, pred=ModerateOUD, modx=IRSEX, jnplot=TRUE)
sim_slopes(mod_2, pred=IRSEX, modx=ModerateOUD, jnplot=TRUE)
design.nsduh$variables$OUDSeverity <- factor(design.nsduh$variables$MildOUD,
levels = c(1, -140/469, -204/469),
labels = c("Mild OUD", "Moderate OUD", "Severe OUD"))
grouped_ggbetweenstats(
data = design.nsduh$variables,
x = sex,
y = KSSLR6MONED,
grouping.var = OUDSeverity,
xlab = "Sex",
ylab = "K6 Distress Total",
results.subtitle = TRUE
)Calculating R-Squared and Adjusted R-Squared
rss <- sum((design.nsduh$variables$KSSLR6MONED - mod$fitted.values)^2)
tss <- sum((design.nsduh$variables$KSSLR6MONED - svymean(~KSSLR6MONED,design.nsduh))^2)
r2 <- 1 - rss/tss
r2## [1] 0.1302239aov(mod)## Call:
## aov(formula = mod)
##
## Terms:
## SevereOUD ModerateOUD opioid30C mhtC IRSEXC
## Sum of Squares 1478.726 1.407 90.854 1764.794 213.600
## Deg. of Freedom 1 1 1 1 1
## SevereOUD:IRSEXC ModerateOUD:IRSEXC Residuals
## Sum of Squares 89.072 1086.962 25427.350
## Deg. of Freedom 1 1 805
##
## Residual standard error: 5.620211
## 1 out of 9 effects not estimable
## Estimated effects may be unbalanced1-(1-r2)*(812/805)## [1] 0.1226606Semi-Partial Correlation Coefficients
design.nsduh$variables$moudsex <- design.nsduh$variables$ModerateOUD * design.nsduh$variables$IRSEX
design.nsduh$variables$soudsex <- design.nsduh$variables$SevereOUD * design.nsduh$variables$IRSEX
spcorsK6 <- spcor(design.nsduh$variables[,c("SevereOUD","ModerateOUD","opioid30","mhtused","IRSEX","moudsex","soudsex","KSSLR6MONED")])$estimate
SOUD <- spcorsK6["KSSLR6MONED","SevereOUD"]
MOUD <- spcorsK6["KSSLR6MONED","ModerateOUD"]
Use <- spcorsK6["KSSLR6MONED","opioid30"]
MHT <- spcorsK6["KSSLR6MONED","mhtused"]
Sex <- spcorsK6["KSSLR6MONED","IRSEX"]
MOUDSex <- spcorsK6["KSSLR6MONED","moudsex"]
SOUDSex <- spcorsK6["KSSLR6MONED","soudsex"]
SOUD^2## [1] 0.02261588MOUD^2## [1] 7.053011e-06Use^2## [1] 0.003675568MHT^2## [1] 0.04324362Sex^2## [1] 0.001242361MOUDSex^2## [1] 0.004743117SOUDSex^2## [1] 0.0001079587SOUD^2 + MOUD^2 + Use^2 + MHT^2 + Sex^2 + MOUDSex^2 + SOUDSex^2## [1] 0.075635550.1341447 - 0.07563555## [1] 0.05850915Testing Assumptions of Linear Regression
plot(mod$fitted.values,mod$residuals,ylab="Residual K6 Points",xlab="Predicted K6 Total",xlim = c(0, 25),main="K6 Model Fit")
abline(0,0)
qqnorm(mod$residuals, main = "Q-Q Plot of K6 Residuals")
qqline(mod$residuals, col = "red")
shapiro.test(mod$residuals)##
## Shapiro-Wilk normality test
##
## data: mod$residuals
## W = 0.98402, p-value = 9.467e-08vif(svyglm(KSSLR6MONED ~ SevereOUD + ModerateOUD + opioid30C + mhtC + IRSEXC,family=gaussian(),design=design.nsduh))## SevereOUD ModerateOUD opioid30C mhtC IRSEXC
## 1.202237 1.163058 1.115007 1.153922 1.0954021/vif(svyglm(KSSLR6MONED ~ SevereOUD + ModerateOUD + opioid30C + mhtC + IRSEXC,family=gaussian(),design=design.nsduh))## SevereOUD ModerateOUD opioid30C mhtC IRSEXC
## 0.8317830 0.8598026 0.8968555 0.8666098 0.9129065Removing and Testing Outliers
cooks_dist <- cooks.distance(mod)
cooks_dist
threshold <- 4 / length(cooks_dist)
outliers <- which(cooks_dist > threshold)
outliers
mod2 <- svyglm(KSSLR6MONED ~ SevereOUD + ModerateOUD + MildOUD + opioid30C + mhtC + IRSEXC + SevereOUD*IRSEXC + ModerateOUD*IRSEXC,family=gaussian(),design=design.nsduh,subset=-c(1461,5411,5753,6582,6749,8038,10186,12213,12649,13272,13797,17549,17847,18310,18949,19839,22436,22823,23878,24253,27421,29262,31974,33832,34185,34544,34744,35049,36437,38685,38879,39068,41482,42007,44449,45505,45763,46607,47952,51304,51977,52002,53394,54066,54843,56502,56550))
summary(mod2)
plot(mod2$fitted.values,mod2$residuals,ylab="Residual K6 Points",xlab="Predicted K6 Total",xlim = c(0, 25),main="K6 Model Fit")
abline(0,0)
qqnorm(mod2$residuals, main = "Q-Q Plot of K6 Residuals")
qqline(mod2$residuals, col = "red")Demographics
nsduh$SEXIDENT22 <- replace(nsduh$SEXIDENT22, nsduh$SEXIDENT22 > 5, NA)
nsduh$IRMARIT <- replace(nsduh$IRMARIT, nsduh$IRMARIT > 4, NA)
nsduh$SPEAKENGL <- replace(nsduh$SPEAKENGL, nsduh$SPEAKENGL > 4, NA)
nsduh$WRKSTATWK2 <- replace(nsduh$WRKSTATWK2, nsduh$WRKSTATWK2 > 9, NA)
nsduh$PNRNMLAS1 <- replace(nsduh$PNRNMLAS1, nsduh$PNRNMLAS1 > 37, 0)
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==2] <- 1 #Hydrocodone
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==3] <- 1
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==4] <- 1
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==5] <- 1
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==6] <- 2 #Oxycodone
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==7] <- 2
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==8] <- 2
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==9] <- 2
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==10] <- 2
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==11] <- 3 #Tramadol
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==13] <- 3
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==14] <- 3
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==15] <- 3
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==16] <- 4 #Codeine
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==17] <- 4
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==20] <- 5 #Morphine
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==21] <- 5
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==22] <- 5
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==23] <- 6 #Fentanyl
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==24] <- 6
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==25] <- 6
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==26] <- 7 #Buprenorphine
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==27] <- 7
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==28] <- 7
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==29] <- 8 #Oxymorphone
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==30] <- 8
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==31] <- 8
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==32] <- 8
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==33] <- 9 #Demerol
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==34] <- 10 #Hydromorphone
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==35] <- 10
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==36] <- 11 #Methadone
nsduh$PNRNMLAS1[nsduh$PNRNMLAS1==37] <- 12 #Other
nsduh$CATAG6 <- as.factor(nsduh$CATAG6)
nsduh$SEXRACE <- as.factor(nsduh$SEXRACE)
nsduh$IRMARIT <- as.factor(nsduh$IRMARIT)
nsduh$IRPINC3 <- as.factor(nsduh$IRPINC3)
nsduh$herused <- ifelse(nsduh$HERREC==1,1,0)
nsduh$herused <- as.factor(nsduh$herused)
nsduh$pnrnused <- ifelse(nsduh$IRPNRNMREC==1,1,0)
nsduh$pnrnused <- as.factor(nsduh$pnrnused)
nsduh$SEXIDENT22 <- as.factor(nsduh$SEXIDENT22)
nsduh$EDUHIGHCAT <- as.factor(nsduh$EDUHIGHCAT)
nsduh$SPEAKENGL <- as.factor(nsduh$SPEAKENGL)
nsduh$WRKSTATWK2 <- as.factor(nsduh$WRKSTATWK2)
nsduh$HEALTH2 <- as.factor(nsduh$HEALTH2)
nsduh$SVYROPIANY <- as.factor(nsduh$SVYROPIANY)
nsduh$SUTOUTOPIPY <- as.factor(nsduh$SUTOUTOPIPY)
nsduh$PNRNMLAS1 <- as.factor(nsduh$PNRNMLAS1)
nsduh$IRMHTOUTDOC <- as.factor(nsduh$IRMHTOUTDOC)
nsduh$IRMHTOUTHOSP <- as.factor(nsduh$IRMHTOUTHOSP)
nsduh$IRMHTOUTMHCR <- as.factor(nsduh$IRMHTOUTMHCR)
nsduh$IRMHTOUTRHAB <- as.factor(nsduh$IRMHTOUTRHAB)
nsduh$IRMHTOUTSCHL <- as.factor(nsduh$IRMHTOUTSCHL)
nsduh$IRMHTOUTTHRP <- as.factor(nsduh$IRMHTOUTTHRP)
nsduh$MHTOUTOTPY <- as.factor(nsduh$MHTOUTOTPY)
nsduh$MHTINPPY <- as.factor(nsduh$MHTINPPY)Demographic Information
DemoTable <- design.nsduh %>%
tbl_svysummary(
by = IRSEX,
include = c(
CATAG6, NEWRACE2, IRMARIT, IRPINC3, SEXIDENT22, EDUHIGHCAT, SPEAKENGL,
WRKSTATWK2, HEALTH2, SVYROPIANY, SUTOUTOPIPY, PNRNMLAS1, herused,
IRHERFM, pnrnused, IRPNRNM30FQ, mhtused, MHTINPPY
),
statistic = list(
all_categorical() ~ "{n} ({p}%)",
all_continuous() ~ "{mean} ({sd})"
),
digits = list(
all_categorical() ~ c(0, 2),
all_continuous() ~ c(1, 2)
),
type = list(
mhtused ~ "continuous"
),
label = list(
CATAG6 ~ "Age group",
NEWRACE2 ~ "Racial group",
IRMARIT ~ "Marital status",
IRPINC3 ~ "Income group",
SEXIDENT22 ~ "Sexual orientation",
EDUHIGHCAT ~ "Education status",
SPEAKENGL ~ "How well do you speak English?",
WRKSTATWK2 ~ "Work status (past week)",
HEALTH2 ~ "Perceived health status",
SVYROPIANY ~ "OUD severity",
SUTOUTOPIPY ~ "Used OUD outpatient treatment (past year)",
PNRNMLAS1 ~ "Pain reliever misused (past year)",
herused ~ "Used heroin in past month (yes/no)",
IRHERFM ~ "Mean days used heroin",
pnrnused ~ "Misused pain reliever in past month (yes/no)",
IRPNRNM30FQ ~ "Mean days misused pain reliever",
mhtused ~ "Mean number of outpatient MHTs used",
MHTINPPY ~ "Used inpatient MHT (past year)"
)
) %>%
modify_header(label = "**Variable**") %>%
modify_caption("Weighted descriptive statistics") %>%
bold_labels()Sex Differences in Demographics
svychisq(~CATAG6 + IRSEX, statistic="adjWald", design = design.nsduh)
svychisq(~NEWRACE2 + IRSEX, statistic="F", design = design.nsduh)
svychisq(~IRMARIT + IRSEX, statistic="F", design = design.nsduh)
svychisq(~IRPINC3 + IRSEX, statistic="F", design = design.nsduh)
svychisq(~SEXIDENT22 + IRSEX, statistic="F", design = design.nsduh)
svychisq(~EDUHIGHCAT + IRSEX, statistic="adjWald", design = design.nsduh)
svychisq(~SPEAKENGL + IRSEX, statistic="F", design = design.nsduh)
svychisq(~WRKSTATWK2 + IRSEX, statistic="F", design = design.nsduh)
svychisq(~HEALTH2 + IRSEX, statistic="F", design = design.nsduh)
svychisq(~SVYROPIANY + IRSEX, statistic="adjWald", design = design.nsduh)
svychisq(~SUTOUTOPIPY + IRSEX, statistic="F", design = design.nsduh)
svychisq(~PNRNMLAS1 + IRSEX, statistic="F", design = design.nsduh)
svychisq(~herused + IRSEX, statistic="F", design = design.nsduh)
svychisq(~IRHERFM + IRSEX, statistic="F", design = design.nsduh)
svychisq(~pnrnused + IRSEX, statistic="F", design = design.nsduh)
svychisq(~IRPNRNM30FQ + IRSEX, statistic="F", design = design.nsduh)
svychisq(~mhtused + IRSEX, statistic="F", design = design.nsduh)
svychisq(~MHTINPPY + IRSEX, statistic="F", design = design.nsduh)