Physical Activity and Low Back Pain in Switzerland

A Cross-Sectional Study based on the 2022 Swiss Health Survey

Authors
Affiliations

Alexander Schurz

University of Bern (Switzerland)

Bern University of Applied Sciences (Switzerland)

Marcel Zwahlen

Jan Taeymans

Nathanael Lutz

Published

October 22, 2025

1 Load packages

Code 
library(pacman)
pacman::p_load(
  tidyverse,
  ggplot2,
  ggdag,
  dagitty,
  dplyr, 
  purrr,
  epitools,
  parameters,
  kableExtra,
  DataExplorer,
  naniar,
  nnet,
  naniar,
  gtsummary,
  srvyr, # develop survey
  survey, # develop survey
  mice, 
  emmeans,
  VIM,
  table1,
  latex2exp,
  svyVGAM, # model analysis
  parameters, # model description and presentation
  flextable) #to extract model outputs

2 Data import

Code 
# Data from telephone interviews
df.tel.raw <- read.delim("../data/tel22_CH.txt", header=TRUE)

# Data from written surveys
df.sfb.raw <- read.delim("../data/sfb22_CH.txt", header=TRUE)

# Calculated indices
df.ind.raw <- read.delim("../data/indic22_CH.txt", header=TRUE)

3 Research questions

  1. What is the association between leisure-time PA and LBP in adults living in Switzerland in 2022?
  2. What other risk factors (e.g. BMI, age, occupational factors, psychosocial factors, and medical history) are associated with LBP in adults living in Switzerland in 2022?

4 Data preparation per survey type

4.1 Data of the telephone-based survey (df.tel.raw)

Code 
# data preparation df.tel.ad
df.tel.ad <- df.tel.raw %>% 
  mutate(
    # ID
    id = IDNO,
    # Sociodemographic: ALTER (int), SEX (int),
    # Kanton (int), TGEZU01 (num, height), TGEZU02 (num, weight), TSOUN24 (int, subj QoL)
    age = as.numeric(ALTER),
    sex = factor(SEX, levels = c(1,2), labels = c("Male", "Female")),
    strata = factor(KANTON, levels = c(1:26), labels = 
                      c("Zurich", "Bern", "Lucerne", "Uri",
                        "Schwyz", "Obwalden", "Nidwalden", "Glarus",
                        "Zug", "Fribourg", "Solothurn", "Basel-City",
                        "Basel-Country", "Schaffhausen", "Appenzell Ausserrhoden",
                        "Appenzell Innerrhoden", "St. Gallen", "Graubünden",
                        "Aargau", "Thurgau", "Tessin", "Vaud",
                        "Wallis", "Neuchâtel", "Geneva", "Jura")),
    height = as.numeric(TGEZU01),
    weight = as.numeric(TGEZU02),
    qol = factor(TSOUN24, levels = c(1:5), labels = c("Very good", "Good", 
                                                      "Neither good nor bad", "Bad", "Very bad")),
    # Age categories (ALTER7)
    age_cat = factor(ALTER7, levels = c(1:7), labels = c("18-24 years", "25-34 years", 
                                                         "35-44 years", "45-54 years", 
                                                         "55-64 years", "65-74 years", 
                                                         "75+ years")),
    
    # Dependent variable: TKRSY01 (int, back pain/low back pain, values 1 to 3, not -2)
    lbp = factor(TKRSY01, levels = c(1:3), labels = c("No", "A little", "Strong")),
    
    # Sedentary behaviour:
    sed_beh = as.numeric(TKOBW17), # Sitting: hours/day 
    # Weight of telephone survey: WGHT (int)
    weighting.tel = as.numeric(WGHT),
    
    # Potential confounder
    # preg: TSCHW01 (int, preg at the moment)
    preg = factor(TSCHW01, levels = c(1,2), labels = c("Yes", "No")),
    
    # Connection between back pain and work: TKRSY35 (int, back/low back)
    lbp_w = factor(TKRSY35, levels = c(1:5), 
                              labels = c("Yes, certainly", 
                                         "Yes, rather", 
                                         "Partly", 
                                         "No, rather not", 
                                         "No, definitely not")),
    
    # Fever last 4 weeks
    fever = factor(TKRSY09, levels = c(1:3), labels = c("No", "A little", "Strong")),
    
    # Shoulder, Neck, Arm pain last 4 weeks
    shoul_neck_pain = factor(TKRSY09, levels = c(1:3), labels = c("No", "A little", "Strong")),
    
    # Medical diagnoses e.g. arthroses, asthma, incontinence
    diag_n = rowSums(across(c(
    TKRAN10l, TKRAN10m, TKRAN10d, TKRAN10a, TKRAN10e, 
    TKRAN10b, TKRAN10i, TKRAN10o, TKRAN10j, TKRAN11l, 
    TKRAN11m, TKRAN11d, TKRAN11a, TKRAN11e, TKRAN11b, 
    TKRAN11i, TKRAN11o, TKRAN12a, TKRAN12b, TKRAN12c), ~ . == "1"), na.rm = TRUE),
    
    # Received care via SPITEX within the last 12 months
    spitex = factor(TINAN24, levels = c(1:2), labels = c("Yes", "No")),
    
    # Hospital stays in last 12 months in days
    hospital = as.numeric(TINAN68))

4.2 Data of the written survey (df.sfb.raw)

Code 
df.sfb.ad <- df.sfb.raw %>%
  rename(id = IDNO) %>%
  mutate(income = factor(case_when(
    SSODE02a %in% 1:2 ~ 1,             # < 4000 Franken
    SSODE02a == 3     ~ 2,             # 4001 bis 6000 Franken
    SSODE02a == 4     ~ 3,             # 6001 bis 8000 Franken
    SSODE02a == 5     ~ 5,             # 8001 bis 10 000 Franken
    SSODE02a %in% 6:11 ~ 6             # > 10 000 Franken
  ), levels = c(1, 2, 3, 5, 6),
     labels = c(
       "< 4000 CHF",
       "4001 bis 6000 CHF",
       "6001 bis 8000 CHF",
       "8001 bis 10 000 CHF",
       "> 10 000 CHF")),
     weighting.sfb = as.numeric(WECRIT))

4.3 Data of the indices (df.ind.raw)

Code 
# data preparation for df.ind.ad
df.ind.ad <- df.ind.raw %>% 
  mutate(
    # ID
    id = IDNO,
    # Sociodemographic: BMI (num), BMI4 (int), civil status (ETATCIV, int), region in CH (SPRACHE, int), rural/urban area (STALA, int)
    bmi = round(as.numeric(BMI),2),
    bmi_cat = factor(BMI4, levels = c(1:4), labels = c("Underweight", "Normal weight", 
                                                       "Overweight", "Obesity")),
    civil_status = factor(ETATCIV, levels = c(1:4), labels = c("Single", "Married", 
                                                               "Divorced", "Widowed")),
    # nat: NATION2 (int, Swiss or foreigner), NATION5 (int, region of origin)
    nat = factor(NATION2, levels = c(1:2), labels = c("Swiss", "Foreigner")),
    nat_cat = factor(NATION5, levels = c(1:5), labels = c("Switzerland", 
                                                                  "Northern & Western Europe", 
                                                                  "South-West Europe", 
                                                                  "Eastern & South-Eastern Europe",
                                                                  "Outside Europe")),
    lang_r = factor(SPRACHE, levels = c(1:3), labels = c("German-speaking Switzerland",
                                                                  "French-speaking Switzerland", 
                                                                  "Italian-speaking Switzerland")),
    region = factor(STALA, levels = c(1:3), 
                    labels = c("Urban", "Intermediate", "Rural")),
    # Independent variable: ACTPHY5
    pa_cat = factor(ACTPHY5, levels = c(1:5), 
                    labels = c("Inactive", "Partially active", 
                               "Irregularly active", 
                               "Regularly active", 
                               "Trained")),
    # Working status
    working_status = factor(ERWERB, levels = c(1:3), 
                            labels = c("Inactive", "Unemployed", "Employed")),
    # Level of employment
    work_per_cat = factor(TXTRAV4, levels = c(1:4), 
                              labels = c("Part-time (<50%)", 
                                         "Part-time (50-69%)", 
                                                                  
                                         "Part-time (70-89%)", 
                                         "Full-time (90-100%)")),
    # Profession
    prof = factor(CH_ISCO_19_MAJOR, levels = c(0:9), 
                        labels = c("members of the reular armed forces", 
                                   "managerial occupations", 
                                                                      
                                   "academic occupation", 
                                   "technicians and related non-technical occupations",
                                   "office workers and related professions", 
                                   "service professions and sales people",
                                   "specialists in agriculture, forestry and fisheries", 
                                   "craft and related trades",
                                   "plant and machine operators and assembly occupations", 
                                   "support staff")),
    # Lifestyle factors
    ## Diet
    diet = factor(FIVEDAY, levels = c(1:4), 
                  labels = c("<5 days per week", 
                             "0-2 portions", 
                             "3-4 portions", 
                             "≥5 portions")),
    ## Sleeping disorder
    sleep = factor(SOMMEIL, levels = c(1:3), 
                   labels = c("none or few", 
                              "moderate", 
                              "pathological")),
    
    # Social relationships - Social support strong as reference value
    soc = factor(OSS3, levels = c(1:3), 
                         labels = c("low", "medium", "strong")), 

    ## Risky substances - Tobacco - TABAC3 (int), DAYSMOKE (int), Kind of tobacco (TABACTYP), number of cigarettes per day (NBCIG5, int)
    tobacco = factor(DAYSMOKE, levels = c(1:3), 
                     labels = c("non-smoker", 
                                "occasional smoker", 
                                "daily smoker")),
    
    ## Risky substances - Alcohol - BIERFREQ (int), POPSFREQ (int), SPIRFREQ (int), WEINFREQ (int), ALCCHRON5 (chronic alcohol consumption, int)
    alcohol = factor(ALCCHRON5, levels = c(1:5), 
                     labels = c("abstinent", 
                                "low risk", 
                                "minimal risk", 
                                "medium risk", 
                                "increased risk")),

    ## Risky substances - Drugs) - total drug consumption (DROGCONS, int)
    drugs = factor(DROGCONS, levels = c(1:5), 
                   labels = c("never used drugs", 
                              "never used cannabis", 
                              ">12 months ago", 
                              "in the last 12 months", 
                              "in the last 30 days")),
    
    # Model by Hartvigsen
    # Depressive symptoms, DEPPHQ5 (int)
    dep = factor(DEPPHQ5, 
                        levels = c(1, 2, 3, 4, 5), 
                        labels = c("No/ minimal symptoms of depression", 
                                   "Mild depression symptoms", 
                                   "Moderate depression symptoms", 
                                   "Rather severe depression symptoms", 
                                   "Severe depression symptoms")),
    
    # Anxiety disorder: ANXGAD4 (int)
    anxiety = factor(ANXGAD4, 
                     levels = c(1, 2, 3, 4), 
                     labels = c("None", "Light", "Moderate", "Severe")),
    
    # Educational level
    education = factor(AUSBILD3, 
                       levels = c(1, 2, 3), 
                       labels = c("Compulsory school", 
                                  "Upper secondary education", 
                                  "Tertiary level")),
    
    # Arterial Hypertension 
    aht = factor(HYPERTENS, 
                       levels = c(2, 1), 
                       labels = c("No", "Yes")),
    
    # Hypercholesterolaemia
    hyperchol = factor(HYPERCHOL, 
                   levels = c(2, 1), 
                   labels = c("No", "Yes")),
    
    # Diabetes
    diabetes = factor(DIABETE, 
                  levels = c(2, 1), 
                  labels = c("No", "Yes")),
    
    # Limitations in functional ADL (sight, hearing, speech, ability to walk)
    func_adl_lim = factor(LIMFONC, levels = c(1:4), 
                          labels = c("no difficulties", 
                                     "slight difficulties", 
                                     "severe difficulties", 
                                     "inability")),
    
    # Psychological stress
    psy_stress = factor(DETPSY3, levels = c(1:3), labels = c("low", "medium", "high")))

5 Variables used for analysis

5.1 Sociodemographic/socioeconomic variables, strata (cantons), weight variable

Code 
# Sociodemographic variables
var.socdem <- c("id", "age_cat", "bmi_cat", "civil_status", "lang_r", "nat", "sex")

# Socioeconomic variables
var.soceco <- c("education")

# Income was not included as it had >15k missing values.
# Stratae variable for cantons
var.strata <- c("strata")

## Weight variable for written survey to be used as income was also considered from written survey
var.weight <- c("weighting.tel")

5.2 Dependent (LBP) & independent variable (PA)

Code 
var.out <- c("lbp", "pa_cat")

5.3 Lifestyle factors and factors by Hartvigsen et al. 2018

Code 
# Lifestyle factors
var.life <- c("diet", "sleep", "sed_beh", "soc",
              "tobacco", "alcohol", "drugs")

# Model by Hartvigsen et al. 2018
var.model <- c(
    # sex + BMI + education already included as confounder
    # Depressive symptoms DEPPHQ5
    "dep",
    # Anxiety ANXGAD4
    "anxiety",
    # Arterial Hypertension, HYPERTENS
    "aht",
    # Hypercholesterolaemia, HYPERCHOL
    "hyperchol",
    # Diabetes, DIABETE
    "diabetes")
Code 
# without confounder
var.all <- c(var.socdem,
             var.soceco,
             var.strata, 
             var.out, 
             var.life, 
             var.model, 
             var.weight)
Code 
# Categorial variables
var.cat <- c("id", "sex", "age_cat", 
             "bmi_cat", "civil_status", "nat", "lang_r", 
             "strata", "weighting.tel", 
             "lbp", "pa_cat", "diet", "sleep", "soc", 
             "tobacco", "alcohol", "drugs", "dep", "anxiety", 
             "education", "aht", "hyperchol", "diabetes")

6 Merge datasets

6.1 Adapt merged datasets

Code 
# 1. Dataset written survey df.sfb.ad, n = 19,137 
df.m.raw <-  df.tel.ad %>% 
  left_join(df.ind.ad, by = "id") %>% 
  select(var.all, age)

write.csv(df.m.raw, "../datasets/df.m.raw.csv", row.names = TRUE)

# 2. Excluding those age < 18 years
df.m.ad <- df.m.raw %>% 
  filter(age>= 18)

write.csv(df.m.ad, "../datasets/df.m.ad.csv", row.names = TRUE)

# 3. Excluding proxy answers (= negative responses)
df.m.ad.2 <- df.m.ad %>% 
  filter(sed_beh > 0) %>% 
  select(-age)
summary(df.m.ad$sed_beh)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -6.000   3.000   5.000   4.976   7.000  16.000
Code 
# 4. complete cases only
df.complete <- df.m.ad.2 %>%
  drop_na()
visdat::vis_miss(df.m.ad.2)

Code 
naniar::miss_var_summary(df.m.ad.2)
# A tibble: 24 × 3
   variable  n_miss pct_miss
   <chr>      <int>    <num>
 1 dep         2522   12.5  
 2 anxiety     2305   11.5  
 3 sleep       2052   10.2  
 4 hyperchol   2017   10.0  
 5 drugs       1426    7.09 
 6 diabetes    1386    6.89 
 7 soc          748    3.72 
 8 aht          487    2.42 
 9 pa_cat       370    1.84 
10 bmi_cat      174    0.865
# ℹ 14 more rows

6.2 Datasets for analysis

Code 
# Df with LBP as binary variable
# Merged dataset with LBP as binary variable, without LBP work related, with confounder variables
df.m.LBP.bin <- df.complete %>% 
  mutate(
    lbp_bin = case_when(
      lbp == "No" ~ "No",
      lbp %in% c("A little", "Strong") ~ "Yes",
      TRUE ~ NA_character_
    ),
    lbp_bin = factor(lbp_bin, levels = c("Yes", "No")))

# Separate df per sex
df.male <- df.m.LBP.bin %>% 
  filter(sex == "Male")
save(df.male, file = "../data/df.male.rda")

df.female <- df.m.LBP.bin %>% 
  filter(sex == "Female")
save(df.female, file = "../data/df.female.rda")

7 Summary (unweighted) sample

Code 
tbl.demo.lbp <- df.m.LBP.bin %>%
  select(
    sex, age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_summary(
    by = lbp,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category",
      sex = "Sex"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.demo.lbp.male <- df.male %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_summary(
    by = lbp,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()


tbl.demo.lbp.fem <- df.female %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_summary(
    by = lbp,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.demo.lbp.sex <- gtsummary::tbl_merge(
  tbls = list(tbl.demo.lbp.male, tbl.demo.lbp.fem, tbl.demo.lbp),
  tab_spanner = c("**Male**", "**Female**", "**Overall**")
)

tbl.demo.lbp.sex %>% 
  as_flex_table() %>% 
  autofit() %>%
  bold(part = "header") %>%
  align(j = 1, align = "left", part = "all")
tbl.demo.lbp.sex %>%
  as_flex_table() %>% 
  set_caption("Table. Descriptive (unweighted) statistics for low back pain, separated by sex.") %>%
  flextable::set_table_properties(width = 1, layout = "autofit") %>% 
  flextable::fontsize(size = 8, part = "all") %>% 
  flextable::font(fontname = "Times New Roman", part = "all") %>% 
  flextable::save_as_docx(path = "../tables/unweighted descriptives/tbl.demo.lbp.sex.docx")
Table 1: Summary (unweighted) sample by lbp - separated by sex

Male

Female

Overall

Characteristic

No
N = 3,7671

A little
N = 2,1091

Strong
N = 3351

No
N = 3,6301

A little
N = 2,9551

Strong
N = 6221

No
N = 7,3971

A little
N = 5,0641

Strong
N = 9571

Age category

18-24 years

328 (67%)

147 (30%)

11 (2.3%)

256 (48%)

235 (44%)

41 (7.7%)

584 (57%)

382 (38%)

52 (5.1%)

25-34 years

408 (59%)

249 (36%)

30 (4.4%)

400 (48%)

361 (43%)

71 (8.5%)

808 (53%)

610 (40%)

101 (6.6%)

35-44 years

564 (61%)

317 (34%)

43 (4.7%)

614 (52%)

471 (40%)

105 (8.8%)

1,178 (56%)

788 (37%)

148 (7.0%)

45-54 years

715 (60%)

403 (34%)

64 (5.4%)

727 (51%)

597 (42%)

113 (7.9%)

1,442 (55%)

1,000 (38%)

177 (6.8%)

55-64 years

846 (61%)

450 (32%)

92 (6.6%)

762 (50%)

613 (40%)

149 (9.8%)

1,608 (55%)

1,063 (37%)

241 (8.3%)

65-74 years

655 (59%)

387 (35%)

69 (6.2%)

646 (53%)

480 (39%)

96 (7.9%)

1,301 (56%)

867 (37%)

165 (7.1%)

75+ years

251 (58%)

156 (36%)

26 (6.0%)

225 (48%)

198 (42%)

47 (10%)

476 (53%)

354 (39%)

73 (8.1%)

BMI category

Underweight

30 (57%)

20 (38%)

3 (5.7%)

212 (54%)

150 (38%)

29 (7.4%)

242 (55%)

170 (38%)

32 (7.2%)

Normal weight

1,663 (60%)

953 (35%)

134 (4.9%)

2,330 (53%)

1,757 (40%)

328 (7.4%)

3,993 (56%)

2,710 (38%)

462 (6.4%)

Overweight

1,604 (63%)

810 (32%)

126 (5.0%)

758 (47%)

714 (44%)

156 (9.6%)

2,362 (57%)

1,524 (37%)

282 (6.8%)

Obesity

470 (54%)

326 (38%)

72 (8.3%)

330 (43%)

334 (43%)

109 (14%)

800 (49%)

660 (40%)

181 (11%)

Civil status

Single

1,131 (61%)

639 (35%)

75 (4.1%)

974 (50%)

813 (42%)

147 (7.6%)

2,105 (56%)

1,452 (38%)

222 (5.9%)

Married

2,253 (60%)

1,250 (33%)

229 (6.1%)

2,094 (51%)

1,670 (41%)

352 (8.6%)

4,347 (55%)

2,920 (37%)

581 (7.4%)

Divorced

326 (61%)

183 (34%)

27 (5.0%)

403 (48%)

354 (42%)

87 (10%)

729 (53%)

537 (39%)

114 (8.3%)

Widowed

57 (58%)

37 (38%)

4 (4.1%)

159 (51%)

118 (38%)

36 (12%)

216 (53%)

155 (38%)

40 (9.7%)

Education

Compulsory school

256 (59%)

157 (36%)

22 (5.1%)

331 (46%)

302 (42%)

93 (13%)

587 (51%)

459 (40%)

115 (9.9%)

Upper secondary education

1,558 (59%)

931 (35%)

168 (6.3%)

1,763 (49%)

1,504 (42%)

344 (9.5%)

3,321 (53%)

2,435 (39%)

512 (8.2%)

Tertiary level

1,953 (63%)

1,021 (33%)

145 (4.6%)

1,536 (54%)

1,149 (40%)

185 (6.4%)

3,489 (58%)

2,170 (36%)

330 (5.5%)

Nationality category

Swiss

3,115 (61%)

1,726 (34%)

262 (5.1%)

3,038 (51%)

2,454 (41%)

491 (8.2%)

6,153 (56%)

4,180 (38%)

753 (6.8%)

Foreigner

652 (59%)

383 (35%)

73 (6.6%)

592 (48%)

501 (41%)

131 (11%)

1,244 (53%)

884 (38%)

204 (8.7%)

Language region

German-speaking Switzerland

2,751 (62%)

1,462 (33%)

233 (5.2%)

2,582 (51%)

2,054 (41%)

427 (8.4%)

5,333 (56%)

3,516 (37%)

660 (6.9%)

French-speaking Switzerland

748 (56%)

514 (38%)

82 (6.1%)

768 (48%)

689 (43%)

153 (9.5%)

1,516 (51%)

1,203 (41%)

235 (8.0%)

Italian-speaking Switzerland

268 (64%)

133 (32%)

20 (4.8%)

280 (52%)

212 (40%)

42 (7.9%)

548 (57%)

345 (36%)

62 (6.5%)

Physical activity category

Inactive

175 (56%)

97 (31%)

42 (13%)

199 (45%)

167 (38%)

78 (18%)

374 (49%)

264 (35%)

120 (16%)

Partially active

505 (58%)

317 (36%)

53 (6.1%)

551 (46%)

533 (44%)

120 (10.0%)

1,056 (51%)

850 (41%)

173 (8.3%)

Irregularly active

1,082 (59%)

676 (37%)

89 (4.8%)

1,035 (48%)

927 (43%)

184 (8.6%)

2,117 (53%)

1,603 (40%)

273 (6.8%)

Regularly active

591 (58%)

351 (35%)

70 (6.9%)

721 (51%)

572 (40%)

125 (8.8%)

1,312 (54%)

923 (38%)

195 (8.0%)

Trained

1,414 (65%)

668 (31%)

81 (3.7%)

1,124 (56%)

756 (38%)

115 (5.8%)

2,538 (61%)

1,424 (34%)

196 (4.7%)

Sex

Male

3,767 (61%)

2,109 (34%)

335 (5.4%)

Female

3,630 (50%)

2,955 (41%)

622 (8.6%)

1n (%)

Code 
tbl.demo.pa <- df.m.LBP.bin %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_summary(
    by = pa_cat,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.demo.pa.male <- df.male %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_summary(
    by = pa_cat,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()


tbl.demo.pa.fem <- df.female %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_summary(
    by = pa_cat,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.demo.pa.sex <- gtsummary::tbl_merge(
  tbls = list(tbl.demo.pa.male, tbl.demo.pa.fem, tbl.demo.pa),
  tab_spanner = c("**Male**", "**Female**", "**Overall**")
)

tbl.demo.pa.sex %>% 
  as_flex_table() %>% 
  autofit() %>%
  bold(part = "header") %>%
  align(j = 1, align = "left", part = "all")
tbl.demo.pa.sex %>%
  as_flex_table() %>% 
  set_caption("Table. Descriptive (unweighted) statistics for physical activity, separated by sex.") %>%
  flextable::set_table_properties(width = 1, layout = "autofit") %>% 
  flextable::fontsize(size = 8, part = "all") %>% 
  flextable::font(fontname = "Times New Roman", part = "all") %>% 
  flextable::save_as_docx(path = "../tables/unweighted descriptives/tbl.demo.pa.sex.docx")
Table 2: Summary (unweighted) sample by pa - separated by sex

Male

Female

Overall

Characteristic

Inactive
N = 3141

Partially active
N = 8751

Irregularly active
N = 1,8471

Regularly active
N = 1,0121

Trained
N = 2,1631

Inactive
N = 4441

Partially active
N = 1,2041

Irregularly active
N = 2,1461

Regularly active
N = 1,4181

Trained
N = 1,9951

Inactive
N = 7581

Partially active
N = 2,0791

Irregularly active
N = 3,9931

Regularly active
N = 2,4301

Trained
N = 4,1581

Age category

18-24 years

11 (2.3%)

49 (10%)

94 (19%)

64 (13%)

268 (55%)

25 (4.7%)

99 (19%)

144 (27%)

60 (11%)

204 (38%)

36 (3.5%)

148 (15%)

238 (23%)

124 (12%)

472 (46%)

25-34 years

21 (3.1%)

103 (15%)

188 (27%)

82 (12%)

293 (43%)

45 (5.4%)

154 (19%)

239 (29%)

152 (18%)

242 (29%)

66 (4.3%)

257 (17%)

427 (28%)

234 (15%)

535 (35%)

35-44 years

46 (5.0%)

169 (18%)

289 (31%)

131 (14%)

289 (31%)

63 (5.3%)

230 (19%)

356 (30%)

228 (19%)

313 (26%)

109 (5.2%)

399 (19%)

645 (31%)

359 (17%)

602 (28%)

45-54 years

79 (6.7%)

201 (17%)

348 (29%)

176 (15%)

378 (32%)

76 (5.3%)

265 (18%)

435 (30%)

267 (19%)

394 (27%)

155 (5.9%)

466 (18%)

783 (30%)

443 (17%)

772 (29%)

55-64 years

68 (4.9%)

212 (15%)

445 (32%)

233 (17%)

430 (31%)

103 (6.8%)

237 (16%)

483 (32%)

296 (19%)

405 (27%)

171 (5.9%)

449 (15%)

928 (32%)

529 (18%)

835 (29%)

65-74 years

60 (5.4%)

104 (9.4%)

345 (31%)

223 (20%)

379 (34%)

79 (6.5%)

149 (12%)

361 (30%)

278 (23%)

355 (29%)

139 (6.0%)

253 (11%)

706 (30%)

501 (21%)

734 (31%)

75+ years

29 (6.7%)

37 (8.5%)

138 (32%)

103 (24%)

126 (29%)

53 (11%)

70 (15%)

128 (27%)

137 (29%)

82 (17%)

82 (9.1%)

107 (12%)

266 (29%)

240 (27%)

208 (23%)

BMI category

Underweight

2 (3.8%)

10 (19%)

14 (26%)

9 (17%)

18 (34%)

24 (6.1%)

60 (15%)

98 (25%)

95 (24%)

114 (29%)

26 (5.9%)

70 (16%)

112 (25%)

104 (23%)

132 (30%)

Normal weight

112 (4.1%)

354 (13%)

726 (26%)

425 (15%)

1,133 (41%)

206 (4.7%)

684 (15%)

1,294 (29%)

880 (20%)

1,351 (31%)

318 (4.4%)

1,038 (14%)

2,020 (28%)

1,305 (18%)

2,484 (35%)

Overweight

124 (4.9%)

359 (14%)

824 (32%)

407 (16%)

826 (33%)

123 (7.6%)

301 (18%)

523 (32%)

287 (18%)

394 (24%)

247 (5.9%)

660 (16%)

1,347 (32%)

694 (17%)

1,220 (29%)

Obesity

76 (8.8%)

152 (18%)

283 (33%)

171 (20%)

186 (21%)

91 (12%)

159 (21%)

231 (30%)

156 (20%)

136 (18%)

167 (10%)

311 (19%)

514 (31%)

327 (20%)

322 (20%)

Civil status

Single

96 (5.2%)

284 (15%)

467 (25%)

247 (13%)

751 (41%)

108 (5.6%)

354 (18%)

539 (28%)

308 (16%)

625 (32%)

204 (5.4%)

638 (17%)

1,006 (27%)

555 (15%)

1,376 (36%)

Married

160 (4.3%)

516 (14%)

1,199 (32%)

656 (18%)

1,201 (32%)

231 (5.6%)

673 (16%)

1,289 (31%)

858 (21%)

1,065 (26%)

391 (5.0%)

1,189 (15%)

2,488 (32%)

1,514 (19%)

2,266 (29%)

Divorced

50 (9.3%)

67 (13%)

151 (28%)

83 (15%)

185 (35%)

81 (9.6%)

137 (16%)

229 (27%)

176 (21%)

221 (26%)

131 (9.5%)

204 (15%)

380 (28%)

259 (19%)

406 (29%)

Widowed

8 (8.2%)

8 (8.2%)

30 (31%)

26 (27%)

26 (27%)

24 (7.7%)

40 (13%)

89 (28%)

76 (24%)

84 (27%)

32 (7.8%)

48 (12%)

119 (29%)

102 (25%)

110 (27%)

Education

Compulsory school

49 (11%)

56 (13%)

101 (23%)

90 (21%)

139 (32%)

102 (14%)

130 (18%)

183 (25%)

160 (22%)

151 (21%)

151 (13%)

186 (16%)

284 (24%)

250 (22%)

290 (25%)

Upper secondary education

164 (6.2%)

391 (15%)

796 (30%)

461 (17%)

845 (32%)

219 (6.1%)

581 (16%)

1,124 (31%)

726 (20%)

961 (27%)

383 (6.1%)

972 (16%)

1,920 (31%)

1,187 (19%)

1,806 (29%)

Tertiary level

101 (3.2%)

428 (14%)

950 (30%)

461 (15%)

1,179 (38%)

123 (4.3%)

493 (17%)

839 (29%)

532 (19%)

883 (31%)

224 (3.7%)

921 (15%)

1,789 (30%)

993 (17%)

2,062 (34%)

Nationality category

Swiss

223 (4.4%)

722 (14%)

1,548 (30%)

834 (16%)

1,776 (35%)

306 (5.1%)

989 (17%)

1,820 (30%)

1,189 (20%)

1,679 (28%)

529 (4.8%)

1,711 (15%)

3,368 (30%)

2,023 (18%)

3,455 (31%)

Foreigner

91 (8.2%)

153 (14%)

299 (27%)

178 (16%)

387 (35%)

138 (11%)

215 (18%)

326 (27%)

229 (19%)

316 (26%)

229 (9.8%)

368 (16%)

625 (27%)

407 (17%)

703 (30%)

Language region

German-speaking Switzerland

203 (4.6%)

572 (13%)

1,372 (31%)

699 (16%)

1,600 (36%)

240 (4.7%)

727 (14%)

1,553 (31%)

1,038 (21%)

1,505 (30%)

443 (4.7%)

1,299 (14%)

2,925 (31%)

1,737 (18%)

3,105 (33%)

French-speaking Switzerland

90 (6.7%)

233 (17%)

354 (26%)

233 (17%)

434 (32%)

152 (9.4%)

374 (23%)

439 (27%)

258 (16%)

387 (24%)

242 (8.2%)

607 (21%)

793 (27%)

491 (17%)

821 (28%)

Italian-speaking Switzerland

21 (5.0%)

70 (17%)

121 (29%)

80 (19%)

129 (31%)

52 (9.7%)

103 (19%)

154 (29%)

122 (23%)

103 (19%)

73 (7.6%)

173 (18%)

275 (29%)

202 (21%)

232 (24%)

Low back pain

No

175 (4.6%)

505 (13%)

1,082 (29%)

591 (16%)

1,414 (38%)

199 (5.5%)

551 (15%)

1,035 (29%)

721 (20%)

1,124 (31%)

374 (5.1%)

1,056 (14%)

2,117 (29%)

1,312 (18%)

2,538 (34%)

A little

97 (4.6%)

317 (15%)

676 (32%)

351 (17%)

668 (32%)

167 (5.7%)

533 (18%)

927 (31%)

572 (19%)

756 (26%)

264 (5.2%)

850 (17%)

1,603 (32%)

923 (18%)

1,424 (28%)

Strong

42 (13%)

53 (16%)

89 (27%)

70 (21%)

81 (24%)

78 (13%)

120 (19%)

184 (30%)

125 (20%)

115 (18%)

120 (13%)

173 (18%)

273 (29%)

195 (20%)

196 (20%)

1n (%)

8 Univariate distribution

Code 
DataExplorer::plot_bar(df.m.LBP.bin, parallel = TRUE, ncol = 2, nrow = 3)

8.1 Correlation analysis

Code 
DataExplorer::plot_correlation(df.m.LBP.bin)

Code 
DataExplorer::plot_correlation(df.m.LBP.bin %>% 
                                 select(lbp, pa_cat, diet, sleep, sed_beh, soc, tobacco, alcohol, drugs, dep, anxiety))

9 Validation of calculated LBP values with reported 2022 data

Code 
df.BFS.data.raw <- rio::import("../reference data/SGB_2022_prepared.xlsx")
df.BFS.data.ad <- df.BFS.data.raw %>% 
  mutate(LBP_pct = round(LBP_pct, 1),
         LBP_CI = round(LBP_CI,1))
save(df.BFS.data.ad, file = "../data/df.BFS.data.ad.rda")
Table 3: Comparison of reported data by BFS and calculated data (mean % and Margin of error: MoE).

Sex

Age category

Prevalence LBP (BFS) (%)

Prevalence LBP (BFS) (MoE)

Calculated LBP data (%)

Calculated LBP data (MoE)

Total

Total

45.1

0.8

45.1

0.8

Male

Total

40.1

1.2

40.1

1.2

Female

Total

50.0

1.1

50.0

1.1

Total

15-24

41.9

2.6

41.9

2.6

Total

25-34

45.6

2.5

45.6

2.5

Total

35-44

43.7

2.1

43.7

2.1

Total

45-54

45.0

1.9

45.1

1.9

Total

55-64

45.6

1.8

45.6

1.8

Total

65-74

45.3

2.0

45.3

2.0

Total

75+

48.8

2.2

48.8

2.2

Male

15-24

32.5

3.7

32.6

3.7

Male

25-34

40.9

3.5

40.9

3.5

Male

35-44

38.9

3.0

38.9

3.0

Male

45-54

41.5

2.8

41.5

2.8

Male

55-64

41.5

2.7

41.5

2.7

Male

65-74

42.4

3.0

42.4

3.0

Male

75+

42.1

3.3

42.1

3.3

Female

15-24

52.0

3.7

52.0

3.7

Female

25-34

50.9

3.5

50.9

3.5

Female

35-44

48.2

2.8

48.2

2.8

Female

45-54

48.5

2.6

48.5

2.6

Female

55-64

49.6

2.5

49.6

2.5

Female

65-74

47.8

2.7

47.8

2.7

Female

75+

54.4

2.9

54.4

2.9

10 Survey development for weighted statistics

Code 
# In case of only one Singleton PSU
options(survey.lonely.psu = "adjust")
df.m.LBP.bin$lbp_bin <- relevel(df.m.LBP.bin$lbp_bin, ref = "No")
df.m.LBP.bin$bmi_cat <- relevel(df.m.LBP.bin$bmi_cat, ref = "Normal weight")
df.m.LBP.bin$age_cat <- relevel(df.m.LBP.bin$age_cat, ref = "45-54 years")

# Srvy design
srvy.LBP.bin.c <- df.m.LBP.bin %>% 
  as_survey_design(weights = weighting.tel, strata = strata, nest = TRUE)
save(srvy.LBP.bin.c, file = "../data/srvy.LBP.bin.c.rda")

# Srvy design: Male
srvy.male <- srvy.LBP.bin.c %>% 
  filter(sex == "Male") %>%
  mutate(bmi_cat = relevel(factor(bmi_cat), ref = "Normal weight"),
         age_cat = relevel(factor(age_cat), ref = "45-54 years"))

save(srvy.male, file = "../data/srvy.male.rda")

# Srvy design: Female
srvy.female <- srvy.LBP.bin.c %>% 
  filter(sex == "Female") %>% 
  mutate(bmi_cat = relevel(factor(bmi_cat), ref = "Normal weight"),
         age_cat = relevel(factor(age_cat), ref = "45-54 years"))

save(srvy.female, file = "../data/srvy.female.rda")
Code 
library(survey)

survey.LBP.bin.c <- svydesign(
  weights = ~weighting.tel,
  strata = ~strata,
  ids = ~1,
  nest = TRUE,
  data = df.m.LBP.bin)

11 Summary (weighted) statistics

Code 
srvy.LBP.bin.c <- srvy.LBP.bin.c %>%
  dplyr::mutate(
    age_cat = factor(
      age_cat,
      levels = c(
        "18-24 years",
        "25-34 years",
        "35-44 years",
        "45-54 years",
        "55-64 years",
        "65-74 years",
        "75+ years"
      )
    )
  )

tbl.LBP.bin.conf.c <- srvy.LBP.bin.c %>%
  select(
    sex, age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, lbp_bin, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = sex,
    percent = "column",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"),
    label = list(
      age = "Age",
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      lbp_bin = "Low back pain - binary",
      pa_cat = "Physical activity category"),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2,
    include = everything()
  ) %>%
  gtsummary::add_overall(last = TRUE) %>%
  bold_labels()

tbl.LBP.bin.conf.c %>% 
  as_flex_table() %>% 
  set_caption("Table. Descriptive (weighted) statistics for low back pain, separated by sex.") %>%
  flextable::set_table_properties(width = 1, layout = "autofit") %>% 
  flextable::fontsize(size = 8, part = "all") %>% 
  flextable::font(fontname = "Times New Roman", part = "all") %>% 
  flextable::save_as_docx(path = "../tables/weighted descriptives/tbl.overall.docx")

tbl.LBP.bin.conf.c %>% 
  as_flex_table() %>% 
  autofit() %>%
  bold(part = "header") %>%
  align(j = 1, align = "left", part = "all")
Table 4: Summary (weighted) sample - overall

Characteristic

Male
N = 2,140,9861

Female
N = 2,234,3051

Overall
N = 4,375,2921

Age category

18-24 years

9.6%

9.1%

9.3%

25-34 years

18%

17%

17%

35-44 years

18%

19%

18%

45-54 years

18%

18%

18%

55-64 years

19%

18%

19%

65-74 years

13%

13%

13%

75+ years

4.7%

5.0%

4.9%

BMI category

Normal weight

46%

62%

54%

Underweight

1.1%

5.5%

3.3%

Overweight

40%

22%

30%

Obesity

13%

11%

12%

Civil status

Single

38%

33%

36%

Married

52%

51%

51%

Divorced

8.6%

12%

10%

Widowed

1.2%

3.7%

2.5%

Education

Compulsory school

7.7%

9.8%

8.8%

Upper secondary education

43%

49%

46%

Tertiary level

49%

42%

45%

Nationality category

Swiss

78%

79%

78%

Foreigner

22%

21%

22%

Language region

German-speaking Switzerland

73%

72%

73%

French-speaking Switzerland

23%

23%

23%

Italian-speaking Switzerland

4.2%

4.8%

4.5%

Low back pain

No

60%

50%

55%

A little

34%

41%

38%

Strong

5.3%

8.7%

7.1%

Low back pain - binary

40%

50%

45%

Physical activity category

Inactive

5.3%

6.3%

5.8%

Partially active

14%

17%

15%

Irregularly active

29%

30%

29%

Regularly active

16%

19%

18%

Trained

36%

28%

32%

1%

Code 
tbl.lbp.overall <- srvy.LBP.bin.c %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = lbp,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.lbp.male <- srvy.male %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = lbp,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()


tbl.lbp.fem <- srvy.female %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = lbp,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.lbp.sex <- gtsummary::tbl_merge(
  tbls = list(tbl.lbp.male, tbl.lbp.fem, tbl.lbp.overall ),
  tab_spanner = c("**Male**", "**Female**", "**Overall**")
)

tbl.lbp.sex %>% 
  as_flex_table() %>% 
  autofit() %>%
  bold(part = "header") %>%
  align(j = 1, align = "left", part = "all")
tbl.lbp.sex %>%
  as_flex_table() %>% 
  set_caption("Table. Descriptive (weighted) statistics for low back pain, separated by sex.") %>%
  flextable::set_table_properties(width = 1, layout = "autofit") %>% 
  flextable::fontsize(size = 8, part = "all") %>% 
  flextable::font(fontname = "Times New Roman", part = "all") %>% 
  flextable::save_as_docx(path = "../tables/weighted descriptives/tbl1.lbp.sex.docx")
Table 5: Summary (weighted) sample by lbp - separated by sex

Male

Female

Overall

Characteristic

No
N = 1,294,6131

A little
N = 732,1221

Strong
N = 114,2521

No
N = 1,127,6541

A little
N = 911,7601

Strong
N = 194,8911

No
N = 2,422,2671

A little
N = 1,643,8821

Strong
N = 309,1431

Age category

45-54 years

59%

35%

6.6%

50%

42%

8.1%

54%

38%

7.4%

18-24 years

69%

29%

1.7%

49%

45%

6.6%

59%

37%

4.2%

25-34 years

59%

38%

3.8%

49%

43%

8.0%

54%

40%

5.9%

35-44 years

62%

33%

4.7%

52%

39%

9.5%

57%

36%

7.2%

55-64 years

60%

33%

7.0%

50%

40%

10%

55%

36%

8.7%

65-74 years

58%

35%

6.8%

54%

39%

7.7%

56%

37%

7.3%

75+ years

59%

36%

5.3%

49%

40%

11%

54%

38%

8.4%

BMI category

Normal weight

61%

35%

4.5%

53%

39%

7.6%

56%

38%

6.3%

Underweight

59%

33%

7.9%

52%

39%

8.3%

53%

38%

8.2%

Overweight

63%

32%

5.0%

46%

45%

9.5%

57%

36%

6.6%

Obesity

52%

39%

9.1%

44%

42%

14%

48%

40%

11%

Civil status

Single

61%

35%

3.6%

51%

41%

7.0%

57%

38%

5.3%

Married

61%

33%

6.3%

51%

40%

8.8%

56%

37%

7.5%

Divorced

58%

35%

7.4%

47%

41%

12%

51%

39%

10%

Widowed

55%

40%

4.5%

50%

40%

11%

51%

40%

9.3%

Education

Compulsory school

61%

35%

3.8%

44%

42%

14%

51%

39%

9.4%

Upper secondary education

57%

36%

7.0%

49%

41%

9.8%

53%

39%

8.5%

Tertiary level

63%

33%

4.1%

53%

40%

6.4%

59%

36%

5.2%

Nationality category

Swiss

61%

34%

5.0%

51%

41%

8.3%

56%

38%

6.7%

Foreigner

60%

33%

6.5%

48%

41%

10%

54%

37%

8.4%

Language region

German-speaking Switzerland

61%

34%

5.3%

51%

40%

8.4%

56%

37%

6.9%

French-speaking Switzerland

58%

36%

5.8%

47%

43%

9.8%

52%

40%

7.9%

Italian-speaking Switzerland

64%

31%

4.3%

52%

41%

7.6%

57%

36%

6.1%

Physical activity category

Inactive

55%

31%

13%

42%

41%

17%

48%

37%

15%

Partially active

59%

35%

6.1%

46%

44%

10%

52%

40%

8.4%

Irregularly active

58%

37%

4.7%

48%

43%

8.4%

53%

40%

6.6%

Regularly active

58%

37%

5.6%

51%

40%

9.3%

54%

38%

7.7%

Trained

65%

31%

4.3%

57%

37%

5.8%

61%

34%

5.0%

1%

Code 
tbl.pa.overall <- srvy.LBP.bin.c %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = pa_cat,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.pa.male <- srvy.male %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = pa_cat,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()


tbl.pa.fem <- srvy.female %>%
  select(
    age_cat, bmi_cat, 
    civil_status, education, nat, lang_r,
    lbp, pa_cat
  ) %>%
  gtsummary::tbl_svysummary(
    by = pa_cat,
    percent = "row",
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{p}%"
    ),
    label = list(
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      education = "Education",
      nat = "Nationality category",
      lang_r = "Language region",
      lbp = "Low back pain",
      pa_cat = "Physical activity category"
    ),
    missing_text = "(Missing)",
    digits = all_continuous() ~ 2
  ) %>%
  bold_labels()

tbl.pa.sex <- gtsummary::tbl_merge(
  tbls = list(tbl.pa.male, tbl.pa.fem, tbl.pa.overall ),
  tab_spanner = c("**Male**", "**Female**", "**Overall**")
)

tbl.pa.sex %>% 
  as_flex_table() %>% 
  autofit() %>%
  bold(part = "header") %>%
  align(j = 1, align = "left", part = "all")
tbl.pa.sex %>% 
  as_flex_table() %>% 
  set_caption("Table. Descriptive (weighted) statistics for physical activity, separated by sex.") %>%
  flextable::set_table_properties(width = 1, layout = "autofit") %>% 
  flextable::fontsize(size = 8, part = "all") %>% 
  flextable::font(fontname = "Times New Roman", part = "all") %>%  
  flextable::save_as_docx(path = "../tables/weighted descriptives/tbl1.pa.sex.docx")
Table 6: Summary (weighted) sample by pa - separated by sex

Male

Female

Overall

Characteristic

Inactive
N = 112,6531

Partially active
N = 296,6961

Irregularly active
N = 613,8911

Regularly active
N = 338,5191

Trained
N = 779,2281

Inactive
N = 140,6271

Partially active
N = 380,2641

Irregularly active
N = 659,5001

Regularly active
N = 432,6811

Trained
N = 621,2341

Inactive
N = 253,2801

Partially active
N = 676,9601

Irregularly active
N = 1,273,3911

Regularly active
N = 771,1991

Trained
N = 1,400,4611

Age category

45-54 years

7.8%

16%

28%

15%

33%

5.5%

18%

32%

18%

26%

6.6%

17%

30%

17%

29%

18-24 years

3.4%

9.1%

19%

14%

55%

5.2%

18%

25%

14%

38%

4.3%

14%

22%

14%

46%

25-34 years

3.2%

14%

27%

12%

45%

6.1%

18%

29%

19%

28%

4.6%

16%

28%

15%

37%

35-44 years

5.9%

18%

31%

16%

30%

5.7%

19%

29%

20%

27%

5.8%

18%

30%

18%

28%

55-64 years

5.0%

15%

31%

17%

32%

6.8%

16%

32%

19%

26%

5.9%

16%

31%

18%

29%

65-74 years

5.1%

9.5%

32%

20%

34%

6.8%

13%

28%

22%

29%

5.9%

11%

30%

21%

32%

75+ years

6.7%

9.1%

32%

23%

28%

11%

14%

27%

29%

17%

9.2%

12%

30%

27%

23%

BMI category

Normal weight

4.4%

13%

25%

15%

42%

4.8%

16%

29%

20%

31%

4.6%

15%

27%

18%

36%

Underweight

4.4%

20%

23%

19%

33%

6.8%

18%

24%

22%

28%

6.5%

19%

24%

22%

29%

Overweight

5.3%

14%

32%

16%

34%

7.9%

19%

32%

17%

24%

6.2%

16%

32%

16%

30%

Obesity

8.2%

17%

33%

19%

23%

11%

19%

31%

21%

17%

9.7%

18%

32%

20%

20%

Civil status

Single

5.2%

14%

25%

13%

42%

5.4%

19%

27%

17%

32%

5.3%

16%

26%

15%

37%

Married

4.6%

14%

32%

17%

32%

5.9%

17%

32%

20%

25%

5.2%

15%

32%

19%

29%

Divorced

9.0%

11%

26%

16%

38%

11%

15%

27%

21%

26%

9.9%

14%

27%

19%

31%

Widowed

7.4%

10%

31%

26%

25%

6.3%

13%

27%

24%

29%

6.6%

13%

28%

25%

28%

Education

Compulsory school

13%

14%

22%

17%

34%

15%

18%

24%

21%

21%

14%

16%

23%

19%

26%

Upper secondary education

6.3%

14%

29%

17%

34%

6.2%

17%

30%

20%

27%

6.3%

15%

30%

19%

30%

Tertiary level

3.2%

14%

30%

15%

39%

4.3%

17%

30%

18%

31%

3.7%

15%

30%

16%

35%

Nationality category

Swiss

4.5%

14%

30%

16%

36%

4.8%

17%

30%

20%

28%

4.6%

16%

30%

18%

32%

Foreigner

8.0%

13%

26%

16%

37%

12%

17%

26%

18%

27%

9.9%

15%

26%

17%

32%

Language region

German-speaking Switzerland

4.9%

13%

30%

15%

37%

5.0%

15%

30%

20%

29%

4.9%

14%

30%

18%

33%

French-speaking Switzerland

6.6%

17%

25%

17%

34%

9.6%

23%

27%

16%

25%

8.2%

20%

26%

16%

30%

Italian-speaking Switzerland

5.1%

17%

29%

17%

32%

9.4%

19%

29%

22%

20%

7.4%

18%

29%

20%

25%

Low back pain

No

4.8%

14%

28%

15%

39%

5.3%

16%

28%

20%

31%

5.0%

14%

28%

17%

35%

A little

4.8%

14%

31%

17%

33%

6.3%

18%

31%

19%

25%

5.7%

16%

31%

18%

29%

Strong

13%

16%

25%

16%

29%

12%

20%

28%

21%

19%

13%

18%

27%

19%

23%

1%

12 Multinomial logistic regression

12.1 Definition of reference levels

Code 
# variables m1
variables.q1 <- c("pa_cat", "sex", "age_cat", "bmi_cat", "civil_status", "nat", "lang_r", "education")

variables.m3 <- c("pa_cat", "diet", "sleep", "soc", "sed_beh", "sex", "age_cat", "bmi_cat", "civil_status", "nat", "lang_r", "education")

variables.m4 <- c("pa_cat", "diet", "sleep", "soc", "sed_beh", "tobacco", "alcohol", "drugs", "dep", "anxiety", "aht", "hyperchol", "diabetes", "sex", "age_cat", "bmi_cat", "civil_status", "nat", "lang_r", "education")

label.names <- c(
  pa_cat   = "PA category",
  age_cat = "Age category",
  civil_status = "Civil status",
  nat = "Nationality",
  lang_r = "Language region",
  diet     = "Diet",
  sleep    = "Sleep",
  sed_beh  = "Sedentary behaviour",
  soc      = "Social relationships",
  tobacco  = "Tobacco",
  alcohol  = "Alcohol",
  drugs    = "Drugs",
  sex      = "Sex",
  bmi_cat  = "BMI category",
  dep      = "Depression",
  anxiety  = "Anxiety",
  education= "Education",
  aht      = "Art. hypertension",
  hyperchol= "Hypercholesterolaemia",
  diabetes = "Diabetes")

variables.conf <- c(var.socdem, var.soceco)

Reference levels:

Type Outcome/Variable Abbreviation Reference Value
Primary LBP lbp No
Physical activity pa_cat Inactive
Secondary Alcohol alc Abstinent
Anxiety anxiety None
Arterial hypertension aht No
Depression dep No/ minimal symptoms of depression
Diabetes mellitus diabetes No
Diet diet <5 days per week
Drugs drugs Never used drugs
Hypercholesterolaemia hyperchol No
Sedentary behaviour (numeric, hours per day) sed_beh 0
Sleeping disorder sleep None or few
Social relationships/ Social support soc Low
Tobacco tobacco Non-smoker
Confounder Age category age_cat 45-54 years
BMI bmi_cat Normal weight
Civil status civil_status Single
Education education Compulsory school
Language region lang_r German-speaking Switzerland
Nationality nat Swiss
Sex sex Male
Code 
# relevel

srvy.LBP.bin.c$variables <- srvy.LBP.bin.c$variables %>% 
  mutate(
    age_cat = fct_relevel(age_cat, "45-54 years"))

# srvy.LBP.bin.c$variables <- srvy.LBP.bin.c$variables %>%
#   mutate(
#     # Primary
#     lbp          = fct_relevel(lbp, "No"),
#     pa_cat       = fct_relevel(pa_cat, "Inactive"),
#     # Secondary
#     alc          = fct_relevel(alc, "Abstinent"),
#     anxiety      = fct_relevel(anxiety, "None"),
#     aht          = fct_relevel(aht, "No"),
#     dep          = fct_relevel(dep, "No/ minimal symptoms of depression"),
#     diabetes     = fct_relevel(diabetes, "No"),
#     diet         = fct_relevel(diet, "<5 days per week"),
#     drugs        = fct_relevel(drugs, "Never used drugs"),
#     hyperchol    = fct_relevel(hyperchol, "No"),
#     # sed beh = numeric
#     sleep        = fct_relevel(sleep, "None or few"),
#     soc          = fct_relevel(soc, "Low"),
#     tobacco      = fct_relevel(tobacco, "Non-smoker"),
#     # Confounder
#     age_cat      = fct_relevel(age_cat, "45-54 years"),
#     bmi_cat      = fct_relevel(bmi_cat, "Normal weight"),
#     civil_status = fct_relevel(civil_status, "Single"),
#     education    = fct_relevel(education, "Compulsory school"),
#     lang_r       = fct_relevel(lang_r, "German-speaking Switzerland"),
#     nat          = fct_relevel(nat, "Swiss"),
#     sex          = fct_relevel(sex, "Male"))
Code 
adapt.models <- function(df, variables) {
  df %>% 
    filter(!str_detect(Parameter, "^\\(Intercept\\)")) %>% 
    mutate(
      `Odds Ratio` = round(Coefficient, 2),
      SE = round(SE, 2),
      CI_low = round(CI_low, 2),
      CI_high = round(CI_high, 2),
      `95% CI` = paste0("[", CI_low, "; ", CI_high, "]"),
      p = ifelse(p == 0, "<0.01", sprintf("%.3f", p)),
      LBP = case_when(
        str_detect(Parameter, ":1$") ~ "A little LBP",
        str_detect(Parameter, ":2$") ~ "Strong LBP",
        TRUE ~ NA_character_
      ),
      Parameter = str_replace(Parameter, ":1$", ""),
      Parameter = str_replace(Parameter, ":2$", "")
    ) %>%
    mutate(
      Variable = str_extract(Parameter, paste0(variables, collapse = "|")),
      Parameter = str_trim(str_remove(Parameter, paste0("^(", paste(variables, collapse = "|"), ")_?")))
    ) %>%
    arrange(LBP) %>%
    mutate(
      LBP = ifelse(duplicated(LBP), NA, LBP)
    ) %>%
    select(LBP, Variable, Parameter, `Odds Ratio`, `95% CI`, p)
}

12.2 Model 1: Simple model q1 with LBP and PA

Code 
vglm.q1 <- svy_vglm(
  lbp ~ pa_cat,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.LBP.bin.c))

# reference categories: PA: inactive, LBP: no LBP, x1 = A little LBP, x2 = Strong LBP
# parameters::parameters(vglm.model.pa, exponentiate = TRUE,  digits = 3)
m1.q1.raw <- model_parameters(vglm.q1,
                              exponentiate = TRUE,
                              drop = "Intercept",
                              digits = 3)

# adjusted for sociodemographic variables
vglm.q1.conf <- svy_vglm(
  lbp ~ pa_cat + sex + age_cat + bmi_cat + civil_status + nat + lang_r + education,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.LBP.bin.c)
)

m1.q1.adj <- parameters::model_parameters(vglm.q1.conf, 
                              exponentiate = TRUE,
                              drop = "Intercept",
                              digits = 3)

\[ \begin{aligned} \log{\left(\frac{P\left(LBP_i=k\mid X_i\right)}{P\left(LBP_i=0\mid X_i\right)}\right)}= \beta_0^{\left(k\right)}+\beta_1^{\left(k\right)}\cdot\mathrm{P}\mathrm{A}_{\mathrm{i,Partially}}+\beta_2^{\left(k\right)}\cdot\mathrm{P}\mathrm{A}_{\mathrm{i,Irregular}}+\beta_3^{\left(k\right)}\cdot\mathrm{P}\mathrm{A}_{\mathrm{i,Regular}}+\beta_4^{\left(k\right)}\cdot\mathrm{P}\mathrm{A}_{\mathrm{i,Trained}} \end{aligned} \]

Code 
# unadjusted model
m1.q1.adapted <- adapt.models(m1.q1.raw, variables.q1)

m1.q1.adapted <- m1.q1.adapted %>%
  mutate(
    `crude OR` = `Odds Ratio`,
    `crude 95% CI` = `95% CI`,
    `crude p value` = p) %>%
  select(LBP, Variable, Parameter, `crude OR`, `crude 95% CI`, `crude p value`) %>% 
  mutate(
    Variable = dplyr::recode(Variable,
      pa_cat = "Physical activity level"))

ft.m1.q1.adapted <- flextable(m1.q1.adapted)
ft.m1.q1.adapted <- ft.m1.q1.adapted %>%
  set_caption("Table. Model 1 on association of other factors with LBP - unadjusted") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all")
ft.m1.q1.adapted
save(m1.q1.adapted, file = "../tables/models/m1.q1.unadjusted.rda")

ft.m1.q1.adapted %>%
  set_caption("Model 1 on association of other factors with LBP - unadjusted") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  flextable::line_spacing(space = 1.0) %>%
  flextable::fontsize(size = 8, part = "all") %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::save_as_docx(path = "../tables/models/model1.q1.unadjusted.docx")
Table 7: Model 1 on association of other factors with LBP.

LBP

Variable

Parameter

crude OR

crude 95% CI

crude p value

A little LBP

Physical activity level

Partially active

1.01

[0.8; 1.27]

0.938

Physical activity level

Irregularly active

0.99

[0.8; 1.23]

0.953

Physical activity level

Regularly active

0.93

[0.74; 1.17]

0.538

Physical activity level

Trained

0.72

[0.58; 0.9]

0.003

Strong LBP

Physical activity level

Partially active

0.51

[0.38; 0.7]

0.000

Physical activity level

Irregularly active

0.39

[0.29; 0.53]

0.000

Physical activity level

Regularly active

0.45

[0.33; 0.61]

0.000

Physical activity level

Trained

0.25

[0.19; 0.35]

0.000

\[ \begin{aligned} \log\left(\frac{P(LBP_i = k \mid X_i)}{P(LBP_i = 0 \mid X_i)}\right) &= \beta_0^{(k)} + \beta_1^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Partially}} + \beta_2^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Irregular}} + \beta_3^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Regular}} \\ &\quad + \beta_4^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Trained}} + \beta_5^{(k)} \cdot \mathrm{Sex}_{i,\mathrm{Female}} + \beta_6^{(k)} \cdot \mathrm{Age}_{i,\mathrm{25\text{-}34}} + \beta_7^{(k)} \cdot \mathrm{Age}_{i,\mathrm{35\text{-}44}} \\ &\quad + \beta_8^{(k)} \cdot \mathrm{Age}_{i,\mathrm{45\text{-}54}} + \beta_9^{(k)} \cdot \mathrm{Age}_{i,\mathrm{55\text{-}64}} + \beta_{10}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{65\text{-}74}} + \beta_{11}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{75+}} \\ &\quad + \beta_{12}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Underweight}} + \beta_{13}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Overweight}} + \beta_{14}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Obese}} \\ &\quad + \beta_{15}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Married}} + \beta_{16}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Divorced}} + \beta_{17}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Widowed}} \\ &\quad + \beta_{18}^{(k)} \cdot \mathrm{Education}_{i,\mathrm{UpperSec}} + \beta_{19}^{(k)} \cdot \mathrm{Education}_{i,\mathrm{Tertiary}} \\ &\quad + \beta_{20}^{(k)} \cdot \mathrm{Region}_{i,\mathrm{French}} + \beta_{21}^{(k)} \cdot \mathrm{Region}_{i,\mathrm{Italian}} \end{aligned} \]

Code 
# adjusted model
m1.q1.adj.adapted <- adapt.models(m1.q1.adj, variables.q1)

m1.q1.adj.adapted <- m1.q1.adj.adapted %>%
  mutate(
    `adjusted OR` = `Odds Ratio`,
    `adjusted 95% CI` = `95% CI`,
    `adjusted p value` = p) %>%
  select(LBP, Variable, Parameter, `adjusted OR`, `adjusted 95% CI`, `adjusted p value`) %>% 
  mutate(
    Variable = dplyr::recode(Variable,
      pa_cat = "Physical activity level",
      sex = "Sex",
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      nat = "Nationality",
      lang_r = "Language region",
      education = "Education",
      income = "Monthly household income"))

ft.m1.q1.adj.adapted <- flextable(m1.q1.adj.adapted)
ft.m1.q1.adj.adapted <- ft.m1.q1.adj.adapted %>%
  set_caption("Table. Model 1 on association of other factors with LBP - adjusted") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all")
ft.m1.q1.adj.adapted

LBP

Variable

Parameter

adjusted OR

adjusted 95% CI

adjusted p value

A little LBP

Physical activity level

Partially active

1.05

[0.84; 1.33]

0.658

Physical activity level

Irregularly active

1.07

[0.86; 1.33]

0.552

Physical activity level

Regularly active

0.98

[0.78; 1.24]

0.873

Physical activity level

Trained

0.80

[0.64; 1]

0.051

Sex

Female

1.40

[1.27; 1.54]

0.000

Age category

18-24 years

0.89

[0.71; 1.11]

0.291

Age category

25-34 years

1.13

[0.95; 1.35]

0.177

Age category

35-44 years

0.92

[0.79; 1.07]

0.269

Age category

55-64 years

0.92

[0.8; 1.06]

0.232

Age category

65-74 years

0.91

[0.79; 1.06]

0.233

Age category

75+ years

0.94

[0.77; 1.14]

0.531

BMI category

Underweight

0.97

[0.75; 1.24]

0.783

BMI category

Overweight

1.01

[0.91; 1.13]

0.803

BMI category

Obesity

1.25

[1.08; 1.45]

0.004

Civil status

Married

0.98

[0.86; 1.12]

0.781

Civil status

Divorced

1.10

[0.92; 1.33]

0.292

Civil status

Widowed

1.07

[0.8; 1.42]

0.652

Nationality

Foreigner

1.00

[0.89; 1.12]

0.976

Language region

French-speaking Switzerland

1.14

[1.03; 1.26]

0.012

Language region

Italian-speaking Switzerland

0.93

[0.8; 1.09]

0.366

Education

Upper secondary education

0.97

[0.82; 1.15]

0.732

Education

Tertiary level

0.83

[0.69; 0.98]

0.033

Strong LBP

Physical activity level

Partially active

0.61

[0.44; 0.83]

0.002

Physical activity level

Irregularly active

0.48

[0.36; 0.65]

0.000

Physical activity level

Regularly active

0.52

[0.38; 0.71]

0.000

Physical activity level

Trained

0.35

[0.26; 0.49]

0.000

Sex

Female

1.81

[1.51; 2.16]

0.000

Age category

18-24 years

0.64

[0.41; 1.01]

0.053

Age category

25-34 years

1.10

[0.77; 1.55]

0.604

Age category

35-44 years

1.01

[0.76; 1.35]

0.921

Age category

55-64 years

1.12

[0.86; 1.46]

0.392

Age category

65-74 years

0.90

[0.68; 1.2]

0.472

Age category

75+ years

0.98

[0.68; 1.42]

0.928

BMI category

Underweight

1.18

[0.74; 1.87]

0.493

BMI category

Overweight

1.03

[0.85; 1.26]

0.744

BMI category

Obesity

1.80

[1.41; 2.29]

0.000

Civil status

Married

1.26

[0.99; 1.62]

0.064

Civil status

Divorced

1.70

[1.19; 2.41]

0.003

Civil status

Widowed

1.43

[0.89; 2.28]

0.135

Nationality

Foreigner

1.24

[0.99; 1.54]

0.056

Language region

French-speaking Switzerland

1.16

[0.96; 1.4]

0.115

Language region

Italian-speaking Switzerland

0.77

[0.57; 1.04]

0.084

Education

Upper secondary education

1.06

[0.8; 1.4]

0.687

Education

Tertiary level

0.61

[0.45; 0.82]

0.001

Code 
save(m1.q1.adj.adapted, file = "../tables/models/m1.q1.adjusted.rda")

ft.m1.q1.adj.adapted %>%
  set_caption("Model 1 on association of other factors with LBP - adjusted") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  flextable::line_spacing(space = 1.0) %>%
  flextable::fontsize(size = 8, part = "all") %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::save_as_docx(path = "../tables/models/model1.q1.adjusted.docx")

12.3 Model 2: Simple model q1, separated by sex

\[ \begin{aligned} \log\left(\frac{P(LBP_i = k \mid X_i)}{P(LBP_i = 0 \mid X_i)}\right) &= \beta_0^{(k)} + \beta_1^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Partially}} + \beta_2^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Irregular}} + \beta_3^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Regular}} \\ &\quad + \beta_4^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Trained}} + \beta_5^{(k)} \cdot \mathrm{Age}_{i,\mathrm{25\text{-}34}} + \beta_6^{(k)} \cdot \mathrm{Age}_{i,\mathrm{35\text{-}44}} \\ &\quad + \beta_7^{(k)} \cdot \mathrm{Age}_{i,\mathrm{45\text{-}54}} + \beta_8^{(k)} \cdot \mathrm{Age}_{i,\mathrm{55\text{-}64}} + \beta_{9}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{65\text{-}74}} + \beta_{10}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{75+}} \\ &\quad + \beta_{11}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Underweight}} + \beta_{12}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Overweight}} + \beta_{13}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Obese}} \\ &\quad + \beta_{14}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Married}} + \beta_{15}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Divorced}} + \beta_{16}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Widowed}} \\ &\quad + \beta_{17}^{(k)} \cdot \mathrm{Education}_{i,\mathrm{UpperSec}} + \beta_{18}^{(k)} \cdot \mathrm{Education}_{i,\mathrm{Tertiary}} \\ &\quad + \beta_{19}^{(k)} \cdot \mathrm{Region}_{i,\mathrm{French}} + \beta_{20}^{(k)} \cdot \mathrm{Region}_{i,\mathrm{Italian}} \end{aligned} \]

Code 
# Model for male
vglm.q1.male <- svy_vglm(
  lbp ~ pa_cat + age_cat + bmi_cat + civil_status + nat + lang_r + education,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.male))

m2.q1.male.raw <- model_parameters(vglm.q1.male, 
                                   exponentiate = TRUE,
                                  drop = "Intercept",
                                  digits = 3)

svytable(~ bmi_cat + lbp, design = as_survey_design(srvy.male))
               lbp
bmi_cat                 No   A little     Strong
  Normal weight 597329.916 343909.174  43896.963
  Underweight    13363.510   7538.563   1804.473
  Overweight    536504.780 270733.214  42663.464
  Obesity       147414.369 109941.053  25886.967
Code 
# Small sample for BMI cat underweight adn strong LBP
Code 
# Model for female
vglm.q1.female <- svy_vglm(
  lbp ~ pa_cat + age_cat + bmi_cat + civil_status + nat + lang_r + education,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.female))

m2.q1.female.raw <- model_parameters(vglm.q1.female,
                                      exponentiate = TRUE,
                                      drop = "Intercept",
                                      digits = 3)
Code 
m2.q1.male.adapted <- adapt.models(m2.q1.male.raw, variables.q1)
m2.q1.female.adapted <- adapt.models(m2.q1.female.raw, variables.q1)

m2.q1.male.adapted <- m2.q1.male.adapted %>% 
  mutate(
    `adjusted OR (male)` = `Odds Ratio`,
    `95% CI (male)` = `95% CI`,
    `p value (male)` = p) %>% 
  select(LBP, Variable, Parameter, `adjusted OR (male)`, `95% CI (male)`, `p value (male)`)

m2.q1.female.adapted <- m2.q1.female.adapted %>% 
  mutate(
    `adjusted OR (female)` = `Odds Ratio`,
    `95% CI (female)` = `95% CI`,
    `p value (female)` = p) %>% 
  select(LBP, Variable, Parameter, `adjusted OR (female)`, `95% CI (female)`, `p value (female)`)

m2.q1.combined <- dplyr::bind_cols(
  m2.q1.male.adapted,
  m2.q1.female.adapted %>%
    select(`adjusted OR (female)`, `95% CI (female)`, `p value (female)`)) %>% 
  mutate(
    Variable = dplyr::recode(Variable,
      pa_cat = "Physical activity level",
      sex = "Sex",
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      nat = "Nationality",
      lang_r = "Language region",
      education = "Education",
      income = "Monthly household income"))

ft.m2.q1.combined <- flextable(m2.q1.combined)

ft.m2.q1.combined <- ft.m2.q1.combined %>%
  set_caption("Table. Adjusted model on association of PA with LBP by sex") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all")
ft.m2.q1.combined
save(m2.q1.combined, file = "../tables/models/m2.q1.combined.rda")

flextable::flextable(m2.q1.combined) %>%
  set_caption("Table. Adjusted model on association of PA with LBP by sex") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  flextable::line_spacing(space = 1.0) %>%
  flextable::fontsize(size = 8, part = "all") %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::save_as_docx(path = "../tables/models/model2.q1.by.sex.docx")
Table 8: Model 2: Adjusted model on association of PA with LBP by sex.

LBP

Variable

Parameter

adjusted OR (male)

95% CI (male)

p value (male)

adjusted OR (female)

95% CI (female)

p value (female)

A little LBP

Physical activity level

Partially active

1.11

[0.77; 1.61]

0.568

1.03

[0.77; 1.38]

0.860

Physical activity level

Irregularly active

1.22

[0.86; 1.74]

0.261

0.99

[0.75; 1.31]

0.960

Physical activity level

Regularly active

1.19

[0.83; 1.73]

0.343

0.88

[0.66; 1.17]

0.375

Physical activity level

Trained

0.93

[0.66; 1.33]

0.705

0.73

[0.55; 0.97]

0.032

Age category

18-24 years

0.61

[0.43; 0.87]

0.007

1.20

[0.89; 1.62]

0.222

Age category

25-34 years

1.06

[0.81; 1.38]

0.679

1.17

[0.93; 1.48]

0.190

Age category

35-44 years

0.90

[0.72; 1.13]

0.373

0.93

[0.77; 1.14]

0.491

Age category

55-64 years

0.92

[0.75; 1.13]

0.419

0.91

[0.76; 1.1]

0.332

Age category

65-74 years

1.01

[0.81; 1.25]

0.964

0.83

[0.68; 1.01]

0.067

Age category

75+ years

0.99

[0.74; 1.32]

0.937

0.88

[0.66; 1.16]

0.351

BMI category

Underweight

0.94

[0.49; 1.79]

0.842

1.01

[0.77; 1.32]

0.933

BMI category

Overweight

0.83

[0.71; 0.97]

0.019

1.28

[1.1; 1.49]

0.001

BMI category

Obesity

1.18

[0.95; 1.46]

0.130

1.25

[1.01; 1.53]

0.037

Civil status

Married

0.89

[0.74; 1.08]

0.236

1.08

[0.91; 1.28]

0.394

Civil status

Divorced

1.00

[0.75; 1.33]

0.996

1.22

[0.96; 1.56]

0.100

Civil status

Widowed

1.17

[0.68; 2.03]

0.569

1.16

[0.82; 1.62]

0.401

Nationality

Foreigner

0.95

[0.8; 1.14]

0.589

1.03

[0.88; 1.22]

0.685

Language region

French-speaking Switzerland

1.15

[0.99; 1.34]

0.076

1.12

[0.98; 1.29]

0.108

Language region

Italian-speaking Switzerland

0.86

[0.68; 1.09]

0.219

0.99

[0.8; 1.21]

0.900

Education

Upper secondary education

0.99

[0.75; 1.31]

0.958

0.93

[0.75; 1.16]

0.511

Education

Tertiary level

0.78

[0.59; 1.03]

0.083

0.84

[0.67; 1.06]

0.145

Strong LBP

Physical activity level

Partially active

0.48

[0.28; 0.81]

0.006

0.70

[0.48; 1.03]

0.073

Physical activity level

Irregularly active

0.38

[0.23; 0.61]

0.000

0.56

[0.39; 0.82]

0.003

Physical activity level

Regularly active

0.43

[0.26; 0.71]

0.001

0.58

[0.39; 0.86]

0.007

Physical activity level

Trained

0.35

[0.21; 0.58]

0.000

0.35

[0.24; 0.53]

0.000

Age category

18-24 years

0.27

[0.11; 0.65]

0.003

1.02

[0.59; 1.75]

0.954

Age category

25-34 years

0.78

[0.44; 1.38]

0.387

1.38

[0.89; 2.15]

0.150

Age category

35-44 years

0.72

[0.44; 1.17]

0.180

1.28

[0.9; 1.83]

0.169

Age category

55-64 years

1.09

[0.72; 1.67]

0.680

1.16

[0.83; 1.61]

0.383

Age category

65-74 years

1.08

[0.69; 1.69]

0.734

0.78

[0.54; 1.13]

0.184

Age category

75+ years

0.84

[0.48; 1.49]

0.555

1.03

[0.64; 1.66]

0.909

BMI category

Underweight

2.20

[0.61; 7.96]

0.230

1.12

[0.69; 1.83]

0.653

BMI category

Overweight

0.86

[0.63; 1.17]

0.340

1.23

[0.95; 1.6]

0.114

BMI category

Obesity

1.68

[1.14; 2.46]

0.009

1.77

[1.29; 2.43]

0.000

Civil status

Married

1.36

[0.93; 1.99]

0.109

1.22

[0.88; 1.68]

0.227

Civil status

Divorced

1.45

[0.77; 2.74]

0.249

1.86

[1.22; 2.84]

0.004

Civil status

Widowed

1.02

[0.3; 3.47]

0.970

1.62

[0.95; 2.77]

0.078

Nationality

Foreigner

1.38

[0.96; 1.98]

0.085

1.18

[0.9; 1.54]

0.229

Language region

French-speaking Switzerland

1.12

[0.82; 1.51]

0.483

1.17

[0.93; 1.48]

0.173

Language region

Italian-speaking Switzerland

0.70

[0.41; 1.19]

0.187

0.82

[0.56; 1.19]

0.292

Education

Upper secondary education

2.27

[1.26; 4.08]

0.006

0.78

[0.56; 1.08]

0.139

Education

Tertiary level

1.17

[0.64; 2.13]

0.605

0.47

[0.32; 0.67]

0.000

12.4 Model 3: Adjusted model q2 including lifestyle factors and confounder, separated by sex

\[ \begin{aligned} \log\left(\frac{P(LBP_i = k \mid X_i)}{P(LBP_i = 0 \mid X_i)}\right) &= \beta_0^{(k)} + \beta_1^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Partially}} + \beta_2^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Irregular}} + \beta_3^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Regular}} + \beta_4^{(k)} \cdot \mathrm{PA}_{i,\mathrm{Trained}} \\ &\quad + \beta_5^{(k)} \cdot \mathrm{Diet}_{i,\mathrm{0\text{-}2\,portions}} + \beta_6^{(k)} \cdot \mathrm{Diet}_{i,\mathrm{3\text{-}4\,portions}} + \beta_7^{(k)} \cdot \mathrm{Diet}_{i,\geq\mathrm{5\,portions}} \\ &\quad + \beta_8^{(k)} \cdot \mathrm{Sleep}_{i,\mathrm{Moderate}} + \beta_9^{(k)} \cdot \mathrm{Sleep}_{i,\mathrm{Pathological}} \\ &\quad + \beta_{10}^{(k)} \cdot \mathrm{SocialRelationships}_{i,\mathrm{Medium}} + \beta_{11}^{(k)} \cdot \mathrm{SocialRelationships}_{i,\mathrm{Strong}} \\ &\quad + \beta_{12}^{(k)} \cdot \mathrm{SedentaryBehaviour}_i \\ &\quad + \beta_{13}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{25\text{-}34}} + \beta_{14}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{35\text{-}44}} + \beta_{15}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{45\text{-}54}} \\ &\quad + \beta_{16}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{55\text{-}64}} + \beta_{17}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{65\text{-}74}} + \beta_{18}^{(k)} \cdot \mathrm{Age}_{i,\mathrm{75+}} \\ &\quad + \beta_{19}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Underweight}} + \beta_{20}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Overweight}} + \beta_{21}^{(k)} \cdot \mathrm{BMI}_{i,\mathrm{Obesity}} \\ &\quad + \beta_{22}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Married}} + \beta_{23}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Divorced}} + \beta_{24}^{(k)} \cdot \mathrm{CivilStatus}_{i,\mathrm{Widowed}} \\ &\quad + \beta_{25}^{(k)} \cdot \mathrm{Nationality}_{i,\mathrm{Foreigner}} \\ &\quad + \beta_{26}^{(k)} \cdot \mathrm{Region}_{i,\mathrm{French}} + \beta_{27}^{(k)} \cdot \mathrm{Region}_{i,\mathrm{Italian}} \\ &\quad + \beta_{28}^{(k)} \cdot \mathrm{Education}_{i,\mathrm{UpperSec}} + \beta_{29}^{(k)} \cdot \mathrm{Education}_{i,\mathrm{Tertiary}} \end{aligned} \]

Code 
vglm.q2.reduced <- svy_vglm(
  lbp ~ pa_cat + diet + sleep + soc + sed_beh + sex + age_cat + bmi_cat + civil_status + nat + lang_r + education, 
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.LBP.bin.c))

m3.q2.raw <- parameters::model_parameters(vglm.q2.reduced, 
                                          exponentiate = TRUE,
                                          drop = "Intercept",
                                          digits = 3)
Code 
vglm.q2.male <- svy_vglm(
  lbp ~ pa_cat + diet + sleep + soc + sed_beh + age_cat + bmi_cat + civil_status + nat + lang_r + education, 
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.male))

m3.q2.male <- parameters::model_parameters(vglm.q2.male, 
                                          exponentiate = TRUE,
                                          drop = "Intercept",
                                          digits = 3)

m3.q2.male.adapted <- adapt.models(m3.q2.male, variables.m3)
Code 
vglm.q2.female <- svy_vglm(
  lbp ~ pa_cat + diet + sleep + soc + sed_beh + age_cat + bmi_cat + civil_status + nat + lang_r + education, 
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.female))

m3.q2.female <- parameters::model_parameters(vglm.q2.female, 
                                          exponentiate = TRUE,
                                          drop = "Intercept",
                                          digits = 3)

m3.q2.female.adapted <- adapt.models(m3.q2.female, variables.m3)
Code 
m3.q2.male.adapted <- m3.q2.male.adapted %>% 
  mutate(
    `Adjusted OR (male)` = `Odds Ratio`,
    `95% CI (male)` = `95% CI`,
    `p (male)` = p) %>% 
  select(LBP, Variable, Parameter, `Adjusted OR (male)`, `95% CI (male)`, `p (male)`)

m3.q2.female.adapted <- m3.q2.female.adapted %>% 
  mutate(
    `Adjusted OR (female)` = `Odds Ratio`,
    `95% CI (female)` = `95% CI`,
    `p (female)` = p) %>% 
  select(LBP, Variable, Parameter, `Adjusted OR (female)`, `95% CI (female)`, `p (female)`)

m3.q2.combined <- dplyr::bind_cols(
  m3.q2.male.adapted,
  m3.q2.female.adapted %>%
    select(`Adjusted OR (female)`, `95% CI (female)`, `p (female)`)) %>% 
  mutate(
    Variable = dplyr::recode(Variable,
      pa_cat = "Physical activity level",
      sex = "Sex",
      diet = "Diet",
      sleep = "Sleep",
      soc = "Social relationships",
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      nat = "Nationality",
      sed_beh = "Sedentary behaviour",
      lang_r = "Language region",
      education = "Education",
      income = "Monthly household income"))



ft.m3.q2.combined <- flextable(m3.q2.combined)

ft.m3.q2.combined <- ft.m3.q2.combined %>%
  set_caption("Table. Adjusted model on association of lifestyle factors with LBP by sex") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all")
ft.m3.q2.combined
save(m3.q2.combined, file = "../tables/models/m3.q2.combined.rda")

flextable::flextable(m3.q2.combined) %>%
  set_caption("Table. Adjusted model on association of lifestyle factors with LBP by sex") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  flextable::line_spacing(space = 1.0) %>%
  flextable::fontsize(size = 8, part = "all") %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::save_as_docx(path = "../tables/models/model3.q2.by.sex.docx")
Table 9: Model 3: Adjusted model on association of lifestyle factors with LBP by sex.

LBP

Variable

Parameter

Adjusted OR (male)

95% CI (male)

p (male)

Adjusted OR (female)

95% CI (female)

p (female)

A little LBP

Physical activity level

Partially active

1.17

[0.81; 1.7]

0.404

1.06

[0.79; 1.43]

0.692

Physical activity level

Irregularly active

1.30

[0.92; 1.85]

0.140

1.05

[0.79; 1.4]

0.712

Physical activity level

Regularly active

1.24

[0.86; 1.8]

0.248

0.93

[0.69; 1.25]

0.624

Physical activity level

Trained

0.96

[0.67; 1.36]

0.802

0.79

[0.6; 1.06]

0.116

Diet

0-2 portions

1.04

[0.85; 1.29]

0.689

0.92

[0.7; 1.22]

0.569

Diet

3-4 portions

1.29

[1.03; 1.61]

0.029

0.88

[0.67; 1.17]

0.383

Diet

≥5 portions

1.12

[0.84; 1.49]

0.440

0.90

[0.67; 1.2]

0.461

Sleep

moderate

1.63

[1.38; 1.92]

0.000

1.31

[1.14; 1.52]

0.000

Sleep

pathological

1.63

[1.21; 2.2]

0.001

1.57

[1.26; 1.96]

0.000

Social relationships

medium

0.93

[0.71; 1.23]

0.620

0.78

[0.61; 0.99]

0.043

Social relationships

strong

0.87

[0.66; 1.14]

0.306

0.70

[0.54; 0.89]

0.004

Sedentary behaviour

0.99

[0.97; 1.02]

0.647

1.01

[0.99; 1.04]

0.253

Age category

18-24 years

0.65

[0.45; 0.93]

0.017

1.22

[0.91; 1.63]

0.193

Age category

25-34 years

1.10

[0.84; 1.44]

0.496

1.19

[0.94; 1.5]

0.150

Age category

35-44 years

0.91

[0.72; 1.14]

0.417

0.96

[0.79; 1.18]

0.722

Age category

55-64 years

0.91

[0.74; 1.12]

0.384

0.90

[0.74; 1.08]

0.262

Age category

65-74 years

0.97

[0.78; 1.21]

0.784

0.84

[0.69; 1.03]

0.094

Age category

75+ years

0.90

[0.67; 1.21]

0.496

0.88

[0.66; 1.16]

0.365

BMI category

Underweight

0.82

[0.42; 1.61]

0.562

1.00

[0.77; 1.31]

0.981

BMI category

Overweight

0.82

[0.7; 0.96]

0.012

1.28

[1.1; 1.49]

0.001

BMI category

Obesity

1.18

[0.95; 1.46]

0.142

1.25

[1.01; 1.53]

0.036

Civil status

Married

0.90

[0.75; 1.1]

0.308

1.10

[0.92; 1.32]

0.279

Civil status

Divorced

1.01

[0.76; 1.34]

0.963

1.24

[0.97; 1.58]

0.087

Civil status

Widowed

1.21

[0.69; 2.12]

0.513

1.18

[0.84; 1.66]

0.343

Nationality

Foreigner

0.91

[0.76; 1.09]

0.323

1.00

[0.85; 1.17]

0.965

Language region

French-speaking Switzerland

1.07

[0.92; 1.26]

0.368

1.07

[0.93; 1.23]

0.342

Language region

Italian-speaking Switzerland

0.83

[0.65; 1.05]

0.126

0.95

[0.77; 1.17]

0.642

Education

Upper secondary education

1.01

[0.77; 1.33]

0.941

0.94

[0.76; 1.17]

0.589

Education

Tertiary level

0.79

[0.59; 1.04]

0.096

0.86

[0.68; 1.08]

0.192

Strong LBP

Physical activity level

Partially active

0.63

[0.36; 1.1]

0.102

0.79

[0.53; 1.18]

0.245

Physical activity level

Irregularly active

0.53

[0.31; 0.9]

0.019

0.69

[0.46; 1.01]

0.057

Physical activity level

Regularly active

0.62

[0.35; 1.07]

0.085

0.68

[0.45; 1.02]

0.065

Physical activity level

Trained

0.48

[0.27; 0.85]

0.012

0.44

[0.29; 0.66]

0.000

Diet

0-2 portions

0.83

[0.54; 1.28]

0.398

0.82

[0.52; 1.29]

0.390

Diet

3-4 portions

0.86

[0.53; 1.39]

0.538

0.76

[0.48; 1.19]

0.232

Diet

≥5 portions

1.16

[0.65; 2.06]

0.615

0.95

[0.59; 1.55]

0.845

Sleep

moderate

2.19

[1.58; 3.05]

0.000

2.43

[1.91; 3.11]

0.000

Sleep

pathological

4.06

[2.64; 6.24]

0.000

5.78

[4.27; 7.82]

0.000

Social relationships

medium

0.50

[0.32; 0.79]

0.003

0.53

[0.37; 0.76]

0.001

Social relationships

strong

0.43

[0.27; 0.69]

0.000

0.43

[0.3; 0.62]

0.000

Sedentary behaviour

1.01

[0.96; 1.06]

0.668

1.02

[0.97; 1.06]

0.454

Age category

18-24 years

0.34

[0.14; 0.84]

0.020

1.15

[0.67; 1.99]

0.610

Age category

25-34 years

0.87

[0.48; 1.57]

0.653

1.51

[0.98; 2.34]

0.064

Age category

35-44 years

0.75

[0.46; 1.22]

0.244

1.44

[1; 2.08]

0.048

Age category

55-64 years

1.09

[0.7; 1.68]

0.705

1.12

[0.8; 1.58]

0.517

Age category

65-74 years

1.07

[0.68; 1.7]

0.759

0.83

[0.57; 1.22]

0.337

Age category

75+ years

0.75

[0.42; 1.37]

0.352

1.11

[0.67; 1.84]

0.688

BMI category

Underweight

1.66

[0.4; 6.89]

0.484

1.08

[0.64; 1.81]

0.773

BMI category

Overweight

0.84

[0.61; 1.15]

0.284

1.23

[0.95; 1.6]

0.123

BMI category

Obesity

1.60

[1.09; 2.36]

0.016

1.71

[1.23; 2.37]

0.001

Civil status

Married

1.54

[1.04; 2.29]

0.030

1.26

[0.92; 1.74]

0.154

Civil status

Divorced

1.42

[0.76; 2.66]

0.273

1.84

[1.21; 2.81]

0.005

Civil status

Widowed

1.15

[0.35; 3.78]

0.817

1.65

[0.96; 2.84]

0.072

Nationality

Foreigner

1.21

[0.83; 1.78]

0.319

1.08

[0.82; 1.42]

0.573

Language region

French-speaking Switzerland

1.00

[0.73; 1.36]

0.991

1.03

[0.81; 1.31]

0.820

Language region

Italian-speaking Switzerland

0.66

[0.39; 1.13]

0.130

0.73

[0.5; 1.07]

0.106

Education

Upper secondary education

2.29

[1.23; 4.26]

0.009

0.81

[0.58; 1.14]

0.223

Education

Tertiary level

1.15

[0.61; 2.18]

0.662

0.49

[0.33; 0.72]

0.000

12.5 Model 4: Model q2 including lifestyle factors and Model Hartvigsen et al. 2018

\[ \begin{aligned} \log\left(\frac{P(LBP_i = k \mid X_i)}{P(LBP_i = 0 \mid X_i)}\right) =\ & \beta_0^{(k)} + \beta_1^{(k)} \cdot \mathrm{PA}_{i,\text{Partially}} + \beta_2^{(k)} \cdot \mathrm{PA}_{i,\text{Irregular}} + \beta_3^{(k)} \cdot \mathrm{PA}_{i,\text{Regular}} + \beta_4^{(k)} \cdot \mathrm{PA}_{i,\text{Trained}} \\ & + \beta_5^{(k)} \cdot \mathrm{Diet}_{i,\text{0--2 portions}} + \beta_6^{(k)} \cdot \mathrm{Diet}_{i,\text{3--4 portions}} + \beta_7^{(k)} \cdot \mathrm{Diet}_{i,\geq5\ \text{portions}} \\ & + \beta_8^{(k)} \cdot \mathrm{Sleep}_{i,\text{moderate}} + \beta_9^{(k)} \cdot \mathrm{Sleep}_{i,\text{pathological}} \\ & + \beta_{10}^{(k)} \cdot \mathrm{SocialRelationships}_{i,\text{medium}} + \beta_{11}^{(k)} \cdot \mathrm{SocialRelationships}_{i,\text{strong}} \\ & + \beta_{12}^{(k)} \cdot \mathrm{SedentaryBehaviour}_i \\ & + \beta_{13}^{(k)} \cdot \mathrm{Tobacco}_{i,\text{occasional}} + \beta_{14}^{(k)} \cdot \mathrm{Tobacco}_{i,\text{daily}} \\ & + \beta_{15}^{(k)} \cdot \mathrm{Alcohol}_{i,\text{low risk}} + \beta_{16}^{(k)} \cdot \mathrm{Alcohol}_{i,\text{minimal risk}} \\ & + \beta_{17}^{(k)} \cdot \mathrm{Alcohol}_{i,\text{medium risk}} + \beta_{18}^{(k)} \cdot \mathrm{Alcohol}_{i,\text{increased risk}} \\ & + \beta_{19}^{(k)} \cdot \mathrm{Drugs}_{i,\text{never used}} + \beta_{20}^{(k)} \cdot \mathrm{Drugs}_{i,>12\text{mo}} \\ & + \beta_{21}^{(k)} \cdot \mathrm{Drugs}_{i,\text{last 12 mo}} + \beta_{22}^{(k)} \cdot \mathrm{Drugs}_{i,\text{last 30 d}} \\ & + \beta_{23}^{(k)} \cdot \mathrm{Depression}_{i,\text{mild}} + \beta_{24}^{(k)} \cdot \mathrm{Depression}_{i,\text{moderate}} \\ & + \beta_{25}^{(k)} \cdot \mathrm{Depression}_{i,\text{rather severe}} + \beta_{26}^{(k)} \cdot \mathrm{Depression}_{i,\text{severe}} \\ & + \beta_{27}^{(k)} \cdot \mathrm{Anxiety}_{i,\text{light}} + \beta_{28}^{(k)} \cdot \mathrm{Anxiety}_{i,\text{moderate}} + \beta_{29}^{(k)} \cdot \mathrm{Anxiety}_{i,\text{severe}} \\ & + \beta_{30}^{(k)} \cdot \mathrm{Hypertension}_{i,\text{yes}} + \beta_{31}^{(k)} \cdot \mathrm{Cholesterol}_{i,\text{yes}} + \beta_{32}^{(k)} \cdot \mathrm{Diabetes}_{i,\text{yes}} \\ & + \beta_{33}^{(k)} \cdot \mathrm{Age}_{i,\text{25--34}} + \beta_{34}^{(k)} \cdot \mathrm{Age}_{i,\text{35--44}} + \beta_{35}^{(k)} \cdot \mathrm{Age}_{i,\text{45--54}} \\ & + \beta_{36}^{(k)} \cdot \mathrm{Age}_{i,\text{55--64}} + \beta_{37}^{(k)} \cdot \mathrm{Age}_{i,\text{65--74}} + \beta_{38}^{(k)} \cdot \mathrm{Age}_{i,\text{75+}} \\ & + \beta_{39}^{(k)} \cdot \mathrm{BMI}_{i,\text{Underweight}} + \beta_{40}^{(k)} \cdot \mathrm{BMI}_{i,\text{Overweight}} + \beta_{41}^{(k)} \cdot \mathrm{BMI}_{i,\text{Obese}} \\ & + \beta_{42}^{(k)} \cdot \mathrm{CivilStatus}_{i,\text{Married}} + \beta_{43}^{(k)} \cdot \mathrm{CivilStatus}_{i,\text{Divorced}} + \beta_{44}^{(k)} \cdot \mathrm{CivilStatus}_{i,\text{Widowed}} \\ & + \beta_{45}^{(k)} \cdot \mathrm{Nationality}_{i,\text{Foreigner}} \\ & + \beta_{46}^{(k)} \cdot \mathrm{Region}_{i,\text{French}} + \beta_{47}^{(k)} \cdot \mathrm{Region}_{i,\text{Italian}} \\ & + \beta_{48}^{(k)} \cdot \mathrm{Education}_{i,\text{UpperSec}} + \beta_{49}^{(k)} \cdot \mathrm{Education}_{i,\text{Tertiary}} \end{aligned} \]

Code 
vglm.q2 <- svy_vglm(
  lbp ~ pa_cat + diet + sleep + soc + sed_beh + tobacco + alcohol + drugs + dep + anxiety + aht + hyperchol + diabetes + sex + age_cat + bmi_cat + civil_status + nat + lang_r + education,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.LBP.bin.c))

m4.q2.raw <- parameters::model_parameters(vglm.q2, 
                                          exponentiate = TRUE,
                                          drop = "Intercept",
                                          digits = 3)

m4.q2.adapted <- adapt.models(m4.q2.raw, variables.m4)
Code 
vglm.m4.q2.male <- svy_vglm(
  lbp ~ pa_cat + diet + sleep + soc + sed_beh + tobacco + alcohol + drugs + dep + anxiety + aht + hyperchol + diabetes + age_cat + bmi_cat + civil_status + nat + lang_r + education, 
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.male))

m4.q2.male <- parameters::model_parameters(vglm.m4.q2.male, 
                                          exponentiate = TRUE,
                                          drop = "Intercept",
                                          digits = 3)

m4.q2.male.adapted <- adapt.models(m4.q2.male, variables.m4)
Code 
vglm.m4.q2.female <- svy_vglm(
  lbp ~ pa_cat + diet + sleep + soc + sed_beh + tobacco + alcohol + drugs + dep + anxiety + aht + hyperchol + diabetes + age_cat + bmi_cat + civil_status + nat + lang_r + education, 
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.female))

m4.q2.female <- parameters::model_parameters(vglm.m4.q2.female, 
                                          exponentiate = TRUE,
                                          drop = "Intercept",
                                          digits = 3)

m4.q2.female.adapted <- adapt.models(m4.q2.female, variables.m4)
Code 
m4.q2.male.adapted <- m4.q2.male.adapted %>% 
  mutate(
    `Adjusted OR (male)` = `Odds Ratio`,
    `95% CI (male)` = `95% CI`,
    `p (male)` = p) %>% 
  select(LBP, Variable, Parameter, `Adjusted OR (male)`, `95% CI (male)`, `p (male)`)

m4.q2.female.adapted <- m4.q2.female.adapted %>% 
  mutate(
    `Adjusted OR (female)` = `Odds Ratio`,
    `95% CI (female)` = `95% CI`,
    `p (female)` = p) %>% 
  select(LBP, Variable, Parameter, `Adjusted OR (female)`, `95% CI (female)`, `p (female)`)

m4.q2.combined <- dplyr::bind_cols(
  m4.q2.male.adapted,
  m4.q2.female.adapted %>%
    select(`Adjusted OR (female)`, `95% CI (female)`, `p (female)`)) %>% 
  mutate(
    Variable = dplyr::recode(Variable,
      pa_cat = "Physical activity level",
      diet = "Diet",
      sleep = "Sleep",
      soc = "Social relationships",
      tobacco = "Tobacco",
      alcohol = "Alcohol",
      drugs = "Drugs",
      dep = "Depression",
      anxiety = "Anxiety",
      aht = "Arterial hypertension",
      sex = "Sex",
      age_cat = "Age category",
      bmi_cat = "BMI category",
      civil_status = "Civil status",
      nat = "Nationality",
      sed_beh = "Sedentary behaviour",
      lang_r = "Language region",
      nat = "Nationality",
      education = "Education",
      hyperchol = "Hypercholesterolaemia",
      diabetes = "Diabetes"))


ft.m4.q2.combined <- flextable(m4.q2.combined)

ft.m4.q2.combined <- ft.m4.q2.combined %>%
  set_caption("Table. Adjusted model on association of lifestyle factors, model by Hartvigsen with LBP by sex") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all")
ft.m4.q2.combined
save(m4.q2.combined, file = "../tables/models/m4.q2.combined.rda")

flextable::flextable(m4.q2.combined) %>%
  set_caption("Table. Adjusted model on association of lifestyle factors, model by Hartvigsen with LBP by sex") %>%
  autofit() %>%
  bold(part = "header") %>%
  align(align = "center", part = "all") %>%
  flextable::line_spacing(space = 1.0) %>%
  flextable::fontsize(size = 8, part = "all") %>%
  flextable::font(fontname = "Times New Roman", part = "all") %>%
  flextable::save_as_docx(path = "../tables/models/model4.q2.by.sex.docx")
Table 10: Model 4: Adjusted model on association of lifestyle factors, model by Hartvigsen with LBP by sex.

LBP

Variable

Parameter

Adjusted OR (male)

95% CI (male)

p (male)

Adjusted OR (female)

95% CI (female)

p (female)

A little LBP

Physical activity level

Partially active

1.14

[0.78; 1.65]

0.492

1.06

[0.79; 1.43]

0.699

Physical activity level

Irregularly active

1.32

[0.93; 1.87]

0.123

1.07

[0.8; 1.42]

0.663

Physical activity level

Regularly active

1.27

[0.88; 1.83]

0.204

0.92

[0.69; 1.24]

0.583

Physical activity level

Trained

0.98

[0.69; 1.4]

0.930

0.82

[0.61; 1.09]

0.176

Diet

0-2 portions

1.08

[0.87; 1.34]

0.468

0.91

[0.69; 1.22]

0.539

Diet

3-4 portions

1.34

[1.06; 1.68]

0.013

0.89

[0.67; 1.18]

0.412

Diet

≥5 portions

1.17

[0.88; 1.55]

0.290

0.88

[0.65; 1.19]

0.417

Sleep

moderate

1.42

[1.2; 1.68]

0.000

1.16

[1; 1.35]

0.043

Sleep

pathological

1.22

[0.88; 1.7]

0.226

1.18

[0.93; 1.49]

0.170

Social relationships

medium

1.03

[0.78; 1.35]

0.851

0.79

[0.62; 1.02]

0.071

Social relationships

strong

0.97

[0.73; 1.28]

0.834

0.74

[0.57; 0.95]

0.018

Sedentary behaviour

0.99

[0.97; 1.01]

0.447

1.01

[0.98; 1.03]

0.595

Tobacco

occasional smoker

0.86

[0.67; 1.11]

0.238

0.85

[0.66; 1.1]

0.217

Tobacco

daily smoker

1.07

[0.87; 1.31]

0.536

1.05

[0.87; 1.27]

0.603

Alcohol

low risk

1.70

[1.32; 2.2]

0.000

1.07

[0.89; 1.28]

0.467

Alcohol

minimal risk

1.88

[1.37; 2.57]

0.000

1.21

[0.96; 1.53]

0.109

Alcohol

medium risk

1.82

[1.12; 2.95]

0.016

0.94

[0.64; 1.39]

0.749

Alcohol

increased risk

2.00

[1.11; 3.6]

0.020

1.54

[0.76; 3.16]

0.234

Drugs

>12 months ago

1.13

[0.96; 1.34]

0.140

1.34

[1.13; 1.58]

0.001

Drugs

in the last 12 months

0.87

[0.6; 1.27]

0.477

1.36

[0.88; 2.1]

0.171

Drugs

in the last 30 days

0.96

[0.65; 1.4]

0.825

1.09

[0.64; 1.86]

0.749

Depression

Mild depression symptoms

1.42

[1.17; 1.73]

0.000

1.25

[1.06; 1.48]

0.007

Depression

Moderate depression symptoms

1.60

[1.07; 2.38]

0.021

1.74

[1.27; 2.39]

0.001

Depression

Rather severe depression symptoms

1.76

[0.95; 3.28]

0.074

1.26

[0.74; 2.14]

0.395

Depression

Severe depression symptoms

2.24

[0.82; 6.11]

0.116

1.34

[0.5; 3.56]

0.562

Anxiety

Light

1.28

[1.03; 1.58]

0.028

1.31

[1.1; 1.56]

0.003

Anxiety

Moderate

1.74

[1.1; 2.74]

0.017

1.48

[1.03; 2.14]

0.035

Anxiety

Severe

0.86

[0.38; 1.91]

0.706

0.83

[0.47; 1.47]

0.529

Arterial hypertension

Yes

1.10

[0.91; 1.32]

0.324

1.17

[0.97; 1.4]

0.106

Hypercholesterolaemia

Yes

0.93

[0.77; 1.12]

0.434

1.00

[0.83; 1.22]

0.971

Diabetes

Yes

1.12

[0.84; 1.49]

0.455

0.84

[0.6; 1.18]

0.313

Age category

18-24 years

0.64

[0.44; 0.92]

0.016

1.10

[0.82; 1.49]

0.527

Age category

25-34 years

1.10

[0.84; 1.45]

0.491

1.10

[0.87; 1.4]

0.433

Age category

35-44 years

0.91

[0.72; 1.14]

0.406

0.92

[0.75; 1.12]

0.401

Age category

55-64 years

0.92

[0.74; 1.13]

0.426

0.92

[0.76; 1.11]

0.381

Age category

65-74 years

1.03

[0.81; 1.3]

0.818

0.87

[0.71; 1.08]

0.217

Age category

75+ years

0.96

[0.7; 1.31]

0.802

0.92

[0.69; 1.24]

0.586

BMI category

Underweight

0.85

[0.44; 1.67]

0.640

0.97

[0.74; 1.28]

0.836

BMI category

Overweight

0.81

[0.69; 0.95]

0.009

1.30

[1.11; 1.51]

0.001

BMI category

Obesity

1.16

[0.92; 1.45]

0.201

1.25

[1.01; 1.55]

0.043

Civil status

Married

0.97

[0.8; 1.18]

0.780

1.17

[0.98; 1.4]

0.080

Civil status

Divorced

1.05

[0.79; 1.4]

0.746

1.27

[1; 1.63]

0.053

Civil status

Widowed

1.29

[0.73; 2.28]

0.373

1.28

[0.91; 1.81]

0.155

Nationality

Foreigner

0.96

[0.8; 1.15]

0.646

1.03

[0.87; 1.22]

0.750

Language region

French-speaking Switzerland

1.06

[0.91; 1.25]

0.461

1.08

[0.93; 1.24]

0.305

Language region

Italian-speaking Switzerland

0.76

[0.59; 0.97]

0.025

0.93

[0.76; 1.15]

0.524

Education

Upper secondary education

1.00

[0.76; 1.32]

0.987

0.92

[0.74; 1.15]

0.490

Education

Tertiary level

0.77

[0.57; 1.03]

0.073

0.85

[0.67; 1.07]

0.165

Strong LBP

Physical activity level

Partially active

0.73

[0.41; 1.32]

0.299

0.81

[0.54; 1.22]

0.316

Physical activity level

Irregularly active

0.66

[0.38; 1.17]

0.155

0.73

[0.49; 1.09]

0.126

Physical activity level

Regularly active

0.77

[0.43; 1.37]

0.373

0.71

[0.47; 1.09]

0.117

Physical activity level

Trained

0.60

[0.33; 1.1]

0.098

0.49

[0.32; 0.75]

0.001

Diet

0-2 portions

1.04

[0.67; 1.6]

0.876

0.85

[0.54; 1.35]

0.487

Diet

3-4 portions

1.12

[0.69; 1.83]

0.649

0.83

[0.52; 1.32]

0.426

Diet

≥5 portions

1.52

[0.84; 2.73]

0.165

0.99

[0.6; 1.62]

0.965

Sleep

moderate

1.46

[1.04; 2.05]

0.027

1.98

[1.54; 2.54]

0.000

Sleep

pathological

1.63

[1.01; 2.63]

0.046

3.53

[2.51; 4.97]

0.000

Social relationships

medium

0.72

[0.44; 1.17]

0.181

0.60

[0.41; 0.87]

0.007

Social relationships

strong

0.66

[0.4; 1.1]

0.110

0.52

[0.35; 0.76]

0.001

Sedentary behaviour

1.01

[0.96; 1.06]

0.832

1.01

[0.97; 1.05]

0.681

Tobacco

occasional smoker

1.89

[1.13; 3.15]

0.014

1.24

[0.78; 1.97]

0.371

Tobacco

daily smoker

2.07

[1.43; 2.99]

0.000

1.39

[1.04; 1.86]

0.024

Alcohol

low risk

0.99

[0.65; 1.51]

0.974

0.89

[0.67; 1.19]

0.440

Alcohol

minimal risk

0.98

[0.53; 1.78]

0.936

0.71

[0.45; 1.12]

0.141

Alcohol

medium risk

0.43

[0.16; 1.15]

0.092

0.70

[0.37; 1.3]

0.257

Alcohol

increased risk

2.17

[0.84; 5.6]

0.109

0.51

[0.12; 2.15]

0.359

Drugs

>12 months ago

1.02

[0.72; 1.44]

0.932

1.37

[1.02; 1.83]

0.034

Drugs

in the last 12 months

1.25

[0.6; 2.6]

0.557

1.60

[0.76; 3.35]

0.212

Drugs

in the last 30 days

1.24

[0.62; 2.51]

0.542

1.28

[0.54; 2.99]

0.575

Depression

Mild depression symptoms

2.26

[1.54; 3.31]

0.000

1.52

[1.12; 2.07]

0.008

Depression

Moderate depression symptoms

4.84

[2.5; 9.35]

0.000

2.25

[1.4; 3.62]

0.001

Depression

Rather severe depression symptoms

4.75

[2.09; 10.8]

0.000

2.77

[1.37; 5.61]

0.005

Depression

Severe depression symptoms

6.06

[1.9; 19.26]

0.002

1.62

[0.51; 5.08]

0.411

Anxiety

Light

1.19

[0.78; 1.82]

0.417

1.52

[1.12; 2.06]

0.007

Anxiety

Moderate

2.70

[1.31; 5.55]

0.007

1.61

[0.97; 2.67]

0.068

Anxiety

Severe

1.10

[0.41; 2.96]

0.846

0.85

[0.39; 1.86]

0.680

Arterial hypertension

Yes

1.48

[1.05; 2.1]

0.026

1.27

[0.92; 1.74]

0.149

Hypercholesterolaemia

Yes

1.23

[0.84; 1.79]

0.280

1.16

[0.83; 1.62]

0.382

Diabetes

Yes

1.23

[0.77; 1.98]

0.383

0.99

[0.56; 1.76]

0.971

Age category

18-24 years

0.33

[0.13; 0.87]

0.025

0.95

[0.54; 1.68]

0.863

Age category

25-34 years

0.90

[0.49; 1.66]

0.730

1.33

[0.85; 2.08]

0.209

Age category

35-44 years

0.77

[0.46; 1.29]

0.319

1.37

[0.94; 1.98]

0.100

Age category

55-64 years

1.09

[0.7; 1.71]

0.704

1.19

[0.83; 1.71]

0.338

Age category

65-74 years

1.30

[0.78; 2.17]

0.308

0.92

[0.61; 1.4]

0.698

Age category

75+ years

0.91

[0.48; 1.73]

0.771

1.24

[0.71; 2.17]

0.459

BMI category

Underweight

1.62

[0.42; 6.19]

0.482

1.01

[0.59; 1.71]

0.979

BMI category

Overweight

0.77

[0.55; 1.07]

0.117

1.21

[0.92; 1.57]

0.169

BMI category

Obesity

1.34

[0.89; 1.99]

0.157

1.62

[1.15; 2.27]

0.005

Civil status

Married

1.98

[1.32; 2.98]

0.001

1.42

[1.03; 1.97]

0.033

Civil status

Divorced

1.57

[0.86; 2.88]

0.145

1.87

[1.21; 2.88]

0.005

Civil status

Widowed

1.32

[0.38; 4.62]

0.660

1.87

[1.08; 3.24]

0.026

Nationality

Foreigner

1.22

[0.83; 1.78]

0.316

1.10

[0.83; 1.45]

0.512

Language region

French-speaking Switzerland

0.95

[0.69; 1.3]

0.737

1.06

[0.83; 1.35]

0.656

Language region

Italian-speaking Switzerland

0.61

[0.35; 1.07]

0.083

0.70

[0.47; 1.03]

0.068

Education

Upper secondary education

2.72

[1.44; 5.13]

0.002

0.84

[0.6; 1.19]

0.325

Education

Tertiary level

1.51

[0.79; 2.89]

0.217

0.53

[0.36; 0.78]

0.001

12.6 Summary table

Code 
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)

clean_model <- function(df, model_name, sex) {
  df %>%
    select(Parameter, Coefficient, CI_low, CI_high) %>%
    rename(OR = Coefficient) %>%
    mutate(Model = model_name, Sex = sex)
  }

# Model 1 (no sex split)
m1_raw <- clean_model(m1.q1.raw, "Model 1 raw", "All")
m1_adj <- clean_model(m1.q1.adj, "Model 1 adj", "All")

# Model 2
m2_male   <- clean_model(m2.q1.male.raw,   "Model 2", "Male")
m2_female <- clean_model(m2.q1.female.raw, "Model 2", "Female")

# Model 3
m3_male   <- clean_model(m3.q2.male,   "Model 3", "Male")
m3_female <- clean_model(m3.q2.female, "Model 3", "Female")

# Model 4
m4_male   <- clean_model(m4.q2.male,   "Model 4", "Male")
m4_female <- clean_model(m4.q2.female, "Model 4", "Female")

# Combine
df.all.models <- bind_rows(
  m1_raw, m1_adj,
  m2_male, m2_female,
  m3_male, m3_female,
  m4_male, m4_female)

# Adapt model
adapt.models.all <- function(df, variables) {
  df %>%
    # Remove intercept rows
    filter(!str_detect(Parameter, "^\\(Intercept\\)")) %>%
    
    # Round numeric values and create 95% CI columns
    mutate(
      `Odds Ratio` = round(OR, 2),
      CI_low = round(CI_low, 2),
      CI_high = round(CI_high, 2),
      `95% CI` = paste0("[", CI_low, "; ", CI_high, "]")  # optional, keeps text version
    ) %>%
    
    # Identify outcome
    mutate(
      LBP = case_when(
        str_detect(Parameter, ":1$") ~ "A little LBP",
        str_detect(Parameter, ":2$") ~ "Strong LBP",
        TRUE ~ NA_character_
      ),
      Parameter = str_replace(Parameter, ":1$", ""),
      Parameter = str_replace(Parameter, ":2$", "")
    ) %>%
    
    # Extract variable group
    mutate(
      Variable = str_extract(Parameter, paste0(variables, collapse = "|")),
      Parameter = str_trim(str_remove(Parameter, paste0("^(", paste(variables, collapse = "|"), ")_?")))
    ) %>%
    
    # Optional: arrange by LBP
    arrange(LBP) %>%
    
    # Keep relevant columns including numeric CI limits
    select(Model, Sex, LBP, Variable, Parameter, `Odds Ratio`, CI_low, CI_high, `95% CI`)
}


# Adapt the model
df.all.models.ad <- adapt.models.all(df.all.models, variables.m4)
df.all.models.ad <- df.all.models.ad %>%
  mutate(Variable = recode(Variable, !!!label.names))

variable_order <- c(
  "PA category", "Diet", "Sleep", "Social relationships", "Sedentary behaviour",
  "Tobacco", "Alcohol", "Drugs", "Depression", "Anxiety",
  "Art. hypertension", "Hypercholesterolaemia", "Diabetes", "Age category", 
  "BMI category", "Civil status", "Nationality", "Language region", "Education"
)

df.all.models.ad <- df.all.models.ad %>%
  mutate(Variable = factor(Variable, levels = variable_order))
Code 
# Custom group order
group_order <- c(
  "PA category", "Diet", "Sleep", "Social relationships", "Sedentary behaviour", 
  "Tobacco", "Alcohol", "Drugs", "Depression", "Anxiety", 
  "Art. hypertension", "Hypercholesterolaemia", "Diabetes", "Sex", 
  "Age category", "BMI category", "Civil status", "Nationality", 
  "Language region", "Education"
)


df.LBP.little <- df.all.models.ad %>% 
  filter(LBP == "A little LBP") %>% 
  mutate(Variable = if_else(is.na(Variable), "Sex", Variable)) %>% 
  mutate(Variable_Parameter = paste(Variable, Parameter, sep = ":")) %>% 
  mutate(Variable_Parameter = factor(Variable_Parameter, levels = unique(Variable_Parameter))) %>% 
  mutate(
    Group = sapply(strsplit(as.character(Variable_Parameter), ":"), `[`, 1),
    Group = factor(Group, levels = rev(group_order)))

df.LBP.little <- df.LBP.little %>%
  arrange(Group, Variable_Parameter) %>%  # Sort within groups alphabetically
  mutate(Variable_Parameter = factor(Variable_Parameter, levels = unique(Variable_Parameter)))


plot.lbp.little <- ggplot(df.LBP.little, 
                          aes(y = Variable_Parameter,
                              x = `Odds Ratio`,
                              xmin = CI_low,
                              xmax = CI_high,
                              color = Model)) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_pointrange(position = position_dodge(width = 0.7)) +
  scale_x_log10() +
  facet_wrap(~Sex) +
  labs(x = "Odds Ratio", y = "",
       title = "a) A little low back pain") +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "none",
    axis.text.y = element_text(hjust = 1),
    plot.title.position = "plot",         # Ensure title aligns to plot area
    plot.title = element_text(hjust = 0)  # hjust = 0 aligns left, 1 aligns right
  )


plot.lbp.little

Code 
df.LBP.strong <- df.all.models.ad %>% 
  filter(LBP == "Strong LBP") %>% 
  mutate(Variable = if_else(is.na(Variable), "Sex", Variable)) %>% 
  mutate(Variable_Parameter = paste(Variable, Parameter, sep = ":")) %>% 
  mutate(Variable_Parameter = factor(Variable_Parameter, levels = unique(Variable_Parameter))) %>% 
  mutate(
    Group = sapply(strsplit(as.character(Variable_Parameter), ":"), `[`, 1),
    Group = factor(Group, levels = rev(group_order)))

df.LBP.strong <- df.LBP.strong %>%
  arrange(Group, Variable_Parameter) %>%  # Sort within groups alphabetically
  mutate(Variable_Parameter = factor(Variable_Parameter, levels = unique(Variable_Parameter)))


plot.lbp.strong <- ggplot(df.LBP.strong, 
                          aes(y = Variable_Parameter,
                              x = `Odds Ratio`,
                              xmin = CI_low,
                              xmax = CI_high,
                              color = Model)) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_pointrange(position = position_dodge(width = 0.7)) +
  scale_x_log10() +
  facet_wrap(~Sex) +
  labs(x = "Odds Ratio", y = "",
       title = "b) Strong low back pain") +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "bottom",
    axis.text.y = element_text(hjust = 1),
    plot.title.position = "plot",         # Ensure title aligns to plot area
    plot.title = element_text(hjust = 0)  # hjust = 0 aligns left, 1 aligns right
  )


plot.lbp.strong

Code 
# combined plot
library(patchwork)
plot.lbp.combined <- plot.lbp.little/plot.lbp.strong +
  plot_layout(heights = c(1, 1))
Code 
ggsave(
  filename = "../figures/plot.models.little.sum.pdf",
  plot = plot.lbp.little,
  width = 10, height = 8, 
  dpi = 1200) # Plot a little LBP

ggsave(
  filename = "../figures/plot.models.strong.sum.pdf",
  plot = plot.lbp.strong,
  width = 10, height = 8, 
  dpi = 1200) # Plot strong LBP

# combined plot
ggsave(filename = "../figures/plot.models.combined.pdf", 
       plot = plot.lbp.combined, 
       width = 10, height = 16, dpi = 1200)

13 Supplementary files

13.1 Selection of package for model analysis

Information

  • nnet: analyses not possible with stratified data.

  • mlogit: analyses not possible with weighted data.

  • survey: analysis only possible with binary variables.

  • svyVGAM: package considers both weights and stratification. THerefore svyVGAM was used for the analyses.

  • No stratification possible
Code 
# weights are considered, stratification is not! - nnet cannot be used
# nnet.model.pa <- nnet::multinom(lbp ~ pa_cat, data = as.data.frame(srvy.LBP.bin.conf.c), weights = weighting.tel)
# par.nnet.model.lbp.pa <- parameters::parameters(nnet.model.pa,  ci = 0.95, digits = 3, exponentiate = TRUE)
# par.nnet.model.lbp.pa

  • weights necessary

  • mlogit(lbp ~ |pa_cat, data, #weights #strata) not used and possible with weighted survey data + stratification

  • Not used

  • svyglm model only applicable for binary outcomes (Only possible for lbp_bin)

  • lbp is not binary

Code 
# survey design

survey.LBP.bin.c <- svydesign(
  weights = ~weighting.tel,
  strata = ~strata,
  ids = ~1,
  nest = TRUE,
  data = df.m.LBP.bin
)

# only possible using svyglm with lbp_bin
survey.model.pa <- svyglm(lbp ~ pa_cat, design = survey.LBP.bin.c, family = quasibinomial())
par.survey.model.lbp.bin.pa <- parameters::parameters(survey.model.pa,  
                                                      ci = 0.95, 
                                                      digits = 3, 
                                                      exponentiate = TRUE) # having no LBP as reference level
summary(survey.model.pa) # it treats lbp as a binary variable - not useful

Call:
svyglm(formula = lbp ~ pa_cat, design = survey.LBP.bin.c, family = quasibinomial())

Survey design:
svydesign(weights = ~weighting.tel, strata = ~strata, ids = ~1, 
    nest = TRUE, data = df.m.LBP.bin)

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               0.07876    0.09176   0.858   0.3907    
pa_catPartially active   -0.14724    0.10607  -1.388   0.1651    
pa_catIrregularly active -0.20206    0.09985  -2.024   0.0430 *  
pa_catRegularly active   -0.23753    0.10486  -2.265   0.0235 *  
pa_catTrained            -0.53877    0.10030  -5.372 7.92e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for quasibinomial family taken to be 1.000075)

Number of Fisher Scoring iterations: 4
Code 
library(see)  # for the plot() method from parameters
plot.survey.model.lbp.bin <-  plot(par.survey.model.lbp.bin.pa, 
     show_labels = TRUE, 
     size_text = 4, 
     show_intercept = TRUE)
# export via ggsave
ggsave("../figures/plot.model.lbp.bin.png", plot = plot.survey.model.lbp.bin, width = 6, height = 4, dpi = 300)
  • Not used
Code 
# It uses survey design object from survey package
library(svyVGAM)
library(parameters)

vglm.m1 <- svy_vglm(
  lbp ~ pa_cat,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.LBP.bin.c)
)
summary(vglm.m1)
svy_vglm.survey.design(lbp ~ pa_cat, family = multinomial(refLevel = "No"), 
    design = as_survey_design(srvy.LBP.bin.c))
Stratified Independent Sampling design (with replacement)
Called via srvyr
Data variables: 
  - id (int), age_cat (fct), bmi_cat (fct), civil_status (fct), lang_r (fct),
    nat (fct), sex (fct), education (fct), strata (fct), lbp (fct), pa_cat
    (fct), diet (fct), sleep (fct), sed_beh (dbl), soc (fct), tobacco (fct),
    alcohol (fct), drugs (fct), dep (fct), anxiety (fct), aht (fct), hyperchol
    (fct), diabetes (fct), weighting.tel (dbl), lbp_bin (fct)
                                 Coef         SE       z         p
(Intercept):1              -0.2694763  0.1029545 -2.6174  0.008859
(Intercept):2              -1.1451685  0.1258706 -9.0980 < 2.2e-16
pa_catPartially active:1    0.0090632  0.1171684  0.0774  0.938344
pa_catPartially active:2   -0.6683198  0.1588354 -4.2076 2.581e-05
pa_catIrregularly active:1 -0.0065476  0.1108548 -0.0591  0.952901
pa_catIrregularly active:2 -0.9325987  0.1491919 -6.2510 4.078e-10
pa_catRegularly active:1   -0.0715576  0.1161328 -0.6162  0.537782
pa_catRegularly active:2   -0.8056273  0.1559499 -5.1659 2.392e-07
pa_catTrained:1            -0.3278969  0.1112864 -2.9464  0.003215
pa_catTrained:2            -1.3678662  0.1568012 -8.7236 < 2.2e-16
Code 
# including confounder age and bmi (psych_stress missing)
vgl.m2 <- svy_vglm(
  lbp ~ pa_cat + age_cat + bmi_cat,
  family = multinomial(refLevel = "No"),
  design = as_survey_design(srvy.LBP.bin.c)
)
# reference categories: PA: inactive, LBP: no LBP, x1 = A little LBP, x2 = Strong LBP
# parameters::parameters(vglm.model.pa, exponentiate = TRUE,  digits = 3)
m1 <- model_parameters(vglm.m1, exponentiate = TRUE)
m1
Parameter                       | Odds Ratio |   SE |       95% CI |     z |      p
-----------------------------------------------------------------------------------
(Intercept) × 1                 |       0.76 | 0.08 | [0.62, 0.93] | -2.62 | 0.009 
(Intercept) × 2                 |       0.32 | 0.04 | [0.25, 0.41] | -9.10 | < .001
pa cat [Partially active] × 1   |       1.01 | 0.12 | [0.80, 1.27] |  0.08 | 0.938 
pa cat [Partially active] × 2   |       0.51 | 0.08 | [0.38, 0.70] | -4.21 | < .001
pa cat [Irregularly active] × 1 |       0.99 | 0.11 | [0.80, 1.23] | -0.06 | 0.953 
pa cat [Irregularly active] × 2 |       0.39 | 0.06 | [0.29, 0.53] | -6.25 | < .001
pa cat [Regularly active] × 1   |       0.93 | 0.11 | [0.74, 1.17] | -0.62 | 0.538 
pa cat [Regularly active] × 2   |       0.45 | 0.07 | [0.33, 0.61] | -5.17 | < .001
pa cat [Trained] × 1            |       0.72 | 0.08 | [0.58, 0.90] | -2.95 | 0.003 
pa cat [Trained] × 2            |       0.25 | 0.04 | [0.19, 0.35] | -8.72 | < .001
Code 
m2 <- model_parameters(vgl.m2, exponentiate = TRUE)
m2
Parameter                       | Odds Ratio |   SE |       95% CI |     z |      p
-----------------------------------------------------------------------------------
(Intercept) × 1                 |       0.77 | 0.09 | [0.61, 0.97] | -2.24 | 0.025 
(Intercept) × 2                 |       0.28 | 0.05 | [0.21, 0.39] | -7.69 | < .001
pa cat [Partially active] × 1   |       1.02 | 0.12 | [0.81, 1.28] |  0.14 | 0.890 
pa cat [Partially active] × 2   |       0.54 | 0.09 | [0.39, 0.73] | -3.93 | < .001
pa cat [Irregularly active] × 1 |       1.00 | 0.11 | [0.81, 1.25] |  0.04 | 0.971 
pa cat [Irregularly active] × 2 |       0.41 | 0.06 | [0.31, 0.55] | -5.95 | < .001
pa cat [Regularly active] × 1   |       0.94 | 0.11 | [0.75, 1.18] | -0.54 | 0.589 
pa cat [Regularly active] × 2   |       0.46 | 0.07 | [0.34, 0.63] | -4.94 | < .001
pa cat [Trained] × 1            |       0.73 | 0.08 | [0.59, 0.91] | -2.80 | 0.005 
pa cat [Trained] × 2            |       0.29 | 0.05 | [0.21, 0.39] | -7.92 | < .001
age cat [18-24 years] × 1       |       0.94 | 0.09 | [0.77, 1.14] | -0.63 | 0.527 
age cat [18-24 years] × 2       |       0.60 | 0.12 | [0.40, 0.88] | -2.59 | 0.010 
age cat [25-34 years] × 1       |       1.09 | 0.09 | [0.93, 1.28] |  1.07 | 0.283 
age cat [25-34 years] × 2       |       0.88 | 0.14 | [0.65, 1.21] | -0.78 | 0.433 
age cat [35-44 years] × 1       |       0.91 | 0.07 | [0.78, 1.05] | -1.29 | 0.197 
age cat [35-44 years] × 2       |       0.97 | 0.14 | [0.73, 1.28] | -0.23 | 0.815 
age cat [55-64 years] × 1       |       0.93 | 0.07 | [0.81, 1.07] | -0.98 | 0.325 
age cat [55-64 years] × 2       |       1.18 | 0.16 | [0.91, 1.53] |  1.26 | 0.209 
age cat [65-74 years] × 1       |       0.95 | 0.07 | [0.82, 1.09] | -0.73 | 0.463 
age cat [65-74 years] × 2       |       0.98 | 0.14 | [0.75, 1.29] | -0.13 | 0.900 
age cat [75+ years] × 1         |       0.97 | 0.09 | [0.80, 1.18] | -0.29 | 0.773 
age cat [75+ years] × 2         |       1.09 | 0.20 | [0.76, 1.56] |  0.45 | 0.655 
bmi cat [Underweight] × 1       |       1.06 | 0.13 | [0.83, 1.36] |  0.50 | 0.618 
bmi cat [Underweight] × 2       |       1.36 | 0.32 | [0.86, 2.16] |  1.31 | 0.191 
bmi cat [Overweight] × 1        |       0.95 | 0.05 | [0.86, 1.05] | -0.96 | 0.339 
bmi cat [Overweight] × 2        |       0.96 | 0.10 | [0.79, 1.17] | -0.41 | 0.679 
bmi cat [Obesity] × 1           |       1.21 | 0.09 | [1.05, 1.40] |  2.60 | 0.009 
bmi cat [Obesity] × 2           |       1.77 | 0.21 | [1.40, 2.24] |  4.72 | < .001

14 Package Versions

Code 
sessionInfo()
R version 4.5.0 (2025-04-11)
Platform: x86_64-apple-darwin20
Running under: macOS Sequoia 15.7.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Zurich
tzcode source: internal

attached base packages:
 [1] splines   stats4    grid      stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] see_0.11.0         patchwork_1.3.0    flextable_0.9.9    svyVGAM_1.2       
 [5] VGAM_1.1-13        latex2exp_0.9.6    table1_1.4.3       VIM_6.2.2         
 [9] colorspace_2.1-1   emmeans_1.11.1     mice_3.18.0        survey_4.4-2      
[13] survival_3.8-3     Matrix_1.7-3       srvyr_1.3.0        gtsummary_2.2.0   
[17] nnet_7.3-20        naniar_1.1.0       DataExplorer_0.8.3 kableExtra_1.4.0  
[21] parameters_0.26.0  epitools_0.5-10.1  dagitty_0.3-4      ggdag_0.2.13      
[25] lubridate_1.9.4    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
[29] purrr_1.0.4        readr_2.1.5        tidyr_1.3.1        tibble_3.2.1      
[33] ggplot2_3.5.2      tidyverse_2.0.0    pacman_0.5.1      

loaded via a namespace (and not attached):
  [1] R.oo_1.27.1             cellranger_1.1.0        datawizard_1.1.0       
  [4] rpart_4.1.24            lifecycle_1.0.4         Rdpack_2.6.4           
  [7] lattice_0.22-7          MASS_7.3-65             insight_1.3.0          
 [10] backports_1.5.0         magrittr_2.0.3          vcd_1.4-13             
 [13] rmarkdown_2.29          yaml_2.3.10             zip_2.3.3              
 [16] askpass_1.2.1           sp_2.2-0                DBI_1.2.3              
 [19] minqa_1.2.8             RColorBrewer_1.1-3      multcomp_1.4-28        
 [22] abind_1.4-8             R.utils_2.13.0          TH.data_1.1-3          
 [25] sandwich_3.1-1          gdtools_0.4.2           data.tree_1.1.0        
 [28] cards_0.6.0             cardx_0.2.4             svglite_2.2.1          
 [31] codetools_0.2-20        xml2_1.3.8              tidyselect_1.2.1       
 [34] shape_1.4.6.1           farver_2.1.2            lme4_1.1-37            
 [37] jsonlite_2.0.0          e1071_1.7-16            mitml_0.4-5            
 [40] tidygraph_1.3.1         Formula_1.2-5           iterators_1.0.14       
 [43] systemfonts_1.2.3       foreach_1.5.2           tools_4.5.0            
 [46] ragg_1.4.0              rio_1.2.3               Rcpp_1.0.14            
 [49] glue_1.8.0              gridExtra_2.3           pan_1.9                
 [52] laeken_0.5.3            xfun_0.52               ranger_0.17.0          
 [55] withr_3.0.2             fastmap_1.2.0           mitools_2.4            
 [58] boot_1.3-31             openssl_2.3.3           digest_0.6.37          
 [61] timechange_0.3.0        R6_2.6.1                estimability_1.5.1     
 [64] visdat_0.6.0            textshaping_1.0.1       networkD3_0.4.1        
 [67] R.methodsS3_1.8.2       utf8_1.2.5              generics_0.1.4         
 [70] fontLiberation_0.1.0    data.table_1.17.4       robustbase_0.99-4-1    
 [73] class_7.3-23            htmlwidgets_1.6.4       pkgconfig_2.0.3        
 [76] gtable_0.3.6            lmtest_0.9-40           htmltools_0.5.8.1      
 [79] fontBitstreamVera_0.1.1 carData_3.0-5           scales_1.4.0           
 [82] reformulas_0.4.1        knitr_1.50              rstudioapi_0.17.1      
 [85] tzdb_0.5.0              reshape2_1.4.4          uuid_1.2-1             
 [88] coda_0.19-4.1           nlme_3.1-168            curl_6.2.3             
 [91] nloptr_2.2.1            proxy_0.4-27            zoo_1.8-14             
 [94] parallel_4.5.0          pillar_1.10.2           vctrs_0.6.5            
 [97] car_3.1-3               jomo_2.7-6              xtable_1.8-4           
[100] evaluate_1.0.3          mvtnorm_1.3-3           cli_3.6.5              
[103] compiler_4.5.0          rlang_1.1.6             labeling_0.4.3         
[106] plyr_1.8.9              stringi_1.8.7           viridisLite_0.4.2      
[109] glmnet_4.1-9            bayestestR_0.16.0       V8_6.0.4               
[112] fontquiver_0.2.1        hms_1.1.3               rbibutils_2.3          
[115] igraph_2.1.4            broom_1.0.8             DEoptimR_1.1-3-1       
[118] readxl_1.4.5            officer_0.6.10