languageserver
languageserver copied to clipboard
Published
20 hours ago •
REditorSupport
REditorSupport
Lsp fails to start in spacemacs ESS(R) mode
trafficstars
Hi,
I am on spacemacs, so I am not sure whether this is an ESS issue or a spacemacs issue. In any case, I am trying to use the lsp server for R files but it fails to start (even after restarts. The lsp--r::stderr file reads:
[2021-07-05 09:22:04.206] error handling json: Error: lexical error: invalid char in json text.
= dplyr::lag(polls_diff_D_R{"jsonrpc":"2.0","method":"textDo
(right here) ------^
Stack trace:
1: parse_string(txt, bigint_as_char)
2: parseJSON(txt, bigint_as_char)
3: parse_and_simplify(txt = txt, simplifyVector = simplifyVector,
simplifyDataFrame = simplifyDataFrame, simplifyMatrix = simplifyMatrix,
flatten = flatten, ...)
4: jsonlite::fromJSON(data, simplifyVector = FALSE)
[2021-07-05 09:22:04.248] Error: Unexpected non-empty line
Call: self$read_header()
Stack trace:
1: stop("Unexpected non-empty line")
2: self$read_header()
3: self$fetch(blocking = FALSE)
[2021-07-05 09:22:04.248] exiting
Process lsp-r stderr finished
I tried installing the latest version of the language server as suggested here in issue 121).
Any help would greatly be appreciated.
Looks like you language client send a corrupted request to the language server. Does it rarely or frequently occur, or could you reliably reproduce it?
I think I may be getting similar problems in Emacs, however I am using ESS in the doom flavour. See below for an excerpt from a stack trace:
Stack trace:
1: parse_string(txt, bigint_as_char)
2: parseJSON(txt, bigint_as_char)
3: parse_and_simplify(txt = txt, simplifyVector = simplifyVector,
simplifyDataFrame = simplifyDataFrame, simplifyMatrix = simplifyMatrix,
flatten = flatten, ...)
4: jsonlite::fromJSON(data, simplifyVector = FALSE)
[2021-12-02 08:17:51.797] Error: Unexpected non-empty line
Call: self$read_header()
Stack trace:
1: stop("Unexpected non-empty line")
2: self$read_header()
3: self$fetch(blocking = FALSE)
[2021-12-02 08:17:51.798] exiting
Other relevant info:
ProductName: macOS
ProductVersion: 11.6.1
BuildVersion: 20G224
GNU Emacs 27.2
> Executing 'doom version' with Emacs 27.2 at 2021-12-02 08:21:11
Doom emacs v21.12.0-alpha grafted, HEAD -> develop, origin/develop, origin/HEAD 0352ade 2021-11-30 14:05:57 +0100
Doom core v3.0.0-alpha grafted, HEAD -> develop, origin/develop, origin/HEAD 0352ade 2021-11-30 14:05:57 +0100
lsp-mode 20211201.1200
ess 20211122.1708
polymode 20211124.913
poly-R 20210930.1921
poly-markdown 20210625.803
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)
languageserver 0.3.12
jsonlite 1.7.2
This is happening regularly on larger files (e.g., 3000+ lines), both .R and .Rmd.
Could you enable debug mode and post the log? One way to do it is to consider the following option in .Rprofile
options(languageserver.debug = TRUE)
You also need to enable debug log in the lsp plugin for your editor.
Not sure what any of it means (lol), but please see below for an example.
[Trace - 03:06:15 PM] Sending request 'textDocument/codeAction - (374)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
},
"range": {
"start": {
"line": 911,
"character": 0
},
"end": {
"line": 911,
"character": 0
}
},
"context": {
"diagnostics": []
}
}
[Trace - 03:06:15 PM] Received response 'textDocument/codeAction - (374)' in 88ms.
Result: []
[Trace - 03:06:15 PM] Sending request 'textDocument/codeAction - (375)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
},
"range": {
"start": {
"line": 911,
"character": 0
},
"end": {
"line": 911,
"character": 0
}
},
"context": {
"diagnostics": []
}
}
[Trace - 03:06:15 PM] Sending request 'textDocument/documentLink - (376)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
}
}
[Trace - 03:06:15 PM] Received response 'textDocument/codeAction - (375)' in 91ms.
Result: []
[Trace - 03:06:16 PM] Sending notification '$/cancelRequest'.
Params: {
"id": 376
}
[Trace - 03:06:18 PM] Sending request 'textDocument/codeAction - (377)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
},
"range": {
"start": {
"line": 911,
"character": 0
},
"end": {
"line": 911,
"character": 0
}
},
"context": {
"diagnostics": []
}
}
[Trace - 03:06:18 PM] Received response 'textDocument/codeAction - (377)' in 176ms.
Result: []
[Trace - 03:06:24 PM] Sending request 'textDocument/codeAction - (378)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
},
"range": {
"start": {
"line": 911,
"character": 0
},
"end": {
"line": 911,
"character": 0
}
},
"context": {
"diagnostics": []
}
}
[Trace - 03:06:24 PM] Received response 'textDocument/codeAction - (378)' in 90ms.
Result: []
[Trace - 03:06:24 PM] Sending request 'textDocument/codeAction - (379)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
},
"range": {
"start": {
"line": 911,
"character": 0
},
"end": {
"line": 911,
"character": 0
}
},
"context": {
"diagnostics": []
}
}
[Trace - 03:06:24 PM] Sending request 'textDocument/documentLink - (380)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
}
}
[Trace - 03:06:25 PM] Sending notification 'textDocument/didClose'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
}
}
[Trace - 03:06:25 PM] Sending notification 'textDocument/didOpen'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R",
"languageId": "r",
"version": 0,
"text": "### ============================================================================\n### SOURCE SCRIPTS\n### ============================================================================\nrm(list = ls(all = TRUE))\nsource(\"./code/func.R\")\nsource(\"./code/load.R\")\n\nload(\"./data/proc/draft.RData\")\nls(all = TRUE)\n\n## save.image(\"./data/tmp/draft.RData\")\n\n### ============================================================================\n### LOAD DATA\n### ============================================================================\n\n\n## mplus missing data analysis\nlpa_miss <-\n readRDS(here(.m1_optseed))$c0_optseed_m_v1.out\n\n## full lpa output\nll_tests_m <-\n readRDS(here(.m2_ll_tests))\n\n## first split half lpa output\nll_tests_m_rv <-\n readRDS(here(.m2_ll_tests_rv))\n\n## three class lpa output\nlpa_c3 <-\n ll_tests_m$c3_2_ll_tests_m_v1.out\n\n## final data set with three class lpa profiles added\nlpa_fin <-\n readRDS(here(.lpa_fin))\n\nlpa_fin %>%\n select(session, final_session, pb, tb, bhs, ac, dep, anx, c, si) %>%\n filter(session == 1 | final_session == 1) %>%\n group_by(final_session) %>%\n .summary_by(c)\n\nlpa_fin %>%\n select(ID, c, session, final_session) %>%\n filter(session == 1) %>%\n describeBy(\"c\")\n\nlpa_fin %>%\n select(ID, c, session, final_session) %>%\n filter(session != 1 & final_session == 1) %>%\n describeBy(\"c\")\n\nlpa_fin %>%\n group_by(c) %>%\n summarise(\n n_final = length(final_session[which(final_session == 1 & session != 1)]),\n n_init = length(session[session %in% 1]),\n ) %>%\n filter(!is.na(c)) %>%\n mutate(\n perc_lost = round((n_init - n_final) / n_init * 100, digits = 1),\n c = as.character(c)\n ) %>%\n janitor::adorn_totals()\n\n## CREATE MI DATASETS\n## =============================================================================\n## SUBSET TIME 1: DEMO\n## -----------------------------------------------------------------------------\nlpa_t0 <-\n lpa_fin %>%\n filter(!is.na(c), session == 1 | final_session == 1) %>%\n group_by(ID) %>%\n mutate(\n total_session = max(session), .after = session,\n planned_final = max(planned_final)\n ) %>%\n ungroup() %>%\n filter(session == 1) %>%\n select(ID:suburb)\n\nlpa_t0 %>% names()\nlpa_t0 %>% summary()\n\n## SUBSET TIME 1: MI\n## -----------------------------------------------------------------------------\nlpa_t0_mi <-\n lpa_fin %>%\n filter(session == 1 & !is.na(c)) %>%\n mutate(across(c(c:sex, atsi:closehxdsh, suicide_attempts), as.factor)) %>%\n select(c, ID, pb, ac, bhs, tb, age, sex, dep, anx, si)\n\nlpa_t0_mi %>% names()\nlpa_t0_mi %>% .summary_by(c)\n\n## SUBSET TIME 1 WITH MATCHING TIME 2\n## -----------------------------------------------------------------------------\nlpa_t0t1 <-\n lpa_fin %>%\n ## filter no class var\n filter(!is.na(c)) %>%\n group_by(ID) %>%\n ## filter first and last session\n filter(session == 1 | final_session == 1) %>%\n ungroup() %>%\n mutate(\n ## replace NA session value\n session = if_else(is.na(session), -999, session),\n ## create conditional variable to drop those with only one session\n drop = if_else(\n session == 1 & final_session == 1,\n 1, 0\n ), .after = ID\n ) %>%\n ## select those with more than one session\n filter(drop == 0) %>%\n ## select vars for mi\n select(\n c, ID, age, sex, session, final_session,\n planned_final, pb, ac, bhs, tb, dep, anx, si\n )\n\nlpa_t0t1_mi <-\n lpa_t0t1 %>%\n panelr::panel_data(\n id = ID,\n wave = final_session,\n ) %>%\n panelr::widen_panel(\n separator = \"_\",\n ignore.attributes = FALSE,\n varying = NULL) %>%\n ungroup() %>%\n ## need to convert to data frame as panelr creates a list object\n as.data.frame() %>%\n ## use age at intake, drop secon ac as not routinely collected\n select(-age_1, -ac_1) %>%\n rename(age = age_0) %>%\n mutate(\n ## remove -999 placehoder session number\n session_1 = if_else(session_1 == -999, NA_real_, session_1)\n )\n\nlpa_t0t1_mi %>% str()\nlpa_t0t1_mi %>% summary()\n\n### ============================================================================\n### MICE MI\n### ============================================================================\n## TODO explore missingness overall\n## TODO explore missingness by profile\n## NOTE need to as a minimum report missingness percentages\n## Nice to have would be to determine if MAR or MNAR instead of just assuming\n## TODO subset vars\n## DONE get time 1 and final session data and bring into imputation model\n## DONE decide if can also bring in the demographic variables to be imputed\n## NOTE Frank suggested this was unusual and so was not done\n## https://www.gerkovink.com/miceVignettes/Multi_level/Multi_level_data.html\n\n## TIME 1\n## =============================================================================\n.summary_by(lpa_t0_mi, c)\n\nmd.pattern(lpa_t0_mi)\n\n.miss_df(lpa_t0_mi)\n\nfluxplot(lpa_t0_mi)\n## NOTE conisder historgrams inspecting missingness here too\n\n## GET AND SET PREDICTOR VARS\n## -----------------------------------------------------------------------------\n## NOTE this was used to create an excel for inspection and manipulation\npred_test <-\n quickpred(lpa_t0_mi) %>%\n .mice_pred_df()\n\n.pred_ <-\n quickpred(lpa_t0_mi)\n.pred_\n\n## import pred via excel manipulation\npred_ <-\n readxl::read_xlsx(\n here(\"data/tmp/mi_working.xlsx\"),\n sheet = \"pred_matrix\",\n range = \"A1:L12\"\n ) %>%\n .mice_pred_matrix()\npred_\n\n## comfirm they are different\nidentical(pred_, .pred_)\n\n## create pred here without exporting to excel\npred <- quickpred(lpa_t0_mi)\n## c\n## pred[,\"c\"] <- c(0, 0, seq(-2, -2, length.out = 9))\npred[, \"c\"] <- c(0, 0, seq(1, 1, length.out = 9))\n## pred[,\"c\"] <-\n## ifelse(pred[,\"c\"] == 0, 0, -2)\n## pred[\"age\", \"c\"] <- 0\npred[, \"ID\"] <- ifelse(pred[,\"ID\"] == 0, 0, 0)\n## pred[,c(\"ID\", \"ac\", \"anx\", \"dep\", \"si\", \"sex\", \"age\")] <-\n## ifelse(pred[,c(\"ID\", \"ac\", \"anx\", \"dep\", \"si\", \"sex\", \"age\")] == 0, 0, 0)\npred\n\n## comfirm they are different\nidentical(pred, pred_)\n\n## IMPUTE DATA\n## =============================================================================\nlpa_t0_mi <-\n lpa_t0_mi %>%\n mutate(c = as.integer(c))\n\n## init mi model\nimp_t0 <-\n mice(\n lpa_t0_mi,\n pred = pred,\n meth = \"pmm\",\n maxit = 0,\n m = 10\n )\n\nimp_t0_5 <-\n mice.mids(\n imp_t0,\n maxit = 5\n )\n\nimp_t0_15 <-\n mice.mids(\n imp_t0_5,\n maxit = 10\n )\n\nimp_t0_25 <-\n mice.mids(\n imp_t0_15,\n maxit = 10\n )\n\nimp_t0_35 <-\n mice.mids(\n imp_t0_25,\n maxit = 10\n )\n\nimp_t0_50 <-\n mice.mids(\n imp_t0_35,\n maxit = 15\n )\n\nimp_t0_65 <-\n mice.mids(\n imp_t0_50,\n maxit = 15\n )\n\nimp_t0_80 <-\n mice.mids(\n imp_t0_65,\n maxit = 15\n )\n\nimp_t0_100 <-\n mice.mids(\n imp_t0_80,\n maxit = 20\n )\n\nimp_t0_115 <-\n mice.mids(\n imp_t0_100,\n maxit = 15\n )\n\nimp_t0_150 <-\n mice.mids(\n imp_t0_115,\n maxit = 35\n )\n\nimp_t0_200 <-\n mice.mids(\n imp_t0_150,\n maxit = 50\n )\n\nimp_t0_300 <-\n mice.mids(\n imp_t0_200,\n maxit = 100\n )\n\nimp_t0_400 <-\n mice.mids(\n imp_t0_300,\n maxit = 100\n )\n\n## inspect each iterations convergence\nsummary(lpa_t0_mi)\nsummary(complete(imp_t0_5))\nsummary(complete(imp_t0_15))\nsummary(complete(imp_t0_25))\nsummary(complete(imp_t0_35))\nsummary(complete(imp_t0_50))\nsummary(complete(imp_t0_65))\nsummary(complete(imp_t0_80))\n\nsummary(complete(imp_t0_100))\nsummary(complete(imp_t0_115))\nsummary(complete(imp_t0_150))\nsummary(complete(imp_t0_200))\nsummary(complete(imp_t0_300))\nsummary(complete(imp_t0_400))\n.summary_by(complete(imp_t0_400), c)\n\nplot(imp_t0_5)\nplot(imp_t0_15)\nplot(imp_t0_25)\nplot(imp_t0_35)\nplot(imp_t0_50)\nplot(imp_t0_65) # this looks okay, perhaps use this one\nplot(imp_t0_80) # this seems to be getting worse\nplot(imp_t0_100) # not too bad\nplot(imp_t0_115) # okay\nplot(imp_t0_150)\nplot(imp_t0_200)\nplot(imp_t0_300)\nplot(imp_t0_400)\n\ndensityplot(imp_t0_5)\ndensityplot(imp_t0_15)\ndensityplot(imp_t0_25)\ndensityplot(imp_t0_35)\ndensityplot(imp_t0_50)\ndensityplot(imp_t0_65)\ndensityplot(imp_t0_80)\ndensityplot(imp_t0_100)\ndensityplot(imp_t0_115)\ndensityplot(imp_t0_150)\ndensityplot(imp_t0_200)\ndensityplot(imp_t0_300)\ndensityplot(imp_t0_400)\n\n\n## NOTE no perfect solution, so will chose to use the imp_t0_400 mids as this\n## has the most iterations and is assumed to be the best as a result.\n\n## TIME 1 AND TIME 2\n## =============================================================================\n## GET AND SET PREDICTOR VARS\n## -----------------------------------------------------------------------------\n## NOTE this was used to create an excel for inspection and manipulation\n## NOTE this was guided by the longitudinal imputation based on the firewords\n## disater in the mice textbook\npred_test_fu <-\n quickpred(lpa_t0t1_mi) %>%\n .mice_pred_df()\npred_test_fu\n\n.pred_fu_ <-\n quickpred(lpa_t0t1_mi)\n.pred_fu_\n\n## import pred via excel manipulation\npred_fu_ <-\n readxl::read_xlsx(\n here(\"data/tmp/mi_working.xlsx\"),\n sheet = \"pred_matrix_fu_simp\",\n range = \"A1:V22\"\n ) %>%\n .mice_pred_matrix()\npred_fu_\n\n## comfirm they are different\nidentical(pred_fu_, .pred_)\n\n## IMPUTE DATA\n## =============================================================================\nlpa_t0t1_mi <-\n lpa_t0t1_mi %>%\n mutate(c = as.integer(c))\n\n## init mi model\nimp_t0t1 <-\n mice(\n lpa_t0t1_mi,\n pred = pred_fu_,\n meth = \"pmm\",\n maxit = 0,\n m = 10\n )\n\nimp_t0t1_5 <-\n mice.mids(\n imp_t0t1,\n maxit = 5\n )\n\nimp_t0t1_15 <-\n mice.mids(\n imp_t0t1_5,\n maxit = 10\n )\n\nimp_t0t1_25 <-\n mice.mids(\n imp_t0t1_15,\n maxit = 10\n )\n\nimp_t0t1_35 <-\n mice.mids(\n imp_t0t1_25,\n maxit = 10\n )\n\nimp_t0t1_50 <-\n mice.mids(\n imp_t0t1_35,\n maxit = 15\n )\n\nimp_t0t1_100 <-\n mice.mids(\n imp_t0t1_50,\n maxit = 50\n )\n\nimp_t0t1_115 <-\n mice.mids(\n imp_t0t1_100,\n maxit = 15\n )\n\n## inspect each iterations convergence\n.summary_by(lpa_t0t1_mi, c)\n.summary_by(complete(imp_t0t1_5), c)\n.summary_by(complete(imp_t0t1_15), c)\n.summary_by(complete(imp_t0t1_25), c)\n.summary_by(complete(imp_t0t1_35), c)\n.summary_by(complete(imp_t0t1_50), c)\n.summary_by(complete(imp_t0t1_100), c)\n.summary_by(complete(imp_t0t1_115), c)\n\n\nplot(imp_t0t1_5)\nplot(imp_t0t1_15)\nplot(imp_t0t1_25)\nplot(imp_t0t1_35) # potential candidate\nplot(imp_t0t1_50) # pretty good with exception to tb_0\nplot(imp_t0t1_100) # not much of an improvement on the above\nplot(imp_t0t1_115) # convergence warnings, will use 50\n\ndensityplot(imp_t0t1_5)\ndensityplot(imp_t0t1_15)\ndensityplot(imp_t0t1_25)\ndensityplot(imp_t0t1_35)\ndensityplot(imp_t0t1_50)\ndensityplot(imp_t0t1_100)\ndensityplot(imp_t0t1_115)\n\n### ============================================================================\n### ANALYSESES\n### ============================================================================\n\n## LPA CLASS TABL\n## =============================================================================\n\nlpa_class_tab <-\n SummaryTable(\n ll_tests_m,\n keepCols = c(\n \"Title\",\n \"LL\",\n \"Entropy\",\n \"BIC\",\n \"T11_KM1LL\",\n \"T11_VLMR_PValue\"\n ),\n sortBy = \"Title\"\n ) %>%\n mutate(\n Classes = str_c(str_extract(Title, \"\\\\d\"), \" Classes\"), .before = 1,\n across(where(is.numeric) & !ends_with(\"Value\"), ~ round(.x, 2)),\n across(ends_with(\"Value\"), ~ round(.x, 3))\n ) %>%\n select(-Title) %>%\n filter(Classes != \"1 Classes\")\nlpa_class_tab\n\n\n## RELIABILITY\n## =============================================================================\n\n.pb_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(inq_1:inq_6) %>%\n psych::alpha()\nround(summary(.pb_alpha)$raw_alpha, 2)\n\n.tb_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(inq_7r:inq_15r, inq_9, inq_11:inq_12) %>%\n psych::alpha()\nsummary(.tb_alpha)\n\n.bhs_20 <-\n lpa_fin %>%\n filter(!is.na(c), session == 1, protocol == 1) %>%\n ## names() %>%\n select(starts_with(\"bhs_\") & ends_with(\"r\")) %>%\n psych::alpha()\nsummary(.bhs_20)\n\n.bhs_sf <-\n lpa_fin %>%\n filter(!is.na(c), session == 1, protocol == 2) %>%\n select(starts_with(\"bhsSF_\") & ends_with(\"lessone\")) %>%\n psych::alpha()\nsummary(.bhs_sf)\n\n.si_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(mssi_1:mssi_4, starts_with(\"rc_mssi\")) %>%\n psych::alpha()\nsummary(.si_alpha)\n\n.ac_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(c(acss_8r:acss_13r, acss_7, acss_11, acss_14, acss_19)) %>%\n psych::alpha()\nsummary(.ac_alpha)\n\n.dep_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"dass_\", c(3, 5, 10, 13, 16, 17, 21))) %>%\n psych::alpha()\nsummary(.dep_alpha)\n\n.anx_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"dass_\", c(2, 4, 7, 9, 15, 19, 20))) %>%\n psych::alpha()\nsummary(.anx_alpha)\n\n.str_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"dass_\", c(1, 6, 8, 11, 12, 14, 18))) %>%\n psych::alpha()\nsummary(.str_alpha)\n\n.is_aff_reg_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(1, 14, 27))) %>%\n psych::alpha()\nsummary(.is_aff_reg_alpha)\n\n.is_int_bnd_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(2, 15, 28))) %>%\n psych::alpha()\nsummary(.is_int_bnd_alpha)\n\n.is_slf_pun_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(3, 16, 29))) %>%\n psych::alpha()\nsummary(.is_slf_pun_alpha)\n\n.is_slf_car_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(4, 17, 30))) %>%\n psych::alpha()\nsummary(.is_slf_car_alpha)\n\n.is_ant_dis_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(5, 18, 31))) %>%\n psych::alpha()\nsummary(.is_ant_dis_alpha)\n\n.is_ant_sui_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(6, 19, 32))) %>%\n psych::alpha()\nsummary(.is_ant_sui_alpha)\n\n.is_sen_see_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(7, 20, 33))) %>%\n psych::alpha()\nsummary(.is_sen_see_alpha)\n\n.is_peer_bnd_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(8, 21, 34))) %>%\n psych::alpha()\nsummary(.is_peer_bnd_alpha)\n\n.is_int_inf_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(9, 22, 35))) %>%\n psych::alpha()\nsummary(.is_int_inf_alpha)\n\n.is_tough_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(10, 23, 36))) %>%\n psych::alpha()\nsummary(.is_tough_alpha)\n\n.is_mrk_dis_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(11, 24, 37))) %>%\n psych::alpha()\nsummary(.is_mrk_dis_alpha)\n\n.is_revenge_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(12, 25, 38))) %>%\n psych::alpha()\nsummary(.is_revenge_alpha)\n\n.is_autnmy_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(13, 36, 39))) %>%\n psych::alpha()\nsummary(.is_autnmy_alpha)\n\n.isas_intra_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(\n 1, 3, 5, 6, 11, 14, 16, 18,\n 19, 24, 27, 29, 31, 32, 37\n ))) %>%\n psych::alpha()\nsummary(.isas_intra_alpha)\n\n.isas_social_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(\n 2, 4, 7, 8, 9, 10, 12, 13, 15, 17,\n 20, 21, 22, 23, 25, 26, 28, 30,\n 33, 34, 35, 36, 38, 39\n ))) %>%\n psych::alpha()\nsummary(.isas_social_alpha)\n\nmeas_rel <-\n list(\n pb_rel = .pb_alpha,\n tb_rel = .tb_alpha,\n bhs_20_rel = .bhs_20,\n bhs_sf_rel = .bhs_sf,\n si_rel = .si_alpha,\n ac_rel = .ac_alpha,\n dep_rel = .dep_alpha,\n anx_rel = .anx_alpha,\n str_rel = .str_alpha,\n is_aff_reg_rel = .is_aff_reg_alpha,\n is_int_bnd_rel = .is_int_bnd_alpha,\n is_slf_pun_rel = .is_slf_pun_alpha,\n is_slf_car_rel = .is_slf_car_alpha,\n is_ant_dis_rel = .is_ant_dis_alpha,\n is_ant_sui_rel = .is_ant_sui_alpha,\n is_sen_see_rel = .is_sen_see_alpha,\n is_peer_bnd_rel = .is_peer_bnd_alpha,\n is_int_inf_rel = .is_int_inf_alpha,\n is_tough_rel = .is_tough_alpha,\n is_mrk_dis_rel = .is_mrk_dis_alpha,\n is_revenge_rel = .is_revenge_alpha,\n is_autnmy_rel = .is_autnmy_alpha,\n isas_intra_rel = .isas_intra_alpha,\n isas_social_rel = .isas_social_alpha\n )\n\n## NOTE e.g., do not run\n## round(summary(meas_rel$pb_rel)$raw_alpha, 2)\n\n## CROSS SECTIONAL: NON-MIDS\n## =============================================================================\n## Demographics\n## ANOVA on c predicted by covariates\n## Pairwise t-tests of covariates and LPA predcitors\n## NOTE optional analysis\n## multilevel regression si ~ all the other vars\n## Multinomial logistic regression instead of AVNOA, this could be interesting\n## to note units of change and therefor more informative from a prediction\n## perspective.\n\n.lpa_ch <-\n complete(imp_t0t1_50, \"long\") %>%\n group_by(.imp) %>%\n mutate(\n ## calculate change scores\n pb_ch = pb_1 - pb_0,\n tb_ch = tb_1 - tb_0,\n bhs_ch = bhs_1 - bhs_0,\n pb_ch = pb_1 - pb_0,\n dep_ch = dep_1 - dep_0,\n anx_ch = anx_1 - anx_0,\n ## str_ch = str_1 - str_0,\n si_ch = si_1 - si_0,\n ) %>%\n ungroup()\n\nAnova(lm(pb_ch + tb_ch + bhs_ch + dep_ch + anx_ch + si_ch ~ c, .lpa_ch), type = 3) %>%\n tidy()\n\n.lpa_ch %>%\n select(ends_with(\"_ch\")) %>%\n ## select(ends_with(\"_1\")) %>%\n map(~ tidy(pairwise.t.test(., .lpa_ch$c, p.adjust.method = \"bonf\")))\n\n.fu_look <-\n .lpa_ch %>%\n mutate(across(everything(), as.numeric)) %>%\n group_by(c) %>%\n summarise(\n across(ends_with(c(\"_0\", \"_1\", \"_ch\")) & !contains(\"session\"),\n list(m = mean, sd = sd),\n .names = \"{.col}_{.fn}\")\n ) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"_0|_1|_ch|_m|_sd\", \"\"),\n time = str_extract(key, \"0|1|ch\"),\n type = str_extract(key, \"m|sd\"),\n ) %>%\n select(-key) %>%\n pivot_wider(\n names_from = c(c, type),\n values_from = value\n ) %>%\n clean_names()\n\n\n.lpa_ch %>%\n mutate(across(everything(), as.numeric)) %>%\n group_by(c) %>%\n t.test(pb_0 ~ pb_1, paried = TRUE)\n\nAnova(lm(pb_ch ~ c, data = .lpa_ch), type = 3) %>%\n\n## ANOVA BY C\n## -----------------------------------------------------------------------------\n## TODO get the standardised residuals of each model\n## https://www.statology.org/standardized-residuals-in-r/\n## get anova stats\n.lpa_anova <-\n Anova(lm(c ~ pb + tb + bhs + ac, data = lpa_t0), type = 3) %>%\n tidy() %>%\n select(term, statistic, p.value) %>%\n filter(term %in% c(\"pb\", \"tb\", \"bhs\", \"ac\")) %>%\n rename(var = term)\n\n.lpa_anova_fx <-\n Anova(lm(c ~ pb + tb + bhs + ac, data = lpa_t0), type = 3) %>%\n effectsize::omega_squared() %>%\n as.data.frame() %>%\n mutate(\n fx_size = interpret_omega_squared(Omega2_partial, rules = \"cohen1992\"),\n omega2 = str_remove(sprintf(\"%.2f\", Omega2_partial), \"^0\")\n ) %>%\n rename(var = Parameter) %>%\n select(var, omega2, fx_size)\n\nlpa_sample <-\n lpa_t0 %>%\n select(c, pb, tb, bhs, ac) %>%\n mutate(\n across(c(pb, tb, bhs, ac),\n ~ scale(.x),\n .names = \"{.col}_z\"\n )\n ) %>%\n group_by(c) %>%\n summarise(\n N = n(),\n across(\n c(pb, tb, bhs, ac),\n list(\n m = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE)\n )\n ),\n across(\n ends_with(\"_z\"),\n ~ round(mean(.x, na.rm = TRUE), 2)\n )\n ) %>%\n mutate(\n prc = N / sum(N) * 100, .after = N\n ) %>%\n select(c, N, prc, ends_with(c(\"_m\", \"_sd\", \"_z\")))\n\n## get mean and sd\n.lpa_m_sd <-\n ## NOTE cannot find original object `lpa_sample`\n lpa_sample %>%\n select(c, ends_with(c(\"_m\", \"_sd\"))) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"_m|_sd\", \"\"),\n type = str_extract(key, \"m|sd\"),\n ) %>%\n select(-key) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value\n )\n\n## make contrasts wide\n.lpa_t <-\n lpa_t0 %>%\n select(pb, tb, bhs, ac) %>%\n map( ~ tidy(pairwise.t.test(., lpa_t0$c, p.adjust.method = \"bonf\"))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n mutate(across(where(is.character), as.numeric)) %>%\n pivot_longer(\n ## col = ends_with(c(\"group1\", \"group2\", \"value\")),\n col = contains(\"value\"),\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n select(contains(c(\"group1\", \"group2\")), everything()) %>%\n arrange(key) %>%\n ## NOTE this was taken from SO\n ## TODO find tidy way to do this\n .[!duplicated(lapply(., summary))] %>%\n rename(group1 = 1, group2 = 2, var = 3, p_value = 4) %>%\n mutate(\n var = str_extract(var, \"ac|bhs|pb|tb\")\n ) %>%\n pivot_wider(\n names_from = c(group1, group2),\n values_from = p_value\n ) %>%\n janitor::clean_names()\n\n## join the lpa var tests\n.lpa_anova_full <-\n .lpa_anova %>%\n full_join(.lpa_m_sd) %>%\n full_join(.lpa_t) %>%\n full_join(.lpa_anova_fx)\n\n## generate the table for the report\nlpa_anova <-\n .lpa_anova_full %>%\n mutate(\n group_diff = case_when(\n x2_1 < 0.05 & x3_1 < 0.05 & x3_2 < 0.05 ~ \"1 < 2,3; 2 < 3\",\n TRUE ~ \"null\"\n ),\n f = case_when(\n ## TODO check that markdown formatted asterix\n p.value < 0.001 ~ paste0(as.character(round(statistic, 2)), \"<sup>***</sup>\"),\n p.value < 0.05 ~ paste0(as.character(round(statistic, 2)), \"<sup>*</sup>\")\n ),\n across(starts_with(c(\"m_\", \"sd_\")), ~ sprintf(\"%.2f\", round(.x, 2))),\n low_si = str_c(m_1, \" (\", sd_1, \")\"),\n mod_si = str_c(m_2, \" (\", sd_2, \")\"),\n hi_si = str_c(m_3, \" (\", sd_3, \")\"),\n var = case_when(\n var == \"pb\" ~ \"Burdensomeness\",\n var == \"tb\" ~ \"Belongingness\",\n var == \"ac\" ~ \"Capability\",\n var == \"bhs\" ~ \"Hopelessness\"\n )\n ) %>%\n select(var, low_si:hi_si, f, omega2, group_diff)\n\n## FOLLOW UP ANALYSIS BY C\n## -----------------------------------------------------------------------------\n## recode values of demo vars\n.fu <-\n ## use full data set here to look at isas items\n lpa_fin %>%\n filter(!is.na(c), session %in% 1 | final_session %in% 1) %>%\n group_by(ID) %>%\n mutate(\n total_session = max(session), .after = session,\n planned_final = max(planned_final)\n ) %>%\n ungroup() %>%\n filter(session %in% 1) %>%\n select(\n ID:suburb, isas_1a:isas_1m, isas_3a,\n isas_3b, isas_4:isas_7, siss\n ) %>%\n mutate(\n siss_n = as.numeric(as.character(siss)), .after = everything(),\n isas_days = floor(\n interval(isas_3b, date) / duration(n = 1, unit = \"days\")),\n ## .before = date,\n ## remove strings from isas vars and make all numeric class\n across(c(isas_1i, isas_1m), ~ str_extract(.x, \"\\\\d*\")),\n across(starts_with(\"isas_1\"), as.numeric),\n ## prorate isas higher-order factors for comparison\n isas_intra = isas_intra / 5,\n isas_social = isas_social / 8\n ) %>%\n relocate(isas_days, .before = date) %>%\n rowwise() %>%\n mutate(\n ## count the number of methods of dsh\n dsh_n = sum(c_across(starts_with(\"isas_1\")) > 0, na.rm = TRUE)\n ) %>%\n ungroup() %>%\n mutate(\n across(isas_4:isas_7, as.factor),\n sex = case_when(\n sex == 0 ~ \"Female\",\n sex == 1 ~ \"Male\",\n sex == 2 ~ \"Other\",\n TRUE ~ NA_character_\n ),\n atsi = case_when(\n ## correct atsi status from old sp_db\n sp_db == \"old\" & atsi == 4 ~ \"Neither\",\n ## collapse categories as cell sizes too small\n atsi == 0 ~ \"Neither\",\n atsi %in% c(1, 2, 3) ~ \"ATSI\",\n TRUE ~ NA_character_\n ),\n living = case_when(\n living == 1 ~ \"family\",\n living == 2 ~ \"other\",\n living == 3 ~ \"other\",\n living == 4 ~ \"other\",\n TRUE ~ NA_character_\n ),\n relation = case_when(\n relation == 1 ~ \"single/separated/divorced\",\n relation == 2 ~ \"relationship/married/defacto\",\n relation == 3 ~ \"relationship/married/defacto\",\n relation == 4 ~ \"single/separated/divorced\",\n TRUE ~ NA_character_\n ),\n income = case_when(\n income == 1 ~ \"working\",\n income == 2 ~ \"working\",\n income == 3 ~ \"not-working\",\n income == 4 ~ \"not-working\",\n income == 5 ~ \"not-working\",\n income == 6 ~ \"student\",\n TRUE ~ NA_character_\n ),\n across(starts_with(\"diag\"), ~ case_when(\n .x == 1 ~ \"no diagnosis warranted\",\n .x == 2 ~ \"depression\",\n .x == 3 ~ \"generalised anxiety\",\n .x == 4 ~ \"phobia\",\n .x == 5 ~ \"adjustment disorder\",\n .x == 6 ~ \"acute stress disorder\",\n .x == 7 ~ \"ptsd\",\n .x == 8 ~ \"obsessive-compulsive\",\n .x == 9 ~ \"eating disorder\",\n .x == 13 ~ \"bipolar\",\n .x == 15 ~ \"personality disorder (e.g. BPD)\",\n .x == 17 ~ \"neurological disorder\",\n .x == 18 ~ \"other\",\n TRUE ~ NA_character_\n )),\n across(\n c(\n c, planned_final, sp_db, location, sex, atsi:diag1,\n starts_with(c(\"past\", \"hx\", \"fxhx\")), suicide_attempts,\n isas_4:isas_7, siss\n ) & -siss_n,\n as.factor),\n )\n.fu\n\n## CONTINUOUS VARS\n## -------------------------------\n## class sizes full only\n.lpa_n_fu_full <-\n .fu %>%\n select(ID) %>%\n summarise(\n n_sd = sum(!is.na(ID))\n ) %>%\n mutate(\n n_m = round(n_sd / sum(n_sd) * 100, 2),\n c = \"full\"\n )\n\n## class sizes all\n.lpa_n_fu <-\n .fu %>%\n select(c, ID) %>%\n group_by(c) %>%\n summarise(\n n_sd = sum(!is.na(ID)),\n ) %>%\n mutate(\n n_m = round(n_sd / sum(n_sd) * 100, 2),\n c = as.character(c)\n ) %>%\n bind_rows(.lpa_n_fu_full) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"_m|_sd\", \"\"),\n type = str_extract(key, \"(m|sd)$\")\n ) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value,\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup() %>%\n mutate(var = str_replace(var, \"n\", \"Class size\")) %>%\n select(var, starts_with(c(\"m\", \"sd\"))) %>%\n relocate(ends_with(\"_full\"), .after = var)\n.lpa_n_fu\n\n## univariate anovas of covariates\n## PARAMETRIC\n.fu_test <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n map(~ tidy(Anova(lm(. ~ .fu$c, na.action = na.exclude), type = 3))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n filter(across(ends_with(\"term\"), ~ .x == \".fu$c\")) %>%\n select(ends_with(c(\"statistic\", \"p.value\"))) %>%\n pivot_longer(\n everything(),\n names_to = c(\"var\", \"stat\"),\n names_sep = \"\\\\.\"\n ) %>%\n pivot_wider(\n names_from = stat,\n values_from = value\n )\n.fu_test %>%\n print(n = Inf)\n\n## individual contrasts of fu covariate data\n.fu_test_cntrst <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n map(~ tidy(pairwise.t.test(., .fu$c, p.adjust.method = \"bonf\"))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n mutate(across(where(is.character), as.numeric)) %>%\n pivot_longer(\n col = contains(\"value\"),\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n select(contains(c(\"group1\", \"group2\")), everything()) %>%\n arrange(key) %>%\n ## TODO find tidy way to do this\n .[!duplicated(lapply(., summary))] %>%\n rename(group1 = 1, group2 = 2, var = 3, p_value = 4) %>%\n mutate(\n var = str_replace(var, \"\\\\.p\\\\.value\", \"\")\n ) %>%\n pivot_wider(\n names_from = c(group1, group2),\n values_from = p_value\n ) %>%\n janitor::clean_names()\n\n.fu_test_cntrst %>%\n print(n = Inf)\n\n## univariate anovas of covariates\n## NON-PARAMETRIC: KRUSKAL WALLIS\n## NOTE used to check that that the results for the isas section 1 is the same\n## NOT used for anlalysis otherwise\n.fu_test_kw <-\n .fu %>%\n ## .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n map(~ kruskal.test(. ~ .fu$c, na.action = na.exclude), type = 3) %>%\n map(`[`, c(\"statistic\", \"p.value\")) %>%\n unlist(recursive = FALSE) %>%\n map(~ as_tibble(.)) %>%\n map(~ janitor::clean_names(.)) %>%\n bind_rows(., .id = \"var\") %>%\n rename(key = var) %>%\n mutate(\n var = str_replace_all(key, \"\\\\.statistic|\\\\.p\\\\.value\", \"\"),\n type = str_extract(key, \"statistic|p\\\\.value\"),\n ) %>%\n pivot_wider(\n names_from = type,\n values_from = value\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup() %>%\n rename(\n statistic_kw = statistic,\n p_kw = p.value\n )\n\n.fu_test_kw %>%\n print(n = Inf)\n\n## individual contrasts of fu covariate data\n.fu_test_kw_cntrst <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n ## map(~ pgirmess::kruskalmc(. ~ .fu$c)$dif.com) %>%\n map(~ agricolae::kruskal(., .fu$c, p.adj = \"bonf\", group = FALSE)$comparison) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n comp = c(\"1v2\", \"1v3\", \"2v3\"),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n as_tibble() %>%\n select(var, comp, everything()) %>%\n arrange(var, comp)\n\n.fu_test_kw_cntrst %>%\n print(n = Inf)\n\n## join the the anova and kw tests to eyeball for consitency\n.fu_test_cont_full <-\n .lpa_m_sd_fu %>%\n full_join(.fu_test_kw ) %>%\n full_join(.fu_test_kw_cntrst) %>%\n full_join(.fu_test) %>%\n full_join(.fu_test_cntrst)\n## NOTE eyebal check finishes here\n\n\n## CONTINOUS MEANS\n.lpa_m_sd_fu_full <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n mutate(across(everything(), as.numeric)) %>%\n summarise(\n across(\n everything(),\n list(\n m = ~ round(mean(.x, na.rm = TRUE), 2),\n sd = ~ round(sd(.x, na.rm = TRUE), 2)\n )\n )\n ) %>%\n mutate(\n c = \"full\", .before = 1\n ) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"(?!_mrk_)(_m|_sd)\", \"\"),\n type = str_extract(key, \"(m|sd)$\"),\n ) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value,\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup()\n\n.lpa_m_sd_fu_full %>%\n print(n = Inf)\n\n## means and sd by class\n.lpa_m_sd_fu <-\n .fu %>%\n select(\n c, si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n mutate(across(everything(), as.numeric)) %>%\n group_by(c) %>%\n summarise(\n across(\n everything(),\n list(\n m = ~ round(mean(.x, na.rm = TRUE), 2),\n sd = ~ round(sd(.x, na.rm = TRUE), 2)\n )\n )\n ) %>%\n ungroup() %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"(?!_mrk_)(_m|_sd)\", \"\"),\n type = str_extract(key, \"(m|sd)$\"),\n ) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value,\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup() %>%\n full_join(.lpa_m_sd_fu_full) %>%\n select(var, ends_with(\"full\"), everything())\n\n.lpa_m_sd_fu %>%\n print(n = Inf)\n\n## FULL CONTINUOUS FU TABLE\nfu_test_cont <-\n .lpa_n_fu %>%\n full_join(.lpa_m_sd_fu) %>%\n full_join(.fu_test) %>%\n full_join(.fu_test_cntrst)\n\nfu_test_cont %>%\n print(n = Inf)\n\n## CATEGORICAL FU VARS\n## --------------------------\n## chisq on fu covariate data\n.fu_test_nonp <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex)\n ) %>%\n select(\n sex,\n location, atsi:diag1,\n starts_with(c(\"past\", \"hx\", \"fxhx\")),\n suicide_attempts, planned_final,\n isas_4:isas_7, siss\n ) %>%\n mutate(location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n map(~ tidy(chisq.test(., .fu$c))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n filter(across(ends_with(\"term\"), ~ .x == \".fu$c\")) %>%\n select(ends_with(c(\"statistic\", \"p.value\"))) %>%\n pivot_longer(\n everything(),\n names_to = c(\"var\", \"stat\"),\n names_sep = \"\\\\.\"\n ) %>%\n pivot_wider(\n names_from = stat,\n values_from = value\n )\n\n.fu_test_nonp %>%\n print(n = Inf)\n\n\n.fu_resd_use <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex),\n hxsa = case_when(\n suicide_attempts %in% \"0\" ~ \"no\",\n suicide_attempts %!in% c(NA, \"0\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n dx1 = case_when(\n diag1 %in% \"no diagnosis warranted\" ~ \"no\",\n diag1 %!in% c(NA, \"no diagnosis warranted\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n across(c(hxsa, dx1), as.factor),\n ) %>%\n select(\n location,\n sex,\n atsi:diag1, dx1,\n starts_with(c(\"past\", \"hx\", \"fxhx\")),\n suicide_attempts, planned_final,\n isas_4:isas_7, siss\n ) %>%\n mutate(\n location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n ## map(~ chisq.test(table(., .fu_bin$c))) %>%\n map(~ chisq.test(table(., .fu$c))) %>%\n map(`[`, c(\"stdres\", \"statistic\", \"p.value\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n rename(\n ## to match cont df\n value = stdres,\n p = p_value,\n x2_1 = x1,\n x3_1 = x2,\n x3_2 = x3\n )\n.fu_resd_use\n\n\n.lpa_m_sd_fu_nonp_full <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex),\n hxsa = case_when(\n suicide_attempts %in% \"0\" ~ \"no\",\n suicide_attempts %!in% c(NA, \"0\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n dx1 = case_when(\n diag1 %in% \"no diagnosis warranted\" ~ \"no\",\n diag1 %!in% c(NA, \"no diagnosis warranted\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n across(c(hxsa, dx1), as.factor),\n ) %>%\n select(where(is.factor)) %>%\n mutate(location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n pivot_longer(cols = -c) %>%\n count(name, value) %>%\n filter(!is.na(value)) %>%\n mutate(c = \"full\", .before = 1) %>%\n group_by(c, name) %>%\n rename(sd = n) %>%\n mutate(across(\n sd,\n ~ round(.x / sum(.x) * 100, 2),\n .names = \"m\"\n )) %>%\n mutate(\n row = row_number(),\n across(where(is.integer), as.double)\n ) %>%\n pivot_wider(\n names_from = c,\n values_from = c(sd, m),\n ) %>%\n select(-row)\n\n.lpa_m_sd_fu_nonp_full %>%\n print(n = Inf)\n\n.lpa_m_sd_fu_nonp <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex),\n hxsa = case_when(\n suicide_attempts %in% \"0\" ~ \"no\",\n suicide_attempts %!in% c(NA, \"0\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n dx1 = case_when(\n diag1 %in% \"no diagnosis warranted\" ~ \"no\",\n diag1 %!in% c(NA, \"no diagnosis warranted\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n across(c(hxsa, dx1), as.factor),\n ) %>%\n select(where(is.factor)) %>%\n mutate(location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n pivot_longer(cols = -c) %>%\n count(c, name, value) %>%\n filter(!is.na(value)) %>%\n group_by(c, name) %>%\n rename(sd = n) %>%\n pivot_wider(\n names_from = c,\n names_prefix = \"sd_\",\n values_from = sd,\n ) %>%\n ## full_join(.lpa_m_sd_fu_nonp_sex) %>%\n full_join(.lpa_m_sd_fu_nonp_full) %>%\n mutate(\n across(\n sd_1:sd_3,\n ~ round(.x / sd_full * 100, 2),\n .names = \"m_{1:3}\"\n )) %>%\n arrange(name, value) %>%\n group_by(name, value) %>%\n ## summarise(across(everything(), .coalesce_by_column)) %>%\n arrange(name, desc(m_full)) %>%\n ungroup() %>%\n select(name, value, ends_with(c(\"full\", \"1\", \"2\", \"3\")), everything()) %>%\n rename(var = name) %>%\n full_join(.fu_resd_use) %>%\n rename(name = var) %>%\n mutate(\n across(where(is.integer), as.double),\n value = case_when(\n str_detect(name, \"^(cald|past|hx|fxhx|planned)\") & value == 0 ~ \"no\",\n str_detect(name, \"^(cald|past|hx|fxhx|planned)\") & value == 1 ~ \"yes\",\n name == \"suicide_attempts\" & value == 0 ~ \"none\",\n name == \"suicide_attempts\" & value == 1 ~ \"one\",\n name == \"suicide_attempts\" & value == 2 ~ \"two\",\n name == \"suicide_attempts\" & value == 3 ~ \"three\",\n name == \"suicide_attempts\" & value == 4 ~ \"four or more\",\n name == \"living\" & value == 1 ~ \"family\",\n name == \"living\" & value == 2 ~ \"friends\",\n name == \"living\" & value == 3 ~ \"other(s)\",\n name == \"living\" & value == 4 ~ \"alone\",\n str_detect(name, \"^(isas_4|isas_5|isas_7)\") & value == 1 ~ \"yes\",\n str_detect(name, \"^(isas_4|isas_5|isas_7)\") & value == 2 ~ \"sometimes\",\n str_detect(name, \"^(isas_4|isas_5|isas_7)\") & value == 3 ~ \"no\",\n name == \"isas_6\" & value == 1 ~ \"00--1 hr\",\n name == \"isas_6\" & value == 2 ~ \"01--3 hrs\",\n name == \"isas_6\" & value == 3 ~ \"03--6 hrs\",\n name == \"isas_6\" & value == 4 ~ \"06--12 hrs\",\n name == \"isas_6\" & value == 5 ~ \"12--24 hrs\",\n name == \"isas_6\" & value == 6 ~ \"24+ hrs\",\n name == \"siss\" & value == 0 ~ \"si_none\",\n name == \"siss\" & value == 1 ~ \"si_vague\",\n name == \"siss\" & value == 2 ~ \"si_clear_no_plan\",\n name == \"siss\" & value == 3 ~ \"si_plans_no_int\",\n name == \"siss\" & value == 4 ~ \"si_plans_int\",\n name == \"siss\" & value == 5 ~ \"si_attempt\",\n TRUE ~ as.character(value)\n )\n ) %>%\n rename(var = name)\n\nfu_test_cont_nonp <-\n .lpa_m_sd_fu_nonp %>%\n mutate(test = \"fact\")\n\nfu_test_cont_nonp %>%\n print(n = Inf)\n\nfu_test_cont <-\n .lpa_n_fu %>%\n full_join(.lpa_m_sd_fu) %>%\n full_join(.fu_test) %>%\n full_join(.fu_test_cntrst) %>%\n mutate(\n test = \"cont\",\n bonf = 0.05\n )\n\n.fu_test_full <-\n fu_test_cont %>%\n rename(value = var) %>%\n mutate(var = \"\") %>%\n full_join(fu_test_cont_nonp) %>%\n relocate(c(var, value), .before = 1) %>%\n filter(var != \"sp_db\") %>%\n mutate(\n value = case_when(\n value == \"isas_1a\" ~ \"dsh_cutting\",\n value == \"isas_1b\" ~ \"dsh_biting\",\n value == \"isas_1c\" ~ \"dsh_burning\",\n value == \"isas_1d\" ~ \"dsh_carving\",\n value == \"isas_1e\" ~ \"dsh_pinching\",\n value == \"isas_1f\" ~ \"dsh_pull_hair\",\n value == \"isas_1g\" ~ \"dsh_scratch\",\n value == \"isas_1h\" ~ \"dsh_hit_self\",\n value == \"isas_1i\" ~ \"dsh_int_wound\",\n value == \"isas_1j\" ~ \"dsh_skin_rub\",\n value == \"isas_1k\" ~ \"dsh_needles\",\n value == \"isas_1l\" ~ \"dsh_swall_dang\",\n value == \"isas_1m\" ~ \"dsh_other\",\n value == \"isas_3a\" ~ \"dsh_age_onset\",\n TRUE ~ as.character(value)\n ),\n var = case_when(\n str_detect(value, \"^(is_|isas_|dsh_)\") ~ \"dsh\",\n var == \"isas_4\" ~ \"dsh_pain\",\n var == \"isas_5\" ~ \"dsh_alone\",\n var == \"isas_6\" ~ \"dsh_time_wait\",\n var == \"isas_7\" ~ \"dsh_want_stop\",\n value %in% c(\"dep\", \"anx\", \"str\") ~ \"mood\",\n value == \"total_session\" ~ \"total_session\",\n value == \"si\" ~ \"si\",\n value == \"siss_n\" ~ \"siss_n\",\n value == \"age\" ~ \"age\",\n value == \"Class size\" ~ \"Class size\",\n TRUE ~ as.character(var)\n ),\n type = case_when(\n var == \"Class size\" ~ \"Class size\",\n str_detect(var, \"^(age|atsi|cald|inc|liv|loc|rel|sex)\") ~ \"demo\",\n str_detect(var, \"^(mood|dsh|si|siss_n)\") ~ \"psych\",\n str_detect(var, \"^(fxhx|hx|suicide_|hxsa)\") ~ \"risk\",\n str_detect(var, \"^(total|diag|dx1|past|planned)\") ~ \"tx\"\n ),\n order = case_when(\n str_detect(var, \"^(dx1|diag|pastmed|planned)\") ~ 0,\n var == \"location\" ~ 0,\n var == \"Class size\" ~ 1,\n var == \"sex\" ~ 2,\n var == \"age\" ~ 3,\n var == \"atsi\" ~ 4,\n var == \"cald\" ~ 5,\n var == \"income\" ~ 6,\n var == \"living\" ~ 7,\n var == \"relation\" ~ 8,\n var == \"mood\" ~ 9,\n value == \"dsh_n\" ~ 10,\n value == \"isas_intra\" ~ 11,\n str_detect(value, \"^(is_)(aff|slf_p|ant|mrk)\") ~ 12,\n value == \"isas_social\" ~ 13,\n str_detect(value, \"^(is_)(slf_c|tou|aut|sen|int|rev|peer)\") ~ 14,\n value == \"isas_days\" ~ 15,\n str_detect(var, \"^(dsh_)(alone|want|time|pain)\") ~ 16,\n value == \"dsh_age_onset\" ~ 17,\n str_detect(value,\n \"^(dsh_)(cu|bi|bu|ca|pi|pu|sc|hi|in|sk|ne|sw)\") ~ 18,\n str_detect(value,\n \"dsh_other\") ~ 19,\n str_detect(var, \"hxdsh|fxhxdsh\") ~ 20,\n str_detect(var, \"^(hx|fxhx)(trauma|suic|sa)\") ~ 21,\n var == \"si\" ~ 24,\n value == \"si_clear_no_plan\" ~ 23,\n value == \"si_vague\" ~ 25,\n value == \"si_plans_no_int\" ~ 22,\n value == \"si_none\" ~ 26,\n value == \"si_plans_int\" ~ 27,\n value == \"si_attempt\" ~ 28,\n var == \"siss_n\" ~ 29,\n var == \"pastpsyc\" ~ 30,\n var == \"total_session\" ~ 31\n )\n ) %>%\n relocate(c(order, type), .before = 1) %>%\n arrange(order, type, var, desc(m_full))\n\n## generate the table for the report\n## load the manual group diff entries\n.fu_test_grp_dif <-\n readxl::read_xlsx(\n \"./data/tmp/anova_results.xlsx\",\n sheet = \"group_diff\"\n ) %>%\n select(order:value, group_diff) %>%\n mutate(order = as.double(order))\n\nfu_test <-\n .fu_test_full %>%\n select(-order) %>%\n mutate(\n across(\n starts_with(\"sd_\"),\n ~ case_when(\n .data$test == \"cont\" && .data$var != \"Class size\" ~\n sprintf(\"%.2f\", round(.x, 2)),\n TRUE ~ as.character(.x)\n )\n ),\n across(\n starts_with(\"m_\"), ~ sprintf(\"%.2f\", round(.x, 2))\n ),\n full = str_c(m_full, \" (\", sd_full, \")\"),\n low_si = str_c(m_1, \" (\", sd_1, \")\"),\n mod_si = str_c(m_2, \" (\", sd_2, \")\"),\n hi_si = str_c(m_3, \" (\", sd_3, \")\"),\n chisq_f = case_when(\n ## html tags for supescript\n p < 0.001 ~ paste0(sprintf(\"%.2f\", statistic), \"<sup>***</sup>\"),\n p < 0.01 ~ paste0(sprintf(\"%.2f\", statistic), \"<sup>**</sup>\"),\n p < 0.05 ~ paste0(sprintf(\"%.2f\", statistic), \"<sup>*</sup>\"),\n p > 0.05 ~ paste0(sprintf(\"%.2f\", statistic), \"\"),\n TRUE ~ NA_character_\n ),\n ) %>%\n full_join(.fu_test_grp_dif) %>%\n ## relocate(starts_with(\"group_diff\"), .after = everything()) %>%\n group_by(order, type, var, value) %>%\n ## arrange(order, type, var, desc(m_3), .by_group = TRUE) %>%\n arrange(order, type, var, .by_group = TRUE) %>%\n ungroup() %>%\n relocate(value2, .before = value) %>%\n relocate(order, .before = 1)\n\n\n## MULTINOMIAL LOGISTIC REGRESSION\n## =============================================================================\n\n.fu_<-\n .fu %>%\n filter(!is.na(c)) %>%\n mutate(\n across(where(is.factor), as.character),\n suicide_attempts = as.numeric(suicide_attempts),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n across(starts_with(c(\"cald\", \"past\", \"hx\", \"fxhx\", \"planned\")),\n ~ case_when(\n .x == 0 ~ \"no\",\n .x == 1 ~ \"yes\",\n TRUE ~ NA_character_\n )\n ),\n hxsa = case_when(\n suicide_attempts < 1 ~ \"no\",\n suicide_attempts >= 1 ~ \"yes\",\n TRUE ~ NA_character_\n ),\n ## living = case_when(\n ## living == 1 ~ \"family\",\n ## living == 2 ~ \"other\",\n ## living == 3 ~ \"other\",\n ## living == 4 ~ \"other\",\n ## TRUE ~ NA_character_\n ## ),\n across(starts_with(\"isas_\") & ends_with(c(\"4\", \"5\", \"7\")),\n ~ case_when(\n .x == 1 ~ \"yes\",\n .x == 2 ~ \"sometimes\",\n .x == 3 ~ \"no\",\n TRUE ~ NA_character_\n )\n ),\n across(where(is.character), as.factor),\n ) %>%\n rowwise() %>%\n mutate(\n dsh_n = sum(c_across(starts_with(\"isas_1\")) > 0, na.rm = TRUE)\n ) %>%\n ungroup() %>%\n rename(\n dsh_pain = isas_4,\n dsh_alone = isas_5,\n dsh_stop = isas_7\n ) %>%\n ## str()\n mutate(\n across(c(living, sex), ~ relevel(.x, ref = 2)),\n across(starts_with(\"dsh\") & -dsh_n, ~ relevel(.x, ref = 2)),\n across(relation, ~ relevel(.x, ref = 2)),\n across(atsi, ~ relevel(.x, ref = 2)),\n )\n\n.fu_ %>%\n select(starts_with(\"fx\")) %>%\n summary()\n\n.fu_log <-\n dfidx::dfidx(.fu_, choice = \"c\", shape = \"wide\")\n\n.fu_log_test_3 <-\n mlogit(\n c ~ 1 |\n dep + anx + str + age + si + sex +\n atsi + cald + income + living + relation +\n dsh_pain + dsh_alone + dsh_stop + dsh_n +\n isas_intra + isas_social +\n hxdsh + hxsa + hxsuicide + hxtrauma +\n ## fxhxdsh + fxhxsuicide + pastmed +\n pastpsyc,\n data = .fu_log,\n reflevel = \"3\",\n na.action = na.exclude\n )\n\nfu_ml_sum_3 <-\n summary(.fu_log_test_3)\n\n.fu_ml_3 <-\n fu_ml_sum_3$CoefTable %>%\n as.data.frame() %>%\n mutate(var = attr(fu_ml_sum_3$coefficients, \"names\"), .before = 1) %>%\n tibble() %>%\n clean_names() %>%\n mutate(\n exp_b = round(exp(estimate), 2), .after = var\n ) %>%\n full_join(\n tibble(clean_names(cbind(\n data.frame(var = attr(exp(confint(.fu_log_test_3)), \"dimnames\")[[1]]),\n data.frame(exp(confint(.fu_log_test_3)))\n )))\n ) %>%\n select(var, exp_b, x2_5, x97_5, pr_z) %>%\n filter(!str_detect(var, \"(Intercept)\")) %>%\n mutate(\n name = str_replace_all(var, \":\\\\d\", \"\"),\n comp = str_extract(var, \":\\\\d\"),\n comp = str_replace(comp, \":\", \"3v\"),\n ) %>%\n select(-var) %>%\n ## NOTE consider transfering below to rmd for table pres\n pivot_wider(\n names_from = comp,\n values_from = c(exp_b, x2_5, x97_5, pr_z)\n ) %>%\n select(name, ends_with(c(\"3v1\", \"3v2\"))) %>%\n mutate(\n ## across(where(is.double), ~ round(.x, 2)),\n across(starts_with(\"x\"), ~ round(.x, 2)),\n across(starts_with(\"pr_\"), ~ round(.x, 3))\n )\n\n.fu_ml_3 %>%\n filter(\n if_any(\n starts_with(\"pr_z\"),\n ~.x < 0.05)\n )\n\n.fu_log_test_2 <-\n mlogit(\n c ~ 1 |\n dep + anx + str + age + si + sex +\n atsi + cald + income + living + relation +\n dsh_pain + dsh_alone + dsh_stop + dsh_n +\n isas_intra + isas_social +\n hxdsh + hxsa + hxsuicide + hxtrauma +\n ## fxhxdsh + fxhxsuicide + pastmed +\n pastpsyc,\n data = .fu_log,\n reflevel = \"2\",\n na.action = na.exclude\n )\n\nfu_ml_sum_2 <-\n summary(.fu_log_test_2)\n\n.fu_ml_2 <-\n fu_ml_sum_2$CoefTable %>%\n as.data.frame() %>%\n mutate(var = attr(fu_ml_sum_2$coefficients, \"names\"), .before = 1) %>%\n tibble() %>%\n clean_names() %>%\n mutate(\n exp_b = round(exp(estimate), 2), .after = var\n ) %>%\n full_join(\n tibble(clean_names(cbind(\n data.frame(var = attr(exp(confint(.fu_log_test_2)), \"dimnames\")[[1]]),\n data.frame(exp(confint(.fu_log_test_2)))\n )))\n ) %>%\n select(var, exp_b, x2_5, x97_5, pr_z) %>%\n filter(!str_detect(var, \"(Intercept)\")) %>%\n mutate(\n name = str_replace_all(var, \":\\\\d\", \"\"),\n comp = str_extract(var, \":\\\\d\"),\n comp = str_replace(comp, \":\", \"2v\"),\n ) %>%\n select(-var) %>%\n ## NOTE consider transfering below to rmd for table pres\n pivot_wider(\n names_from = comp,\n values_from = c(exp_b, x2_5, x97_5, pr_z)\n ) %>%\n select(name, ends_with(c(\"2v1\"))) %>%\n mutate(\n across(starts_with(\"x\"), ~ round(.x, 2)),\n across(starts_with(\"pr_\"), ~ round(.x, 3))\n )\n\n.fu_ml_2 %>%\n filter(\n if_any(\n starts_with(\"pr_z\"),\n ~ .x < 0.05)\n )\n\n\n## table for presentation\n.fu_ml_order <-\n readxl::read_xlsx(\n \"./data/tmp/anova_results.xlsx\",\n sheet = \"ml_table\"\n ) %>%\n select(order:name2) %>%\n mutate(order = as.double(order))\n\nfu_ml <-\n .fu_ml_3 %>%\n full_join(.fu_ml_2) %>%\n select(name, ends_with(c(\"2v1\", \"3v1\", \"3v2\"))) %>%\n full_join(.fu_ml_order) %>%\n arrange(order) %>%\n relocate(c(order, type, name2), .before = name)\n\nfu_ml %>%\n print(n = Inf)\n\n\n.fu_log_test_1 <-\n mlogit(\n c ~ 1 |\n dep + anx + str + age + si + sex +\n atsi + cald + income + living + relation +\n dsh_pain + dsh_alone + dsh_stop + dsh_n +\n isas_intra + isas_social +\n hxdsh + hxsa + hxsuicide + hxtrauma +\n ## fxhxdsh + fxhxsuicide + pastmed +\n pastpsyc,\n data = .fu_log,\n reflevel = \"1\",\n )\n\nfu_ml_sum_1 <-\n summary(.fu_log_test_1)\n\n.fu_ml_1 <-\n fu_ml_sum_1$CoefTable %>%\n as.data.frame() %>%\n mutate(var = attr(fu_ml_sum_1$coefficients, \"names\"), .before = 1) %>%\n tibble() %>%\n clean_names() %>%\n mutate(\n exp_b = round(exp(estimate), 2), .after = var\n ) %>%\n full_join(\n tibble(clean_names(cbind(\n data.frame(var = attr(exp(confint(.fu_log_test_1)), \"dimnames\")[[1]]),\n data.frame(exp(confint(.fu_log_test_1)))\n )))\n ) %>%\n select(var, exp_b, x2_5, x97_5, pr_z) %>%\n filter(!str_detect(var, \"(Intercept)\")) %>%\n mutate(\n name = str_replace_all(var, \":\\\\d\", \"\"),\n comp = str_extract(var, \":\\\\d\"),\n comp = str_replace(comp, \":\", \"1v\"),\n ) %>%\n select(-var) %>%\n ## NOTE consider transfering below to rmd for table pres\n pivot_wider(\n names_from = comp,\n values_from = c(exp_b, x2_5, x97_5, pr_z)\n ) %>%\n ## select(name, ends_with(c(\"1v1\"))) %>%\n mutate(\n across(starts_with(\"x\"), ~ round(.x, 2)),\n across(starts_with(\"pr_\"), ~ round(.x, 3))\n )\n\n.fu_ml_1 %>%\n print(n = Inf)\n\n## CREATE FU GRAPHS\n## -----------------------------------------------------------------------------\n## continuous\n.fu_z <-\n .fu_ %>%\n mutate(\n across(\n c(\n pb, tb, bhs, ac,\n dep, anx, str,\n isas_days,\n is_aff_reg, is_ant_dis, is_ant_sui, is_slf_pun,\n dsh_n,\n isas_intra, isas_social,\n si,\n total_session\n ),\n ## where(is.numeric),\n ~ scale(.x),\n .names = \"{.col}_z\"\n )\n ) %>%\n select(c, contains(\"_z\")) %>%\n mutate(\n across(\n everything(),\n ~ str_remove(.x, \"_z\")\n ),\n across(\n everything(), as.double\n )\n ) %>%\n rename_with(~str_remove(., \"_z\"))\n\n## cat\n.fu_plot_cat_names <-\n c(\"DSH[hx]\", \"SA[hx]\", \"SI[hx]\", \"TRM[hx]\", \"FAM\", \"REL\", \"FEM\")\n\nfu_plot_cat <-\n .fu_ %>%\n select(\n c, sex,\n living, relation, dsh_stop, dsh_alone,\n hxsa, hxdsh, hxsuicide, hxtrauma,\n ) %>%\n group_by(c) %>%\n summarise(\n n = n(),\n sex_fem =\n length(sex[sex == \"Female\"]) /\n length(sex[!is.na(sex)]) * 100,\n liv_fam =\n length(living[living == \"family\"]) /\n length(living[!is.na(living)]) * 100,\n rel_yes =\n length(relation[str_detect(relation, \"relation\")]) /\n length(relation[!is.na(relation)]) * 100,\n across(starts_with(c(\"hx\", \"past\")),\n ~ length(.x[str_detect(.x, \"yes\")]) /\n length(.x[!is.na(.x)]) * 100\n )\n ) %>%\n pivot_longer(\n -c\n ) %>%\n filter(name != \"n\") %>%\n ggplot(\n aes(\n x = name,\n y = value,\n linetype = factor(c),\n group = factor(c),\n fill = factor(c),\n shape = factor(c)\n )\n ) +\n geom_bar(position = position_dodge(), stat=\"identity\") +\n theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ggtitle(\"Significant caterogical measures\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n scale_x_discrete(label = parse(text = .fu_plot_cat_names))\n\n## http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/\n\npd <-\n position_dodge(0.01) # move them .05 to the left and right\n\n.lpa_plot_cont_names <-\n c(\n \"PB\", \"TB\", \"BHS\", \"AC\"\n )\n\nlpa_plot_cont <-\n .fu_z %>%\n pivot_longer(\n cols = -c\n ) %>%\n filter(name %in% c(\"pb\", \"tb\", \"bhs\", \"ac\")) %>%\n group_by(name, c) %>%\n summarise(\n across(\n everything(),\n list(\n N = ~ length(.x),\n mean = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n se = ~ sd(.x, na.rm = TRUE) / sqrt(length(.x)),\n ci = ~ (sd(.x, na.rm = TRUE) / sqrt(length(.x))) * qt(0.95/2 + .5, length(.x) - 1)\n ),\n .names = \"{.col}.{.fn}\")\n ) %>%\n ungroup() %>%\n group_by(name) %>%\n mutate(\n ID = case_when(\n name == \"pb\" ~ 1,\n name == \"tb\" ~ 2,\n name == \"bhs\" ~ 3,\n name == \"ac\" ~ 4\n ),\n .before = name,\n name = factor(\n name,\n levels = c(\n \"pb\", \"tb\", \"bhs\", \"ac\"\n )\n )\n ) %>%\n ungroup() %>%\n ggplot(\n aes(\n x = name,\n y = value.mean,\n colour = factor(c),\n shape = factor(c),\n linetype = factor(c),\n )\n ) +\n geom_errorbar(\n aes(\n ymin = value.mean - value.se,\n ymax = value.mean + value.se\n ), width = .1, position = pd\n ) +\n geom_point(position = pd) +\n geom_line(\n aes(\n x = ID,\n y = value.mean,\n colour = as.factor(c)\n ),\n position = pd) +\n theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ## ggtitle(\"\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n scale_x_discrete(label = .lpa_plot_cont_names)\nlpa_plot_cont\n\nggsave(\n \"./figures/lpa_plot_cont.png\",\n lpa_plot_cont\n)\n\nggsave(\n \"./figures/lpa_plot_cont.pdf\",\n lpa_plot_cont\n)\n\n.fu_plot_cont_names <-\n c(\n ## \"PB\", \"TB\", \"BHS\", \"AC\",\n \"DEP\", \"ANX\", \"STR\", \"SI\",\n ## \"DAYS\",\n \"METH\", \"INTRA\",\n \"AR\", \"AD\", \"AS\", \"SP\",\n \"INTER\"\n ## \"SESSION\"\n )\n\nfu_plot_cont <-\n .fu_z %>%\n pivot_longer(\n cols = -c\n ) %>%\n filter(!name %in% c(\"siss_n\", \"pb\", \"tb\", \"bhs\", \"ac\", \"isas_days\", \"total_session\")) %>%\n group_by(name, c) %>%\n summarise(\n across(\n everything(),\n list(\n N = ~ length(.x),\n mean = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n se = ~ sd(.x, na.rm = TRUE) / sqrt(length(.x)),\n ci = ~ (sd(.x, na.rm = TRUE) / sqrt(length(.x))) * qt(0.95/2 + .5, length(.x) - 1)\n ),\n .names = \"{.col}.{.fn}\")\n ) %>%\n ## print(n = Inf) %>%\n ungroup() %>%\n group_by(name) %>%\n mutate(\n ID = case_when(\n name == \"dep\" ~ 1,\n name == \"anx\" ~ 2,\n name == \"str\" ~ 3,\n name == \"si\" ~ 4,\n name == \"dsh_n\" ~ 5,\n name == \"isas_intra\" ~ 6,\n name == \"is_aff_reg\" ~ 7,\n name == \"is_ant_dis\" ~ 8,\n name == \"is_ant_sui\" ~ 9,\n name == \"is_slf_pun\" ~ 10,\n name == \"isas_social\" ~ 11\n ),\n .before = name,\n name = factor(\n name,\n levels = c(\n ## \"pb\", \"tb\", \"bhs\", \"ac\",\n \"dep\", \"anx\", \"str\", \"si\",\n ## \"siss_n\",\n ## \"isas_days\",\n \"dsh_n\", \"isas_intra\",\n \"is_aff_reg\", \"is_ant_dis\", \"is_ant_sui\", \"is_slf_pun\",\n \"isas_social\"\n ## \"total_session\"\n )\n )\n ) %>%\n ungroup() %>%\n ## filter(name %!in% c(\"pb\", \"tb\", \"bhs\", \"ac\")) %>%\n ggplot(\n aes(\n x = name,\n y = value.mean,\n colour = factor(c),\n shape = factor(c),\n linetype = factor(c)\n )\n ) +\n ## geom_bar(position=position_dodge(), stat=\"identity\") +\n geom_errorbar(\n aes(\n ymin = value.mean - value.se,\n ymax = value.mean + value.se\n ), width = .1, position = pd\n ) +\n geom_point(\n position = pd,\n ) +\n geom_line(\n aes(\n x = ID,\n y = value.mean,\n ),\n position = pd\n ) +\n labs(colour = \"Class\", shape = \"Class\", linetype = \"Class\") +\n ## scale_colour_viridis(\n ## discrete = TRUE,\n ## alpha = 0.7,\n ## )\n ## theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ggtitle(\"Significant continous measures\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n scale_x_discrete(label = .fu_plot_cont_names)\n\ndesign = \"\n112\n112\n112\n332\n332\n\"\n\n.lca_plot_names <-\n ## c(\"Capability\", \"Hopelessness\", \"Burdensomeness\", \"Belongingness\")\n c(\"AC\", \"BHS\", \"PB\", \"TB\")\n\nlpa_fig_plot <-\n lpa_fig +\n theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ## ggtitle(\"\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n ## ylim(\"\") +\n scale_x_discrete(label = .lca_plot_names) +\n theme(\n ## legend.position = \"bottom\",\n text = element_text(size = 7),\n title = element_text(size = 8),\n )\n\nlpa_fig_plot +\n theme(\n ## legend.position = \"bottom\",\n text = element_text(size = 7),\n title = element_text(size = 8),\n )\n\nggsave(\n \"./figures/lpa_fig_plot.png\",\n lpa_fig_plot\n)\n\nlpa_combi_plot <-\n ## lpa_plot_cont +\n fu_plot_cont +\n fu_plot_cat +\n plot_layout(\n ## guides = \"collect\",\n ## design = design\n ## heights = c(2, 1)\n widths = c(2, 1)\n ) &\n theme(\n legend.position = \"bottom\",\n text = element_text(size = 7),\n title = element_text(size = 8),\n )\n\n\nggsave(\n \"./figures/lpa_combi_plot.png\",\n lpa_combi_plot\n)\n\n### ============================================================================\n### OUTCOMES\n### ============================================================================\nlpa_t0t1_outc <-\n lpa_fin %>%\n ## filter no class var\n filter(!is.na(c)) %>%\n group_by(ID) %>%\n ## filter first and last session\n filter(session == 1 | final_session == 1) %>%\n ungroup() %>%\n mutate(\n ## replace NA session value as 5 as only 5 sessions of data available\n session = if_else(is.na(session), 5, session),\n ## create conditional variable to drop those with only one session\n drop = if_else(\n session == 1 & final_session == 1,\n 1, 0\n ), .before = ID,\n hxsa = case_when(\n suicide_attempts <= 0 ~ 0,\n suicide_attempts >= 1 ~ 1,\n TRUE ~ NA_real_\n ),\n living = factor(case_when(\n living == 1 ~ \"family\",\n living == 2 ~ \"other\",\n living == 3 ~ \"other\",\n living == 4 ~ \"other\",\n TRUE ~ NA_character_\n )),\n relation = factor(case_when(\n relation == 1 ~ \"single/separated/divorced\",\n relation == 2 ~ \"relationship/married/defacto\",\n relation == 3 ~ \"relationship/married/defacto\",\n relation == 4 ~ \"single/separated/divorced\",\n TRUE ~ NA_character_\n )),\n sex = factor(if_else(sex == 2, NA_character_, as.character(sex)))\n ) %>%\n ## select those with more than one session\n filter(drop == 0) %>%\n group_by(ID) %>%\n mutate(age = min(age)) %>%\n ungroup() %>%\n ## names() %>%\n ## select vars for mi\n dplyr::select(\n c, ID, age, sex, session, final_session,\n planned_final, pb, ac, bhs, tb, dep, anx, str,\n isas_intra, isas_social,\n relation, living,\n hxsa, hxsuicide, hxdsh, hxtrauma,\n si\n ) %>%\n group_by(ID) %>%\n mutate(\n cov_sess = max(session), .after = final_session,\n ## hold AC as constant at time 1 to include as a covariate\n ac = ac[session == 1]\n ) %>%\n ungroup() %>%\n group_by(final_session) %>%\n mutate(\n across(\n c(\n si, pb, tb, bhs, dep, anx, str, ac, isas_social, isas_intra,\n age, cov_sess, hxsa, hxsuicide, hxdsh, hxtrauma\n ),\n ~ scale(.x),\n .names = \"{.col}_z\"\n ),\n ) %>%\n ungroup() %>%\n mutate(\n fin_ses_c = with(\n .,\n factor(final_session):factor(c)\n ),\n .after = c,\n across(\n c(c, sex, ID, final_session, planned_final),\n ~ as.factor(.x)\n ),\n )\n\nlpa_t0t1_outc_w <-\n lpa_t0t1_outc %>%\n dplyr::select(-fin_ses_c) %>%\n panelr::panel_data(\n id = ID,\n wave = final_session,\n ) %>%\n panelr::widen_panel(\n separator = \"_\",\n ignore.attributes = FALSE,\n varying = NULL) %>%\n ungroup() %>%\n ## need to convert to data frame as panelr creates a list object\n as.data.frame() %>%\n ## use age at intake, drop secon ac as not routinely collected\n dplyr::select(-planned_final_0)\n\n\nlpa_t0t1_outc %>% str()\nlpa_t0t1_outc_w %>% str()\n\n\nlpa_t0t1_outc_hisi <-\n lpa_t0t1_outc_w_hisi %>%\n select(!contains(\"_z\")) %>%\n panelr::long_panel(\n prefix = \"_\",\n periods = c(\"0\", \"1\"),\n id = \"ID\",\n wave = \"final_session\",\n label_location = \"end\"\n ) %>%\n ungroup() %>%\n mutate(\n c2 = case_when(\n c == \"2\" ~ \"lomi\",\n c == \"1\" ~ \"lomi\",\n c == \"3\" ~ \"high\",\n TRUE ~ NA_character_\n ), .after = c,\n c2 = as.factor(c2),\n ) %>%\n mutate(\n c3 = case_when(\n c == \"2\" ~ \"2\",\n c == \"3\" ~ \"3\",\n TRUE ~ NA_character_\n ), .after = c2,\n c3 = as.factor(c3),\n ) %>%\n mutate(fin_ses_c = with(., c:final_session), .after = c3) %>%\n mutate(fin_ses_c2 = with(., c2:final_session), .after = fin_ses_c) %>%\n mutate(fin_ses_c3 = with(., c3:final_session), .after = fin_ses_c2) %>%\n mutate(\n across(c(pb, ac, tb, bhs, dep, anx, str, si),\n ~scale(.x),\n .names = \"{.col}_z\"\n )\n )\n\nlpa_t0t1_outc %>%\n filter(final_session == 0) %>%\n group_by(c) %>%\n summarise(\n full = n_distinct(ID),\n hisi = n_distinct(ID[si > 8])\n ) %>%\n adorn_totals(\"row\")\n\nlpa_t0t1_outc_hisi %>%\n group_by(c) %>%\n summarise(\n n = n_distinct(ID)\n ) %>%\n adorn_totals(\"row\")\n\nlpa_t0t1_outc_w_hisi %>%\n select(ID, age, c, starts_with(c(\"pb\", \"si\", \"ac\", \"tb\", \"dep\"))) %>%\n .summary_by(c)\n\n## DISTRIBUTION PLOTS OF ANOVA VARS\n## -----------------------------------------------------------------------------\nplot_iter <-\n lpa_t0t1_outc %>%\n dplyr::select(\n ## ends_with(\"_z\") & !starts_with(\"isas\"),\n si, pb, tb, bhs, dep, anx,\n str, age, sex, cov_sess\n ) %>%\n names() %>%\n purrr::set_names()\n\ndens_fun <-\n function(dat, x, filly, facet) {\n ggplot(dat, aes(x = .data[[x]], fill = .data[[filly]])) +\n geom_density(alpha = 0.3) +\n facet_grid(.data[[facet]] ~ .) +\n theme_bw()\n }\n\n.plots_lpa_long <-\n map(plot_iter, ~ dens_fun(lpa_t0t1_outc, .x, \"c\", \"final_session\"))\n\nplots_lpa_long <-\n wrap_plots(.plots_lpa_long, guides = \"collect\")\n\nggsave(\n \"./figures/plots_lpa_long.png\",\n plots_lpa_long\n\n\n## PLOT CHANGE OVER TIME\n## =============================================================================\n## TODO get plots\n\n## LME4 PACKAGE\n## =============================================================================\n## Mixed Effects ANOVA\nstep_data <-\n lpa_t0t1_outc\n\nno_na_step_data <-\n na.omit(\n step_data[c(\n \"si\", \"pb\", \"bhs\", \"tb\", \"ac\", \"sex\",\n \"age\",\n \"dep\", \"anx\", \"str\",\n \"isas_social\", \"isas_intra\", \"hxsa\", \"hxdsh\", \"hxtrauma\", \"hxsuicide\",\n \"cov_sess\", \"c\",\n \"si_z\", \"pb_z\", \"bhs_z\", \"tb_z\", \"ac_z\", \"sex_z\",\n \"final_session\", \"age_z\",\n \"dep_z\", \"anx_z\", \"str_z\",\n \"isas_social_z\", \"isas_intra_z\", \"hxsa_z\", \"hxdsh_z\", \"hxtrauma_z\", \"hxsuicide_z\",\n \"cov_sess_z\"\n )])\n\n\n\n.tt0 <-\n no_na_step_data %>%\n filter(final_session == 0)\n\n.tt1 <-\n no_na_step_data %>%\n filter(final_session == 1)\n\n.tt0_lmer <-\n lmer(\n si ~\n pb + tb + bhs + dep + anx + str +\n ac + cov_sess + age + sex +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n isas_social + isas_intra +\n (1 | c),\n data = .tt0,\n na.action = na.exclude,\n REML = FALSE\n )\n\n.tt1_lmer <-\n lmer(\n si ~\n pb + tb + bhs + dep + anx + str +\n ac + cov_sess + age + sex +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n isas_social + isas_intra +\n (1 | c),\n data = .tt1,\n na.action = na.exclude,\n REML = FALSE\n )\n\n\ntt0_fixmodel <-\n lm(formula(.tt0_lmer, fixed.only=TRUE),\n data=eval(getCall(.tt0_lmer)$data))\n\nstep_model_tt0 <-\n stats::step(\n tt0_fixmodel,\n direction = \"both\",\n trace = 0,\n test = \"F\"\n )\n\nstep_model_tt0$anova\nstep_model_tt0$coefficients\nsummary(step_model_tt0)\n\ntt1_fixmodel <-\n lm(formula(.tt1_lmer, fixed.only=TRUE),\n data=eval(getCall(.tt1_lmer)$data))\n\nstep_model_tt1 <-\n stats::step(\n tt1_fixmodel,\n direction = \"both\",\n trace = 0,\n test = \"F\"\n )\n\nstep_model_tt1$anova\nstep_model_tt1$coefficients\nsummary(step_model_tt1)\n\n\nlpa_out_ni_step <-\n lmer(\n si ~\n pb + tb + bhs +\n dep + anx + str +\n ac + cov_sess + age + sex +\n isas_social + isas_intra +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n (1 | final_session/c),\n ## (1 | final_session/c),\n data = no_na_step_data,\n ## na.action = na.omit,\n ## control = lmerControl(optimizer = \"Nelder_Mead\"),\n REML = FALSE # use ML estimation\n )\n\nfixmodel <-\n lm(formula(lpa_out_ni_step, fixed.only=TRUE),\n data=eval(getCall(lpa_out_ni_step)$data))\n\nstep_model <-\n stats::step(\n fixmodel,\n direction = \"both\",\n trace = 0,\n test = \"F\"\n )\n\nstep_model$anova\nstep_model$coefficients\nsummary(step_model)\n\nlpa_out_step <-\n lmer(\n si_z ~\n pb_z + dep_z + ac_z + cov_sess_z + age_z +\n isas_social_z + hxsa_z +\n (1 | final_session/c),\n ## (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\n## NOTE lmerTest used to add t and p values\ncoef(summary(lpa_out_step, correlation = FALSE))\ncoef(summary(as(lpa_out_step, \"merModLmerTest\")))\nconfint(lpa_out_step, method = \"Wald\", oldNames = FALSE)\nsummary(as(lpa_out_step, \"merModLmerTest\"))\n\nlpa_out_ni <-\n lmer(\n si ~\n pb + tb + bhs + dep +\n ac + cov_sess + age + sex +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(lpa_out_ni)\n\n## NOTE lmerTest used to add t and p values\ncoef(summary(lpa_out_ni, correlation = FALSE))\ncoef(summary(as(lpa_out_ni, \"merModLmerTest\")))\nconfint(lpa_out_ni, method = \"Wald\", oldNames = FALSE)\nsummary(as(lpa_out_ni, \"merModLmerTest\"))\nvcov(lpa_out_ni)\n\n## calculate ICC for lmer model object\nicc_lmer <- function(m){\n vc <- as.data.frame((VarCorr(m)))\n l <- vc$vcov\n tibble(grp=vc$grp, icc=sapply(l, function(x){x/sum(l)}))\n}\n\nicc_lmer(lpa_out_ipt)\nVarCorr(lpa_out_ipt)\n\nrcompanion::nagelkerke(\n lpa_out_ni,\n lpa_out_step\n)\n\n## fit with all optimizers\nall_lpa <-\n allFit(lpa_out_ni)\nss <-\n summary(all_lpa)\nss$which.OK ## logical vector: which optimizers worked?\n## the other components only contain values for the optimizers that worked\nss$llik ## vector of log-likelihoods\nss$fixef ## table of fixed effects\nss$sdcor ## table of random effect SDs and correlations\nss$theta ## table of random effects parameters, Cholesky scale\n\nlpa_out_ni_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ni, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goews with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_out_ni_tidy %>%\n print(n = Inf)\n\n## by Ben Bolker\n## stackoverflow.com/questions/25142901/standardized-coefficients-for-lmer-model\nstdCoef.merMod <-\n function(object) {\n sdy <- sd(getME(object,\"y\"))\n sdx <- apply(getME(object,\"X\"), 2, sd)\n sc <- fixef(object)*sdx/sdy\n se.fixef <- coef(summary(object))[,\"Std. Error\"]\n se <- se.fixef*sdx/sdy\n return(data.frame(stdcoef=sc, stdse=se))\n}\n\nlpa_out_ni %>%\n stdCoef.merMod()\n\nlpa_out_ipt_hisi <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session), # does not provide singularity error\n ## (1 | final_session/c), # same results as above, with singularity error\n ## (1 | final_session/c2), # same results as above, with singularity error\n ## (1 | final_session/c3), # only pb sig, with ac almost sig, signul error\n ## data = lpa_t0t1_outc_hisi[!is.na(lpa_t0t1_outc_hisi$c3), ], # c3 data\n data = lpa_t0t1_outc_hisi,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\n\nsummary(as(lpa_out_ipt_hisi, \"merModLmerTest\"))\n## singular fit potentially becuase of low cell numbers in c2\n## recommendation when singular fit to simplify the model, therefore nesting removed\nVarCorr(lpa_out_ipt_hisi)\n\n\nlpa_out_ipt_hisi_plot <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n ac_z +\n ## cov_sess + age + sex +\n (1 | final_session), # does not provide singularity error\n ## (1 | final_session/c), # same results as above, with singularity error\n ## (1 | final_session/c2), # same results as above, with singularity error\n ## (1 | final_session/c3), # only pb sig, with ac almost sig, signul error\n ## data = lpa_t0t1_outc_hisi[!is.na(lpa_t0t1_outc_hisi$c3), ], # c3 data\n data = lpa_t0t1_outc_hisi,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(as(lpa_out_ipt_hisi_plot, \"merModLmerTest\"))\n\nlpa_int_plot_hihsi <-\n interact_plot(\n lpa_out_ipt_hisi,\n pred = pb_z,\n modx = ac_z,\n mod2 = tb_z,\n plot.points = FALSE\n ) +\n theme(\n text = element_text(size = 9),\n ) +\n labs(\n title = \"Meets SI Cut-off\"\n )\n\n\n\nlpa_int_plot2 <-\n interact_plot(\n lpa_out_ipt,\n pred = pb_z,\n modx = ac_z,\n plot.points = FALSE\n ) +\n theme(\n text = element_text(size = 9),\n )+\n labs(\n title = \"SI ~ PB * AC\"\n )\n\nlpa_int_plot <-\n interact_plot(\n lpa_out_ipt,\n pred = pb_z,\n modx = ac_z,\n plot.points = FALSE\n ) +\n theme(\n legend.position = \"none\",\n text = element_text(size = 9),\n )+\n labs(\n title = \"Full Sample\"\n )\n\nlpa_int_plot_combi <-\n lpa_int_plot +\n lpa_int_plot_hihsi &\n plot_layout(\n nrow = 2,\n ) &\n theme(\n text = element_text(size = 9),\n title = element_text(size = 9),\n )\n\n\nggsave(\n \"./figures/lpa_int_plot.png\",\n lpa_int_plot2\n)\n\nggsave(\n \"./figures/lpa_int_plot_combi.png\",\n lpa_int_plot_combi\n)\n\nlpa_t0t1_outc_hisi %>%\n group_by(final_session) %>%\n summarise(\n full = n_distinct(ID),\n hisi = n_distinct(ID[si > 8 & final_session %in% \"0\"]),\n ) %>%\n adorn_totals(\"row\") %>%\n mutate(\n prc = round(hisi/full * 100, 1)\n )\n\nlpa_main_fx_c_ipt_hisi <-\n lpa_t0t1_outc_hisi %>%\n ## split(.$final_session) %>%\n split(.$c) %>%\n map(~\n ## lm(\n lmer(\n si ~\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n ## map(~ summary(.)) %>%\n map(~ tidy(as(., \"merModLmerTest\"), conf.method = \"profile\", conf.int = TRUE)) %>%\n bind_rows(., .id = \"var\")\n\nlpa_main_fx_c_ipt_hisi %>%\n filter(p.value < 0.5) %>%\n print(n = Inf)\n\nlpa_t0t1_outc %>%\n group_by(c) %>%\n summarise(\n n = n_distinct(ID)\n )\n\n## START MULTILEVEL MODEL STUFF\n## -----------------------------------------------------------------------------\nlpa_out_ipt <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(as(lpa_out_ipt, \"merModLmerTest\"))\n\nlpa_out_ipt_std <-\n lpa_out_ipt %>%\n stdCoef.merMod()\n\nlpa_out_ipt_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ipt, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goes with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_out_ipt_tidy %>%\n print(n = Inf)\n\nlpa_out_ipt_rs <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (c | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\n## complare with added slope, random intercept enough\nanova(lpa_out_ipt_rs, lpa_out_ipt)\n\nsummary(as(lpa_out_ipt_rs, \"merModLmerTest\"))\n\nlpa_out_ipt_rs_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ipt_rs, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goes with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_main_fx_c_ipt <-\n lpa_t0t1_outc %>%\n split(.$c) %>%\n map(~\n lmer(\n si ~\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n map(~ tidy(as(., \"merModLmerTest\"), conf.method = \"profile\", conf.int = TRUE)) %>%\n bind_rows(., .id = \"var\")\n\nlpa_main_fx_c_ipt %>%\n filter(p.value < 0.05) %>%\n print(n = Inf)\n\nprint(VarCorr(lpa_out_ipt), comp = c(\"Variance\", \"Std.Dev\"), digits = 2)\nVarCorr(lpa_out_ipt)\n\nlpa_out_ipt_all <-\n bind_rows(\n lpa_out_ipt_tidy,\n lpa_main_fx_c_ipt\n ) %>%\n clean_names()\n\n\nlpa_int_pb_ac <-\n interact_plot(\n lpa_out_ipt,\n pred = pb_z,\n modx = ac_z,\n ## cluster = \"final_session/c\",\n modx.values = \"terciles\",\n plot.points = FALSE\n )\n\nlpa_int_tb_ac <-\n interact_plot(\n lpa_out_ipt,\n pred = tb_z,\n modx = ac_z,\n plot.points = FALSE\n )\n\n## ENDS MULTILEVEL MODEL STUFF\n## -----------------------------------------------------------------------------\n\nlpa_main_fx_c <-\n lpa_t0t1_outc %>%\n split(.$c) %>%\n map(~\n lmer(\n si ~\n pb + tb + bhs + dep +\n ac + cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n map(~ tidy(as(., \"merModLmerTest\"), conf.method = \"profile\", conf.int = TRUE)) %>%\n bind_rows(., .id = \"var\")\n\nlpa_out <-\n bind_rows(\n lpa_out_ni_tidy,\n lpa_main_fx_c\n )\n\nlpa_out_tbl <-\n lpa_out %>%\n select(var, effect, group, term, estimate, std.error, starts_with(\"conf\"), p.value) %>%\n pivot_wider(\n ## -c(effect:group),\n names_from = var,\n values_from = c(estimate, std.error, conf.low, conf.high, p.value)\n ) %>%\n filter(is.na(group)) %>%\n select(term, ends_with(c(\"full\", \"_1\", \"_2\", \"_3\"))) %>%\n clean_names()\n\nlpa_out_tbl %>%\n filter(if_any(starts_with(\"p_value\"), ~.x < 0.05)) %>%\n select(term, starts_with(\"p_value\"))\n\nlpa_main_fx_c %>%\n filter(p.value < 0.05) %>%\n print(n = Inf)\n\nlpa_t0t1_outc %>%\n split(.$c) %>%\n map(~\n lmer(\n si ~\n pb + tb + bhs + dep +\n ac + cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n map(~ stdCoef.merMod(.))\n\n## test for interactions\nlpa_out_i <-\n lmer(\n si ~\n pb + tb + bhs + dep + anx + str +\n ac + cov_sess + age + sex +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n isas_social + isas_intra +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n control = lmerControl(optimizer = \"Nelder_Mead\"),\n REML = FALSE # use ML estimation\n )\n\n## singular fit arises\nall_lpa_i <- allFit(lpa_out_i)\nsummary(all_lpa_i)\n\nanova(\n lpa_out_i,\n lpa_out_step\n)\n## null model with random effects\nlpa_out_null_lmer <-\n lmer(\n si_z ~ 1 +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE\n )\n\n## null model with no random effects\nlpa_out_null_gls <-\n nlme::gls(\n si_z ~ 1,\n data = lpa_t0t1_outc,\n method = \"ML\",\n na.action = na.exclude\n )\n\nrcompanion::nagelkerke(\n lpa_out_step,\n lpa_out_null_lmer,\n )\n\nrcompanion::nagelkerke(\n lpa_out_step,\n lpa_out_null_gls,\n )\n\n.lpa_t_outc <-\n lpa_t0t1_outc %>%\n dplyr::select(si, pb, tb, bhs, dep) %>%\n map( ~ tidy(pairwise.t.test(., lpa_t0t1_outc$fin_ses_c, p.adjust.method = \"bonf\"))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n dplyr::select(pb.group1, pb.group2, ends_with(\"value\")) %>%\n rename(group_1 = 1, group_2 = 2) %>%\n janitor::clean_names()\n.lpa_t_outc\n\n### ============================================================================\n### RELIABLE CHANGE\n### ============================================================================\n\n## RELIABLE CHANGE MSSI CASENESS\n## =============================================================================\n## reliable change calculations\n## pb rci\n\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"pb_0\",\n post = \"pb_1\",\n group = \"c\",\n reliability = meas_rel$pb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\npb_rci_hisi <- JTRCIdf\npb_rci_hisi_fig <-\n plotRCI(\n pb_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Burdensomeness\"\n )\npb_rci_hisi_tbl <-\n tableRCI(\n pb_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"pb\", .before = everything())\n\n\n## tb rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"tb_0\",\n post = \"tb_1\",\n group = \"c\",\n reliability = meas_rel$tb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ntb_rci_hisi <- JTRCIdf\ntb_rci_hisi_fig <-\n plotRCI(\n tb_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Belongingness\"\n )\ntb_rci_hisi_tbl <-\n tableRCI(\n tb_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"tb\", .before = everything())\n\n## bhs rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"bhs_0\",\n post = \"bhs_1\",\n group = \"c\",\n reliability = (meas_rel$bhs_sf_rel$total$raw_alpha + meas_rel$bhs_20_rel$total$raw_alpha) / 2,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nbhs_rci_hisi <- JTRCIdf\nbhs_rci_hisi_fig <-\n plotRCI(\n bhs_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Hopelessness\"\n )\nbhs_rci_hisi_tbl <-\n tableRCI(\n bhs_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"bhs\", .before = everything())\n\n## dep rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"dep_0\",\n post = \"dep_1\",\n group = \"c\",\n reliability = meas_rel$dep_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ndep_rci_hisi <- JTRCIdf\ndep_rci_hisi_fig <-\n plotRCI(\n dep_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Depression\"\n )\ndep_rci_hisi_tbl <-\n tableRCI(\n dep_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"dep\", .before = everything())\n\n## anx rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"anx_0\",\n post = \"anx_1\",\n group = \"c\",\n reliability = meas_rel$anx_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nanx_rci_hisi <- JTRCIdf\nanx_rci_hisi_fig <-\n plotRCI(\n anx_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Anxiety\"\n )\nanx_rci_hisi_tbl <-\n tableRCI(\n anx_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"anx\", .before = everything())\n\n## str rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"str_0\",\n post = \"str_1\",\n group = \"c\",\n reliability = meas_rel$str_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nstr_rci_hisi <- JTRCIdf\nstr_rci_hisi_fig <-\n plotRCI(\n str_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Stress\"\n )\nstr_rci_hisi_tbl <-\n tableRCI(\n str_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"str\", .before = everything())\n\n## si rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"si_0\",\n post = \"si_1\",\n group = \"c\",\n reliability = meas_rel$si_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nsi_rci_hisi <- JTRCIdf\nsi_rci_hisi_fig <-\n plotRCI(\n si_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n plottitle = \"Suicidal Ideation\"\n )\nsi_rci_hisi_tbl <-\n tableRCI(\n si_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"si\", .before = everything())\n\n.rci_hisi_tbl <-\n bind_rows(\n si_rci_hisi_tbl,\n pb_rci_hisi_tbl,\n tb_rci_hisi_tbl,\n bhs_rci_hisi_tbl,\n dep_rci_hisi_tbl\n ) %>%\n rename(rci_hisi = 2) %>%\n mutate(rci_hisi = case_when(\n str_detect(rci_hisi, \"deteriorated\") ~ \"det\",\n str_detect(rci_hisi, \"no reliable\") ~ \"noch\",\n str_detect(rci_hisi, \"improved\") ~ \"imp\"\n ),\n across(\n where(is.character) & -c(var, rci_hisi),\n as.numeric\n ),\n ) %>%\n group_by(var) %>%\n mutate(\n across(\n where(is.numeric),\n ~ round(.x / sum(.x) * 100, 1),\n .names = \"{.col}_pr\"\n )\n ) %>%\n ungroup()\n\n.rci_hisi_tbl_chisq <-\n list(\n si = si_rci_hisi_tbl,\n pb = pb_rci_hisi_tbl,\n tb = tb_rci_hisi_tbl,\n bhs = bhs_rci_hisi_tbl,\n dep = dep_rci_hisi_tbl\n ) %>%\n map(`[`, 3:5) %>%\n map(~ .x %>% mutate(across(where(is.character), as.numeric))) %>%\n map(~ .x %>% stats::chisq.test(table(.))) %>%\n map(`[`, c(\"p.value\", \"statistic\", \"stdres\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ .x %>%\n rename(\n stdres_1 = 3,\n stdres_2 = 4,\n stdres_3 = 5\n )) %>%\n ## map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n mutate(\n across(\n starts_with(\"stdres_\"),\n ~ if_else(abs(.x) > bonf, \"yes\", \"no\"),\n ## .names = \"{.col}_sig\"\n .names = \"{.col}_sig\"\n )\n ) %>%\n mutate(\n rci_hisi =\n rep(c(\"det\", \"noch\", \"imp\"), 5),\n .after = var\n )\n\nrci_hisi_tbl <-\n .rci_hisi_tbl %>%\n full_join(.rci_hisi_tbl_chisq)\n\n.rci_full_kbl <-\n .rci_tbl %>%\n select(var:x3) %>%\n pivot_wider(\n names_from = var,\n values_from = x1:x3\n ) %>%\n ## rownames_to_column() %>%\n pivot_longer(\n -rci,\n 'variable',\n 'value'\n ) %>%\n pivot_wider(variable, rci) %>%\n separate(variable, into = c(\"c\", \"var\"), sep = \"_\") %>%\n arrange(var, c) %>%\n rowwise() %>%\n mutate(\n N = sum(c_across(det:imp)), .after = var,\n across(\n det:imp,\n ~ round(.x / N * 100, 1)\n )\n ) %>%\n ungroup() %>%\n mutate(c = str_remove(c, \"x\")) %>%\n select(var, c, N, imp, noch, det)\n\n\nrci_full_kbl <-\n .rci_tbl_chisq %>%\n group_by(var) %>%\n slice(1) %>%\n ungroup() %>%\n select(var, statistic, p_value) %>%\n right_join(.rci_full_kbl) %>%\n select(var, c:det, statistic, p_value)\n\n.rci_hisi_kbl <-\n .rci_hisi_tbl %>%\n select(var:x3) %>%\n pivot_wider(\n names_from = var,\n values_from = x1:x3\n ) %>%\n ## rownames_to_column() %>%\n pivot_longer(\n -rci_hisi,\n 'variable',\n 'value'\n ) %>%\n pivot_wider(variable, rci_hisi) %>%\n separate(variable, into = c(\"c\", \"var\"), sep = \"_\") %>%\n arrange(var, c) %>%\n rowwise() %>%\n mutate(\n N = sum(c_across(det:imp)), .after = var,\n across(\n det:imp,\n ~ round(.x / N * 100, 1)\n )\n ) %>%\n ungroup() %>%\n mutate(c = str_remove(c, \"x\")) %>%\n select(var, c, N, imp, noch, det)\n\nrci_hisi_kbl <-\n .rci_hisi_tbl_chisq %>%\n group_by(var) %>%\n slice(1) %>%\n ungroup() %>%\n select(var, statistic, p_value) %>%\n right_join(.rci_hisi_kbl) %>%\n select(var, c:det, statistic, p_value)\n\n\n.rci_kbl_group_diff <-\n readxl::read_xlsx(\n \"./data/tmp/anova_results.xlsx\",\n sheet = \"rci_tbl\"\n ) %>%\n select(sample:rci, stdres_1_sig:group_diff) %>%\n pivot_longer(\n -c(sample, rci, var, group_diff)\n ) %>%\n group_by(sample, var, rci, name) %>%\n slice(1) %>%\n ungroup() %>%\n mutate(\n c = str_extract(name, \"\\\\d\")\n ) %>%\n select(-name) %>%\n pivot_wider(\n names_from = c(sample, rci),\n values_from = group_diff\n ) %>%\n unite(\"group_diff_full\", c(full_imp, full_noch, full_det), na.rm = TRUE, sep = \";<br>\") %>%\n unite(\"group_diff_hisi\", c(hisi_imp, hisi_noch, hisi_det), na.rm = TRUE, sep = \";<br>\") %>%\n mutate(\n across(\n starts_with(\"group_\"),\n ~ if_else(.x == \"\", NA_character_, .x)\n )\n ) %>%\n select(-value) %>%\n group_by(var, c) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup()\n\n\nrci_kbl <-\n list(\n full = rci_full_kbl,\n hisi = rci_hisi_kbl\n ) %>%\n bind_rows(.id = \"vers\") %>%\n pivot_wider(\n names_from = vers,\n values_from = N:p_value\n ) %>%\n full_join(.rci_kbl_group_diff) %>%\n select(\n var, c,\n ends_with(c(\"full\", \"hisi\"))\n ) %>%\n mutate(\n order = case_when(\n var == \"si\" ~ 1,\n var == \"pb\" ~ 2,\n var == \"tb\" ~ 3,\n var == \"bhs\" ~ 4,\n var == \"dep\" ~ 5\n ), .before = 1\n ) %>%\n arrange(order)\n\n\nrci_hisi_combi_plot <-\n pb_rci_hisi_fig +\n tb_rci_hisi_fig +\n bhs_rci_hisi_fig +\n dep_rci_hisi_fig +\n anx_rci_hisi_fig +\n str_rci_hisi_fig +\n si_rci_hisi_fig +\n guide_area() +\n plot_layout(\n nrow = 2,\n guides = \"collect\",\n ) &\n theme(\n plot.title = element_text(size = 10),\n legend.title = element_text(size = 10),\n legend.text = element_text(size = 9),\n axis.title.x = element_text(size = 9),\n axis.title.y = element_text(size = 9),\n axis.text = element_text(size = 7),\n ) &\n labs(\n colour = \"RCI_HISI classification\",\n shape = \"Class\"\n ) &\n scale_fill_brewer(palette = \"Dark2\") &\n scale_colour_brewer(palette = \"Dark2\")\n\nggsave(\n \"./figures/rci_hisi_combi_plot.png\",\n rci_hisi_combi_plot\n)\n\n\n\n## RELIABLE CHANGE\n## =============================================================================\n## reliable change calculations\n## pb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"pb_0\",\n post = \"pb_1\",\n group = \"c\",\n reliability = meas_rel$pb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\npb_rci <- JTRCIdf\npb_rci_fig <-\n plotRCI(\n pb_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Burdensomeness\"\n )\npb_rci_tbl <-\n tableRCI(\n pb_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"pb\", .before = everything())\n\n\n## tb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"tb_0\",\n post = \"tb_1\",\n group = \"c\",\n reliability = meas_rel$tb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ntb_rci <- JTRCIdf\ntb_rci_fig <-\n plotRCI(\n tb_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Belongingness\"\n )\ntb_rci_tbl <-\n tableRCI(\n tb_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"tb\", .before = everything())\n\n## bhs rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"bhs_0\",\n post = \"bhs_1\",\n group = \"c\",\n reliability = (meas_rel$bhs_sf_rel$total$raw_alpha + meas_rel$bhs_20_rel$total$raw_alpha) / 2,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nbhs_rci <- JTRCIdf\nbhs_rci_fig <-\n plotRCI(\n bhs_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Hopelessness\"\n )\nbhs_rci_tbl <-\n tableRCI(\n bhs_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"bhs\", .before = everything())\n\n## dep rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"dep_0\",\n post = \"dep_1\",\n group = \"c\",\n reliability = meas_rel$dep_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ndep_rci <- JTRCIdf\ndep_rci_fig <-\n plotRCI(\n dep_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Depression\"\n )\ndep_rci_tbl <-\n tableRCI(\n dep_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"dep\", .before = everything())\n\n## anx rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"anx_0\",\n post = \"anx_1\",\n group = \"c\",\n reliability = meas_rel$anx_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nanx_rci <- JTRCIdf\nanx_rci_fig <-\n plotRCI(\n anx_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Anxiety\"\n )\nanx_rci_tbl <-\n tableRCI(\n anx_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"anx\", .before = everything())\n\n## str rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"str_0\",\n post = \"str_1\",\n group = \"c\",\n reliability = meas_rel$str_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nstr_rci <- JTRCIdf\nstr_rci_fig <-\n plotRCI(\n str_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Stress\"\n )\nstr_rci_tbl <-\n tableRCI(\n str_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"str\", .before = everything())\n\n## si rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"si_0\",\n post = \"si_1\",\n group = \"c\",\n reliability = meas_rel$si_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nsi_rci <- JTRCIdf\nsi_rci_fig <-\n plotRCI(\n si_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n plottitle = \"Suicidal Ideation\"\n )\nsi_rci_tbl <-\n tableRCI(\n si_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"si\", .before = everything())\n\n.rci_tbl <-\n bind_rows(\n si_rci_tbl,\n pb_rci_tbl,\n tb_rci_tbl,\n bhs_rci_tbl,\n dep_rci_tbl\n ) %>%\n rename(rci = 2) %>%\n mutate(rci = case_when(\n str_detect(rci, \"deteriorated\") ~ \"det\",\n str_detect(rci, \"no reliable\") ~ \"noch\",\n str_detect(rci, \"improved\") ~ \"imp\"\n ),\n across(\n where(is.character) & -c(var, rci),\n as.numeric\n ),\n ) %>%\n group_by(var) %>%\n mutate(\n across(\n where(is.numeric),\n ~ round(.x / sum(.x) * 100, 1),\n .names = \"{.col}_pr\"\n )\n ) %>%\n ungroup()\n\n.rci_tbl_chisq <-\n list(\n si = si_rci_tbl,\n pb = pb_rci_tbl,\n tb = tb_rci_tbl,\n bhs = bhs_rci_tbl,\n dep = dep_rci_tbl\n ) %>%\n map(`[`, 3:5) %>%\n map(~ .x %>% mutate(across(where(is.character), as.numeric))) %>%\n map(~ .x %>% stats::chisq.test(table(.))) %>%\n map(`[`, c(\"p.value\", \"statistic\", \"stdres\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ .x %>%\n rename(\n stdres_1 = 3,\n stdres_2 = 4,\n stdres_3 = 5\n )) %>%\n ## map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n mutate(\n across(\n starts_with(\"stdres_\"),\n ~ if_else(abs(.x) > bonf, \"yes\", \"no\"),\n ## .names = \"{.col}_sig\"\n .names = \"{.col}_sig\"\n )\n ) %>%\n mutate(\n rci =\n rep(c(\"det\", \"noch\", \"imp\"), 5),\n .after = var\n )\n\nrci_tbl <-\n .rci_tbl %>%\n full_join(.rci_tbl_chisq)\n\nrci_combi_plot <-\n pb_rci_fig +\n tb_rci_fig +\n bhs_rci_fig +\n dep_rci_fig +\n anx_rci_fig +\n str_rci_fig +\n si_rci_fig +\n guide_area() +\n plot_layout(\n nrow = 2,\n guides = \"collect\",\n ) &\n theme(\n plot.title = element_text(size = 10),\n legend.title = element_text(size = 10),\n legend.text = element_text(size = 9),\n axis.title.x = element_text(size = 9),\n axis.title.y = element_text(size = 9),\n axis.text = element_text(size = 7),\n ) &\n labs(\n colour = \"RCI classification\",\n shape = \"Class\"\n ) &\n scale_fill_brewer(palette = \"Dark2\") &\n scale_colour_brewer(palette = \"Dark2\")\n\nggsave(\n \"./figures/rci_combi_plot.png\",\n rci_combi_plot\n)\n\n\n## JT CLIN SIG CH\n## =============================================================================\n## reliable change calculations\n## pb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"pb_0\",\n post = \"pb_1\",\n group = \"c\",\n reliability = meas_rel$pb_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\npb_jt <- JTRCIdf\npb_jt_fig <-\n plotJT(\n pb_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Burdensomeness\"\n )\npb_jt_tbl <-\n tableJT(\n pb_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"pb\", .before = everything())\n\n## tb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"tb_0\",\n post = \"tb_1\",\n group = \"c\",\n reliability = meas_rel$tb_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ntb_jt <- JTRCIdf\ntb_jt_fig <-\n plotJT(\n tb_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Belongingness\"\n )\ntb_jt_tbl <-\n tableJT(\n tb_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"tb\", .before = everything())\n\n## bhs rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"bhs_0\",\n post = \"bhs_1\",\n group = \"c\",\n reliability = (meas_rel$bhs_sf_rel$total$raw_alpha + meas_rel$bhs_20_rel$total$raw_alpha) / 2,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nbhs_jt <- JTRCIdf\nbhs_jt_fig <-\n plotJT(\n bhs_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Hopelessness\"\n )\nbhs_jt_tbl <-\n tableJT(\n bhs_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"bhs\", .before = everything())\n\n## dep rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"dep_0\",\n post = \"dep_1\",\n group = \"c\",\n reliability = meas_rel$dep_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ndep_jt <- JTRCIdf\ndep_jt_fig <-\n plotJT(\n dep_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Depression\"\n )\ndep_jt_tbl <-\n tableJT(\n dep_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"dep\", .before = everything())\n\n## anx rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"anx_0\",\n post = \"anx_1\",\n group = \"c\",\n reliability = meas_rel$anx_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nanx_jt <- JTRCIdf\nanx_jt_fig <-\n plotJT(\n anx_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Anxiety\"\n )\nanx_jt_tbl <-\n tableJT(\n anx_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"anx\", .before = everything())\n\n## str rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"str_0\",\n post = \"str_1\",\n group = \"c\",\n reliability = meas_rel$str_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nstr_jt <- JTRCIdf\nstr_jt_fig <-\n plotJT(\n str_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Stress\"\n )\nstr_jt_tbl <-\n tableJT(\n str_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"str\", .before = everything())\n\n## si rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"si_0\",\n post = \"si_1\",\n group = \"c\",\n reliability = meas_rel$si_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nsi_jt <- JTRCIdf\nsi_jt_fig <-\n plotJT(\n si_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n plottitle = \"Suicidal Ideation\"\n )\nsi_jt_tbl <-\n tableJT(\n si_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"si\", .before = everything())\n\njt_tbl <-\n pb_jt_tbl %>%\n full_join(tb_jt_tbl) %>%\n full_join(bhs_jt_tbl) %>%\n full_join(dep_jt_tbl) %>%\n full_join(anx_jt_tbl) %>%\n full_join(str_jt_tbl) %>%\n full_join(si_jt_tbl)\n\n\n.jt_tbl <-\n bind_rows(\n si_jt_tbl,\n pb_jt_tbl,\n tb_jt_tbl,\n bhs_jt_tbl,\n dep_jt_tbl\n ) %>%\n rename(jt = 2) %>%\n mutate(jt = case_when(\n str_detect(jt, \"deteriorated\") ~ \"det\",\n str_detect(jt, \"unchanged\") ~ \"unch\",\n str_detect(jt, \"improved\") ~ \"imp\",\n str_detect(jt, \"non reliab\") ~ \"nr_rec\",\n TRUE ~ \"rec\"\n ),\n across(\n where(is.character) & -c(var, jt),\n as.numeric\n ),\n ) %>%\n group_by(var) %>%\n mutate(\n across(\n where(is.numeric),\n ~ round(.x / sum(.x) * 100, 1),\n .names = \"{.col}_pr\"\n )\n ) %>%\n ungroup()\n\n.jt_tbl_chisq <-\n list(\n si = si_jt_tbl,\n pb = pb_jt_tbl,\n tb = tb_jt_tbl,\n bhs = bhs_jt_tbl,\n dep = dep_jt_tbl\n ) %>%\n map(`[`, 3:5) %>%\n map(~ .x %>% mutate(across(where(is.character), as.numeric))) %>%\n map(~ .x %>% stats::chisq.test(table(.))) %>%\n map(`[`, c(\"p.value\", \"statistic\", \"stdres\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ .x %>%\n rename(\n stdres_1 = 3,\n stdres_2 = 4,\n stdres_3 = 5\n )) %>%\n ## map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n mutate(\n across(\n starts_with(\"stdres_\"),\n ~ if_else(abs(.x) > bonf, \"yes\", \"no\"),\n ## .names = \"{.col}_sig\"\n .names = \"{.col}_sig\"\n )\n ) %>%\n mutate(\n jt =\n rep(c(\"det\", \"unch\", \"imp\", \"nr_rec\", \"rec\"), 5),\n .after = var\n )\n\njt_tbl <-\n .jt_tbl %>%\n full_join(.jt_tbl_chisq)\n\n\njt_combi_plot <-\n pb_jt_fig +\n tb_jt_fig +\n bhs_jt_fig +\n dep_jt_fig +\n anx_jt_fig +\n str_jt_fig +\n si_jt_fig +\n guide_area() +\n plot_layout(\n nrow = 2,\n guides = \"collect\",\n ) &\n theme(\n plot.title = element_text(size = 10),\n legend.title = element_text(size = 10),\n legend.text = element_text(size = 9),\n axis.title.x = element_text(size = 9),\n axis.title.y = element_text(size = 9),\n axis.text = element_text(size = 7),\n ) &\n labs(\n colour = \"JT Classification\",\n shape = \"Class\"\n ) &\n scale_fill_brewer(palette = \"Dark2\") &\n scale_colour_brewer(palette = \"Dark2\")\n\nggsave(\n \"./figures/jt_combi_plot.png\",\n jt_combi_plot\n)\n\n### ============================================================================\n### ABSTRACT DF\n### ============================================================================\n\nabst_df <-\n lpa_t0 %>%\n summarise(\n N = sum(!is.na(c)),\n age_m = round(mean(age, na.rm = TRUE), 2),\n age_sd = round(sd(age, na.rm = TRUE), 2),\n prc_1 = round(sum(!is.na(c[c == 1])) / N * 100, 2),\n prc_2 = round(sum(!is.na(c[c == 2])) / N * 100, 2),\n prc_3 = round(sum(!is.na(c[c == 3])) / N * 100, 2),\n )\n\nlpa_t0 %>%\n group_by(c, protocol) %>%\n summarise(\n N = n()\n )\n\n### ============================================================================\n### END SCRIPT\n### ============================================================================\n\n## CORRELATION TABLES\n## =============================================================================\ncorr_des <-\n lpa_t0t1_outc %>%\n filter(final_session == 0) %>%\n select(\n age, sex, hxsa, hxdsh, hxsuicide, hxtrauma,\n pb, si, tb, dep, anx, str, ac,\n isas_social, isas_intra,\n ) %>%\n mutate(across(everything(), as.numeric))\n\ncts <-\n corr_des %>%\n ## mutate_all(as.numeric) %>%\n corr.test()\n\ncts_p <-\n print(corr.p(cts$r, n = 500), short = FALSE) %>%\n ## as_tibble()\n arrange(p, desc(r))\n\ncts_p %>%\n filter(\n p <= 0.05\n )\n\n## correlation plots\n## =============================================================================\ncorr_des %>%\n cor.ci() %>%\n cor.plot.upperLowerCi()\n\n## histograms\n## =============================================================================\ntest_des %>%\n ## colnames() %>%\n mutate_all(as.numeric) %>%\n select(\n age, sex, suicide_attempts,\n contains(\"hx\"),\n ends_with(\"mean\")\n ) %>%\n ## review bindwidth\n multi.hist()\n\n### ============================================================================\n### PARTCIPANTS FOR METHODS SECTION\n### ============================================================================\n## consenters\nlpa_consenters <-\n read_csv(here(.phd_sp_full_df)) %>%\n filter(session %in% 1, sp_episode %in% 1) %>%\n group_by(consent) %>%\n summarise(\n n = n()\n )\n\n## MULTIPLE EPISODES\n## -----------------------------------------------------------------------------\n## extract the uci of those used for the lpa\n.lpa_target <-\n lpa_fin %>%\n filter(!is.na(c)) %>%\n select(uci) %>%\n distinct(uci) %>%\n deframe()\n\n## get the uci and class id\n.lpa_c_uci <-\n lpa_fin %>%\n filter(!is.na(c)) %>%\n select(uci, c) %>%\n group_by(uci) %>%\n slice(1) %>%\n ungroup()\n\n## get the number of episodes by class\n.lpa_sp_eps <-\n read_csv(here(.phd_sp_full_df)) %>%\n full_join(.lpa_c_uci, by = \"uci\") %>%\n select(c, uci, sp_episode, session) %>%\n filter(!is.na(c)) %>%\n group_by(uci, sp_episode) %>%\n slice(1) %>%\n group_by(uci) %>%\n filter(sp_episode == max(sp_episode)) %>%\n group_by(sp_episode, c) %>%\n summarise(\n n = n(),\n ) %>%\n group_by(c) %>%\n mutate(prc_ep = paste0(round(n/sum(n)*100,1), \"%\"))\n\n## SEX AND AGE DIFFERENCE\n## -----------------------------------------------------------------------------\n## sex diff with age\n.lpa_sex_diff <-\n lpa_fin %>%\n filter(!is.na(c)) %>%\n group_by(ID) %>%\n slice(1) %>%\n ungroup() %>%\n select(ID, c, sex, age)\n\nAnova(lm(age ~ sex, data = .lpa_sex_diff , na.action = na.exclude), type = 3) %>%\n tidy()\n\n## parametric\nAnova(lm(age ~ sex, data = .lpa_sex_diff, na.action = na.exclude), type = 3) %>%\n tidy() %>%\n filter(term %in% \"sex\") %>%\n select(term, p.value)\n\n.lpa_sex_diff %>%\n select(age) %>%\n map( ~ tidy(pairwise.t.test(., .lpa_sex_diff$sex, p.adjust.method = \"bonf\")))\n\n## nonparametric\n.lpa_sex_diff %>%\n select(age) %>%\n map( ~ tidy(pairwise.wilcox.test(., .lpa_sex_diff$sex, p.adjust.method = \"bonf\")))\n\n.lpa_sex_diff %>%\n split(.$c) %>%\n map(\n ~ lm(age ~ sex, data = .x, na.action = na.exclude) %>%\n summary() %>%\n tidy() %>%\n filter(term %in% \"sex\") %>%\n select(term, p.value)\n ) %>%\n bind_rows(.id = \"c\")\n\nlpa_sex_age_diff <-\n .lpa_sex_diff %>%\n group_by(sex) %>%\n summarise(\n n = n(),\n across(\n age,\n list(\n mean = ~ mean(.x, na.rm = TRUE),\n median = ~ median(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n min = ~ min(.x, na.rm = TRUE),\n max = ~ max(.x, na.rm = TRUE)\n ))\n ) %>%\n filter(sex %in% c(0:2)) %>%\n mutate(\n across(c(age_mean, age_sd), ~ sprintf(\"%.2f\", .x))\n )\n\n##\n## PLANNED FINAL DIFFERENCES\n## -----------------------------------------------------------------------------\n## planned final\nlpa_planned_fin_df <-\n lpa_t0t1_outc %>%\n select(c:planned_final, si, pb:dep, -fin_ses_c, -session) %>%\n pivot_wider(\n names_from = final_session,\n values_from = c(planned_final:dep)\n ) %>%\n ## filter(final_session %in% 1) %>%\n group_by(planned_final_1, c) %>%\n summarise(\n n = n(),\n across(\n starts_with(c(\"pb\", \"ac\", \"tb\", \"bhs\", \"dep\", \"si\", \"age\", \"cov\")),\n list(\n mean = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n median = ~ median(.x, na.rm = TRUE)\n )\n )) %>%\n group_by(c) %>%\n mutate(\n prc = paste0(round(n/sum(n)*100,1),\"%\"), .after = n,\n across(is.numeric, ~ round(.x, 2))\n ) %>%\n arrange(c)\n\n\nlpa_dropout <-\n lmer(\n cov_sess ~\n ## predictors scaled due to lmer warning\n planned_final +\n (1 | c),\n data = lpa_t0t1_outc[lpa_t0t1_outc$final_session %in% 1, ],\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(as(lpa_dropout, \"merModLmerTest\"))\n\nlpa_t0t1_outc <-\n lpa_t0t1_outc %>%\n group_by(ID) %>%\n mutate(\n planned_final = max(as.numeric(as.character(planned_final))),\n ) %>%\n ungroup()\n\n\nlm(\n dep ~\n planned_final,\n data = lpa_t0t1_outc[lpa_t0t1_outc$c %in% 3 & lpa_t0t1_outc$final_session %in% 0, ],\n na.action = na.exclude\n) %>%\n summary()\n\n.lpa_plan_fin <-\n lpa_fin %>%\n filter(!is.na(c), final_session %in% 1) %>%\n group_by(c, planned_final) %>%\n summarise(\n n = n()\n ) %>%\n group_by(c) %>%\n mutate(prc_fin = paste0(round(n/sum(n)*100,1), \"%\"))\n\n## is there a difference in dropout by class?\nlpa_fin %>%\n filter(!is.na(c), final_session %in% 1) %>%\n ## glm(planned_final ~ c, data = ., na.action = na.exclude) %>%\n lm(planned_final ~ c, data = ., na.action = na.exclude) %>%\n summary() %>%\n tidy()\n\n## regression with palnned final\nlpa_out_ipt_r2 <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex + planned_final +\n ## cov_sess + age + sex +\n ## (cov_sess + planned_final | final_session/c),\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n ## control = lmerControl(optimizer = \"Nelder_Mead\"),\n REML = FALSE # use ML estimation\n )\n\nall_lpa_r2 <- allFit(lpa_out_ipt_r2)\n\nsummary(all_lpa_r2)\n\nsummary(as(lpa_out_ipt_r2, \"merModLmerTest\"))\n\nlpa_out_ipt_r2_std <-\n lpa_out_ipt_r2 %>%\n stdCoef.merMod()\n\nlpa_out_ipt_r2_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ipt_r2, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goes with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_out_ipt_r2_tidy %>%\n filter(p.value < 0.05) %>%\n print(n = Inf)\n\n### ============================================================================\n### WRITE IMAGE AND DATA\n### ============================================================================\n\n## TODO remove P from env\nsave.image(\"./data/proc/draft.RData\")\n\nrep_df <-\n mget(ls())\n\nsaveRDS(\n rep_df,\n \"./data/proc/lpa_df.rds\"\n)\n\n### ============================================================================\n### ENDS\n### ============================================================================\n"
}
}
[Trace - 03:06:25 PM] Received response 'textDocument/codeAction - (379)' in 120ms.
Result: []
[Trace - 03:06:25 PM] Received notification 'textDocument/publishDiagnostics'.
Params: {
"diagnostics": [],
"version": null,
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
}
[Trace - 03:06:26 PM] Sending request 'textDocument/documentSymbol - (381)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
}
}
[Trace - 03:06:26 PM] Sending notification 'textDocument/didChange'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R",
"version": 1
},
"contentChanges": [
{
"text": "### ============================================================================\n### SOURCE SCRIPTS\n### ============================================================================\nrm(list = ls(all = TRUE))\nsource(\"./code/func.R\")\nsource(\"./code/load.R\")\n\nload(\"./data/proc/draft.RData\")\nls(all = TRUE)\n\n## save.image(\"./data/tmp/draft.RData\")\n\n### ============================================================================\n### LOAD DATA\n### ============================================================================\n\n\n## mplus missing data analysis\nlpa_miss <-\n readRDS(here(.m1_optseed))$c0_optseed_m_v1.out\n\n## full lpa output\nll_tests_m <-\n readRDS(here(.m2_ll_tests))\n\n## first split half lpa output\nll_tests_m_rv <-\n readRDS(here(.m2_ll_tests_rv))\n\n## three class lpa output\nlpa_c3 <-\n ll_tests_m$c3_2_ll_tests_m_v1.out\n\n## final data set with three class lpa profiles added\nlpa_fin <-\n readRDS(here(.lpa_fin))\n\nlpa_fin %>%\n select(session, final_session, pb, tb, bhs, ac, dep, anx, c, si) %>%\n filter(session == 1 | final_session == 1) %>%\n group_by(final_session) %>%\n .summary_by(c)\n\nlpa_fin %>%\n select(ID, c, session, final_session) %>%\n filter(session == 1) %>%\n describeBy(\"c\")\n\nlpa_fin %>%\n select(ID, c, session, final_session) %>%\n filter(session != 1 & final_session == 1) %>%\n describeBy(\"c\")\n\nlpa_fin %>%\n group_by(c) %>%\n summarise(\n n_final = length(final_session[which(final_session == 1 & session != 1)]),\n n_init = length(session[session %in% 1]),\n ) %>%\n filter(!is.na(c)) %>%\n mutate(\n perc_lost = round((n_init - n_final) / n_init * 100, digits = 1),\n c = as.character(c)\n ) %>%\n janitor::adorn_totals()\n\n## CREATE MI DATASETS\n## =============================================================================\n## SUBSET TIME 1: DEMO\n## -----------------------------------------------------------------------------\nlpa_t0 <-\n lpa_fin %>%\n filter(!is.na(c), session == 1 | final_session == 1) %>%\n group_by(ID) %>%\n mutate(\n total_session = max(session), .after = session,\n planned_final = max(planned_final)\n ) %>%\n ungroup() %>%\n filter(session == 1) %>%\n select(ID:suburb)\n\nlpa_t0 %>% names()\nlpa_t0 %>% summary()\n\n## SUBSET TIME 1: MI\n## -----------------------------------------------------------------------------\nlpa_t0_mi <-\n lpa_fin %>%\n filter(session == 1 & !is.na(c)) %>%\n mutate(across(c(c:sex, atsi:closehxdsh, suicide_attempts), as.factor)) %>%\n select(c, ID, pb, ac, bhs, tb, age, sex, dep, anx, si)\n\nlpa_t0_mi %>% names()\nlpa_t0_mi %>% .summary_by(c)\n\n## SUBSET TIME 1 WITH MATCHING TIME 2\n## -----------------------------------------------------------------------------\nlpa_t0t1 <-\n lpa_fin %>%\n ## filter no class var\n filter(!is.na(c)) %>%\n group_by(ID) %>%\n ## filter first and last session\n filter(session == 1 | final_session == 1) %>%\n ungroup() %>%\n mutate(\n ## replace NA session value\n session = if_else(is.na(session), -999, session),\n ## create conditional variable to drop those with only one session\n drop = if_else(\n session == 1 & final_session == 1,\n 1, 0\n ), .after = ID\n ) %>%\n ## select those with more than one session\n filter(drop == 0) %>%\n ## select vars for mi\n select(\n c, ID, age, sex, session, final_session,\n planned_final, pb, ac, bhs, tb, dep, anx, si\n )\n\nlpa_t0t1_mi <-\n lpa_t0t1 %>%\n panelr::panel_data(\n id = ID,\n wave = final_session,\n ) %>%\n panelr::widen_panel(\n separator = \"_\",\n ignore.attributes = FALSE,\n varying = NULL) %>%\n ungroup() %>%\n ## need to convert to data frame as panelr creates a list object\n as.data.frame() %>%\n ## use age at intake, drop secon ac as not routinely collected\n select(-age_1, -ac_1) %>%\n rename(age = age_0) %>%\n mutate(\n ## remove -999 placehoder session number\n session_1 = if_else(session_1 == -999, NA_real_, session_1)\n )\n\nlpa_t0t1_mi %>% str()\nlpa_t0t1_mi %>% summary()\n\n### ============================================================================\n### MICE MI\n### ============================================================================\n## TODO explore missingness overall\n## TODO explore missingness by profile\n## NOTE need to as a minimum report missingness percentages\n## Nice to have would be to determine if MAR or MNAR instead of just assuming\n## TODO subset vars\n## DONE get time 1 and final session data and bring into imputation model\n## DONE decide if can also bring in the demographic variables to be imputed\n## NOTE Frank suggested this was unusual and so was not done\n## https://www.gerkovink.com/miceVignettes/Multi_level/Multi_level_data.html\n\n## TIME 1\n## =============================================================================\n.summary_by(lpa_t0_mi, c)\n\nmd.pattern(lpa_t0_mi)\n\n.miss_df(lpa_t0_mi)\n\nfluxplot(lpa_t0_mi)\n## NOTE conisder historgrams inspecting missingness here too\n\n## GET AND SET PREDICTOR VARS\n## -----------------------------------------------------------------------------\n## NOTE this was used to create an excel for inspection and manipulation\npred_test <-\n quickpred(lpa_t0_mi) %>%\n .mice_pred_df()\n\n.pred_ <-\n quickpred(lpa_t0_mi)\n.pred_\n\n## import pred via excel manipulation\npred_ <-\n readxl::read_xlsx(\n here(\"data/tmp/mi_working.xlsx\"),\n sheet = \"pred_matrix\",\n range = \"A1:L12\"\n ) %>%\n .mice_pred_matrix()\npred_\n\n## comfirm they are different\nidentical(pred_, .pred_)\n\n## create pred here without exporting to excel\npred <- quickpred(lpa_t0_mi)\n## c\n## pred[,\"c\"] <- c(0, 0, seq(-2, -2, length.out = 9))\npred[, \"c\"] <- c(0, 0, seq(1, 1, length.out = 9))\n## pred[,\"c\"] <-\n## ifelse(pred[,\"c\"] == 0, 0, -2)\n## pred[\"age\", \"c\"] <- 0\npred[, \"ID\"] <- ifelse(pred[,\"ID\"] == 0, 0, 0)\n## pred[,c(\"ID\", \"ac\", \"anx\", \"dep\", \"si\", \"sex\", \"age\")] <-\n## ifelse(pred[,c(\"ID\", \"ac\", \"anx\", \"dep\", \"si\", \"sex\", \"age\")] == 0, 0, 0)\npred\n\n## comfirm they are different\nidentical(pred, pred_)\n\n## IMPUTE DATA\n## =============================================================================\nlpa_t0_mi <-\n lpa_t0_mi %>%\n mutate(c = as.integer(c))\n\n## init mi model\nimp_t0 <-\n mice(\n lpa_t0_mi,\n pred = pred,\n meth = \"pmm\",\n maxit = 0,\n m = 10\n )\n\nimp_t0_5 <-\n mice.mids(\n imp_t0,\n maxit = 5\n )\n\nimp_t0_15 <-\n mice.mids(\n imp_t0_5,\n maxit = 10\n )\n\nimp_t0_25 <-\n mice.mids(\n imp_t0_15,\n maxit = 10\n )\n\nimp_t0_35 <-\n mice.mids(\n imp_t0_25,\n maxit = 10\n )\n\nimp_t0_50 <-\n mice.mids(\n imp_t0_35,\n maxit = 15\n )\n\nimp_t0_65 <-\n mice.mids(\n imp_t0_50,\n maxit = 15\n )\n\nimp_t0_80 <-\n mice.mids(\n imp_t0_65,\n maxit = 15\n )\n\nimp_t0_100 <-\n mice.mids(\n imp_t0_80,\n maxit = 20\n )\n\nimp_t0_115 <-\n mice.mids(\n imp_t0_100,\n maxit = 15\n )\n\nimp_t0_150 <-\n mice.mids(\n imp_t0_115,\n maxit = 35\n )\n\nimp_t0_200 <-\n mice.mids(\n imp_t0_150,\n maxit = 50\n )\n\nimp_t0_300 <-\n mice.mids(\n imp_t0_200,\n maxit = 100\n )\n\nimp_t0_400 <-\n mice.mids(\n imp_t0_300,\n maxit = 100\n )\n\n## inspect each iterations convergence\nsummary(lpa_t0_mi)\nsummary(complete(imp_t0_5))\nsummary(complete(imp_t0_15))\nsummary(complete(imp_t0_25))\nsummary(complete(imp_t0_35))\nsummary(complete(imp_t0_50))\nsummary(complete(imp_t0_65))\nsummary(complete(imp_t0_80))\n\nsummary(complete(imp_t0_100))\nsummary(complete(imp_t0_115))\nsummary(complete(imp_t0_150))\nsummary(complete(imp_t0_200))\nsummary(complete(imp_t0_300))\nsummary(complete(imp_t0_400))\n.summary_by(complete(imp_t0_400), c)\n\nplot(imp_t0_5)\nplot(imp_t0_15)\nplot(imp_t0_25)\nplot(imp_t0_35)\nplot(imp_t0_50)\nplot(imp_t0_65) # this looks okay, perhaps use this one\nplot(imp_t0_80) # this seems to be getting worse\nplot(imp_t0_100) # not too bad\nplot(imp_t0_115) # okay\nplot(imp_t0_150)\nplot(imp_t0_200)\nplot(imp_t0_300)\nplot(imp_t0_400)\n\ndensityplot(imp_t0_5)\ndensityplot(imp_t0_15)\ndensityplot(imp_t0_25)\ndensityplot(imp_t0_35)\ndensityplot(imp_t0_50)\ndensityplot(imp_t0_65)\ndensityplot(imp_t0_80)\ndensityplot(imp_t0_100)\ndensityplot(imp_t0_115)\ndensityplot(imp_t0_150)\ndensityplot(imp_t0_200)\ndensityplot(imp_t0_300)\ndensityplot(imp_t0_400)\n\n\n## NOTE no perfect solution, so will chose to use the imp_t0_400 mids as this\n## has the most iterations and is assumed to be the best as a result.\n\n## TIME 1 AND TIME 2\n## =============================================================================\n## GET AND SET PREDICTOR VARS\n## -----------------------------------------------------------------------------\n## NOTE this was used to create an excel for inspection and manipulation\n## NOTE this was guided by the longitudinal imputation based on the firewords\n## disater in the mice textbook\npred_test_fu <-\n quickpred(lpa_t0t1_mi) %>%\n .mice_pred_df()\npred_test_fu\n\n.pred_fu_ <-\n quickpred(lpa_t0t1_mi)\n.pred_fu_\n\n## import pred via excel manipulation\npred_fu_ <-\n readxl::read_xlsx(\n here(\"data/tmp/mi_working.xlsx\"),\n sheet = \"pred_matrix_fu_simp\",\n range = \"A1:V22\"\n ) %>%\n .mice_pred_matrix()\npred_fu_\n\n## comfirm they are different\nidentical(pred_fu_, .pred_)\n\n## IMPUTE DATA\n## =============================================================================\nlpa_t0t1_mi <-\n lpa_t0t1_mi %>%\n mutate(c = as.integer(c))\n\n## init mi model\nimp_t0t1 <-\n mice(\n lpa_t0t1_mi,\n pred = pred_fu_,\n meth = \"pmm\",\n maxit = 0,\n m = 10\n )\n\nimp_t0t1_5 <-\n mice.mids(\n imp_t0t1,\n maxit = 5\n )\n\nimp_t0t1_15 <-\n mice.mids(\n imp_t0t1_5,\n maxit = 10\n )\n\nimp_t0t1_25 <-\n mice.mids(\n imp_t0t1_15,\n maxit = 10\n )\n\nimp_t0t1_35 <-\n mice.mids(\n imp_t0t1_25,\n maxit = 10\n )\n\nimp_t0t1_50 <-\n mice.mids(\n imp_t0t1_35,\n maxit = 15\n )\n\nimp_t0t1_100 <-\n mice.mids(\n imp_t0t1_50,\n maxit = 50\n )\n\nimp_t0t1_115 <-\n mice.mids(\n imp_t0t1_100,\n maxit = 15\n )\n\n## inspect each iterations convergence\n.summary_by(lpa_t0t1_mi, c)\n.summary_by(complete(imp_t0t1_5), c)\n.summary_by(complete(imp_t0t1_15), c)\n.summary_by(complete(imp_t0t1_25), c)\n.summary_by(complete(imp_t0t1_35), c)\n.summary_by(complete(imp_t0t1_50), c)\n.summary_by(complete(imp_t0t1_100), c)\n.summary_by(complete(imp_t0t1_115), c)\n\n\nplot(imp_t0t1_5)\nplot(imp_t0t1_15)\nplot(imp_t0t1_25)\nplot(imp_t0t1_35) # potential candidate\nplot(imp_t0t1_50) # pretty good with exception to tb_0\nplot(imp_t0t1_100) # not much of an improvement on the above\nplot(imp_t0t1_115) # convergence warnings, will use 50\n\ndensityplot(imp_t0t1_5)\ndensityplot(imp_t0t1_15)\ndensityplot(imp_t0t1_25)\ndensityplot(imp_t0t1_35)\ndensityplot(imp_t0t1_50)\ndensityplot(imp_t0t1_100)\ndensityplot(imp_t0t1_115)\n\n### ============================================================================\n### ANALYSESES\n### ============================================================================\n\n## LPA CLASS TABL\n## =============================================================================\n\nlpa_class_tab <-\n SummaryTable(\n ll_tests_m,\n keepCols = c(\n \"Title\",\n \"LL\",\n \"Entropy\",\n \"BIC\",\n \"T11_KM1LL\",\n \"T11_VLMR_PValue\"\n ),\n sortBy = \"Title\"\n ) %>%\n mutate(\n Classes = str_c(str_extract(Title, \"\\\\d\"), \" Classes\"), .before = 1,\n across(where(is.numeric) & !ends_with(\"Value\"), ~ round(.x, 2)),\n across(ends_with(\"Value\"), ~ round(.x, 3))\n ) %>%\n select(-Title) %>%\n filter(Classes != \"1 Classes\")\nlpa_class_tab\n\n\n## RELIABILITY\n## =============================================================================\n\n.pb_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(inq_1:inq_6) %>%\n psych::alpha()\nround(summary(.pb_alpha)$raw_alpha, 2)\n\n.tb_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(inq_7r:inq_15r, inq_9, inq_11:inq_12) %>%\n psych::alpha()\nsummary(.tb_alpha)\n\n.bhs_20 <-\n lpa_fin %>%\n filter(!is.na(c), session == 1, protocol == 1) %>%\n ## names() %>%\n select(starts_with(\"bhs_\") & ends_with(\"r\")) %>%\n psych::alpha()\nsummary(.bhs_20)\n\n.bhs_sf <-\n lpa_fin %>%\n filter(!is.na(c), session == 1, protocol == 2) %>%\n select(starts_with(\"bhsSF_\") & ends_with(\"lessone\")) %>%\n psych::alpha()\nsummary(.bhs_sf)\n\n.si_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(mssi_1:mssi_4, starts_with(\"rc_mssi\")) %>%\n psych::alpha()\nsummary(.si_alpha)\n\n.ac_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(c(acss_8r:acss_13r, acss_7, acss_11, acss_14, acss_19)) %>%\n psych::alpha()\nsummary(.ac_alpha)\n\n.dep_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"dass_\", c(3, 5, 10, 13, 16, 17, 21))) %>%\n psych::alpha()\nsummary(.dep_alpha)\n\n.anx_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"dass_\", c(2, 4, 7, 9, 15, 19, 20))) %>%\n psych::alpha()\nsummary(.anx_alpha)\n\n.str_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"dass_\", c(1, 6, 8, 11, 12, 14, 18))) %>%\n psych::alpha()\nsummary(.str_alpha)\n\n.is_aff_reg_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(1, 14, 27))) %>%\n psych::alpha()\nsummary(.is_aff_reg_alpha)\n\n.is_int_bnd_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(2, 15, 28))) %>%\n psych::alpha()\nsummary(.is_int_bnd_alpha)\n\n.is_slf_pun_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(3, 16, 29))) %>%\n psych::alpha()\nsummary(.is_slf_pun_alpha)\n\n.is_slf_car_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(4, 17, 30))) %>%\n psych::alpha()\nsummary(.is_slf_car_alpha)\n\n.is_ant_dis_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(5, 18, 31))) %>%\n psych::alpha()\nsummary(.is_ant_dis_alpha)\n\n.is_ant_sui_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(6, 19, 32))) %>%\n psych::alpha()\nsummary(.is_ant_sui_alpha)\n\n.is_sen_see_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(7, 20, 33))) %>%\n psych::alpha()\nsummary(.is_sen_see_alpha)\n\n.is_peer_bnd_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(8, 21, 34))) %>%\n psych::alpha()\nsummary(.is_peer_bnd_alpha)\n\n.is_int_inf_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(9, 22, 35))) %>%\n psych::alpha()\nsummary(.is_int_inf_alpha)\n\n.is_tough_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(10, 23, 36))) %>%\n psych::alpha()\nsummary(.is_tough_alpha)\n\n.is_mrk_dis_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(11, 24, 37))) %>%\n psych::alpha()\nsummary(.is_mrk_dis_alpha)\n\n.is_revenge_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(12, 25, 38))) %>%\n psych::alpha()\nsummary(.is_revenge_alpha)\n\n.is_autnmy_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(13, 36, 39))) %>%\n psych::alpha()\nsummary(.is_autnmy_alpha)\n\n.isas_intra_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(\n 1, 3, 5, 6, 11, 14, 16, 18,\n 19, 24, 27, 29, 31, 32, 37\n ))) %>%\n psych::alpha()\nsummary(.isas_intra_alpha)\n\n.isas_social_alpha <-\n lpa_fin %>%\n filter(!is.na(c), session == 1) %>%\n select(num_range(\"isas_8.\", c(\n 2, 4, 7, 8, 9, 10, 12, 13, 15, 17,\n 20, 21, 22, 23, 25, 26, 28, 30,\n 33, 34, 35, 36, 38, 39\n ))) %>%\n psych::alpha()\nsummary(.isas_social_alpha)\n\nmeas_rel <-\n list(\n pb_rel = .pb_alpha,\n tb_rel = .tb_alpha,\n bhs_20_rel = .bhs_20,\n bhs_sf_rel = .bhs_sf,\n si_rel = .si_alpha,\n ac_rel = .ac_alpha,\n dep_rel = .dep_alpha,\n anx_rel = .anx_alpha,\n str_rel = .str_alpha,\n is_aff_reg_rel = .is_aff_reg_alpha,\n is_int_bnd_rel = .is_int_bnd_alpha,\n is_slf_pun_rel = .is_slf_pun_alpha,\n is_slf_car_rel = .is_slf_car_alpha,\n is_ant_dis_rel = .is_ant_dis_alpha,\n is_ant_sui_rel = .is_ant_sui_alpha,\n is_sen_see_rel = .is_sen_see_alpha,\n is_peer_bnd_rel = .is_peer_bnd_alpha,\n is_int_inf_rel = .is_int_inf_alpha,\n is_tough_rel = .is_tough_alpha,\n is_mrk_dis_rel = .is_mrk_dis_alpha,\n is_revenge_rel = .is_revenge_alpha,\n is_autnmy_rel = .is_autnmy_alpha,\n isas_intra_rel = .isas_intra_alpha,\n isas_social_rel = .isas_social_alpha\n )\n\n## NOTE e.g., do not run\n## round(summary(meas_rel$pb_rel)$raw_alpha, 2)\n\n## CROSS SECTIONAL: NON-MIDS\n## =============================================================================\n## Demographics\n## ANOVA on c predicted by covariates\n## Pairwise t-tests of covariates and LPA predcitors\n## NOTE optional analysis\n## multilevel regression si ~ all the other vars\n## Multinomial logistic regression instead of AVNOA, this could be interesting\n## to note units of change and therefor more informative from a prediction\n## perspective.\n\n.lpa_ch <-\n complete(imp_t0t1_50, \"long\") %>%\n group_by(.imp) %>%\n mutate(\n ## calculate change scores\n pb_ch = pb_1 - pb_0,\n tb_ch = tb_1 - tb_0,\n bhs_ch = bhs_1 - bhs_0,\n pb_ch = pb_1 - pb_0,\n dep_ch = dep_1 - dep_0,\n anx_ch = anx_1 - anx_0,\n ## str_ch = str_1 - str_0,\n si_ch = si_1 - si_0,\n ) %>%\n ungroup()\n\nAnova(lm(pb_ch + tb_ch + bhs_ch + dep_ch + anx_ch + si_ch ~ c, .lpa_ch), type = 3) %>%\n tidy()\n\n.lpa_ch %>%\n select(ends_with(\"_ch\")) %>%\n ## select(ends_with(\"_1\")) %>%\n map(~ tidy(pairwise.t.test(., .lpa_ch$c, p.adjust.method = \"bonf\")))\n\n.fu_look <-\n .lpa_ch %>%\n mutate(across(everything(), as.numeric)) %>%\n group_by(c) %>%\n summarise(\n across(ends_with(c(\"_0\", \"_1\", \"_ch\")) & !contains(\"session\"),\n list(m = mean, sd = sd),\n .names = \"{.col}_{.fn}\")\n ) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"_0|_1|_ch|_m|_sd\", \"\"),\n time = str_extract(key, \"0|1|ch\"),\n type = str_extract(key, \"m|sd\"),\n ) %>%\n select(-key) %>%\n pivot_wider(\n names_from = c(c, type),\n values_from = value\n ) %>%\n clean_names()\n\n\n.lpa_ch %>%\n mutate(across(everything(), as.numeric)) %>%\n group_by(c) %>%\n t.test(pb_0 ~ pb_1, paried = TRUE)\n\nAnova(lm(pb_ch ~ c, data = .lpa_ch), type = 3) %>%\n\n## ANOVA BY C\n## -----------------------------------------------------------------------------\n## TODO get the standardised residuals of each model\n## https://www.statology.org/standardized-residuals-in-r/\n## get anova stats\n.lpa_anova <-\n Anova(lm(c ~ pb + tb + bhs + ac, data = lpa_t0), type = 3) %>%\n tidy() %>%\n select(term, statistic, p.value) %>%\n filter(term %in% c(\"pb\", \"tb\", \"bhs\", \"ac\")) %>%\n rename(var = term)\n\n.lpa_anova_fx <-\n Anova(lm(c ~ pb + tb + bhs + ac, data = lpa_t0), type = 3) %>%\n effectsize::omega_squared() %>%\n as.data.frame() %>%\n mutate(\n fx_size = interpret_omega_squared(Omega2_partial, rules = \"cohen1992\"),\n omega2 = str_remove(sprintf(\"%.2f\", Omega2_partial), \"^0\")\n ) %>%\n rename(var = Parameter) %>%\n select(var, omega2, fx_size)\n\nlpa_sample <-\n lpa_t0 %>%\n select(c, pb, tb, bhs, ac) %>%\n mutate(\n across(c(pb, tb, bhs, ac),\n ~ scale(.x),\n .names = \"{.col}_z\"\n )\n ) %>%\n group_by(c) %>%\n summarise(\n N = n(),\n across(\n c(pb, tb, bhs, ac),\n list(\n m = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE)\n )\n ),\n across(\n ends_with(\"_z\"),\n ~ round(mean(.x, na.rm = TRUE), 2)\n )\n ) %>%\n mutate(\n prc = N / sum(N) * 100, .after = N\n ) %>%\n select(c, N, prc, ends_with(c(\"_m\", \"_sd\", \"_z\")))\n\n## get mean and sd\n.lpa_m_sd <-\n ## NOTE cannot find original object `lpa_sample`\n lpa_sample %>%\n select(c, ends_with(c(\"_m\", \"_sd\"))) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"_m|_sd\", \"\"),\n type = str_extract(key, \"m|sd\"),\n ) %>%\n select(-key) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value\n )\n\n## make contrasts wide\n.lpa_t <-\n lpa_t0 %>%\n select(pb, tb, bhs, ac) %>%\n map( ~ tidy(pairwise.t.test(., lpa_t0$c, p.adjust.method = \"bonf\"))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n mutate(across(where(is.character), as.numeric)) %>%\n pivot_longer(\n ## col = ends_with(c(\"group1\", \"group2\", \"value\")),\n col = contains(\"value\"),\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n select(contains(c(\"group1\", \"group2\")), everything()) %>%\n arrange(key) %>%\n ## NOTE this was taken from SO\n ## TODO find tidy way to do this\n .[!duplicated(lapply(., summary))] %>%\n rename(group1 = 1, group2 = 2, var = 3, p_value = 4) %>%\n mutate(\n var = str_extract(var, \"ac|bhs|pb|tb\")\n ) %>%\n pivot_wider(\n names_from = c(group1, group2),\n values_from = p_value\n ) %>%\n janitor::clean_names()\n\n## join the lpa var tests\n.lpa_anova_full <-\n .lpa_anova %>%\n full_join(.lpa_m_sd) %>%\n full_join(.lpa_t) %>%\n full_join(.lpa_anova_fx)\n\n## generate the table for the report\nlpa_anova <-\n .lpa_anova_full %>%\n mutate(\n group_diff = case_when(\n x2_1 < 0.05 & x3_1 < 0.05 & x3_2 < 0.05 ~ \"1 < 2,3; 2 < 3\",\n TRUE ~ \"null\"\n ),\n f = case_when(\n ## TODO check that markdown formatted asterix\n p.value < 0.001 ~ paste0(as.character(round(statistic, 2)), \"<sup>***</sup>\"),\n p.value < 0.05 ~ paste0(as.character(round(statistic, 2)), \"<sup>*</sup>\")\n ),\n across(starts_with(c(\"m_\", \"sd_\")), ~ sprintf(\"%.2f\", round(.x, 2))),\n low_si = str_c(m_1, \" (\", sd_1, \")\"),\n mod_si = str_c(m_2, \" (\", sd_2, \")\"),\n hi_si = str_c(m_3, \" (\", sd_3, \")\"),\n var = case_when(\n var == \"pb\" ~ \"Burdensomeness\",\n var == \"tb\" ~ \"Belongingness\",\n var == \"ac\" ~ \"Capability\",\n var == \"bhs\" ~ \"Hopelessness\"\n )\n ) %>%\n select(var, low_si:hi_si, f, omega2, group_diff)\n\n## FOLLOW UP ANALYSIS BY C\n## -----------------------------------------------------------------------------\n## recode values of demo vars\n.fu <-\n ## use full data set here to look at isas items\n lpa_fin %>%\n filter(!is.na(c), session %in% 1 | final_session %in% 1) %>%\n group_by(ID) %>%\n mutate(\n total_session = max(session), .after = session,\n planned_final = max(planned_final)\n ) %>%\n ungroup() %>%\n filter(session %in% 1) %>%\n select(\n ID:suburb, isas_1a:isas_1m, isas_3a,\n isas_3b, isas_4:isas_7, siss\n ) %>%\n mutate(\n siss_n = as.numeric(as.character(siss)), .after = everything(),\n isas_days = floor(\n interval(isas_3b, date) / duration(n = 1, unit = \"days\")),\n ## .before = date,\n ## remove strings from isas vars and make all numeric class\n across(c(isas_1i, isas_1m), ~ str_extract(.x, \"\\\\d*\")),\n across(starts_with(\"isas_1\"), as.numeric),\n ## prorate isas higher-order factors for comparison\n isas_intra = isas_intra / 5,\n isas_social = isas_social / 8\n ) %>%\n relocate(isas_days, .before = date) %>%\n rowwise() %>%\n mutate(\n ## count the number of methods of dsh\n dsh_n = sum(c_across(starts_with(\"isas_1\")) > 0, na.rm = TRUE)\n ) %>%\n ungroup() %>%\n mutate(\n across(isas_4:isas_7, as.factor),\n sex = case_when(\n sex == 0 ~ \"Female\",\n sex == 1 ~ \"Male\",\n sex == 2 ~ \"Other\",\n TRUE ~ NA_character_\n ),\n atsi = case_when(\n ## correct atsi status from old sp_db\n sp_db == \"old\" & atsi == 4 ~ \"Neither\",\n ## collapse categories as cell sizes too small\n atsi == 0 ~ \"Neither\",\n atsi %in% c(1, 2, 3) ~ \"ATSI\",\n TRUE ~ NA_character_\n ),\n living = case_when(\n living == 1 ~ \"family\",\n living == 2 ~ \"other\",\n living == 3 ~ \"other\",\n living == 4 ~ \"other\",\n TRUE ~ NA_character_\n ),\n relation = case_when(\n relation == 1 ~ \"single/separated/divorced\",\n relation == 2 ~ \"relationship/married/defacto\",\n relation == 3 ~ \"relationship/married/defacto\",\n relation == 4 ~ \"single/separated/divorced\",\n TRUE ~ NA_character_\n ),\n income = case_when(\n income == 1 ~ \"working\",\n income == 2 ~ \"working\",\n income == 3 ~ \"not-working\",\n income == 4 ~ \"not-working\",\n income == 5 ~ \"not-working\",\n income == 6 ~ \"student\",\n TRUE ~ NA_character_\n ),\n across(starts_with(\"diag\"), ~ case_when(\n .x == 1 ~ \"no diagnosis warranted\",\n .x == 2 ~ \"depression\",\n .x == 3 ~ \"generalised anxiety\",\n .x == 4 ~ \"phobia\",\n .x == 5 ~ \"adjustment disorder\",\n .x == 6 ~ \"acute stress disorder\",\n .x == 7 ~ \"ptsd\",\n .x == 8 ~ \"obsessive-compulsive\",\n .x == 9 ~ \"eating disorder\",\n .x == 13 ~ \"bipolar\",\n .x == 15 ~ \"personality disorder (e.g. BPD)\",\n .x == 17 ~ \"neurological disorder\",\n .x == 18 ~ \"other\",\n TRUE ~ NA_character_\n )),\n across(\n c(\n c, planned_final, sp_db, location, sex, atsi:diag1,\n starts_with(c(\"past\", \"hx\", \"fxhx\")), suicide_attempts,\n isas_4:isas_7, siss\n ) & -siss_n,\n as.factor),\n )\n.fu\n\n## CONTINUOUS VARS\n## -------------------------------\n## class sizes full only\n.lpa_n_fu_full <-\n .fu %>%\n select(ID) %>%\n summarise(\n n_sd = sum(!is.na(ID))\n ) %>%\n mutate(\n n_m = round(n_sd / sum(n_sd) * 100, 2),\n c = \"full\"\n )\n\n## class sizes all\n.lpa_n_fu <-\n .fu %>%\n select(c, ID) %>%\n group_by(c) %>%\n summarise(\n n_sd = sum(!is.na(ID)),\n ) %>%\n mutate(\n n_m = round(n_sd / sum(n_sd) * 100, 2),\n c = as.character(c)\n ) %>%\n bind_rows(.lpa_n_fu_full) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"_m|_sd\", \"\"),\n type = str_extract(key, \"(m|sd)$\")\n ) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value,\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup() %>%\n mutate(var = str_replace(var, \"n\", \"Class size\")) %>%\n select(var, starts_with(c(\"m\", \"sd\"))) %>%\n relocate(ends_with(\"_full\"), .after = var)\n.lpa_n_fu\n\n## univariate anovas of covariates\n## PARAMETRIC\n.fu_test <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n map(~ tidy(Anova(lm(. ~ .fu$c, na.action = na.exclude), type = 3))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n filter(across(ends_with(\"term\"), ~ .x == \".fu$c\")) %>%\n select(ends_with(c(\"statistic\", \"p.value\"))) %>%\n pivot_longer(\n everything(),\n names_to = c(\"var\", \"stat\"),\n names_sep = \"\\\\.\"\n ) %>%\n pivot_wider(\n names_from = stat,\n values_from = value\n )\n.fu_test %>%\n print(n = Inf)\n\n## individual contrasts of fu covariate data\n.fu_test_cntrst <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n map(~ tidy(pairwise.t.test(., .fu$c, p.adjust.method = \"bonf\"))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n mutate(across(where(is.character), as.numeric)) %>%\n pivot_longer(\n col = contains(\"value\"),\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n select(contains(c(\"group1\", \"group2\")), everything()) %>%\n arrange(key) %>%\n ## TODO find tidy way to do this\n .[!duplicated(lapply(., summary))] %>%\n rename(group1 = 1, group2 = 2, var = 3, p_value = 4) %>%\n mutate(\n var = str_replace(var, \"\\\\.p\\\\.value\", \"\")\n ) %>%\n pivot_wider(\n names_from = c(group1, group2),\n values_from = p_value\n ) %>%\n janitor::clean_names()\n\n.fu_test_cntrst %>%\n print(n = Inf)\n\n## univariate anovas of covariates\n## NON-PARAMETRIC: KRUSKAL WALLIS\n## NOTE used to check that that the results for the isas section 1 is the same\n## NOT used for anlalysis otherwise\n.fu_test_kw <-\n .fu %>%\n ## .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n map(~ kruskal.test(. ~ .fu$c, na.action = na.exclude), type = 3) %>%\n map(`[`, c(\"statistic\", \"p.value\")) %>%\n unlist(recursive = FALSE) %>%\n map(~ as_tibble(.)) %>%\n map(~ janitor::clean_names(.)) %>%\n bind_rows(., .id = \"var\") %>%\n rename(key = var) %>%\n mutate(\n var = str_replace_all(key, \"\\\\.statistic|\\\\.p\\\\.value\", \"\"),\n type = str_extract(key, \"statistic|p\\\\.value\"),\n ) %>%\n pivot_wider(\n names_from = type,\n values_from = value\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup() %>%\n rename(\n statistic_kw = statistic,\n p_kw = p.value\n )\n\n.fu_test_kw %>%\n print(n = Inf)\n\n## individual contrasts of fu covariate data\n.fu_test_kw_cntrst <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n ## map(~ pgirmess::kruskalmc(. ~ .fu$c)$dif.com) %>%\n map(~ agricolae::kruskal(., .fu$c, p.adj = \"bonf\", group = FALSE)$comparison) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n comp = c(\"1v2\", \"1v3\", \"2v3\"),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n as_tibble() %>%\n select(var, comp, everything()) %>%\n arrange(var, comp)\n\n.fu_test_kw_cntrst %>%\n print(n = Inf)\n\n## join the the anova and kw tests to eyeball for consitency\n.fu_test_cont_full <-\n .lpa_m_sd_fu %>%\n full_join(.fu_test_kw ) %>%\n full_join(.fu_test_kw_cntrst) %>%\n full_join(.fu_test) %>%\n full_join(.fu_test_cntrst)\n## NOTE eyebal check finishes here\n\n\n## CONTINOUS MEANS\n.lpa_m_sd_fu_full <-\n .fu %>%\n select(\n si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n mutate(across(everything(), as.numeric)) %>%\n summarise(\n across(\n everything(),\n list(\n m = ~ round(mean(.x, na.rm = TRUE), 2),\n sd = ~ round(sd(.x, na.rm = TRUE), 2)\n )\n )\n ) %>%\n mutate(\n c = \"full\", .before = 1\n ) %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"(?!_mrk_)(_m|_sd)\", \"\"),\n type = str_extract(key, \"(m|sd)$\"),\n ) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value,\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup()\n\n.lpa_m_sd_fu_full %>%\n print(n = Inf)\n\n## means and sd by class\n.lpa_m_sd_fu <-\n .fu %>%\n select(\n c, si, dep, anx, str, age, starts_with(\"is_\"),\n isas_intra, isas_social, total_session,\n isas_days, isas_1a:isas_1m, dsh_n, isas_3a, siss_n\n ) %>%\n mutate(across(everything(), as.numeric)) %>%\n group_by(c) %>%\n summarise(\n across(\n everything(),\n list(\n m = ~ round(mean(.x, na.rm = TRUE), 2),\n sd = ~ round(sd(.x, na.rm = TRUE), 2)\n )\n )\n ) %>%\n ungroup() %>%\n pivot_longer(\n cols = !c,\n names_to = \"key\",\n values_to = \"value\"\n ) %>%\n mutate(\n var = str_replace_all(key, \"(?!_mrk_)(_m|_sd)\", \"\"),\n type = str_extract(key, \"(m|sd)$\"),\n ) %>%\n pivot_wider(\n names_from = c(type, c),\n values_from = value,\n ) %>%\n select(-key) %>%\n group_by(var) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup() %>%\n full_join(.lpa_m_sd_fu_full) %>%\n select(var, ends_with(\"full\"), everything())\n\n.lpa_m_sd_fu %>%\n print(n = Inf)\n\n## FULL CONTINUOUS FU TABLE\nfu_test_cont <-\n .lpa_n_fu %>%\n full_join(.lpa_m_sd_fu) %>%\n full_join(.fu_test) %>%\n full_join(.fu_test_cntrst)\n\nfu_test_cont %>%\n print(n = Inf)\n\n## CATEGORICAL FU VARS\n## --------------------------\n## chisq on fu covariate data\n.fu_test_nonp <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex)\n ) %>%\n select(\n sex,\n location, atsi:diag1,\n starts_with(c(\"past\", \"hx\", \"fxhx\")),\n suicide_attempts, planned_final,\n isas_4:isas_7, siss\n ) %>%\n mutate(location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n map(~ tidy(chisq.test(., .fu$c))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n filter(across(ends_with(\"term\"), ~ .x == \".fu$c\")) %>%\n select(ends_with(c(\"statistic\", \"p.value\"))) %>%\n pivot_longer(\n everything(),\n names_to = c(\"var\", \"stat\"),\n names_sep = \"\\\\.\"\n ) %>%\n pivot_wider(\n names_from = stat,\n values_from = value\n )\n\n.fu_test_nonp %>%\n print(n = Inf)\n\n\n.fu_resd_use <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex),\n hxsa = case_when(\n suicide_attempts %in% \"0\" ~ \"no\",\n suicide_attempts %!in% c(NA, \"0\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n dx1 = case_when(\n diag1 %in% \"no diagnosis warranted\" ~ \"no\",\n diag1 %!in% c(NA, \"no diagnosis warranted\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n across(c(hxsa, dx1), as.factor),\n ) %>%\n select(\n location,\n sex,\n atsi:diag1, dx1,\n starts_with(c(\"past\", \"hx\", \"fxhx\")),\n suicide_attempts, planned_final,\n isas_4:isas_7, siss\n ) %>%\n mutate(\n location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n ## map(~ chisq.test(table(., .fu_bin$c))) %>%\n map(~ chisq.test(table(., .fu$c))) %>%\n map(`[`, c(\"stdres\", \"statistic\", \"p.value\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n rename(\n ## to match cont df\n value = stdres,\n p = p_value,\n x2_1 = x1,\n x3_1 = x2,\n x3_2 = x3\n )\n.fu_resd_use\n\n\n.lpa_m_sd_fu_nonp_full <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex),\n hxsa = case_when(\n suicide_attempts %in% \"0\" ~ \"no\",\n suicide_attempts %!in% c(NA, \"0\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n dx1 = case_when(\n diag1 %in% \"no diagnosis warranted\" ~ \"no\",\n diag1 %!in% c(NA, \"no diagnosis warranted\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n across(c(hxsa, dx1), as.factor),\n ) %>%\n select(where(is.factor)) %>%\n mutate(location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n pivot_longer(cols = -c) %>%\n count(name, value) %>%\n filter(!is.na(value)) %>%\n mutate(c = \"full\", .before = 1) %>%\n group_by(c, name) %>%\n rename(sd = n) %>%\n mutate(across(\n sd,\n ~ round(.x / sum(.x) * 100, 2),\n .names = \"m\"\n )) %>%\n mutate(\n row = row_number(),\n across(where(is.integer), as.double)\n ) %>%\n pivot_wider(\n names_from = c,\n values_from = c(sd, m),\n ) %>%\n select(-row)\n\n.lpa_m_sd_fu_nonp_full %>%\n print(n = Inf)\n\n.lpa_m_sd_fu_nonp <-\n .fu %>%\n mutate(\n sex = as.character(sex),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n sex = as.factor(sex),\n hxsa = case_when(\n suicide_attempts %in% \"0\" ~ \"no\",\n suicide_attempts %!in% c(NA, \"0\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n dx1 = case_when(\n diag1 %in% \"no diagnosis warranted\" ~ \"no\",\n diag1 %!in% c(NA, \"no diagnosis warranted\") ~ \"yes\",\n TRUE ~ NA_character_\n ),\n across(c(hxsa, dx1), as.factor),\n ) %>%\n select(where(is.factor)) %>%\n mutate(location = case_when(\n location == \"Wollongong\" ~ \"metro\",\n location %in% c(\"Nowra\", \"Bega\") ~ \"regional\",\n TRUE ~ NA_character_\n )) %>%\n pivot_longer(cols = -c) %>%\n count(c, name, value) %>%\n filter(!is.na(value)) %>%\n group_by(c, name) %>%\n rename(sd = n) %>%\n pivot_wider(\n names_from = c,\n names_prefix = \"sd_\",\n values_from = sd,\n ) %>%\n ## full_join(.lpa_m_sd_fu_nonp_sex) %>%\n full_join(.lpa_m_sd_fu_nonp_full) %>%\n mutate(\n across(\n sd_1:sd_3,\n ~ round(.x / sd_full * 100, 2),\n .names = \"m_{1:3}\"\n )) %>%\n arrange(name, value) %>%\n group_by(name, value) %>%\n ## summarise(across(everything(), .coalesce_by_column)) %>%\n arrange(name, desc(m_full)) %>%\n ungroup() %>%\n select(name, value, ends_with(c(\"full\", \"1\", \"2\", \"3\")), everything()) %>%\n rename(var = name) %>%\n full_join(.fu_resd_use) %>%\n rename(name = var) %>%\n mutate(\n across(where(is.integer), as.double),\n value = case_when(\n str_detect(name, \"^(cald|past|hx|fxhx|planned)\") & value == 0 ~ \"no\",\n str_detect(name, \"^(cald|past|hx|fxhx|planned)\") & value == 1 ~ \"yes\",\n name == \"suicide_attempts\" & value == 0 ~ \"none\",\n name == \"suicide_attempts\" & value == 1 ~ \"one\",\n name == \"suicide_attempts\" & value == 2 ~ \"two\",\n name == \"suicide_attempts\" & value == 3 ~ \"three\",\n name == \"suicide_attempts\" & value == 4 ~ \"four or more\",\n name == \"living\" & value == 1 ~ \"family\",\n name == \"living\" & value == 2 ~ \"friends\",\n name == \"living\" & value == 3 ~ \"other(s)\",\n name == \"living\" & value == 4 ~ \"alone\",\n str_detect(name, \"^(isas_4|isas_5|isas_7)\") & value == 1 ~ \"yes\",\n str_detect(name, \"^(isas_4|isas_5|isas_7)\") & value == 2 ~ \"sometimes\",\n str_detect(name, \"^(isas_4|isas_5|isas_7)\") & value == 3 ~ \"no\",\n name == \"isas_6\" & value == 1 ~ \"00--1 hr\",\n name == \"isas_6\" & value == 2 ~ \"01--3 hrs\",\n name == \"isas_6\" & value == 3 ~ \"03--6 hrs\",\n name == \"isas_6\" & value == 4 ~ \"06--12 hrs\",\n name == \"isas_6\" & value == 5 ~ \"12--24 hrs\",\n name == \"isas_6\" & value == 6 ~ \"24+ hrs\",\n name == \"siss\" & value == 0 ~ \"si_none\",\n name == \"siss\" & value == 1 ~ \"si_vague\",\n name == \"siss\" & value == 2 ~ \"si_clear_no_plan\",\n name == \"siss\" & value == 3 ~ \"si_plans_no_int\",\n name == \"siss\" & value == 4 ~ \"si_plans_int\",\n name == \"siss\" & value == 5 ~ \"si_attempt\",\n TRUE ~ as.character(value)\n )\n ) %>%\n rename(var = name)\n\nfu_test_cont_nonp <-\n .lpa_m_sd_fu_nonp %>%\n mutate(test = \"fact\")\n\nfu_test_cont_nonp %>%\n print(n = Inf)\n\nfu_test_cont <-\n .lpa_n_fu %>%\n full_join(.lpa_m_sd_fu) %>%\n full_join(.fu_test) %>%\n full_join(.fu_test_cntrst) %>%\n mutate(\n test = \"cont\",\n bonf = 0.05\n )\n\n.fu_test_full <-\n fu_test_cont %>%\n rename(value = var) %>%\n mutate(var = \"\") %>%\n full_join(fu_test_cont_nonp) %>%\n relocate(c(var, value), .before = 1) %>%\n filter(var != \"sp_db\") %>%\n mutate(\n value = case_when(\n value == \"isas_1a\" ~ \"dsh_cutting\",\n value == \"isas_1b\" ~ \"dsh_biting\",\n value == \"isas_1c\" ~ \"dsh_burning\",\n value == \"isas_1d\" ~ \"dsh_carving\",\n value == \"isas_1e\" ~ \"dsh_pinching\",\n value == \"isas_1f\" ~ \"dsh_pull_hair\",\n value == \"isas_1g\" ~ \"dsh_scratch\",\n value == \"isas_1h\" ~ \"dsh_hit_self\",\n value == \"isas_1i\" ~ \"dsh_int_wound\",\n value == \"isas_1j\" ~ \"dsh_skin_rub\",\n value == \"isas_1k\" ~ \"dsh_needles\",\n value == \"isas_1l\" ~ \"dsh_swall_dang\",\n value == \"isas_1m\" ~ \"dsh_other\",\n value == \"isas_3a\" ~ \"dsh_age_onset\",\n TRUE ~ as.character(value)\n ),\n var = case_when(\n str_detect(value, \"^(is_|isas_|dsh_)\") ~ \"dsh\",\n var == \"isas_4\" ~ \"dsh_pain\",\n var == \"isas_5\" ~ \"dsh_alone\",\n var == \"isas_6\" ~ \"dsh_time_wait\",\n var == \"isas_7\" ~ \"dsh_want_stop\",\n value %in% c(\"dep\", \"anx\", \"str\") ~ \"mood\",\n value == \"total_session\" ~ \"total_session\",\n value == \"si\" ~ \"si\",\n value == \"siss_n\" ~ \"siss_n\",\n value == \"age\" ~ \"age\",\n value == \"Class size\" ~ \"Class size\",\n TRUE ~ as.character(var)\n ),\n type = case_when(\n var == \"Class size\" ~ \"Class size\",\n str_detect(var, \"^(age|atsi|cald|inc|liv|loc|rel|sex)\") ~ \"demo\",\n str_detect(var, \"^(mood|dsh|si|siss_n)\") ~ \"psych\",\n str_detect(var, \"^(fxhx|hx|suicide_|hxsa)\") ~ \"risk\",\n str_detect(var, \"^(total|diag|dx1|past|planned)\") ~ \"tx\"\n ),\n order = case_when(\n str_detect(var, \"^(dx1|diag|pastmed|planned)\") ~ 0,\n var == \"location\" ~ 0,\n var == \"Class size\" ~ 1,\n var == \"sex\" ~ 2,\n var == \"age\" ~ 3,\n var == \"atsi\" ~ 4,\n var == \"cald\" ~ 5,\n var == \"income\" ~ 6,\n var == \"living\" ~ 7,\n var == \"relation\" ~ 8,\n var == \"mood\" ~ 9,\n value == \"dsh_n\" ~ 10,\n value == \"isas_intra\" ~ 11,\n str_detect(value, \"^(is_)(aff|slf_p|ant|mrk)\") ~ 12,\n value == \"isas_social\" ~ 13,\n str_detect(value, \"^(is_)(slf_c|tou|aut|sen|int|rev|peer)\") ~ 14,\n value == \"isas_days\" ~ 15,\n str_detect(var, \"^(dsh_)(alone|want|time|pain)\") ~ 16,\n value == \"dsh_age_onset\" ~ 17,\n str_detect(value,\n \"^(dsh_)(cu|bi|bu|ca|pi|pu|sc|hi|in|sk|ne|sw)\") ~ 18,\n str_detect(value,\n \"dsh_other\") ~ 19,\n str_detect(var, \"hxdsh|fxhxdsh\") ~ 20,\n str_detect(var, \"^(hx|fxhx)(trauma|suic|sa)\") ~ 21,\n var == \"si\" ~ 24,\n value == \"si_clear_no_plan\" ~ 23,\n value == \"si_vague\" ~ 25,\n value == \"si_plans_no_int\" ~ 22,\n value == \"si_none\" ~ 26,\n value == \"si_plans_int\" ~ 27,\n value == \"si_attempt\" ~ 28,\n var == \"siss_n\" ~ 29,\n var == \"pastpsyc\" ~ 30,\n var == \"total_session\" ~ 31\n )\n ) %>%\n relocate(c(order, type), .before = 1) %>%\n arrange(order, type, var, desc(m_full))\n\n## generate the table for the report\n## load the manual group diff entries\n.fu_test_grp_dif <-\n readxl::read_xlsx(\n \"./data/tmp/anova_results.xlsx\",\n sheet = \"group_diff\"\n ) %>%\n select(order:value, group_diff) %>%\n mutate(order = as.double(order))\n\nfu_test <-\n .fu_test_full %>%\n select(-order) %>%\n mutate(\n across(\n starts_with(\"sd_\"),\n ~ case_when(\n .data$test == \"cont\" && .data$var != \"Class size\" ~\n sprintf(\"%.2f\", round(.x, 2)),\n TRUE ~ as.character(.x)\n )\n ),\n across(\n starts_with(\"m_\"), ~ sprintf(\"%.2f\", round(.x, 2))\n ),\n full = str_c(m_full, \" (\", sd_full, \")\"),\n low_si = str_c(m_1, \" (\", sd_1, \")\"),\n mod_si = str_c(m_2, \" (\", sd_2, \")\"),\n hi_si = str_c(m_3, \" (\", sd_3, \")\"),\n chisq_f = case_when(\n ## html tags for supescript\n p < 0.001 ~ paste0(sprintf(\"%.2f\", statistic), \"<sup>***</sup>\"),\n p < 0.01 ~ paste0(sprintf(\"%.2f\", statistic), \"<sup>**</sup>\"),\n p < 0.05 ~ paste0(sprintf(\"%.2f\", statistic), \"<sup>*</sup>\"),\n p > 0.05 ~ paste0(sprintf(\"%.2f\", statistic), \"\"),\n TRUE ~ NA_character_\n ),\n ) %>%\n full_join(.fu_test_grp_dif) %>%\n ## relocate(starts_with(\"group_diff\"), .after = everything()) %>%\n group_by(order, type, var, value) %>%\n ## arrange(order, type, var, desc(m_3), .by_group = TRUE) %>%\n arrange(order, type, var, .by_group = TRUE) %>%\n ungroup() %>%\n relocate(value2, .before = value) %>%\n relocate(order, .before = 1)\n\n\n## MULTINOMIAL LOGISTIC REGRESSION\n## =============================================================================\n\n.fu_<-\n .fu %>%\n filter(!is.na(c)) %>%\n mutate(\n across(where(is.factor), as.character),\n suicide_attempts = as.numeric(suicide_attempts),\n sex = case_when(\n sex == \"Other\" ~ NA_character_,\n TRUE ~ as.character(sex)\n ),\n across(starts_with(c(\"cald\", \"past\", \"hx\", \"fxhx\", \"planned\")),\n ~ case_when(\n .x == 0 ~ \"no\",\n .x == 1 ~ \"yes\",\n TRUE ~ NA_character_\n )\n ),\n hxsa = case_when(\n suicide_attempts < 1 ~ \"no\",\n suicide_attempts >= 1 ~ \"yes\",\n TRUE ~ NA_character_\n ),\n ## living = case_when(\n ## living == 1 ~ \"family\",\n ## living == 2 ~ \"other\",\n ## living == 3 ~ \"other\",\n ## living == 4 ~ \"other\",\n ## TRUE ~ NA_character_\n ## ),\n across(starts_with(\"isas_\") & ends_with(c(\"4\", \"5\", \"7\")),\n ~ case_when(\n .x == 1 ~ \"yes\",\n .x == 2 ~ \"sometimes\",\n .x == 3 ~ \"no\",\n TRUE ~ NA_character_\n )\n ),\n across(where(is.character), as.factor),\n ) %>%\n rowwise() %>%\n mutate(\n dsh_n = sum(c_across(starts_with(\"isas_1\")) > 0, na.rm = TRUE)\n ) %>%\n ungroup() %>%\n rename(\n dsh_pain = isas_4,\n dsh_alone = isas_5,\n dsh_stop = isas_7\n ) %>%\n ## str()\n mutate(\n across(c(living, sex), ~ relevel(.x, ref = 2)),\n across(starts_with(\"dsh\") & -dsh_n, ~ relevel(.x, ref = 2)),\n across(relation, ~ relevel(.x, ref = 2)),\n across(atsi, ~ relevel(.x, ref = 2)),\n )\n\n.fu_ %>%\n select(starts_with(\"fx\")) %>%\n summary()\n\n.fu_log <-\n dfidx::dfidx(.fu_, choice = \"c\", shape = \"wide\")\n\n.fu_log_test_3 <-\n mlogit(\n c ~ 1 |\n dep + anx + str + age + si + sex +\n atsi + cald + income + living + relation +\n dsh_pain + dsh_alone + dsh_stop + dsh_n +\n isas_intra + isas_social +\n hxdsh + hxsa + hxsuicide + hxtrauma +\n ## fxhxdsh + fxhxsuicide + pastmed +\n pastpsyc,\n data = .fu_log,\n reflevel = \"3\",\n na.action = na.exclude\n )\n\nfu_ml_sum_3 <-\n summary(.fu_log_test_3)\n\n.fu_ml_3 <-\n fu_ml_sum_3$CoefTable %>%\n as.data.frame() %>%\n mutate(var = attr(fu_ml_sum_3$coefficients, \"names\"), .before = 1) %>%\n tibble() %>%\n clean_names() %>%\n mutate(\n exp_b = round(exp(estimate), 2), .after = var\n ) %>%\n full_join(\n tibble(clean_names(cbind(\n data.frame(var = attr(exp(confint(.fu_log_test_3)), \"dimnames\")[[1]]),\n data.frame(exp(confint(.fu_log_test_3)))\n )))\n ) %>%\n select(var, exp_b, x2_5, x97_5, pr_z) %>%\n filter(!str_detect(var, \"(Intercept)\")) %>%\n mutate(\n name = str_replace_all(var, \":\\\\d\", \"\"),\n comp = str_extract(var, \":\\\\d\"),\n comp = str_replace(comp, \":\", \"3v\"),\n ) %>%\n select(-var) %>%\n ## NOTE consider transfering below to rmd for table pres\n pivot_wider(\n names_from = comp,\n values_from = c(exp_b, x2_5, x97_5, pr_z)\n ) %>%\n select(name, ends_with(c(\"3v1\", \"3v2\"))) %>%\n mutate(\n ## across(where(is.double), ~ round(.x, 2)),\n across(starts_with(\"x\"), ~ round(.x, 2)),\n across(starts_with(\"pr_\"), ~ round(.x, 3))\n )\n\n.fu_ml_3 %>%\n filter(\n if_any(\n starts_with(\"pr_z\"),\n ~.x < 0.05)\n )\n\n.fu_log_test_2 <-\n mlogit(\n c ~ 1 |\n dep + anx + str + age + si + sex +\n atsi + cald + income + living + relation +\n dsh_pain + dsh_alone + dsh_stop + dsh_n +\n isas_intra + isas_social +\n hxdsh + hxsa + hxsuicide + hxtrauma +\n ## fxhxdsh + fxhxsuicide + pastmed +\n pastpsyc,\n data = .fu_log,\n reflevel = \"2\",\n na.action = na.exclude\n )\n\nfu_ml_sum_2 <-\n summary(.fu_log_test_2)\n\n.fu_ml_2 <-\n fu_ml_sum_2$CoefTable %>%\n as.data.frame() %>%\n mutate(var = attr(fu_ml_sum_2$coefficients, \"names\"), .before = 1) %>%\n tibble() %>%\n clean_names() %>%\n mutate(\n exp_b = round(exp(estimate), 2), .after = var\n ) %>%\n full_join(\n tibble(clean_names(cbind(\n data.frame(var = attr(exp(confint(.fu_log_test_2)), \"dimnames\")[[1]]),\n data.frame(exp(confint(.fu_log_test_2)))\n )))\n ) %>%\n select(var, exp_b, x2_5, x97_5, pr_z) %>%\n filter(!str_detect(var, \"(Intercept)\")) %>%\n mutate(\n name = str_replace_all(var, \":\\\\d\", \"\"),\n comp = str_extract(var, \":\\\\d\"),\n comp = str_replace(comp, \":\", \"2v\"),\n ) %>%\n select(-var) %>%\n ## NOTE consider transfering below to rmd for table pres\n pivot_wider(\n names_from = comp,\n values_from = c(exp_b, x2_5, x97_5, pr_z)\n ) %>%\n select(name, ends_with(c(\"2v1\"))) %>%\n mutate(\n across(starts_with(\"x\"), ~ round(.x, 2)),\n across(starts_with(\"pr_\"), ~ round(.x, 3))\n )\n\n.fu_ml_2 %>%\n filter(\n if_any(\n starts_with(\"pr_z\"),\n ~ .x < 0.05)\n )\n\n\n## table for presentation\n.fu_ml_order <-\n readxl::read_xlsx(\n \"./data/tmp/anova_results.xlsx\",\n sheet = \"ml_table\"\n ) %>%\n select(order:name2) %>%\n mutate(order = as.double(order))\n\nfu_ml <-\n .fu_ml_3 %>%\n full_join(.fu_ml_2) %>%\n select(name, ends_with(c(\"2v1\", \"3v1\", \"3v2\"))) %>%\n full_join(.fu_ml_order) %>%\n arrange(order) %>%\n relocate(c(order, type, name2), .before = name)\n\nfu_ml %>%\n print(n = Inf)\n\n\n.fu_log_test_1 <-\n mlogit(\n c ~ 1 |\n dep + anx + str + age + si + sex +\n atsi + cald + income + living + relation +\n dsh_pain + dsh_alone + dsh_stop + dsh_n +\n isas_intra + isas_social +\n hxdsh + hxsa + hxsuicide + hxtrauma +\n ## fxhxdsh + fxhxsuicide + pastmed +\n pastpsyc,\n data = .fu_log,\n reflevel = \"1\",\n )\n\nfu_ml_sum_1 <-\n summary(.fu_log_test_1)\n\n.fu_ml_1 <-\n fu_ml_sum_1$CoefTable %>%\n as.data.frame() %>%\n mutate(var = attr(fu_ml_sum_1$coefficients, \"names\"), .before = 1) %>%\n tibble() %>%\n clean_names() %>%\n mutate(\n exp_b = round(exp(estimate), 2), .after = var\n ) %>%\n full_join(\n tibble(clean_names(cbind(\n data.frame(var = attr(exp(confint(.fu_log_test_1)), \"dimnames\")[[1]]),\n data.frame(exp(confint(.fu_log_test_1)))\n )))\n ) %>%\n select(var, exp_b, x2_5, x97_5, pr_z) %>%\n filter(!str_detect(var, \"(Intercept)\")) %>%\n mutate(\n name = str_replace_all(var, \":\\\\d\", \"\"),\n comp = str_extract(var, \":\\\\d\"),\n comp = str_replace(comp, \":\", \"1v\"),\n ) %>%\n select(-var) %>%\n ## NOTE consider transfering below to rmd for table pres\n pivot_wider(\n names_from = comp,\n values_from = c(exp_b, x2_5, x97_5, pr_z)\n ) %>%\n ## select(name, ends_with(c(\"1v1\"))) %>%\n mutate(\n across(starts_with(\"x\"), ~ round(.x, 2)),\n across(starts_with(\"pr_\"), ~ round(.x, 3))\n )\n\n.fu_ml_1 %>%\n print(n = Inf)\n\n## CREATE FU GRAPHS\n## -----------------------------------------------------------------------------\n## continuous\n.fu_z <-\n .fu_ %>%\n mutate(\n across(\n c(\n pb, tb, bhs, ac,\n dep, anx, str,\n isas_days,\n is_aff_reg, is_ant_dis, is_ant_sui, is_slf_pun,\n dsh_n,\n isas_intra, isas_social,\n si,\n total_session\n ),\n ## where(is.numeric),\n ~ scale(.x),\n .names = \"{.col}_z\"\n )\n ) %>%\n select(c, contains(\"_z\")) %>%\n mutate(\n across(\n everything(),\n ~ str_remove(.x, \"_z\")\n ),\n across(\n everything(), as.double\n )\n ) %>%\n rename_with(~str_remove(., \"_z\"))\n\n## cat\n.fu_plot_cat_names <-\n c(\"DSH[hx]\", \"SA[hx]\", \"SI[hx]\", \"TRM[hx]\", \"FAM\", \"REL\", \"FEM\")\n\nfu_plot_cat <-\n .fu_ %>%\n select(\n c, sex,\n living, relation, dsh_stop, dsh_alone,\n hxsa, hxdsh, hxsuicide, hxtrauma,\n ) %>%\n group_by(c) %>%\n summarise(\n n = n(),\n sex_fem =\n length(sex[sex == \"Female\"]) /\n length(sex[!is.na(sex)]) * 100,\n liv_fam =\n length(living[living == \"family\"]) /\n length(living[!is.na(living)]) * 100,\n rel_yes =\n length(relation[str_detect(relation, \"relation\")]) /\n length(relation[!is.na(relation)]) * 100,\n across(starts_with(c(\"hx\", \"past\")),\n ~ length(.x[str_detect(.x, \"yes\")]) /\n length(.x[!is.na(.x)]) * 100\n )\n ) %>%\n pivot_longer(\n -c\n ) %>%\n filter(name != \"n\") %>%\n ggplot(\n aes(\n x = name,\n y = value,\n linetype = factor(c),\n group = factor(c),\n fill = factor(c),\n shape = factor(c)\n )\n ) +\n geom_bar(position = position_dodge(), stat=\"identity\") +\n theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ggtitle(\"Significant caterogical measures\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n scale_x_discrete(label = parse(text = .fu_plot_cat_names))\n\n## http://www.cookbook-r.com/Graphs/Plotting_means_and_error_bars_(ggplot2)/\n\npd <-\n position_dodge(0.01) # move them .05 to the left and right\n\n.lpa_plot_cont_names <-\n c(\n \"PB\", \"TB\", \"BHS\", \"AC\"\n )\n\nlpa_plot_cont <-\n .fu_z %>%\n pivot_longer(\n cols = -c\n ) %>%\n filter(name %in% c(\"pb\", \"tb\", \"bhs\", \"ac\")) %>%\n group_by(name, c) %>%\n summarise(\n across(\n everything(),\n list(\n N = ~ length(.x),\n mean = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n se = ~ sd(.x, na.rm = TRUE) / sqrt(length(.x)),\n ci = ~ (sd(.x, na.rm = TRUE) / sqrt(length(.x))) * qt(0.95/2 + .5, length(.x) - 1)\n ),\n .names = \"{.col}.{.fn}\")\n ) %>%\n ungroup() %>%\n group_by(name) %>%\n mutate(\n ID = case_when(\n name == \"pb\" ~ 1,\n name == \"tb\" ~ 2,\n name == \"bhs\" ~ 3,\n name == \"ac\" ~ 4\n ),\n .before = name,\n name = factor(\n name,\n levels = c(\n \"pb\", \"tb\", \"bhs\", \"ac\"\n )\n )\n ) %>%\n ungroup() %>%\n ggplot(\n aes(\n x = name,\n y = value.mean,\n colour = factor(c),\n shape = factor(c),\n linetype = factor(c),\n )\n ) +\n geom_errorbar(\n aes(\n ymin = value.mean - value.se,\n ymax = value.mean + value.se\n ), width = .1, position = pd\n ) +\n geom_point(position = pd) +\n geom_line(\n aes(\n x = ID,\n y = value.mean,\n colour = as.factor(c)\n ),\n position = pd) +\n theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ## ggtitle(\"\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n scale_x_discrete(label = .lpa_plot_cont_names)\nlpa_plot_cont\n\nggsave(\n \"./figures/lpa_plot_cont.png\",\n lpa_plot_cont\n)\n\nggsave(\n \"./figures/lpa_plot_cont.pdf\",\n lpa_plot_cont\n)\n\n.fu_plot_cont_names <-\n c(\n ## \"PB\", \"TB\", \"BHS\", \"AC\",\n \"DEP\", \"ANX\", \"STR\", \"SI\",\n ## \"DAYS\",\n \"METH\", \"INTRA\",\n \"AR\", \"AD\", \"AS\", \"SP\",\n \"INTER\"\n ## \"SESSION\"\n )\n\nfu_plot_cont <-\n .fu_z %>%\n pivot_longer(\n cols = -c\n ) %>%\n filter(!name %in% c(\"siss_n\", \"pb\", \"tb\", \"bhs\", \"ac\", \"isas_days\", \"total_session\")) %>%\n group_by(name, c) %>%\n summarise(\n across(\n everything(),\n list(\n N = ~ length(.x),\n mean = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n se = ~ sd(.x, na.rm = TRUE) / sqrt(length(.x)),\n ci = ~ (sd(.x, na.rm = TRUE) / sqrt(length(.x))) * qt(0.95/2 + .5, length(.x) - 1)\n ),\n .names = \"{.col}.{.fn}\")\n ) %>%\n ## print(n = Inf) %>%\n ungroup() %>%\n group_by(name) %>%\n mutate(\n ID = case_when(\n name == \"dep\" ~ 1,\n name == \"anx\" ~ 2,\n name == \"str\" ~ 3,\n name == \"si\" ~ 4,\n name == \"dsh_n\" ~ 5,\n name == \"isas_intra\" ~ 6,\n name == \"is_aff_reg\" ~ 7,\n name == \"is_ant_dis\" ~ 8,\n name == \"is_ant_sui\" ~ 9,\n name == \"is_slf_pun\" ~ 10,\n name == \"isas_social\" ~ 11\n ),\n .before = name,\n name = factor(\n name,\n levels = c(\n ## \"pb\", \"tb\", \"bhs\", \"ac\",\n \"dep\", \"anx\", \"str\", \"si\",\n ## \"siss_n\",\n ## \"isas_days\",\n \"dsh_n\", \"isas_intra\",\n \"is_aff_reg\", \"is_ant_dis\", \"is_ant_sui\", \"is_slf_pun\",\n \"isas_social\"\n ## \"total_session\"\n )\n )\n ) %>%\n ungroup() %>%\n ## filter(name %!in% c(\"pb\", \"tb\", \"bhs\", \"ac\")) %>%\n ggplot(\n aes(\n x = name,\n y = value.mean,\n colour = factor(c),\n shape = factor(c),\n linetype = factor(c)\n )\n ) +\n ## geom_bar(position=position_dodge(), stat=\"identity\") +\n geom_errorbar(\n aes(\n ymin = value.mean - value.se,\n ymax = value.mean + value.se\n ), width = .1, position = pd\n ) +\n geom_point(\n position = pd,\n ) +\n geom_line(\n aes(\n x = ID,\n y = value.mean,\n ),\n position = pd\n ) +\n labs(colour = \"Class\", shape = \"Class\", linetype = \"Class\") +\n ## scale_colour_viridis(\n ## discrete = TRUE,\n ## alpha = 0.7,\n ## )\n ## theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ggtitle(\"Significant continous measures\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n scale_x_discrete(label = .fu_plot_cont_names)\n\ndesign = \"\n112\n112\n112\n332\n332\n\"\n\n.lca_plot_names <-\n ## c(\"Capability\", \"Hopelessness\", \"Burdensomeness\", \"Belongingness\")\n c(\"AC\", \"BHS\", \"PB\", \"TB\")\n\nlpa_fig_plot <-\n lpa_fig +\n theme_ipsum() +\n labs(\n colour = \"Class\",\n shape = \"Class\",\n linetype = \"Class\",\n fill = \"Class\"\n ) +\n ## ggtitle(\"\") +\n xlab(\"\") +\n ## TODO correct for z scale these limits and labels\n ylab(\"\") +\n ## ylim(\"\") +\n scale_x_discrete(label = .lca_plot_names) +\n theme(\n ## legend.position = \"bottom\",\n text = element_text(size = 7),\n title = element_text(size = 8),\n )\n\nlpa_fig_plot +\n theme(\n ## legend.position = \"bottom\",\n text = element_text(size = 7),\n title = element_text(size = 8),\n )\n\nggsave(\n \"./figures/lpa_fig_plot.png\",\n lpa_fig_plot\n)\n\nlpa_combi_plot <-\n ## lpa_plot_cont +\n fu_plot_cont +\n fu_plot_cat +\n plot_layout(\n ## guides = \"collect\",\n ## design = design\n ## heights = c(2, 1)\n widths = c(2, 1)\n ) &\n theme(\n legend.position = \"bottom\",\n text = element_text(size = 7),\n title = element_text(size = 8),\n )\n\n\nggsave(\n \"./figures/lpa_combi_plot.png\",\n lpa_combi_plot\n)\n\n### ============================================================================\n### OUTCOMES\n### ============================================================================\nlpa_t0t1_outc <-\n lpa_fin %>%\n ## filter no class var\n filter(!is.na(c)) %>%\n group_by(ID) %>%\n ## filter first and last session\n filter(session == 1 | final_session == 1) %>%\n ungroup() %>%\n mutate(\n ## replace NA session value as 5 as only 5 sessions of data available\n session = if_else(is.na(session), 5, session),\n ## create conditional variable to drop those with only one session\n drop = if_else(\n session == 1 & final_session == 1,\n 1, 0\n ), .before = ID,\n hxsa = case_when(\n suicide_attempts <= 0 ~ 0,\n suicide_attempts >= 1 ~ 1,\n TRUE ~ NA_real_\n ),\n living = factor(case_when(\n living == 1 ~ \"family\",\n living == 2 ~ \"other\",\n living == 3 ~ \"other\",\n living == 4 ~ \"other\",\n TRUE ~ NA_character_\n )),\n relation = factor(case_when(\n relation == 1 ~ \"single/separated/divorced\",\n relation == 2 ~ \"relationship/married/defacto\",\n relation == 3 ~ \"relationship/married/defacto\",\n relation == 4 ~ \"single/separated/divorced\",\n TRUE ~ NA_character_\n )),\n sex = factor(if_else(sex == 2, NA_character_, as.character(sex)))\n ) %>%\n ## select those with more than one session\n filter(drop == 0) %>%\n group_by(ID) %>%\n mutate(age = min(age)) %>%\n ungroup() %>%\n ## names() %>%\n ## select vars for mi\n dplyr::select(\n c, ID, age, sex, session, final_session,\n planned_final, pb, ac, bhs, tb, dep, anx, str,\n isas_intra, isas_social,\n relation, living,\n hxsa, hxsuicide, hxdsh, hxtrauma,\n si\n ) %>%\n group_by(ID) %>%\n mutate(\n cov_sess = max(session), .after = final_session,\n ## hold AC as constant at time 1 to include as a covariate\n ac = ac[session == 1]\n ) %>%\n ungroup() %>%\n group_by(final_session) %>%\n mutate(\n across(\n c(\n si, pb, tb, bhs, dep, anx, str, ac, isas_social, isas_intra,\n age, cov_sess, hxsa, hxsuicide, hxdsh, hxtrauma\n ),\n ~ scale(.x),\n .names = \"{.col}_z\"\n ),\n ) %>%\n ungroup() %>%\n mutate(\n fin_ses_c = with(\n .,\n factor(final_session):factor(c)\n ),\n .after = c,\n across(\n c(c, sex, ID, final_session, planned_final),\n ~ as.factor(.x)\n ),\n )\n\nlpa_t0t1_outc_w <-\n lpa_t0t1_outc %>%\n dplyr::select(-fin_ses_c) %>%\n panelr::panel_data(\n id = ID,\n wave = final_session,\n ) %>%\n panelr::widen_panel(\n separator = \"_\",\n ignore.attributes = FALSE,\n varying = NULL) %>%\n ungroup() %>%\n ## need to convert to data frame as panelr creates a list object\n as.data.frame() %>%\n ## use age at intake, drop secon ac as not routinely collected\n dplyr::select(-planned_final_0)\n\n\nlpa_t0t1_outc %>% str()\nlpa_t0t1_outc_w %>% str()\n\n\nlpa_t0t1_outc_hisi <-\n lpa_t0t1_outc_w_hisi %>%\n select(!contains(\"_z\")) %>%\n panelr::long_panel(\n prefix = \"_\",\n periods = c(\"0\", \"1\"),\n id = \"ID\",\n wave = \"final_session\",\n label_location = \"end\"\n ) %>%\n ungroup() %>%\n mutate(\n c2 = case_when(\n c == \"2\" ~ \"lomi\",\n c == \"1\" ~ \"lomi\",\n c == \"3\" ~ \"high\",\n TRUE ~ NA_character_\n ), .after = c,\n c2 = as.factor(c2),\n ) %>%\n mutate(\n c3 = case_when(\n c == \"2\" ~ \"2\",\n c == \"3\" ~ \"3\",\n TRUE ~ NA_character_\n ), .after = c2,\n c3 = as.factor(c3),\n ) %>%\n mutate(fin_ses_c = with(., c:final_session), .after = c3) %>%\n mutate(fin_ses_c2 = with(., c2:final_session), .after = fin_ses_c) %>%\n mutate(fin_ses_c3 = with(., c3:final_session), .after = fin_ses_c2) %>%\n mutate(\n across(c(pb, ac, tb, bhs, dep, anx, str, si),\n ~scale(.x),\n .names = \"{.col}_z\"\n )\n )\n\nlpa_t0t1_outc %>%\n filter(final_session == 0) %>%\n group_by(c) %>%\n summarise(\n full = n_distinct(ID),\n hisi = n_distinct(ID[si > 8])\n ) %>%\n adorn_totals(\"row\")\n\nlpa_t0t1_outc_hisi %>%\n group_by(c) %>%\n summarise(\n n = n_distinct(ID)\n ) %>%\n adorn_totals(\"row\")\n\nlpa_t0t1_outc_w_hisi %>%\n select(ID, age, c, starts_with(c(\"pb\", \"si\", \"ac\", \"tb\", \"dep\"))) %>%\n .summary_by(c)\n\n## DISTRIBUTION PLOTS OF ANOVA VARS\n## -----------------------------------------------------------------------------\nplot_iter <-\n lpa_t0t1_outc %>%\n dplyr::select(\n ## ends_with(\"_z\") & !starts_with(\"isas\"),\n si, pb, tb, bhs, dep, anx,\n str, age, sex, cov_sess\n ) %>%\n names() %>%\n purrr::set_names()\n\ndens_fun <-\n function(dat, x, filly, facet) {\n ggplot(dat, aes(x = .data[[x]], fill = .data[[filly]])) +\n geom_density(alpha = 0.3) +\n facet_grid(.data[[facet]] ~ .) +\n theme_bw()\n }\n\n.plots_lpa_long <-\n map(plot_iter, ~ dens_fun(lpa_t0t1_outc, .x, \"c\", \"final_session\"))\n\nplots_lpa_long <-\n wrap_plots(.plots_lpa_long, guides = \"collect\")\n\nggsave(\n \"./figures/plots_lpa_long.png\",\n plots_lpa_long\n\n\n## PLOT CHANGE OVER TIME\n## =============================================================================\n## TODO get plots\n\n## LME4 PACKAGE\n## =============================================================================\n## Mixed Effects ANOVA\nstep_data <-\n lpa_t0t1_outc\n\nno_na_step_data <-\n na.omit(\n step_data[c(\n \"si\", \"pb\", \"bhs\", \"tb\", \"ac\", \"sex\",\n \"age\",\n \"dep\", \"anx\", \"str\",\n \"isas_social\", \"isas_intra\", \"hxsa\", \"hxdsh\", \"hxtrauma\", \"hxsuicide\",\n \"cov_sess\", \"c\",\n \"si_z\", \"pb_z\", \"bhs_z\", \"tb_z\", \"ac_z\", \"sex_z\",\n \"final_session\", \"age_z\",\n \"dep_z\", \"anx_z\", \"str_z\",\n \"isas_social_z\", \"isas_intra_z\", \"hxsa_z\", \"hxdsh_z\", \"hxtrauma_z\", \"hxsuicide_z\",\n \"cov_sess_z\"\n )])\n\n\n\n.tt0 <-\n no_na_step_data %>%\n filter(final_session == 0)\n\n.tt1 <-\n no_na_step_data %>%\n filter(final_session == 1)\n\n.tt0_lmer <-\n lmer(\n si ~\n pb + tb + bhs + dep + anx + str +\n ac + cov_sess + age + sex +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n isas_social + isas_intra +\n (1 | c),\n data = .tt0,\n na.action = na.exclude,\n REML = FALSE\n )\n\n.tt1_lmer <-\n lmer(\n si ~\n pb + tb + bhs + dep + anx + str +\n ac + cov_sess + age + sex +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n isas_social + isas_intra +\n (1 | c),\n data = .tt1,\n na.action = na.exclude,\n REML = FALSE\n )\n\n\ntt0_fixmodel <-\n lm(formula(.tt0_lmer, fixed.only=TRUE),\n data=eval(getCall(.tt0_lmer)$data))\n\nstep_model_tt0 <-\n stats::step(\n tt0_fixmodel,\n direction = \"both\",\n trace = 0,\n test = \"F\"\n )\n\nstep_model_tt0$anova\nstep_model_tt0$coefficients\nsummary(step_model_tt0)\n\ntt1_fixmodel <-\n lm(formula(.tt1_lmer, fixed.only=TRUE),\n data=eval(getCall(.tt1_lmer)$data))\n\nstep_model_tt1 <-\n stats::step(\n tt1_fixmodel,\n direction = \"both\",\n trace = 0,\n test = \"F\"\n )\n\nstep_model_tt1$anova\nstep_model_tt1$coefficients\nsummary(step_model_tt1)\n\n\nlpa_out_ni_step <-\n lmer(\n si ~\n pb + tb + bhs +\n dep + anx + str +\n ac + cov_sess + age + sex +\n isas_social + isas_intra +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n (1 | final_session/c),\n ## (1 | final_session/c),\n data = no_na_step_data,\n ## na.action = na.omit,\n ## control = lmerControl(optimizer = \"Nelder_Mead\"),\n REML = FALSE # use ML estimation\n )\n\nfixmodel <-\n lm(formula(lpa_out_ni_step, fixed.only=TRUE),\n data=eval(getCall(lpa_out_ni_step)$data))\n\nstep_model <-\n stats::step(\n fixmodel,\n direction = \"both\",\n trace = 0,\n test = \"F\"\n )\n\nstep_model$anova\nstep_model$coefficients\nsummary(step_model)\n\nlpa_out_step <-\n lmer(\n si_z ~\n pb_z + dep_z + ac_z + cov_sess_z + age_z +\n isas_social_z + hxsa_z +\n (1 | final_session/c),\n ## (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\n## NOTE lmerTest used to add t and p values\ncoef(summary(lpa_out_step, correlation = FALSE))\ncoef(summary(as(lpa_out_step, \"merModLmerTest\")))\nconfint(lpa_out_step, method = \"Wald\", oldNames = FALSE)\nsummary(as(lpa_out_step, \"merModLmerTest\"))\n\nlpa_out_ni <-\n lmer(\n si ~\n pb + tb + bhs + dep +\n ac + cov_sess + age + sex +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(lpa_out_ni)\n\n## NOTE lmerTest used to add t and p values\ncoef(summary(lpa_out_ni, correlation = FALSE))\ncoef(summary(as(lpa_out_ni, \"merModLmerTest\")))\nconfint(lpa_out_ni, method = \"Wald\", oldNames = FALSE)\nsummary(as(lpa_out_ni, \"merModLmerTest\"))\nvcov(lpa_out_ni)\n\n## calculate ICC for lmer model object\nicc_lmer <- function(m){\n vc <- as.data.frame((VarCorr(m)))\n l <- vc$vcov\n tibble(grp=vc$grp, icc=sapply(l, function(x){x/sum(l)}))\n}\n\nicc_lmer(lpa_out_ipt)\nVarCorr(lpa_out_ipt)\n\nrcompanion::nagelkerke(\n lpa_out_ni,\n lpa_out_step\n)\n\n## fit with all optimizers\nall_lpa <-\n allFit(lpa_out_ni)\nss <-\n summary(all_lpa)\nss$which.OK ## logical vector: which optimizers worked?\n## the other components only contain values for the optimizers that worked\nss$llik ## vector of log-likelihoods\nss$fixef ## table of fixed effects\nss$sdcor ## table of random effect SDs and correlations\nss$theta ## table of random effects parameters, Cholesky scale\n\nlpa_out_ni_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ni, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goews with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_out_ni_tidy %>%\n print(n = Inf)\n\n## by Ben Bolker\n## stackoverflow.com/questions/25142901/standardized-coefficients-for-lmer-model\nstdCoef.merMod <-\n function(object) {\n sdy <- sd(getME(object,\"y\"))\n sdx <- apply(getME(object,\"X\"), 2, sd)\n sc <- fixef(object)*sdx/sdy\n se.fixef <- coef(summary(object))[,\"Std. Error\"]\n se <- se.fixef*sdx/sdy\n return(data.frame(stdcoef=sc, stdse=se))\n}\n\nlpa_out_ni %>%\n stdCoef.merMod()\n\nlpa_out_ipt_hisi <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session), # does not provide singularity error\n ## (1 | final_session/c), # same results as above, with singularity error\n ## (1 | final_session/c2), # same results as above, with singularity error\n ## (1 | final_session/c3), # only pb sig, with ac almost sig, signul error\n ## data = lpa_t0t1_outc_hisi[!is.na(lpa_t0t1_outc_hisi$c3), ], # c3 data\n data = lpa_t0t1_outc_hisi,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\n\nsummary(as(lpa_out_ipt_hisi, \"merModLmerTest\"))\n## singular fit potentially becuase of low cell numbers in c2\n## recommendation when singular fit to simplify the model, therefore nesting removed\nVarCorr(lpa_out_ipt_hisi)\n\n\nlpa_out_ipt_hisi_plot <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n ac_z +\n ## cov_sess + age + sex +\n (1 | final_session), # does not provide singularity error\n ## (1 | final_session/c), # same results as above, with singularity error\n ## (1 | final_session/c2), # same results as above, with singularity error\n ## (1 | final_session/c3), # only pb sig, with ac almost sig, signul error\n ## data = lpa_t0t1_outc_hisi[!is.na(lpa_t0t1_outc_hisi$c3), ], # c3 data\n data = lpa_t0t1_outc_hisi,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(as(lpa_out_ipt_hisi_plot, \"merModLmerTest\"))\n\nlpa_int_plot_hihsi <-\n interact_plot(\n lpa_out_ipt_hisi,\n pred = pb_z,\n modx = ac_z,\n mod2 = tb_z,\n plot.points = FALSE\n ) +\n theme(\n text = element_text(size = 9),\n ) +\n labs(\n title = \"Meets SI Cut-off\"\n )\n\n\n\nlpa_int_plot2 <-\n interact_plot(\n lpa_out_ipt,\n pred = pb_z,\n modx = ac_z,\n plot.points = FALSE\n ) +\n theme(\n text = element_text(size = 9),\n )+\n labs(\n title = \"SI ~ PB * AC\"\n )\n\nlpa_int_plot <-\n interact_plot(\n lpa_out_ipt,\n pred = pb_z,\n modx = ac_z,\n plot.points = FALSE\n ) +\n theme(\n legend.position = \"none\",\n text = element_text(size = 9),\n )+\n labs(\n title = \"Full Sample\"\n )\n\nlpa_int_plot_combi <-\n lpa_int_plot +\n lpa_int_plot_hihsi &\n plot_layout(\n nrow = 2,\n ) &\n theme(\n text = element_text(size = 9),\n title = element_text(size = 9),\n )\n\n\nggsave(\n \"./figures/lpa_int_plot.png\",\n lpa_int_plot2\n)\n\nggsave(\n \"./figures/lpa_int_plot_combi.png\",\n lpa_int_plot_combi\n)\n\nlpa_t0t1_outc_hisi %>%\n group_by(final_session) %>%\n summarise(\n full = n_distinct(ID),\n hisi = n_distinct(ID[si > 8 & final_session %in% \"0\"]),\n ) %>%\n adorn_totals(\"row\") %>%\n mutate(\n prc = round(hisi/full * 100, 1)\n )\n\nlpa_main_fx_c_ipt_hisi <-\n lpa_t0t1_outc_hisi %>%\n ## split(.$final_session) %>%\n split(.$c) %>%\n map(~\n ## lm(\n lmer(\n si ~\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n ## map(~ summary(.)) %>%\n map(~ tidy(as(., \"merModLmerTest\"), conf.method = \"profile\", conf.int = TRUE)) %>%\n bind_rows(., .id = \"var\")\n\nlpa_main_fx_c_ipt_hisi %>%\n filter(p.value < 0.5) %>%\n print(n = Inf)\n\nlpa_t0t1_outc %>%\n group_by(c) %>%\n summarise(\n n = n_distinct(ID)\n )\n\n## START MULTILEVEL MODEL STUFF\n## -----------------------------------------------------------------------------\nlpa_out_ipt <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(as(lpa_out_ipt, \"merModLmerTest\"))\n\nlpa_out_ipt_std <-\n lpa_out_ipt %>%\n stdCoef.merMod()\n\nlpa_out_ipt_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ipt, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goes with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_out_ipt_tidy %>%\n print(n = Inf)\n\nlpa_out_ipt_rs <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (c | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\n## complare with added slope, random intercept enough\nanova(lpa_out_ipt_rs, lpa_out_ipt)\n\nsummary(as(lpa_out_ipt_rs, \"merModLmerTest\"))\n\nlpa_out_ipt_rs_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ipt_rs, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goes with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_main_fx_c_ipt <-\n lpa_t0t1_outc %>%\n split(.$c) %>%\n map(~\n lmer(\n si ~\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n map(~ tidy(as(., \"merModLmerTest\"), conf.method = \"profile\", conf.int = TRUE)) %>%\n bind_rows(., .id = \"var\")\n\nlpa_main_fx_c_ipt %>%\n filter(p.value < 0.05) %>%\n print(n = Inf)\n\nprint(VarCorr(lpa_out_ipt), comp = c(\"Variance\", \"Std.Dev\"), digits = 2)\nVarCorr(lpa_out_ipt)\n\nlpa_out_ipt_all <-\n bind_rows(\n lpa_out_ipt_tidy,\n lpa_main_fx_c_ipt\n ) %>%\n clean_names()\n\n\nlpa_int_pb_ac <-\n interact_plot(\n lpa_out_ipt,\n pred = pb_z,\n modx = ac_z,\n ## cluster = \"final_session/c\",\n modx.values = \"terciles\",\n plot.points = FALSE\n )\n\nlpa_int_tb_ac <-\n interact_plot(\n lpa_out_ipt,\n pred = tb_z,\n modx = ac_z,\n plot.points = FALSE\n )\n\n## ENDS MULTILEVEL MODEL STUFF\n## -----------------------------------------------------------------------------\n\nlpa_main_fx_c <-\n lpa_t0t1_outc %>%\n split(.$c) %>%\n map(~\n lmer(\n si ~\n pb + tb + bhs + dep +\n ac + cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n map(~ tidy(as(., \"merModLmerTest\"), conf.method = \"profile\", conf.int = TRUE)) %>%\n bind_rows(., .id = \"var\")\n\nlpa_out <-\n bind_rows(\n lpa_out_ni_tidy,\n lpa_main_fx_c\n )\n\nlpa_out_tbl <-\n lpa_out %>%\n select(var, effect, group, term, estimate, std.error, starts_with(\"conf\"), p.value) %>%\n pivot_wider(\n ## -c(effect:group),\n names_from = var,\n values_from = c(estimate, std.error, conf.low, conf.high, p.value)\n ) %>%\n filter(is.na(group)) %>%\n select(term, ends_with(c(\"full\", \"_1\", \"_2\", \"_3\"))) %>%\n clean_names()\n\nlpa_out_tbl %>%\n filter(if_any(starts_with(\"p_value\"), ~.x < 0.05)) %>%\n select(term, starts_with(\"p_value\"))\n\nlpa_main_fx_c %>%\n filter(p.value < 0.05) %>%\n print(n = Inf)\n\nlpa_t0t1_outc %>%\n split(.$c) %>%\n map(~\n lmer(\n si ~\n pb + tb + bhs + dep +\n ac + cov_sess + age + sex +\n (1 | final_session),\n data = .x,\n na.action = na.exclude,\n REML = FALSE\n )\n ) %>%\n map(~ stdCoef.merMod(.))\n\n## test for interactions\nlpa_out_i <-\n lmer(\n si ~\n pb + tb + bhs + dep + anx + str +\n ac + cov_sess + age + sex +\n hxsa + hxsuicide + hxdsh + hxtrauma +\n isas_social + isas_intra +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n control = lmerControl(optimizer = \"Nelder_Mead\"),\n REML = FALSE # use ML estimation\n )\n\n## singular fit arises\nall_lpa_i <- allFit(lpa_out_i)\nsummary(all_lpa_i)\n\nanova(\n lpa_out_i,\n lpa_out_step\n)\n## null model with random effects\nlpa_out_null_lmer <-\n lmer(\n si_z ~ 1 +\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n REML = FALSE\n )\n\n## null model with no random effects\nlpa_out_null_gls <-\n nlme::gls(\n si_z ~ 1,\n data = lpa_t0t1_outc,\n method = \"ML\",\n na.action = na.exclude\n )\n\nrcompanion::nagelkerke(\n lpa_out_step,\n lpa_out_null_lmer,\n )\n\nrcompanion::nagelkerke(\n lpa_out_step,\n lpa_out_null_gls,\n )\n\n.lpa_t_outc <-\n lpa_t0t1_outc %>%\n dplyr::select(si, pb, tb, bhs, dep) %>%\n map( ~ tidy(pairwise.t.test(., lpa_t0t1_outc$fin_ses_c, p.adjust.method = \"bonf\"))) %>%\n unlist(recursive = FALSE) %>%\n as_tibble() %>%\n dplyr::select(pb.group1, pb.group2, ends_with(\"value\")) %>%\n rename(group_1 = 1, group_2 = 2) %>%\n janitor::clean_names()\n.lpa_t_outc\n\n### ============================================================================\n### RELIABLE CHANGE\n### ============================================================================\n\n## RELIABLE CHANGE MSSI CASENESS\n## =============================================================================\n## reliable change calculations\n## pb rci\n\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"pb_0\",\n post = \"pb_1\",\n group = \"c\",\n reliability = meas_rel$pb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\npb_rci_hisi <- JTRCIdf\npb_rci_hisi_fig <-\n plotRCI(\n pb_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Burdensomeness\"\n )\npb_rci_hisi_tbl <-\n tableRCI(\n pb_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"pb\", .before = everything())\n\n\n## tb rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"tb_0\",\n post = \"tb_1\",\n group = \"c\",\n reliability = meas_rel$tb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ntb_rci_hisi <- JTRCIdf\ntb_rci_hisi_fig <-\n plotRCI(\n tb_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Belongingness\"\n )\ntb_rci_hisi_tbl <-\n tableRCI(\n tb_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"tb\", .before = everything())\n\n## bhs rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"bhs_0\",\n post = \"bhs_1\",\n group = \"c\",\n reliability = (meas_rel$bhs_sf_rel$total$raw_alpha + meas_rel$bhs_20_rel$total$raw_alpha) / 2,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nbhs_rci_hisi <- JTRCIdf\nbhs_rci_hisi_fig <-\n plotRCI(\n bhs_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Hopelessness\"\n )\nbhs_rci_hisi_tbl <-\n tableRCI(\n bhs_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"bhs\", .before = everything())\n\n## dep rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"dep_0\",\n post = \"dep_1\",\n group = \"c\",\n reliability = meas_rel$dep_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ndep_rci_hisi <- JTRCIdf\ndep_rci_hisi_fig <-\n plotRCI(\n dep_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Depression\"\n )\ndep_rci_hisi_tbl <-\n tableRCI(\n dep_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"dep\", .before = everything())\n\n## anx rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"anx_0\",\n post = \"anx_1\",\n group = \"c\",\n reliability = meas_rel$anx_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nanx_rci_hisi <- JTRCIdf\nanx_rci_hisi_fig <-\n plotRCI(\n anx_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Anxiety\"\n )\nanx_rci_hisi_tbl <-\n tableRCI(\n anx_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"anx\", .before = everything())\n\n## str rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"str_0\",\n post = \"str_1\",\n group = \"c\",\n reliability = meas_rel$str_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nstr_rci_hisi <- JTRCIdf\nstr_rci_hisi_fig <-\n plotRCI(\n str_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Stress\"\n )\nstr_rci_hisi_tbl <-\n tableRCI(\n str_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"str\", .before = everything())\n\n## si rci_hisi\nJTRCI(\n lpa_t0t1_outc_w_hisi,\n ppid = \"ID\",\n pre = \"si_0\",\n post = \"si_1\",\n group = \"c\",\n reliability = meas_rel$si_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nsi_rci_hisi <- JTRCIdf\nsi_rci_hisi_fig <-\n plotRCI(\n si_rci_hisi,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n plottitle = \"Suicidal Ideation\"\n )\nsi_rci_hisi_tbl <-\n tableRCI(\n si_rci_hisi,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"si\", .before = everything())\n\n.rci_hisi_tbl <-\n bind_rows(\n si_rci_hisi_tbl,\n pb_rci_hisi_tbl,\n tb_rci_hisi_tbl,\n bhs_rci_hisi_tbl,\n dep_rci_hisi_tbl\n ) %>%\n rename(rci_hisi = 2) %>%\n mutate(rci_hisi = case_when(\n str_detect(rci_hisi, \"deteriorated\") ~ \"det\",\n str_detect(rci_hisi, \"no reliable\") ~ \"noch\",\n str_detect(rci_hisi, \"improved\") ~ \"imp\"\n ),\n across(\n where(is.character) & -c(var, rci_hisi),\n as.numeric\n ),\n ) %>%\n group_by(var) %>%\n mutate(\n across(\n where(is.numeric),\n ~ round(.x / sum(.x) * 100, 1),\n .names = \"{.col}_pr\"\n )\n ) %>%\n ungroup()\n\n.rci_hisi_tbl_chisq <-\n list(\n si = si_rci_hisi_tbl,\n pb = pb_rci_hisi_tbl,\n tb = tb_rci_hisi_tbl,\n bhs = bhs_rci_hisi_tbl,\n dep = dep_rci_hisi_tbl\n ) %>%\n map(`[`, 3:5) %>%\n map(~ .x %>% mutate(across(where(is.character), as.numeric))) %>%\n map(~ .x %>% stats::chisq.test(table(.))) %>%\n map(`[`, c(\"p.value\", \"statistic\", \"stdres\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ .x %>%\n rename(\n stdres_1 = 3,\n stdres_2 = 4,\n stdres_3 = 5\n )) %>%\n ## map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n mutate(\n across(\n starts_with(\"stdres_\"),\n ~ if_else(abs(.x) > bonf, \"yes\", \"no\"),\n ## .names = \"{.col}_sig\"\n .names = \"{.col}_sig\"\n )\n ) %>%\n mutate(\n rci_hisi =\n rep(c(\"det\", \"noch\", \"imp\"), 5),\n .after = var\n )\n\nrci_hisi_tbl <-\n .rci_hisi_tbl %>%\n full_join(.rci_hisi_tbl_chisq)\n\n.rci_full_kbl <-\n .rci_tbl %>%\n select(var:x3) %>%\n pivot_wider(\n names_from = var,\n values_from = x1:x3\n ) %>%\n ## rownames_to_column() %>%\n pivot_longer(\n -rci,\n 'variable',\n 'value'\n ) %>%\n pivot_wider(variable, rci) %>%\n separate(variable, into = c(\"c\", \"var\"), sep = \"_\") %>%\n arrange(var, c) %>%\n rowwise() %>%\n mutate(\n N = sum(c_across(det:imp)), .after = var,\n across(\n det:imp,\n ~ round(.x / N * 100, 1)\n )\n ) %>%\n ungroup() %>%\n mutate(c = str_remove(c, \"x\")) %>%\n select(var, c, N, imp, noch, det)\n\n\nrci_full_kbl <-\n .rci_tbl_chisq %>%\n group_by(var) %>%\n slice(1) %>%\n ungroup() %>%\n select(var, statistic, p_value) %>%\n right_join(.rci_full_kbl) %>%\n select(var, c:det, statistic, p_value)\n\n.rci_hisi_kbl <-\n .rci_hisi_tbl %>%\n select(var:x3) %>%\n pivot_wider(\n names_from = var,\n values_from = x1:x3\n ) %>%\n ## rownames_to_column() %>%\n pivot_longer(\n -rci_hisi,\n 'variable',\n 'value'\n ) %>%\n pivot_wider(variable, rci_hisi) %>%\n separate(variable, into = c(\"c\", \"var\"), sep = \"_\") %>%\n arrange(var, c) %>%\n rowwise() %>%\n mutate(\n N = sum(c_across(det:imp)), .after = var,\n across(\n det:imp,\n ~ round(.x / N * 100, 1)\n )\n ) %>%\n ungroup() %>%\n mutate(c = str_remove(c, \"x\")) %>%\n select(var, c, N, imp, noch, det)\n\nrci_hisi_kbl <-\n .rci_hisi_tbl_chisq %>%\n group_by(var) %>%\n slice(1) %>%\n ungroup() %>%\n select(var, statistic, p_value) %>%\n right_join(.rci_hisi_kbl) %>%\n select(var, c:det, statistic, p_value)\n\n\n.rci_kbl_group_diff <-\n readxl::read_xlsx(\n \"./data/tmp/anova_results.xlsx\",\n sheet = \"rci_tbl\"\n ) %>%\n select(sample:rci, stdres_1_sig:group_diff) %>%\n pivot_longer(\n -c(sample, rci, var, group_diff)\n ) %>%\n group_by(sample, var, rci, name) %>%\n slice(1) %>%\n ungroup() %>%\n mutate(\n c = str_extract(name, \"\\\\d\")\n ) %>%\n select(-name) %>%\n pivot_wider(\n names_from = c(sample, rci),\n values_from = group_diff\n ) %>%\n unite(\"group_diff_full\", c(full_imp, full_noch, full_det), na.rm = TRUE, sep = \";<br>\") %>%\n unite(\"group_diff_hisi\", c(hisi_imp, hisi_noch, hisi_det), na.rm = TRUE, sep = \";<br>\") %>%\n mutate(\n across(\n starts_with(\"group_\"),\n ~ if_else(.x == \"\", NA_character_, .x)\n )\n ) %>%\n select(-value) %>%\n group_by(var, c) %>%\n summarise(across(everything(), .coalesce_by_column)) %>%\n ungroup()\n\n\nrci_kbl <-\n list(\n full = rci_full_kbl,\n hisi = rci_hisi_kbl\n ) %>%\n bind_rows(.id = \"vers\") %>%\n pivot_wider(\n names_from = vers,\n values_from = N:p_value\n ) %>%\n full_join(.rci_kbl_group_diff) %>%\n select(\n var, c,\n ends_with(c(\"full\", \"hisi\"))\n ) %>%\n mutate(\n order = case_when(\n var == \"si\" ~ 1,\n var == \"pb\" ~ 2,\n var == \"tb\" ~ 3,\n var == \"bhs\" ~ 4,\n var == \"dep\" ~ 5\n ), .before = 1\n ) %>%\n arrange(order)\n\n\nrci_hisi_combi_plot <-\n pb_rci_hisi_fig +\n tb_rci_hisi_fig +\n bhs_rci_hisi_fig +\n dep_rci_hisi_fig +\n anx_rci_hisi_fig +\n str_rci_hisi_fig +\n si_rci_hisi_fig +\n guide_area() +\n plot_layout(\n nrow = 2,\n guides = \"collect\",\n ) &\n theme(\n plot.title = element_text(size = 10),\n legend.title = element_text(size = 10),\n legend.text = element_text(size = 9),\n axis.title.x = element_text(size = 9),\n axis.title.y = element_text(size = 9),\n axis.text = element_text(size = 7),\n ) &\n labs(\n colour = \"RCI_HISI classification\",\n shape = \"Class\"\n ) &\n scale_fill_brewer(palette = \"Dark2\") &\n scale_colour_brewer(palette = \"Dark2\")\n\nggsave(\n \"./figures/rci_hisi_combi_plot.png\",\n rci_hisi_combi_plot\n)\n\n\n\n## RELIABLE CHANGE\n## =============================================================================\n## reliable change calculations\n## pb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"pb_0\",\n post = \"pb_1\",\n group = \"c\",\n reliability = meas_rel$pb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\npb_rci <- JTRCIdf\npb_rci_fig <-\n plotRCI(\n pb_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Burdensomeness\"\n )\npb_rci_tbl <-\n tableRCI(\n pb_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"pb\", .before = everything())\n\n\n## tb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"tb_0\",\n post = \"tb_1\",\n group = \"c\",\n reliability = meas_rel$tb_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ntb_rci <- JTRCIdf\ntb_rci_fig <-\n plotRCI(\n tb_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Belongingness\"\n )\ntb_rci_tbl <-\n tableRCI(\n tb_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"tb\", .before = everything())\n\n## bhs rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"bhs_0\",\n post = \"bhs_1\",\n group = \"c\",\n reliability = (meas_rel$bhs_sf_rel$total$raw_alpha + meas_rel$bhs_20_rel$total$raw_alpha) / 2,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nbhs_rci <- JTRCIdf\nbhs_rci_fig <-\n plotRCI(\n bhs_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Hopelessness\"\n )\nbhs_rci_tbl <-\n tableRCI(\n bhs_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"bhs\", .before = everything())\n\n## dep rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"dep_0\",\n post = \"dep_1\",\n group = \"c\",\n reliability = meas_rel$dep_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ndep_rci <- JTRCIdf\ndep_rci_fig <-\n plotRCI(\n dep_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Depression\"\n )\ndep_rci_tbl <-\n tableRCI(\n dep_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"dep\", .before = everything())\n\n## anx rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"anx_0\",\n post = \"anx_1\",\n group = \"c\",\n reliability = meas_rel$anx_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nanx_rci <- JTRCIdf\nanx_rci_fig <-\n plotRCI(\n anx_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Anxiety\"\n )\nanx_rci_tbl <-\n tableRCI(\n anx_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"anx\", .before = everything())\n\n## str rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"str_0\",\n post = \"str_1\",\n group = \"c\",\n reliability = meas_rel$str_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nstr_rci <- JTRCIdf\nstr_rci_fig <-\n plotRCI(\n str_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Stress\"\n )\nstr_rci_tbl <-\n tableRCI(\n str_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"str\", .before = everything())\n\n## si rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"si_0\",\n post = \"si_1\",\n group = \"c\",\n reliability = meas_rel$si_rel$total$raw_alpha,\n indextype = \"RCI\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nsi_rci <- JTRCIdf\nsi_rci_fig <-\n plotRCI(\n si_rci,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n plottitle = \"Suicidal Ideation\"\n )\nsi_rci_tbl <-\n tableRCI(\n si_rci,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"si\", .before = everything())\n\n.rci_tbl <-\n bind_rows(\n si_rci_tbl,\n pb_rci_tbl,\n tb_rci_tbl,\n bhs_rci_tbl,\n dep_rci_tbl\n ) %>%\n rename(rci = 2) %>%\n mutate(rci = case_when(\n str_detect(rci, \"deteriorated\") ~ \"det\",\n str_detect(rci, \"no reliable\") ~ \"noch\",\n str_detect(rci, \"improved\") ~ \"imp\"\n ),\n across(\n where(is.character) & -c(var, rci),\n as.numeric\n ),\n ) %>%\n group_by(var) %>%\n mutate(\n across(\n where(is.numeric),\n ~ round(.x / sum(.x) * 100, 1),\n .names = \"{.col}_pr\"\n )\n ) %>%\n ungroup()\n\n.rci_tbl_chisq <-\n list(\n si = si_rci_tbl,\n pb = pb_rci_tbl,\n tb = tb_rci_tbl,\n bhs = bhs_rci_tbl,\n dep = dep_rci_tbl\n ) %>%\n map(`[`, 3:5) %>%\n map(~ .x %>% mutate(across(where(is.character), as.numeric))) %>%\n map(~ .x %>% stats::chisq.test(table(.))) %>%\n map(`[`, c(\"p.value\", \"statistic\", \"stdres\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ .x %>%\n rename(\n stdres_1 = 3,\n stdres_2 = 4,\n stdres_3 = 5\n )) %>%\n ## map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n mutate(\n across(\n starts_with(\"stdres_\"),\n ~ if_else(abs(.x) > bonf, \"yes\", \"no\"),\n ## .names = \"{.col}_sig\"\n .names = \"{.col}_sig\"\n )\n ) %>%\n mutate(\n rci =\n rep(c(\"det\", \"noch\", \"imp\"), 5),\n .after = var\n )\n\nrci_tbl <-\n .rci_tbl %>%\n full_join(.rci_tbl_chisq)\n\nrci_combi_plot <-\n pb_rci_fig +\n tb_rci_fig +\n bhs_rci_fig +\n dep_rci_fig +\n anx_rci_fig +\n str_rci_fig +\n si_rci_fig +\n guide_area() +\n plot_layout(\n nrow = 2,\n guides = \"collect\",\n ) &\n theme(\n plot.title = element_text(size = 10),\n legend.title = element_text(size = 10),\n legend.text = element_text(size = 9),\n axis.title.x = element_text(size = 9),\n axis.title.y = element_text(size = 9),\n axis.text = element_text(size = 7),\n ) &\n labs(\n colour = \"RCI classification\",\n shape = \"Class\"\n ) &\n scale_fill_brewer(palette = \"Dark2\") &\n scale_colour_brewer(palette = \"Dark2\")\n\nggsave(\n \"./figures/rci_combi_plot.png\",\n rci_combi_plot\n)\n\n\n## JT CLIN SIG CH\n## =============================================================================\n## reliable change calculations\n## pb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"pb_0\",\n post = \"pb_1\",\n group = \"c\",\n reliability = meas_rel$pb_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\npb_jt <- JTRCIdf\npb_jt_fig <-\n plotJT(\n pb_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Burdensomeness\"\n )\npb_jt_tbl <-\n tableJT(\n pb_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"pb\", .before = everything())\n\n## tb rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"tb_0\",\n post = \"tb_1\",\n group = \"c\",\n reliability = meas_rel$tb_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ntb_jt <- JTRCIdf\ntb_jt_fig <-\n plotJT(\n tb_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Belongingness\"\n )\ntb_jt_tbl <-\n tableJT(\n tb_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"tb\", .before = everything())\n\n## bhs rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"bhs_0\",\n post = \"bhs_1\",\n group = \"c\",\n reliability = (meas_rel$bhs_sf_rel$total$raw_alpha + meas_rel$bhs_20_rel$total$raw_alpha) / 2,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nbhs_jt <- JTRCIdf\nbhs_jt_fig <-\n plotJT(\n bhs_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Hopelessness\"\n )\nbhs_jt_tbl <-\n tableJT(\n bhs_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"bhs\", .before = everything())\n\n## dep rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"dep_0\",\n post = \"dep_1\",\n group = \"c\",\n reliability = meas_rel$dep_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\ndep_jt <- JTRCIdf\ndep_jt_fig <-\n plotJT(\n dep_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Depression\"\n )\ndep_jt_tbl <-\n tableJT(\n dep_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"dep\", .before = everything())\n\n## anx rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"anx_0\",\n post = \"anx_1\",\n group = \"c\",\n reliability = meas_rel$anx_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nanx_jt <- JTRCIdf\nanx_jt_fig <-\n plotJT(\n anx_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Anxiety\"\n )\nanx_jt_tbl <-\n tableJT(\n anx_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"anx\", .before = everything())\n\n## str rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"str_0\",\n post = \"str_1\",\n group = \"c\",\n reliability = meas_rel$str_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nstr_jt <- JTRCIdf\nstr_jt_fig <-\n plotJT(\n str_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n ## facetplot = TRUE,\n plottitle = \"Stress\"\n )\nstr_jt_tbl <-\n tableJT(\n str_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"str\", .before = everything())\n\n## si rci\nJTRCI(\n lpa_t0t1_outc_w,\n ppid = \"ID\",\n pre = \"si_0\",\n post = \"si_1\",\n group = \"c\",\n reliability = meas_rel$si_rel$total$raw_alpha,\n indextype = \"JT\",\n JTcrit = \"auto\",\n plot = FALSE,\n table = FALSE\n)\nsi_jt <- JTRCIdf\nsi_jt_fig <-\n plotJT(\n si_jt,\n addInfoLegend = \"no\",\n useGroups = TRUE,\n addJitter = FALSE,\n plottitle = \"Suicidal Ideation\"\n )\nsi_jt_tbl <-\n tableJT(\n si_jt,\n useGroups = TRUE\n ) %>%\n clean_names() %>%\n mutate(var = \"si\", .before = everything())\n\njt_tbl <-\n pb_jt_tbl %>%\n full_join(tb_jt_tbl) %>%\n full_join(bhs_jt_tbl) %>%\n full_join(dep_jt_tbl) %>%\n full_join(anx_jt_tbl) %>%\n full_join(str_jt_tbl) %>%\n full_join(si_jt_tbl)\n\n\n.jt_tbl <-\n bind_rows(\n si_jt_tbl,\n pb_jt_tbl,\n tb_jt_tbl,\n bhs_jt_tbl,\n dep_jt_tbl\n ) %>%\n rename(jt = 2) %>%\n mutate(jt = case_when(\n str_detect(jt, \"deteriorated\") ~ \"det\",\n str_detect(jt, \"unchanged\") ~ \"unch\",\n str_detect(jt, \"improved\") ~ \"imp\",\n str_detect(jt, \"non reliab\") ~ \"nr_rec\",\n TRUE ~ \"rec\"\n ),\n across(\n where(is.character) & -c(var, jt),\n as.numeric\n ),\n ) %>%\n group_by(var) %>%\n mutate(\n across(\n where(is.numeric),\n ~ round(.x / sum(.x) * 100, 1),\n .names = \"{.col}_pr\"\n )\n ) %>%\n ungroup()\n\n.jt_tbl_chisq <-\n list(\n si = si_jt_tbl,\n pb = pb_jt_tbl,\n tb = tb_jt_tbl,\n bhs = bhs_jt_tbl,\n dep = dep_jt_tbl\n ) %>%\n map(`[`, 3:5) %>%\n map(~ .x %>% mutate(across(where(is.character), as.numeric))) %>%\n map(~ .x %>% stats::chisq.test(table(.))) %>%\n map(`[`, c(\"p.value\", \"statistic\", \"stdres\")) %>%\n map(~ as.data.frame(.)) %>%\n map(~ .x %>%\n rename(\n stdres_1 = 3,\n stdres_2 = 4,\n stdres_3 = 5\n )) %>%\n ## map(~ pivot_wider(., names_from = 2, values_from = 3)) %>%\n map(~ janitor::clean_names(.)) %>%\n map(.,\n ~ as.data.frame(bind_cols(\n .x,\n ## create bonf corrected cutoff\n bonf = qnorm(1-0.05/(3 * nrow(.x))/2),\n ))\n ) %>%\n bind_rows(., .id = \"var\") %>%\n mutate(\n across(\n starts_with(\"stdres_\"),\n ~ if_else(abs(.x) > bonf, \"yes\", \"no\"),\n ## .names = \"{.col}_sig\"\n .names = \"{.col}_sig\"\n )\n ) %>%\n mutate(\n jt =\n rep(c(\"det\", \"unch\", \"imp\", \"nr_rec\", \"rec\"), 5),\n .after = var\n )\n\njt_tbl <-\n .jt_tbl %>%\n full_join(.jt_tbl_chisq)\n\n\njt_combi_plot <-\n pb_jt_fig +\n tb_jt_fig +\n bhs_jt_fig +\n dep_jt_fig +\n anx_jt_fig +\n str_jt_fig +\n si_jt_fig +\n guide_area() +\n plot_layout(\n nrow = 2,\n guides = \"collect\",\n ) &\n theme(\n plot.title = element_text(size = 10),\n legend.title = element_text(size = 10),\n legend.text = element_text(size = 9),\n axis.title.x = element_text(size = 9),\n axis.title.y = element_text(size = 9),\n axis.text = element_text(size = 7),\n ) &\n labs(\n colour = \"JT Classification\",\n shape = \"Class\"\n ) &\n scale_fill_brewer(palette = \"Dark2\") &\n scale_colour_brewer(palette = \"Dark2\")\n\nggsave(\n \"./figures/jt_combi_plot.png\",\n jt_combi_plot\n)\n\n### ============================================================================\n### ABSTRACT DF\n### ============================================================================\n\nabst_df <-\n lpa_t0 %>%\n summarise(\n N = sum(!is.na(c)),\n age_m = round(mean(age, na.rm = TRUE), 2),\n age_sd = round(sd(age, na.rm = TRUE), 2),\n prc_1 = round(sum(!is.na(c[c == 1])) / N * 100, 2),\n prc_2 = round(sum(!is.na(c[c == 2])) / N * 100, 2),\n prc_3 = round(sum(!is.na(c[c == 3])) / N * 100, 2),\n )\n\nlpa_t0 %>%\n group_by(c, protocol) %>%\n summarise(\n N = n()\n )\n\n### ============================================================================\n### END SCRIPT\n### ============================================================================\n\n## CORRELATION TABLES\n## =============================================================================\ncorr_des <-\n lpa_t0t1_outc %>%\n filter(final_session == 0) %>%\n select(\n age, sex, hxsa, hxdsh, hxsuicide, hxtrauma,\n pb, si, tb, dep, anx, str, ac,\n isas_social, isas_intra,\n ) %>%\n mutate(across(everything(), as.numeric))\n\ncts <-\n corr_des %>%\n ## mutate_all(as.numeric) %>%\n corr.test()\n\ncts_p <-\n print(corr.p(cts$r, n = 500), short = FALSE) %>%\n ## as_tibble()\n arrange(p, desc(r))\n\ncts_p %>%\n filter(\n p <= 0.05\n )\n\n## correlation plots\n## =============================================================================\ncorr_des %>%\n cor.ci() %>%\n cor.plot.upperLowerCi()\n\n## histograms\n## =============================================================================\ntest_des %>%\n ## colnames() %>%\n mutate_all(as.numeric) %>%\n select(\n age, sex, suicide_attempts,\n contains(\"hx\"),\n ends_with(\"mean\")\n ) %>%\n ## review bindwidth\n multi.hist()\n\n### ============================================================================\n### PARTCIPANTS FOR METHODS SECTION\n### ============================================================================\n## consenters\nlpa_consenters <-\n read_csv(here(.phd_sp_full_df)) %>%\n filter(session %in% 1, sp_episode %in% 1) %>%\n group_by(consent) %>%\n summarise(\n n = n()\n )\n\n## MULTIPLE EPISODES\n## -----------------------------------------------------------------------------\n## extract the uci of those used for the lpa\n.lpa_target <-\n lpa_fin %>%\n filter(!is.na(c)) %>%\n select(uci) %>%\n distinct(uci) %>%\n deframe()\n\n## get the uci and class id\n.lpa_c_uci <-\n lpa_fin %>%\n filter(!is.na(c)) %>%\n select(uci, c) %>%\n group_by(uci) %>%\n slice(1) %>%\n ungroup()\n\n## get the number of episodes by class\n.lpa_sp_eps <-\n read_csv(here(.phd_sp_full_df)) %>%\n full_join(.lpa_c_uci, by = \"uci\") %>%\n select(c, uci, sp_episode, session) %>%\n filter(!is.na(c)) %>%\n group_by(uci, sp_episode) %>%\n slice(1) %>%\n group_by(uci) %>%\n filter(sp_episode == max(sp_episode)) %>%\n group_by(sp_episode, c) %>%\n summarise(\n n = n(),\n ) %>%\n group_by(c) %>%\n mutate(prc_ep = paste0(round(n/sum(n)*100,1), \"%\"))\n\n## SEX AND AGE DIFFERENCE\n## -----------------------------------------------------------------------------\n## sex diff with age\n.lpa_sex_diff <-\n lpa_fin %>%\n filter(!is.na(c)) %>%\n group_by(ID) %>%\n slice(1) %>%\n ungroup() %>%\n select(ID, c, sex, age)\n\nAnova(lm(age ~ sex, data = .lpa_sex_diff , na.action = na.exclude), type = 3) %>%\n tidy()\n\n## parametric\nAnova(lm(age ~ sex, data = .lpa_sex_diff, na.action = na.exclude), type = 3) %>%\n tidy() %>%\n filter(term %in% \"sex\") %>%\n select(term, p.value)\n\n.lpa_sex_diff %>%\n select(age) %>%\n map( ~ tidy(pairwise.t.test(., .lpa_sex_diff$sex, p.adjust.method = \"bonf\")))\n\n## nonparametric\n.lpa_sex_diff %>%\n select(age) %>%\n map( ~ tidy(pairwise.wilcox.test(., .lpa_sex_diff$sex, p.adjust.method = \"bonf\")))\n\n.lpa_sex_diff %>%\n split(.$c) %>%\n map(\n ~ lm(age ~ sex, data = .x, na.action = na.exclude) %>%\n summary() %>%\n tidy() %>%\n filter(term %in% \"sex\") %>%\n select(term, p.value)\n ) %>%\n bind_rows(.id = \"c\")\n\nlpa_sex_age_diff <-\n .lpa_sex_diff %>%\n group_by(sex) %>%\n summarise(\n n = n(),\n across(\n age,\n list(\n mean = ~ mean(.x, na.rm = TRUE),\n median = ~ median(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n min = ~ min(.x, na.rm = TRUE),\n max = ~ max(.x, na.rm = TRUE)\n ))\n ) %>%\n filter(sex %in% c(0:2)) %>%\n mutate(\n across(c(age_mean, age_sd), ~ sprintf(\"%.2f\", .x))\n )\n\n##\n## PLANNED FINAL DIFFERENCES\n## -----------------------------------------------------------------------------\n## planned final\nlpa_planned_fin_df <-\n lpa_t0t1_outc %>%\n select(c:planned_final, si, pb:dep, -fin_ses_c, -session) %>%\n pivot_wider(\n names_from = final_session,\n values_from = c(planned_final:dep)\n ) %>%\n ## filter(final_session %in% 1) %>%\n group_by(planned_final_1, c) %>%\n summarise(\n n = n(),\n across(\n starts_with(c(\"pb\", \"ac\", \"tb\", \"bhs\", \"dep\", \"si\", \"age\", \"cov\")),\n list(\n mean = ~ mean(.x, na.rm = TRUE),\n sd = ~ sd(.x, na.rm = TRUE),\n median = ~ median(.x, na.rm = TRUE)\n )\n )) %>%\n group_by(c) %>%\n mutate(\n prc = paste0(round(n/sum(n)*100,1),\"%\"), .after = n,\n across(is.numeric, ~ round(.x, 2))\n ) %>%\n arrange(c)\n\n\nlpa_dropout <-\n lmer(\n cov_sess ~\n ## predictors scaled due to lmer warning\n planned_final +\n (1 | c),\n data = lpa_t0t1_outc[lpa_t0t1_outc$final_session %in% 1, ],\n na.action = na.exclude,\n REML = FALSE # use ML estimation\n )\n\nsummary(as(lpa_dropout, \"merModLmerTest\"))\n\nlpa_t0t1_outc <-\n lpa_t0t1_outc %>%\n group_by(ID) %>%\n mutate(\n planned_final = max(as.numeric(as.character(planned_final))),\n ) %>%\n ungroup()\n\n\nlm(\n dep ~\n planned_final,\n data = lpa_t0t1_outc[lpa_t0t1_outc$c %in% 3 & lpa_t0t1_outc$final_session %in% 0, ],\n na.action = na.exclude\n) %>%\n summary()\n\n.lpa_plan_fin <-\n lpa_fin %>%\n filter(!is.na(c), final_session %in% 1) %>%\n group_by(c, planned_final) %>%\n summarise(\n n = n()\n ) %>%\n group_by(c) %>%\n mutate(prc_fin = paste0(round(n/sum(n)*100,1), \"%\"))\n\n## is there a difference in dropout by class?\nlpa_fin %>%\n filter(!is.na(c), final_session %in% 1) %>%\n ## glm(planned_final ~ c, data = ., na.action = na.exclude) %>%\n lm(planned_final ~ c, data = ., na.action = na.exclude) %>%\n summary() %>%\n tidy()\n\n## regression with palnned final\nlpa_out_ipt_r2 <-\n lmer(\n si ~\n ## predictors scaled due to lmer warning\n pb_z *\n tb_z *\n bhs_z *\n ac_z +\n dep +\n cov_sess + age + sex + planned_final +\n ## cov_sess + age + sex +\n ## (cov_sess + planned_final | final_session/c),\n (1 | final_session/c),\n data = lpa_t0t1_outc,\n na.action = na.exclude,\n ## control = lmerControl(optimizer = \"Nelder_Mead\"),\n REML = FALSE # use ML estimation\n )\n\nall_lpa_r2 <- allFit(lpa_out_ipt_r2)\n\nsummary(all_lpa_r2)\n\nsummary(as(lpa_out_ipt_r2, \"merModLmerTest\"))\n\nlpa_out_ipt_r2_std <-\n lpa_out_ipt_r2 %>%\n stdCoef.merMod()\n\nlpa_out_ipt_r2_tidy <-\n tidy(\n ## NOTE needs to have lme4 specify the model\n as(lpa_out_ipt_r2, \"merModLmerTest\"), # from lmerTest\n ## effects = \"fixed\",\n conf.method = \"profile\",\n ## conf.method = \"boot\", # slower but more accurate\n ## nsim = 1000, # goes with \"boot\" method\n conf.int = TRUE\n ) %>%\n mutate(var = \"full\", .before = everything())\n\nlpa_out_ipt_r2_tidy %>%\n filter(p.value < 0.05) %>%\n print(n = Inf)\n\n### ============================================================================\n### WRITE IMAGE AND DATA\n### ============================================================================\n\n## TODO remove P from env\nsave.image(\"./data/proc/draft.RData\")\n\nrep_df <-\n mget(ls())\n\nsaveRDS(\n rep_df,\n \"./data/proc/lpa_df.rds\"\n)\n\n### ============================================================================\n### ENDS\n### ============================================================================\n"
}
]
}
[Trace - 03:06:26 PM] Sending request 'textDocument/codeAction - (382)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
},
"range": {
"start": {
"line": 911,
"character": 0
},
"end": {
"line": 911,
"character": 0
}
},
"context": {
"diagnostics": []
}
}
[Trace - 03:06:26 PM] Sending request 'textDocument/documentLink - (383)'.
Params: {
"textDocument": {
"uri": "file:///Users/atanas/git/PR_quant_LPA/code/2_analyses.R"
}
}
[Trace - 03:06:26 PM] Sending notification '$/cancelRequest'.
Params: {
"id": 383
}
[Trace - 03:06:26 PM] Sending notification '$/cancelRequest'.
Params: {
"id": 382
}
[Trace - 03:06:26 PM] Sending notification '$/cancelRequest'.
Params: {
"id": 380
}
You are still missing the debug log from languageserver. Make sure to capture stderr of the server in the client. Additionally, I don't see any error in the log.
Sorry about that. This may be an Emacs specific question, but how do I access the debug log from languageserver?
Regarding no error, that's interesting because languageserver has stopped working…
Also, the here is a gist with the lsp-r::stderr output, not sure if it is useful or not (or if this is what you were after).
This is what I wanted.
[2021-12-02 15:10:04.056] received: Content-Length: 116294
[2021-12-02 15:10:04.096] error handling json: Error: lexical error: invalid char in json text.
~ scale(.x),\n .n{"jsonrpc":"2.0","method":"textDo
(right here) ------^
Stack trace:
1: parse_string(txt, bigint_as_char)
2: parseJSON(txt, bigint_as_char)
3: parse_and_simplify(txt = txt, simplifyVector = simplifyVector,
simplifyDataFrame = simplifyDataFrame, simplifyMatrix = simplifyMatrix,
flatten = flatten, ...)
4: jsonlite::fromJSON(data, simplifyVector = FALSE)
For some reasons, there are some random text ~ scale(.x),\n .n prepended in the json text. Any idea where it is from? is it from your file?
Sorry for the late reply to this. The only thing that I can find in that file is something like this:
...
mutate(
across(c(pb, tb, bhs, ac),
~ scale(.x),
.names = "{.col}_z"
)
) %>%
...
I have something similar to this four times. Any thoughts on why that would be happening?
I have updated the options in by .Rprofile and things seem more stable. I'm not too sure why, but posting here just in case it helps resolves this issues:
## ## options to improve performance
options(
## =========================================
## solution 1: low CPU
## =========================================
languageserver.did_change.run_lintr = FALSE,
languageserver.did_change.parse = FALSE,
languageserver.task_manager.cpu_load = 0.5
## languageserver.task_manager.delay = 0.5
## languageserver.r.source.files = c("*.R", "R/*.R"),
## languageserver.excluded.source.files = "renv/**"
)
options(
languageserver.debug = TRUE
)