The following code will show you how I created functions to run my multinomial regressions. This page includes:

Create function to run univariate multinomial regressions.
Run block 1 of the step wise multivariate multinomial regression.
Run block 2 of the step wise multivariate multinomial regression.
Create function for block 3 of step wise multivariate multinomial regression.
Obtain R^2 values for each varaible in regressio model, and R^2 change between blocks of the stepwise regression
Run function for block 3 of the step wise multivariate multinomial regression.
Create table of output from stepwise regression model.

Create Data frames

This first set of code is used to create the dataframes needed for our analysis. Once the data frame is made, then we can use the function below to add values to the data frames.

variable    = 0 
  statistic = 0
  estimate  = 0
  conf.low  = 0
  conf.high = 0
  ci        = 0
  p.value   = 0
  r2_change = 0
  
# Univariate Multinomial 

  Ins_Mod_Uni <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Ins_Act_Uni <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Act_Mod_Uni <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  
# Block 1 Covars - SF8

  Ins_Mod_Cov <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Ins_Act_Cov <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Act_Mod_Cov <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  
# Block 2 Covars + SF8
  
  Ins_Mod_SF8 <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Ins_Act_SF8 <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Act_Mod_SF8 <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]

# Block 3 Full Model
  Ins_Mod  <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Ins_Act  <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]
  Act_Mod  <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))[1,]

Preliminary Regression: Univariate Regression function

The first step in my analysis was to run the multinomial regressions without covariates so I can see how the health conditions are associated with physical activity level independent of covariates. The following UniThesis() function does the following:

Runs the multinomial models using multinom()
Manually calculates the coefficients B, the odds ratios (OR) + 95% Confideince interval (95% CI), and p values
Extracts respective values for each type of comparison (e.g., sufficiently active vs moderatly active)
Returns respective values as a row in a dataframe which can be compiled into a larger data fram per comparison

Regression Code

Create `UniThesis()` Function

# Univariate Logistic Multinomial Regression function

UniThesis<-function(Data,condition,string,contrast){

#1. run model 1
  df <- multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ condition, data = Data)

#2. run model 2
  df2 <- multinom(relevel(Ex_3cat,ref = "Moderate") ~ condition, data = Data)

#4 Coefficients 
  coefs_df <- abs(coef(df))
  coefs_df2 <- abs(coef(df2))
#5. Odds ratio
  OR_df <- exp(coef(df))
  OR_df2 <- exp(coef(df2))
#6. OR 95% CI
  CI_df <- data.frame(exp(confint(df)))
  CI_df2 <- data.frame(exp(confint(df2)))
#7. P value function
  P_value <- function(model) {
  z <- summary(model)$coefficients/summary(model)$standard.errors
  p <- (1 - pnorm(abs(z), 0, 1)) * 2
  return(p[,2])
  }
#8 P value 
  p_val1<- P_value(df)
  p_val2<- P_value(df2)

#10 Insufficient vs Moderate data frame
  if(contrast == "Ins_Mod"){
    variable  = string 
    statistic = coefs_df[[3]]
    estimate  = OR_df[[3]]
    conf.low = CI_df[2,1]
    conf.high = CI_df[2,2]
    ci <- paste(round(OR_df[[3]],2), "(",round(CI_df[2,1],2),"-", round(CI_df[2,2],2),")")
    p.value = p_val1[[1]]
    r2_change = 0
    output <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))
  }
#11 Insufficient vs Active datafraame
  if(contrast == "Ins_Act"){
    variable  = string
    statistic = coefs_df[[4]]
    estimate  = OR_df[[4]]
    conf.low = CI_df[2,3]
    conf.high = CI_df[2,4]
    ci <- paste(round(OR_df[[4]],2), "(",round(CI_df[2,3],2),"-", round(CI_df[2,4],2),")")
    p.value = p_val1[[2]]
    r2_change = 0
    output <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))
    }
#12 Active vs Moderate table 
  if(contrast == "Act_Mod"){
    variable  = string  
    statistic = coefs_df2[[4]]
    estimate  = OR_df2[[4]]
    conf.low = CI_df2[2,3]
    conf.high = CI_df2[2,4]
    ci <- paste(round(OR_df2[[4]],2),"(",round(CI_df2[2,3],2),"-", round(CI_df2[2,4],2),")")
    p.value = p_val2[[2]]
    r2_change = 0
    output <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))
  }
  return(output)
}

Running `UniThesis()`

Now that we have a function, we need to run it. The UniThesis() function is for each health condition and compiles it into one larger data frame Ins_Mod_Uni. The code below is used to extract values for the Insufficient vs Sufficiently active comparisons.

Insufficient vs moderate

### Physical Health Conditions 
#1. Any Disability
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Any_Disability,"Any Disability", "Ins_Mod"))
#2. Arthritis
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Arthritis,"Arthritis", "Ins_Mod"))
#3. Asthma
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Athsma,"Athsma", "Ins_Mod"))
#4. Cancer
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Cancer,"Cancer", "Ins_Mod"))
#5. Chronic Pain
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Chron_pain,"Chronic Pain", "Ins_Mod"))
#6. Diabetes
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Diabetes,"Diabetes", "Ins_Mod"))
#7. Heart Attack
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Hrt_atk,"Heart Attack", "Ins_Mod"))
#8. Heart Disease
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Hrt_dis,"Heart Disease", "Ins_Mod"))
#9. High Cholesterol
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$High_chol,"High Cholesterol", "Ins_Mod"))
#10. Hypertension
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$High_bld_press,"Hypertension", "Ins_Mod"))
#11. Kidney Disease
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Kid_dis,"Kidney Disease", "Ins_Mod"))
#12. Liver Disease
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Liv_dis,"Liver Disease", "Ins_Mod"))
#13. Migrane
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Migrane,"Migrane", "Ins_Mod"))
#14. MS
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$MS,"MS", "Ins_Mod"))
#15. Osteoporosis
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Osteoporosis,"Osteoporosis", "Ins_Mod"))
#16. Rheumatoid Arthritis
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Rhum_arth,"RA", "Ins_Mod"))
#17. Stroke
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$Stroke,"Stroke", "Ins_Mod"))
### Mental Health Conditions
#18. Alcohol Use Disorder
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$LT_AUD,"AUD", "Ins_Mod"))
#19. Drug Use Disorder
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$LT_DUD,"DUD", "Ins_Mod"))
#20. Major Depressive Disorder 
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$LT_MDD,"MDD", "Ins_Mod"))
#21 Nicotine Dependence 
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$LT_NicotineDep,"Nicotine Dep", "Ins_Mod"))
#22 Posttraumatic Stress disorder
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$LT_PTSD,"PTSD", "Ins_Mod"))
#23. Suicidality
Ins_Mod_Uni<-rbind(Ins_Mod_Uni,UniThesis(thesis_NHRVS,thesis_NHRVS$LT_Suicidal,"Suicidality", "Ins_Mod"))
Ins_Mod_Uni <- Ins_Mod_Uni[-(1),]

Lets look at the table!

knitr::kable(Ins_Mod_Uni[1:3,]) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

	variable	statistic	estimate	conf.low	conf.high	ci	p.value
2	Any Disability	0.641804870673004	0.526341587737195	0.406587931217258	0.681366675474878	0.53 ( 0.41 - 0.68 )	1.10027191291273e-06
3	Arthritis	0.131641124556636	0.876655548781944	0.730159209754747	1.05254435052365	0.88 ( 0.73 - 1.05 )	0.158231184859979
4	Athsma	0.439914446936515	0.644091522728985	0.490118648374674	0.84643563558962	0.64 ( 0.49 - 0.85 )	0.00159912089713932

We will repeat this step for the Insufficient vs active comparison (Ins_Act) and Active vs Moderate Comparison (Act_Mod). As most of the analysis in this markdown are repeated I will only be showing the first one, and not the 2 other comparisons to save space.

Export as CSV

write.csv(Ins_Act_Uni, "~/Desktop/Coding/data/Ins_Act_Uni.csv")
write.csv(Ins_Mod_Uni, "~/Desktop/Coding/data/Ins_Mod_Uni.csv")
write.csv(Act_Mod_Uni, "~/Desktop/Coding/data/Act_Mod_Uni.csv")

Step Wise Multivariate Logistic Multinomial Regression

That’s a mouthful! For the following section we will be running code for the full regression model, extracting values like we did for the univariate function, and putting these values into tables for the 3 comparisons at each of the 3 blocks. The blocks of this regession are as follows:

Covariates: Sociodemograpahic and health variables
Physical and Mental Health Functioning: The short form 8 - MCS and PCS
Health Conditions: Each health condition will be ran as a seperate regression and compiled into one dataframe

Regession Code

Block 1: Covariates

This first block of the stepwise regression will do the same thing as we did for the univariate analysis, being that this was not replicated we did not make a function for it.

### Part 1 Create values for table
#1. Run the models
Cov_Reg<-multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ Age + College + Income_60k_plus + YN_work + BMI, data = thesis_NHRVS)
Cov_Reg2<-multinom(relevel(Ex_3cat,ref = "Active")  ~ Age + College + Income_60k_plus + YN_work + BMI, data = thesis_NHRVS)

#2. Coefficients 
  coefs_Cov_Reg <- abs(coef(Cov_Reg))
  coefs_Cov_Reg2 <- abs(coef(Cov_Reg2))
#3. Odds ratio
  OR_Cov_Reg <- exp(coef(Cov_Reg))
  OR_Cov_Reg2 <- exp(coef(Cov_Reg2))

#4. OR 95% CI
  CI_Cov_Reg <- data.frame(exp(confint(Cov_Reg)))
  CI_Cov_Reg2 <- data.frame(exp(confint(Cov_Reg2)))
#5. P value 
  P_value <- function(model) {
  z <- summary(model)$coefficients/summary(model)$standard.errors
  p <- (1 - pnorm(abs(z), 0, 1)) * 2
  return(p)
  }
  p_val1<- P_value(Cov_Reg)
  p_val2<- P_value(Cov_Reg2)
#6. R^2 for model
  R2_Cov_Reg <- performance::r2(Cov_Reg)
  R2_Cov_Reg2 <- performance::r2(Cov_Reg2)

  OR.CI_IvA <- paste(round(OR_Cov_Reg[2,2:6],2), "(",round(CI_Cov_Reg[2:6,3],2),"-", round(CI_Cov_Reg[2:6,4],2),")")
  OR.CI_IvM <- paste(round(OR_Cov_Reg[1,2:6],2), "(",round(CI_Cov_Reg[2:6,1],2),"-", round(CI_Cov_Reg[2:6,2],2),")")
  OR.CI_AvM <- paste(round(OR_Cov_Reg2[2,2:6],2),"(",round(CI_Cov_Reg2[2:6,3],2),"-", round(CI_Cov_Reg2[2:6,4],2),")")

#### Part 2. Make the table and export 

#1. Create data frame
cov_OR_tab <-  data.frame(cbind(
  variable = Cov_Reg$coefnames[2:length(Cov_Reg$coefnames)],
  IvA_ODDS = OR.CI_IvA,
  IvA_p = round(p_val1[2,2:5],4),
  IvM_ODDS = OR.CI_IvM,
  IvM_p = round(p_val1[1,2:5],4),
  AvM_ODDS = OR.CI_AvM, 
  AvM_p = round(p_val2[2,2:5],4)
))
#2. create the row for the R2
  paste("Block 1 R^2 =", round((R2_Cov_Reg[[1]]*100),3),"%")
  variable = paste("Block 1 R^2 =", round((R2_Cov_Reg[[1]]*100),3),"%")
  R2TAB<-data.frame(cbind(
                    variable, IvA_ODDS = NA, IvA_p = NA, IvM_ODDS = NA,
                    IvM_p = NA,AvM_ODDS = NA, AvM_p = NA))

#3. Bind the stuff together
  cov_OR_tab <- rbind(R2TAB, cov_OR_tab)
   
#4. Export
  write_csv(cov_OR_tab,"~/Desktop/Coding/data/covSF8_OR_tab.csv")

  knitr::kable(cov_OR_tab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

variable	IvA_ODDS	IvA_p	IvM_ODDS	IvM_p	AvM_ODDS	AvM_p
Block 1 R^2 = 2.642 %	NA	NA	NA	NA	NA	NA
Age	0.99 ( 0.98 - 0.99 )	3e-04	1 ( 0.99 - 1.01 )	0.4526	1.02 ( 1.01 - 1.02 )	1e-04
CollegeSome_Colege	1.8 ( 1.41 - 2.29 )	0	1.11 ( 0.85 - 1.43 )	0.4483	0.62 ( 0.46 - 0.82 )	9e-04
Income_60k_plusOVer_60k	1.28 ( 1.09 - 1.5 )	0.002	1.25 ( 1.04 - 1.51 )	0.0184	0.98 ( 0.81 - 1.18 )	0.8173
YN_workWorking	0.78 ( 0.65 - 0.94 )	0.0084	0.76 ( 0.61 - 0.95 )	0.0134	0.97 ( 0.78 - 1.2 )	0.7914
BMI	0.93 ( 0.91 - 0.94 )	3e-04	0.97 ( 0.95 - 0.98 )	0.4526	1.04 ( 1.02 - 1.06 )	1e-04

Block 2: Covariates with SF8

Block to is idential to block one with the addition of the MCS and PCS to the regression model. We are running two different regession models (i.e., one with the MCS and one with the PCS) so you can see that below.

### Part 1 SF8 values for table
#1. Run the models

  # Mental health composite score (MCS)
  Cov.MCS_Reg<-multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ Age + College + Income_60k_plus + YN_work + BMI + MCS, data = thesis_NHRVS)
  Cov.MCS_Reg2<-multinom(relevel(Ex_3cat,ref = "Active")  ~ Age + College + Income_60k_plus + YN_work + BMI + MCS, data = thesis_NHRVS)

  # Physical health composite score (PCS)
  Cov.PCS_Reg<-multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ Age + College + Income_60k_plus + YN_work + BMI + PCS, data = thesis_NHRVS)
  Cov.PCS_Reg2<-multinom(relevel(Ex_3cat,ref = "Active")  ~ Age + College + Income_60k_plus + YN_work + BMI + PCS, data = thesis_NHRVS)

#2. Coefficients 
  #MCS
    coefs_Cov.MCS_Reg <- abs(coef(Cov.MCS_Reg))
    coefs_Cov.MCS_Reg2 <- abs(coef(Cov.MCS_Reg2))
  #PCS
    coefs_Cov.PCS_Reg <- abs(coef(Cov.PCS_Reg))
    coefs_Cov.PCS_Reg2 <- abs(coef(Cov.PCS_Reg2))
#3. Odds ratio
  #MCS
    OR_Cov.MCS_Reg <- exp(coef(Cov.MCS_Reg))
    OR_Cov.MCS_Reg2 <- exp(coef(Cov.MCS_Reg2))
  #PCS
    OR_Cov.PCS_Reg <- exp(coef(Cov.PCS_Reg))
    OR_Cov.PCS_Reg2 <- exp(coef(Cov.PCS_Reg2))
#4. OR 95% CI
  #MCS
    CI_Cov.MCS_Reg <- data.frame(exp(confint(Cov.MCS_Reg)))
    CI_Cov.MCS_Reg2 <- data.frame(exp(confint(Cov.MCS_Reg2)))
  #PCS
    CI_Cov.PCS_Reg <- data.frame(exp(confint(Cov.PCS_Reg)))
    CI_Cov.PCS_Reg2 <- data.frame(exp(confint(Cov.PCS_Reg2)))
#5. P value 
  P_value <- function(model) {
  z <- summary(model)$coefficients/summary(model)$standard.errors
  p <- (1 - pnorm(abs(z), 0, 1)) * 2
  return(p)
  }
  MCS.p_val1<- P_value(Cov.MCS_Reg)
  MCS.p_val2<- P_value(Cov.MCS_Reg2)
  PCS.p_val1<- P_value(Cov.PCS_Reg)
  PCS.p_val2<- P_value(Cov.PCS_Reg2)
#6. R^2 for model
  R2_Cov.MCS_Reg <- performance::r2(Cov.MCS_Reg)
  R2_Cov.PCS_Reg <- performance::r2(Cov.PCS_Reg)

#7. Confidence intervals
  
  MCS.OR.CI_IvA <- paste(round(OR_Cov.MCS_Reg[2,7],2), "(",round(CI_Cov.MCS_Reg[7,3],2),"-", round(CI_Cov.MCS_Reg[7,4],2),")")
  MCS.OR.CI_IvM <- paste(round(OR_Cov.MCS_Reg[1,7],2), "(",round(CI_Cov.MCS_Reg[7,1],2),"-", round(CI_Cov.MCS_Reg[7,2],2),")")
  MCS.OR.CI_AvM <- paste(round(OR_Cov.MCS_Reg2[2,7],2),"(",round(CI_Cov.MCS_Reg2[7,3],2),"-", round(CI_Cov.MCS_Reg2[7,4],2),")")
  
    PCS.OR.CI_IvA <- paste(round(OR_Cov.PCS_Reg[2,7],2), "(",round(CI_Cov.PCS_Reg[7,3],2),"-", round(CI_Cov.PCS_Reg[7,4],2),")")
  PCS.OR.CI_IvM <- paste(round(OR_Cov.PCS_Reg[1,7],2), "(",round(CI_Cov.PCS_Reg[7,1],2),"-", round(CI_Cov.PCS_Reg[7,2],2),")")
  PCS.OR.CI_AvM <- paste(round(OR_Cov.PCS_Reg2[2,7],2),"(",round(CI_Cov.PCS_Reg2[7,3],2),"-", round(CI_Cov.PCS_Reg2[7,4],2),")")

#### Part 2. Make the table and export 

#1. Create data frame
SF8_OR_tab <-  data.frame(cbind(
  variable = c("MCS","PCS"),
  IvA_ODDS =c(MCS.OR.CI_IvA,PCS.OR.CI_IvA),
  IvA_p = c(round(MCS.p_val1[2,6],4),round(PCS.p_val1[2,6],4)),
  IvM_ODDS = c(MCS.OR.CI_IvM,PCS.OR.CI_IvM),
  IvM_p =  c(round(MCS.p_val1[1,6],4),round(PCS.p_val1[1,6],4)),
  AvM_ODDS = c(MCS.OR.CI_AvM,PCS.OR.CI_AvM),
  AvM_p =  c(round(MCS.p_val2[2,6],4),round(PCS.p_val2[2,6],4))
))

#2. create the row for the R2
  variable = paste("Block 2 MCS R^2 =", round((R2_Cov.MCS_Reg[[1]]*100),3),"%")
  R2TAB2<-data.frame(cbind(
                    variable, IvA_ODDS = NA, IvA_p = NA, IvM_ODDS = NA,
                    IvM_p = NA,AvM_ODDS = NA, AvM_p = NA))

  variable = paste("Block 2 PCS R^2 =", round((R2_Cov.PCS_Reg[[1]]*100),3),"%")
  R2TAB3<-data.frame(cbind(
                    variable, IvA_ODDS = NA, IvA_p = NA, IvM_ODDS = NA,
                    IvM_p = NA,AvM_ODDS = NA, AvM_p = NA))
#3. Bind the stuff together
  SF8_OR_tab <- rbind(R2TAB2,SF8_OR_tab[1,],R2TAB3,SF8_OR_tab[2,])
  Full_OR_tab <- rbind(cov_OR_tab,SF8_OR_tab)
  
#4. Export
  write_csv(SF8_OR_tab,"~/Desktop/Coding/data/SF8.cov_OR_tab.csv")
  write_csv(Full_OR_tab, "~/Desktop/Coding/data/SF8.cov_OR_tab.csv")

  knitr::kable(Full_OR_tab) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

	variable	IvA_ODDS	IvA_p	IvM_ODDS	IvM_p	AvM_ODDS	AvM_p
1	Block 1 R^2 = 2.642 %	NA	NA	NA	NA	NA	NA
2	Age	0.99 ( 0.98 - 0.99 )	3e-04	1 ( 0.99 - 1.01 )	0.4526	1.02 ( 1.01 - 1.02 )	1e-04
3	CollegeSome_Colege	1.8 ( 1.41 - 2.29 )	0	1.11 ( 0.85 - 1.43 )	0.4483	0.62 ( 0.46 - 0.82 )	9e-04
4	Income_60k_plusOVer_60k	1.28 ( 1.09 - 1.5 )	0.002	1.25 ( 1.04 - 1.51 )	0.0184	0.98 ( 0.81 - 1.18 )	0.8173
5	YN_workWorking	0.78 ( 0.65 - 0.94 )	0.0084	0.76 ( 0.61 - 0.95 )	0.0134	0.97 ( 0.78 - 1.2 )	0.7914
6	BMI	0.93 ( 0.91 - 0.94 )	3e-04	0.97 ( 0.95 - 0.98 )	0.4526	1.04 ( 1.02 - 1.06 )	1e-04
11	Block 2 MCS R^2 = 3.159 %	NA	NA	NA	NA	NA	NA
22	MCS	1.03 ( 1.02 - 1.05 )	0	1.02 ( 1.01 - 1.03 )	1e-04	0.99 ( 0.97 - 1 )	0
31	Block 2 PCS R^2 = 5.217 %	NA	NA	NA	NA	NA	NA
21	PCS	1.06 ( 1.05 - 1.07 )	0	1.04 ( 1.03 - 1.05 )	0.037	0.98 ( 0.97 - 0.99 )	1e-04

Full Model Function

The final step in my analysis was to run the full multinomial model so I can see how the health conditions are associated with physical activity level after controling for covariates. The following MultiThesis() function does the following 1. Runs the 2 different multinomial models using multinom() + Physical health condtions include the MCS + Mental health conditions include the PCS 2. Manually calculates the coefficients B, the odds ratios (OR) + 95% Confideince interval (95% CI), and p values 3. Extracts respective values for each type of comparison (e.g., sufficiently active vs moderatly active) 4. Returns respective values as a row in a dataframe which can be compiled into a larger data fram per comparison

# Multivariate Logistic Multinomial Regression Function

MultiThesis<-function(Data,condition,string,type,contrast){

#1. run model 1
  if(type == "PHC"){
  df <- multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ condition + Age + College + Income_60k_plus + YN_work + BMI + MCS, data = Data)
  }
  if(type == "MHC"){ df <- multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ condition + Age + College + Income_60k_plus + YN_work + BMI + PCS, data = Data)
  }
#2 Specify the reference group2
  Data$Ex_3cat_act <- relevel(Data$Ex_3cat, ref = "Moderate")
#3. run model 2
  if(type == "PHC"){
   df2 <- multinom(relevel(Ex_3cat,ref = "Moderate") ~ condition + Age + College + Income_60k_plus + YN_work + BMI + MCS, data = Data)
  }
  if(type == "MHC"){ 
    df2 <- multinom(relevel(Ex_3cat,ref = "Moderate") ~ condition + Age + College + Income_60k_plus + YN_work + BMI + PCS, data = Data)
  }
#4 Coefficients 
  coefs_df <- abs(coef(df))
  coefs_df2 <- abs(coef(df2))
#5. Odds ratio
  OR_df <- exp(coef(df))
  OR_df2 <- exp(coef(df2))
#6. OR 95% CI
  CI_df <- data.frame(exp(confint(df)))
  CI_df2 <- data.frame(exp(confint(df2)))
#7. P value function
  P_value <- function(model) {
  z <- summary(model)$coefficients/summary(model)$standard.errors
  p <- (1 - pnorm(abs(z), 0, 1)) * 2
  return(p[,2])
  }
#8 P value 
  p_val1<- P_value(df)
  p_val2<- P_value(df2)

#10 Insufficient vs Moderate data frame
  if(contrast == "Ins_Mod"){
    variable  = string 
    statistic = coefs_df[[3]]
    estimate  = OR_df[[3]]
    conf.low = CI_df[2,1]
    conf.high = CI_df[2,2]
    ci <- paste(round(OR_df[[3]],2), "(",round(CI_df[2,1],2),"-", round(CI_df[2,2],2),")")
    p.value = p_val1[[1]]
    r2_change = 0
    output <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))
  }
#11 Insufficient vs Active datafraame
  if(contrast == "Ins_Act"){
    variable  = string
    statistic = coefs_df[[4]]
    estimate  = OR_df[[4]]
    conf.low = CI_df[2,3]
    conf.high = CI_df[2,4]
    ci <- paste(round(OR_df[[4]],2), "(",round(CI_df[2,3],2),"-", round(CI_df[2,4],2),")")
    p.value = p_val1[[2]]
    r2_change = 0
    output <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))
    }
#12 Active vs Moderate table 
  if(contrast == "Act_Mod"){
    variable  = string  
    statistic = coefs_df2[[4]]
    estimate  = OR_df2[[4]]
    conf.low = CI_df2[2,3]
    conf.high = CI_df2[2,4]
    ci <- paste(round(OR_df2[[4]],2),"(",round(CI_df2[2,3],2),"-", round(CI_df2[2,4],2),")")
    p.value = p_val2[[2]]
    r2_change = 0
    output <- data.frame(cbind(variable,statistic,estimate,conf.low,conf.high,ci,p.value,r2_change))
  }
  return(output)
}

#Its done!

R^2 Change

Before we can create the dataframes with the multinomial regession output, I need to find the R^2 change for all of the variables. The following function allows us to look at the R^2 change between the covariates and the IV within the multinomial model. As this is a multinomial model, the R^2 change is only observed for the model as a whole. I created the r2_chan() function to do the following: 1. Run a multinomial model without the predictor (i.e., just covariates) and with the predictor (i.e., the full model) 2. Calculate the R^2 for both models, and find the difference between the covariates and full model. 3. turn this R^2 change into a percentage, and save it as a string with the % sign.

#assign the variable to 0
r2_change = 0

#create the function 
  r2_chan<-function(condition){
  multi_cov<-multinom(relevel(Ex_3cat,ref = "Insufficient") ~ Age + College + Income_60k_plus + YN_work + BMI, data = thesis_NHRVS)
  df <-multinom(relevel(Ex_3cat,ref = "Insufficient") ~ condition + Age + College + Income_60k_plus + YN_work + BMI, data = thesis_NHRVS)
Delta <-  (r2(df)[[1]] -  r2(multi_cov)[[1]])*100
if(Delta >= 0.01){
output <- paste(round(Delta,2), "%")
}
if(Delta < 0.01){
output <- "< 0.01 %"
}
return(output)
  }

Now we get to run the r2_chan() function for all of the variables. The code below compiles a data frame of the R2CHANGE for each variable within each block of the step wise regression.

# Block 1: Covariates
    R2CHANGE <- NA
    R2CHANGE  <- rbind(R2CHANGE,NA,NA,NA,NA)
# Block 2:  MCS / PCS 
  #MCS
    MCS_R2 = 
    R2CHANGE  <- rbind(R2CHANGE,NA,
                       paste(round(((r2(Cov.MCS_Reg)[[1]] - r2(Cov_Reg)[[1]])*100),2), "%"))
  #PCS
    R2CHANGE  <-rbind(R2CHANGE,NA,
                      paste(round(((r2(Cov.PCS_Reg)[[1]] - r2(Cov_Reg)[[1]])*100),2), "%"),
                      NA)
    
# Block 3: run the function for all predictors
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Any_Disability)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Arthritis)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Athsma)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Cancer)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Chron_pain )[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Diabetes )[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Hrt_atk )[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Hrt_dis)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$High_chol)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$High_bld_press)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Kid_dis)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Liv_dis)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Migrane)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$MS)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Osteoporosis)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Rhum_arth)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$Stroke)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$LT_AUD)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$LT_DUD)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$LT_MDD)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$LT_NicotineDep)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$LT_PTSD)[[1]])
    R2CHANGE <- rbind(R2CHANGE,r2_chan(thesis_NHRVS$LT_Suicidal)[[1]])

Lastly, we convert R2CHANGE to a dataframe and export it as a CSV. You can see the dataframe below.

R2CHANGE<-data.frame(R2CHANGE)
  knitr::kable(R2CHANGE) %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

	R2CHANGE
R2CHANGE	NA
X	NA
X.1	NA
X.2	NA
X.3	NA
X.4	NA
X.5	0.52 %
X.6	NA
X.7	2.58 %
X.8	NA
X.9	0.65 %
X.10	0.03 %
X.11	0.16 %
X.12	0.04 %
X.13	0.25 %
X.14	0.14 %
X.15	0.09 %
X.16	0.03 %
X.17	0.11 %
X.18	0.09 %
X.19	0.06 %
X.20	< 0.01 %
X.21	< 0.01 %
X.22	< 0.01 %
X.23	0.08 %
X.24	< 0.01 %
X.25	0.12 %
X.26	0.05 %
X.27	0.11 %
X.28	0.12 %
X.29	0.21 %
X.30	0.02 %
X.31	0.07 %

write_csv(R2CHANGE,"~/Desktop/Coding/data/R2change.csv")

We have the R^2 change values for each value, but for the table. I want to know what the the range of R2 for the full model within block 3. The highest and lowest values are for multiple socrosis (lowest R^2) and Lifetime nictotine dependence (highest R^2)

# Lowest R2 = MS
r2(multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ MS + Age + College + Income_60k_plus + YN_work + BMI + MCS, data = thesis_NHRVS))

## # weights:  27 (16 variable)
## initial  value 3897.876400 
## iter  10 value 3676.884042
## iter  20 value 3633.802638
## final  value 3633.801616 
## converged

## # R2 for Generalized Linear Regression
##   McFadden's R2: 0.032

# Highest R2 =  LT_NicotineDep

r2(multinom(relevel(Ex_3cat,ref = "Insufficient")  ~ LT_NicotineDep + Age + College + Income_60k_plus + YN_work + BMI + PCS, data = thesis_NHRVS))

## # weights:  27 (16 variable)
## initial  value 3897.876400 
## iter  10 value 3583.526084
## iter  20 value 3552.954623
## final  value 3552.953989 
## converged

## # R2 for Generalized Linear Regression
##   McFadden's R2: 0.053

Block 3: The Multivaraite Multinomial

Run the `MultiThesis()` function

Now that we have created our function, we can use it to make our data frames for each contrast. Being that the function just creates the values needed for each condition, we still need to add it to the exsisting data frame (i.e.,rbind())

Insufficient vs Moderate Dataframe

### Physical Health Conditions 
Ins_Mod = 0
#1. Any Disability
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Any_Disability,"Any Disability", "PHC", "Ins_Mod"))
#2. Arthritis
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Arthritis,"Arthritis", "PHC", "Ins_Mod"))
#3. Asthma
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Athsma,"Athsma", "PHC", "Ins_Mod"))
#4. Cancer
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Cancer,"Cancer", "PHC", "Ins_Mod"))
#5. Chronic Pain
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Chron_pain,"Chronic Pain", "PHC", "Ins_Mod"))
#6. Diabetes
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Diabetes,"Diabetes", "PHC", "Ins_Mod"))
#7. Heart Attack
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Hrt_atk,"Heart Attack", "PHC", "Ins_Mod"))
#8. Heart Disease
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Hrt_dis,"Heart Disease", "PHC", "Ins_Mod"))
#9. High Cholesterol
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$High_chol,"High Cholesterol", "PHC", "Ins_Mod"))
#10. Hypertension
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$High_bld_press,"Hypertension", "PHC", "Ins_Mod"))
#11. Kidney Disease
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Kid_dis,"Kidney Disease","PHC", "Ins_Mod"))
#12. Liver Disease
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Liv_dis,"Liver Disease", "PHC", "Ins_Mod"))
#13. Migrane
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Migrane,"Migrane", "PHC", "Ins_Mod"))
#14. MS
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$MS,"MS", "PHC", "Ins_Mod"))
#15. Osteoporosis
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Osteoporosis,"Osteoporosis", "PHC", "Ins_Mod"))
#16. Rheumatoid Arthritis
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Rhum_arth,"RA", "PHC", "Ins_Mod"))
#17. Stroke
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$Stroke,"Stroke", "PHC", "Ins_Mod"))
### Mental Health Conditions
#18. Alcohol Use Disorder
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$LT_AUD,"AUD", "MHC", "Ins_Mod"))
#19. Drug Use Disorder
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$LT_DUD,"DUD", "MHC", "Ins_Mod"))
#20. Major Depressive Disorder 
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$LT_MDD,"MDD", "MHC", "Ins_Mod"))
#21 Nicotine Dependence 
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$LT_NicotineDep,"Nicotine Dep", "MHC", "Ins_Mod"))
#22 Posttraumatic Stress disorder
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$LT_PTSD,"PTSD", "MHC", "Ins_Mod"))
#23. Suicidality
Ins_Mod<-rbind(Ins_Mod,MultiThesis(thesis_NHRVS,thesis_NHRVS$LT_Suicidal,"Suicidality", "MHC", "Ins_Mod"))
Ins_Mod

We’ve now compiled all of our helath conditions and added them into the dataframe. However, we still have that pesky row of zeros. We can use negative indexing to remove this row from the dataframe leaving us with just the stuff we care about.

Ins_Mod <- Ins_Mod[-(1),]

We can now add the R^2 change to the dataframe.

R2_Conditions<- R2CHANGE[11:33,]
Ins_Mod$r2_change <- R2_Conditions
knitr::kable(Ins_Mod[1:3,]) %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

	variable	statistic	estimate	conf.low	conf.high	ci	p.value	r2_change
2	Any Disability	0.471580584343876	0.624015179776385	0.475788279165653	0.818420632963471	0.62 ( 0.48 - 0.82 )	0.000654183824832177	0.65 %
3	Arthritis	0.0253984833697793	0.974921344671965	0.807302237988491	1.17734298701466	0.97 ( 0.81 - 1.18 )	0.791885208764693	0.03 %
4	Athsma	0.340659436074763	0.711301110459597	0.539141100312293	0.938435725727437	0.71 ( 0.54 - 0.94 )	0.015980219803454	0.16 %

Now repeate that step for the other 2 comparisons!

Export Tables to CSV

write_csv(Ins_Mod, "~/Desktop/Coding/data/Ins_Mod.csv")
write_csv(Ins_Act, "~/Desktop/Coding/data/Ins_Act.csv")
write_csv(Act_Mod, "~/Desktop/Coding/data/Act_Mod.csv")

Create Tables for Multinomial Logistic Regression Output

We have everything we need for our final table output!

#1. Import tables
   r2_change <- read_csv("~/Desktop/Coding/data/R2change.csv")
   cov.SF8_OR_tab <- read.csv("~/Desktop/Coding/data/SF8.cov_OR_tab.csv")
   Ins_Mod <- read_csv("~/Desktop/Coding/data/Ins_Mod.csv")
   Ins_Act <- read_csv("~/Desktop/Coding/data/Ins_Act.csv")
   Act_Mod <- read_csv( "~/Desktop/Coding/data/Act_Mod.csv")

#2. Create dataframes to be merged.
    OR_Table = data.frame(
      Ins_Act$variable, Ins_Act$ci, Ins_Mod$ci, Act_Mod$ci,Act_Mod$r2_change
      )
    colnames(OR_Table) <- c("Variable","IvA_ODDS",
                          "IvM_ODDS", "AvM_ODDS",
                          "R2_Change")
    
#3. Create a block 3 row
  Block_3 = cbind(
    Variable = "Block 3: Health Conditions R^2 range = 0.032 - 0.052",
    IvA_ODDS = NA,
    IvM_ODDS = NA,
    AvM_ODDS = NA,
    R2_Change = NA
  )

#4. pull out values from covariates table
  Cov_table = data.frame(
    Full_OR_tab$variable,
    Full_OR_tab$IvA_ODDS, 
    Full_OR_tab$IvM_ODDS,
    Full_OR_tab$AvM_ODDS
  )
  Cov_table$R2_Change = NA
  colnames(Cov_table) <- c("Variable","IvA_ODDS",
                          "IvM_ODDS", "AvM_ODDS",
                          "R2_Change")

#5. Merge all the data frames
  OR_Table <- rbind(Cov_table,Block_3,OR_Table)

#6. update R2 change df and add it to OR tab
  r2_change = rbind(r2_change[1:6,],NA,r2_change[7:33,])
  OR_Table <- cbind(OR_Table[,1:4],r2_change)
  
#6. Rename columns for table 
  colnames(OR_Table) <- c("Variable","Insufficient vs Sufficient ORs",
                          "Insufficient vs Moderate ORs", "Active vs Moderate ORs", "R^2 Change")
  
  
#7. View Table
    knitr::kable(OR_Table) %>%
          kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Variable	Insufficient vs Sufficient ORs	Insufficient vs Moderate ORs	Active vs Moderate ORs	R^2 Change
Block 1 R^2 = 2.642 %	NA	NA	NA	NA
Age	0.99 ( 0.98 - 0.99 )	1 ( 0.99 - 1.01 )	1.02 ( 1.01 - 1.02 )	NA
CollegeSome_Colege	1.8 ( 1.41 - 2.29 )	1.11 ( 0.85 - 1.43 )	0.62 ( 0.46 - 0.82 )	NA
Income_60k_plusOVer_60k	1.28 ( 1.09 - 1.5 )	1.25 ( 1.04 - 1.51 )	0.98 ( 0.81 - 1.18 )	NA
YN_workWorking	0.78 ( 0.65 - 0.94 )	0.76 ( 0.61 - 0.95 )	0.97 ( 0.78 - 1.2 )	NA
BMI	0.93 ( 0.91 - 0.94 )	0.97 ( 0.95 - 0.98 )	1.04 ( 1.02 - 1.06 )	NA
Block 2 MCS R^2 = 3.159 %	NA	NA	NA	NA
MCS	1.03 ( 1.02 - 1.05 )	1.02 ( 1.01 - 1.03 )	0.99 ( 0.97 - 1 )	0.52 %
Block 2 PCS R^2 = 5.217 %	NA	NA	NA	NA
PCS	1.06 ( 1.05 - 1.07 )	1.04 ( 1.03 - 1.05 )	0.98 ( 0.97 - 0.99 )	2.58 %
Block 3: Health Conditions R^2 range = 0.032 - 0.052	NA	NA	NA	NA
Any Disability	0.53 ( 0.42 - 0.67 )	0.62 ( 0.48 - 0.82 )	0.85 ( 0.63 - 1.15 )	0.65 %
Arthritis	1.14 ( 0.97 - 1.34 )	0.97 ( 0.81 - 1.18 )	1.17 ( 0.97 - 1.41 )	0.03 %
Athsma	0.75 ( 0.59 - 0.94 )	0.71 ( 0.54 - 0.94 )	1.05 ( 0.78 - 1.4 )	0.16 %
Cancer	1.18 ( 0.97 - 1.44 )	0.98 ( 0.78 - 1.24 )	1.21 ( 0.96 - 1.52 )	0.04 %
Chronic Pain	0.79 ( 0.66 - 0.94 )	0.73 ( 0.59 - 0.9 )	1.08 ( 0.87 - 1.35 )	0.25 %
Diabetes	0.75 ( 0.61 - 0.91 )	0.93 ( 0.74 - 1.16 )	0.81 ( 0.64 - 1.02 )	0.14 %
Heart Attack	0.7 ( 0.52 - 0.94 )	0.84 ( 0.61 - 1.17 )	0.83 ( 0.58 - 1.19 )	0.09 %
Heart Disease	0.89 ( 0.71 - 1.11 )	1.01 ( 0.79 - 1.29 )	0.88 ( 0.68 - 1.15 )	0.03 %
High Cholesterol	0.82 ( 0.7 - 0.95 )	0.95 ( 0.79 - 1.15 )	0.86 ( 0.71 - 1.03 )	0.11 %
Hypertension	0.82 ( 0.7 - 0.97 )	0.83 ( 0.68 - 1 )	1 ( 0.82 - 1.21 )	0.09 %
Kidney Disease	0.99 ( 0.72 - 1.35 )	0.67 ( 0.45 - 1.01 )	1.47 ( 0.96 - 2.24 )	0.06 %
Liver Disease	1.14 ( 0.66 - 1.98 )	1.26 ( 0.68 - 2.35 )	0.91 ( 0.48 - 1.71 )	< 0.01 %
Migrane	1.16 ( 0.88 - 1.55 )	1.1 ( 0.78 - 1.54 )	1.06 ( 0.75 - 1.5 )	< 0.01 %
MS	1.42 ( 0.43 - 4.65 )	1.25 ( 0.29 - 5.3 )	1.13 ( 0.29 - 4.45 )	< 0.01 %
Osteoporosis	0.69 ( 0.49 - 0.96 )	0.89 ( 0.62 - 1.3 )	0.77 ( 0.52 - 1.14 )	0.08 %
RA	1.01 ( 0.73 - 1.41 )	1.15 ( 0.79 - 1.66 )	0.88 ( 0.6 - 1.3 )	< 0.01 %
Stroke	0.58 ( 0.38 - 0.88 )	0.83 ( 0.53 - 1.3 )	0.7 ( 0.42 - 1.15 )	0.12 %
AUD	0.97 ( 0.82 - 1.13 )	1.1 ( 0.91 - 1.33 )	0.88 ( 0.73 - 1.06 )	0.05 %
DUD	0.84 ( 0.65 - 1.07 )	0.83 ( 0.62 - 1.11 )	1.01 ( 0.75 - 1.36 )	0.11 %
MDD	1.03 ( 0.84 - 1.25 )	1.06 ( 0.84 - 1.35 )	0.96 ( 0.76 - 1.22 )	0.12 %
Nicotine Dep	0.89 ( 0.72 - 1.09 )	0.71 ( 0.55 - 0.91 )	1.25 ( 0.96 - 1.62 )	0.21 %
PTSD	1.38 ( 1.06 - 1.79 )	1.09 ( 0.79 - 1.49 )	1.27 ( 0.93 - 1.74 )	0.02 %
Suicidality	0.94 ( 0.72 - 1.24 )	1.04 ( 0.76 - 1.43 )	0.91 ( 0.66 - 1.25 )	0.07 %

#8 Export Table
    write_csv(OR_Table,"~/Desktop/Coding/data/Thesis_OR_Table.csv")

Logistic Multinomial Regresion