## setup chunk knitr::opts_chunk$set(echo = TRUE, fig.cap = TRUE, message=FALSE, warning=FALSE) ##read in packages library(sessioninfo) library(officedown) library(officer) library(flextable) library(sessioninfo) library(ggplot2) library(tidyverse) library(dplyr) library(gtsummary) library(gt) library(gamlss) library(haven) library(logbin) library(risks) my_theme <- list( "pkgwide-fn:pvalue_fun" = function(x) style_pvalue(x, digits = 2), "pkgwide-fn:prependpvalue_fun" = function(x) style_pvalue(x, digits = 2, prepend_p = TRUE), "tbl_summary-arg:statistic" = list( all_continuous() ~ "{mean} ({sd})", all_categorical() ~ "{n} ({p}%)"), "tbl_summary-str:header-withby" = "**{level}**\n (n={n})", "tbl_summary-arg:digits" = list(all_continuous() ~ 1), "tbl_summary-arg:missing_text" = "Missing", "tbl_summary-fn:addnl-fn-to-run" = bold_labels ) set_gtsummary_theme(my_theme) session_info(pkgs="loaded") # ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── # setting value # version R version 4.4.2 (2024-10-31 ucrt) # os Windows 11 x64 (build 26100) # system x86_64, mingw32 # ui RStudio # language (EN) # collate English_Australia.utf8 # ctype English_Australia.utf8 # tz Australia/Sydney # date 2026-01-13 # rstudio 2025.09.2+418 Cucumberleaf Sunflower (desktop) # pandoc 3.6.3 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown) # # ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── # package * version date (UTC) lib source # askpass 1.2.0 2023-09-03 [1] CRAN (R 4.4.1) # cachem 1.1.0 2024-05-16 [1] CRAN (R 4.4.1) # cli 3.6.3 2024-06-21 [1] CRAN (R 4.4.1) # colorspace 2.1-1 2024-07-26 [1] CRAN (R 4.4.1) # data.table 1.16.0 2024-08-27 [1] CRAN (R 4.4.1) # digest 0.6.37 2024-08-19 [1] CRAN (R 4.4.1) # dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.4.1) # evaluate 1.0.3 2025-01-10 [1] CRAN (R 4.4.3) # fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.4.1) # flextable * 0.9.6 2024-05-05 [1] CRAN (R 4.4.1) # fontBitstreamVera 0.1.1 2017-02-01 [1] CRAN (R 4.4.0) # fontLiberation 0.1.0 2016-10-15 [1] CRAN (R 4.4.0) # fontquiver 0.2.1 2017-02-01 [1] CRAN (R 4.4.1) # forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.4.1) # gamlss * 5.4-22 2024-03-20 [1] CRAN (R 4.4.1) # gamlss.data * 6.0-7 2025-09-04 [1] CRAN (R 4.4.3) # gamlss.dist * 6.1-1 2023-08-23 [1] CRAN (R 4.4.1) # gdtools 0.4.0 2024-08-28 [1] CRAN (R 4.4.1) # generics 0.1.3 2022-07-05 [1] CRAN (R 4.4.1) # ggplot2 * 3.5.1 2024-04-23 [1] CRAN (R 4.4.1) # glm2 1.2.1 2018-08-11 [1] CRAN (R 4.4.0) # glue 1.7.0 2024-01-09 [1] CRAN (R 4.4.1) # gt * 0.11.0 2024-07-09 [1] CRAN (R 4.4.1) # gtable 0.3.6 2024-10-25 [1] CRAN (R 4.4.3) # gtsummary * 2.0.2 2024-09-05 [1] CRAN (R 4.4.1) # haven * 2.5.4 2023-11-30 [1] CRAN (R 4.4.1) # hms 1.1.3 2023-03-21 [1] CRAN (R 4.4.1) # htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.4.1) # iterators 1.0.14 2022-02-05 [1] CRAN (R 4.4.1) # itertools2 0.1.1 2014-08-08 [1] CRAN (R 4.4.2) # knitr 1.48 2024-07-07 [1] CRAN (R 4.4.1) # lattice 0.22-6 2024-03-20 [2] CRAN (R 4.4.2) # lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.4.1) # logbin * 2.0.5 2021-08-09 [1] CRAN (R 4.4.2) # lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.4.1) # magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.4.1) # MASS 7.3-61 2024-06-13 [2] CRAN (R 4.4.2) # Matrix 1.7-1 2024-10-18 [2] CRAN (R 4.4.2) # memoise 2.0.1 2021-11-26 [1] CRAN (R 4.4.1) # munsell 0.5.1 2024-04-01 [1] CRAN (R 4.4.1) # nlme * 3.1-166 2024-08-14 [2] CRAN (R 4.4.2) # officedown * 0.3.3 2024-07-16 [1] CRAN (R 4.4.1) # officer * 0.6.6 2024-05-05 [1] CRAN (R 4.4.1) # openssl 2.2.1 2024-08-16 [1] CRAN (R 4.4.1) # pillar 1.10.2 2025-04-05 [1] CRAN (R 4.4.3) # pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.4.1) # pkgload 1.4.0 2024-06-28 [1] CRAN (R 4.4.3) # purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.4.1) # R6 2.6.1 2025-02-15 [1] CRAN (R 4.4.3) # ragg 1.3.2 2024-05-15 [1] CRAN (R 4.4.1) # Rcpp 1.0.13 2024-07-17 [1] CRAN (R 4.4.1) # readr * 2.1.5 2024-01-10 [1] CRAN (R 4.4.1) # rlang 1.1.4 2024-06-04 [1] CRAN (R 4.4.1) # rmarkdown 2.30 2025-09-28 [1] CRAN (R 4.4.3) # rstudioapi 0.16.0 2024-03-24 [1] CRAN (R 4.4.1) # rvg 0.3.4 2024-08-27 [1] CRAN (R 4.4.1) # scales 1.3.0 2023-11-28 [1] CRAN (R 4.4.1) # sessioninfo * 1.2.2 2021-12-06 [1] CRAN (R 4.4.1) # stringi 1.8.4 2024-05-06 [1] CRAN (R 4.4.0) # stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.4.1) # survival 3.7-0 2024-06-05 [2] CRAN (R 4.4.2) # systemfonts 1.1.0 2024-05-15 [1] CRAN (R 4.4.1) # textshaping 0.4.0 2024-05-24 [1] CRAN (R 4.4.1) # tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.4.1) # tidyr * 1.3.1 2024-01-24 [1] CRAN (R 4.4.1) # tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.4.1) # tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.4.1) # timechange 0.3.0 2024-01-18 [1] CRAN (R 4.4.1) # tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.4.1) # uuid 1.2-1 2024-07-29 [1] CRAN (R 4.4.1) # vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.4.1) # withr 3.0.2 2024-10-28 [1] CRAN (R 4.4.3) # xfun 0.50 2025-01-07 [1] CRAN (R 4.4.2) # xml2 1.3.6 2023-12-04 [1] CRAN (R 4.4.1) # yaml 2.3.10 2024-07-26 [1] CRAN (R 4.4.1) # zip 2.3.1 2024-01-27 [1] CRAN (R 4.4.1) # # [1] C:/Users/ghel3647/AppData/Local/R/win-library/4.4 # [2] C:/Program Files/R/R-4.4.2/library ## ----------------------------------------------------------------------------------------------------------------------------------------------- read_spss("3. FASTStudy-Falls SPSS -Table 2 data.sav") |> filter(!is.na(ID)) |> mutate(FallerYesNo = haven::as_factor(FallerYesNo), RandomisationGroup = factor(RandomisationGroup, labels=c("Control", "Intervention")) ) -> FASTdata ## ----------------------------------------------------------------------------------------------------------------------------------------------- FASTdata |> dplyr::select(RandomisationGroup, Totaldaysinstudy, FallerYesNo, TotalFalls12month) |> tbl_summary(by=RandomisationGroup, statistic = all_continuous() ~ "{mean} ({sd})") |> as_flex_table() ## 12-month falls -------------------------------------------------------------------------------------- FASTdata |> dplyr::select(RandomisationGroup, TotalFalls12month, Totaldaysinstudy) |> na.omit() -> FASTsubset m1 <- gamlss(TotalFalls12month ~ RandomisationGroup + offset(log(Totaldaysinstudy)), data=FASTsubset, family=NBI, trace = FALSE) a = summary(m1) m2 <- gamlss(TotalFalls12month ~ -1 + RandomisationGroup + offset(log(Totaldaysinstudy)), data=FASTsubset, family=NBI, trace = FALSE) newdata = FASTsubset[c(3,1),] p2=predict(m2, newdata=newdata, data=FASTsubset, type="response") se = sqrt(diag(vcov(m2)))[1:2] lowerControl = round(exp(log(p2)[1]-2*se[1]),2) upperControl = round(exp(log(p2)[1]+2*se[1]),2) lowerIntervention = round(exp(log(p2)[2]-2*se[2]),2) upperIntervention = round(exp(log(p2)[2]+2*se[2]),2) ## 12-month rate of falls (mean (95% confidence interval))-------------------------------------------------- # Intervention round(p2[2],2) # [1] 1.85 c(lowerIntervention,upperIntervention) # RandomisationGroupIntervention RandomisationGroupIntervention # 1.45 2.36 # Control round(p2[1],2) # [1] 2.76 c(lowerControl,upperControl) # RandomisationGroupControl RandomisationGroupControl # 2.19 3.48 # p-value for difference in rate of falls between Intervention and Control: round(a[2,4],3) # [1] 0.017 # Note that the above mean rates of falls are slightly higher than what is in the table above, # because they are adjusting for the fact that not all subjects were in the study for 365 days. ## Rate ratio ------------------------------------------------------------------------------------------------------------------------- LT <- round(exp(a[2,1]),2) LT.lower <- round(exp(a[2,1] -2*a[2,2]),2) LT.upper <- round(exp(a[2,1] +2*a[2,2]),2) # Rate ratio Intervention/Control LT c(LT.lower,LT.upper) # p-value round(a[2,4],3) # > LT # [1] 0.67 # > c(LT.lower,LT.upper) # [1] 0.48 0.94 # > # p-value # > round(a[2,4],3) # [1] 0.017 ## Proportion of participants experiencing a fall ------------------------------------------------------------------------------------- FASTdata |> mutate(faller = as.numeric(TotalFalls12month>0)) |> dplyr::select(RandomisationGroup, faller) |> na.omit() -> FASTsubset # log-binomial regression m7 <- logbin(faller ~ RandomisationGroup, data=FASTsubset) a = summary(m7) b = confint(m7) # Relative risk for fall round(exp(a$coefficients[2,1]),2) # 95% confidence interval c(round(exp(b[2,]),2)) # p-value a$coefficients[2,4] # > round(exp(a$coefficients[2,1]),2) # [1] 0.94 # > # 95% confidence interval # > c(round(exp(b[2,]),2)) # 2.5 % 97.5 % # 0.79 1.12 # > # p-value # > a$coefficients[2,4] # [1] 0.5178049 ## Risk difference ------------------------------------------------------------------------------------------------------------------------ fit <- riskdiff(formula = faller ~ RandomisationGroup, data = FASTsubset) a = summary(fit) # RD (Risk reduction) in Table 2 is negative of risk difference round(-fit$coefficients[2],2) round(c(-confint(fit)[2,2], -confint(fit)[2,1]),2) # p-value a$coefficients[2,4] # > round(-fit$coefficients[2],2) # RandomisationGroupIntervention # 0.03 # > round(c(-confint(fit)[2,2], -confint(fit)[2,1]),2) # [1] -0.07 0.13 # > # p-value # > a$coefficients[2,4] # [1] 0.5173388