This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
library(survey)## Warning: package 'survey' was built under R version 4.3.3## Loading required package: grid## Loading required package: Matrix## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchartlibrary(tidyverse)## Warning: package 'tidyverse' was built under R version 4.3.3## Warning: package 'ggplot2' was built under R version 4.3.3## Warning: package 'tidyr' was built under R version 4.3.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ tidyr::pack() masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(mice)## Warning: package 'mice' was built under R version 4.3.3##
## Attaching package: 'mice'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## cbind, rbindlibrary(dplyr)
library(tidyr)
library(stringr)
library(magrittr)##
## Attaching package: 'magrittr'
##
## The following object is masked from 'package:purrr':
##
## set_names
##
## The following object is masked from 'package:tidyr':
##
## extractlibrary(foreign)
library(svrep)## Warning: package 'svrep' was built under R version 4.3.3library(psych)##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(mitools)## Warning: package 'mitools' was built under R version 4.3.3library(srvyr)##
## Attaching package: 'srvyr'
##
## The following object is masked from 'package:stats':
##
## filterlibrary(gt)
library(scales)##
## Attaching package: 'scales'
##
## The following objects are masked from 'package:psych':
##
## alpha, rescale
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor# --- Step 0: Setup and Data Import ---
set.seed(123)
options(survey.adjust.domain.lonely=TRUE)
options(survey.lonely.psu="adjust")
# Load your BRFSS dataset (example for Cameron county)
TX_BRFSS_23 <- read.dta("E:/TX_DSHS_BRFSS/2023/PUDF/state_23_working_pudf.dta")
# Define variables of interest and predictors clearly:
variables <- c("c09q02", "cholch1", "educat3a", "pacat", "rfsmok",
"rfdrhv2", "highbp", "bmicat5", "diabetes", "c09q01", "c07q09", "c08q04", "agegr3")
predictors <- c("fairpoor", "sex", "c08q02a", "c08q03a", "nattmpts",
"llcpwt", "geostr", "ststr", "agegr3")
#build the key
resp_key <- data.frame(
var = c("c09q02", "cholch1", "educat3a", "pacat", "rfsmok", "rfdrhv2", "highbp",
"bmicat5", "diabetes", "c09q01", "c07q09", "c08q04"),
key = c(1, 1, 1, "3-4", 2, 2, 1, "4-5", 1, 1, 1, "2-5"), lancet_modifiable = c("VISION LOSS", "HIGH CHOLESTEROL", "LESS EDUCATION", "PHYSICAL INACTIVITY", "SMOKING", "EXCESSIVE ALCOHOL", "HYPERTENSION", "OBESITY", "DIABETES", "HEARING LOSS", "DEPRESSION", "SOCIAL ISOLATION"), lancet24_RR = c(1.5, 1.3, 1.6, 1.2, 1.3, 1.2, 1.2, 1.3, 1.7, 1.4, 2.2, 1.6), upperRR = c(1.6, 1.4, 2.0, 1.3, 1.4, 1.5, 1.4, 1.7, 1.8, 1.9, 3.0, 1.8), lowerRR = c(1.4, 1.3, 1.3, 1.2, 1.2, 1.0, 1.1, 1.0, 1.6, 1.0, 1.7, 1.3))
# Filter data for Cameron county
cam_data <- TX_BRFSS_23 %>%
filter(countyname == "Cameron") %>%
select(all_of(c(predictors, variables))) %>%
mutate(across(where(is.character), as.factor)) # convert characters to factors
# Convert survey weights and strata properly
cam_data$llcpwt <- as.numeric(cam_data$llcpwt)
cam_data$ststr <- as.factor(cam_data$ststr)
# --- Step 1: Build Predictor Matrix for MICE ---
# Include survey design variables explicitly in imputation
# Initialize mice with maxit=0 to get default predictor matrix
imp_init <- mice(cam_data, maxit=0, print=FALSE)
pred_matrix <- imp_init$predictorMatrix
# Set predictors: include weights (llcpwt) and strata (ststr) as predictors but not imputed
pred_matrix[, "llcpwt"] <- 1 # use weights as predictors
pred_matrix[, "agegr3"] <- 1 # use age as predictors
pred_matrix["llcpwt", ] <- 0 # don't impute weights
pred_matrix["ststr", ] <- 0 # don't impute strata
pred_matrix["diabetes", ] <- 1 # DO impute diabetes, I wonder why it wasn't!
pred_matrix[,"ststr" ] <- 0 # don't use strata for prediction
# You can customize further, e.g., limit predictors to only relevant variables
# --- Step 2: Define imputation methods ---
meth <- imp_init$method
meth["llcpwt"] <- "" # no imputation on weights
meth["ststr"] <- "" # no imputation on strata
# Using predictive mean matching, "pmm", for categorical variables
## THIS CAN BE *REALLY* CONTENTIOUS, but the 'pmm' approach captures tails fairly well compared to 'rf'
## some lit to consider: https://arxiv.org/html/2401.09602v1 and https://doi.org/10.1186/s12874-018-0615-6
meth[variables] <- "rf"
# --- Step 3: Run Multiple Imputation ---
## When dealing with categorical variables/factors, you've gotta do it visually by inspecting the distribution and ensuring that each level/factor that's been imputed is a) representative of all possible choices (i.e. that there aren't any choices or options missing (you should cut those)), and b) that the imputed values are distributed fairly evenly across the categories available
#thanks to Dr. Nerler (https://nerler.github.io/EP16_Multiple_Imputation/slide/07_convergence_and_diagnostics.pdf)
#code below adapted from hers here: https://gist.github.com/NErler/0d00375da460dd33839b98faeee2fdab
#DIAGNOSTIC FUNCTION FOR VISUALIZATION OF IMPUTED VALUE DISTRIBUTIONS across ONE type of imputation method across QUESTIONS (not method v. method unless called twice)
plot_categorical_convergence <- function(imputed_data, variables) {
# Ensure variables are a character vector
if (!is.character(variables)) {
stop("`variables` must be a character vector.")
}
if (inherits(imputed_data, "mids")) {
imputed_data <- complete(imputed_data, action = "long", include = TRUE) %>%
mutate(across(where(is.character), as.factor))
}
imputed_data %>%
filter(.imp != 0) %>%
pivot_longer(cols = all_of(variables), names_to = "variable", values_to = "value") %>%
filter(!is.na(value)) %>%
group_by(.imp, variable, value) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(.imp, variable) %>%
mutate(prop = n / sum(n)) %>%
ggplot(aes(x = value, y = prop, fill = factor(.imp))) +
geom_bar(stat = "identity", position = "dodge") +
facet_wrap(~ variable, scales = "free_x") +
labs(
title = "Convergence of Categorical Variables Across Imputations",
x = "Category",
y = "Proportion",
fill = "Imputation"
) +
theme_minimal()
}
#####
##
get_prop_table <- function(data, label = "original") {
data %>%
pivot_longer(cols = all_of(variables), names_to = "variable", values_to = "value") %>%
filter(!is.na(value)) %>%
group_by(variable, value) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(variable) %>%
mutate(
prop = n / sum(n),
method = label
) %>%
select(method, variable, value, prop)
}
get_imputation_props <- function(mids_obj, method_label) {
complete(mids_obj, action = "long", include = TRUE) %>%
pivot_longer(cols = all_of(variables), names_to = "variable", values_to = "value") %>%
filter(!is.na(value)) %>%
group_by(.imp, variable, value) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(.imp, variable) %>%
mutate(prop = n / sum(n)) %>%
group_by(variable, value) %>%
summarise(
mean_prop = mean(prop),
sd_prop = sd(prop),
method = method_label,
.groups = "drop"
)
}
#
# #########
# #want to compare convergences? they're all really similar
# install.packages(patchwork)
# library(patchwork)
# #
# imp_test4 <- mice(cam_data, m=20, method = "rf", predictorMatrix = pred_matrix, print=FALSE)
# imp_test3 <- mice(cam_data, m=20, method = "cart", predictorMatrix = pred_matrix, print=FALSE)
# imp_test2 <- mice(cam_data, m=20, method = "sample", predictorMatrix = pred_matrix, print=FALSE)
# imp_test1 <- mice(cam_data, m=20, method = "midastouch", predictorMatrix = pred_matrix, print=FALSE)
# imp <- mice(cam_data, m=20, method = "pmm", predictorMatrix = pred_matrix, print=FALSE)
# # Original data (complete cases only)
# original_props <- get_prop_table(TX_BRFSS_23, label = "original")
#
# # Example for two methods
# pmm_props <- get_imputation_props(imp, "pmm")
# cart_props <- get_imputation_props(imp_test3, "cart")
# midas_props <- get_imputation_props(imp_test2, "midastouch")
#
#
# # Combine all
# all_props <- bind_rows(
# original_props %>% rename(mean_prop = prop),
# pmm_props,
# cart_props,
# midas_props
# )
#
#
# diffs <- all_props %>%
# filter(method != "original") %>%
# left_join(original_props, by = c("variable", "value"), suffix = c("_imp", "_orig")) %>%
# mutate(abs_diff = abs(mean_prop - prop))
#
# # Summarise per method and variable
# summary_diff <- diffs %>%
# group_by(method_imp, variable) %>%
# summarise(
# mean_abs_diff = mean(abs_diff, na.rm = TRUE),
# max_abs_diff = max(abs_diff, na.rm = TRUE),
# .groups = "drop"
# )
#
# ggplot(all_props, aes(x = value, y = mean_prop, color = method, group = method, fill = method)) +
# geom_col(position = "dodge") +
# facet_wrap(~ variable, scales = "free_x") +
# labs(title = "Proportion Comparison Across Methods", y = "Proportion", x = "Category")
# plot_categorical_convergence(imp_test1, variables) -> imp_catdistro
# plot_categorical_convergence(imp_test1, variables) -> midastouch_catdistro
# plot_categorical_convergence(imp_test2, variables) -> sample_catdistro
# plot_categorical_convergence(imp_test3, variables) -> cart_catdistro
# plot_categorical_convergence(imp_test4, variables) -> rf_catdistro
#
# imp_catdistro + midastouch_catdistro + sample_catdistro + cart_catdistro + rf_catdistro
#
#
#
# ###
#begin imputation
imp <- mice(cam_data, m=20, method = meth, predictorMatrix = pred_matrix, print=TRUE)##
## iter imp variable
## 1 1 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 2 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 3 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 4 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 5 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 6 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 7 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 8 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 9 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 10 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 11 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 12 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 13 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 14 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 15 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 16 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 17 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 18 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 19 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 1 20 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 1 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 2 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 3 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 4 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 5 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 6 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 7 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 8 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 9 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 10 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 11 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 12 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 13 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 14 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 15 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 16 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 17 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 18 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 19 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 2 20 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 1 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 2 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 3 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 4 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 5 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 6 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 7 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 8 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 9 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 10 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 11 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 12 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 13 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 14 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 15 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 16 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 17 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 18 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 19 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 3 20 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 1 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 2 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 3 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 4 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 5 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 6 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 7 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 8 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 9 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 10 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 11 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 12 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 13 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 14 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 15 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 16 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 17 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 18 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 19 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 4 20 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 1 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 2 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 3 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 4 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 5 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 6 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 7 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 8 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 9 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 10 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 11 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 12 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 13 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 14 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 15 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 16 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 17 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 18 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 19 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04
## 5 20 agegr3 c09q02 cholch1 educat3a pacat rfsmok rfdrhv2 highbp bmicat5 c09q01 c07q09 c08q04cam_BRFSS_23_core_imp <- complete(imp, action = "long", include = TRUE) %>% mutate(across(where(is.character), as.factor))
#cam_BRFSS_23_core_imp_svy<- svydesign(id = ~1, strata = ~ststr, weights = ~llcpwt, data = cam_BRFSS_23_core_imp)
## Make the different survey_design objects necessary to feed into 'survey' and 'svytotal'
svydesign_objs<- split(cam_BRFSS_23_core_imp %>% mutate(llcpwt = as.numeric(llcpwt)), cam_BRFSS_23_core_imp$.imp) %>% imputationList(.) %>% svydesign(id = ~1, strata = ~ststr, weights = ~llcpwt, data = .)
### Compute estimate/point totals for all Lancet-related variables of interest and, while we're at it, separate the variable of interest + the response to make keying, comparison, and evaluation easier
lapply( variables , function( z ) with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",z,"),","factor(agegr3))")), na.rm = TRUE)) %>% MIcombine() %>% summary(.)) %>% do.call(rbind, .) %>% as.data.frame(.) %>% mutate(var = rownames(.), var = var %>% str_remove(., paste0(" factor\\(agegr3\\)\\)")) %>% str_remove(., paste0("interaction\\(factor\\(")) %>% str_remove(., "\\)")) %>% separate_wider_delim(., cols = "var", delim = ",", names_sep="", too_few = "align_start") %>% separate_wider_delim(., cols = "var2", delim = ".", names_sep="", too_few = "align_start") %>% dplyr::rename(., resp=7, agegr3_resp = 8) %>% filter(var1!= "agegr3") %>% left_join(., resp_key, by=c("var1"="var")) %>% group_by(var1, agegr3_resp) %>% mutate(prop = 100 * (results/sum(results))) %>% ungroup(.) %>% arrange(., desc(prop)) -> pooled_imp_lancetvars_23camBRFSS## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(c09q02), factor(agegr3))1.1 12621.14 6926.384 -954.8277
## interaction(factor(c09q02), factor(agegr3))2.1 152244.53 13525.404 125732.1005
## interaction(factor(c09q02), factor(agegr3))1.2 15046.59 3816.293 7566.0429
## interaction(factor(c09q02), factor(agegr3))2.2 108916.36 11581.147 86216.4672
## interaction(factor(c09q02), factor(agegr3))1.3 15564.89 3289.225 9117.8359
## interaction(factor(c09q02), factor(agegr3))2.3 41513.81 5020.040 31673.0810
## upper) missInfo
## interaction(factor(c09q02), factor(agegr3))1.1 26197.10 2 %
## interaction(factor(c09q02), factor(agegr3))2.1 178756.96 4 %
## interaction(factor(c09q02), factor(agegr3))1.2 22527.14 4 %
## interaction(factor(c09q02), factor(agegr3))2.2 131616.25 3 %
## interaction(factor(c09q02), factor(agegr3))1.3 22011.95 3 %
## interaction(factor(c09q02), factor(agegr3))2.3 51354.54 5 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(cholch1), factor(agegr3))1.1 30784.89 7901.738 15207.80
## interaction(factor(cholch1), factor(agegr3))2.1 133039.68 15718.055 102147.60
## interaction(factor(cholch1), factor(agegr3))1.2 48028.41 7733.445 32869.06
## interaction(factor(cholch1), factor(agegr3))2.2 75939.42 10007.674 56323.89
## interaction(factor(cholch1), factor(agegr3))1.3 32927.01 4623.788 23860.24
## interaction(factor(cholch1), factor(agegr3))2.3 24053.47 4061.200 16092.51
## upper) missInfo
## interaction(factor(cholch1), factor(agegr3))1.1 46361.99 32 %
## interaction(factor(cholch1), factor(agegr3))2.1 163931.76 22 %
## interaction(factor(cholch1), factor(agegr3))1.2 63187.76 5 %
## interaction(factor(cholch1), factor(agegr3))2.2 95554.96 3 %
## interaction(factor(cholch1), factor(agegr3))1.3 41993.78 9 %
## interaction(factor(cholch1), factor(agegr3))2.3 32014.42 5 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(educat3a), factor(agegr3))1.1 28356.68 8455.844 11783.524
## interaction(factor(educat3a), factor(agegr3))2.1 72151.44 9965.326 52619.757
## interaction(factor(educat3a), factor(agegr3))3.1 64921.78 9531.574 46240.237
## interaction(factor(educat3a), factor(agegr3))1.2 40612.65 8532.761 23888.707
## interaction(factor(educat3a), factor(agegr3))2.2 26066.18 5120.666 16029.856
## interaction(factor(educat3a), factor(agegr3))3.2 57632.46 7980.051 41991.849
## interaction(factor(educat3a), factor(agegr3))1.3 36153.88 5308.713 25748.747
## interaction(factor(educat3a), factor(agegr3))2.3 6065.62 1555.498 3016.838
## interaction(factor(educat3a), factor(agegr3))3.3 14955.15 2305.604 10436.227
## upper) missInfo
## interaction(factor(educat3a), factor(agegr3))1.1 44929.830 0 %
## interaction(factor(educat3a), factor(agegr3))2.1 91683.116 0 %
## interaction(factor(educat3a), factor(agegr3))3.1 83603.321 0 %
## interaction(factor(educat3a), factor(agegr3))1.2 57336.600 1 %
## interaction(factor(educat3a), factor(agegr3))2.2 36102.499 0 %
## interaction(factor(educat3a), factor(agegr3))3.2 73273.076 0 %
## interaction(factor(educat3a), factor(agegr3))1.3 46559.019 2 %
## interaction(factor(educat3a), factor(agegr3))2.3 9114.402 2 %
## interaction(factor(educat3a), factor(agegr3))3.3 19474.077 1 %## Warning in onestrat(`attr<-`(x[index, , drop = FALSE], "recentering",
## recentering), : Stratum (482112) has only one PSU at stage 1## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(pacat), factor(agegr3))1.1 1843.8975 1855.3107 -1804.5855
## interaction(factor(pacat), factor(agegr3))4.1 157976.5185 29977.0824 97052.3475
## interaction(factor(pacat), factor(agegr3))1.2 636.2548 1174.0579 -1716.7734
## interaction(factor(pacat), factor(agegr3))4.2 119757.3318 21721.2530 75810.0189
## interaction(factor(pacat), factor(agegr3))1.3 369.0328 621.5677 -873.6892
## interaction(factor(pacat), factor(agegr3))4.3 55331.5626 9160.4655 36921.6131
## upper) missInfo
## interaction(factor(pacat), factor(agegr3))1.1 5492.381 24 %
## interaction(factor(pacat), factor(agegr3))4.1 218900.690 78 %
## interaction(factor(pacat), factor(agegr3))1.2 2989.283 62 %
## interaction(factor(pacat), factor(agegr3))4.2 163704.645 73 %
## interaction(factor(pacat), factor(agegr3))1.3 1611.755 58 %
## interaction(factor(pacat), factor(agegr3))4.3 73741.512 65 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se
## interaction(factor(rfsmok), factor(agegr3))1.1 150338.650 14739.7220
## interaction(factor(rfsmok), factor(agegr3))2.1 14430.226 4680.4356
## interaction(factor(rfsmok), factor(agegr3))1.2 111756.811 11496.0893
## interaction(factor(rfsmok), factor(agegr3))2.2 12175.644 4000.7756
## interaction(factor(rfsmok), factor(agegr3))1.3 55815.461 5780.2729
## interaction(factor(rfsmok), factor(agegr3))2.3 1131.953 577.0819
## (lower upper) missInfo
## interaction(factor(rfsmok), factor(agegr3))1.1 1.214456e+05 179231.674 5 %
## interaction(factor(rfsmok), factor(agegr3))2.1 5.256083e+03 23604.369 3 %
## interaction(factor(rfsmok), factor(agegr3))1.2 8.922340e+04 134290.220 3 %
## interaction(factor(rfsmok), factor(agegr3))2.2 4.333835e+03 20017.453 3 %
## interaction(factor(rfsmok), factor(agegr3))1.3 4.448209e+04 67148.831 8 %
## interaction(factor(rfsmok), factor(agegr3))2.3 7.251827e-01 2263.181 5 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se
## interaction(factor(rfdrhv2), factor(agegr3))1.1 147384.002 14601.697
## interaction(factor(rfdrhv2), factor(agegr3))2.1 17201.563 5808.462
## interaction(factor(rfdrhv2), factor(agegr3))1.2 112480.126 12153.991
## interaction(factor(rfdrhv2), factor(agegr3))2.2 10837.429 4408.949
## interaction(factor(rfdrhv2), factor(agegr3))1.3 53825.937 5778.385
## interaction(factor(rfdrhv2), factor(agegr3))2.3 3053.445 1456.558
## (lower upper) missInfo
## interaction(factor(rfdrhv2), factor(agegr3))1.1 118756.5522 176011.451 7 %
## interaction(factor(rfdrhv2), factor(agegr3))2.1 5816.7725 28586.354 2 %
## interaction(factor(rfdrhv2), factor(agegr3))1.2 88622.6538 136337.598 16 %
## interaction(factor(rfdrhv2), factor(agegr3))2.2 2193.9407 19480.917 6 %
## interaction(factor(rfdrhv2), factor(agegr3))1.3 42491.0580 65160.816 12 %
## interaction(factor(rfdrhv2), factor(agegr3))2.3 196.2337 5910.656 12 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(highbp), factor(agegr3))1.1 15483.83 4505.249 6653.683
## interaction(factor(highbp), factor(agegr3))2.1 149946.07 14435.043 121653.902
## interaction(factor(highbp), factor(agegr3))1.2 47577.88 8636.641 30650.298
## interaction(factor(highbp), factor(agegr3))2.2 76718.18 9062.454 58956.094
## interaction(factor(highbp), factor(agegr3))1.3 42758.28 5128.022 32706.731
## interaction(factor(highbp), factor(agegr3))2.3 14369.51 2908.272 8669.372
## upper) missInfo
## interaction(factor(highbp), factor(agegr3))1.1 24313.97 1 %
## interaction(factor(highbp), factor(agegr3))2.1 178238.23 0 %
## interaction(factor(highbp), factor(agegr3))1.2 64505.46 1 %
## interaction(factor(highbp), factor(agegr3))2.2 94480.26 0 %
## interaction(factor(highbp), factor(agegr3))1.3 52809.83 4 %
## interaction(factor(highbp), factor(agegr3))2.3 20069.66 1 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(bmicat5), factor(agegr3))1.1 3958.642 2320.129 -594.2060
## interaction(factor(bmicat5), factor(agegr3))2.1 53267.669 10536.783 32609.6674
## interaction(factor(bmicat5), factor(agegr3))3.1 53393.323 9120.036 35514.5420
## interaction(factor(bmicat5), factor(agegr3))4.1 41498.374 8346.552 25135.8411
## interaction(factor(bmicat5), factor(agegr3))5.1 12478.553 4400.885 3850.2630
## interaction(factor(bmicat5), factor(agegr3))1.2 3634.637 2814.006 -1881.0565
## interaction(factor(bmicat5), factor(agegr3))2.2 17755.951 4550.714 8826.8042
## interaction(factor(bmicat5), factor(agegr3))3.2 42011.959 8108.472 26111.5842
## interaction(factor(bmicat5), factor(agegr3))4.2 46574.776 7818.510 31241.4126
## interaction(factor(bmicat5), factor(agegr3))5.2 13671.372 5101.770 3668.1445
## interaction(factor(bmicat5), factor(agegr3))1.3 1483.193 1046.619 -568.8483
## interaction(factor(bmicat5), factor(agegr3))2.3 11482.700 3129.362 5348.2898
## interaction(factor(bmicat5), factor(agegr3))3.3 22691.480 3565.259 15699.6089
## interaction(factor(bmicat5), factor(agegr3))4.3 19962.261 4079.923 11955.3157
## interaction(factor(bmicat5), factor(agegr3))5.3 1243.440 813.013 -357.1752
## upper) missInfo
## interaction(factor(bmicat5), factor(agegr3))1.1 8511.490 14 %
## interaction(factor(bmicat5), factor(agegr3))2.1 73925.670 7 %
## interaction(factor(bmicat5), factor(agegr3))3.1 71272.104 6 %
## interaction(factor(bmicat5), factor(agegr3))4.1 57860.908 6 %
## interaction(factor(bmicat5), factor(agegr3))5.1 21106.843 7 %
## interaction(factor(bmicat5), factor(agegr3))1.2 9150.331 3 %
## interaction(factor(bmicat5), factor(agegr3))2.2 26685.097 14 %
## interaction(factor(bmicat5), factor(agegr3))3.2 57912.333 9 %
## interaction(factor(bmicat5), factor(agegr3))4.2 61908.140 10 %
## interaction(factor(bmicat5), factor(agegr3))5.2 23674.600 8 %
## interaction(factor(bmicat5), factor(agegr3))1.3 3535.235 8 %
## interaction(factor(bmicat5), factor(agegr3))2.3 17617.111 5 %
## interaction(factor(bmicat5), factor(agegr3))3.3 29683.352 10 %
## interaction(factor(bmicat5), factor(agegr3))4.3 27969.207 15 %
## interaction(factor(bmicat5), factor(agegr3))5.3 2844.055 28 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se
## interaction(factor(diabetes), factor(agegr3))1.1 5018.128 3004.469
## interaction(factor(diabetes), factor(agegr3))2.1 160411.764 14535.148
## interaction(factor(diabetes), factor(agegr3))1.2 33294.963 6958.318
## interaction(factor(diabetes), factor(agegr3))2.2 91016.331 10171.111
## interaction(factor(diabetes), factor(agegr3))1.3 28754.362 4503.243
## interaction(factor(diabetes), factor(agegr3))2.3 28424.407 3959.290
## (lower upper) missInfo
## interaction(factor(diabetes), factor(agegr3))1.1 -870.5822 10906.84 1 %
## interaction(factor(diabetes), factor(agegr3))2.1 131923.3979 188900.13 0 %
## interaction(factor(diabetes), factor(agegr3))1.2 19656.8233 46933.10 1 %
## interaction(factor(diabetes), factor(agegr3))2.2 71081.3183 110951.34 0 %
## interaction(factor(diabetes), factor(agegr3))1.3 19927.7312 37580.99 3 %
## interaction(factor(diabetes), factor(agegr3))2.3 20664.3296 36184.48 1 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se
## interaction(factor(c09q01), factor(agegr3))1.1 3744.22396 2691.5377
## interaction(factor(c09q01), factor(agegr3))2.1 161075.26892 14823.3378
## interaction(factor(c09q01), factor(agegr3))7.1 46.17389 302.7849
## interaction(factor(c09q01), factor(agegr3))1.2 5645.68608 2456.8568
## interaction(factor(c09q01), factor(agegr3))2.2 116216.29198 11691.7216
## interaction(factor(c09q01), factor(agegr3))7.2 2155.70284 1528.4755
## interaction(factor(c09q01), factor(agegr3))1.3 8300.79909 2104.4819
## interaction(factor(c09q01), factor(agegr3))2.3 48056.49853 5461.7895
## interaction(factor(c09q01), factor(agegr3))7.3 721.40290 721.4150
## (lower upper) missInfo
## interaction(factor(c09q01), factor(agegr3))1.1 -1531.6315 9020.0795 4 %
## interaction(factor(c09q01), factor(agegr3))2.1 132020.0013 190130.5365 3 %
## interaction(factor(c09q01), factor(agegr3))7.1 -556.8225 649.1703 52 %
## interaction(factor(c09q01), factor(agegr3))1.2 828.9972 10462.3750 7 %
## interaction(factor(c09q01), factor(agegr3))2.2 93300.1107 139132.4732 2 %
## interaction(factor(c09q01), factor(agegr3))7.2 -840.0541 5151.4598 0 %
## interaction(factor(c09q01), factor(agegr3))1.3 4175.9802 12425.6180 2 %
## interaction(factor(c09q01), factor(agegr3))2.3 37349.9564 58763.0406 5 %
## interaction(factor(c09q01), factor(agegr3))7.3 -692.5445 2135.3503 0 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(c07q09), factor(agegr3))1.1 35499.167 9302.729 17266.002
## interaction(factor(c07q09), factor(agegr3))2.1 129748.497 12501.913 105245.002
## interaction(factor(c07q09), factor(agegr3))1.2 27873.706 6984.403 14184.528
## interaction(factor(c07q09), factor(agegr3))2.2 96437.588 10130.705 76581.741
## interaction(factor(c07q09), factor(agegr3))1.3 6655.962 2026.595 2683.909
## interaction(factor(c07q09), factor(agegr3))2.3 50515.247 5416.538 39898.662
## upper) missInfo
## interaction(factor(c07q09), factor(agegr3))1.1 53732.33 1 %
## interaction(factor(c07q09), factor(agegr3))2.1 154251.99 1 %
## interaction(factor(c07q09), factor(agegr3))1.2 41562.88 0 %
## interaction(factor(c07q09), factor(agegr3))2.2 116293.43 1 %
## interaction(factor(c07q09), factor(agegr3))1.3 10628.02 0 %
## interaction(factor(c07q09), factor(agegr3))2.3 61131.83 2 %## Warning in onestrat(`attr<-`(x[index, , drop = FALSE], "recentering",
## recentering), : Stratum (482111) has only one PSU at stage 1## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(c08q04), factor(agegr3))1.1 72698.5765 11721.9385 49717.3109
## interaction(factor(c08q04), factor(agegr3))2.1 2950.1431 2545.2479 -2059.9733
## interaction(factor(c08q04), factor(agegr3))3.1 0.0000 0.0000 NaN
## interaction(factor(c08q04), factor(agegr3))4.1 4865.4274 3388.3094 -1807.9226
## interaction(factor(c08q04), factor(agegr3))5.1 77160.7118 10444.0990 56689.5561
## interaction(factor(c08q04), factor(agegr3))6.1 7474.9326 3936.1950 -260.8414
## interaction(factor(c08q04), factor(agegr3))1.2 69851.8114 9599.4591 51037.1614
## interaction(factor(c08q04), factor(agegr3))2.2 20342.9289 4104.0255 12299.1825
## interaction(factor(c08q04), factor(agegr3))3.2 6404.2285 3983.7153 -1403.7113
## interaction(factor(c08q04), factor(agegr3))4.2 11728.0506 3396.0127 5071.9878
## interaction(factor(c08q04), factor(agegr3))5.2 12054.9407 4614.0283 3011.6023
## interaction(factor(c08q04), factor(agegr3))6.2 3894.9453 2126.0023 -271.9426
## interaction(factor(c08q04), factor(agegr3))1.3 37776.9940 5140.0271 27702.3729
## interaction(factor(c08q04), factor(agegr3))2.3 6085.9659 1837.0657 2485.3687
## interaction(factor(c08q04), factor(agegr3))3.3 9376.4476 1978.1298 5499.3825
## interaction(factor(c08q04), factor(agegr3))4.3 481.3567 299.3145 -106.2971
## interaction(factor(c08q04), factor(agegr3))5.3 2329.5143 1121.1035 132.1088
## interaction(factor(c08q04), factor(agegr3))6.3 1101.8866 1081.3965 -1017.6180
## upper) missInfo
## interaction(factor(c08q04), factor(agegr3))1.1 95679.842 7 %
## interaction(factor(c08q04), factor(agegr3))2.1 7960.259 27 %
## interaction(factor(c08q04), factor(agegr3))3.1 NaN NaN %
## interaction(factor(c08q04), factor(agegr3))4.1 11538.777 29 %
## interaction(factor(c08q04), factor(agegr3))5.1 97631.868 3 %
## interaction(factor(c08q04), factor(agegr3))6.1 15210.707 22 %
## interaction(factor(c08q04), factor(agegr3))1.2 88666.461 1 %
## interaction(factor(c08q04), factor(agegr3))2.2 28386.675 0 %
## interaction(factor(c08q04), factor(agegr3))3.2 14212.168 0 %
## interaction(factor(c08q04), factor(agegr3))4.2 18384.113 0 %
## interaction(factor(c08q04), factor(agegr3))5.2 21098.279 0 %
## interaction(factor(c08q04), factor(agegr3))6.2 8061.833 0 %
## interaction(factor(c08q04), factor(agegr3))1.3 47851.615 2 %
## interaction(factor(c08q04), factor(agegr3))2.3 9686.563 1 %
## interaction(factor(c08q04), factor(agegr3))3.3 13253.513 0 %
## interaction(factor(c08q04), factor(agegr3))4.3 1069.011 17 %
## interaction(factor(c08q04), factor(agegr3))5.3 4526.920 3 %
## interaction(factor(c08q04), factor(agegr3))6.3 3221.391 1 %
## Multiple imputation results:
## with(svydesign_objs, svytotal(as.formula(paste0("~interaction(factor(",
## z, "),", "factor(agegr3))")), na.rm = TRUE))
## MIcombine.default(.)
## results se (lower
## interaction(factor(agegr3), factor(agegr3))1.1 165429.89 14705.292 136608.0
## interaction(factor(agegr3), factor(agegr3))2.1 0.00 0.000 NaN
## interaction(factor(agegr3), factor(agegr3))3.1 0.00 0.000 NaN
## interaction(factor(agegr3), factor(agegr3))1.2 0.00 0.000 NaN
## interaction(factor(agegr3), factor(agegr3))2.2 124311.29 11723.610 101333.4
## interaction(factor(agegr3), factor(agegr3))3.2 0.00 0.000 NaN
## interaction(factor(agegr3), factor(agegr3))1.3 0.00 0.000 NaN
## interaction(factor(agegr3), factor(agegr3))2.3 0.00 0.000 NaN
## interaction(factor(agegr3), factor(agegr3))3.3 57178.77 5621.779 46160.0
## upper) missInfo
## interaction(factor(agegr3), factor(agegr3))1.1 194251.74 0 %
## interaction(factor(agegr3), factor(agegr3))2.1 NaN NaN %
## interaction(factor(agegr3), factor(agegr3))3.1 NaN NaN %
## interaction(factor(agegr3), factor(agegr3))1.2 NaN NaN %
## interaction(factor(agegr3), factor(agegr3))2.2 147289.17 0 %
## interaction(factor(agegr3), factor(agegr3))3.2 NaN NaN %
## interaction(factor(agegr3), factor(agegr3))1.3 NaN NaN %
## interaction(factor(agegr3), factor(agegr3))2.3 NaN NaN %
## interaction(factor(agegr3), factor(agegr3))3.3 68197.54 2 %rbind(pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "educat3a") %>% filter(resp == "1" & agegr3_resp==1),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "c09q02") %>% filter(resp == key & agegr3_resp == 3),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "c08q04") %>% filter(resp != "1" & resp != "6" & agegr3_resp == 3),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "bmicat5") %>% filter(resp == "4" | resp == "5") %>% filter(agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "c07q09") %>% filter(resp == "1" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "c09q01") %>% filter(resp == "1" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "cholch1") %>% filter(resp == "1" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "diabetes") %>% filter(resp == "1" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "highbp") %>% filter(resp == "1" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "pacat") %>% filter(resp == "3" | resp == "4") %>% filter(agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "rfdrhv2") %>% filter(resp == "2" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "rfsmok") %>% filter(resp == "2" & agegr3_resp == 2),
pooled_imp_lancetvars_23camBRFSS %>% filter(var1 == "c07q09") %>% filter(resp == "1" & agegr3_resp == 2)
) %>% unique(.) %>% group_by(lancet_modifiable) %>% mutate(prop = sum(prop)) %>% arrange(., desc(.$prop)) -> pooled_imp_lancetvars_23_camBRFSS_keyed
### IMPUTED VALUES BUILT: get to building values needed for PIF/PAF like the communality and W; load in the risk reduction factors from 2024 Lancet
## build imputation set for tetrachoric correlation matrix for *11* Lancet risk factors
## we're dropping physical inactivity for Cameron County calculations because it's missing 86% of its damn values
cam_BRFSS_23_core_imp %>% dplyr::rename(., imp=1) %>% filter(imp != 0) %>% mutate(
educat3a_bin = ifelse(educat3a == "1", 1, 0),
c09q02_bin = ifelse(c09q02 == "1", 1, 0),
c08q04_bin = ifelse(!(c08q04 %in% c("1", "6")), 1, 0),
bmicat5_bin = ifelse(bmicat5 %in% c("4", "5"), 1, 0),
c07q09_bin = ifelse(c07q09 == "1", 1, 0),
c09q01_bin = ifelse(c09q01 == "1", 1, 0),
cholch1_bin = ifelse(cholch1 == "1", 1, 0),
diabetes_bin = ifelse(diabetes == "1", 1, 0),
highbp_bin = ifelse(highbp == "1", 1, 0),
rfdrhv2_bin = ifelse(rfdrhv2 == "2", 1, 0),
rfsmok_bin = ifelse(rfsmok == "2", 1, 0)
) -> imps_for_tetrachor
# split them all, we're gonna wanna pool or average our work eventually
imps_for_tetra_list<- split(imps_for_tetrachor, imps_for_tetrachor$imp)
#build the tetrachor matrices
lapply(1:length(imps_for_tetra_list), function(i) imps_for_tetra_list[[i]] %>% .[,24:34] %>% tetrachoric(.) %>% .$rho) -> tetrachor_rhoslapply(1:length(imps_for_tetra_list), function(i) imps_for_tetra_list[[i]] %>% filter(agegr3==1) %>% .[,24:34] %>% tetrachoric(.) %>% .$rho) -> tetrachor_rhos_young## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 4 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 6 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 6 j = 3 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 7 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 4 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 1 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 9 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 8 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 11 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!lapply(1:length(imps_for_tetra_list), function(i) imps_for_tetra_list[[i]] %>% filter(agegr3==2) %>% .[,24:34] %>% tetrachoric(.) %>% .$rho) -> tetrachor_rhos_mid## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 8 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 2 A cell entry of 0 was replaced with correct = 0.5. Check your data!
## For i = 10 j = 6 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## donelapply(1:length(imps_for_tetra_list), function(i) imps_for_tetra_list[[i]] %>% filter(agegr3==3) %>% .[,24:34] %>% tetrachoric(.) %>% .$rho) -> tetrachor_rhos_late## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done## For i = 10 j = 5 A cell entry of 0 was replaced with correct = 0.5. Check your data!#find the number of factors/components we're gonnna reduce down to
lapply(1:length(tetrachor_rhos), function(i) eigen(tetrachor_rhos[[i]]) %>% .$values %>% .[.> 1] %>% length(.)) ->num_factors_tetra_imps
##NOTE: this is '4' factors for nine datasets, and '5' for one (set #8) -- we're gonna code out the dynamic handling and hardcode '4' to build the communalities manually for transparency/defensibility; it's still '4' if we drop 'pacat'-- and also it's still MOSTLY 4 across the age groups (5 a few times w/ 'young')
#perform the actual PCA and extract the H2 (communality as defined by Supplement to: Livingston G, Huntley J, Liu KY, et al. Dementia prevention, intervention, and care: 2024 report of the Lancet standing Commission. Lancet 2024; published online July 31. https://doi.org/10.1016/S0140-6736(24)01296-0),
#see the end of page 9
#the h2 -- communalities -- we're pulling here SHOULD be the same as the 'scaled' ones they build manually through stata, already multiplied by eigenvalues, because the loadings themselves are scaled
## in fact, I'm positive they are: check out the actual code here: https://raw.githubusercontent.com/cran/psych/refs/heads/master/R/principal.R
# loadings <- eigens$vectors %*% sqrt(diag(eigens$values,nrow=length(eigens$values))) #added May 2, 2016 for the weird case of a single variable with covariance > 1
#
# if(nfactors > 0) {loadings <- loadings[,1:nfactors]} else {nfactors <- n}
# if (nfactors > 1) {communalities <- rowSums(loadings^2)} else {communalities <- loadings^2 }
##BUT we're gonna build the communalities manually by pooling the loadings instead of pooling or averaging the communalities -- hence why we're hardcoding the number of components -- because we're dealing with imputed datasets, and this seems more appropriate than just pooling the communalities; I don't want a fuzzy picture of a description of fuzzy pictures, though maybe it'll be the same
#We selected 4 components for all imputations to ensure consistency in pooling and interpretation, given that 90% of imputations retained 4 components based on the eigenvalue ≥ 1 rule.
# produces list of PCAs, one for each imputation
lapply(1:length(tetrachor_rhos), function(i) principal(tetrachor_rhos[[i]], nfactors = 4, rotate = "none", scores = FALSE)) -> tetrachor_PCA
lapply(1:length(tetrachor_rhos_young), function(i) principal(tetrachor_rhos_young[[i]], nfactors = 4, rotate = "none", scores = FALSE)) -> tetrachor_PCA_young
lapply(1:length(tetrachor_rhos_mid), function(i) principal(tetrachor_rhos_mid[[i]], nfactors = 4, rotate = "none", scores = FALSE)) -> tetrachor_PCA_mid
lapply(1:length(tetrachor_rhos_late), function(i) principal(tetrachor_rhos_late[[i]], nfactors = 4, rotate = "none", scores = FALSE)) -> tetrachor_PCA_late
#get the SCALED loading matrices from each imputation/list
loading_matrices<- lapply(tetrachor_PCA, function(x) as.matrix(x$loadings))
loading_matrices_young<- lapply(tetrachor_PCA_young, function(x) as.matrix(x$loadings))
loading_matrices_mid<- lapply(tetrachor_PCA_mid, function(x) as.matrix(x$loadings))
loading_matrices_late<- lapply(tetrachor_PCA_late, function(x) as.matrix(x$loadings))
#slap them together and average/pool them across the values from our imputed datasets
pooled_loadings<- Reduce("+", loading_matrices) / length(loading_matrices)
pooled_loadings_young<- Reduce("+", loading_matrices_young) / length(loading_matrices_young)
pooled_loadings_mid<- Reduce("+", loading_matrices_mid) / length(loading_matrices_mid)
pooled_loadings_late<- Reduce("+", loading_matrices_late) / length(loading_matrices_late)
#build the communalities; go crazy
pooled_communalities <- rowSums(pooled_loadings^2) %>% as.data.frame(.) %>% cbind(., row.names(.)) %>% as.data.frame(.) %>% dplyr::rename(communality=1, var = 2) %>% mutate(var = var %>% str_remove(., "\\_bin"), communality = communality %>% as.numeric(.))
pooled_communalities_young <- rowSums(pooled_loadings_young^2) %>% as.data.frame(.) %>% cbind(., row.names(.)) %>% as.data.frame(.) %>% dplyr::rename(communality=1, var = 2) %>% mutate(var = var %>% str_remove(., "\\_bin"), communality = communality %>% as.numeric(.))
pooled_communalities_mid <- rowSums(pooled_loadings_mid^2) %>% as.data.frame(.) %>% cbind(., row.names(.)) %>% as.data.frame(.) %>% dplyr::rename(communality=1, var = 2) %>% mutate(var = var %>% str_remove(., "\\_bin"), communality = communality %>% as.numeric(.))
pooled_communalities_late <- rowSums(pooled_loadings_late^2) %>% as.data.frame(.) %>% cbind(., row.names(.)) %>% as.data.frame(.) %>% dplyr::rename(communality=1, var = 2) %>% mutate(var = var %>% str_remove(., "\\_bin"), communality = communality %>% as.numeric(.))
##and go ahead and build the "W" to build the adjusted PIF; the "var" for both is the same, I promise, feel free to check
W_df<- {1 - pooled_communalities$communality} %>% as.data.frame(.) %>% cbind(., var = pooled_communalities$var) %>% as.data.frame(.) %>% dplyr::rename(W=1) %>% mutate(W = W %>% as.numeric(.))
# IMPORTANT!!!
# It's worth noting that pooling loadings without sign alignment can underestimate communalities, potentially leading to overoptimistic adjustments to our Potential Impact Fractions (PIF) (i.e. we risk overestimating the adjusted impact of reducing that risk factor). Given our use of multiple imputed datasets to capture data uncertainty, we accepted this limitation as a reasonable trade-off in this context. To account for overly optimistic impact fractions, we introduced a 10% reduction in prevalence rather than the 15% reduction suggested in the existing literature.
#
###### interpreting top communalities:
# ### in essence, the top communalities indicate what should be the highest yield targets because they explain the most "variance" among other risk factors; essentially, they're clustered, and acting on one should help act on others!
# > pooled_communalities %>% arrange(., desc(.$communality))
# communality var
# diabetes_bin 0.7957997 diabetes
# rfdrhv2_bin 0.7872580 rfdrhv2
# c09q02_bin 0.7107840 c09q02
# highbp_bin 0.6725299 highbp
# c07q09_bin 0.6631723 c07q09
# c08q04_bin 0.6547181 c08q04
# c09q01_bin 0.5966956 c09q01
# rfsmok_bin 0.5951654 rfsmok
# bmicat5_bin 0.5739486 bmicat5
# educat3a_bin 0.5281733 educat3a
# cholch1_bin 0.5180080 cholch1
###
### DON'T USE THIS
###
# DYNAMIC FACTOR LOADING AND AUTO COMMUNALITY EXTRACTION W/ VIZ
# lapply(1:length(tetrachor_rhos), function(i) principal(tetrachor_rhos[[i]], nfactors = num_factors_tetra_imps[[i]], rotate = "none", scores = FALSE)) -> tetrachor_PCA_test
#
# lapply(tetrachor_PCA_test, function(x) x$communality) %>% do.call(rbind, .) -> h2s_imps
#
# h2s_imps %>% as.data.frame(.) %>% mutate(id= row.names(.) %>% as.character(.)) %>% reshape2::melt(., id.vars = "id") %>% ggplot(., aes(x=variable, y=value)) + geom_boxplot() + geom_point()
##add in the population values for each group based off of the ACS-5 2023 pop values; here:
library(tidycensus)
rbind(load_variables(2023, "acs5", cache = TRUE) %>% .[,1:3], load_variables(2023, "acs5/profile", cache = TRUE), load_variables(2023, "acs5/subject", cache = TRUE), load_variables(2023, "acs5/cprofile", cache = TRUE)) %>% as.data.frame(.) ->all_acs_vars## Warning: `cache` is deprecated and ignored. tidycensus no longer writes
## variable metadata to a local cache.## Warning: `cache` is deprecated and ignored. tidycensus no longer writes
## variable metadata to a local cache.
## Warning: `cache` is deprecated and ignored. tidycensus no longer writes
## variable metadata to a local cache.
## Warning: `cache` is deprecated and ignored. tidycensus no longer writes
## variable metadata to a local cache.lrgv_23_sexpop <- get_acs(geography = "county", county = c("Hidalgo", "Cameron", "Willacy", "Starr"), table = "B01001", summary_var = "B01001_001", state = "TX", year = 2023, survey = "acs5", geometry = TRUE) %>% left_join(., all_acs_vars, by = c("variable"="name")) %>% separate_wider_delim(., cols = "label", delim = "!!", names_sep="", too_few = "align_start") %>% separate_wider_delim(., cols = "variable", delim = "_", names_sep="", too_few = "align_start") %>% mutate(variable2 = as.numeric(variable2))## Getting data from the 2019-2023 5-year ACS
## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%lrgv_23_sexpop %>% filter(variable2 >= 7 & variable2<=14 | variable2 >= 31 & variable2<=38) %>% filter(grepl("Cameron", NAME)) %>% transform(total_18to44_pop = ave(.$estimate, .$NAME, FUN=sum)) %>% .$total_18to44_pop %>% unique(.) -> cam18to44_23poptotal #149,957
lrgv_23_sexpop %>% filter(variable2 >= 15 & variable2<=19 | variable2 >= 39 & variable2<=43) %>% filter(grepl("Cameron", NAME)) %>% transform(total_45to64_pop = ave(.$estimate, .$NAME, FUN=sum)) %>% .$total_45to64_pop %>% unique(.) -> cam45to64_23poptotal #90,513
lrgv_23_sexpop %>% filter(variable2 >=20 & variable2 <= 25 | variable2 >=44 & variable2 <= 49) %>% filter(grepl("Cameron", NAME)) %>% transform(total_65plus_pop = ave(.$estimate, .$NAME, FUN=sum)) %>% .$total_65plus_pop %>% unique(.) -> cam65plus_23poptotal #59,431
age_margins <- data.frame(
agegr3 = c("1", "2", "3"),
total_pop = c(cam18to44_23poptotal, cam45to64_23poptotal, cam65plus_23poptotal) # real values from ACS/Census
)
######
# the rise of the LHD as front-line defense: modest gain everywhere is significant gain nationally
#
# we're going for a highly conservative estimate of 5% improvement
#
##
target_reduction<- 0.10
#FUNCTION DEFINITON
##calculates the Population Attributable Fraction (PAF)
calc_PAF <- function(P, RR) {
(P * (RR - 1)) / (P * (RR - 1) + 1)
}
## calculates the Potential Impact Fraction (PIF)
calc_PIF <- function(P, P_prime, RR) {
((P - P_prime) * (RR - 1)) / (P * (RR - 1) + 1)
}
## calculates the Adjusted PIF taking the communality into account (AdjPAF)
calc_AdjPIF <- function(PIF, communality) {
(1 - communality) * PIF
}
#join all your values to the key for quality of life; we're dropping "Physical Inactivity" or 'pacat' because it's missing too many values AND entries (not all answer choices represented)
prevs<- left_join(pooled_imp_lancetvars_23_camBRFSS_keyed, age_margins, by=c("agegr3_resp"="agegr3")) %>% mutate(prev = results/total_pop) %>% .[-grepl("pacat", .$var1),] %>% dplyr::select(-key, -lancet_modifiable) %>% mutate(VARIANCE_SUMMED = sum(se^2)) %>% transform(prev = ave(.$prev, .$lancet_modifiable, FUN=sum),
results = ave(.$results, .$lancet_modifiable, FUN=sum),
se = ave(.$VARIANCE_SUMMED, .$lancet_modifiable, FUN=sqrt)) %>% .[,-14] %>% dplyr::rename(., lower=4, upper=5) %>% mutate(upper = results + 2*se, lower = results - 2*se) %>% dplyr::select(-prop, -missInfo, -resp) %>% unique(.)## Adding missing grouping variables: `lancet_modifiable`full_resp_key<- left_join(W_df, pooled_communalities, by=c("var"="var")) %>% left_join(., resp_key %>% dplyr::select(-lancet24_RR, -upperRR, -lowerRR), by=c("var"="var")) %>% left_join(., prevs, by=c("var"="var1")) %>% mutate(prev_prime = prev * (1 - target_reduction), PAF = calc_PAF(prev, lancet24_RR), upperPAF = calc_PAF(prev, upperRR), lowerPAF = calc_PAF(prev,lowerRR), PIF = calc_PIF(prev, prev_prime, lancet24_RR), upperPIF = calc_PIF(prev, prev_prime, upperRR), lowerPIF = calc_PIF(prev, prev_prime, lowerRR), AdjPIF = calc_AdjPIF(PIF, communality), adjPIF_usingW = W*PIF, upperAdjPIF = calc_AdjPIF(upperPIF, communality), lowerAdjPIF = calc_AdjPIF(lowerPIF, communality))
##compare the one using W and the one using the function; they're the same ;)
#homestretch: create the total weighted PIF and GO
total_weighted_PIF <- 1 - prod(1 - full_resp_key$AdjPIF)
total_weighted_upperPIF <- 1 - prod(1 - full_resp_key$upperAdjPIF)
total_weighted_lowerPIF <- 1 - prod(1 - full_resp_key$lowerAdjPIF)
#8800 is the estimated number of dementia cases in those 65 and older in Cameron County in 2020
#source: https://www.alz.org/getmedia/bc3091e7-db20-4293-a919-d2d02f795d2f/tx_prevalence-map-and-spreadsheet.pdf
total_possible_dementiacases_prevented<- 8800 * total_weighted_PIF
upper_possible_dementiacases_prevented<- 8800 * total_weighted_upperPIF
lower_possible_dementiacases_prevented<- 8800 * total_weighted_lowerPIF
print(total_possible_dementiacases_prevented)## [1] 402.4823cat("Max Number Possible:", upper_possible_dementiacases_prevented)## Max Number Possible: 587.6583cat("Lowest Number Possible:", lower_possible_dementiacases_prevented)## Lowest Number Possible: 236.4046## pretty tables
full_resp_key %>% mutate(`Age Group` = case_when(
agegr3_resp == 1 ~ "18 to 44",
agegr3_resp == 2 ~ "45 to 64",
agegr3_resp == 3 ~ "65+",
TRUE ~ NA_character_
), PIF = PIF*100, upperPIF = upperPIF*100, lowerPIF = lowerPIF*100, AdjPIF = AdjPIF*100, upperAdjPIF = upperAdjPIF*100, lowerAdjPIF = lowerAdjPIF*100, PAF = PAF*100, upperPAF = upperPAF*100, lowerPAF = lowerPAF*100, results_ci = sprintf("%.0f (%.0f - %.0f)", results, lower, upper), RR_ci = sprintf("%.1f (%.1f-%.1f)", lancet24_RR, lowerRR, upperRR), PAF_ci = sprintf("%.2f (%.2f - %.2f)", PAF, lowerPAF, upperPAF), PIF_ci = sprintf("%.4f (%.4f - %.4f)", PIF, lowerPIF, upperPIF), AdjPIF_ci = sprintf("%.4f (%.4f - %.4f)", AdjPIF, lowerAdjPIF, upperAdjPIF), attrib = round((1-AdjPIF) * 8800), attribLower = round((1-lowerAdjPIF) * 8800), attribUpper = round((1-upperAdjPIF) * 8800), prop_impact = AdjPIF / sum(AdjPIF), # relative contribution
cases_prevented = prop_impact * total_possible_dementiacases_prevented,
upper_cases_prevented = (upperAdjPIF / sum(upperAdjPIF)) * upper_possible_dementiacases_prevented,
lower_cases_prevented = (lowerAdjPIF / sum(lowerAdjPIF)) * lower_possible_dementiacases_prevented, PossPreventable_ci = sprintf("%.0f (%.0f-%.0f)", cases_prevented, lower_cases_prevented,upper_cases_prevented)) %>% .[,c(5, 30, 3, 28, 29, 15, 31:33, 41)] %>% dplyr::rename(., `Lancet Modifiable Risk Factor`=1, `RR (95% CI)`=2, `Communality, %`=3, `Weighted Estimates (95% CI)`=5, `Age Group Prevalence, %`=6, `PAF, % (95% CI)`=7, `PIF, % (95% CI)`=8, `Adjusted PIF, % (95% CI)`=9, `Possible Number of Dementia Cases Preventable`=10) %>% arrange(., desc(`Age Group Prevalence, %`)) %>% gt(.) %>% opt_stylize(4) %>% cols_align("center") %>% fmt_percent(columns = c(3, 6:9), decimals = 1) %>% fmt_number(columns = c(5), sep_mark = ",", decimals = 0)| Lancet Modifiable Risk Factor | RR (95% CI) | Communality, % | Age Group | Weighted Estimates (95% CI) | Age Group Prevalence, % | PAF, % (95% CI) | PIF, % (95% CI) | Adjusted PIF, % (95% CI) | Possible Number of Dementia Cases Preventable |
|---|---|---|---|---|---|---|---|---|---|
| OBESITY | 1.3 (1.0-1.7) | 57.4% | 45 to 64 | 60246 (41575 - 78918) | 66.6% | 16.64 (0.00 - 31.78) | 1.6645 (0.0000 - 3.1784) | 0.7091 (0.0000 - 1.3541) | 61 (0-116) |
| HIGH CHOLESTEROL | 1.3 (1.3-1.4) | 51.8% | 45 to 64 | 48028 (32562 - 63495) | 53.1% | 13.73 (13.73 - 17.51) | 1.3733 (1.3733 - 1.7509) | 0.6619 (0.6619 - 0.8439) | 57 (58-72) |
| HYPERTENSION | 1.2 (1.1-1.4) | 67.3% | 45 to 64 | 47578 (30305 - 64851) | 52.6% | 9.51 (4.99 - 17.37) | 0.9513 (0.4994 - 1.7373) | 0.3115 (0.1635 - 0.5689) | 27 (14-49) |
| DIABETES | 1.7 (1.6-1.8) | 79.6% | 45 to 64 | 33295 (19378 - 47212) | 36.8% | 20.48 (18.08 - 22.74) | 2.0477 (1.8080 - 2.2737) | 0.4181 (0.3692 - 0.4643) | 36 (32-40) |
| DEPRESSION | 2.2 (1.7-3.0) | 66.3% | 45 to 64 | 27874 (13905 - 41843) | 30.8% | 26.98 (17.73 - 38.12) | 2.6983 (1.7734 - 3.8115) | 0.9089 (0.5973 - 1.2838) | 78 (52-110) |
| SOCIAL ISOLATION | 1.6 (1.3-1.8) | 65.5% | 65+ | 18273 (12396 - 24150) | 30.7% | 15.57 (8.45 - 19.74) | 1.5575 (0.8445 - 1.9742) | 0.5378 (0.2916 - 0.6816) | 46 (25-58) |
| VISION LOSS | 1.5 (1.4-1.6) | 71.1% | 65+ | 15565 (8986 - 22143) | 26.2% | 11.58 (9.48 - 13.58) | 1.1579 (0.9483 - 1.3580) | 0.3349 (0.2743 - 0.3928) | 29 (24-34) |
| LESS EDUCATION | 1.6 (1.3-2.0) | 52.8% | 18 to 44 | 28357 (11445 - 45268) | 18.9% | 10.19 (5.37 - 15.90) | 1.0190 (0.5368 - 1.5903) | 0.4808 (0.2533 - 0.7503) | 41 (22-64) |
| SMOKING | 1.3 (1.2-1.4) | 59.5% | 45 to 64 | 12176 (4174 - 20177) | 13.5% | 3.88 (2.62 - 5.11) | 0.3879 (0.2620 - 0.5106) | 0.1570 (0.1061 - 0.2067) | 14 (9-18) |
| EXCESSIVE ALCOHOL | 1.2 (1.0-1.5) | 78.7% | 45 to 64 | 10837 (2020 - 19655) | 12.0% | 2.34 (0.00 - 5.65) | 0.2339 (0.0000 - 0.5649) | 0.0498 (0.0000 - 0.1202) | 4 (0-10) |
| HEARING LOSS | 1.4 (1.0-1.9) | 59.7% | 45 to 64 | 5646 (732 - 10559) | 6.2% | 2.43 (0.00 - 5.32) | 0.2434 (0.0000 - 0.5315) | 0.0982 (0.0000 - 0.2144) | 8 (0-18) |
### TO-DO: MARKOV CHAIN MODELINGAdd a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.