partitions vs ordinal metadata

enterotypes = readr::read_csv2("enterotypes_prediction_outliers.csv")[,-1] %>% filter(!is.na(Enterotypes_id))

ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
New names:Rows: 17600 Columns: 5── Column specification ───────────────────────────────────────────────────────────────────────────────────────
Delimiter: ";"
chr (2): sample_name, set
dbl (3): ...1, Enterotypes_id, value
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

metadata=read.csv2(system.file("data-raw/Metadata_10317_20191022-112414_curatedv4_VSv1.csv", package = "gutzilla"), stringsAsFactors = FALSE, ) %>% mutate_all(na_if,"")

metadata %>% .[1:5,1:5]

data("UNSD_countries")

metadata <- metadata %>% 
  #select(SAMPLE_NAME, COUNTRY_OF_BIRTH) %>%
  merge(UNSD_countries %>% 
          select(`Country or Area`, `Sub-region Name`, `Region Name`) %>% 
          mutate(`Country or Area` =  gsub("United Kingdom of Great Britain and Northern Ireland", "United Kingdom", `Country or Area`)) %>% 
          mutate(`Country or Area` =  gsub("United States of America", "USA", `Country or Area`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",`Sub-region Name`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Melanesia", "Micronesia","Polynesia"), "Oceania (others)",`Sub-region Name`)), 
        by.x="COUNTRY_OF_BIRTH", by.y="Country or Area", all.x = TRUE)




variable_to_remove= c("HEIGHT_CM", 
"WEIGHT_KG",
"SAMPLE_TYPE",
"BMI_CAT",
"HOST_SUBJECT_ID",
"ANONYMIZED_NAME",
"ASSIGNED_FROM_GEO",
"BIRTH_YEAR",
"DESCRIPTION",
"DNA_EXTRACTED",
"HEIGHT_UNITS",
"WEIGHT_UNITS",
"SCIENTIFIC_NAME",
"SURVEY_ID",
"TAXON_ID",
"WEIGHT_UNITS",
"VIOSCREEN_STATUS","VIOSCREEN_SURVEYID","PUBLIC","TITLE","GEO_LOC_NAME")




health_metadata = c("BMI_CAT",
"ANTIBIOTIC_HISTORY" ,
"ALZHEIMERS",
"AUTOIMMUNE",
"CANCER",
"DIABETES",
"IBD",
"IBS",
"KIDNEY_DISEASE",
"ACID_REFLUX",
"ALLERGIC_TO_I_HAVE_NO_FOOD_ALLERGIES_THAT_I_KNOW_OF",
"CARDIOVASCULAR_DISEASE",
"CDIFF",
"DEPRESSION_BIPOLAR_SCHIZOPHRENIA",
"EPILEPSY_OR_SEIZURE_DISORDER",
"FUNGAL_OVERGROWTH",
"GLUTEN",
"LACTOSE",
"LIVER_DISEASE",
"LUNG_DISEASE",
"MENTAL_ILLNESS",
"MIGRAINE",
"PKU",
"SIBO",
"THYROID")

diet_variables =
  c("TYPES_OF_PLANTS",
  "DIET_TYPE",
  "ALCOHOL_FREQUENCY",
  "RED_MEAT_FREQUENCY",
  "WHOLE_GRAIN_FREQUENCY",
  "VEGETABLE_FREQUENCY",
  "MEAT_EGGS_FREQUENCY",
  "PREPARED_MEALS_FREQUENCY",
  "FRUIT_FREQUENCY",
  "ONE_LITER_OF_WATER_A_DAY_FREQUENCY",
  "SEAFOOD_FREQUENCY",
  "HOMECOOKED_MEALS_FREQUENCY",
  "HIGH_FAT_RED_MEAT_FREQUENCY",
  "READY_TO_EAT_MEALS_FREQUENCY",
  "SUGAR_SWEETENED_DRINK_FREQUENCY",
  "SUGARY_SWEETS_FREQUENCY",
  "SALTED_SNACKS_FREQUENCY",
  "SPECIALIZED_DIET_FODMAP",
  "SPECIALIZED_DIET_EXCLUDE_DAIRY",
  "SPECIALIZED_DIET_EXCLUDE_NIGHTSHADES",
  "SPECIALIZED_DIET_MODIFIED_PALEO_DIET",
  "SPECIALIZED_DIET_PALEODIET_OR_PRIMAL_DIET",
  "OLIVE_OIL")

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,ANTIBIOTIC_HISTORY),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(ANTIBIOTIC_HISTORY = dplyr::recode_factor(ANTIBIOTIC_HISTORY,   
                             `I have not taken antibiotics in the past year.`="I have not taken antibiotics in the past year.",
                             `Year`="Year", `6 months`="6 months",  
                             `Month`="Month", `Week` = "Week",  
                             .ordered = TRUE)  ) %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = MASS::polr(ANTIBIOTIC_HISTORY ~ ., data = df_test, Hess=TRUE)

summary(polr_model)

Call:
MASS::polr(formula = ANTIBIOTIC_HISTORY ~ ., data = df_test, 
    Hess = TRUE)

Coefficients:
       Value Std. Error t value
`1`  -1.3385     0.1255 -10.667
`2`  -0.9025     0.1180  -7.649
`3`  -1.0740     0.1198  -8.966
`4`  -0.8686     0.1207  -7.197
`5`  -1.3373     0.1234 -10.838
`6`  -0.7667     0.1188  -6.452
`7`  -1.2953     0.1295 -10.005
`8`  -0.6436     0.1216  -5.295
`9`  -0.6512     0.1220  -5.338
`10` -0.2498     0.1181  -2.115
`11` -0.2414     0.1246  -1.938
`12` -0.9527     0.1318  -7.228
`13` -0.6256     0.1247  -5.017
`14` -1.4279     0.1598  -8.938
`15`  0.4592     0.1242   3.698
`16` -0.6848     0.1524  -4.494
`17` -0.9662     0.1637  -5.903
`18` -0.9273     0.1786  -5.191

Intercepts:
                                                    Value    Std. Error t value 
I have not taken antibiotics in the past year.|Year  -0.0475   0.1025    -0.4633
Year|6 months                                         0.7979   0.1027     7.7659
6 months|Month                                        2.1292   0.1062    20.0430
Month|Week                                            3.1388   0.1147    27.3689

Residual Deviance: 30339.40 
AIC: 30383.40 

exp(coef(polr_model))

      `1`       `2`       `3`       `4`       `5`       `6`       `7`       `8`       `9`      `10`      `11` 
0.2622394 0.4055686 0.3416509 0.4195224 0.2625527 0.4645654 0.2738110 0.5253882 0.5214126 0.7789634 0.7854937 
     `12`      `13`      `14`      `15`      `16`      `17`      `18` 
0.3857055 0.5349344 0.2398110 1.5828795 0.5041748 0.3805119 0.3956386 

ctable <- coef(summary(polr_model))

p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

Attente de la réalisation du profilage...

metadata$TYPES_OF_PLANTS %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,TYPES_OF_PLANTS),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(TYPES_OF_PLANTS = dplyr::recode_factor(TYPES_OF_PLANTS,   
                             `Less than 5`="Less than 5",
                             `6 to 10`="6 to 10", `11 to 20`="11 to 20",  
                             `21 to 30`="21 to 30", `More than 30` = "More than 30",  
                             .ordered = TRUE)  ) %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = MASS::polr(TYPES_OF_PLANTS ~ ., data = df_test, Hess=TRUE)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

metadata$ONE_LITER_OF_WATER_A_DAY_FREQUENCY %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,ONE_LITER_OF_WATER_A_DAY_FREQUENCY),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(ONE_LITER_OF_WATER_A_DAY_FREQUENCY = dplyr::recode_factor(ONE_LITER_OF_WATER_A_DAY_FREQUENCY,   
                             `Never`="Never",
                             `Rarely`="Rarely (less than once/week)",
                             `Occasionally`="Occasionally (1-2 times/week)",  
                             `Regularly`="Regularly (3-5 times/week)", 
                             `Daily` = "Daily",  
                             .ordered = TRUE)  ) %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = MASS::polr(ONE_LITER_OF_WATER_A_DAY_FREQUENCY ~ ., data = df_test, Hess=TRUE)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

diet_variables =
  c("TYPES_OF_PLANTS",
  "DIET_TYPE",
  "ALCOHOL_FREQUENCY",
  "RED_MEAT_FREQUENCY",
  "WHOLE_GRAIN_FREQUENCY",
  "VEGETABLE_FREQUENCY",
  "MEAT_EGGS_FREQUENCY",
  "PREPARED_MEALS_FREQUENCY",
  "FRUIT_FREQUENCY",
  "ONE_LITER_OF_WATER_A_DAY_FREQUENCY",
  "SEAFOOD_FREQUENCY",
  "HOMECOOKED_MEALS_FREQUENCY",
  "HIGH_FAT_RED_MEAT_FREQUENCY",
  "READY_TO_EAT_MEALS_FREQUENCY",
  "SUGAR_SWEETENED_DRINK_FREQUENCY",
  "SUGARY_SWEETS_FREQUENCY",
  "SALTED_SNACKS_FREQUENCY",
  "SPECIALIZED_DIET_FODMAP",
  "SPECIALIZED_DIET_EXCLUDE_DAIRY",
  "SPECIALIZED_DIET_EXCLUDE_NIGHTSHADES",
  "SPECIALIZED_DIET_MODIFIED_PALEO_DIET",
  "SPECIALIZED_DIET_PALEODIET_OR_PRIMAL_DIET",
  "OLIVE_OIL")



prepare_data_for_modeling = function(target, levels_set, order_levels) {

  target <- enquo(target)
  #levels_set <- enquo(levels_set)
  
  levels_set=levels_set[levels_set %in% (metadata %>% select(SAMPLE_NAME,!!target) %>% pull(2) %>% unique())]
  
  
  enterotypes %>%
    reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
    merge(metadata %>% select(SAMPLE_NAME,!!target),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
    na.omit() %>%
    dplyr::rename(x := !!target) %>%
    mutate(x = forcats::fct_recode(x, !!!levels_set) %>% 
             factor(levels = order_levels , ordered = TRUE) )  %>%
    select(-SAMPLE_NAME,-`19`) -> df_test
  
return(df_test)
    
}



  
polr_modeling = function(data) {
  
  df_test=data
  
  polr_model = MASS::polr(x ~ ., data = df_test, Hess=TRUE)
  
  #summary(polr_model)
  
  #exp(coef(polr_model))
  
  ctable <- coef(summary(polr_model))
  p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
  ctable <- cbind(ctable, "p value" = p)
  suppressMessages(ci <- confint(polr_model))
  #exp(cbind(OR = coef(polr_model), ci))
  #ctable
  
  res=merge(exp(cbind(OR = coef(polr_model), ci)),ctable,by="row.names", all=TRUE)
  return(res)
}
    

target =
  c("RED_MEAT_FREQUENCY",
  "WHOLE_GRAIN_FREQUENCY",
  "VEGETABLE_FREQUENCY",
  "MEAT_EGGS_FREQUENCY",
  "PREPARED_MEALS_FREQUENCY",
  "FRUIT_FREQUENCY",
  "ONE_LITER_OF_WATER_A_DAY_FREQUENCY",
  "SEAFOOD_FREQUENCY",
  "HOMECOOKED_MEALS_FREQUENCY",
  "HIGH_FAT_RED_MEAT_FREQUENCY",
  "READY_TO_EAT_MEALS_FREQUENCY",
  "SUGAR_SWEETENED_DRINK_FREQUENCY",
  "SUGARY_SWEETS_FREQUENCY",
  "SALTED_SNACKS_FREQUENCY",
  "OLIVE_OIL")


#target = c("RED_MEAT_FREQUENCY","WHOLE_GRAIN_FREQUENCY")




levels_set = c(`Never`="Never",
                             `Rarely`="Rarely (less than once/week)",
                             `Occasionally`="Occasionally (1-2 times/week)",  
                             `Regularly`="Regularly (3-5 times/week)", 
                             `Daily` = "Daily")

order_levels = c("Never","Rarely","Occasionally","Regularly","Daily")


combine_polr_model_result = function(target,levels_set, order_levels){

  final_res = NULL
  
for(i in seq_along(target)) {
  
  cat("modeling",target[i],"...")
  
  res= prepare_data_for_modeling(target = target[i], levels_set=levels_set, order_levels=order_levels) %>% 
  polr_modeling()
  
  tmp=data.frame(target=target[i],res)
  
  final_res = rbind(final_res,tmp)
  
  cat(". done\n") }

  return(final_res) 
    
}


diet_final_res = combine_polr_model_result(target = target, levels_set=levels_set, order_levels=order_levels)

save(diet_final_res, file="diet_final_res.rda")

diet_final_res %>%
  mutate(Row.names = Row.names %>% as.character()) %>% 
  mutate(OR = cut(OR, breaks = c(0,1/4,1/2,1,2,4,100)) %>% as.character()) %>%
  filter(!is.na(OR)) %>% #pull(OR)
  mutate(OR = ifelse(p.value>0.05,"NS",OR)) %>% #pull(OR) %>% factor()
  mutate(OR = factor(OR, levels=c("(0.25,0.5]", "(0.5,1]", "NS", "(1,2]", "(2,4]", "(4,100]" )  )) %>% 
  filter(Row.names %in% paste0("`",as.character(1:20),"`")   ) %>%
  ggplot() + geom_tile(aes(x=Row.names%>%as.character, y=target, fill=OR)) + 
  scale_fill_manual("OR", labels=c("0.25-0.5","0.5-1","NS","1-2","2-4",">4"), 
                    values = c("#d8b365","#f6e8c3","#f5f5f5","#c7eae5","#5ab4ac","#01665e")) + cowplot::theme_cowplot()


diet_final_res %>%
  mutate(Row.names = Row.names %>% as.character()) %>%
  filter(!is.na(OR)) %>%
  mutate(OR=log(OR)) %>%
  reshape2::dcast(Row.names~target, value.var = "OR") %>%
  tibble::column_to_rownames("Row.names") %>%
  as.matrix() %>%
  t() %>%
  heatmap()



diet_final_res %>%
  mutate(Row.names = Row.names %>% as.character()) %>%
  filter(!is.na(OR)) %>%
  mutate(OR=log(OR)) %>%
  group_by(target) %>%
  mutate(OR=scale(OR)) %>%
  ggplot() + geom_tile(aes(x=Row.names%>%as.character, y=target, fill=OR)) + scale_fill_distiller("log\nOdds",type="div")

metadata %>% 
  select(contains("_FREQUENCY"), -OTHER_SUPPLEMENT_FREQUENCY, -FERMENTED_FREQUENCY, OLIVE_OIL) %>%
  colnames() -> target

levels_set = c(A = "Never",
               B = "Rarely (less than once/week)",
               B = "Rarely (less than once/week)",              
               C = "Occasionally (1-2 times/week)",  
               D = "Regularly (3-5 times/week)", 
               E = "Daily")

order_levels=LETTERS[1:5]

#order_levels = c("Never","Rarely","Occasionally","Regularly","Daily")

global_final_res = combine_polr_model_result(target = target[1:3], levels_set=levels_set, order_levels=order_levels)


#prepare_data_for_modeling(target = target[1], levels_set=levels_set, order_levels=order_levels)

metadata$`Sub-region Name` %>% table
metadata$`Region Name` %>% table

enterotypes %>%
  mutate(Enterotypes_id=paste0("et",Enterotypes_id)) %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% filter(COUNTRY %in% c("USA","Canada")) %>% select(SAMPLE_NAME,`Sub-region Name`),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(`Sub-region Name` = forcats::fct_relevel(`Sub-region Name`, ref = "Northern America")) %>%
  dplyr::rename(region=`Sub-region Name`) %>%
  select(-SAMPLE_NAME,-et19) -> df_region
  nnet::multinom(region ~ ., data = df_region) -> region_birth_model


region_birth_model_effects = ggeffects::ggeffect(region_birth_model)

save(region_birth_model, file="region_birth_model.rda")

region_birth_model_summary = summary(region_birth_model)

ci = confint(region_birth_model)

zvalues <- region_birth_model_summary$coefficients / region_birth_model_summary$standard.errors
pnorm(abs(zvalues), lower.tail=FALSE)*2


(100*(exp(coef(region_birth_model))-1)) %>% .[,-1] %>%  heatmap()


coef(region_birth_model) %>% .[,-1] %>% heatmap()

zvalues %>% .[,-1] %>% .[,-13]%>%  heatmap()

coef(region_birth_model) %>% .[,-1] %>% as.data.frame() %>%
  tibble::rownames_to_column("region") %>%
  reshape2::melt(id.vars="region") %>%
  group_by(region) %>%
  mutate(value=scale(value)) %>%
  ggplot() + geom_tile(aes(y=region,x=variable,fill=value)) + scale_fill_distiller(type="div")

rbind(coef(region_birth_model) %>% .[,-1] %>% as.data.frame() %>%
  tibble::rownames_to_column("region") %>%
  reshape2::melt(id.vars="region") %>%
  group_by(region) %>%
  mutate(value=scale(value)) %>% dplyr::rename(target=1,partition=2,OR=3) %>% select(1:3),


diet_final_res %>%
  mutate(Row.names = Row.names %>% as.character()) %>%
  filter(!is.na(OR)) %>%
  mutate(OR=log(OR)) %>%
  group_by(target) %>%
  mutate(OR=scale(OR))  %>% dplyr::rename(target=1,partition=2,OR=3) %>% select(1:3) ) %>%
  
  reshape2::dcast(partition~target, value.var = "OR") %>%
  tibble::column_to_rownames("partition") %>%
  as.matrix() %>%
  t() %>%
  heatmap(col = hcl.colors(12, "BrBG", rev = TRUE))

metadata$GLUTEN %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,GLUTEN),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(GLUTEN = forcats::fct_recode(GLUTEN,   
                             `No`="No",
                             `Yes`="I do not eat gluten because it makes me feel bad",
                             `Yes`="I was diagnosed with celiac disease",  
                             `medium`="I was diagnosed with gluten allergy (anti-gluten IgG), but not celiac disease") %>% 
           factor(levels = c("No","medium","Yes") , ordered = TRUE) )  %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = MASS::polr(GLUTEN ~ ., data = df_test, Hess=TRUE)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

metadata$GLUTEN %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,GLUTEN),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(GLUTEN = forcats::fct_recode(GLUTEN,   
                             `No`="No",
                             `Yes`="I do not eat gluten because it makes me feel bad",
                             `Yes`="I was diagnosed with celiac disease",  
                             `Yes`="I was diagnosed with gluten allergy (anti-gluten IgG), but not celiac disease") %>% 
           factor(levels = c("No","Yes") , ordered = TRUE) )  %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = glm(GLUTEN ~ ., data = df_test, family=binomial)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

#p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
#ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

metadata$SPECIALIZED_DIET_EXCLUDE_NIGHTSHADES %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,SPECIALIZED_DIET_EXCLUDE_DAIRY),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(SPECIALIZED_DIET_EXCLUDE_DAIRY = forcats::fct_recode(SPECIALIZED_DIET_EXCLUDE_DAIRY,   
                             `No`="No",
                             `Yes`="Yes",
                             `Yes`="true",  
                             `No`="false") %>% 
           factor(levels = c("No","Yes") , ordered = TRUE) )  %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = glm(SPECIALIZED_DIET_EXCLUDE_DAIRY ~ ., data = df_test, family=binomial)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

#p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
#ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

metadata$DIABETES %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,DIABETES),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(DIABETES = forcats::fct_recode(DIABETES,   
                             `No`="I do not have this condition",
                             `Yes`="Diagnosed by a medical professional (doctor, physician assistant)",
                             `medium`="Diagnosed by an alternative medicine practitioner",  
                             `medium`="Self-diagnosed") %>% 
           factor(levels = c("No","medium","Yes") , ordered = TRUE) )  %>%
  filter(DIABETES != "medium") %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = glm(DIABETES ~ ., data = df_test, family=binomial)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

#p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
#ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

metadata$BMI_CAT %>% table

enterotypes %>%
  reshape2::dcast(`sample_name` ~ `Enterotypes_id`, length) %>%
  merge(metadata %>% select(SAMPLE_NAME,BMI_CAT),., by.y="sample_name", by.x="SAMPLE_NAME") %>%
  na.omit() %>%
  mutate(BMI_CAT = forcats::fct_recode(BMI_CAT,   
                             `Underweight`="Underweight",
                             `Normal`="Normal",
                             `Overweight`="Overweight",  
                             `Obese`="Obese") %>% 
           factor(levels = c("Underweight","Normal","Overweight","Obese") , ordered = TRUE) )  %>%
  filter(BMI_CAT != "Underweight") %>%
  mutate(BMI_CAT = factor(BMI_CAT, levels = c("Normal","Overweight","Obese") , ordered = TRUE)) %>%
  select(-SAMPLE_NAME,-`19`) -> df_test
  
polr_model = MASS::polr(BMI_CAT ~ ., data = df_test, Hess=TRUE)

summary(polr_model)

exp(coef(polr_model))

ctable <- coef(summary(polr_model))

p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
ctable <- cbind(ctable, "p value" = p)

ci <- confint(polr_model)

exp(cbind(OR = coef(polr_model), ci))
ctable

multinomial regression

we aimed here to predict microbiome partition from metadata. Partition #1, the most centroid in branch structure and highest diversity, will be used as reference

prepare data

import metadata and merge with UNSD regions

enterotypes = readr::read_csv2("enterotypes_prediction_outliers.csv")[,-1] %>% filter(!is.na(Enterotypes_id))
metadata=read.csv2(system.file("data-raw/Metadata_10317_20191022-112414_curatedv4_VSv1.csv", package = "gutzilla"), stringsAsFactors = FALSE, ) %>% mutate_all(na_if,"")

metadata %>% .[1:5,1:5]


data("UNSD_countries")

metadata <- metadata %>% 
  #select(SAMPLE_NAME, COUNTRY_OF_BIRTH) %>%
  merge(UNSD_countries %>% 
          select(`Country or Area`, `Sub-region Name`, `Region Name`) %>% 
          mutate(`Country or Area` =  gsub("United Kingdom of Great Britain and Northern Ireland", "United Kingdom", `Country or Area`)) %>% 
          mutate(`Country or Area` =  gsub("United States of America", "USA", `Country or Area`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Central Asia", "Southern Asia"), "Central and Southern Asia",`Sub-region Name`)) %>%
          mutate(`Sub-region Name` =  ifelse(`Sub-region Name` %in% c("Melanesia", "Micronesia","Polynesia"), "Oceania (others)",`Sub-region Name`)), 
        by.x="COUNTRY_OF_BIRTH", by.y="Country or Area", all.x = TRUE)

clean metadata enterotypes

metadata %>%
  select(SAMPLE_NAME, contains("_FREQUENCY"), 
         TYPES_OF_PLANTS, OLIVE_OIL, DIET_TYPE, 
         contains("SPECIALIZED_DIET_"),
         all_of(health_metadata), contains("ALCOHOL_TYPES_"), 
         AGE_YEARS, SEX, COLLECTION_SEASON, `Sub-region Name`) %>%
  select(-SPECIALIZED_DIET_WESTENPRICE_OR_OTHER_LOWGRAIN_LOW_PROCESSED_FO, -ALCOHOL_TYPES_UNSPECIFIED) %>%
  dplyr::rename(REGION_BIRTH=`Sub-region Name`) %>%
  reshape2::melt(id.vars="SAMPLE_NAME") %>% 
  filter(value!="LabControl test") %>%
  filter(value != "Self-diagnosed") %>%
  filter(value != "Diagnosed by an alternative medicine practitioner") %>%
  mutate(value = forcats::fct_recode(value, Never = "Never",
                                              Rarely = "Rarely (less than once/week)",
                                              Rarely = "Rarely (a few times/month)",              
                                              Occasionally = "Occasionally (1-2 times/week)",                 
                                              Regularly = "Regularly (3-5 times/week)", 
                                              Daily = "Daily",
                                              Yes = "true",
                                              No  = "false",
                                              Yes = "Diagnosed by a medical professional (doctor, physician assistant)",
                                              No  = "I do not have this condition" )) %>%
  mutate(value = forcats::fct_explicit_na(value, "Missing") ) %>%
  filter(!is.na(value)) %>%
  mutate(value= value %>% as.character) %>%
  reshape2::dcast(SAMPLE_NAME~variable, value.var = "value", fill="Missing") %>%
  mutate(AGE_YEARS = AGE_YEARS %>% as.character() %>% as.numeric) %>%
  mutate(BOWEL_MOVEMENT_FREQUENCY = forcats::fct_relevel(BOWEL_MOVEMENT_FREQUENCY,ref="One")) %>%
  mutate(REGION_BIRTH             = forcats::fct_relevel(REGION_BIRTH, ref = "Northern America")) %>%
  mutate(ANTIBIOTIC_HISTORY       = forcats::fct_relevel(ANTIBIOTIC_HISTORY, ref = "I have not taken antibiotics in the past year.")) %>%
  merge(enterotypes %>% select(sample_name, Enterotypes_id), by.y="sample_name", by.x="SAMPLE_NAME") %>%
  mutate(Enterotypes_id = paste0("m",Enterotypes_id)) %>% select(-SAMPLE_NAME) -> metadata_enterotypes

Avis : NAs introduced by coercion

clean metadata enterotypes numerical

range0_1 = function(x) {(x-min(x, na.rm = TRUE))/(max(x, na.rm="TRUE")-min(x,na.rm = TRUE))}

metadata %>%
  select(SAMPLE_NAME, contains("_FREQUENCY"), 
         TYPES_OF_PLANTS, OLIVE_OIL, DIET_TYPE, 
         contains("SPECIALIZED_DIET_"),
         all_of(health_metadata),BMI, contains("ALCOHOL_TYPES_"),
         AGE_YEARS, SEX, COLLECTION_SEASON, `Sub-region Name`) %>%
  select(-SPECIALIZED_DIET_WESTENPRICE_OR_OTHER_LOWGRAIN_LOW_PROCESSED_FO, -BMI_CAT, -ALCOHOL_TYPES_UNSPECIFIED) %>%
  dplyr::rename(REGION_BIRTH=`Sub-region Name`) %>% #select(-SAMPLE_NAME) %>% apply(2,table)
  tidyext::onehot(var = c("REGION_BIRTH","GLUTEN","COLLECTION_SEASON")) %>%
  select(-GLUTEN_No) %>%
  reshape2::melt(id.vars=c("SAMPLE_NAME")) %>% 
  mutate(value= ifelse(value %in% c("LabControl test",
                                    "other",
                                    "Self-diagnosed",
                                    "Diagnosed by an alternative medicine practitioner"),NA,value)) %>%
  mutate(value = forcats::fct_recode(value,   `0` = "Never",
                                              `0.25` = "Rarely (less than once/week)",
                                              `0.25` = "Rarely (a few times/month)",              
                                              `0.5` = "Occasionally (1-2 times/week)",                 
                                              `0.75` = "Regularly (3-5 times/week)", 
                                              `1` = "Daily",
                                              `1` = "true",
                                              `0` = "false",
                                              `1` = "Yes",
                                              `0` = "No",
                                              `1` = "Diagnosed by a medical professional (doctor, physician assistant)",
                                              `0` = "I do not have this condition",
                                              `0` = "female",
                                              `1` = "male",
                                              `0` = "I have not taken antibiotics in the past year.",
                                              `0.25` = "Year",
                                              `0.5` = "6 months",
                                              `0.75` = "Month",
                                              `1` = "Week",
                                              `0` = "Less than 5",
                                              `0.25` = "6 to 10",
                                              `0.5` = "11 to 20",              
                                              `0.75` = "21 to 30",                 
                                              `1` = "More than 30", 
                                              `0` = "Less than one",
                                              `0.2` = "One",
                                              `0.4` = "Two",
                                              `0.6` = "Three",              
                                              `0.6` = "Four",                 
                                              `1` = "Five or more",
                                              `0` = "Omnivore" ,
                                              `0.25` = "Omnivore but do not eat red meat",
                                              `0.5` = "Vegetarian but eat seafood",              
                                              `0.75` = "Vegetarian",                 
                                              `1` = "Vegan"
                                    
                                       
                                       
                                       )) %>%
  mutate(value= value %>% as.character %>% as.numeric) %>%
  group_by(variable) %>%
  mutate(value = ifelse(is.na(value),mean(value,na.rm = TRUE),value)) %>% #slice_sample(prop=0.1)
  reshape2::dcast(SAMPLE_NAME~variable, value.var = "value", fill="Missing") %>%
  mutate(AGE_YEARS = AGE_YEARS %>% as.character() %>% as.numeric %>% range0_1() ) %>%
  mutate(BMI = BMI %>% as.character() %>% as.numeric %>% range0_1() ) %>%
  #mutate(AGE_YEARS = ifelse(is.na(AGE_YEARS),mean(AGE_YEARS,na.rm = TRUE),value)) %>%
  #mutate(REGION_BIRTH             = forcats::fct_relevel(REGION_BIRTH, ref = "Northern America")) %>%
  merge(enterotypes %>% select(sample_name, Enterotypes_id), by.y="sample_name", by.x="SAMPLE_NAME") %>%
  mutate(Enterotypes_id = paste0("m",Enterotypes_id)) %>% tibble::column_to_rownames("SAMPLE_NAME") -> metadata_enterotypes_num

metadata_enterotypes_num %>%
  select(-Enterotypes_id) %>%
  mutate_all(as.numeric) %>% cor(method="spearman") %>%
  .^2 %>% sqrt() %>% heatmap

NA
NA

metadata enterotypes correspondance analysis

metadata_enterotypes_num %>%
  select(-Enterotypes_id) %>%
  mutate_all(as.numeric) %>%
  ade4::dudi.pca(scannf=F,nf=30) -> metadata_enterotypes_num_pca


cumsum(metadata_enterotypes_num_pca$eig/sum(metadata_enterotypes_num_pca$eig))

metadata_enterotypes_num_pca$co %>% .[1:30] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  ggplot() + geom_bar(aes(x=metadata,y=value), stat = "identity") + 
  facet_wrap(~variable, ncol = 10) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))



metadata_enterotypes_num_pca$co %>% .[1:30] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata") %>%
  group_by(variable) %>%
  mutate(value_w=abs(value)) %>%
  group_by(variable) %>%
  mutate(value_tot = sum(value_w)) %>%
  group_by(variable) %>%
  mutate(value_w=value_w/value_tot) %>%
  group_by(variable) %>%
  mutate(value_tot = sum(value_tot)) %>%
  arrange(variable,desc(value_w)) %>%
  group_by(variable) %>%
  mutate(value_cum = cumsum(value_w)) %>%
  filter(value_cum < 0.25) %>%
  pull(metadata) %>% unique()
 
metadata_enterotypes_num_pca$co %>% .[1:30] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata") %>%
  group_by(variable) %>%
  mutate(value_w=abs(value)) %>%
  group_by(variable) %>%
  mutate(value_tot = sum(value_w)) %>%
  group_by(variable) %>%
  mutate(value_w=value_w/value_tot) %>%
  group_by(variable) %>%
  mutate(value_tot = sum(value_tot)) %>%
  arrange(variable,desc(value_w)) %>%
  group_by(variable) %>%
  mutate(value_cum = cumsum(value_w)) %>%
  group_by(variable) %>%
  slice_max(abs(value), n = 2) %>%
  pull(metadata) %>% unique()




  metadata_enterotypes_num_bca= ade4::bca(metadata_enterotypes_num_pca, fac=metadata_enterotypes_num$Enterotypes_id %>% as.character() %>% as.factor(), scannf=F, n=30) 

  #ade4::randtest(metadata_enterotypes_num_bca)
  
  
  metadata_enterotypes_num_bca$co %>%
    tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata") %>%
  group_by(variable) %>%
  mutate(value_w=abs(value)) %>%
  group_by(variable) %>%
  mutate(value_tot = sum(value_w)) %>%
  group_by(variable) %>%
  mutate(value_w=value_w/value_tot) %>%
  group_by(variable) %>%
  mutate(value_tot = sum(value_tot)) %>%
  arrange(variable,desc(value_w)) %>%
  group_by(variable) %>%
  mutate(value_cum = cumsum(value_w)) %>%
  group_by(variable) %>%
  slice_max(abs(value), n = 5) %>%
  #filter(value_cum < 0.05) %>%   
  #filter(abs(value)> 0.1) %>%  
  pull(metadata) %>% unique()
    
  
dim(metadata_enterotypes_num)  

cbind(metadata_enterotypes_num %>%
  select(-DIET_TYPE,-COLLECTION_SEASON,-REGION_BIRTH,-Enterotypes_id,-GLUTEN) %>%
  mutate_all(as.numeric),

metadata_enterotypes_num %>%
  select(DIET_TYPE,COLLECTION_SEASON,REGION_BIRTH,Enterotypes_id,GLUTEN) %>%
  mutate_all(as.character) %>%
  mutate_all(function(x){ifelse(is.na(x),"Missing",x)}) %>%
  mutate_all(as.factor)) %>%
  select(-Enterotypes_id) %>%
  ade4::dudi.mix(scannf=FALSE, nf=32) -> metadata_enterotypes_coa

metadata_enterotypes_num %>%
  select(-DIET_TYPE,-COLLECTION_SEASON,-REGION_BIRTH,-Enterotypes_id,-GLUTEN) %>%
  mutate_all(as.numeric) %>%
  ade4::dudi.mix(scannf=FALSE, nf=30) -> metadata_enterotypes_pca


ade4::scatter(metadata_enterotypes_coa)


cumsum(metadata_enterotypes_coa$eig/sum(metadata_enterotypes_coa$eig))

ade4::s.corcircle(metadata_enterotypes_coa$co)


metadata_enterotypes_coa$co %>% .[1:4] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  ggplot() + geom_bar(aes(x=metadata,y=value), stat = "identity") + 
  facet_wrap(~variable, ncol = 4) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))



metadata_enterotypes_coa$li %>%
  cbind(metadata_enterotypes_num) %>%
  ggplot() + geom_point(aes(x=Axis1, y=Axis2, col=Enterotypes_id))


metadata_enterotypes_coa$co %>% .[1:32] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  filter(abs(value)>0.5)


cumsum(metadata_enterotypes_pca$eig/sum(metadata_enterotypes_pca$eig))

metadata_enterotypes_pca$c1 %>% .[1:24] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  mutate(type=NA) %>%
  mutate(type=ifelse(grepl("REGIO", metadata),"Region",type )) %>%
  mutate(type=ifelse(grepl("GLUTE", metadata),"Gluten",type )) %>%
  mutate(type=ifelse(grepl("COLLE", metadata),"Collection",type )) %>%
  mutate(type=ifelse(grepl("DIET", metadata),"Diet",type )) %>%
  filter(is.na(type)) %>%
  filter(abs(value)>0.25) %>%
  ggplot() + geom_bar(aes(x=metadata,y=value), stat = "identity") + 
  facet_wrap(~variable, ncol = 6, scales = "free_x") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))


metadata_enterotypes_pca$li %>%
  cbind(metadata_enterotypes_num) %>%
  ggplot() + geom_point(aes(x=Axis1, y=Axis2, col=Enterotypes_id))

metadata_enterotypes_pca$c1 %>% .[1:23] %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  filter(abs(value)>0.25)


metadata_enterotypes_bca = ade4::bca(metadata_enterotypes_coa, fac=as.factor(metadata_enterotypes_num$Enterotypes_id), nf=18, scannf=FALSE)

metadata_enterotypes_bca_randtest = ade4::randtest(metadata_enterotypes_bca)

cumsum(metadata_enterotypes_bca$eig/sum(metadata_enterotypes_bca$eig))


metadata_enterotypes_bca$c1 %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  mutate(type=NA) %>%
  mutate(type=ifelse(grepl("REGIO", metadata),"Region",type )) %>%
  mutate(type=ifelse(grepl("GLUTE", metadata),"Gluten",type )) %>%
  mutate(type=ifelse(grepl("COLLE", metadata),"Collection",type )) %>%
  mutate(type=ifelse(grepl("DIET", metadata),"Diet",type )) %>%
  filter(is.na(type)) %>%
  filter(abs(value)>0.25) %>%
  ggplot() + geom_bar(aes(x=metadata,y=value), stat = "identity") + 
  facet_wrap(~variable, ncol = 6, scales = "free_x") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

metadata_enterotypes_bca$c1 %>%
  tibble::rownames_to_column("metadata") %>%
  reshape2::melt(id.vars="metadata")  %>%
  mutate(type=NA) %>%
  mutate(type=ifelse(grepl("REGIO", metadata),"Region",type )) %>%
  mutate(type=ifelse(grepl("GLUTE", metadata),"Gluten",type )) %>%
  mutate(type=ifelse(grepl("COLLE", metadata),"Collection",type )) %>%
  mutate(type=ifelse(grepl("DIET", metadata),"Diet",type )) %>%
  filter(!is.na(type)) %>%
  filter(abs(value)>1)

ade4::s.class(metadata_enterotypes_bca$ls, fac=as.factor(metadata_enterotypes_num$Enterotypes_id))

ade4::s.label(metadata_enterotypes_bca$li)

alcool types explo

library(UpSetR)

alcohol_types_list = list(BEERCIDER   = metadata %>% filter(ALCOHOL_TYPES_BEERCIDER == "Yes") %>% pull(SAMPLE_NAME),
                          SPIRITSHARD = metadata %>% filter(ALCOHOL_TYPES_SPIRITSHARD_ALCOHOL == "Yes") %>% pull(SAMPLE_NAME),
                          RED_WINE    = metadata %>% filter(ALCOHOL_TYPES_RED_WINE == "Yes") %>% pull(SAMPLE_NAME),
                          SOUR_BEERS  = metadata %>% filter(ALCOHOL_TYPES_SOUR_BEERS == "Yes") %>% pull(SAMPLE_NAME),
                          WHITE_WINE  = metadata %>% filter(ALCOHOL_TYPES_WHITE_WINE == "Yes") %>% pull(SAMPLE_NAME),
                          UNSPECIFIED = metadata %>% filter(ALCOHOL_TYPES_UNSPECIFIED  == "Yes") %>% pull(SAMPLE_NAME))

upset(alcohol_types_list)

metadata %>%
  select(contains("ALCOHOL_TYPES_"), -ALCOHOL_TYPES_UNSPECIFIED) %>%
  mutate_all(function(x) ifelse(x=="Yes",1,0)) %>% 
upset(group.by = "sets",cutoff = 3)

build model

library(nnet)
library(GGally)
library(broom.helpers)

metadata_enterotypes_tmp =  
  metadata_enterotypes_num %>%
  mutate(ANTIBIOTIC_HISTORY = ANTIBIOTIC_HISTORY %>% as.numeric,
         BOWEL_MOVEMENT_FREQUENCY = BOWEL_MOVEMENT_FREQUENCY %>% as.numeric,
         AGE_YEARS = AGE_YEARS %>% as.numeric,
         SEX = SEX %>% as.numeric,
         TYPES_OF_PLANTS = TYPES_OF_PLANTS %>% as.numeric,
         VEGETABLE_FREQUENCY = VEGETABLE_FREQUENCY %>% as.numeric,
         IBD = IBD %>% as.numeric,
         FROZEN_DESSERT_FREQUENCY = FROZEN_DESSERT_FREQUENCY %>% as.numeric,
         TEETHBRUSHING_FREQUENCY = TEETHBRUSHING_FREQUENCY %>% as.numeric) %>%
  #slice_sample(n = 5000, replace = FALSE) %>%
  na.omit()
  


regm  <- multinom(Enterotypes_id ~ ANTIBIOTIC_HISTORY + BOWEL_MOVEMENT_FREQUENCY + 
                    AGE_YEARS + SEX + TYPES_OF_PLANTS + VEGETABLE_FREQUENCY + 
                    IBD + FROZEN_DESSERT_FREQUENCY +  TEETHBRUSHING_FREQUENCY + REGION_BIRTH + GLUTEN , data = metadata_enterotypes_tmp) #MaxNWts = 1026
regm2 <- step(regm)


save(regm, file="regm.rda")
save(regm2, file="regm2.rda")

#regm2_OR   <- Rfast::odds.ratio(regm2)
#save(regm2_OR, file="regm2_OR.rda")

regm2_tidy <- tidy_plus_plus(regm2, exponentiate = TRUE)
regm2_plot <- ggcoef_multinom(regm2,exponentiate = TRUE, type = "faceted")

save(regm2_tidy, file="regm2_tidy.rda")
save(regm2_plot, file="regm2_plot.rda")


#regm2_plot$data %>% ggcoef_plot(facet_col="y.level", exponentiate = TRUE)

library(nnet)
library(GGally)
library(broom.helpers)

metadata_enterotypes_tmp =  
  metadata_enterotypes_num %>%
  mutate(Enterotypes_id = Enterotypes_id %>% forcats::fct_relevel(ref = "m1")) %>%
  mutate_if(is.character, as.numeric) %>%
  na.omit()
  


regm  <- multinom(Enterotypes_id ~ . , data = metadata_enterotypes_tmp, MaxNWts = 1938) 

regm_tidy <- tidy_plus_plus(regm, exponentiate = TRUE)

save(regm_tidy, file="regm_tidy.rda")

# regm2 <- step(regm)
# 
# save(regm, file="regm.rda")
# save(regm2, file="regm2.rda")
# 
# 
# regm2_tidy <- tidy_plus_plus(regm2, exponentiate = TRUE)
# 
# save(regm2_tidy, file="regm2_tidy.rda")
# 
# 
# 
# exp(coef(regm))
# 
# ctable <- coef(summary(regm))
# 
# p <- pnorm(abs(ctable[, "t value"]), lower.tail = FALSE) * 2
# ctable <- cbind(ctable, "p value" = p)

shannon_path = system.file("data-raw/qiime/generated-files-20190512/alpha/shannon.qza", package = "gutzilla")

shannon = qiime2R::read_qza(shannon_path)$data %>% as.data.frame

enterotypes = readr::read_csv2("enterotypes_prediction_outliers.csv")[,-1] %>% filter(!is.na(Enterotypes_id))

ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
New names:Rows: 17600 Columns: 5── Column specification ───────────────────────────────────────────────────────────────────────────────────────
Delimiter: ";"
chr (2): sample_name, set
dbl (3): ...1, Enterotypes_id, value
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

merge(enterotypes,shannon,by.x="sample_name", by.y="row.names") %>%
  group_by(Enterotypes_id) %>%
  summarise(shannon_median=median(shannon)) %>%
  mutate(Enterotypes_id = paste0("m",Enterotypes_id))-> enterotypes_shannon_median

coef(regm) %>% 
  reshape2::melt() %>% 
  filter(Var2 != "(Intercept)") %>% 
  merge(enterotypes_shannon_median, by.x="Var1", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(Var1 %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(Var1 %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(Var1 %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  mutate(Var1 = Var1 %>% forcats::fct_reorder(shannon_median)) %>% 
  ggplot() + geom_point(aes(x=Var1,y=10^value, col=branch)) + scale_y_log10() + facet_wrap(~Var2, scale="free_y")

regm_tidy %>% 
  merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  mutate(y.level = y.level %>% forcats::fct_reorder(shannon_median)) %>%
  ggplot() + geom_point(aes(x=y.level,y=estimate, col=branch, alpha=p.value<0.05)) + scale_y_continuous(trans = "log") + facet_wrap(~term, scale="free_y")

regm_tidy %>% 
  merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  mutate(y.level = y.level %>% forcats::fct_reorder(shannon_median)) %>%
  filter(p.value < 0.05) %>%
  ggplot() + geom_point(aes(x=term,y=log(estimate), col=branch)) + 
  #scale_y_continuous(trans = "log") + 
  facet_wrap(~y.level, scale="free_x", nrow=2) + coord_flip() + theme_classic() + geom_hline(yintercept = 0)

NA
NA

data viz

# regm_tidy %>% 
#   merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
#   mutate(branch ="Clostridiales DMM types") %>%
#   mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
#   mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
#   mutate(branch = ifelse(y.level %in% c("m13"), "Akkermansia DMM types",branch)) %>%
#   mutate(estimate = log(estimate)) %>%
#   #mutate(estimate = ifelse(p.value>0.05,0,estimate)) %>%
#   group_by(variable) %>%
#   mutate(s=mean(abs(estimate))) %>% 
#   filter(s>0) %>%
#   mutate(y.level = y.level %>% forcats::fct_reorder(shannon_median)) %>%
#   mutate(variable = variable %>% forcats::fct_reorder(s)) %>%
#   ggplot() + geom_tile(aes(x=y.level, y=variable, fill=estimate)) + scale_fill_gradient2() + facet_grid(.~branch, scale="free_x", space = "free_x")

load("regm_tidy.rda")  
  
  regm_tidy %>%
    mutate(variable_group = "disease") %>%
    mutate(variable_group = 
             case_when(grepl("FREQUENCY", variable) ~ "Diet, demographic and lifestyle",
                       grepl("ALCOHOL", variable) ~ "Diet, demographic and lifestyle",
                       variable %in% c("AGE_YEARS","SEX","DIET_TYPE","TYPES_OF_PLANTS") ~ "Diet, demographic and lifestyle",
                       grepl("SPECIALIZED", variable) ~"Specialized Diet",
                       grepl("REGION", variable) ~ "Region of Birth",
                       TRUE ~ "Health associated"
                            )) %>%
    merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM partitions") %>%
  mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM partitions",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM partitions",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m13"), "Akk.",branch)) %>%
  mutate(estimate = log(estimate)) %>%
  mutate(estimate2 = ifelse(p.value>0.05,0,estimate)) %>%
  group_by(variable) %>%
  mutate(s=mean(abs(estimate2))) %>% 
  filter(s>0) %>%
  ungroup() %>%
  mutate(y.level = y.level %>% forcats::fct_reorder(shannon_median, .desc = FALSE)) %>%
  mutate(variable = variable %>% as.character %>% as.factor() %>% forcats::fct_reorder(s, .desc = FALSE)) %>% 
  ggplot() + 
    geom_tile(aes(x=y.level, y=variable, fill=estimate)) + 
    geom_point(aes(x=y.level, y=variable, alpha=p.value<0.05)) + 
    scale_fill_gradient2("Log OR", low=scales::muted("blue"), mid="white", high=scales::muted("red")) +
    facet_grid(variable_group~branch, scale="free", space = "free") + 
    cowplot::theme_cowplot() +
    theme(axis.title.x=element_blank(),
        axis.text.x=element_blank()) +
  ylab("Microbiome partition\nordered by branch and by alpha-diversity") + ylab("AGP questionnaire variable") -> fig4A
    fig4A

  ggsave("partitions_multinomial_logOR.pdf", height = 12, width = 14)
  
  
  # reshape2::dcast(y.level~variable, value.var = "estimate") %>%
  # tibble::column_to_rownames("y.level") %>%
  # as.matrix() %>%
  # t() %>%
  # heatmap(scale = "none", col = cm.colors(256))

  
regm_tidy %>%
    mutate(variable_group = "disease") %>%
    mutate(variable_group = 
             case_when(grepl("FREQUENCY", variable) ~ "Diet, demographic and lifestyle",
                       grepl("ALCOHOL", variable) ~ "Diet, demographic and lifestyle",
                       variable %in% c("AGE_YEARS","SEX","DIET_TYPE","TYPES_OF_PLANTS") ~ "Diet, demographic and lifestyle",
                       grepl("SPECIALIZED", variable) ~"Specialized Diet",
                       grepl("REGION", variable) ~ "Region of Birth",
                       TRUE ~ "Health associated"
                            )) %>%
    merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m13"), "Akk.",branch)) %>%
  mutate(estimate = log(estimate)) %>%
  mutate(estimate2 = ifelse(p.value>0.05,0,estimate)) %>%
  group_by(variable) %>%
  mutate(s=mean(abs(estimate2))) %>% 
  filter(s>0) %>%
  ungroup() %>%
  mutate(y.level = y.level %>% forcats::fct_reorder(shannon_median, .desc = FALSE)) %>%
  mutate(variable = variable %>% as.character %>% as.factor() %>% forcats::fct_reorder(s, .desc = FALSE)) %>%  
  group_by(variable,variable_group,branch) %>%
  summarise(estimate=sum(abs(estimate2))) %>%
  ggplot() +
  geom_bar(aes(x=variable,y=estimate,fill=branch), stat="identity", position = "stack") +
  facet_grid(variable_group~., scale="free", space = "free") + coord_flip() + scale_fill_brewer(type="qual") +
  xlab("") + ylab("absolute cummulated Log OR") +
  cowplot::theme_cowplot() +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank()) -> fig4B

`summarise()` has grouped output by 'variable', 'variable_group'. You can override using the `.groups` argument.

  





 regm_tidy %>%
    mutate(variable_group = "disease") %>%
    mutate(variable_group = 
             case_when(grepl("FREQUENCY", variable) ~ "Diet, demographic and lifestyle",
                       grepl("ALCOHOL", variable) ~ "Diet, demographic and lifestyle",
                       variable %in% c("AGE_YEARS","SEX","DIET_TYPE","TYPES_OF_PLANTS") ~ "Diet, demographic and lifestyle",
                       grepl("SPECIALIZED", variable) ~"Specialized Diet",
                       grepl("REGION", variable) ~ "Region of Birth",
                       TRUE ~ "Health associated"
                            )) %>%
    merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM partitions") %>%
  mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM partitions",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM partitions",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m13"), "Akk.",branch)) %>%
  mutate(estimate = log(estimate)) %>%
  mutate(estimate2 = ifelse(p.value>0.05,0,estimate)) %>%
  group_by(variable) %>%
  mutate(s=mean(abs(estimate2))) %>% 
  filter(s>0) %>%
  ungroup() %>%
  mutate(y.level = y.level %>% forcats::fct_reorder(shannon_median, .desc = FALSE)) %>%
  mutate(variable = variable %>% as.character %>% as.factor() %>% forcats::fct_reorder(s, .desc = FALSE)) %>%  
  group_by(y.level, variable_group, branch) %>%
  summarise(estimate=sum(abs(estimate2))) %>%
  ggplot() +
  geom_bar(aes(x=y.level,y=estimate,fill=variable_group), stat="identity", position = "stack") +
  facet_grid(.~branch, scale="free", space = "free") + scale_fill_brewer("",type="qual", palette = "Set2") +
  xlab("") + ylab("absolute cummulated Log OR") + 
   cowplot::theme_cowplot() -> fig4C

`summarise()` has grouped output by 'y.level', 'variable_group'. You can override using the `.groups` argument.

focus on Bacteroides branches

  regm_tidy %>%
  group_by(variable) %>% mutate(p.value = p.adjust(p.value, method="fdr")) %>%
    mutate(variable_group = "disease") %>%
   mutate(variable_group = 
             case_when(variable %in% c("AGE_YEARS","SEX","BMI","BOWEL_MOVEMENT_FREQUENCY") ~ "Intrinsic\nHost",
                       grepl("MEALS", variable) ~ "Lifestyle\n& hygiene",
                       variable %in% c("EXERCISE_FREQUENCY","POOL_FREQUENCY","COSMETICS_FREQUENCY","PREPARED_MEALS_FREQUENCY","TEETHBRUSHING_FREQUENCY","SMOKING_FREQUENCY") ~ "Lifestyle\n& hygiene",
                       grepl("FREQUENCY", variable) ~ "Diet",
                       grepl("ALCOHOL", variable) ~ "Diet",
                       variable %in% c("DIET_TYPE","TYPES_OF_PLANTS") ~ "Diet",
                       grepl("SPECIALIZED", variable) ~"Spec.\nDiet",
                       grepl("REGION", variable) ~ "Reg. of\nBirth",
                       TRUE ~ "Medical history"
                            )) %>%
    merge(enterotypes_shannon_median, by.x="y.level", by.y="Enterotypes_id") %>%
  mutate(branch ="Clostridiales DMM partitions") %>%
  mutate(branch = ifelse(y.level %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM partitions",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM partitions",branch)) %>%
  mutate(branch = ifelse(y.level %in% c("m13"), "Akk.",branch)) %>%
  mutate(substates = "High") %>%
  mutate(substates = ifelse(y.level %in% c("m11","m6","m10","m3"), "Medium",substates)) %>%
  mutate(substates = ifelse(y.level %in% c("m15"), "Low",substates)) %>%
  filter(branch == "Bacteroides DMM partitions") %>%
  mutate(substates = substates %>% forcats::fct_relevel("High",after=2 )) %>%
  mutate(estimate = log(estimate)) %>%
  mutate(estimate2 = ifelse(p.value>p_cut_off,0,estimate)) %>%
  group_by(variable) %>%
  mutate(s=mean(abs(estimate2))) %>% 
  filter(s>0) %>%
  ungroup() %>%
  mutate(y.level = y.level %>% stringr::str_to_upper() %>% forcats::fct_reorder(shannon_median, .desc = TRUE)) %>%
  
  mutate(branch = branch %>% gsub("DMM types", "\nDMM types",.)) %>% #pull(variable) %>% grep("COLL",.)
  filter(!(variable %in% grep("COLL", variable, value = TRUE) )) %>%
  
  mutate(variable = gsub("_FREQUENCY|SPECIALIZED_DIET_|ALCOHOL_TYPES_|REGION_BIRTH_","", variable)) %>%
  mutate(variable = 
           case_when(grepl("NO_FOOD", variable)~"No food allergy",
                     grepl("not.eat.gluten", variable)~"Gluten intolerance",
                     grepl("gluten.IgG", variable)~"Gluten allergy (IgG)",
                     grepl("WESTENPRICE", variable)~"Low processed",
                     grepl("OTHER_RESTRICTIONS", variable)~"Other restrictions",
                     grepl("SEX", variable)~"Sex (male)",
                     TRUE ~ variable)) %>%
  mutate(variable = gsub("_|\\."," ", variable)) %>%
  mutate(variable = variable %>% as.character %>% as.factor() %>% forcats::fct_reorder(s, .desc = TRUE)) %>% 
     mutate(branch = branch %>% gsub("DMM partitions", "\nDMM partitions",.)) %>%
  ggplot() + 
    geom_tile(aes(y=y.level, x=variable, fill=estimate)) + 
    geom_point(aes(y=y.level, x=variable, alpha=p.value<p_cut_off)) + 
    #scale_fill_gradient2("Log OR", low=scales::muted("blue"), mid="white", high=scales::muted("red")) +
    scale_fill_gradient("Log OR", low="white", high="#1a1a1a") +  
    facet_grid(variable_group~substates, scale="free", space = "free") + 
    cowplot::theme_cowplot() +
    scale_alpha_discrete("FDR < 10%") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    ylab("Microbiome Bacteroides partitions") + xlab("AGP questionnaire data") + coord_flip() 

Avis : Using alpha for a discrete variable is not advised.

NA
NA
NA
NA
NA
NA

to do: plot shannon vs cumulated model weight and make an id card for partition 1 (ref)

fig4A + coord_flip()


(fig4A+fig4B+ plot_layout(widths = c(2,1))) / (fig4C) + plot_layout(heights = c(3,1))

  
  Z = model.matrix(object)
    probs <- object$fitted # avoid napredict from fitted.default
    coefs <- coef(object)
    if (is.vector(coefs)){ # ie there are only 2 response categories
        coefs <- t(as.matrix(coefs))
        probs <- cbind(1 - probs, probs)
    }
    coefdim <- dim(coefs)
    p <- coefdim[2L]
    k <- coefdim[1L]
    ncoefs <- k * p
    kpees <- rep(p, k)
    n <- dim(Z)[1L]
##  Now compute the observed (= expected, in this case) information,
##  e.g. as in T Amemiya "Advanced Econometrics" (1985) pp295-6.
##  Here i and j are as in Amemiya, and x, xbar are vectors
##  specific to (i,j) and to i respectively.
    info <- matrix(0, ncoefs, ncoefs)
    Names <- dimnames(coefs)
    if (is.null(Names[[1L]])) {Names <- Names[[2L]] } else {
      Names <- as.vector(outer(Names[[2L]], Names[[1L]],
                             function(name2, name1)
                                 paste(name1, name2, sep = ":")))}
    dimnames(info) <- list(Names, Names)
    x0 <- matrix(0, p, k+1L)
    row.totals <- object$weights
    pb <- txtProgressBar(min=1,max=n)
    for (i in seq_len(n)) {
      
        Zi <- Z[i, ]
        xbar <- Zi * rep(probs[i, -1, drop = FALSE], kpees)
        for (j in seq_len(k+1)){
            x <- x0
            x[, j] <- Zi
            x <- x[, -1, drop = FALSE]
            x <- x - xbar
            dim(x) <- c(1, ncoefs)
            info <- info + (row.totals[i] * probs[i, j] * crossprod(x))
        }
        setTxtProgressBar(pb, i)
    }
    info

hysteresis exploration

metadata_enterotypes_num %>%
  mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  filter(branch=="Bacteroides DMM types" | Enterotypes_id =="m5") %>%
  select(SALTED_SNACKS_FREQUENCY,VEGETABLE_FREQUENCY,SUGARY_SWEETS_FREQUENCY,TYPES_OF_PLANTS,FRUIT_FREQUENCY,ANTIBIOTIC_HISTORY) %>%
  mutate_all(as.numeric) %>%
  ade4::dudi.pca(scannf=FALSE, nf=2) -> metadata_enterotypes_pca_branch
  


metadata_enterotypes_num_shannon = merge(metadata_enterotypes_num,shannon, by="row.names")
  

metadata_enterotypes_num_shannon %>%
  mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  filter(branch=="Bacteroides DMM types" | Enterotypes_id =="m5") %>%
  select(TYPES_OF_PLANTS, shannon, Enterotypes_id) %>%
  mutate(Enterotypes_id= forcats::fct_reorder(Enterotypes_id, shannon)) %>%
  #mutate(TYPES_OF_PLANTS=TYPES_OF_PLANTS%>% as.numeric()) %>%
  ggplot() + 
  geom_histogram(aes(x=TYPES_OF_PLANTS,fill=Enterotypes_id), stat="count", position="fill")


 merge(metadata_enterotypes_num_shannon,metadata_enterotypes_pca_branch$li, by.x="Row.names", by.y="row.names") %>% 
   mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  filter(branch=="Bacteroides DMM types" | Enterotypes_id =="m5") %>%
  mutate(Enterotypes_id= forcats::fct_reorder(Enterotypes_id, shannon)) %>%
   select(Axis1,shannon,Enterotypes_id) %>%
   ggplot() + 
  geom_point(aes(x=Axis1,y=shannon,col=Enterotypes_id), size=1,alpha=0.5)

merge(metadata_enterotypes_num_shannon,metadata_enterotypes_pca_branch$li, by.x="Row.names", by.y="row.names") %>% 
   mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  filter(branch=="Bacteroides DMM types" | Enterotypes_id =="m5") %>%
  mutate(Enterotypes_id= forcats::fct_reorder(Enterotypes_id, shannon)) %>%
   select(Axis1,shannon,Enterotypes_id) %>%
  ggplot() + 
  geom_histogram(aes(x=rank((-Axis1))+1,fill=Enterotypes_id), position="fill", bins=30) +
  scale_fill_viridis_d()  + scale_x_log10()


merge(metadata_enterotypes_num_shannon,metadata_enterotypes_pca_branch$li, by.x="Row.names", by.y="row.names") %>% 
   mutate(branch ="Clostridiales DMM types") %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m18","m14","m16","m17","m7"), "Prevotella DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m11","m15","m6","m8","m10","m3","m4","m2"), "Bacteroides DMM types",branch)) %>%
  mutate(branch = ifelse(Enterotypes_id %in% c("m13"), "Akkermansia DMM types",branch)) %>%
  filter(branch=="Bacteroides DMM types" | Enterotypes_id =="m5") %>%
  mutate(Enterotypes_id= forcats::fct_reorder(Enterotypes_id, shannon)) %>%
   select(Axis1,shannon,Enterotypes_id) %>%
  ggplot() + 
  geom_point(aes(x=rank((-Axis1))+1,y=shannon)) +
  scale_fill_viridis_d()