library(tidyverse)
library(haven)
library(survey)
library(gt)
## Data Wrangling - BRFSS Weight Calculations
brfss_23 <- read_sav("~/Downloads/GU23FINLvSplit.SAV")
brfss_24 <- read_sav("~/Downloads/GU24FINLvSplit.SAV")
## Calculate sample sizes
n23 <- nrow(brfss_23)
n24 <- nrow(brfss_24)
total_n <- n23 + n24
## Compute adjusted weights that will be used with the combined dataset (brfss_23+brfss_24)
brfss_23 <- brfss_23 %>%
mutate(combine_final_weight = v74 * (n23 / total_n))
brfss_24 <- brfss_24 %>%
mutate(combine_final_weight = v74 * (n24 / total_n))We adjusted the BRFSS sampling weights for the pooled 2023-2024 BRFSS guam dataset by taking the product each respondent’s original survey weight by the proportion of the combined sample contributed by their respective survey year.
This is the CDC link in which we got the approach:
Re-coding of variables from BRFSS 2023
## Recoding variables
brfss_23 <- brfss_23 %>%
mutate(
# Ethnicity
ethnicity = case_when(
Ethnic1 %in% c(1) ~ "CHamoru",
Ethnic1 %in% 9:14 ~ "Other Micronesian",
Ethnic1 == 2 ~ "Filipino",
TRUE ~ "Other"
),
# Ethnicity for Regression only
ethnic_regression = case_when(
Ethnic1 %in% c(1,9,10,11,12) ~ "Pacific",
Ethnic1 %in% c(2,4,5,6,7) ~ "Asian",
TRUE ~ "Other"
),
# Sex
sex = case_when(
v145 == 1 ~ "Male",
v145 == 2 ~ "Female",
TRUE ~ NA_character_
),
# Age categories
age_cat = case_when(
v146 %in% 1:1 ~ "18 to 24",
v146 %in% 2:3 ~ "25 to 34",
v146 %in% 4:5 ~ "35 to 44",
v146 %in% 6:7 ~ "45 to 54",
v146 %in% 8:9 ~ "55 to 64",
v146 >= 10 ~ "65 or older",
TRUE ~ NA_character_
),
# Numeric age midpoint
age_num = case_when(
AGE %in% c(8,9) ~ NA_integer_,
TRUE ~ as.numeric(AGE)
),
# Education
edu_status = case_when(
EDUCA %in% 1:4 ~ "High School or Below",
EDUCA %in% 5:6 ~ "College",
EDUCA == 9 ~ NA_character_,
TRUE ~ NA_character_
),
# Gout
gout_status = case_when(
is.na(GoutDiag) ~ NA_character_,
GoutDiag == 97 ~ NA_character_,
GoutDiag == 96 ~ "No",
TRUE ~ "Yes"
),
early_gout_onset = case_when(
is.na(GoutDiag) ~ NA_character_,
GoutDiag == 97 ~ NA_character_,
GoutDiag == 96 ~ "No",
GoutDiag < 40 ~ "Yes",
TRUE ~ NA_character_
),
# Diabetes
diabetes = case_when(
DIABETE4 %in% c(1,2,4) ~ "Yes",
DIABETE4 == 2 ~ "No",
DIABETE4 %in% c(7,9) ~ NA_character_,
TRUE ~ NA_character_
),
# Kidney disease
kidney_disease = case_when(
CHCKDNY2 == 1 ~ "Yes",
CHCKDNY2 == 2 ~ "No",
CHCKDNY2 == 7 ~ NA_character_,
TRUE ~ NA_character_
),
# BMI scaled
bmi = ifelse(is.na(v150), NA_real_, v150 / 100),
# Alcohol
alcohol = case_when(
ALCDAY4 %in% c(0, 888) ~ "No",
ALCDAY4 >= 1 ~ "Yes",
ALCDAY4 %in% c(77,777,99,999) ~ NA_character_,
TRUE ~ NA_character_
),
# Current smoker
smoke = case_when(
v157 == 1 ~ "No",
v157 == 2 ~ "Yes",
v157 == 9 ~ NA_character_,
TRUE ~ NA_character_
),
# Smokeless tobacco (USENOW3)
smokeless_tobacco = case_when(
USENOW3 %in% c(2) ~ "Yes",
USENOW3 == 3 ~ "No",
USENOW3 == 9 ~ NA_character_,
TRUE ~ NA_character_
),
# Betel nut ever
betel_nut = case_when(
BetChw1 == 1 ~ "Yes",
BetChw1 == 2 ~ "No",
BetChw1 %in% c(7,9) ~ "No",
TRUE ~ "No"
),
# Betel nut current
betel_nut_current = case_when(
BetChw2 %in% c(1, 2) ~ "Yes",
BetChw2 %in% c(3,7,9) ~ "No",
TRUE ~ "No"
)
) %>%
select(
stratum = v2,
sampling_unit = v4,
final_weight = v74,
combine_final_weight,
ethnicity,
sex,
age_cat,
age_num,
edu_status,
GoutDiag,
gout_status,
early_gout_onset,
diabetes,
kidney_disease,
bmi,
alcohol,
smoke,
smokeless_tobacco,
betel_nut,
betel_nut_current,
ethnic_regression
)For most variables, missing data was excluded from population estimates/ preavalences except for areca nut. Including missing data as “No” responses was adopted to match what is done in the literature for this specific behavior as seen in Paulino et al. 2017.
Re-coding of variables from BRFSS 2024
brfss_24 <- brfss_24 %>%
mutate(
# Ethnicity
ethnicity = case_when(
Ethnic1 %in% c(1) ~ "CHamoru",
Ethnic1 %in% 9:14 ~ "Other Micronesian",
Ethnic1 == 2 ~ "Filipino",
TRUE ~ "Other"
),
# Ethnicity for regression only
ethnic_regression = case_when(
Ethnic1 %in% c(1, 9, 10, 11, 12) ~ "Pacific",
Ethnic1 %in% c(2, 4, 5, 6, 7) ~ "Asian",
TRUE ~ "Other"
),
# Sex
sex = case_when(
SEXVAR == 1 ~ "Male",
SEXVAR == 2 ~ "Female",
TRUE ~ NA_character_
),
# Numeric age
age_num = case_when(
AGE %in% c(7, 9) ~ NA_real_,
TRUE ~ as.numeric(AGE)
),
# Age categories
age_cat = case_when(
age_num >= 18 & age_num <= 24 ~ "18 to 24",
age_num >= 25 & age_num <= 34 ~ "25 to 34",
age_num >= 35 & age_num <= 44 ~ "35 to 44",
age_num >= 45 & age_num <= 54 ~ "45 to 54",
age_num >= 55 & age_num <= 64 ~ "55 to 64",
age_num >= 65 ~ "65 or older",
TRUE ~ NA_character_
),
# Education
edu_status = case_when(
EDUCA %in% 1:4 ~ "High School or Below",
EDUCA %in% 5:6 ~ "College",
TRUE ~ NA_character_
),
# Gout
gout_status = case_when(
is.na(GoutDiag) ~ NA_character_,
GoutDiag == 97 ~ NA_character_,
GoutDiag == 96 ~ "No",
TRUE ~ "Yes"
),
# Early-onset gout (diagnosis before age 40)
early_gout_onset = case_when(
is.na(GoutDiag) ~ NA_character_,
GoutDiag == 97 ~ NA_character_,
GoutDiag == 96 ~ "No",
GoutDiag < 40 ~ "Yes",
TRUE ~ NA_character_
),
# Diabetes
diabetes = case_when(
DIABETE4 %in% c(1, 2,4) ~ "Yes",
DIABETE4 == 3 ~ "No",
DIABETE4 %in% c(7) ~ NA_character_,
TRUE ~ NA_character_
),
# Kidney disease
kidney_disease = case_when(
CHCKDNY2 == 1 ~ "Yes",
CHCKDNY2 == 2 ~ "No",
CHCKDNY2 %in% c(7, 9) ~ NA_character_,
TRUE ~ NA_character_
),
# BMI scaled
bmi = ifelse(
is.na(v140),
NA_real_,
as.numeric(v140) / 100
),
# Alcohol
alcohol = case_when(
ALCDAY4 %in% c(0, 888) ~ "No",
ALCDAY4 %in% c(77, 777, 99, 999) ~ NA_character_,
ALCDAY4 >= 1 & ALCDAY4 <= 231 ~ "Yes",
TRUE ~ NA_character_
),
# Current smoker
smoke = case_when(
v163 == 1 ~ "No",
v163 == 2 ~ "Yes",
v163 == 9 ~ NA_character_,
TRUE ~ NA_character_
),
# Smokeless tobacco (USENOW3)
smokeless_tobacco = case_when(
USENOW3 %in% c(2) ~ "Yes",
USENOW3 == 3 ~ "No",
TRUE ~ NA_character_
),
# Betel nut ever
betel_nut = case_when(
BetChw1 == 1 ~ "Yes",
BetChw1 %in% c(2,7,9) ~ "No",
BetChw1 %in% c(7, 9) ~ "No",
TRUE ~ "No"
),
# Betel nut current
betel_nut_current = case_when(
BetChw2 %in% c(1, 2) ~ "Yes",
BetChw2 %in% c(3,7,9) ~ "No",
TRUE ~ "No"
)
) %>%
select(
stratum = v2,
sampling_unit = v4,
final_weight = v74,
combine_final_weight,
ethnicity,
sex,
age_cat,
age_num,
edu_status,
GoutDiag,
gout_status,
early_gout_onset,
diabetes,
kidney_disease,
bmi,
alcohol,
smoke,
smokeless_tobacco,
betel_nut,
betel_nut_current,
ethnic_regression
)Similarly described above - for most variables, missing data was excluded from population estimates/ preavalences except for areca nut.
Merge the two BRFSS datasets
# ============================================================
# MERGE BRFSS 2023 & 2024 DATA
# ============================================================
combined_years <- rbind(brfss_23, brfss_24)
# ============================================================
# CONVERT VARIABLES TO FACTORS
# ============================================================
combined_years <- combined_years %>%
mutate(
age_cat = factor(
age_cat,
levels = c(
"18 to 24",
"25 to 34",
"35 to 44",
"45 to 54",
"55 to 64",
"65 or older"
)
),
ethnicity = factor(
ethnicity,
levels = c(
"CHamoru",
"Filipino",
"Other Micronesian",
"Other"
)
),
ethnic_regression = factor(
ethnic_regression,
levels = c("Pacific", "Asian", "Other")
),
sex = factor(
sex,
levels = c("Male", "Female")
),
edu_status = factor(
edu_status,
levels = c("High School or Below", "College")
),
gout_status = factor(
gout_status,
levels = c("No", "Yes")
),
early_gout_onset = factor(
early_gout_onset,
levels = c("No", "Yes")
),
kidney_disease = factor(
kidney_disease,
levels = c("No", "Yes")
),
diabetes = factor(
diabetes,
levels = c("No", "Yes")
),
alcohol = factor(
alcohol,
levels = c("No", "Yes")
),
smoke = factor(
smoke,
levels = c("No", "Yes")
),
# FIX: Use smokeless_tobacco, NOT smoke
smokeless_tobacco = factor(
smokeless_tobacco,
levels = c("No", "Yes")
),
betel_nut = factor(
betel_nut,
levels = c("No", "Yes")
),
betel_nut_current = factor(
betel_nut_current,
levels = c("No", "Yes")
),
)
table(brfss_23$betel_nut_current)
table(brfss_24$betel_nut_current)
table(combined_years$betel_nut_current)The analysis used the pooled Guam BRFSS 2023–2024 dataset. Descriptive estimates were calculated with the survey design variables and adjusted sampling weights defined above. Unweighted respondent counts are reported alongside survey-weighted percentages and 95% confidence intervals. Regression analyses were conducted using survey-weighted logistic regression with a quasibinomial family to account for the complex survey design.
Formatting of the subsequent tables may not entirely match what is seen in the published manuscript - but all estimates and values remain the same.
library(dplyr)
library(survey)
library(gt)
options(survey.lonely.psu = "adjust")
# ============================================================
# TABLE 1: DEMOGRAPHIC DESCRIPTION
# Guam BRFSS 2023–2024
# ============================================================
# Prepare data
table_1 <- combined_years %>%
select(
combine_final_weight,
stratum,
sampling_unit,
sex,
age_cat,
ethnicity,
edu_status
) %>%
mutate(
sex = factor(
sex,
levels = c("Female", "Male")
),
age_cat = factor(
age_cat,
levels = c(
"18 to 24",
"25 to 34",
"35 to 44",
"45 to 54",
"55 to 64",
"65 or older"
)
),
ethnicity = factor(
ethnicity,
levels = c(
"CHamoru",
"Filipino",
"Other Micronesian",
"Other"
)
),
edu_status = factor(
edu_status,
levels = c(
"College",
"High School or Below"
)
)
)
# Create survey design
table1_survey <- svydesign(
ids = ~sampling_unit,
strata = ~stratum,
weights = ~combine_final_weight,
data = table_1,
nest = TRUE
)
# ============================================================
# FUNCTION: UNWEIGHTED N + WEIGHTED % AND 95% CI
# ============================================================
summarize_variable <- function(var_name, design) {
values <- design$variables[[var_name]]
categories <- levels(values)
valid <- !is.na(values)
design_valid <- design[valid, ]
bind_rows(lapply(categories, function(category) {
indicator <- as.numeric(
as.character(
design_valid$variables[[var_name]]
) == category
)
category_design <- update(
design_valid,
category_indicator = indicator
)
estimate <- svyciprop(
~category_indicator,
design = category_design,
method = "logit",
na.rm = TRUE
)
ci <- confint(estimate)
data.frame(
Variable = category,
Frequency = sum(
values == category,
na.rm = TRUE
),
Proportion_CI = sprintf(
"%.1f%% (%.1f%% – %.1f%%)",
100 * as.numeric(coef(estimate)),
100 * ci[1],
100 * ci[2]
)
)
}))
}
# ============================================================
# GENERATE TABLE
# ============================================================
table_1_table <- bind_rows(
summarize_variable("sex", table1_survey) %>%
mutate(Group = "Sex"),
summarize_variable("age_cat", table1_survey) %>%
mutate(Group = "Age Category"),
summarize_variable("ethnicity", table1_survey) %>%
mutate(Group = "Ethnicity"),
summarize_variable("edu_status", table1_survey) %>%
mutate(Group = "Education Level")
) %>%
mutate(
Group = factor(
Group,
levels = c(
"Sex",
"Age Category",
"Ethnicity",
"Education Level"
)
)
) %>%
select(
Group,
Variable,
Frequency,
Proportion_CI
)
# ============================================================
# FORMAT TABLE USING GT
# ============================================================
table_1_table_object <- table_1_table %>%
gt(groupname_col = "Group") %>%
tab_header(
title = "Table 1. Demographic Description",
subtitle = "Guam BRFSS 2023–2024"
) %>%
cols_label(
Variable = "Variable",
Frequency = paste0(
"Frequency (N = ",
format(
nrow(combined_years),
big.mark = ","
),
")"
),
Proportion_CI = "Proportion (95% CI)"
) %>%
fmt_number(
columns = Frequency,
decimals = 0,
use_seps = TRUE
) %>%
cols_align(
align = "left",
columns = Variable
) %>%
cols_align(
align = "center",
columns = c(Frequency, Proportion_CI)
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_row_groups()
) %>%
tab_source_note(
source_note = paste(
"Frequencies are unweighted.",
"Proportions and 95% confidence intervals",
"are survey-weighted.",
"Missing responses are excluded",
"from variable-specific denominators."
)
)
# Display table
# The formatted table is displayed in the final Tables section.
invisible(table_1_table_object)Formatting of the subsequent tables may not entirely match what is seen in the published manuscript - but all estimates and values remain the same.
# ============================================================
# TABLE 2: PREVALENCE OF SELECTED HEALTH CONDITIONS
# AND BEHAVIORS FROM GUAM BRFSS 2023–2024
# ============================================================
library(dplyr)
library(survey)
library(gt)
options(survey.lonely.psu = "adjust")
# ------------------------------------------------------------
# 1. PREPARE VARIABLES
# ------------------------------------------------------------
table_2 <- combined_years %>%
mutate(
# Correct early-onset gout definition
early_gout_onset = case_when(
GoutDiag == 96 ~ "No",
GoutDiag >= 17 & GoutDiag < 40 ~ "Yes",
GoutDiag >= 40 & GoutDiag <= 95 ~ "No",
TRUE ~ NA_character_
),
# Obesity: BMI >= 30 kg/m2
obesity = case_when(
is.na(bmi) ~ NA_character_,
bmi >= 30 ~ "Yes",
bmi < 30 ~ "No",
TRUE ~ NA_character_
),
# Standardize binary variables
across(
c(
gout_status,
early_gout_onset,
obesity,
diabetes,
kidney_disease,
alcohol,
smoke,
smokeless_tobacco,
betel_nut,
betel_nut_current
),
~factor(.x, levels = c("No", "Yes"))
)
) %>%
select(
combine_final_weight,
stratum,
sampling_unit,
gout_status,
early_gout_onset,
obesity,
diabetes,
kidney_disease,
alcohol,
smoke,
smokeless_tobacco,
betel_nut,
betel_nut_current
)
# ------------------------------------------------------------
# 2. CREATE SURVEY DESIGN
# ------------------------------------------------------------
table2_survey <- svydesign(
ids = ~sampling_unit,
strata = ~stratum,
weights = ~combine_final_weight,
data = table_2,
nest = TRUE
)
# ------------------------------------------------------------
# 3. FUNCTION: FREQUENCY, WEIGHTED % AND 95% CI
# ------------------------------------------------------------
summarize_condition <- function(
var_name,
label,
survey_design
) {
# Extract variable
values <- survey_design$variables[[var_name]]
# Exclude missing responses from denominator
valid <- !is.na(values)
design_valid <- survey_design[valid, ]
# Number of respondents with valid responses
valid_n <- sum(valid)
# Unweighted number reporting Yes
n_yes <- sum(values == "Yes", na.rm = TRUE)
# Create indicator for Yes
indicator <- as.numeric(
as.character(
design_valid$variables[[var_name]]
) == "Yes"
)
# Add indicator to survey design
design_valid <- update(
design_valid,
condition_yes = indicator
)
# Survey-weighted prevalence and CI
estimate <- svyciprop(
~condition_yes,
design = design_valid,
method = "logit",
na.rm = TRUE
)
ci <- confint(estimate)
data.frame(
Variable = label,
Frequency = n_yes,
Valid_N = valid_n,
Proportion = 100 * as.numeric(coef(estimate)),
CI_Lower = 100 * ci[1],
CI_Upper = 100 * ci[2]
)
}
# ------------------------------------------------------------
# 4. GENERATE HEALTH CONDITION SUMMARIES
# ------------------------------------------------------------
health_conditions <- bind_rows(
summarize_condition(
"gout_status",
"Gout",
table2_survey
),
summarize_condition(
"early_gout_onset",
"Early-onset gout",
table2_survey
),
summarize_condition(
"obesity",
"Obesity",
table2_survey
),
summarize_condition(
"diabetes",
"Diabetes",
table2_survey
),
summarize_condition(
"kidney_disease",
"Kidney disease",
table2_survey
)
)
# ------------------------------------------------------------
# 5. GENERATE HEALTH BEHAVIOR SUMMARIES
# ------------------------------------------------------------
health_behaviors <- bind_rows(
summarize_condition(
"alcohol",
"Alcohol consumption",
table2_survey
),
summarize_condition(
"smoke",
"Smoking",
table2_survey
),
summarize_condition(
"smokeless_tobacco",
"Smokeless tobacco use",
table2_survey
),
summarize_condition(
"betel_nut",
"Areca nut use (ever used)",
table2_survey
),
summarize_condition(
"betel_nut_current",
"Areca nut use (currently use)",
table2_survey
)
)
# ------------------------------------------------------------
# 6. COMBINE TABLE 2 RESULTS
# ------------------------------------------------------------
table_2_table <- bind_rows(
health_conditions %>%
mutate(Group = "Health Conditions"),
health_behaviors %>%
mutate(Group = "Health Behaviors")
) %>%
mutate(
Group = factor(
Group,
levels = c(
"Health Conditions",
"Health Behaviors"
)
),
percent_ci = sprintf(
"%.1f%% (%.1f%% – %.1f%%)",
Proportion,
CI_Lower,
CI_Upper
)
) %>%
select(
Group,
Variable,
Frequency,
Valid_N,
percent_ci
)
# ------------------------------------------------------------
# 7. FORMAT MANUSCRIPT TABLE WITH GT
# ------------------------------------------------------------
table_2_table_object <- table_2_table %>%
select(
Group,
Variable,
Frequency,
percent_ci
) %>%
gt(groupname_col = "Group") %>%
tab_header(
title = paste(
"Table 2. Prevalence of Selected Health",
"Conditions and Behaviors"
),
subtitle = "Guam BRFSS 2023–2024"
) %>%
cols_label(
Variable = "Variable",
Frequency = paste0(
"Frequency (N = ",
format(
nrow(combined_years),
big.mark = ","
),
")"
),
percent_ci = "Proportion (95% CI)"
) %>%
fmt_number(
columns = Frequency,
decimals = 0,
use_seps = TRUE
) %>%
cols_align(
align = "left",
columns = Variable
) %>%
cols_align(
align = "center",
columns = c(
Frequency,
percent_ci
)
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_row_groups()
) %>%
tab_source_note(
source_note = paste(
"1. Proportions are survey-weighted using",
"pooled Guam BRFSS 2023 and 2024 weights.",
"2. Parentheses indicate 95% confidence intervals.",
"3. Frequencies are unweighted.",
"Missing responses are excluded from",
"variable-specific denominators.",
"4. Early-onset gout is defined as gout",
"diagnosed before age 40."
)
)
# ------------------------------------------------------------
# 8. DISPLAY TABLE
# ------------------------------------------------------------
# The formatted table is displayed in the final Tables section.
invisible(table_2_table_object)Four exploratory survey-weighted logistic regression models were specified to examine associations between areca nut use and gout outcomes. The models were organized by exposure (ever use versus current use) and outcome (early-onset gout versus gout overall). The reported coefficients are exponentiated to obtain odds ratios (ORs) and are presented with 95% confidence intervals and P values. These cross-sectional associations do not establish causality.
The regression dataset contains the gout outcomes, areca nut exposure variables, respondent age, sex, and ethnicity, along with survey strata, sampling units, and adjusted sampling weights. Binary outcomes and exposures were coded with No as the reference category. Male was the reference category for sex, and Pacific Islander was the reference category for the ethnicity variable used in regression.
library(dplyr)
library(survey)
library(gt)
options(survey.lonely.psu = "adjust")
# ============================================================
# TABLE 3: SURVEY-WEIGHTED LOGISTIC REGRESSION
# Guam BRFSS 2023–2024
# ============================================================
# ------------------------------------------------------------
# 1. PREPARE REGRESSION DATA
# ------------------------------------------------------------
regression_data <- combined_years %>%
mutate(
gout_status = factor(
gout_status,
levels = c("No", "Yes")
),
early_gout_onset = factor(
early_gout_onset,
levels = c("No", "Yes")
),
betel_nut = factor(
betel_nut,
levels = c("No", "Yes")
),
betel_nut_current = factor(
betel_nut_current,
levels = c("No", "Yes")
),
sex = factor(
sex,
levels = c("Male", "Female")
),
ethnic_regression = factor(
ethnic_regression,
levels = c("Pacific", "Asian", "Other")
),
age_num = as.numeric(age_num)
) %>%
select(
combine_final_weight,
stratum,
sampling_unit,
gout_status,
early_gout_onset,
betel_nut,
betel_nut_current,
age_num,
sex,
ethnic_regression
)
# ------------------------------------------------------------
# 2. CREATE SURVEY DESIGN
# ------------------------------------------------------------
regression_survey <- svydesign(
ids = ~sampling_unit,
strata = ~stratum,
weights = ~combine_final_weight,
data = regression_data,
nest = TRUE
)Model I — Early-onset gout and ever areca nut use. The outcome was early-onset gout, and the primary exposure was ever having used areca nut. The model adjusted for sex and ethnicity, without age as an additional covariate.
Model II — Overall gout and ever areca nut use. The outcome was gout status, and the primary exposure was ever having used areca nut. The model adjusted for age (continuous, per year), sex, and ethnicity.
Model III — Early-onset gout and current areca nut use. The outcome was early-onset gout, and the primary exposure was current areca nut use. The model adjusted for sex and ethnicity, without age as an additional covariate.
Model IV — Overall gout and current areca nut use. The outcome was gout status, and the primary exposure was current areca nut use. The model adjusted for age (continuous, per year), sex, and ethnicity.
The No category for the respective areca nut exposure is the comparison group in each model. The early-onset outcome used here follows the recoding in the pooled regression dataset; its handling of later-onset gout should be reconciled with the descriptive prevalence definition before publication.
# ------------------------------------------------------------
# 3. DEFINE THE FOUR REGRESSION MODELS
# ------------------------------------------------------------
# Model I: Early-onset gout and ever areca nut use
model_1 <- svyglm(
early_gout_onset ~
betel_nut +
sex +
ethnic_regression,
design = regression_survey,
family = quasibinomial()
)
# Model II: Overall gout and ever areca nut use
model_2 <- svyglm(
gout_status ~
betel_nut +
age_num +
sex +
ethnic_regression,
design = regression_survey,
family = quasibinomial()
)
# Model III: Early-onset gout and current areca nut use
model_3 <- svyglm(
early_gout_onset ~
betel_nut_current +
sex +
ethnic_regression,
design = regression_survey,
family = quasibinomial()
)
# Model IV: Overall gout and current areca nut use
model_4 <- svyglm(
gout_status ~
betel_nut_current +
age_num +
sex +
ethnic_regression,
design = regression_survey,
family = quasibinomial()
)For each fitted model, regression coefficients were exponentiated to obtain adjusted odds ratios. Confidence intervals were calculated from model-based coefficient estimates and survey-design-adjusted standard errors using the residual degrees of freedom. P values were obtained from the survey-weighted model summaries. Values are rounded for presentation only; calculations use the unrounded estimates.
# ------------------------------------------------------------
# 4. EXTRACT ODDS RATIOS, 95% CI AND P VALUES
# ------------------------------------------------------------
extract_regression <- function(model, model_name) {
beta <- coef(model)
se <- sqrt(diag(vcov(model)))
# Survey-adjusted confidence intervals
df <- model$df.residual
critical_value <- qt(0.975, df = df)
ci_lower <- beta - critical_value * se
ci_upper <- beta + critical_value * se
p_values <- summary(model)$coefficients[, 4]
data.frame(
Term = names(beta),
Model = model_name,
OR = exp(beta),
CI_Lower = exp(ci_lower),
CI_Upper = exp(ci_upper),
P = p_values,
row.names = NULL
)
}
# Extract all models
results_1 <- extract_regression(model_1, "Model I")
results_2 <- extract_regression(model_2, "Model II")
results_3 <- extract_regression(model_3, "Model III")
results_4 <- extract_regression(model_4, "Model IV")
# ------------------------------------------------------------
# 5. FORMAT REGRESSION RESULTS
# ------------------------------------------------------------
regression_results <- bind_rows(
results_1,
results_2,
results_3,
results_4
) %>%
filter(Term != "(Intercept)") %>%
mutate(
Predictor = recode(
Term,
"age_num" = "Age",
"betel_nutYes" = "Areca nut (Ever)",
"betel_nut_currentYes" = "Areca nut (Current)",
"ethnic_regressionAsian" = "Asian",
"ethnic_regressionOther" = "Other",
"sexFemale" = "Female"
),
OR_CI = sprintf(
"%.2f (%.2f–%.2f)",
OR,
CI_Lower,
CI_Upper
),
P_formatted = case_when(
is.na(P) ~ "—",
P < 0.01 ~ "<0.01",
TRUE ~ sprintf("%.2f", P)
)
)Results from the four models were combined into a single table. Ethnicity and sex are shown as grouped predictors, with explicit reference rows for Pacific Islander and Male. Dashes indicate reference categories or terms that were not included in a particular model. The finished table is presented in the final section alongside the descriptive tables.
# ============================================================
# 6. CREATE COMBINED TABLE 3 WITH REFERENCE ROWS
# ============================================================
predictor_order <- c(
"Age",
"Areca nut (Ever)",
"Areca nut (Current)",
"Ethnicity",
"Pacific Islander (Reference)",
"Asian",
"Other",
"Sex",
"Male (Reference)",
"Female"
)
# Convert regression results into wide format
table_3_data <- regression_results %>%
select(
Predictor,
Model,
OR_CI,
P_formatted
) %>%
tidyr::pivot_wider(
names_from = Model,
values_from = c(OR_CI, P_formatted),
names_glue = "{Model}_{.value}"
)
# Add explicit reference and section-header rows
reference_rows <- tibble(
Predictor = c(
"Ethnicity",
"Pacific Islander (Reference)",
"Sex",
"Male (Reference)"
)
)
table_3_data <- reference_rows %>%
bind_rows(table_3_data)
# Make sure all model columns exist
model_names <- c(
"Model I",
"Model II",
"Model III",
"Model IV"
)
for (m in model_names) {
for (suffix in c("OR_CI", "P_formatted")) {
col <- paste0(m, "_", suffix)
if (!col %in% names(table_3_data)) {
table_3_data[[col]] <- NA_character_
}
}
}
# Arrange predictors in manuscript order
table_3_data <- table_3_data %>%
mutate(
Predictor = factor(
Predictor,
levels = predictor_order
)
) %>%
arrange(Predictor) %>%
mutate(
Predictor = as.character(Predictor)
) %>%
select(
Predictor,
`Model I_OR_CI`,
`Model I_P_formatted`,
`Model II_OR_CI`,
`Model II_P_formatted`,
`Model III_OR_CI`,
`Model III_P_formatted`,
`Model IV_OR_CI`,
`Model IV_P_formatted`
) %>%
mutate(
across(
-Predictor,
~tidyr::replace_na(.x, "—")
)
)
# Section-header rows should have blank cells,
# while reference-category rows should show dashes.
table_3_data <- table_3_data %>%
mutate(
across(
-Predictor,
~if_else(
Predictor %in% c("Ethnicity", "Sex"),
"",
.x
)
)
)
# ============================================================
# 7. FORMAT TABLE 3 WITH GT
# ============================================================
table_3_table_object <- table_3_data %>%
gt() %>%
tab_header(
title = paste(
"Table 3. Multiple Regression Analysis",
"of Areca Nut Use and Gout"
),
subtitle = "Guam BRFSS 2023–2024"
) %>%
cols_label(
Predictor = "Predictor",
`Model I_OR_CI` = "OR (95% CI)",
`Model I_P_formatted` = "P",
`Model II_OR_CI` = "OR (95% CI)",
`Model II_P_formatted` = "P",
`Model III_OR_CI` = "OR (95% CI)",
`Model III_P_formatted` = "P",
`Model IV_OR_CI` = "OR (95% CI)",
`Model IV_P_formatted` = "P"
) %>%
tab_spanner(
label = "Model I: Early-onset gout / Ever use",
columns = c(
`Model I_OR_CI`,
`Model I_P_formatted`
)
) %>%
tab_spanner(
label = "Model II: Overall gout / Ever use",
columns = c(
`Model II_OR_CI`,
`Model II_P_formatted`
)
) %>%
tab_spanner(
label = "Model III: Early-onset gout / Current use",
columns = c(
`Model III_OR_CI`,
`Model III_P_formatted`
)
) %>%
tab_spanner(
label = "Model IV: Overall gout / Current use",
columns = c(
`Model IV_OR_CI`,
`Model IV_P_formatted`
)
) %>%
cols_align(
align = "left",
columns = Predictor
) %>%
cols_align(
align = "center",
columns = -Predictor
) %>%
# Bold the ethnicity and sex section headers
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
rows = Predictor %in% c("Ethnicity", "Sex")
)
) %>%
# Italicize reference-category labels
tab_style(
style = cell_text(style = "italic"),
locations = cells_body(
columns = Predictor,
rows = Predictor %in% c(
"Pacific Islander (Reference)",
"Male (Reference)"
)
)
) %>%
tab_source_note(
source_note = paste(
"Survey-weighted logistic regression.",
"OR = odds ratio; CI = confidence interval.",
"Reference categories: Pacific Islander",
"ethnicity, male sex, and no areca nut use.",
"Models I and III adjust for sex and ethnicity.",
"Models II and IV additionally adjust for age.",
"Dashes indicate reference categories",
"or predictors not included in a model."
)
)
# Display Table 3
# The formatted table is displayed in the final Tables section.
invisible(table_3_table_object)This section presents the final formatted tables for the pooled Guam BRFSS 2023–2024 analyses. Frequencies are unweighted; prevalence estimates and regression results account for the specified survey design and weights. The R code used to generate each table can be expanded using the HTML Code controls in the sections above.
| Table 1. Demographic Description | ||
| Guam BRFSS 2023–2024 | ||
| Variable | Frequency (N = 3,219) | Proportion (95% CI) |
|---|---|---|
| Sex | ||
| Female | 1,546 | 49.4% (45.1% – 53.8%) |
| Male | 1,673 | 50.6% (46.2% – 54.9%) |
| Age Category | ||
| 18 to 24 | 261 | 11.6% (9.4% – 14.2%) |
| 25 to 34 | 435 | 22.3% (19.2% – 25.8%) |
| 35 to 44 | 477 | 16.7% (13.4% – 20.6%) |
| 45 to 54 | 523 | 17.1% (14.1% – 20.5%) |
| 55 to 64 | 636 | 16.1% (13.0% – 19.9%) |
| 65 or older | 865 | 16.2% (13.0% – 19.9%) |
| Ethnicity | ||
| CHamoru | 1,365 | 43.6% (39.2% – 48.0%) |
| Filipino | 858 | 30.1% (26.1% – 34.4%) |
| Other Micronesian | 219 | 7.6% (5.9% – 9.9%) |
| Other | 777 | 18.7% (16.2% – 21.5%) |
| Education Level | ||
| College | 1,944 | 48.2% (43.8% – 52.6%) |
| High School or Below | 1,261 | 51.8% (47.4% – 56.2%) |
| Frequencies are unweighted. Proportions and 95% confidence intervals are survey-weighted. Missing responses are excluded from variable-specific denominators. | ||
| Table 2. Prevalence of Selected Health Conditions and Behaviors | ||
| Guam BRFSS 2023–2024 | ||
| Variable | Frequency (N = 3,219) | Proportion (95% CI) |
|---|---|---|
| Health Conditions | ||
| Gout | 336 | 13.5% (9.9% – 18.0%) |
| Early-onset gout | 131 | 5.1% (3.6% – 7.3%) |
| Obesity | 1,019 | 37.8% (33.6% – 42.3%) |
| Diabetes | 788 | 31.4% (26.7% – 36.5%) |
| Kidney disease | 150 | 3.7% (2.8% – 5.0%) |
| Health Behaviors | ||
| Alcohol consumption | 1,257 | 42.9% (38.6% – 47.3%) |
| Smoking | 468 | 18.2% (14.0% – 23.2%) |
| Smokeless tobacco use | 85 | 3.9% (2.7% – 5.8%) |
| Areca nut use (ever used) | 1,162 | 29.2% (25.9% – 32.6%) |
| Areca nut use (currently use) | 214 | 7.4% (5.7% – 9.6%) |
| 1. Proportions are survey-weighted using pooled Guam BRFSS 2023 and 2024 weights. 2. Parentheses indicate 95% confidence intervals. 3. Frequencies are unweighted. Missing responses are excluded from variable-specific denominators. 4. Early-onset gout is defined as gout diagnosed before age 40. | ||
Models I and III assess early-onset gout, while Models II and IV assess overall gout. Models I and II use ever areca nut use as the exposure; Models III and IV use current use. Reference categories for ethnicity and sex are shown explicitly in the table.
| Table 3. Multiple Regression Analysis of Areca Nut Use and Gout | ||||||||
| Guam BRFSS 2023–2024 | ||||||||
| Predictor |
Model I: Early-onset gout / Ever use
|
Model II: Overall gout / Ever use
|
Model III: Early-onset gout / Current use
|
Model IV: Overall gout / Current use
|
||||
|---|---|---|---|---|---|---|---|---|
| OR (95% CI) | P | OR (95% CI) | P | OR (95% CI) | P | OR (95% CI) | P | |
| Age | — | — | 1.05 (1.03–1.07) | <0.01 | — | — | 1.05 (1.03–1.07) | <0.01 |
| Areca nut (Ever) | 2.24 (1.15–4.37) | 0.02 | 0.93 (0.45–1.93) | 0.84 | — | — | — | — |
| Areca nut (Current) | — | — | — | — | 2.06 (0.67–6.30) | 0.20 | 1.46 (0.55–3.87) | 0.45 |
| Ethnicity | ||||||||
| Pacific Islander (Reference) | — | — | — | — | — | — | — | — |
| Asian | 1.42 (0.66–3.05) | 0.37 | 1.72 (0.78–3.80) | 0.18 | 1.10 (0.46–2.62) | 0.83 | 1.89 (0.92–3.87) | 0.08 |
| Other | 0.23 (0.09–0.58) | <0.01 | 0.39 (0.18–0.84) | 0.02 | 0.20 (0.08–0.52) | <0.01 | 0.40 (0.19–0.85) | 0.02 |
| Sex | ||||||||
| Male (Reference) | — | — | — | — | — | — | — | — |
| Female | 0.09 (0.05–0.19) | <0.01 | 0.16 (0.09–0.31) | <0.01 | 0.09 (0.04–0.18) | <0.01 | 0.17 (0.09–0.31) | <0.01 |
| Survey-weighted logistic regression. OR = odds ratio; CI = confidence interval. Reference categories: Pacific Islander ethnicity, male sex, and no areca nut use. Models I and III adjust for sex and ethnicity. Models II and IV additionally adjust for age. Dashes indicate reference categories or predictors not included in a model. | ||||||||
The BRFSS datasets are not listed here as the investigators are not approved to disseminate the data. A formal request can be made to the Guam BRFSS state coordinator. Once the form has been submitted and processed, arrangements can be made in terms of access to the specified data. Request form: https://dphss.guam.gov/sites/default/files/OER%20Data%20Request%20Form.pdf