Untitled
Dawoon Lee
2025-11-28
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
setwd("~/Desktop/FDU LAB")
TLFB2 <- read.csv("TLFB_9-22-22(TLFB).csv")
TLFB <-TLFB2[-c(1,2, 618),]
# summary(TLFB)
# require(psych)
# describe(TLFB)
# getOption("max.print")
# options(max.print = 100000)
# describeTLFB <- describe(TLFB)
# describeTLFB
# TLFB$Condom_Access
WIDE <- read.csv("WIDE_FINAL(WIDE_FINAL).csv")
#describe(WIDE)
#describeWIDE <- describe(WIDE)
#describeWIDE
#There are two datasets from my dissertation project, one is in wide format and
#the other is in long format. The datasets are the same,
#the data structure differs based on the type of analyses that were conducted -
#long for multilevel models and wide for everything else.
###### variables
#####demographics
#demo <- WIDE[,1:23]
######## condomaccess
#CDAC <-TLFB$Condom_Access
#######partner STI
#PSTI <-TLFB$STD
require(stringr)## Loading required package: stringrTLFB$STD <- gsub("2", "1", TLFB$STD)
#TLFB$STD
######partner types
#ptype <- TLFB$Partner_Type
#TLFB$Partner_Type
#########Motivation
#moti <-TLFB$UCLA_Total
#########information
#info <-TLFB$STD_KQ_Total
#########behavior skill
# beha <- TLFB$CISQ_W_Total
# beha <- TLFB$CISQ_S_Total
# beha <- TLFB$CISQ_DR_Total
# beha <- TLFB$CISQ_RC_Total
# beha <- TLFB$CISQ_R_Total
# beha <- TLFB$CISQ_D_Total
# beha <- TLFB$CISQ_P_Total
BSsum <- data.frame(TLFB$CISQ_W_Total, TLFB$CISQ_S_Total, TLFB$CISQ_DR_Total, TLFB$CISQ_RC_Total,
TLFB$CISQ_R_Total, TLFB$CISQ_D_Total, TLFB$CISQ_P_Total)
TLFB$BSsum <-rowSums(BSsum)
#TLFB$BSsum
##########substance use
#drinkuse <- TLFB$Drink_YN
#drug <-TLFB$Drugs_YN
######??lme4
#library(dplyr)
#install.packages("pacman")
#pacman::p_load(Matrix)
library(lme4)## Loading required package: Matrix###### as numeric
TLFB$UCLA_Total <- as.numeric(as.character(TLFB$UCLA_Total))
TLFB$STD_KQ_Total <- as.numeric(as.character(TLFB$STD_KQ_Total))
TLFB$BSsum <- as.numeric(as.character(TLFB$BSsum))
TLFB$Condom_Access <- as.numeric(as.character(TLFB$Condom_Access))
TLFB$ID <- as.numeric(as.character(TLFB$ID))
TLFB$STD <- as.numeric(as.character(TLFB$STD))
TLFB$Gender <- as.numeric(as.character(TLFB$Gender))
#TLFB$Drink_YN
#TLFB$Drugs_YN
# combine the drug and alcohol use
subs <- data.frame(as.numeric(TLFB$Drink_YN), as.numeric(TLFB$Drugs_YN))
require(stringr)
#subs$as.numeric.TLFB.Drugs_YN. <- as.numeric(gsub("1", "2", subs$as.numeric.TLFB.Drugs_YN.)) #drug = 2
#TLFB$substanceuse <-rowSums(subs) # 0 = no use, 1 = drink, 2 = drug, 3 = combined
#TLFB$substanceuse
TLFB$substanceuse <-rowSums(subs)
TLFB$substanceuse <- gsub("2", "1", TLFB$substanceuse)
#TLFB$substanceuse
#TLFB$Partner_Type
#TLFB$Partner_Type <- gsub("1", "0", TLFB$Partner_Type)
#TLFB$Partner_Type <- gsub("2", "1", TLFB$Partner_Type) # 0 New + casual, 1 regular
#TLFB$Partner_Type
###scaled numeric values
TLFB$UCLA_Total_scaled <- scale(TLFB$UCLA_Total)
TLFB$STD_KQ_Total_scaled <- scale(TLFB$STD_KQ_Total)
TLFB$BSsum_scaled <- scale(TLFB$BSsum)
TLFB$ID_scaled <- scale(TLFB$ID)
#Test whether the categorical variables are factor variables
#Change to factor if the results of the above were "FALSE"
TLFB$substanceuse <- factor(TLFB$substanceuse, levels = c(0,1))
TLFB$Gender <- factor(TLFB$Gender, levels = c(0,1))
TLFB$Condom_Access <- factor(TLFB$Condom_Access, levels = c(0,1))
TLFB$STD <- factor(TLFB$STD, levels = c(0,1))
TLFB$Partner_Type <- factor(TLFB$Partner_Type, levels = c(0,1,2))
TLFB$Drugs_YN <- factor(TLFB$Drugs_YN, levels = c(0,1))
TLFB$Drink_YN <- factor(TLFB$Drink_YN, levels = c(0,1))
is.factor(TLFB$substanceuse)## [1] TRUEis.factor(TLFB$Gender)## [1] TRUEis.factor(TLFB$Condom_Access)## [1] TRUEis.factor(TLFB$STD)## [1] TRUEis.factor(TLFB$Partner_Type)## [1] TRUEis.factor(TLFB$Drink_YN)## [1] TRUEis.factor(TLFB$Drugs_YN)## [1] TRUE###Fit the full model
model_condom_accessibility <- glmer(Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender + substanceuse + Partner_Type +
(substanceuse + Partner_Type | ID), # Complex random effects
data = TLFB, family = binomial(link = "logit"),
glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000))) # Set a higher maxfun## boundary (singular) fit: see help('isSingular')summary(model_condom_accessibility)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled +
## Gender + substanceuse + Partner_Type + (substanceuse + Partner_Type |
## ID)
## Data: TLFB
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000))
##
## AIC BIC logLik -2*log(L) df.resid
## 337.7 417.3 -150.8 301.7 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1845 0.0069 0.0292 0.1047 2.0561
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 91.96 9.590
## substanceuse1 87.21 9.338 -0.89
## Partner_Type1 73.65 8.582 -0.59 0.27
## Partner_Type2 16.49 4.060 -0.24 -0.21 0.83
## Number of obs: 615, groups: ID, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.9107 2.7692 3.579 0.000345 ***
## UCLA_Total_scaled 1.1220 0.4670 2.402 0.016286 *
## STD_KQ_Total_scaled 0.1875 0.4275 0.439 0.660881
## BSsum_scaled 0.7729 0.5651 1.368 0.171380
## Gender1 -2.0545 1.2306 -1.670 0.095011 .
## substanceuse1 -4.2369 2.0399 -2.077 0.037802 *
## Partner_Type1 2.5878 2.2555 1.147 0.251248
## Partner_Type2 -0.9259 1.3669 -0.677 0.498188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) UCLA_T STD_KQ BSsm_s Gendr1 sbstn1 Prt_T1
## UCLA_Ttl_sc 0.033
## STD_KQ_Ttl_ 0.064 0.041
## BSsum_scald 0.038 0.457 -0.254
## Gender1 -0.414 0.344 -0.134 0.403
## substances1 -0.777 -0.093 -0.022 -0.098 0.044
## Partnr_Typ1 -0.298 0.100 0.090 0.084 0.025 0.029
## Partnr_Typ2 -0.541 0.043 -0.054 0.017 0.174 0.048 0.544
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')isSingular(model_condom_accessibility)## [1] TRUE####Overall Fit & Uniformity: Generate simulated residuals and plot them
require(DHARMa)## Loading required package: DHARMa## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')simulationOutput <- DHARMa::simulateResiduals(fittedModel = model_condom_accessibility, plot = TRUE)
####The points in the Q-Q plot should align with the diagonal line,
####and the p-value for the Kolmogorov-Smirnov test should be above 0.05, indicating uniformity.
DHARMa::testDispersion(simulationOutput)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9949, p-value = 1
## alternative hypothesis: two.sided###A p-value below 0.05 usually indicates a problem that needs addressing (e.g., using a negative binomial family instead of Poisson).
qqnorm(ranef(model_condom_accessibility)[[1]][,1])
qqline(ranef(model_condom_accessibility)[[1]][,1])
###Normality of Random Effects: Visually inspect the distribution of the extracted random effects using Q-Q plots to ensure they are approximately normal:
###The Q-Q plot compares the distribution of your sample data (the extracted random effects) against a theoretical normal distribution.
###If the data are normally distributed, the points should fall closely along a straight diagonal line.
##The visual inspection of the provided Q-Q plot suggests a deviation from the normal distribution, particularly in the lower (left) tail.
summary(model_condom_accessibility)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled +
## Gender + substanceuse + Partner_Type + (substanceuse + Partner_Type |
## ID)
## Data: TLFB
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000))
##
## AIC BIC logLik -2*log(L) df.resid
## 337.7 417.3 -150.8 301.7 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.1845 0.0069 0.0292 0.1047 2.0561
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 91.96 9.590
## substanceuse1 87.21 9.338 -0.89
## Partner_Type1 73.65 8.582 -0.59 0.27
## Partner_Type2 16.49 4.060 -0.24 -0.21 0.83
## Number of obs: 615, groups: ID, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.9107 2.7692 3.579 0.000345 ***
## UCLA_Total_scaled 1.1220 0.4670 2.402 0.016286 *
## STD_KQ_Total_scaled 0.1875 0.4275 0.439 0.660881
## BSsum_scaled 0.7729 0.5651 1.368 0.171380
## Gender1 -2.0545 1.2306 -1.670 0.095011 .
## substanceuse1 -4.2369 2.0399 -2.077 0.037802 *
## Partner_Type1 2.5878 2.2555 1.147 0.251248
## Partner_Type2 -0.9259 1.3669 -0.677 0.498188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) UCLA_T STD_KQ BSsm_s Gendr1 sbstn1 Prt_T1
## UCLA_Ttl_sc 0.033
## STD_KQ_Ttl_ 0.064 0.041
## BSsum_scald 0.038 0.457 -0.254
## Gender1 -0.414 0.344 -0.134 0.403
## substances1 -0.777 -0.093 -0.022 -0.098 0.044
## Partnr_Typ1 -0.298 0.100 0.090 0.084 0.025 0.029
## Partnr_Typ2 -0.541 0.043 -0.054 0.017 0.174 0.048 0.544
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')# install.packages("sjPlot")
library(sjPlot)
# Plot the fixed effects (estimates transformed from log-odds to odds ratios by default)
plot_model(model_condom_accessibility, type = "est")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the sjPlot package.
## Please report the issue at <https://github.com/strengejacke/sjPlot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
library(ggeffects)
library(ggplot2) # Optional, for further customization##
## Attaching package: 'ggplot2'## The following object is masked from 'package:sjPlot':
##
## set_theme# Use sjPlot to quickly generate Q-Q plots for all random effects
#plot_model(model_condom_accessibility, type = "pred")
#first_re_plot <- plot_model(model_condom_accessibility, type = "re")[[1]]
#se_re_plot <- plot_model(model_condom_accessibility, type = "re")[[2]]
#th_re_plot <- plot_model(model_condom_accessibility, type = "re")[[3]]
#modified_plot <- first_re_plot +
# labs(
# title = "Participant ID Random Effects",
# x = "ID",
# y = "Condom Access",
# caption = "Model: Condom_Access"
# )
#
# print(modified_plot)
#
# modified_plot2 <- se_re_plot +
# labs(
# title = "Partner Type Random Effects",
# subtitle = "0 = New, 1 = Causal, 2 = Regular",
# x = "Partner Type",
# y = "Condom Access",
# caption = "Model: Condom_Access"
# )
#
# print(modified_plot2)
#
#
#
# modified_plot3 <- th_re_plot +
# labs(
# title = "Substance UseRandom Effects",
# subtitle = "0 = no use, 1 = use",
# x = "Substance use",
# y = "Condom Access",
# caption = "Model: Condom_Access"
# )
#
# print(modified_plot3)
sjPlot::tab_model(model_condom_accessibility)## boundary (singular) fit: see help('isSingular')| Condom_Access | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 20144.51 | 88.51 – 4584926.51 | <0.001 |
| UCLA Total scaled | 3.07 | 1.23 – 7.67 | 0.016 |
| STD KQ Total scaled | 1.21 | 0.52 – 2.79 | 0.661 |
| BSsum scaled | 2.17 | 0.72 – 6.56 | 0.171 |
| Gender [1] | 0.13 | 0.01 – 1.43 | 0.095 |
| substanceuse [1] | 0.01 | 0.00 – 0.79 | 0.038 |
| Partner Type [1] | 13.30 | 0.16 – 1106.00 | 0.251 |
| Partner Type [2] | 0.40 | 0.03 – 5.77 | 0.498 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 91.96 | ||
| τ11 ID.substanceuse1 | 87.21 | ||
| τ11 ID.Partner_Type1 | 73.65 | ||
| τ11 ID.Partner_Type2 | 16.49 | ||
| ρ01 | -0.89 | ||
| -0.59 | |||
| -0.24 | |||
| N ID | 53 | ||
| Observations | 615 | ||
| Marginal R2 / Conditional R2 | 0.691 / NA | ||
#library(performance)
# This will return the correct ICC calculation for your model
#icc(model_condom_accessibility)
##Warning message:
#Can't compute random effect variances. Some variance components equal zero. Your model may suffer from singularity (see
# `?lme4::isSingular` and `?performance::check_singularity`).
# Decrease the `tolerance` level to force the calculation of random effect variances, or impose priors on your random effects
# parameters (using packages like `brms` or `glmmTMB`).
## Random effects - P VALUE
# Make sure to use the 'model_condom_accessibility' you defined earlier.
# The original model already used maximum likelihood (Laplace Approximation)
# 1. Test the significance of the 'ID' random effect
# Create a reduced model that removes (1 | ID) but keeps other random effects
# model_no_ID <- glmer(Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender +
# (substanceuse | ID) + (Partner_Type | ID),
# data = TLFB, family = binomial(link = "logit"),
# glmerControl(optimizer = "bobyqa"))
#
# Compare the full model to the model without ID
#anova(model_no_ID, model_condom_accessibility)
# The output will provide a Chi-squared statistic and a p-value for the difference in fit.
# 2. Test the significance of the 'Partner_Type' random effect
# Create a reduced model that removes (1 | Partner_Type)
# model_no_PT <- glmer(Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender +
# (1 | ID) + (substanceuse | ID),
# data = TLFB, family = binomial(link = "logit"),
# glmerControl(optimizer = "bobyqa"))
# Compare the full model to the model without Partner_Type
#anova(model_no_PT, model_condom_accessibility)
# 3. Test the significance of the 'substanceuse' random effect
# Create a reduced model that removes (1 | substanceuse)
#model_no_SU <- glmer(Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender +
# (1 | ID) + (Partner_Type | ID),
# data = TLFB, family = binomial(link = "logit"),
# glmerControl(optimizer = "bobyqa"))
# Compare the full model to the model without substanceuse
#anova(model_no_SU, model_condom_accessibility)
#since there is no variance on substance use,
#model_condom_accessibility_NoSub <-glmer(Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender +
# (1 |ID)+
# (1 |Partner_Type),
# data = TLFB, family = binomial(link = "logit"),
# glmerControl(optimizer = "bobyqa"))
#summary(model_condom_accessibility_NoSub)
### to get rid of error
TLFB$Partner_Type## [1] 0 0 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 0 1 0 2 2
## [38] 0 1 1 1 1 1 0 1 1 1 0 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 1
## [75] 1 1 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 1 1 1 1 2 2 2 2 2 2 2 1 1 0 1 1 1 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 2 0 1 0 1 0 0 0 0 1 1 0
## [149] 1 0 1 1 1 0 0 1 0 0 0 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 1 1 1 1 2 2
## [223] 2 0 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 0 1 2 0 0 0
## [260] 1 0 1 1 0 2 2 2 1 1 1 1 1 0 2 2 2 2 2 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2
## [297] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [334] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 0 0 0 0 1 1 1
## [371] 1 1 1 1 1 1 1 2 2 2 2 0 2 0 0 1 1 1 0 0 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [408] 2 0 1 1 1 0 1 1 1 1 1 1 1 1 1 2 2 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 2 2
## [445] 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 1 1 0 0 0
## [482] 1 0 1 1 0 1 0 1 1 2 2 1 2 2 2 1 2 1 2 2 2 2 2 2 2 2 0 0 0 1 0 0 1 1 1 0 1
## [519] 1 0 2 2 2 0 0 2 2 2 2 2 2 2 2 0 1 1 1 0 1 1 1 1 2 2 2 2 2 2 2 2 0 1 0 1 1
## [556] 1 0 1 1 0 0 1 2 2 2 2 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 1 1 1 0 0 0
## [593] 0 1 2 2 2 2 2 2 2 0 2 2 0 0 0 1 0 0 0 1 0 1 1
## Levels: 0 1 2TLFB$Partner_Type_modified <- gsub("1", "0", TLFB$Partner_Type)
TLFB$Partner_Type_modified <- gsub("2", "1", TLFB$Partner_Type) # 0 New + casual, 1 regular
TLFB$Partner_Type_modified## [1] "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [19] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1" "0" "1"
## [37] "1" "0" "1" "1" "1" "1" "1" "0" "1" "1" "1" "0" "1" "1" "1" "1" "1" "1"
## [55] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1" "1"
## [73] "1" "1" "1" "1" "1" "1" "0" "0" "0" "0" "1" "0" "0" "0" "1" "1" "0" "0"
## [91] "0" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1"
## [109] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [127] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1" "0" "1" "0" "0" "0"
## [145] "0" "1" "1" "0" "1" "0" "1" "1" "1" "0" "0" "1" "0" "0" "0" "1" "1" "1"
## [163] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [181] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [199] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0"
## [217] "1" "1" "1" "1" "1" "1" "1" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [235] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [253] "1" "0" "1" "1" "0" "0" "0" "1" "0" "1" "1" "0" "1" "1" "1" "1" "1" "1"
## [271] "1" "1" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [289] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [307] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [325] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [343] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [361] "1" "1" "1" "0" "0" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [379] "1" "1" "1" "0" "1" "0" "0" "1" "1" "1" "0" "0" "0" "1" "1" "1" "1" "1"
## [397] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1" "1" "1" "0" "1"
## [415] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "0" "0" "0" "0" "1" "1" "1"
## [433] "1" "1" "1" "1" "1" "1" "1" "0" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1"
## [451] "1" "1" "1" "1" "1" "1" "1" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [469] "1" "1" "1" "1" "1" "0" "1" "1" "1" "1" "0" "0" "0" "1" "0" "1" "1" "0"
## [487] "1" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [505] "1" "1" "1" "0" "0" "0" "1" "0" "0" "1" "1" "1" "0" "1" "1" "0" "1" "1"
## [523] "1" "0" "0" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1" "1" "1" "0" "1" "1"
## [541] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1" "0" "1" "0" "1" "1" "1" "0" "1"
## [559] "1" "0" "0" "1" "1" "1" "1" "1" "0" "1" "1" "1" "1" "1" "1" "1" "1" "1"
## [577] "1" "1" "1" "1" "1" "0" "1" "1" "1" "0" "1" "1" "1" "0" "0" "0" "0" "1"
## [595] "1" "1" "1" "1" "1" "1" "1" "0" "1" "1" "0" "0" "0" "1" "0" "0" "0" "1"
## [613] "0" "1" "1"TLFB$Partner_Type_modified <- factor(TLFB$Partner_Type_modified, levels = c(0,1))
model_condom_accessibility_reduced <- glmer(Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender + substanceuse + Partner_Type_modified +
(substanceuse + Partner_Type_modified | ID), # correlated random effects
data = TLFB, family = binomial(link = "logit"),
glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000))) # Set a higher maxfun
summary(model_condom_accessibility_reduced)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Condom_Access ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled +
## Gender + substanceuse + Partner_Type_modified + (substanceuse +
## Partner_Type_modified | ID)
## Data: TLFB
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 10000))
##
## AIC BIC logLik -2*log(L) df.resid
## 341.2 398.7 -157.6 315.2 602
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.0957 0.0243 0.0579 0.1431 2.7340
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 46.30 6.804
## substanceuse1 36.91 6.075 -0.91
## Partner_Type_modified1 15.26 3.906 -0.48 0.15
## Number of obs: 615, groups: ID, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.83133 2.88078 2.371 0.0177 *
## UCLA_Total_scaled 0.94277 0.43515 2.167 0.0303 *
## STD_KQ_Total_scaled 0.09541 0.37887 0.252 0.8012
## BSsum_scaled 0.69304 0.51499 1.346 0.1784
## Gender1 -1.19850 1.10360 -1.086 0.2775
## substanceuse1 -2.73873 2.35983 -1.161 0.2458
## Partner_Type_modified1 0.67910 1.17467 0.578 0.5632
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) UCLA_T STD_KQ BSsm_s Gendr1 sbstn1
## UCLA_Ttl_sc 0.044
## STD_KQ_Ttl_ 0.103 -0.056
## BSsum_scald 0.028 0.402 -0.323
## Gender1 -0.386 0.246 -0.212 0.438
## substances1 -0.873 -0.041 -0.050 -0.066 0.087
## Prtnr_Typ_1 -0.445 0.029 -0.096 0.026 0.238 0.126isSingular(model_condom_accessibility_reduced)## [1] FALSE########STD
###Fit the full model
model_STD <-glmer(STD ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender +substanceuse + Partner_Type+
(substanceuse + Partner_Type |ID),
data = TLFB, family = binomial(link = "logit"),
glmerControl(optimizer = "bobyqa"))## boundary (singular) fit: see help('isSingular')summary(model_STD)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: STD ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled +
## Gender + substanceuse + Partner_Type + (substanceuse + Partner_Type |
## ID)
## Data: TLFB
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik -2*log(L) df.resid
## 289.8 369.4 -126.9 253.8 597
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8892 -0.0275 -0.0022 0.0552 6.0797
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 12.72 3.566
## substanceuse1 8.54 2.922 0.31
## Partner_Type1 135.47 11.639 0.89 -0.16
## Partner_Type2 109.73 10.475 1.00 0.38 0.85
## Number of obs: 615, groups: ID, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.5939 2.0744 1.250 0.211143
## UCLA_Total_scaled -0.3665 0.8200 -0.447 0.654940
## STD_KQ_Total_scaled 1.6951 0.8435 2.010 0.044474 *
## BSsum_scaled -0.9877 1.0742 -0.919 0.357862
## Gender1 -5.0833 2.5871 -1.965 0.049434 *
## substanceuse1 1.7461 1.3172 1.326 0.184967
## Partner_Type1 -3.9770 2.5291 -1.573 0.115828
## Partner_Type2 -11.1080 3.0748 -3.613 0.000303 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) UCLA_T STD_KQ BSsm_s Gendr1 sbstn1 Prt_T1
## UCLA_Ttl_sc -0.086
## STD_KQ_Ttl_ 0.299 -0.001
## BSsum_scald 0.088 -0.015 -0.223
## Gender1 -0.792 0.059 -0.420 0.347
## substances1 -0.108 0.071 0.309 -0.016 -0.024
## Partnr_Typ1 0.191 0.174 0.133 0.025 0.031 -0.260
## Partnr_Typ2 -0.079 -0.201 -0.267 0.278 0.185 -0.162 0.146
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')isSingular(model_STD)## [1] TRUE####Overall Fit & Uniformity: Generate simulated residuals and plot them
require(DHARMa)
simulationOutput_STD <- DHARMa::simulateResiduals(fittedModel = model_STD, plot = TRUE)## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
####The points in the Q-Q plot should align with the diagonal line,
####and the p-value for the Kolmogorov-Smirnov test should be above 0.05, indicating uniformity.
DHARMa::testDispersion(simulationOutput_STD)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.0229, p-value = 0.904
## alternative hypothesis: two.sided###A p-value below 0.05 usually indicates a problem that needs addressing (e.g., using a negative binomial family instead of Poisson).
qqnorm(ranef(model_STD)[[1]][,1])
qqline(ranef(model_STD)[[1]][,1])
###Normality of Random Effects: Visually inspect the distribution of the extracted random effects using Q-Q plots to ensure they are approximately normal:
###The Q-Q plot compares the distribution of your sample data (the extracted random effects) against a theoretical normal distribution.
###If the data are normally distributed, the points should fall closely along a straight diagonal line.
##The visual inspection of the provided Q-Q plot suggests a deviation from the normal distribution, particularly in the lower (left) tail.
#summary(model_STD)
# install.packages("sjPlot")
library(sjPlot)
# Plot the fixed effects (estimates transformed from log-odds to odds ratios by default)
plot_model(model_STD, type = "est")
#library(ggeffects)
#library(ggplot2) # Optional, for further customization
# Use sjPlot to quickly generate Q-Q plots for all random effects
# plot_model(model_STD, type = "pred")
# first_re_plot <- plot_model(model_STD, type = "re")[[1]]
# se_re_plot <- plot_model(model_STD, type = "re")[[2]]
# th_re_plot <- plot_model(model_STD, type = "re")[[3]]
#
# modified_plot <- first_re_plot +
# labs(
# title = "Participant ID Random Effects",
# x = "ID",
# y = "STD",
# caption = "Model: STD"
# )
#
# print(modified_plot)
#
# modified_plot2 <- se_re_plot +
# labs(
# title = "Partner Type Random Effects",
# subtitle = "0 = New, 1 = Causal, 2 = Regular",
# x = "Partner Type",
# y = "STD",
# caption = "Model: STD"
# )
#
# print(modified_plot2)
#
#
#
# modified_plot3 <- th_re_plot +
# labs(
# title = "Substance UseRandom Effects",
# subtitle = "0 = no use, 1 = use",
# x = "Substance use",
# y = "STD",
# caption = "Model: STD"
# )
#
# print(modified_plot3)
#
tab_model(model_STD)## boundary (singular) fit: see help('isSingular')| STD | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 13.38 | 0.23 – 780.30 | 0.211 |
| UCLA Total scaled | 0.69 | 0.14 – 3.46 | 0.655 |
| STD KQ Total scaled | 5.45 | 1.04 – 28.46 | 0.044 |
| BSsum scaled | 0.37 | 0.05 – 3.06 | 0.358 |
| Gender [1] | 0.01 | 0.00 – 0.99 | 0.049 |
| substanceuse [1] | 5.73 | 0.43 – 75.77 | 0.185 |
| Partner Type [1] | 0.02 | 0.00 – 2.66 | 0.116 |
| Partner Type [2] | 0.00 | 0.00 – 0.01 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 ID | 12.72 | ||
| τ11 ID.substanceuse1 | 8.54 | ||
| τ11 ID.Partner_Type1 | 135.47 | ||
| τ11 ID.Partner_Type2 | 109.73 | ||
| ρ01 | 0.31 | ||
| 0.89 | |||
| 1.00 | |||
| N ID | 53 | ||
| Observations | 615 | ||
| Marginal R2 / Conditional R2 | 0.891 / NA | ||
########reduced version ##### I need to drop one of random slope.. ## still not fitting.
model_STD_reduced <-glmer(STD ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled + Gender +substanceuse + Partner_Type_modified+
(substanceuse + Partner_Type_modified || ID), ## uncorrelated random slope
data = TLFB, family = binomial(link = "logit"),
glmerControl(optimizer = "bobyqa"))## boundary (singular) fit: see help('isSingular')summary(model_STD_reduced)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: STD ~ UCLA_Total_scaled + STD_KQ_Total_scaled + BSsum_scaled +
## Gender + substanceuse + Partner_Type_modified + (substanceuse +
## Partner_Type_modified || ID)
## Data: TLFB
## Control: glmerControl(optimizer = "bobyqa")
##
## AIC BIC logLik -2*log(L) df.resid
## 338.0 399.9 -155.0 310.0 601
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.7554 -0.1370 -0.0197 0.0857 3.6909
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 8.956e-10 2.993e-05
## ID.1 substanceuse0 2.018e+01 4.492e+00
## substanceuse1 1.012e+01 3.182e+00 1.00
## ID.2 Partner_Type_modified0 6.204e+00 2.491e+00
## Partner_Type_modified1 5.617e+01 7.495e+00 1.00
## Number of obs: 615, groups: ID, 53
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4978 1.9733 0.759 0.44783
## UCLA_Total_scaled -2.1690 1.0439 -2.078 0.03772 *
## STD_KQ_Total_scaled 1.1272 0.7989 1.411 0.15828
## BSsum_scaled -0.1104 1.4680 -0.075 0.94005
## Gender1 -3.1709 2.3639 -1.341 0.17981
## substanceuse1 1.9730 0.7978 2.473 0.01340 *
## Partner_Type_modified1 -4.3270 1.6198 -2.671 0.00755 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) UCLA_T STD_KQ BSsm_s Gendr1 sbstn1
## UCLA_Ttl_sc -0.197
## STD_KQ_Ttl_ 0.181 -0.295
## BSsum_scald 0.090 0.322 -0.339
## Gender1 -0.717 0.285 -0.301 0.353
## substances1 -0.322 -0.021 0.097 0.236 0.168
## Prtnr_Typ_1 0.042 0.120 -0.184 0.163 0.003 -0.013
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')isSingular(model_STD_reduced)## [1] TRUE