#factors on full dataset (three waves)
mydata$sex.f <- factor(mydata$sex,levels = c(1,2),labels = c("male", "female"))
mydata$intersex.f <- factor(mydata$intersex,levels = c(1,2,3),labels = c("yes", "no", "unsure"))
mydata$gender.f <- factor(mydata$gender,levels = c(1,2, 4),labels = c("man/boy", "woman/girl", "not listed"))
mydata$gender.everyday.f <- factor(mydata$gender.everyday,levels = c(1,2,3,4),labels = c("man/boy", "woman/girl", "changes", "other"))
mydata$indigenous.f <- factor(mydata$indigenous,levels = c(1,2),labels = c("indigenous", "not indigenous"))
mydata$black.identity.f <- factor(mydata$black.identity,levels = c(1,2,3,4),labels = c("black african", "black canadian", "african-american", "black caribbean"))
mydata$poc.f <- factor(mydata$poc,levels = c(1,2),labels = c("person of color", "not a person of color"))
mydata$disability.f <- factor(mydata$disability,levels = c(1,2),labels = c("disability", "no disability"))
mydata$disability.visible.f <- factor(mydata$dis.visible,levels = c(1,2,3),labels = c("visible always", "visible sometimes", "non-visible"))
mydata$gsa.f <- factor(mydata$gsa,levels = c(1,2),labels = c("part of a gsa", "not part of a gsa"))
mydata$time.f <- factor(mydata$time,levels = c(1,2,3),labels = c("pre", "post", "followup"))
mydata$dated.f <- factor(mydata$dated,levels = c(1,2),labels = c("relationship", "single"))
mydata1 %>% 
  select(age, sex.f, gender.f, indigenous.f, disability.f, disability.visible.f, gsa.f, dated.f) %>%  
  tbl_summary(statistic = list(
      all_categorical() ~ "{n} ({p}%)"),
            missing_text = "(Missing, prefer not to answer, or not applicable)",
            label = list(
      age ~ "Age",
      sex.f ~ "Sex at Birth",
      gender.f ~ "Gender Identity",
      indigenous.f ~ "Indigenous Identity",
      disability.f ~ "Disability status",
      disability.visible.f ~ "Disability visibility",
      gsa.f ~ "GSA membership",
      dated.f ~ "Dated in past year"))
Characteristic N = 211
Age
    13 2 (9.5%)
    14 5 (24%)
    15 7 (33%)
    16 5 (24%)
    17 2 (9.5%)
Sex at Birth
    female 19 (90%)
    male 2 (9.5%)
Gender Identity
    man/boy 8 (38%)
    not listed 5 (24%)
    woman/girl 8 (38%)
Indigenous Identity
    indigenous 1 (4.8%)
    not indigenous 20 (95%)
Disability status
    disability 10 (48%)
    no disability 11 (52%)
Disability visibility
    non-visible 6 (43%)
    visible always 1 (7.1%)
    visible sometimes 7 (50%)
    (Missing, prefer not to answer, or not applicable) 7
GSA membership
    not part of a gsa 1 (4.8%)
    part of a gsa 20 (95%)
Dated in past year
    relationship 14 (67%)
    single 7 (33%)
1 n (%) 
mydata1 %>% 
  select(bisexual, gay, lesbian, asexual, pansexual, queer, straight, two.spirit) %>%  
  tbl_summary(statistic = list(
      all_categorical() ~ "{n}"),
            missing_text = "(Missing, prefer not to answer, or not applicable)")
Characteristic N = 211
bisexual 1
    (Missing, prefer not to answer, or not applicable) 20
gay 1
    (Missing, prefer not to answer, or not applicable) 20
lesbian 2
    (Missing, prefer not to answer, or not applicable) 19
asexual 3
    (Missing, prefer not to answer, or not applicable) 18
pansexual 1
    (Missing, prefer not to answer, or not applicable) 20
queer 2
    (Missing, prefer not to answer, or not applicable) 19
straight 13
    (Missing, prefer not to answer, or not applicable) 8
two.spirit 0
    (Missing, prefer not to answer, or not applicable) 21
1 n 
mydata1 %>% 
  select(sears, self.efficacy, wellbeing, depression) %>%  
  tbl_summary(statistic = list(
      all_continuous() ~ "{mean} ({sd})   [{min}, {max}]",
      all_categorical() ~ "{n} ({p}%)"),
            missing_text = "(Missing, prefer not to answer, or not applicable)",
            label = list(
      sears ~ "Social Emotional Assets and Resilience Scale",
      self.efficacy ~ "Self Efficacy",
      wellbeing ~ "Positive Mental Health",
      depression ~ "Depression"))
Characteristic N = 211
Social Emotional Assets and Resilience Scale 2.49 (0.44) [1.66, 3.06]
    (Missing, prefer not to answer, or not applicable) 3
Self Efficacy 3.83 (0.57) [2.77, 4.92]
    (Missing, prefer not to answer, or not applicable) 3
Positive Mental Health 2.74 (1.04) [1.00, 4.85]
    (Missing, prefer not to answer, or not applicable) 5
Depression 2.71 (0.91) [1.00, 4.00]
    (Missing, prefer not to answer, or not applicable) 5
1 Mean (SD) [Range] 
LGBIS <- dplyr::select(mydata1, lgbis1,lgbis2R, lgbis3R, lgbis4R, lgbis5, lgbis6,
                                                             lgbis7R, lgbis8R, lgbis9R, lgbis10R, lgbis11R,
                                                             lgbis12, lgbis13R)
psych::alpha(LGBIS, check.keys=TRUE) #alpha = .82

QPIAS <- dplyr::select(mydata1, qpias1R, qpias2, qpias3, qpias4R, qpias5, qpias6,
                                                        qpias7, qpias8, qpias9R, qpias10, qpias11, qpias12)
psych::alpha(QPIAS, check.keys=TRUE) #alpha = .87

TIS <- dplyr::select(mydata1, tis1R, tis2, tis3, tis4R, tis5R, tis6R, tis7R, tis7R,
                                                                                                                     tis8, tis9R, tis10R, tis11, tis12R, tis13, tis14R,
                                                                                                                     tis15R, tis16R, tis17, tis18R, tis19R, tis20R, tis21,
                                                                                                                     tis22, tis23, tis24, tis25R, tis26, tis27, tis28, tis29R,
                                                                                                                     tis30R, tis31)
psych::alpha(TIS,check.keys=TRUE) #alpha = .92

SEARS <- dplyr::select(mydata1, sears1:sears35)
psych::alpha(SEARS,check.keys=TRUE) #alpha = .93

SCBS <- dplyr::select(mydata1, scbs1:scbs10)
psych::alpha(SCBS,check.keys=TRUE) #alpha = .81

CADR <- dplyr::select(mydata1, cadr2, cadr4, cadr6, cadr8, cadr10, cadr12, cadr14, cadr16, cadr18, cadr20,
                                            cadr1, cadr3, cadr5, cadr7, cadr9, cadr11, cadr13, cadr15, cadr17,cadr19)
psych::alpha(CADR,check.keys=TRUE) #alpha = .79

EIS <- dplyr::select(mydata1, eis1:eis6)
psych::alpha(EIS,check.keys=TRUE) #alpha = .94

YRBS <- dplyr::select(mydata1, yrbs1, yrbs2, yrbs3, yrbs4)
psych::alpha(YRBS,check.keys=TRUE) #alpha = .81

SE <- dplyr::select(mydata1, hrp.se1:hrp.se13)
psych::alpha(SE,check.keys=TRUE) #alpha = .87

BEHAV <- dplyr::select(mydata1, hrp.behav1:hrp.behav22)
psych::alpha(BEHAV,check.keys=TRUE) #alpha = .91

MHC <- dplyr::select(mydata1, mhc1:mhc14)
psych::alpha(MHC,check.keys=TRUE) #alpha = .94
tapply(mydata$sears, mydata$time.f, summary) 
## $pre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.657   2.151   2.500   2.485   2.871   3.057       3 
## 
## $post
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.029   2.307   2.657   2.639   3.007   3.200       3 
## 
## $followup
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.057   2.193   2.400   2.429   2.629   2.886       1
tapply(mydata$self.efficacy, mydata$time.f, summary) 
## $pre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.769   3.538   3.731   3.833   4.269   4.923       3 
## 
## $post
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   3.154   3.923   4.115   4.019   4.250   4.538       3 
## 
## $followup
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.923   3.769   4.000   3.808   4.000   4.231       1
tapply(mydata$wellbeing, mydata$time.f, summary) 
## $pre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   2.286   2.643   2.745   3.304   4.846       5 
## 
## $post
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      NA      NA      NA     NaN      NA      NA      11 
## 
## $followup
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.857   2.357   2.500   2.614   3.071   3.286       2
tapply(mydata$depression, mydata$time.f, summary) 
## $pre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   1.839   2.893   2.705   3.429   4.000       5 
## 
## $post
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      NA      NA      NA     NaN      NA      NA      11 
## 
## $followup
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.714   3.357   3.679   3.482   3.804   3.857       3
tapply(mydata$wellbeing, mydata$time.f, summary) 
## $pre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   2.286   2.643   2.745   3.304   4.846       5 
## 
## $post
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      NA      NA      NA     NaN      NA      NA      11 
## 
## $followup
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.857   2.357   2.500   2.614   3.071   3.286       2
# slight decrease in average wellbeing from pre to followup

# subset wb data before program
pre <- subset(mydata,  time.f == "pre", wellbeing,
                 drop = TRUE)

# subset wb data after program
followup <- subset(mydata,  time.f == "followup", wellbeing,
                 drop = TRUE)

pd <- paired(pre, followup)
plot(pd, type = "profile") + theme_bw()

tapply(mydata$depression, mydata$time.f, summary) 
## $pre
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.000   1.839   2.893   2.705   3.429   4.000       5 
## 
## $post
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      NA      NA      NA     NaN      NA      NA      11 
## 
## $followup
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   2.714   3.357   3.679   3.482   3.804   3.857       3
# increase  in average depression from pre to followup

# subset wb data before program
pre <- subset(mydata,  time.f == "pre", depression,
                 drop = TRUE)

# subset wb data after program
followup <- subset(mydata,  time.f == "followup", depression,
                 drop = TRUE)

pd2 <- paired(pre, followup)
plot(pd2, type = "profile") + theme_bw()

corr.df<- mydata[,c("age", "sex", "sears", "bullying", "self.efficacy", "wellbeing", "depression" )]

#calculating correlations and CIs
cor1 <- cor.mtest(corr.df, use="pairwise.complete.obs", conf.level = 0.95)
cor1
## $p
##                     age        sex      sears   bullying self.efficacy
## age           0.0000000 0.52831998 0.53074056 0.66441785    0.73362154
## sex           0.5283200 0.00000000 0.08631929 0.50822224    0.29823093
## sears         0.5307406 0.08631929 0.00000000 0.29716295    0.01150881
## bullying      0.6644178 0.50822224 0.29716295 0.00000000    0.90056102
## self.efficacy 0.7336215 0.29823093 0.01150881 0.90056102    0.00000000
## wellbeing     0.2826385 0.20360322 0.30089400 0.03125455    0.19427076
## depression    0.1045814 0.50656168 0.20300659 0.84011615    0.45716573
##                wellbeing depression
## age           0.28263846 0.10458142
## sex           0.20360322 0.50656168
## sears         0.30089400 0.20300659
## bullying      0.03125455 0.84011615
## self.efficacy 0.19427076 0.45716573
## wellbeing     0.00000000 0.07497734
## depression    0.07497734 0.00000000
## 
## $lowCI
##                       age        sex      sears   bullying self.efficacy
## age            1.00000000 -0.5432856 -0.5820628 -0.5634249    -0.5317130
## sex           -0.54328565  1.0000000 -0.7390552 -0.3361073    -0.2358695
## sears         -0.58206277 -0.7390552  1.0000000 -0.5846447     0.1091642
## bullying      -0.56342490 -0.3361073 -0.5846447  1.0000000    -0.4348452
## self.efficacy -0.53171304 -0.2358695  0.1091642 -0.4348452     1.0000000
## wellbeing     -0.68472737 -0.7128181 -0.2168483 -0.7621742    -0.6445775
## depression    -0.09465741 -0.3473283 -0.6538132 -0.4139148    -0.5740361
##                wellbeing  depression
## age           -0.6847274 -0.09465741
## sex           -0.7128181 -0.34732832
## sears         -0.2168483 -0.65381317
## bullying      -0.7621742 -0.41391476
## self.efficacy -0.6445775 -0.57403610
## wellbeing      1.0000000 -0.71982170
## depression    -0.7198217  1.00000000
## 
## $uppCI
##                     age        sex      sears    bullying self.efficacy
## age           1.0000000 0.30509726 0.33330727  0.38830960     0.3965855
## sex           0.3050973 1.00000000 0.06364069  0.60306320     0.6479316
## sears         0.3333073 0.06364069 1.00000000  0.20412387     0.6844897
## bullying      0.3883096 0.60306320 0.20412387  1.00000000     0.3890382
## self.efficacy 0.3965855 0.64793156 0.68448968  0.38903821     1.0000000
## wellbeing     0.2441964 0.19189383 0.60664115 -0.05061704     0.1566690
## depression    0.7583204 0.61986762 0.16720707  0.49272740     0.2887484
##                 wellbeing depression
## age            0.24419637  0.7583204
## sex            0.19189383  0.6198676
## sears          0.60664115  0.1672071
## bullying      -0.05061704  0.4927274
## self.efficacy  0.15666900  0.2887484
## wellbeing      1.00000000  0.0434200
## depression     0.04342000  1.0000000
cor1b <- cor(corr.df, use="pairwise.complete.obs")

rownames(cor1b) <- c("age",
                                         "female",
                                        "SEARS",
                                         "SCBS",
                                         "SE",
                                         "MHC",
                                         "DASS")


colnames(cor1b) <- c("age",
                                         "female",
                                         "SEARS",
                                         "SCBS",
                                         "SE",
                                         "MHC",
                                         "DASS")


#correlation Matrix
corrplot(cor1b, method = "color",type = "upper", 
                   p.mat = cor1$p,  addCoef.col = 'black', sig.level = 0.05,  
                                    insig = "pch", number.cex=0.85, tl.col="black",tl.cex = 1,
                     col=colorRampPalette(c("hotpink", "white", "#728393"))(50), cl.pos = 'n')

ggviolin(data=subset(mydata1, !is.na(dated.f)), x = "dated.f", y = "sears", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(dated.f)), x = "dated.f", y = "bullying", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(dated.f)), x = "dated.f", y = "self.efficacy", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(dated.f)), x = "dated.f", y = "wellbeing", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(dated.f)), x = "dated.f", y = "depression", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(gender.f)), x = "gender.f", y = "sears", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(gender.f)), x = "gender.f", y = "bullying", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(gender.f)), x = "gender.f", y = "self.efficacy", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(gender.f)), x = "gender.f", y = "wellbeing", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata1, !is.na(gender.f)), x = "gender.f", y = "depression", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata, !is.na(disability.f)), x = "disability.f", y = "sears", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata, !is.na(disability.f)), x = "disability.f", y = "bullying", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata, !is.na(disability.f)), x = "disability.f", y = "self.efficacy", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata, !is.na(disability.f)), x = "disability.f", y = "wellbeing", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()

ggviolin(data=subset(mydata, !is.na(disability.f)), x = "disability.f", y = "depression", fill = "#94AEBC",
         add = "boxplot", add.params = list(fill = "hotpink")) +
  coord_flip()