3.1 UNWEIGHTED regression
3.1a Regresion Function
regre_table_uw <- function(fomular,data,imputed_data,outcome,table_labels){
model <- glmmTMB(as.formula(fomular),data,
family = binomial(link = "logit"))
print(summary(model))
Coef1 <- data.frame(coefficients(summary(model))$cond)%>%
mutate(variables = rownames(.)) %>%
filter(variables != "(Intercept)") %>%
mutate(Coef = format(round(exp(Estimate), digits = 2), nsmall = 2)) %>%
mutate('95% CI' = paste0(format(round( exp(Estimate - 1.96*Std..Error) , digits = 2),nsmall = 2),", " ,
format(round( exp(Estimate + 1.96*Std..Error) , digits = 2),nsmall = 2))) %>%
mutate('p-value' = format(round(2*(1 - pnorm(abs(z.value))), digits = 3 ), nsmall = 3)) %>%
mutate('p-value' = ifelse(`p-value` == "0.000","<0.001",`p-value`))
imputed_model <- lapply(imputed_data, function(x) glmmTMB(as.formula(fomular),x,
family = binomial(link = "logit")))
pooled_fit <- pool(imputed_model)
Coef2 <- data.frame(summary(pooled_fit)) %>%
filter(term != "(Intercept)") %>%
mutate('Imputed Coef ' = format(round(exp(estimate), digits = 2), nsmall = 2)) %>%
mutate('Imputed 95% CI ' = paste0(format(round( exp(estimate - 1.96*std.error) , digits = 2),nsmall = 2)," , " ,
format(round( exp(estimate + 1.96*std.error) , digits = 2),nsmall = 2))) %>%
mutate('p-value ' = format(round(p.value, digits = 3 ), nsmall = 3)) %>%
mutate('Imputed p-value ' = ifelse(`p-value ` == "0.000","<0.001",`p-value `))
finaltable <- cbind(dplyr::select(Coef1,variables,Coef,'95% CI','p-value'),
dplyr::select(Coef2,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `))
finaltable <- table_labels %>%
dplyr::left_join(finaltable, by = "variables") %>%
dplyr::select(-variables) %>%
mutate(Outcome=outcome)
list(finaltable)
}
3.1b Run unweighted regression models
results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Variable = c("PM2.5 per 10, centred(10)","Round 2","Round 3","Round 4","Age per 10yrs",
"Female", "Post-secondary education","In paid employment",
"IRSAD Score per 100","Work exposed","Pre-mine fire Asthma or COPD",
"Ex-smoker", "Current smoker","ln(cigarette packyears + 1)",
"Pm2.5:Round 2","Pm2.5:Round 3",
"Pm2.5:Round 4"),
variables = c("firePM_cent10per10","TimeR2","TimeR3","TimeR4","age_per10",
"genderFemale", "educationPost-secondary","employmentstatusPaid employment",
"Base_IRSAD_Score_per100","Base_Work_Exposed_2grpExposed","Base_Asthma_COPDyes",
"Base_smokestatusEx-smoker",
"Base_smokestatusCurrent smoker","Base_lnpackyears",
"firePM_cent10per10:TimeR2","firePM_cent10per10:TimeR3",
"firePM_cent10per10:TimeR4"))
## Results
model <- regre_table_uw(fomular,dta_long,imputed_dta,outcome,table_labels)
results[[outcome]] <- model
}
results_combined <- map_dfr(results, bind_rows)
results<-results_combined %>%
dplyr::select(Outcome,Variable,Coef,'95% CI','p-value',`Imputed Coef `,
`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(results,here::here("Results","All519_unweighted.xlsx"))
# saveRDS( results, here::here("Results","All519_unweighted.rds"))
results<-readRDS(here::here("Results","All519_unweighted.rds"))
# for (i in 1:length(outcomes)) {
#
# print(knitr::kable(results[results$Outcome==outcomes[i], ],format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))
# }
3.1c Linear function
lincom_table_uw <- function(fomular,imputed_data,outcome,labels){
mod<- lapply(imputed_data, function(x) multcomp::glht(glmmTMB(as.formula(fomular),x,
family = binomial(link = "logit")),
linfct =
c("firePM_cent10per10 = 0",
"firePM_cent10per10 + firePM_cent10per10:TimeR2 = 0",
"firePM_cent10per10 + firePM_cent10per10:TimeR3 = 0",
"firePM_cent10per10 + firePM_cent10per10:TimeR4 = 0"
)
)
)
pooled_lfct <- pool(mod)
Coef <- data.frame(summary(pooled_lfct)) %>%
mutate('Imputed Coef ' = format(round(exp(estimate), digits = 2), nsmall = 2)) %>%
mutate('Imputed 95% CI ' = paste0(format(round( exp(estimate - 1.96*std.error) , digits = 2),nsmall = 2)," , " ,
format(round( exp(estimate + 1.96*std.error) , digits = 2),nsmall = 2))) %>%
mutate('p-value ' = format(round(p.value, digits = 3 ), nsmall = 3)) %>%
mutate('Imputed p-value ' = ifelse(`p-value ` == "0.000","<0.001",`p-value `))
finaltable <- Coef %>% dplyr::select(contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
finaltable <- table_labels %>%
dplyr::left_join(finaltable, by = "contrast") %>%
dplyr::select(-contrast) %>%
mutate(Outcome=outcome)
list(finaltable)
}
3.1d Run linear combinations
lin_results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Contrast = c("PM2.5 per 10, R1","PM2.5 R1 + Pm2.5:Round 2",
"PM2.5 R1 + Pm2.5:Round 3","PM2.5 R1 + Pm2.5:Round 4"),
contrast = c("firePM_cent10per10","firePM_cent10per10 + firePM_cent10per10:TimeR2",
"firePM_cent10per10 + firePM_cent10per10:TimeR3",
"firePM_cent10per10 + firePM_cent10per10:TimeR4"))
## Results
model <- lincom_table_uw(fomular,imputed_dta,outcome,table_labels)
lin_results[[outcome]] <- model
}
lin_combined <- map_dfr(lin_results, bind_rows)
lin_res<-lin_combined %>%
dplyr::select(Outcome,Contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(lin_res,here::here("Results","lin_all519_unweighted.xlsx"))
# saveRDS( lin_res, here::here("Results","lin_all519_unweighted.rds"))
lin_results<-readRDS(here::here("Results","lin_all519_unweighted.rds"))
# for (i in 1:length(outcomes)) {
#
# print(knitr::kable(results[results$Outcome==outcomes[i], ],format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))
# }
3.1e Sensitivity analyses: Observed imputed UNweighted linear
combinations
lin_results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Contrast = c("PM2.5 per 10, R1","PM2.5 R1 + Pm2.5:Round 2",
"PM2.5 R1 + Pm2.5:Round 3","PM2.5 R1 + Pm2.5:Round 4"),
contrast = c("firePM_cent10per10","firePM_cent10per10 + firePM_cent10per10:TimeR2",
"firePM_cent10per10 + firePM_cent10per10:TimeR3",
"firePM_cent10per10 + firePM_cent10per10:TimeR4"))
## Results
model <- lincom_table_uw(fomular,imputed_observed,outcome,table_labels)
lin_results[[outcome]] <- model
}
lin_combined <- map_dfr(lin_results, bind_rows)
lin_res<-lin_combined %>%
dplyr::select(Outcome,Contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(lin_res,here::here("Results","lin_observed_unweighted.xlsx"))
# saveRDS( lin_res, here::here("Results","lin_observed_unweighted.rds"))
lin_results<-readRDS(here::here("Results","lin_observed_unweighted.rds"))
# for (i in 1:length(outcomes)) {
#
# print(knitr::kable(results[results$Outcome==outcomes[i], ],format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))
# }
3.1f UNWeighted Stratified
3.1f.1 By Pre-fire Asthma
imputed_asthma <- list()
for (i in 1:50) {
imputed_asthma[[i]] <-imputed_dta[[i]] %>%
filter(Base_Asthma=="Yes")
}
lin_results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Contrast = c("PM2.5 per 10, R1","PM2.5 R1 + Pm2.5:Round 2",
"PM2.5 R1 + Pm2.5:Round 3","PM2.5 R1 + Pm2.5:Round 4"),
contrast = c("firePM_cent10per10","firePM_cent10per10 + firePM_cent10per10:TimeR2",
"firePM_cent10per10 + firePM_cent10per10:TimeR3",
"firePM_cent10per10 + firePM_cent10per10:TimeR4"))
## Results
model <- lincom_table_uw(fomular,imputed_asthma,outcome,table_labels)
lin_results[[outcome]] <- model
}
lin_combined <- map_dfr(lin_results, bind_rows)
lin_res<-lin_combined %>%
dplyr::select(Outcome,Contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(lin_res,here::here("Results","lin_all519_unw_strata_Asthma.xlsx"))
# saveRDS( lin_res, here::here("Results","lin_all519_unw_strata_Asthma.rds"))
3.1f.2 No Pre-fire Asthma
imputed_NONasthma <- list()
for (i in 1:50) {
imputed_NONasthma[[i]] <-imputed_dta[[i]] %>%
filter(Base_Asthma=="No")
}
lin_results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Contrast = c("PM2.5 per 10, R1","PM2.5 R1 + Pm2.5:Round 2",
"PM2.5 R1 + Pm2.5:Round 3","PM2.5 R1 + Pm2.5:Round 4"),
contrast = c("firePM_cent10per10","firePM_cent10per10 + firePM_cent10per10:TimeR2",
"firePM_cent10per10 + firePM_cent10per10:TimeR3",
"firePM_cent10per10 + firePM_cent10per10:TimeR4"))
## Results
model <- lincom_table_uw(fomular,imputed_NONasthma,outcome,table_labels)
lin_results[[outcome]] <- model
}
lin_combined <- map_dfr(lin_results, bind_rows)
lin_res<-lin_combined %>%
dplyr::select(Outcome,Contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(lin_res,here::here("Results","lin_all519_unw_strata_NONAsthma.xlsx"))
# saveRDS( lin_res, here::here("Results","lin_all519_unw_strata_NONAsthma.rds"))
4 WEIGHTED regression
4.1 Regresion Function WEIGHTED
regre_table <- function(fomular,data,imputed_data,outcome,table_labels){
model <- glmmTMB(as.formula(fomular),data,weights=Base_sample_weight,
family = binomial(link = "logit"))
print(summary(model))
Coef1 <- data.frame(coefficients(summary(model))$cond)%>%
mutate(variables = rownames(.)) %>%
filter(variables != "(Intercept)") %>%
mutate(Coef = format(round(exp(Estimate), digits = 2), nsmall = 2)) %>%
mutate('95% CI' = paste0(format(round( exp(Estimate - 1.96*Std..Error) , digits = 2),nsmall = 2),", " ,
format(round( exp(Estimate + 1.96*Std..Error) , digits = 2),nsmall = 2))) %>%
mutate('p-value' = format(round(2*(1 - pnorm(abs(z.value))), digits = 3 ), nsmall = 3)) %>%
mutate('p-value' = ifelse(`p-value` == "0.000","<0.001",`p-value`))
imputed_model <- lapply(imputed_data, function(x) glmmTMB(as.formula(fomular),x,weights=Base_sample_weight,
family = binomial(link = "logit")))
pooled_fit <- pool(imputed_model)
Coef2 <- data.frame(summary(pooled_fit)) %>%
filter(term != "(Intercept)") %>%
mutate('Imputed Coef ' = format(round(exp(estimate), digits = 2), nsmall = 2)) %>%
mutate('Imputed 95% CI ' = paste0(format(round( exp(estimate - 1.96*std.error) , digits = 2),nsmall = 2)," , " ,
format(round( exp(estimate + 1.96*std.error) , digits = 2),nsmall = 2))) %>%
mutate('p-value ' = format(round(p.value, digits = 3 ), nsmall = 3)) %>%
mutate('Imputed p-value ' = ifelse(`p-value ` == "0.000","<0.001",`p-value `))
finaltable <- cbind(dplyr::select(Coef1,variables,Coef,'95% CI','p-value'),
dplyr::select(Coef2,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `))
finaltable <- table_labels %>%
dplyr::left_join(finaltable, by = "variables") %>%
dplyr::select(-variables) %>%
mutate(Outcome=outcome)
list(finaltable)
}
4.2 Weighted Linear function
lincom_table <- function(fomular,imputed_data,outcome,labels){
mod<- lapply(imputed_data, function(x) multcomp::glht(glmmTMB(as.formula(fomular),x, weights=Base_sample_weight,
family = binomial(link = "logit")),
linfct =
c("firePM_cent10per10 = 0",
"firePM_cent10per10 + firePM_cent10per10:TimeR2 = 0",
"firePM_cent10per10 + firePM_cent10per10:TimeR3 = 0",
"firePM_cent10per10 + firePM_cent10per10:TimeR4 = 0"
)
)
)
pooled_lfct <- pool(mod)
Coef <- data.frame(summary(pooled_lfct)) %>%
mutate('Imputed Coef ' = format(round(exp(estimate), digits = 2), nsmall = 2)) %>%
mutate('Imputed 95% CI ' = paste0(format(round( exp(estimate - 1.96*std.error) , digits = 2),nsmall = 2)," , " ,
format(round( exp(estimate + 1.96*std.error) , digits = 2),nsmall = 2))) %>%
mutate('p-value ' = format(round(p.value, digits = 3 ), nsmall = 3)) %>%
mutate('Imputed p-value ' = ifelse(`p-value ` == "0.000","<0.001",`p-value `))
finaltable <- Coef %>% dplyr::select(contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
finaltable <- table_labels %>%
dplyr::left_join(finaltable, by = "contrast") %>%
dplyr::select(-contrast) %>%
mutate(Outcome=outcome)
list(finaltable)
}
4.5c Run weighted regression
results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Variable = c("PM2.5 per 10, centred(10)","Round 2","Round 3","Round 4","Age per 10yrs",
"Female", "Post-secondary education","In paid employment",
"IRSAD Score per 100","Work exposed","Pre-mine fire Asthma or COPD",
"Ex-smoker", "Current smoker","ln(cigarette packyears + 1)",
"Pm2.5:Round 2","Pm2.5:Round 3",
"Pm2.5:Round 4"),
variables = c("firePM_cent10per10","TimeR2","TimeR3","TimeR4","age_per10",
"genderFemale", "educationPost-secondary","employmentstatusPaid employment",
"Base_IRSAD_Score_per100","Base_Work_Exposed_2grpExposed","Base_Asthma_COPDyes",
"Base_smokestatusEx-smoker",
"Base_smokestatusCurrent smoker","Base_lnpackyears",
"firePM_cent10per10:TimeR2","firePM_cent10per10:TimeR3",
"firePM_cent10per10:TimeR4"))
## Results
model <- regre_table(fomular,dta_long,imputed_dta,outcome,table_labels)
results[[outcome]] <- model
}
results_combined <- map_dfr(results, bind_rows)
results<-results_combined %>%
dplyr::select(Outcome,Variable,Coef,'95% CI','p-value',`Imputed Coef `,
`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(results,here::here("Results","All519_weighted.xlsx"))
# saveRDS( results, here::here("Results","All519_weighted.rds"))
#
results<-readRDS(here::here("Results","All519_weighted.rds"))
# for (i in 1:length(outcomes)) {
#
# print(knitr::kable(results[results$Outcome==outcomes[i], ],format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))
# }
#
4.5d Run weighted linear combinations
lin_results <- list()
#loop though all outcomes,
for (i in 1:length(outcomes)) {
#outcome variable
outcome <- outcomes[i]
#regression table
fomular <- paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")
table_labels <- tibble(Contrast = c("PM2.5 per 10, R1","PM2.5 R1 + Pm2.5:Round 2",
"PM2.5 R1 + Pm2.5:Round 3","PM2.5 R1 + Pm2.5:Round 4"),
contrast = c("firePM_cent10per10","firePM_cent10per10 + firePM_cent10per10:TimeR2",
"firePM_cent10per10 + firePM_cent10per10:TimeR3",
"firePM_cent10per10 + firePM_cent10per10:TimeR4"))
## Results
model <- lincom_table(fomular,imputed_dta,outcome,table_labels)
lin_results[[outcome]] <- model
}
lin_combined <- map_dfr(lin_results, bind_rows)
lin_res<-lin_combined %>%
dplyr::select(Outcome,Contrast,`Imputed Coef `,`Imputed 95% CI `,`Imputed p-value `)
# write_xlsx(lin_res,here::here("Results","lin_all519_weighted.xlsx"))
# saveRDS( lin_res, here::here("Results","lin_all519_weighted.rds"))
lin_results<-readRDS(here::here("Results","lin_all519_weighted.rds"))
# for (i in 1:length(outcomes)) {
#
# print(knitr::kable(results[results$Outcome==outcomes[i], ],format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))
# }
5. Interaction p-values
5.1 Unweighted Interaction tests
p_load("mitml")
int_test <- list()
for (i in 1:length(outcomes)) {
outcome <- outcomes[i]
mod1 <- lapply(imputed_dta, function(x) glmmTMB(as.formula(paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")),x,family = binomial(link = "logit")))
mod2 <- lapply(imputed_dta, function(x) glmmTMB(as.formula(paste0( outcome, " ~ firePM_cent10per10 + Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")),x,family = binomial(link = "logit")))
interaction_p <- testModels(mod1,mod2)
tst <- as.tibble(interaction_p$test) %>%
mutate(Outcome = outcome )
int_test[[outcome]] <- tst
}
test_combined <- map_dfr(int_test, bind_rows)
p_load("writexl")
# write_xlsx(test_combined,here::here("Results","All519_unweighted_interactiontests.xlsx"))
# saveRDS( test_combined, here::here("Results","All519_unweighted_interactiontests.rds"))
interaction <- readRDS(here::here("Results","All519_unweighted_interactiontests.rds")) %>%
dplyr::select(Outcome, F.value,df1,df2,`P(>F)`) %>%
mutate(Fvalue = format(round(F.value, digits = 2),nsmall=2),
pvalue = format(round(`P(>F)`,digits=3),nsmall=3)) %>%
mutate(pvalue = ifelse(pvalue == "0.000","<0.001",pvalue)) %>%
dplyr::select(Outcome, Fvalue,df1,df2,pvalue)
#print(knitr::kable(interaction,format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))
5.3 Weighted Interaction tests
int_test <- list()
for (i in 1:length(outcomes)) {
outcome <- outcomes[i]
mod1 <- lapply(imputed_dta, function(x) glmmTMB(as.formula(paste0( outcome, " ~ firePM_cent10per10*Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")),x,weights=Base_sample_weight,
family = binomial(link = "logit")))
mod2 <- lapply(imputed_dta, function(x) glmmTMB(as.formula(paste0( outcome, " ~ firePM_cent10per10 + Time + age_per10 + gender +
education + employmentstatus + Base_IRSAD_Score_per100 + Base_Work_Exposed_2grp +
Base_Asthma_COPD + Base_smokestatus + Base_lnpackyears +
(1 | RESIDENT_ID)")),x,weights=Base_sample_weight,
family = binomial(link = "logit")))
interaction_p <- testModels(mod1,mod2)
tst <- as.tibble(interaction_p$test) %>%
mutate(Outcome = outcome )
int_test[[outcome]] <- tst
}
test_combined <- map_dfr(int_test, bind_rows)
# write_xlsx(test_combined,here::here("Results","All519_weighted_interactiontests.xlsx"))
# saveRDS( test_combined, here::here("Results","All519_weighted_interactiontests.rds"))
interaction <- readRDS(here::here("Results","All519_weighted_interactiontests.rds")) %>%
dplyr::select(Outcome, F.value,df1,df2,`P(>F)`) %>%
mutate(Fvalue = format(round(F.value, digits = 2),nsmall=2),
pvalue = format(round(`P(>F)`,digits=3),nsmall=3)) %>%
mutate(pvalue = ifelse(pvalue == "0.000","<0.001",pvalue)) %>%
dplyr::select(Outcome, Fvalue,df1,df2,pvalue)
#print(knitr::kable(interaction,format="html") %>% kable_styling(full_width = TRUE, bootstrap_options = c("striped")))