Disclosure Analyses
2024-01-30
pacman::p_load(tidyverse, MASS,
magrittr, robustHD,
tibble, psych,
kableExtra, e1071,
knitr, tidyr,
lavaan, semPlot,
jtools, car,
lmtest, ggpubr,
FSA, rstatix,
writexl, readxl,
rcompanion, coin, lm.beta,
rstatix, Hmisc, pwr, nnet,
gtsummary, marginaleffects,
ggeffects, margins, VGAM, brant,
splines, effects, ggeffects, foreign, nnet, reshape2)Data Cleaning and Variable Creation
SAdata <- read_excel("~/Downloads/SAdata.xlsx")
SAdata$SV_1_disclosure<- as.numeric(SAdata$SV_1_disclosure) #making these variable numeric
SAdata$SV_2_disclosure<- as.numeric(SAdata$SV_2_disclosure)
SAdata$SV_3_disclosure<- as.numeric(SAdata$SV_3_disclosure)
SAdata$SV_4_disclosure<- as.numeric(SAdata$SV_4_disclosure)
SA.data.1<- subset(SAdata, is.na(SV_perp_1)) #getting rid of people who were assaulted by their current partners
disclosedata<- SA.data.1 %>%
select(SubjectID, relLength, Age_1, Relationship_type, Education ,Income, Race_1, Race_2, Race_3, Race_4, Race_5, Race_6, Race_7, SV_1_disclosure, SV_2_disclosure, SV_3_disclosure, SV_4_disclosure, PRQC, NSSS, PRI_total,Trust, AAQ.ANX, AAQ.AV)
disclosedata <- disclosedata %>% #two people said not assaulted by current partners but also listed current partner as perp for disclosure question
dplyr::filter(SubjectID != 90,
SubjectID !=325)
disclosedata <- disclosedata %>% #recoding so higher scores mean you disclosed more
mutate(
SV1.disclose = dplyr::recode(`SV_1_disclosure`, "1"=3, "2"=2, "3"=1))
disclosedata <- disclosedata %>%
mutate(
SV2.disclose = dplyr::recode(`SV_2_disclosure`, "1"=3, "2"=2, "3"=1))
disclosedata <- disclosedata %>%
mutate(
SV3.disclose = dplyr::recode(`SV_3_disclosure`, "1"=3, "2"=2, "3"=1))
disclosedata <- disclosedata %>%
mutate(
SV4.disclose = dplyr::recode(`SV_4_disclosure`, "1"=3, "2"=2, "3"=1))
disclosedata <- disclosedata %>% #creating disclosure ratio variable for amount of experiences you disclosed to current partner
rowwise() %>%
mutate(disclosure.ratio = mean(c(SV1.disclose, SV2.disclose, SV3.disclose, SV4.disclose), na.rm = TRUE))%>%
ungroup()
hist(disclosedata$disclosure.ratio)
mean(disclosedata$disclosure.ratio)## [1] 1.938131sd(disclosedata$disclosure.ratio)## [1] 0.8121303sum(disclosedata$disclosure.ratio>1 & disclosedata$disclosure.ratio<3)/length(disclosedata$disclosure.ratio) #proportion that disclosed some of experiences## [1] 0.3787879sum(disclosedata$disclosure.ratio==1)/length(disclosedata$disclosure.ratio) #proportion that disclosed none of experience## [1] 0.3363636sum(disclosedata$disclosure.ratio==3)/length(disclosedata$disclosure.ratio) #prop that disclosed all ## [1] 0.2848485# Merging sexual orientation and gender recoded data
OrientationGender <- read_excel("~/Downloads/RECODED Sexual Orientation_Gender Complete.xlsx")
OrientationGender <- OrientationGender %>% #renaming ID variable to merge two dataframes
rename(SubjectID = ID)
disclosedata<-merge(x=disclosedata, y=OrientationGender, by='SubjectID', all.x=T) #merging data frames
# Creating a variable that reflects whether participants are LGBTQ (so either not heterosexual or #not cisgender) OR straight and cisgender
disclosedata<- disclosedata %>%
mutate(LGBTQ = case_when(
disclosedata$Questioning==1 ~ "LGBTQ",
disclosedata$Orientation=="Sexual or heterosexual" ~ "StraightCis",
disclosedata$Orientation=="Bisexual or pansexual" ~ "LGBTQ",
disclosedata$Orientation=="Lesbian or gay" ~ "LGBTQ",
disclosedata$Orientation=="other" ~ "LGBTQ",
disclosedata$GenderQueer==1 ~ "LGBTQ",
disclosedata$TransMan==1 ~ "LGBTQ",
disclosedata$TransWoman==1 ~ "LGBTQ",
disclosedata$Orientation=="Sexual or heterosexual" & disclosedata$Questioning==1 ~ "LGBTQ"
))###### Disclosure frequencies, 4 indicates people who said that current partner was perp (these people were removed for correlations)
kableExtra::kable(table(SAdata$SV_1_disclosure), booktabs = TRUE, col.names = c("SV 1 Disclosure", "Frequency")) %>%
kableExtra::kable_styling(bootstrap_options = "striped", full_width = TRUE) | SV 1 Disclosure | Frequency |
|---|---|
| 1 | 71 |
| 2 | 88 |
| 3 | 62 |
| 4 | 5 |
kableExtra::kable(table(SAdata$SV_2_disclosure), booktabs = TRUE, col.names = c("SV 2 Disclosure", "Frequency")) %>%
kableExtra::kable_styling(bootstrap_options = "striped", full_width = TRUE) | SV 2 Disclosure | Frequency |
|---|---|
| 1 | 72 |
| 2 | 62 |
| 3 | 72 |
| 4 | 2 |
kableExtra::kable(table(SAdata$SV_3_disclosure), booktabs = TRUE, col.names = c("SV 3 Disclosure", "Frequency")) %>%
kableExtra::kable_styling(bootstrap_options = "striped", full_width = TRUE)| SV 3 Disclosure | Frequency |
|---|---|
| 1 | 56 |
| 2 | 55 |
| 3 | 75 |
| 4 | 6 |
kableExtra::kable(table(SAdata$SV_4_disclosure), booktabs = TRUE, col.names = c("SV 4 Disclosure", "Frequency")) %>%
kableExtra::kable_styling(bootstrap_options = "striped", full_width = TRUE) | SV 4 Disclosure | Frequency |
|---|---|
| 1 | 45 |
| 2 | 38 |
| 3 | 66 |
| 4 | 5 |
disclosedata<- disclosedata %>% #Creating variable for that represents whether participants are people of color
mutate(POC = case_when(
disclosedata$Race_1==1 ~ "POC", #Alaskan Native, Native American, or Native Hawaiian/Pacific Islander
disclosedata$Race_2==1 ~ "POC", #Asian American or Asian
disclosedata$Race_3==1 ~ "POC", #Black or African American
disclosedata$Race_4==1 ~ "POC", #Latinx or Hispanic
disclosedata$Race_5==1 ~ "POC", #Middle Eastern or North African
disclosedata$Race_6==1 ~ "POC" #Native American, Hawaiian Native, or Alaskan Native
)) #NA variables mean that they weren't Asian, Black, etc.
disclosedata$POC<- disclosedata$POC %>% replace_na("White") # Replacing NA variables with "White"
disclosedata<- disclosedata %>% #Created categorical disclosure variable
mutate(disclose.cat = case_when(
disclosedata$disclosure.ratio==1 ~ 1, #None
disclosedata$disclosure.ratio>1 & disclosedata$disclosure.ratio<3 ~ 2, #Some
disclosedata$disclosure.ratio==3 ~ 3 #All
))
disclosedata$disclose.cat.fac<- as.factor(disclosedata$disclose.cat)
disclosedata$POC <- as.factor(disclosedata$POC)
disclosedata$LGBTQ <- as.factor(disclosedata$LGBTQ)Scale Reliabilities
disclose.reliabilities<- disclosedata %>%
select(SubjectID)
items<- SAdata %>%
select(SubjectID, PRQC_1:PRQC_18, AAQ_1.r, AAQ_2, AAQ_3.r, AAQ_5:AAQ_9, AAQ_4.r, AAQ_10, AAQ_11, AAQ_12.r, AAQ_13, AAQ_14.r, AAQ_15, AAQ_16.r, AAQ_17.r, PRI_1:PRI_4,PRI_5.r: PRI_8.r)
disclose.reliabilities<-merge(x=disclose.reliabilities, y=items, by='SubjectID', all.x=T) #merging data frames
#PRQC
disclose.reliabilities %>%
select(PRQC_1:PRQC_18) %>%
psych::alpha( ,check.keys = F, na.rm = TRUE)##
## Reliability analysis
## Call: psych::alpha(x = ., na.rm = TRUE, check.keys = F)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.96 0.97 0.99 0.62 30 0.0031 5.7 1.1 0.63
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.96 0.96 0.97
## Duhachek 0.96 0.96 0.97
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## PRQC_1 0.96 0.96 0.98 0.61 27 0.0034 0.019 0.63
## PRQC_2 0.96 0.96 0.98 0.61 27 0.0033 0.019 0.63
## PRQC_3 0.96 0.96 0.98 0.61 27 0.0034 0.019 0.63
## PRQC_4 0.96 0.97 0.98 0.63 28 0.0032 0.019 0.64
## PRQC_5 0.96 0.97 0.98 0.62 28 0.0032 0.019 0.63
## PRQC_6 0.96 0.97 0.98 0.62 28 0.0032 0.019 0.64
## PRQC_7 0.96 0.97 0.98 0.62 28 0.0033 0.020 0.64
## PRQC_8 0.96 0.96 0.98 0.62 27 0.0033 0.020 0.63
## PRQC_9 0.96 0.96 0.98 0.61 27 0.0034 0.020 0.63
## PRQC_10 0.96 0.97 0.98 0.62 28 0.0032 0.019 0.64
## PRQC_11 0.96 0.97 0.98 0.62 28 0.0033 0.019 0.64
## PRQC_12 0.96 0.97 0.98 0.62 28 0.0032 0.019 0.64
## PRQC_13 0.96 0.97 0.98 0.62 28 0.0033 0.019 0.64
## PRQC_14 0.96 0.97 0.98 0.64 30 0.0030 0.016 0.65
## PRQC_15 0.96 0.97 0.98 0.64 30 0.0029 0.015 0.65
## PRQC_16 0.96 0.97 0.98 0.62 28 0.0032 0.019 0.64
## PRQC_17 0.96 0.96 0.98 0.62 27 0.0033 0.020 0.63
## PRQC_18 0.96 0.96 0.98 0.62 27 0.0033 0.019 0.63
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PRQC_1 330 0.87 0.87 0.87 0.85 5.6 1.33
## PRQC_2 330 0.87 0.87 0.87 0.85 5.7 1.32
## PRQC_3 330 0.88 0.89 0.89 0.86 5.7 1.35
## PRQC_4 330 0.73 0.76 0.75 0.70 6.4 0.97
## PRQC_5 330 0.78 0.81 0.81 0.75 6.4 1.02
## PRQC_6 330 0.78 0.81 0.80 0.75 6.4 1.04
## PRQC_7 330 0.80 0.77 0.76 0.77 5.2 1.57
## PRQC_8 330 0.85 0.84 0.84 0.82 5.9 1.35
## PRQC_9 330 0.88 0.88 0.87 0.86 5.8 1.40
## PRQC_10 330 0.77 0.77 0.76 0.73 5.8 1.58
## PRQC_11 330 0.81 0.81 0.81 0.78 5.9 1.42
## PRQC_12 330 0.76 0.77 0.77 0.73 5.9 1.42
## PRQC_13 330 0.82 0.78 0.78 0.79 4.9 1.76
## PRQC_14 330 0.71 0.66 0.66 0.66 4.4 1.93
## PRQC_15 330 0.70 0.65 0.64 0.64 4.2 1.98
## PRQC_16 330 0.76 0.78 0.77 0.74 6.5 0.97
## PRQC_17 330 0.83 0.84 0.84 0.80 6.0 1.41
## PRQC_18 330 0.82 0.84 0.84 0.79 6.2 1.30
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 7 miss
## PRQC_1 0.01 0.02 0.05 0.09 0.20 0.34 0.28 0
## PRQC_2 0.01 0.02 0.05 0.08 0.17 0.33 0.34 0
## PRQC_3 0.01 0.02 0.05 0.10 0.18 0.32 0.32 0
## PRQC_4 0.00 0.00 0.02 0.05 0.08 0.17 0.68 0
## PRQC_5 0.00 0.00 0.02 0.05 0.08 0.19 0.65 0
## PRQC_6 0.00 0.00 0.02 0.06 0.08 0.19 0.65 0
## PRQC_7 0.02 0.05 0.09 0.15 0.21 0.22 0.26 0
## PRQC_8 0.01 0.02 0.05 0.07 0.17 0.24 0.44 0
## PRQC_9 0.01 0.02 0.05 0.10 0.14 0.26 0.42 0
## PRQC_10 0.03 0.04 0.02 0.08 0.14 0.21 0.48 0
## PRQC_11 0.01 0.03 0.04 0.07 0.15 0.20 0.50 0
## PRQC_12 0.01 0.03 0.04 0.06 0.14 0.24 0.47 0
## PRQC_13 0.05 0.07 0.11 0.11 0.23 0.22 0.20 0
## PRQC_14 0.10 0.11 0.13 0.16 0.17 0.16 0.18 0
## PRQC_15 0.12 0.12 0.13 0.13 0.18 0.15 0.17 0
## PRQC_16 0.00 0.00 0.02 0.03 0.09 0.16 0.69 0
## PRQC_17 0.02 0.02 0.04 0.08 0.11 0.19 0.55 0
## PRQC_18 0.02 0.01 0.03 0.04 0.11 0.17 0.62 0# Trust
disclose.reliabilities %>%
select(PRQC_10:PRQC_12) %>%
psych::alpha( ,check.keys = F, na.rm = TRUE)##
## Reliability analysis
## Call: psych::alpha(x = ., na.rm = TRUE, check.keys = F)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.92 0.93 0.91 0.81 13 0.0076 5.9 1.4 0.78
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.91 0.92 0.94
## Duhachek 0.91 0.92 0.94
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## PRQC_10 0.95 0.95 0.91 0.91 19.9 0.0053 NA 0.91
## PRQC_11 0.85 0.85 0.74 0.74 5.6 0.0167 NA 0.74
## PRQC_12 0.87 0.88 0.78 0.78 7.1 0.0137 NA 0.78
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PRQC_10 330 0.91 0.90 0.80 0.78 5.8 1.6
## PRQC_11 330 0.96 0.96 0.95 0.90 5.9 1.4
## PRQC_12 330 0.94 0.94 0.93 0.87 5.9 1.4
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 7 miss
## PRQC_10 0.03 0.04 0.02 0.08 0.14 0.21 0.48 0
## PRQC_11 0.01 0.03 0.04 0.07 0.15 0.20 0.50 0
## PRQC_12 0.01 0.03 0.04 0.06 0.14 0.24 0.47 0# Attachment Avoidance
disclose.reliabilities %>%
select(AAQ_1.r, AAQ_2, AAQ_3.r, AAQ_5:AAQ_9) %>%
psych::alpha( ,check.keys = F, na.rm = TRUE)##
## Reliability analysis
## Call: psych::alpha(x = ., na.rm = TRUE, check.keys = F)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.88 0.88 0.89 0.47 7.2 0.0099 3.9 1.3 0.49
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.86 0.88 0.9
## Duhachek 0.86 0.88 0.9
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## AAQ_1.r 0.87 0.87 0.87 0.48 6.5 0.011 0.036 0.50
## AAQ_2 0.88 0.87 0.88 0.50 7.0 0.010 0.030 0.51
## AAQ_3.r 0.89 0.89 0.89 0.54 8.3 0.009 0.017 0.53
## AAQ_5 0.85 0.85 0.86 0.45 5.8 0.012 0.026 0.49
## AAQ_6 0.85 0.84 0.85 0.44 5.4 0.013 0.023 0.46
## AAQ_7 0.86 0.85 0.86 0.46 5.9 0.012 0.031 0.46
## AAQ_8 0.85 0.84 0.85 0.44 5.4 0.013 0.023 0.46
## AAQ_9 0.87 0.87 0.88 0.49 6.6 0.011 0.032 0.50
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## AAQ_1.r 330 0.71 0.71 0.65 0.61 4.0 1.8
## AAQ_2 330 0.64 0.65 0.57 0.53 5.0 1.8
## AAQ_3.r 330 0.49 0.49 0.38 0.35 3.4 1.8
## AAQ_5 330 0.81 0.81 0.80 0.74 3.5 1.8
## AAQ_6 330 0.87 0.87 0.88 0.82 3.6 1.8
## AAQ_7 330 0.80 0.80 0.77 0.72 4.7 1.9
## AAQ_8 330 0.87 0.86 0.87 0.81 3.8 1.9
## AAQ_9 330 0.69 0.69 0.62 0.58 3.5 1.8
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 7 miss
## AAQ_1.r 0.11 0.15 0.16 0.16 0.16 0.16 0.09 0
## AAQ_2 0.06 0.06 0.09 0.10 0.22 0.25 0.22 0
## AAQ_3.r 0.15 0.19 0.22 0.19 0.09 0.09 0.07 0
## AAQ_5 0.15 0.18 0.16 0.19 0.16 0.08 0.06 0
## AAQ_6 0.15 0.19 0.12 0.18 0.21 0.08 0.07 0
## AAQ_7 0.09 0.08 0.09 0.12 0.20 0.23 0.18 0
## AAQ_8 0.16 0.15 0.12 0.14 0.20 0.13 0.09 0
## AAQ_9 0.18 0.17 0.16 0.19 0.15 0.08 0.07 0# Attachment Anxiety
disclose.reliabilities %>%
select(AAQ_4.r, AAQ_10, AAQ_11, AAQ_12.r, AAQ_13, AAQ_14.r, AAQ_15, AAQ_16.r, AAQ_17.r) %>%
psych::alpha( ,check.keys = F, na.rm = TRUE)##
## Reliability analysis
## Call: psych::alpha(x = ., na.rm = TRUE, check.keys = F)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.82 0.81 0.85 0.33 4.4 0.015 3.5 1.2 0.35
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.79 0.82 0.85
## Duhachek 0.79 0.82 0.85
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## AAQ_4.r 0.80 0.80 0.83 0.33 3.9 0.016 0.036 0.35
## AAQ_10 0.82 0.82 0.84 0.36 4.4 0.015 0.032 0.38
## AAQ_11 0.78 0.77 0.81 0.30 3.4 0.019 0.035 0.23
## AAQ_12.r 0.80 0.80 0.84 0.33 4.0 0.016 0.035 0.35
## AAQ_13 0.82 0.81 0.83 0.35 4.3 0.015 0.032 0.38
## AAQ_14.r 0.80 0.79 0.83 0.33 3.9 0.017 0.031 0.35
## AAQ_15 0.81 0.80 0.83 0.34 4.1 0.015 0.035 0.38
## AAQ_16.r 0.78 0.78 0.81 0.31 3.5 0.018 0.031 0.35
## AAQ_17.r 0.79 0.78 0.81 0.31 3.6 0.018 0.032 0.35
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## AAQ_4.r 330 0.64 0.62 0.56 0.51 4.3 2.0
## AAQ_10 330 0.45 0.49 0.40 0.33 2.8 1.6
## AAQ_11 330 0.79 0.78 0.77 0.70 3.1 2.1
## AAQ_12.r 330 0.64 0.62 0.54 0.50 3.7 2.2
## AAQ_13 330 0.49 0.53 0.46 0.36 2.5 1.6
## AAQ_14.r 330 0.66 0.64 0.59 0.54 4.6 1.9
## AAQ_15 330 0.54 0.58 0.52 0.41 3.2 1.8
## AAQ_16.r 330 0.76 0.73 0.72 0.66 4.3 2.0
## AAQ_17.r 330 0.73 0.72 0.69 0.63 2.8 1.8
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 7 miss
## AAQ_4.r 0.11 0.12 0.12 0.12 0.18 0.20 0.15 0
## AAQ_10 0.24 0.28 0.21 0.13 0.08 0.04 0.03 0
## AAQ_11 0.32 0.18 0.11 0.12 0.08 0.09 0.10 0
## AAQ_12.r 0.20 0.19 0.14 0.09 0.11 0.10 0.17 0
## AAQ_13 0.40 0.22 0.12 0.14 0.07 0.04 0.02 0
## AAQ_14.r 0.06 0.13 0.11 0.17 0.13 0.18 0.23 0
## AAQ_15 0.19 0.22 0.20 0.12 0.15 0.07 0.05 0
## AAQ_16.r 0.12 0.14 0.11 0.12 0.14 0.18 0.19 0
## AAQ_17.r 0.34 0.18 0.15 0.12 0.09 0.05 0.05 0#PRI Total
disclose.reliabilities %>%
select(PRI_1:PRI_4,PRI_5.r: PRI_8.r) %>%
psych::alpha( ,check.keys = F, na.rm = TRUE)##
## Reliability analysis
## Call: psych::alpha(x = ., na.rm = TRUE, check.keys = F)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.96 0.96 0.97 0.74 22 0.0037 3.9 1.2 0.71
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.95 0.96 0.96
## Duhachek 0.95 0.96 0.96
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## PRI_1 0.95 0.95 0.96 0.74 20 0.0041 0.0066 0.71
## PRI_2 0.95 0.95 0.96 0.74 20 0.0041 0.0066 0.71
## PRI_3 0.95 0.95 0.96 0.73 19 0.0043 0.0080 0.69
## PRI_4 0.95 0.95 0.96 0.74 20 0.0042 0.0073 0.71
## PRI_5.r 0.95 0.95 0.96 0.74 20 0.0042 0.0062 0.71
## PRI_6.r 0.95 0.95 0.96 0.74 19 0.0044 0.0071 0.71
## PRI_7.r 0.95 0.95 0.96 0.73 19 0.0044 0.0072 0.71
## PRI_8.r 0.95 0.95 0.96 0.74 19 0.0044 0.0077 0.71
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PRI_1 330 0.85 0.86 0.85 0.81 3.7 1.2
## PRI_2 330 0.86 0.87 0.85 0.82 3.8 1.3
## PRI_3 330 0.89 0.90 0.88 0.86 3.9 1.3
## PRI_4 330 0.87 0.88 0.87 0.83 3.8 1.3
## PRI_5.r 330 0.87 0.86 0.84 0.82 4.0 1.5
## PRI_6.r 330 0.89 0.88 0.87 0.85 4.0 1.4
## PRI_7.r 330 0.90 0.89 0.88 0.86 3.8 1.5
## PRI_8.r 330 0.89 0.89 0.87 0.85 3.8 1.5
##
## Non missing response frequency for each item
## 0 1 2 3 4 5 miss
## PRI_1 0.01 0.06 0.09 0.20 0.34 0.30 0
## PRI_2 0.02 0.06 0.08 0.16 0.30 0.38 0
## PRI_3 0.02 0.04 0.07 0.15 0.27 0.45 0
## PRI_4 0.02 0.05 0.08 0.17 0.26 0.42 0
## PRI_5.r 0.04 0.07 0.06 0.07 0.20 0.56 0
## PRI_6.r 0.03 0.05 0.08 0.10 0.18 0.56 0
## PRI_7.r 0.04 0.08 0.08 0.12 0.20 0.48 0
## PRI_8.r 0.04 0.08 0.09 0.10 0.19 0.50 0Demographics
# relationship length, age, and study variable descriptives
disclosedata %>%
dplyr::select(disclose.cat, relLength, Age_1, PRQC, Trust, PRI_total, AAQ.ANX, AAQ.AV)%>%
psych::describe(na.rm=TRUE) %>%
as.data.frame() %>%
dplyr::select("n", "mean", "sd", "median", "min", "max", "range", "skew", "kurtosis") %>%
kableExtra::kable(caption= "Age and Relationship Length Descriptives", digits = 2) %>%
kable_styling(bootstrap_options = "striped", full_width = TRUE)| n | mean | sd | median | min | max | range | skew | kurtosis | |
|---|---|---|---|---|---|---|---|---|---|
| disclose.cat | 330 | 1.95 | 0.79 | 2.00 | 1.00 | 3.00 | 2.00 | 0.09 | -1.39 |
| relLength | 330 | 7.73 | 7.14 | 5.50 | 1.00 | 45.92 | 44.92 | 2.32 | 6.98 |
| Age_1 | 330 | 36.98 | 12.10 | 34.00 | 19.00 | 78.00 | 59.00 | 0.87 | 0.09 |
| PRQC | 330 | 5.72 | 1.11 | 6.00 | 1.78 | 7.00 | 5.22 | -1.19 | 1.01 |
| Trust | 330 | 5.88 | 1.38 | 6.33 | 1.00 | 7.00 | 6.00 | -1.47 | 1.63 |
| PRI_total | 330 | 3.86 | 1.21 | 4.25 | 0.00 | 5.00 | 5.00 | -1.14 | 0.53 |
| AAQ.ANX | 330 | 3.49 | 1.21 | 3.44 | 1.00 | 7.00 | 6.00 | 0.25 | -0.26 |
| AAQ.AV | 330 | 3.94 | 1.35 | 4.00 | 1.00 | 7.00 | 6.00 | -0.28 | -0.54 |
# Relationship Type
table(disclosedata$Relationship_type)##
## Casually dating Engaged Friends with benefits
## 7 15 2
## Living together Married Seriously dating
## 157 72 77prop.table(table(disclosedata$Relationship_type))##
## Casually dating Engaged Friends with benefits
## 0.021212121 0.045454545 0.006060606
## Living together Married Seriously dating
## 0.475757576 0.218181818 0.233333333# Sexual Orientation
table(disclosedata$Orientation)##
## Bisexual or pansexual Lesbian or gay other
## 110 22 7
## Sexual or heterosexual
## 191prop.table(table(disclosedata$Orientation))##
## Bisexual or pansexual Lesbian or gay other
## 0.33333333 0.06666667 0.02121212
## Sexual or heterosexual
## 0.57878788# Gender
length(which(disclosedata$CisWoman==1))## [1] 308length(which(disclosedata$GenderQueer==1))## [1] 18length(which(disclosedata$TransMan==1))## [1] 0length(which(disclosedata$TransWoman==1))## [1] 2length(which(disclosedata$Questioning==1))## [1] 5# Education
table(disclosedata$Education)##
## Associate's degree Bachelor's degree Doctorate degree
## 28 118 12
## High school degree Master's degree Some college
## 40 40 83
## Some graduate school Some high school
## 7 2prop.table(table(disclosedata$Education))##
## Associate's degree Bachelor's degree Doctorate degree
## 0.084848485 0.357575758 0.036363636
## High school degree Master's degree Some college
## 0.121212121 0.121212121 0.251515152
## Some graduate school Some high school
## 0.021212121 0.006060606# Income
table(disclosedata$Income)##
## $15,000-34,999 $35,000-49,999 $50,000-74,999 $75,000 or more
## 50 58 82 114
## Less than $15,000
## 26prop.table(table(disclosedata$Income))##
## $15,000-34,999 $35,000-49,999 $50,000-74,999 $75,000 or more
## 0.15151515 0.17575758 0.24848485 0.34545455
## Less than $15,000
## 0.07878788# Race
length(which(disclosedata$Race_1==1)) # Alaskan Native, Native American, or Native Hawaiian/Pacific Isander## [1] 3length(which(disclosedata$Race_2==1)) # Asian or Asian American## [1] 11length(which(disclosedata$Race_3==1)) # Black or African American## [1] 33length(which(disclosedata$Race_4==1)) # Latinx or Hispanic## [1] 36length(which(disclosedata$Race_5==1)) # Middle Eastern or North African## [1] 2length(which(disclosedata$Race_6==1)) # Native American, Hawaiian Native, or Alaskan Native## [1] 5length(which(disclosedata$Race_7==1)) # White or European American## [1] 277#1 and 6 accidentally repeated, fixing it
length(which(disclosedata$Race_1==1 & disclosedata$Race_6==1)) #so only 1 person aid yes to option 1 and 6## [1] 1# 5 people said yes to 6 and 3 to q1, so 5+3=8-1=7 because removing repeat
# Disclosure
# 1 = none; 2 = some; 3 = all
table(disclosedata$disclose.cat.fac)##
## 1 2 3
## 111 125 94prop.table(table(disclosedata$disclose.cat.fac))##
## 1 2 3
## 0.3363636 0.3787879 0.2848485Correlations between Study Variables
disclosecor<- disclosedata %>%
select(disclose.cat ,PRQC, Trust, PRI_total, AAQ.ANX, AAQ.AV)
(rcorr<- rcorr(as.matrix(disclosecor, type = c("spearman"))))## disclose.cat PRQC Trust PRI_total AAQ.ANX AAQ.AV
## disclose.cat 1.00 0.25 0.17 0.21 -0.04 -0.15
## PRQC 0.25 1.00 0.84 0.75 -0.34 -0.23
## Trust 0.17 0.84 1.00 0.71 -0.40 -0.26
## PRI_total 0.21 0.75 0.71 1.00 -0.39 -0.27
## AAQ.ANX -0.04 -0.34 -0.40 -0.39 1.00 0.34
## AAQ.AV -0.15 -0.23 -0.26 -0.27 0.34 1.00
##
## n= 330
##
##
## P
## disclose.cat PRQC Trust PRI_total AAQ.ANX AAQ.AV
## disclose.cat 0.0000 0.0026 0.0001 0.5145 0.0053
## PRQC 0.0000 0.0000 0.0000 0.0000 0.0000
## Trust 0.0026 0.0000 0.0000 0.0000 0.0000
## PRI_total 0.0001 0.0000 0.0000 0.0000 0.0000
## AAQ.ANX 0.5145 0.0000 0.0000 0.0000 0.0000
## AAQ.AV 0.0053 0.0000 0.0000 0.0000 0.0000# Exporting to data frame so that I can export to xl
cor.table<- (as.data.frame(rcorr$r))
cor.table<- round(cor.table, 2)
write_xlsx(cor.table, "~/Downloads/cor.tables.xlsx") #function that exports this dataframe to a xcel fileData Visualization
by(data=disclosedata$Trust, INDICES=disclosedata$disclose.cat.fac, FUN=mean,na.rm=TRUE) ## disclosedata$disclose.cat.fac: 1
## [1] 5.708709
## ------------------------------------------------------------
## disclosedata$disclose.cat.fac: 2
## [1] 5.722667
## ------------------------------------------------------------
## disclosedata$disclose.cat.fac: 3
## [1] 6.304965by(data=disclosedata$PRQC, INDICES=disclosedata$disclose.cat, FUN=mean,na.rm=TRUE) ## disclosedata$disclose.cat: 1
## [1] 5.400901
## ------------------------------------------------------------
## disclosedata$disclose.cat: 2
## [1] 5.707556
## ------------------------------------------------------------
## disclosedata$disclose.cat: 3
## [1] 6.101064by(data=disclosedata$PRI_total, INDICES=disclosedata$disclose.cat, FUN=mean,na.rm=TRUE) ## disclosedata$disclose.cat: 1
## [1] 3.599099
## ------------------------------------------------------------
## disclosedata$disclose.cat: 2
## [1] 3.809
## ------------------------------------------------------------
## disclosedata$disclose.cat: 3
## [1] 4.240691by(data=disclosedata$AAQ.ANX, INDICES=disclosedata$disclose.cat, FUN=mean,na.rm=TRUE) ## disclosedata$disclose.cat: 1
## [1] 3.483483
## ------------------------------------------------------------
## disclosedata$disclose.cat: 2
## [1] 3.583111
## ------------------------------------------------------------
## disclosedata$disclose.cat: 3
## [1] 3.362884by(data=disclosedata$AAQ.AV, INDICES=disclosedata$disclose.cat, FUN=mean,na.rm=TRUE) ## disclosedata$disclose.cat: 1
## [1] 4.112613
## ------------------------------------------------------------
## disclosedata$disclose.cat: 2
## [1] 4.062
## ------------------------------------------------------------
## disclosedata$disclose.cat: 3
## [1] 3.575798disclosedata$disclose.cat<- as.numeric(disclosedata$disclose.cat)
disclosedata %>%
group_by(POC) %>%
summarise(mean = mean(disclose.cat)) disclosedata %>%
group_by(LGBTQ) %>%
summarise(mean = mean(disclose.cat)) lapply(disclosedata[, c("disclose.cat.fac", "POC", "LGBTQ")], table) # Frequencies of factor variables## $disclose.cat.fac
##
## 1 2 3
## 111 125 94
##
## $POC
##
## POC White
## 80 250
##
## $LGBTQ
##
## LGBTQ StraightCis
## 140 190ftable(xtabs(~ POC + disclose.cat.fac, data = disclosedata)) # Frequencies POC x disclosure## disclose.cat.fac 1 2 3
## POC
## POC 21 36 23
## White 90 89 71ftable(xtabs(~ LGBTQ + disclose.cat.fac, data = disclosedata)) # Frequencies LGBTQ x disclosure## disclose.cat.fac 1 2 3
## LGBTQ
## LGBTQ 34 52 54
## StraightCis 77 73 40ggplot(disclosedata, aes(x = disclose.cat.fac, y = PRQC)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) 
ggplot(disclosedata, aes(x = disclose.cat.fac, y = Trust)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) 
ggplot(disclosedata, aes(x = disclose.cat.fac, y = PRI_total)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) 
ggplot(disclosedata, aes(x = disclose.cat.fac, y = AAQ.ANX)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
ggplot(disclosedata, aes(x = disclose.cat.fac, y = AAQ.AV)) +
geom_boxplot(size = .75) +
geom_jitter(alpha = .5) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
Detecting Outliers
PRQCmin <- mean(disclosedata$PRQC, na.rm=T) - (3*(sd(disclosedata$PRQC, na.rm=T)))
PRQCmax <- mean(disclosedata$PRQC, na.rm=T) + (3*(sd(disclosedata$PRQC, na.rm=T)))
disclosedata$PRQC[which(disclosedata$PRQC < PRQCmin | disclosedata$PRQC > PRQCmax)] # 3 outliers## [1] 2.277778 1.944444 1.777778PRImin <- mean(disclosedata$PRI_total, na.rm=T) - (3*(sd(disclosedata$PRI_total, na.rm=T)))
PRImax <- mean(disclosedata$PRI_total, na.rm=T) + (3*(sd(disclosedata$PRI_total, na.rm=T)))
disclosedata$PRI_total[which(disclosedata$PRI_total < PRImin | disclosedata$PRI_total > PRImax)] # 3 outliers## [1] 0 0 0Trustmin <- mean(disclosedata$Trust, na.rm=T) - (3*(sd(disclosedata$Trust, na.rm=T)))
Trustmax <- mean(disclosedata$Trust, na.rm=T) + (3*(sd(disclosedata$Trust, na.rm=T)))
disclosedata$Trust[which(disclosedata$Trust < Trustmin | disclosedata$Trust > Trustmax)] # 5 outliers## [1] 1.000000 1.666667 1.666667 1.333333 1.000000AAQ.AVmin <- mean(disclosedata$AAQ.AV, na.rm=T) - (3*(sd(disclosedata$AAQ.AV, na.rm=T)))
AAQ.AVmax <- mean(disclosedata$AAQ.AV, na.rm=T) + (3*(sd(disclosedata$AAQ.AV, na.rm=T)))
disclosedata$AAQ.AV[which(disclosedata$AAQ.AV < AAQ.AVmin | disclosedata$AAQ.AV > AAQ.AVmax)] # 0 outliers## numeric(0)AAQ.ANXmin <- mean(disclosedata$AAQ.ANX, na.rm=T) - (3*(sd(disclosedata$AAQ.ANX, na.rm=T)))
AAQ.ANXmax <- mean(disclosedata$AAQ.ANX, na.rm=T) + (3*(sd(disclosedata$AAQ.ANX, na.rm=T)))
disclosedata$AAQ.ANX[which(disclosedata$AAQ.ANX < AAQ.ANXmin | disclosedata$AAQ.ANX > AAQ.ANXmax)] # 0 outliers## numeric(0)# Winsoring Outliers to +- 3 sds from mean
disclosedata$PRQC.win<- winsorize(disclosedata$PRQC, method="zscore",
threshold=3, robust=FALSE)
disclosedata$PRI.win<- winsorize(disclosedata$PRI_total, method="zscore",
threshold=3, robust=FALSE)
disclosedata$Trust.win<- winsorize(disclosedata$Trust, method="zscore",
threshold=3, robust=FALSE) Mean Centering
disclosedata$cPRQC<- disclosedata$PRQC.win - mean(disclosedata$PRQC.win, na.rm=TRUE)
disclosedata$cPRI<- disclosedata$PRI.win - mean(disclosedata$PRI.win, na.rm=TRUE)
disclosedata$cTrust<- disclosedata$Trust.win - mean(disclosedata$Trust.win, na.rm=TRUE)
disclosedata$cANX<- disclosedata$AAQ.ANX - mean(disclosedata$AAQ.ANX, na.rm=TRUE)
disclosedata$cAV<- disclosedata$AAQ.AV - mean(disclosedata$AAQ.AV, na.rm=TRUE)Linear Regression
lin.reg<- lm(disclosure.ratio ~ POC + LGBTQ + cTrust + cPRQC + cPRI + cAV + cANX, data=disclosedata)
summary(lin.reg) #Looks like a spurious coefficient##
## Call:
## lm(formula = disclosure.ratio ~ POC + LGBTQ + cTrust + cPRQC +
## cPRI + cAV + cANX, data = disclosedata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.43624 -0.65106 -0.08813 0.69214 1.75922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.16742 0.09774 22.175 < 2e-16 ***
## POCWhite -0.08245 0.09912 -0.832 0.40611
## LGBTQStraightCis -0.28975 0.08856 -3.272 0.00118 **
## cTrust -0.15633 0.07492 -2.087 0.03771 *
## cPRQC 0.27056 0.08405 3.219 0.00142 **
## cPRI 0.08639 0.06362 1.358 0.17543
## cAV -0.07498 0.03408 -2.200 0.02851 *
## cANX 0.02135 0.04092 0.522 0.60220
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7649 on 322 degrees of freedom
## Multiple R-squared: 0.1318, Adjusted R-squared: 0.1129
## F-statistic: 6.983 on 7 and 322 DF, p-value: 9.323e-08# Does not seem to meet assumptions. Ratios as DVs can act strangely. This might be because of the distribution of disclosure ratios (U curve). Proportional Odds Logistic Regression (Ordinal Logistic Regression)
# https://www.bookdown.org/rwnahhas/RMPH/blr-ordinal.html
# https://stats.oarc.ucla.edu/r/dae/ordinal-logistic-regression/
# Hess = TRUE produces a matrix necessary to calculate the standard errors
# polr function is from MASS package
# Make SURE that variables are factors and not characters, this will change results
disclosedata$LGBTQ <- relevel(disclosedata$LGBTQ, ref = "StraightCis") # Setting straightcis as #reference level
disclosedata$POC <- relevel(disclosedata$POC, ref = "POC") # Setting POC as reference level
polr<- MASS:: polr(disclose.cat.fac ~ POC + LGBTQ + cTrust + cPRQC + cPRI + cAV + cANX,
data=disclosedata,
Hess=TRUE)
summary(polr)## Call:
## MASS::polr(formula = disclose.cat.fac ~ POC + LGBTQ + cTrust +
## cPRQC + cPRI + cAV + cANX, data = disclosedata, Hess = TRUE)
##
## Coefficients:
## Value Std. Error t value
## POCWhite -0.14249 0.24039 -0.5928
## LGBTQLGBTQ 0.80084 0.21955 3.6477
## cTrust -0.33989 0.18888 -1.7995
## cPRQC 0.58621 0.21080 2.7809
## cPRI 0.19640 0.15506 1.2666
## cAV -0.20866 0.08464 -2.4653
## cANX 0.09677 0.10360 0.9341
##
## Intercepts:
## Value Std. Error t value
## 1|2 -0.5220 0.2393 -2.1815
## 2|3 1.2556 0.2478 5.0678
##
## Residual Deviance: 677.4246
## AIC: 695.4246tbl_regression(polr, exp=TRUE) # Odds ratios| Characteristic | OR1 | 95% CI1 |
|---|---|---|
| POC | ||
| POC | — | — |
| White | 0.87 | 0.54, 1.39 |
| LGBTQ | ||
| StraightCis | — | — |
| LGBTQ | 2.23 | 1.45, 3.44 |
| cTrust | 0.71 | 0.49, 1.03 |
| cPRQC | 1.80 | 1.19, 2.73 |
| cPRI | 1.22 | 0.90, 1.65 |
| cAV | 0.81 | 0.69, 0.96 |
| cANX | 1.10 | 0.90, 1.35 |
| 1 OR = Odds Ratio, CI = Confidence Interval | ||
## Store coefficients
(coef<- coef(summary(polr)))## Value Std. Error t value
## POCWhite -0.14249397 0.24039190 -0.5927569
## LGBTQLGBTQ 0.80084180 0.21954689 3.6477028
## cTrust -0.33989149 0.18887660 -1.7995426
## cPRQC 0.58620661 0.21079701 2.7809057
## cPRI 0.19640037 0.15506041 1.2666055
## cAV -0.20866450 0.08463893 -2.4653489
## cANX 0.09677394 0.10360001 0.9341113
## 1|2 -0.52199020 0.23928492 -2.1814588
## 2|3 1.25564536 0.24777029 5.0677802## Calculate and store p-values
p <- pnorm(abs(coef[, "t value"]), lower.tail = FALSE) * 2
## Table of coefficients and pvalues
(ctable <- cbind(coef, "p value" = p))## Value Std. Error t value p value
## POCWhite -0.14249397 0.24039190 -0.5927569 5.533438e-01
## LGBTQLGBTQ 0.80084180 0.21954689 3.6477028 2.645955e-04
## cTrust -0.33989149 0.18887660 -1.7995426 7.193289e-02
## cPRQC 0.58620661 0.21079701 2.7809057 5.420748e-03
## cPRI 0.19640037 0.15506041 1.2666055 2.052964e-01
## cAV -0.20866450 0.08463893 -2.4653489 1.368799e-02
## cANX 0.09677394 0.10360001 0.9341113 3.502465e-01
## 1|2 -0.52199020 0.23928492 -2.1814588 2.914950e-02
## 2|3 1.25564536 0.24777029 5.0677802 4.024818e-07# Confidence intervals
(ci <- confint(polr))## Waiting for profiling to be done...## 2.5 % 97.5 %
## POCWhite -0.6145584 0.32912178
## LGBTQLGBTQ 0.3729102 1.23437953
## cTrust -0.7144500 0.02762468
## cPRQC 0.1754660 1.00332142
## cPRI -0.1072886 0.50192950
## cAV -0.3759839 -0.04368108
## cANX -0.1057590 0.30110670# Proportional Odds Assumption
# Whether odds of moving from one category to the next is the same across all (DV) categroes
# Ho: Parallel assumption holds
# Ha: Parallel assumption does not hold
poTest(polr)##
## Tests for Proportional Odds
## MASS::polr(formula = disclose.cat.fac ~ POC + LGBTQ + cTrust +
## cPRQC + cPRI + cAV + cANX, data = disclosedata, Hess = TRUE)
##
## b[polr] b[>1] b[>2] Chisquare df Pr(>Chisq)
## Overall 11.13 7 0.133
## POCWhite -0.14249 -0.34346 0.00377 1.07 1 0.301
## LGBTQLGBTQ 0.80084 0.69396 0.86043 0.34 1 0.562
## cTrust -0.33989 -0.58529 -0.00513 5.52 1 0.019 *
## cPRQC 0.58621 0.77355 0.36401 2.32 1 0.128
## cPRI 0.19640 0.16627 0.20415 0.03 1 0.861
## cAV -0.20866 -0.16743 -0.26164 0.74 1 0.389
## cANX 0.09677 0.09678 0.06985 0.04 1 0.838
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1brant(polr) # Also trying Brant-Wald test (Brant, 1990), according to both assumption of proportional odds was not met## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 11.13 7 0.13
## POCWhite 1.07 1 0.3
## LGBTQLGBTQ 0.34 1 0.56
## cTrust 5.52 1 0.02
## cPRQC 2.32 1 0.13
## cPRI 0.03 1 0.86
## cAV 0.74 1 0.39
## cANX 0.04 1 0.84
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds# These results suggest a violation of the proportional odds assumption for trustPartial Proportional Odds Logistic Regression
# As we could see, the trust variable did not meet the assumption of proportional odds so we will try partial proportional odds which relaxes the assumption of proportional odds for variables that did not meet that assumption
part_polr <- vglm(disclose.cat.fac ~ POC + LGBTQ + cTrust + cPRQC + cPRI + cAV + cANX,
data = disclosedata,
family = cumulative(parallel = FALSE ~ cTrust, # relaxing parallel odds assumption for trust
reverse = TRUE))## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered()summary(part_polr)##
## Call:
## vglm(formula = disclose.cat.fac ~ POC + LGBTQ + cTrust + cPRQC +
## cPRI + cAV + cANX, family = cumulative(parallel = FALSE ~
## cTrust, reverse = TRUE), data = disclosedata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.51062 0.24125 2.117 0.034296 *
## (Intercept):2 -1.29475 0.25260 -5.126 2.96e-07 ***
## POCWhite -0.15832 0.24395 -0.649 0.516326
## LGBTQLGBTQ 0.81006 0.22144 3.658 0.000254 ***
## cTrust:1 -0.49014 0.20079 -2.441 0.014643 *
## cTrust:2 -0.12110 0.20801 -0.582 0.560433
## cPRQC 0.58562 0.21103 2.775 0.005519 **
## cPRI 0.20769 0.15870 1.309 0.190642
## cAV -0.20105 0.08493 -2.367 0.017919 *
## cANX 0.08938 0.10152 0.880 0.378610
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 671.524 on 650 degrees of freedom
##
## Log-likelihood: -335.762 on 650 degrees of freedom
##
## Number of Fisher scoring iterations: 5
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## POCWhite LGBTQLGBTQ cTrust:1 cTrust:2 cPRQC cPRI cAV
## 0.8535724 2.2480446 0.6125407 0.8859438 1.7961022 1.2308306 0.8178729
## cANX
## 1.0934986round(confint(part_polr), 2)## 2.5 % 97.5 %
## (Intercept):1 0.04 0.98
## (Intercept):2 -1.79 -0.80
## POCWhite -0.64 0.32
## LGBTQLGBTQ 0.38 1.24
## cTrust:1 -0.88 -0.10
## cTrust:2 -0.53 0.29
## cPRQC 0.17 1.00
## cPRI -0.10 0.52
## cAV -0.37 -0.03
## cANX -0.11 0.29round(coef(part_polr),2)## (Intercept):1 (Intercept):2 POCWhite LGBTQLGBTQ cTrust:1
## 0.51 -1.29 -0.16 0.81 -0.49
## cTrust:2 cPRQC cPRI cAV cANX
## -0.12 0.59 0.21 -0.20 0.09round(exp(coef(part_polr)),2)## (Intercept):1 (Intercept):2 POCWhite LGBTQLGBTQ cTrust:1
## 1.67 0.27 0.85 2.25 0.61
## cTrust:2 cPRQC cPRI cAV cANX
## 0.89 1.80 1.23 0.82 1.09confint(part_polr)## 2.5 % 97.5 %
## (Intercept):1 0.03778278 0.98345183
## (Intercept):2 -1.78983278 -0.79965743
## POCWhite -0.63644868 0.31979884
## LGBTQLGBTQ 0.37605264 1.24406891
## cTrust:1 -0.88367415 -0.09660570
## cTrust:2 -0.52878913 0.28658557
## cPRQC 0.17201104 0.99922672
## cPRI -0.10335883 0.51873722
## cAV -0.36750300 -0.03459361
## cANX -0.10958884 0.28835344### INTERPRETATION
# Coefficient #1: Cat 1 vs Cat 2 + 3
# Coefficient #2: Cat 1 + Cat 2 vs Cat 3
# Predicted Probabilities
newdat <- data.frame(
LGBTQ = rep(0:1, each = 200),
PRQC.win = rep(seq(from = 1, to = 7, length.out = 100), 4))
Effect(focal.predictors = c("LGBTQ"), polr)##
## LGBTQ effect (probability) for 1
## LGBTQ
## StraightCis LGBTQ
## 0.3979437 0.2288384
##
## LGBTQ effect (probability) for 2
## LGBTQ
## StraightCis LGBTQ
## 0.3983997 0.4082510
##
## LGBTQ effect (probability) for 3
## LGBTQ
## StraightCis LGBTQ
## 0.2036566 0.3629106plot(Effect(focal.predictors = c("cPRQC"), polr), rug = F)
plot(Effect(focal.predictors = c("cTrust"), polr), rug = F)
plot(Effect(focal.predictors = c("cAV"), polr), rug = F)
plot(Effect(focal.predictors = c("LGBTQ"), polr), rug = F)
# Assumption of multicollinearity
vif(polr) # All values are under 10 so multicollinearity does not seem to be an issue## Warning in vif.default(polr): No intercept: vifs may not be sensible.## POC LGBTQ cTrust cPRQC cPRI cAV cANX
## 1.911551 3.142458 3.321454 2.336421 1.212377 1.403616 4.590137# calculate the percent change in odds by taking (OR - 1) * 100.
(0.4448310 - 1) * 100 # LGBTQ = -55.5169## [1] -55.5169(0.6125407 - 1) * 100 # Trust1 = -38.75## [1] -38.74593(1.7961022 - 1) * 100 # PRQC## [1] 79.61022(0.8178729 - 1) * 100 # Attachment Avoidance## [1] -18.21271# Running separate model with only trust. When trust is not conditioned on other variables in the model, more trust is assocated with greater likelihood of disclosing.
trust.part_polr <- vglm(disclose.cat.fac ~ cTrust,
data = disclosedata,
family = cumulative(parallel = FALSE,
reverse = TRUE))## Warning in eval(slot(family, "initialize")): response should be ordinal---see
## ordered() summary(trust.part_polr)##
## Call:
## vglm(formula = disclose.cat.fac ~ cTrust, family = cumulative(parallel = FALSE,
## reverse = TRUE), data = disclosedata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.6821 0.1169 5.837 5.33e-09 ***
## (Intercept):2 -0.9695 0.1277 -7.594 3.11e-14 ***
## cTrust:1 0.1408 0.1185 1.188 0.234822
## cTrust:2 0.4882 0.1403 3.479 0.000502 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3])
##
## Residual deviance: 707.6957 on 656 degrees of freedom
##
## Log-likelihood: -353.8479 on 656 degrees of freedom
##
## Number of Fisher scoring iterations: 4
##
## No Hauck-Donner effect found in any of the estimates
##
##
## Exponentiated coefficients:
## cTrust:1 cTrust:2
## 1.151190 1.629302AIC(trust.part_polr)## [1] 715.6957Interpretation: For every one unit increase in relationship quality, we would expect 1.80 greater odds that they disclose to their romantic partner. For every unit increase in attachment avoidance, we would expect a .82 odds decrease in the likelihood that they disclose to their romantic partner. LGBTQ individuals are 2.25 greater odds of disclosing their sexual assault to their romantic partners compared to straight and cisgender individuals, while holding other variables constant.
Multinomial Logistic Regression
# Reference category is "All"
# https://stats.oarc.ucla.edu/r/dae/multinomial-logistic-regression/
# multinom is from the nnet package
mult.log <- multinom(disclose.cat ~ POC + LGBTQ + cTrust + cPRQC + cPRI + cAV + cANX, data = disclosedata)## # weights: 27 (16 variable)
## initial value 362.542055
## iter 10 value 339.183228
## iter 20 value 333.474869
## final value 333.447359
## convergedsummary(mult.log)## Call:
## multinom(formula = disclose.cat ~ POC + LGBTQ + cTrust + cPRQC +
## cPRI + cAV + cANX, data = disclosedata)
##
## Coefficients:
## (Intercept) POCWhite LGBTQLGBTQ cTrust cPRQC cPRI cAV
## 2 0.3338050 -0.4240286 0.4091446 -0.6955341 0.7788146 0.1123150 -0.0717473
## 3 -0.4979546 -0.2444012 1.0920859 -0.3819082 0.7737413 0.2622391 -0.3026483
## cANX
## 2 0.0836727
## 3 0.1171934
##
## Std. Errors:
## (Intercept) POCWhite LGBTQLGBTQ cTrust cPRQC cPRI cAV
## 2 0.3069511 0.3213405 0.2913154 0.2454753 0.2708426 0.1993778 0.1114797
## 3 0.3519340 0.3598937 0.3171866 0.2767686 0.3035812 0.2351878 0.1215819
## cANX
## 2 0.132131
## 3 0.145844
##
## Residual Deviance: 666.8947
## AIC: 698.8947# Coefficients
z <- summary(mult.log)$coefficients/summary(mult.log)$standard.errors
z## (Intercept) POCWhite LGBTQLGBTQ cTrust cPRQC cPRI cAV
## 2 1.087486 -1.3195618 1.404473 -2.833418 2.875525 0.5633273 -0.6435909
## 3 -1.414909 -0.6790929 3.443039 -1.379883 2.548713 1.1150197 -2.4892541
## cANX
## 2 0.6332557
## 3 0.8035532# 2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1)) * 2 # wald tests
p## (Intercept) POCWhite LGBTQLGBTQ cTrust cPRQC cPRI cAV
## 2 0.2768221 0.1869814 0.160178070 0.004605307 0.004033559 0.573212 0.51984074
## 3 0.1570952 0.4970790 0.000575216 0.167622799 0.010812139 0.264842 0.01280114
## cANX
## 2 0.5265667
## 3 0.4216551## Risk ratios: extract the coefficients from the model and exponentiate
exp(coef(mult.log))## (Intercept) POCWhite LGBTQLGBTQ cTrust cPRQC cPRI cAV
## 2 1.3962709 0.6544052 1.505529 0.4988080 2.178888 1.118865 0.9307661
## 3 0.6077725 0.7831733 2.980485 0.6825577 2.167862 1.299837 0.7388589
## cANX
## 2 1.087273
## 3 1.124337tbl_regression(mult.log, exp = TRUE) # ORs## ℹ Multinomial models have a different underlying structure than the models
## gtsummary was designed for. Other gtsummary functions designed to work with
## tbl_regression objects may yield unexpected results.| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| 2 | |||
| POC | |||
| POC | — | — | |
| White | 0.65 | 0.35, 1.23 | 0.2 |
| LGBTQ | |||
| StraightCis | — | — | |
| LGBTQ | 1.51 | 0.85, 2.66 | 0.2 |
| cTrust | 0.50 | 0.31, 0.81 | 0.005 |
| cPRQC | 2.18 | 1.28, 3.70 | 0.004 |
| cPRI | 1.12 | 0.76, 1.65 | 0.6 |
| cAV | 0.93 | 0.75, 1.16 | 0.5 |
| cANX | 1.09 | 0.84, 1.41 | 0.5 |
| 3 | |||
| POC | |||
| POC | — | — | |
| White | 0.78 | 0.39, 1.59 | 0.5 |
| LGBTQ | |||
| StraightCis | — | — | |
| LGBTQ | 2.98 | 1.60, 5.55 | <0.001 |
| cTrust | 0.68 | 0.40, 1.17 | 0.2 |
| cPRQC | 2.17 | 1.20, 3.93 | 0.011 |
| cPRI | 1.30 | 0.82, 2.06 | 0.3 |
| cAV | 0.74 | 0.58, 0.94 | 0.013 |
| cANX | 1.12 | 0.84, 1.50 | 0.4 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
tidy(mult.log, conf.int = TRUE, exponentiate = TRUE) %>%
kable() %>%
kable_styling("basic", full_width = FALSE)| y.level | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|
| 2 | (Intercept) | 1.3962709 | 0.3069511 | 1.0874859 | 0.2768221 | 0.7650545 | 2.5482789 |
| 2 | POCWhite | 0.6544052 | 0.3213405 | -1.3195618 | 0.1869814 | 0.3485951 | 1.2284917 |
| 2 | LGBTQLGBTQ | 1.5055295 | 0.2913154 | 1.4044729 | 0.1601781 | 0.8505916 | 2.6647557 |
| 2 | cTrust | 0.4988080 | 0.2454753 | -2.8334184 | 0.0046053 | 0.3083078 | 0.8070163 |
| 2 | cPRQC | 2.1788879 | 0.2708426 | 2.8755252 | 0.0040336 | 1.2814253 | 3.7049000 |
| 2 | cPRI | 1.1188652 | 0.1993778 | 0.5633273 | 0.5732120 | 0.7569498 | 1.6538209 |
| 2 | cAV | 0.9307661 | 0.1114797 | -0.6435909 | 0.5198407 | 0.7480814 | 1.1580631 |
| 2 | cANX | 1.0872730 | 0.1321310 | 0.6332557 | 0.5265667 | 0.8392058 | 1.4086681 |
| 3 | (Intercept) | 0.6077725 | 0.3519340 | -1.4149091 | 0.1570952 | 0.3049118 | 1.2114565 |
| 3 | POCWhite | 0.7831733 | 0.3598937 | -0.6790929 | 0.4970790 | 0.3868262 | 1.5856230 |
| 3 | LGBTQLGBTQ | 2.9804845 | 0.3171866 | 3.4430391 | 0.0005752 | 1.6006528 | 5.5497906 |
| 3 | cTrust | 0.6825577 | 0.2767686 | -1.3798826 | 0.1676228 | 0.3967834 | 1.1741545 |
| 3 | cPRQC | 2.1678618 | 0.3035812 | 2.5487125 | 0.0108121 | 1.1957015 | 3.9304331 |
| 3 | cPRI | 1.2998373 | 0.2351878 | 1.1150197 | 0.2648420 | 0.8197789 | 2.0610154 |
| 3 | cAV | 0.7388589 | 0.1215819 | -2.4892541 | 0.0128011 | 0.5821981 | 0.9376748 |
| 3 | cANX | 1.1243368 | 0.1458440 | 0.8035532 | 0.4216551 | 0.8447999 | 1.4963701 |
# Now using info from https://bookdown.org/sarahwerth2024/CategoricalBook/multinomial-logit-regression-r.html#running-a-mlr-in-r
# OR do table with all info
tidy(mult.log, conf.int = TRUE)# relative risk ratios
tbl_regression(mult.log, exp = TRUE)## ℹ Multinomial models have a different underlying structure than the models
## gtsummary was designed for. Other gtsummary functions designed to work with
## tbl_regression objects may yield unexpected results.| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| 2 | |||
| POC | |||
| POC | — | — | |
| White | 0.65 | 0.35, 1.23 | 0.2 |
| LGBTQ | |||
| StraightCis | — | — | |
| LGBTQ | 1.51 | 0.85, 2.66 | 0.2 |
| cTrust | 0.50 | 0.31, 0.81 | 0.005 |
| cPRQC | 2.18 | 1.28, 3.70 | 0.004 |
| cPRI | 1.12 | 0.76, 1.65 | 0.6 |
| cAV | 0.93 | 0.75, 1.16 | 0.5 |
| cANX | 1.09 | 0.84, 1.41 | 0.5 |
| 3 | |||
| POC | |||
| POC | — | — | |
| White | 0.78 | 0.39, 1.59 | 0.5 |
| LGBTQ | |||
| StraightCis | — | — | |
| LGBTQ | 2.98 | 1.60, 5.55 | <0.001 |
| cTrust | 0.68 | 0.40, 1.17 | 0.2 |
| cPRQC | 2.17 | 1.20, 3.93 | 0.011 |
| cPRI | 1.30 | 0.82, 2.06 | 0.3 |
| cAV | 0.74 | 0.58, 0.94 | 0.013 |
| cANX | 1.12 | 0.84, 1.50 | 0.4 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
vif(mult.log)## Warning in vif.default(mult.log): No intercept: vifs may not be sensible.## POC LGBTQ cTrust cPRQC cPRI cAV cANX
## 5.727115 2.660958 4.362215 4.393048 3.135980 1.639429 1.901750# Model comparison
AIC(part_polr) # Best model## [1] 691.524AIC(polr)## [1] 695.4246AIC(mult.log)## [1] 698.8947AIC(lin.reg) # Linear regression is worst fitting model## [1] 769.5168## So best fitting model is the partial proportional odds modelQualitative Data Frequencies
# Merging sexual orientation and qualitative data data
qual_dat <- read_excel("~/Downloads/qual_dat.xlsx")
qual_dat.complete<-merge(x=disclosedata, y=qual_dat, by='SubjectID', all.x=T) #merging data frames
# Frequencies
length(which(qual_dat.complete$Emotional.support==1))## [1] 149length(which(qual_dat.complete$Treat.differently==1))## [1] 1length(which(qual_dat.complete$Distraction.Dismissive==1))## [1] 13length(which(qual_dat.complete$Take.control==1))## [1] 0length(which(qual_dat.complete$Tangible.aid==1))## [1] 4length(which(qual_dat.complete$Victim.blame==1))## [1] 2length(which(qual_dat.complete$angry.sad.upset==1))## [1] 55length(which(qual_dat.complete$Partner.protective==1))## [1] 12length(which(qual_dat.complete$topic.for.therapist==1))## [1] 4length(which(qual_dat.complete$Too.hard.talk.about==1))## [1] 34length(which(qual_dat.complete$Partner.unsupportive==1))## [1] 10length(which(qual_dat.complete$Worried.what.partner.think==1))## [1] 9length(which(qual_dat.complete$keep.things.to.myself==1))## [1] 12length(which(qual_dat.complete$upset.partner==1))## [1] 9length(which(qual_dat.complete$embarrassed.ashamed==1))## [1] 19length(which(qual_dat.complete$not.important==1))## [1] 34length(which(qual_dat.complete$Other==1))## [1] 2length(which(qual_dat.complete$Partner.did.this==1))## [1] 0length(which(qual_dat.complete$Not.enough.info==1))## [1] 44# Percentages (frequency/330)
length(which(qual_dat.complete$Emotional.support==1))/330## [1] 0.4515152length(which(qual_dat.complete$Treat.differently==1))/330## [1] 0.003030303length(which(qual_dat.complete$Distraction.Dismissive==1))/330## [1] 0.03939394length(which(qual_dat.complete$Take.control==1))/330## [1] 0length(which(qual_dat.complete$Tangible.aid==1))/330## [1] 0.01212121length(which(qual_dat.complete$Victim.blame==1))/330## [1] 0.006060606length(which(qual_dat.complete$angry.sad.upset==1))/330## [1] 0.1666667length(which(qual_dat.complete$Partner.protective==1))/330## [1] 0.03636364length(which(qual_dat.complete$topic.for.therapist==1))/330## [1] 0.01212121length(which(qual_dat.complete$Too.hard.talk.about==1))/330## [1] 0.1030303length(which(qual_dat.complete$Partner.unsupportive==1))/330## [1] 0.03030303length(which(qual_dat.complete$Worried.what.partner.think==1))/330## [1] 0.02727273length(which(qual_dat.complete$keep.things.to.myself==1))/330## [1] 0.03636364length(which(qual_dat.complete$upset.partner==1))/330## [1] 0.02727273length(which(qual_dat.complete$embarrassed.ashamed==1))/330## [1] 0.05757576length(which(qual_dat.complete$not.important==1))/330## [1] 0.1030303length(which(qual_dat.complete$Other==1))/330## [1] 0.006060606length(which(qual_dat.complete$Partner.did.this==1))/330## [1] 0length(which(qual_dat.complete$Not.enough.info==1))/330## [1] 0.1333333# exporting to de-identified data file
qual.describe.disclosure<- qual_dat.complete %>%
select(SubjectID, Relationship_type, Race_1, Race_2, Race_3, Race_4, Race_5, Race_6, Race_7, CisWoman, GenderQueer, TransMan, TransWoman, Questioning, Orientation, disclose.cat)
write_xlsx(qual.describe.disclosure, "~/Downloads/SA_disclose_qualitative.xlsx")