## full dataset = data_new
library(haven)
library(dplyr)
library(ggplot2)
library(broom)
library(caTools)
library(pROC)
library(ROCR)
library(patchwork)
install.packages("gtsummary")
install.packages("flextable")
install.packages("officer")
install.packages("broom")
install.packages("knitr")
library(broom)
library(knitr)
library(officer)
```
head(data_new)
tail(data_new)
str(data_new)
library(readxl)
library(tidyr)
install.packages("tidyr", tpye = "binary")
install.packages("readr", tpye = "binary")
library(dplyr)


# Keep only rows in b that have PAT_ID also in a
ToothLossModelingPopulation_09112025_filtered <- semi_join(ToothLossModelingPopulation_09112025, ToothLossModelingPopulation_07102025_1_, by = "PAT_ID")

# Check dimensions
dim(ToothLossModelingPopulation_09112025_filtered)

# IDs in a but not in b
missing_ids <- setdiff(ToothLossModelingPopulation_07102025_1_$PAT_ID, ToothLossModelingPopulation_09112025$PAT_ID)

length(missing_ids)  # how many are missing
head(missing_ids)    # look at first few
# check for duplicates
sum(duplicated(ToothLossModelingPopulation_07102025_1_$PAT_ID))  # duplicates in a
sum(duplicated(ToothLossModelingPopulation_09112025$PAT_ID))  # duplicates in b

# Identify the missing patients
missing_ids <- setdiff(ToothLossModelingPopulation_07102025_1_$PAT_ID, ToothLossModelingPopulation_09112025$PAT_ID)

# Subset those patients from dataset a
ToothLossModelingPopulation_07102025_1_missing <- ToothLossModelingPopulation_07102025_1_[ToothLossModelingPopulation_07102025_1_$PAT_ID %in% missing_ids, ]

#combine a & b
# First filter b to keep only patients also in a
ToothLossModelingPopulation_09112025_filtered <- semi_join(ToothLossModelingPopulation_09112025, ToothLossModelingPopulation_07102025_1_, by = "PAT_ID")

# make sure keys are characters
ToothLossModelingPopulation_07102025_1_$PAT_ID <- as.character(ToothLossModelingPopulation_07102025_1_$PAT_ID)
ToothLossModelingPopulation_09112025$PAT_ID <- as.character(ToothLossModelingPopulation_09112025$PAT_ID)

# ZIPs should be character (not numeric)
ToothLossModelingPopulation_07102025_1_$ZIP_CODE <- as.character(ToothLossModelingPopulation_07102025_1_$ZIP_CODE)
ToothLossModelingPopulation_09112025$ZIP_CODE <- as.character(ToothLossModelingPopulation_09112025$ZIP_CODE)

# IDs in a but not in b
missing_ids <- setdiff(ToothLossModelingPopulation_07102025_1_$PAT_ID, ToothLossModelingPopulation_09112025$PAT_ID)

ToothLossModelingPopulation_07102025_1_missing   <- ToothLossModelingPopulation_07102025_1_ %>% filter(PAT_ID %in% missing_ids)
ToothLossModelingPopulation_09112025_filtered  <- semi_join(ToothLossModelingPopulation_09112025, ToothLossModelingPopulation_07102025_1_, by = "PAT_ID")

# ensure both have the same columns and types before binding
ToothLossModelingPopulation_07102025_1_missing_for_ToothLossModelingPopulation_09112025 <- ToothLossModelingPopulation_07102025_1_missing %>% select(any_of(names(ToothLossModelingPopulation_09112025_filtered)))
ToothLossModelingPopulation_09112025_completed     <- bind_rows(ToothLossModelingPopulation_09112025_filtered, ToothLossModelingPopulation_07102025_1_missing_for_ToothLossModelingPopulation_09112025)
nrow(ToothLossModelingPopulation_09112025_completed)   # should be 14,683

data_new <- ToothLossModelingPopulation_09112025_completed

# recode for area-based poverty level
data_new <- ToothLossModelingPopulation_07102025_1_
data_new$abpl <- NULL
data_new$ABPLnum <- NULL
data_new$ABPL <- NULL
data_new$Zipcode <- data_new$ZIP_CODE
data_new$ABPL[data_new$ZIP_CODE %in% c("10954", "10024", "10023","07043","10536","10011", "10021","10282","10803","11791","10028","41042","11375","12732","10128","10543","08553",
                                       "10708","07052","10471","11215","10005","07006","10309","11520","10703","11518","10603","10604","12763","10014","10541","11415","10016","10014","07649","11969","10314",
                                       "07733","19454","10312","13492","07030","11412","10553","10552","11109","11003","10008","10113","12518","11742","06883","10591","10523","10583","10010","11357","10065",
                                       "11234","10805","07642","11356","10579","34711","08690","07666","07621","12590","11710","11429","11361","07632","11516","12533","07024","10970","11552","10562","10994","07024",
                                       "10538","33647","10930","11379","19605","07670","10804","11580","10022","10709","10989","03768","11379","10528","10590","07410","11716","10980","11776","11782",
                                       "10528","10305","10017","11954","10710","10965","11727","06830","07020","17111","11422","10580","10516","10530","07656","11946","06461","01027","11363",
                                       "11222","07605","12553","11557","12302","10917","10007","10962","10598","12603","10706","11598","07079","10307","07110","11360","94619",
                                       "22310","11507","11590","10549","12528","07646","12520","10982","06811","10918","11793","11561","07860","11358","07645","11510","11717",
                                       "10707","11787","11411","10464","11717","11040","10901","10920","33029","11364","07603","08530","07663","20782","11768","10986",
                                       "07628","10522","10533","07838","11362","11004","12775","10547","08060","11753","07641","10924","07028","04107","32903","10567","10514",
                                       "12570","11563","11581","11963","06776","10960","68507","22508","27513","11414","07920","11545","12582","10921","10510","11746",
                                       "11756","10597","12508","33470","07031","11023","10308","08831","08837","11566","11783","10520","10956","07003","07620","07661","11731",
                                       "11783","11530","10983","08759","07648","07 643","93720","07607","07960","10992","33131","07640","07086","11777","07647","87508",
                                       "07094","07624","07012","11021","06807","10607","10006","12473","06880","12051","12205","13413","07071","11706","11413","07981","06840","08520",
                                       "11050","12549","11559","18235","17025","12524","98102","07065","10512","14221","18901","28273","11714","95037","06905","12548","16112","06824",
                                       "11572","12766","07013","11570","11968","12572","10156","11509","06612")] <- "low"
data_new$ABPL[data_new$ZIP_CODE %in% c("10701", "10025", "10704", "10034", "10040","10033", "10463","10036", "11368", "10301","11101","11421","11106","10461","10001","11421","11105","11374","10019","11102","11230","11436","10705","10462","
                                       11367","11373","32129","11372","11210","12723","11691","11236","27406","10566","07631","11201","10470","06902","14883","10306","12491",
                                       "10993","12401","10465","11423","11238","11435","11218","10702","10475","11103","11693","11238","11435","11226","10601","10475","10469","11434","10044","11385","11416","11209","11370","11216","11204",
                                       "11550","11202","10310","21740","11377","10310","07010","11417","10012","11104","11354","11432","11203","11978","07307","12601","07601","33896","10018","10075","10003",
                                       "11229","07650","07660","11225","11231","07601","10940","11203","11235","07047","11419","32666","11217","92614","11428","11369","12563",
                                       "10013","03431","11433","10801","10303","11937","11232","10927","10118","12586","12721","10304","12498","02657","12538","33056","10941",
                                       "12487","10573","12534","10573","11418","97850","13468","11378","12208","11365","11951","10595","12440","11427","11420","11427",
                                       "12303","11365","11420","10987","07022","77019","07306","11553","11722","12446","10302","12203","12550","13126","33442","07514","12701",
                                       "33138","30213","29486","10925","07029","12758","07420","07407","07712","06810","17512","20002","12602","10101","32807","11694","22172",
                                       "10276","18103","10923","11228","07008","11228","11694","28371","06451","06901","07644","10606","11247","11366","95834","18853",
                                       "33637","11120","07305","12404","07305")] <- "medium"
data_new$ABPL[data_new$ZIP_CODE %in% c("10032", "10039","10031", "10037","10458","10027","10467","10468","11221","10038","11213","10026","11207","10002","10950","11249","11220",
                                       "10455","10009","11208","11355","11239","10466","28311","10473","11211","11223","11237","10550","11233","11224","11214","10069","11692","07093",
                                       "07304","07093","33311","08105","10977","75220","13502","10949","07522","12592","06511","06515","12202","14213","07087","27405",
                                       "07106","07104","06604","02905","14760","18202","12461","07206","06704","14850","33060","30344","08010","07055")] <- "high"
data_new$ABPL[data_new$ZIP_CODE %in% c("10029", "10030", "10453", "10452","11219","10457","10459", "10035",
                                       "10460","10451","10454","11212","10456","11205","10472","12201","11206","10474","07114","10952","12733","07108","07102","13208",
                                       "10545")] <- "very high"

## numerize abpl
## 0 = low, 1 = medium, 2 = high, 3 = very high
data_new$ABPLnum[data_new$ABPL == "low"] = "0"
data_new$ABPLnum[data_new$ABPL == "medium"] = "1"
data_new$ABPLnum[data_new$ABPL == "high"] = "2"
data_new$ABPLnum[data_new$ABPL == "very high"] = "3"
data_new$ABPLnum <- factor(data_new$ABPLnum, levels = c("0", "1", "2", "3"))
str(data_new$ABPLnum)

##recode Rheumatoid arthritis
## 0 = no, 1 = yes
data_new$ARTHRITIS[data_new$ARTHRITIS == "no"] = "0"
data_new$ARTHRITIS[data_new$ARTHRITIS == "yes"] = "1"
data_new$ARTHRITIS <- as.numeric(as.character(data_new$ARTHRITIS))
str(data_new$ARTHRITIS)

#recode depression
## 0 = no, 1 = yes
str(data_new$DEPRESSION)

##recode ethnicity
data_new$Ethnicity_race_recode[data_new$ETHNICITY == "HISPANIC OR LATINO OR SPANISH ORIGIN"] = "Hispanic"
data_new$Ethnicity_race_recode[data_new$ETHNICITY == "NOT HISPANIC OR LATINO OR SPANISH ORIGIN" & data_new$RACE == "WHITE"] = "Non-Hispanic white"
data_new$Ethnicity_race_recode[data_new$ETHNICITY == "NOT HISPANIC OR LATINO OR SPANISH ORIGIN" & data_new$RACE == "BLACK OR AFRICAN AMERICAN"] = "Non-Hispanic black"
data_new$ABPLnum <- as.numeric(as.character(data_new$ABPLnum))

## numerize ethnicity recode
## 0 = white, 1 = black, 2 = hispanic
data_new$Ethnicity_race_recodeNum[data_new$Ethnicity_race_recode == "Non-Hispanic white"] = "0"
data_new$Ethnicity_race_recodeNum[data_new$Ethnicity_race_recode == "Non-Hispanic black"] = "1"
data_new$Ethnicity_race_recodeNum[data_new$Ethnicity_race_recode == "Hispanic"] = "2"
data_new$Ethnicity_race_recodeNum <- factor(data_new$Ethnicity_race_recodeNum, levels = c("0", "1", "2"))

## recode functional dentition
## 0 = 21 or Over, 1 = 20 or Less
data_new$TEETH_BUCKET_binary[data_new$TEETH_BUCKET == "21 or Over"] = "0"
data_new$TEETH_BUCKET_binary[data_new$TEETH_BUCKET == "20 or Less"] = "1"
data_new$TEETH_BUCKET_binary <- as.numeric(as.character(data_new$TEETH_BUCKET_binary))

## recode gender 
## 0 = male, 1 = female
data_new$GENDER_binary[data_new$GENDER == "Male"] = "0"
data_new$GENDER_binary[data_new$GENDER == "Female"] = "1"
data_new$GENDER_binary <- as.numeric(as.character(data_new$GENDER_binary))

## recode employment status
## 0 = Full Time, 1 = Part Time, 2 = unemployed & disabled, 3 = Retired
data_new$EMPLOYMENT_STATUS_cat <- NULL
data_new$EMPLOYMENT_STATUS_cat[data_new$EMPLOYMENT_STATUS == "Full Time"] = "0"
data_new$EMPLOYMENT_STATUS_cat[data_new$EMPLOYMENT_STATUS == "Part Time"] = "1"
data_new$EMPLOYMENT_STATUS_cat[data_new$EMPLOYMENT_STATUS == "Not Employed"] = "2"
data_new$EMPLOYMENT_STATUS_cat[data_new$EMPLOYMENT_STATUS == "Not Disabled"] = "2"
data_new$EMPLOYMENT_STATUS_cat[data_new$EMPLOYMENT_STATUS == "Retired"] = "3"
data_new$EMPLOYMENT_STATUS_cat <- factor(data_new$EMPLOYMENT_STATUS_cat, levels = c("0", "1", "2", "3"))

## recode interpreter needs
data_new$NEEDS_INTERPRETER_binary[data_new$NEEDS_INTERPRETER == "N"] = "0"
data_new$NEEDS_INTERPRETER_binary[data_new$NEEDS_INTERPRETER == "Y"] = "1"
str(data_new$NEEDS_INTERPRETER_binary)
data_new$NEEDS_INTERPRETER_binary <- factor(data_new$NEEDS_INTERPRETER_binary, levels = c("0", "1"))

## dental insurance recode
## 0 = self pay, 1 = Medicaid, 2 = Medicare
data_new$INSURANCE_cat[data_new$DENTAL_FC == "Self-Pay"] = "0"
data_new$INSURANCE_cat[data_new$DENTAL_FC == "Medicaid"] = "1"
data_new$INSURANCE_cat[data_new$DENTAL_FC == "Managed Medicaid"] = "1"
data_new$INSURANCE_cat[data_new$DENTAL_FC == "Managed Medicare"] = "2"
data_new$INSURANCE_cat[data_new$DENTAL_FC == "Medicare"] = "2"
str(data_new$INSURANCE_cat)
data_new$INSURANCE_cat <- factor(data_new$INSURANCE_cat, levels = c("0", "1", "2"))

## tobacco use recode
## 0 = never, 1 = smoker at some point, . = never assessed/unknown
data_new$SMOKING_binary[data_new$SMOKING_TOB_USE == "Never"] = "0"
data_new$SMOKING_binary[data_new$SMOKING_TOB_USE == "Every Day"] = "1"
data_new$SMOKING_binary[data_new$SMOKING_TOB_USE == "Some Days"] = "1"
data_new$SMOKING_binary[data_new$SMOKING_TOB_USE == "Former"] = "1"
data_new$SMOKING_binary[data_new$SMOKING_TOB_USE == "Light Smoker"] = "1"
str(data_new$SMOKING_binary)
data_new$SMOKING_binary <- factor(data_new$SMOKING_binary, levels = c("0", "1"))

## HIV
str(data_new$HIV)
## Diabetes
str(data_new$DIABETES)

## 3. age cat with all
data_new$AGE_ALL <- ifelse(data_new$AGE_LAST_VISIT >= 30, 
                           data_new$AGE_LAST_VISIT, 
                           NA)
## categorized - 0 = <40, 1 = 40-49, 2 = 50-59, 3 = 60-69, 4 = 70-79, 5 = 79+
data_new$AGE_ALL_CAT <- cut(
  data_new$AGE_ALL,
  breaks = c(29, 39, 49, 59, 69, 79, Inf),  # defines the intervals
  labels = c(0, 1, 2, 3, 4, 5),          # assigns your custom labels
  right = TRUE                     # upper limit is inclusive
)
table(data_new$AGE_ALL_CAT, useNA = "ifany")


## Clinical attachment loss recode

## 1. number of CAL with 1-2mm
data_new$NUM_CAL_1_2 <- data_new$NUM_SITES_CAL_1 + data_new$NUM_SITES_CAL_2
summary(data_new$NUM_CAL_1_2)

## Proportion of 1_2 >30% total probed
data_new$Proportion_CAL_1_2 = data_new$NUM_CAL_1_2/data_new$NUM_TOTAL_SITES_PROBED
summary(data_new$Proportion_CAL_1_2)

## 2. number of probing sites with 3-4mm
data_new$NUM_CAL_3_4 <- data_new$NUM_SITES_CAL_3 + data_new$NUM_SITES_CAL_4
summary(data_new$NUM_CAL_3_4)
## Proportion of 3_4 >30% total probed
data_new$Proportion_CAL_3_4 = data_new$NUM_CAL_3_4/data_new$NUM_TOTAL_SITES_PROBED
summary(data_new$Proportion_CAL_3_4)

## 3. combine 1&2 (anything below 5mm)
data_new$NUM_CAL_less_5 <- data_new$NUM_SITES_CAL_1 + data_new$NUM_SITES_CAL_2 + data_new$NUM_SITES_CAL_3 + data_new$NUM_SITES_CAL_4
summary(data_new$NUM_CAL_less_5)
data_new$Proportion_CAL_less_5 = data_new$NUM_CAL_less_5/data_new$NUM_TOTAL_SITES_PROBED
summary(data_new$Proportion_CAL_less_5)

## 4. number of probing sites with 5&+mm
data_new$NUM_CAL_5_PLUS <- data_new$NUM_SITES_CAL_5 + data_new$NUM_SITES_CAL_6_PLUS
summary(data_new$NUM_CAL_5_PLUS)
## Proportion of 5&+ >30% total probed
data_new$Proportion_CAL_5_PLUS = data_new$NUM_CAL_5_PLUS/data_new$NUM_TOTAL_SITES_PROBED
summary(data_new$Proportion_CAL_5_PLUS)

## reflect the highest level of disease
## 0 = highest 1-2mm, 1 = highest 3-4mm, 2 = highest 5+mm
data_new$Highest_PD_Category <- with(data_new, ifelse(
  data_new$NUM_CAL_5_PLUS > 0, 2,
  ifelse(data_new$NUM_CAL_3_4 > 0, 1, 0)
))
table(data_new$Highest_PD_Category, useNA = "ifany")

data_new$Highest_PD_Category <- factor(
  data_new$Highest_PD_Category,
  levels = c(0, 1, 2),
  ordered = TRUE
)

## Presence of generalized chronic periodontitis - <30%
data_new$Generalized_Perio_1_2_binary <- ifelse(
  is.na(data_new$Proportion_CAL_1_2) | data_new$Proportion_CAL_1_2 < 0.3,
  0,
  1
)
table(data_new$Generalized_Perio_1_2_binary)

data_new$Generalized_Perio_3_4_binary <- ifelse(
  is.na(data_new$Proportion_CAL_3_4) | data_new$Proportion_CAL_3_4 < 0.3,
  0,
  1
)
table(data_new$Generalized_Perio_3_4_binary)

data_new$Generalized_Perio_5_PLUS_binary <- ifelse(
  is.na(data_new$NUM_CAL_5_PLUS) | data_new$Proportion_CAL_5_PLUS < 0.3,
  0,
  1
)
table(data_new$Generalized_Perio_5_PLUS_binary)

## CAL ordinal
## 0 = highest 1-2mm, 1 = highest 3-4mm, 2 = highest 5+mm
data_new$Highest_Generalized_Perio_Category_CAL <- with(data_new, ifelse(
  Generalized_Perio_5_PLUS_binary > 0, 2,
  ifelse(data_new$Generalized_Perio_3_4_binary > 0, 1, 0)
))
table(data_new$Highest_Generalized_Perio_Category_CAL, useNA = "ifany")

data_new$Highest_Generalized_Perio_Category_CAL <- 
  as.character(data_new$Highest_Generalized_Perio_Category_CAL)

# Replace NA with "0"
data_new$Highest_Generalized_Perio_Category_CAL[is.na(data_new$Highest_Generalized_Perio_Category_CAL)] <- "0"

# Convert back to factor (optional, keeps levels clean)
data_new$Highest_Generalized_Perio_Category_CAL <- factor(
  data_new$Highest_Generalized_Perio_Category_CAL,
  levels = c("0", "1", "2")
)
## 0=13477, 1=477, 2=729 


## 3/4 & 5+ combined
data_new$Generalized_Perio_Binary_CAL <- with(data_new, ifelse(
  Generalized_Perio_5_PLUS_binary > 0, 1,
  ifelse(data_new$Generalized_Perio_3_4_binary > 0, 1, 0)
))
table(data_new$Generalized_Perio_Binary_CAL, useNA = "ifany")

# Replace NA with "0"
data_new$Generalized_Perio_Binary_CAL[is.na(data_new$Generalized_Perio_Binary_CAL)] <- "0"

# Convert back to factor (optional, keeps levels clean)
data_new$Generalized_Perio_Binary_CAL <- factor(
  data_new$Generalized_Perio_Binary_CAL,
  levels = c("0", "1")
)
## 0 = 13477, 1 = 1206 

table(data_new$TEETH_BUCKET_binary) ## 21 or Over = 0 (11217), 20 or Less = 1 (3466)

#cross tabs
## 1. functional dentition
xtabs(~ Ethnicity_race_recode + TEETH_BUCKET_binary, data=data_new)
## With row percentages - "What % of people within each ethnicity group have/do not have functional dentition?"
prop.table(xtabs(~ Ethnicity_race_recode + TEETH_BUCKET_binary, data=data_new), margin = 1)*100

## distribution with function dentition vs. without
data_new %>%
  
  table_age <- table(data_new$AGE_ALL_CAT, data_new$TEETH_BUCKET_binary)
prop.table(table_age, margin = 2)  # column-wise %
table_gender <- table(data_new$GENDER_binary, data_new$TEETH_BUCKET_binary)
prop.table(table_gender, margin = 2)
table_marital <- table(data_new$MARITAL_STATUS_binary, data_new$TEETH_BUCKET_binary)
prop.table(table_marital, margin = 2)
xtabs(~ AGE_ALL_CAT + Highest_Generalized_Perio_Category_CAL, data=data_new)
tab1 <- xtabs(~ AGE_ALL_CAT + Highest_Generalized_Perio_Category_CAL, data = data_new)
row_pct <- prop.table(tab1, margin = 1) * 100
counts <- as.data.frame(tab1)
percent <- as.data.frame(row_pct)
merge(counts, percent, by = c("AGE_ALL_CAT", "Highest_Generalized_Perio_Category_CAL"),
      suffixes = c("_count", "_pct"))

xtabs(~ AGE_ALL_CAT + TEETH_BUCKET_binary, data=data_new) 
tab2 <- xtabs(~ AGE_ALL_CAT + TEETH_BUCKET_binary, data = data_new)
row_pct <- prop.table(tab2, margin = 1) * 100
counts <- as.data.frame(tab2)
percent <- as.data.frame(row_pct)
merge(counts, percent, by = c("AGE_ALL_CAT", "TEETH_BUCKET_binary"),
      suffixes = c("_count", "_pct"))

xtabs(~ GENDER_binary + Highest_Generalized_Perio_Category_CAL, data=data_new)
tab3 <- xtabs(~ GENDER_binary + Highest_Generalized_Perio_Category_CAL, data = data_new)
row_pct <- prop.table(tab3, margin = 1) * 100
counts <- as.data.frame(tab3)
percent <- as.data.frame(row_pct)
merge(counts, percent, by = c("GENDER_binary", "Highest_Generalized_Perio_Category_CAL"),
      suffixes = c("_count", "_pct"))

xtabs(~ GENDER_binary + TEETH_BUCKET_binary, data=data_new)
tab4 <- xtabs(~ GENDER_binary + TEETH_BUCKET_binary, data = data_new)
row_pct <- prop.table(tab4, margin = 1) * 100
counts <- as.data.frame(tab4)
percent <- as.data.frame(row_pct)
merge(counts, percent, by = c("GENDER_binary", "TEETH_BUCKET_binary"),
      suffixes = c("_count", "_pct"))

xtabs(~ GENDER_binary + TEETH_BUCKET_binary, data=data_new)
xtabs(~ AGE_ALL_CAT + Highest_Generalized_Perio_Category, data=data_new)

table_hiv <- table(data_new$HIV, data_new$TEETH_BUCKET_binary)
prop.table(table_hiv, margin = 2)
table_diabetes <- table(data_new$DIABETES, data_new$TEETH_BUCKET_binary)
prop.table(table_diabetes, margin = 2)
table_smoking <- table(data_new$SMOKING_binary, data_new$TEETH_BUCKET_binary)
prop.table(table_smoking, margin = 2)
table_insurance <- table(data_new$INSURANCE_cat, data_new$TEETH_BUCKET_binary)
prop.table(table_insurance, margin = 2)
table_depression <- table(data_new$DEPRESSION, data_new$TEETH_BUCKET_binary)
prop.table(table_depression, margin = 2)
table_arthritis <- table(data_new$ARTHRITIS, data_new$TEETH_BUCKET_binary)
prop.table(table_arthritis, margin = 2)
table_ethnicity <- table(data_new$Ethnicity_race_recodeNum, data_new$TEETH_BUCKET_binary)
prop.table(table_ethnicity, margin = 2)
chisq.test(table_age)
chisq.test(table_gender)
chisq.test(table_ethnicity)
chisq.test(table_hiv)
chisq.test(table_diabetes)
chisq.test(table_smoking)
chisq.test(table_insurance)
chisq.test(table_depression)
chisq.test(table_arthritis)
table_genperioCAL <- table(data_new$Highest_Generalized_Perio_Category_CAL, data_new$TEETH_BUCKET_binary)
prop.table(table_genperioCAL, margin = 2)
chisq.test(table_genperioCAL)
table_abpl <- table(data_new$ABPLnum, data_new$TEETH_BUCKET_binary)
prop.table(table_abpl, margin = 2)  # column-wise %
table_interpreter <- table(data_new$NEEDS_INTERPRETER_binary, data_new$TEETH_BUCKET_binary)
prop.table(table_interpreter, margin = 2)
table_employment <- table(data_new$EMPLOYMENT_STATUS_cat, data_new$TEETH_BUCKET_binary)
prop.table(table_employment, margin = 2)
chisq.test(table_abpl)
chisq.test(table_interpreter)
chisq.test(table_employment)

## generate table functional dentition
make_prevalence_table <- function(data, group_var, outcome_var = "TEETH_BUCKET_binary") {
  # Create table of counts
  tab <- table(data[[group_var]], data[[outcome_var]])
  tab_df <- as.data.frame.matrix(tab)
  colnames(tab_df) <- c("Functional Dentition", "Functional Tooth Loss")
  tab_df$Category <- rownames(tab_df)
  tab_df$Total <- rowSums(tab_df[,1:2])
  
  # Calculate percentages
  tab_df$`Functional Dentition %` <- round(100 * tab_df$`Functional Dentition` / tab_df$Total, 1)
  tab_df$`Functional Tooth Loss %` <- round(100 * tab_df$`Functional Tooth Loss` / tab_df$Total, 1)
  
  # Reorder and format
  tab_df <- tab_df %>%
    select(Category, `Functional Dentition`, `Functional Dentition %`,
           `Functional Tooth Loss`, `Functional Tooth Loss %`, Total)
  
  # Add p-value
  pval <- tryCatch({
    chisq.test(table(data[[group_var]], data[[outcome_var]]))$p.value
  }, error = function(e) NA)
  
  tab_df$p_value <- ifelse(1:nrow(tab_df) == 1, round(pval, 4), "")
  return(tab_df)
}
vars <- c("AGE_ALL_CAT", "GENDER_binary", "Ethnicity_race_recodeNum", 
          "HIV", "DIABETES", "ARTHRITIS", "DEPRESSION", "SMOKING_binary", "INSURANCE_cat", "ABPLnum", "NEEDS_INTERPRETER_binary", "EMPLOYMENT_STATUS_cat", "Highest_Generalized_Perio_Category_CAL")

table_list <- lapply(vars, function(v) {
  df <- make_prevalence_table(data_new, v)
  df$Variable <- v
  df
})

# Combine all tables
prevalence_table <- bind_rows(table_list) %>%
  relocate(Variable, .before = Category)

write.csv(prevalence_table, "Functional_Tooth_Loss_Prevalence_Table.csv", row.names = FALSE)

count_percent <- function(var, group) {
  tab <- table(var, group)
  pct <- prop.table(tab, margin = 2) * 100  # column-wise %
  result <- matrix(paste0(tab, " (", round(pct, 1), "%)"), 
                   nrow = nrow(tab), 
                   dimnames = dimnames(tab))
  return(result)
}
count_percent(data_new$AGE_50_OLDER_CAT, data_new$TEETH_BUCKET_binary)
count_percent(data_new$GENDER_binary, data_new$TEETH_BUCKET_binary)
count_percent(data_new$Ethnicity_race_recodeNum, data_new$TEETH_BUCKET_binary)
count_percent(data_new$HIV, data_new$TEETH_BUCKET_binary)
count_percent(data_new$DIABETES, data_new$TEETH_BUCKET_binary)
count_percent(data_new$SMOKING_binary, data_new$TEETH_BUCKET_binary)
count_percent(data_new$INSURANCE_cat, data_new$TEETH_BUCKET_binary)
count_percent(data_new$ABPLnum, data_new$TEETH_BUCKET_binary)
count_percent(data_new$NEEDS_INTERPRETER_binary, data_new$TEETH_BUCKET_binary)
count_percent(data_new$EMPLOYMENT_STATUS_cat, data_new$TEETH_BUCKET_binary)
count_percent(data_new$Highest_Generalized_Perio_Category_CAL, data_new$TEETH_BUCKET_binary)
table_genperioCAL <- table(data_new$Highest_Generalized_Perio_Category_CAL, data_new$TEETH_BUCKET_binary)
prop.table(table_genperioCAL, margin = 2)
xtabs(~ Highest_Generalized_Perio_Category_CAL + TEETH_BUCKET_binary, data=data_new)

table_genperioCAL <- table(data_new$Highest_Generalized_Perio_Category_CAL, data_new$TEETH_BUCKET_binary)
prop.table(table_genperioCAL, margin = 2)
xtabs(~ Highest_Generalized_Perio_Category_CAL + TEETH_BUCKET_binary, data=data_new)

chisq.test(table_age)
chisq.test(table_gender)
chisq.test(table_ethnicity)
chisq.test(table_hiv)
chisq.test(table_diabetes)
chisq.test(table_smoking)
chisq.test(table_insurance)
table(data_new$Highest_Generalized_Perio_Category_CAL)

##interaction

model_interaction <- glm(
  TEETH_BUCKET_binary ~ AGE_ALL_CAT * INSURANCE_cat + AGE_ALL_CAT * EMPLOYMENT_STATUS_cat,
  data = data_new,
  family = binomial
)

summary(model_interaction)

model_interaction2 <- glm(
  TEETH_BUCKET_binary ~ EMPLOYMENT_STATUS_cat * INSURANCE_cat,
  data = data_new,
  family = binomial
)

summary(model_interaction2)


#logistic regression - functional dentition
## model - age, sex & ethnicity/race, HIV, diabetes, smoking, dental insurance, ABPL, EMPLOYMENT, interpreter
logistic_model_122225 <- glm(data_new$TEETH_BUCKET_binary ~ data_new$AGE_ALL_CAT + data_new$GENDER_binary + data_new$Ethnicity_race_recodeNum + 
                               data_new$HIV + data_new$DIABETES + data_new$SMOKING_binary + data_new$INSURANCE_cat + data_new$ABPLnum + data_new$EMPLOYMENT_STATUS_cat + data_new$NEEDS_INTERPRETER_binary + 
                               data_new$ARTHRITIS + data_new$DEPRESSION,
                             data = data_new, 
                             family = binomial)
summary(logistic_model_122225)
exp(confint(logistic_model_122225))
tidy(logistic_model_122225,exponentiate = TRUE, conf.int = TRUE) %>%
  print(n = Inf)
nobs(logistic_model_122225)

## adjusted model with generalized periodontitis
logistic_model_122225_adjustedforperio <- glm(data_new$TEETH_BUCKET_binary ~ data_new$AGE_ALL_CAT + data_new$GENDER_binary + data_new$Ethnicity_race_recodeNum + 
                                                data_new$HIV + data_new$DIABETES + data_new$SMOKING_binary + data_new$INSURANCE_cat + data_new$ABPLnum + data_new$EMPLOYMENT_STATUS_cat + data_new$NEEDS_INTERPRETER_binary + data_new$ARTHRITIS + data_new$DEPRESSION +
                                                data_new$Highest_Generalized_Perio_Category_CAL,
                                              data = data_new, 
                                              family = binomial)
summary(logistic_model_122225_adjustedforperio)
exp(confint(logistic_model_122225_adjustedforperio))
tidy(logistic_model_122225_adjustedforperio,exponentiate = TRUE, conf.int = TRUE) %>%
  print(n = Inf)

nobs(logistic_model_122225_adjustedforperio)

## variance (deviance/model variability)
summary(logistic_model_122225)$deviance
summary(logistic_model_122225_adjustedforperio)$deviance 
nobs(logistic_model_122225)
nobs(logistic_model_122225_adjustedforperio)
##likelihood ratio test
anova(logistic_model_122225, logistic_model_122225_adjustedforperio, test = "Chisq")