<jupyter_start><jupyter_text>Species extinction probability<jupyter_code># Housekeeping

library(car)
library(ggplot2)
library(InformationValue)
library(MASS)
library(plyr)
library(pscl)
# Read in data

df = read.table("../../../data/amplicon/extinction.txt",
                sep = "\t",
                header = T)

df = df[order(df$Strain, df$AB, df$imigration),]
head(df)
# Keep only strains that went extinct in at least one treatment

meta = as.data.frame(table(df$Strain, by = df$extinction))
never_extinct_strains = meta$Var1[meta$by == 1 & meta$Freq == 0]
df = df[!(df$Strain %in% never_extinct_strains),]
head(df)<jupyter_output><empty_output><jupyter_text>Model specification<jupyter_code># Binary response --> logistic regression (= binomial glm)

# Let's start from model including all experimental treatments

M1 = glm(extinction ~ factor(imigration) + 
           AB +
           Strain +
         MIC_sm,
         family = binomial, data = df)

# Let's use AIC to select best model among all interactions

M2 = stepAIC(M1, 
             scope=list(upper = ~factor(imigration)*
                            AB*
                            Strain*MIC_sm,
                            lower = ~1),
            trace = FALSE)

anova(M2, test = "Chi")

# Antibiotic, immigration and species all affect extinction probability
# The effect of antibiotic depends on the species
pR2(M2)

# Maximum likelihood pseudo r-squared (r2ML) indicates that model
# explains 75 % of variation in extinction probability (very good)
# Sketch the model fit

MyData = ddply(df, 
                .(imigration, Strain, AB),
                summarize,
               AB = c(0,4,16,128))
ilink = family(M2)$linkinv
MyData = cbind(MyData, predict(M2, MyData, type = "link", se.fit = TRUE)[1:2])
MyData = transform(MyData, Fitted = ilink(fit), Upper = ilink(fit + (2 * se.fit)),
                Lower = ilink(fit - (2 * se.fit)))

colnames(MyData)[1] = "Immigration"
MyData$Immigration = factor(MyData$Immigration, 
                            levels = c("0", "1"),
                            labels = c("italic('Immigration'*' absent')",
                                       "italic('Immigration'*' present')"))
MyData$AB = factor(MyData$AB, 
                            levels = c("0", "4", "16", "128"),
                            labels = c("0~mu*g*' mL'^-1",
                                       "4~mu*g*' mL'^-1",
                                       "16~mu*g*' mL'^-1",
                                       "128~mu*g*' mL'^-1"))

MyData$Strain = factor(MyData$Strain, 
                       levels = unique(MyData$Strain))

mylabels = c(expression(paste(italic("Bordetella avium"), " HAMBI 2160")),
           expression(paste(italic("Comamonas testosteroni"), " HAMBI 403")),
           expression(paste(italic("Cupriavidus necator"), " HAMBI 2164")),
           expression(paste(italic("Escherichia coli"), " JE2571(RP4)")),
           expression(paste(italic("Hafnia alvei"), " HAMBI 1279")),
           expression(paste(italic("Kluyvera intermedia"), " HAMBI 1299")),
           expression(paste(italic("M<NAME>"), " HAMBI 1292")),
           expression(paste(italic("Myroides odoratus"), " HAMBI 1923")),
           expression(paste(italic("Niabella yanshanensis"), " HAMBI 3031")),
           expression(paste(italic("Paracoccus denitrificans"), " HAMBI 2443")),
           expression(paste(italic("Pseudomonas chlororaphis"), " HAMBI 1977")),
           expression(paste(italic("Pseudomonas putida"), " HAMBI 6")),
           expression(paste(italic("Roseomonas gilardii"), " HAMBI 2470")),
           expression(paste(italic("Sphingobacterium multivorum"), " HAMBI 1874")),
           expression(paste(italic("Sphingobacterium spiritivorum"), " HAMBI 1896")),
           expression(paste(italic("Sphingobium yanoikuyae"), " HAMBI 1842")),
           expression(paste(italic("Stenotrophomonas maltophilia"), " HAMBI 2659")))
MyData2 = merge(MyData, df[, c("Strain", "MIC_sm")], by="Strain")
colnames(MyData2) = c('Strain', 'Immigration', 'AB', 'fit', 'se.fit', 'Fitted', 'Upper', 'Lower', 'MIC')
head(MyData2)
MyData2$MIC2 = ifelse(MyData2$MIC < 4, "<4", 
                      ifelse(MyData2$MIC >= 4 & MyData2$MIC < 16, "<16", 
                            ifelse(MyData2$MIC >= 16 & MyData2$MIC < 128, "<128", ">128")))
MyData2$MIC2 = factor(MyData2$MIC2, levels = c(">128","<128","<16", "<4"))
ggplot(MyData2, aes(x=Strain, y=Fitted, colour = MIC2)) + 
    geom_errorbar(aes(ymin=Lower, ymax=Upper), width=.1, colour = "grey") +
    geom_point() +
    facet_grid(AB~Immigration, labeller=label_parsed) +
    ylab("Extinction probability") +
    theme_classic() +
    scale_x_discrete(labels = parse(text=mylabels)) +
    scale_color_manual(values = c("black", "#522639","#8f4364","#cd6090")) + 
    labs(colour = "Intrinsic resistance level") +
    theme(axis.text.x  = element_text(angle = 90, vjust = 0.5, hjust = 1),
          strip.background = element_rect(colour = "white", fill = "white"),
          strip.text = element_text(face = "italic"),
          axis.title.x = element_blank())

ggsave("../../../manuscript/figures/extinction.pdf")
# ### Model validation
# 
# #Decide on optimal prediction probability cutoff for the model
# #The default cutoff prediction probability score is 0.5 or the ratio of 1's
# #and 0's in the training data. But sometimes, tuning the probability cutoff
# #can improve the accuracy in both the development and validation samples.
# #The InformationValue::optimalCutoff function provides ways to find the optimal
# #cutoff to improve the prediction of 1's, 0's, both 1's and 0's and to reduce
# #the misclassification error. Lets compute the optimal score that minimizes
# #the misclassification error for the above model.
predicted <- plogis(predict(M2, df))
optCutOff <- optimalCutoff(df$extinction, predicted)[1] 
# 

# #VIF - Variance inflation factor
# #Like in case of linear regression, we should check for multicollinearity
# #in the model. As seen below, all X variables in the model have VIF well below 4. --> OK
vif(M2)
# 
# #Misclassification Error
# #Misclassification error is the percentage mismatch of predcited vs actuals,
# #irrespective of 1's or 0's. The lower the misclassification error, the better is your model.
misClassError(df$extinction, predicted, threshold = optCutOff)
# # 0.0956 --> quite good
# 
# #ROC
# #Receiver Operating Characteristics Curve traces the percentage of true positives accurately
# #predicted by a given logit model as the prediction probability cutoff is lowered from 1 to 0.
# #For a good model, as the cutoff is lowered, it should mark more of actual 1's as positives
# # and lesser of actual 0's as 1's. So for a good model, the curve should rise steeply, indicating
# #that the TPR (Y-Axis) increases faster than the FPR (X-Axis) as the cutoff score decreases.
# #Greater the area under the ROC curve, better the predictive ability of the model.
plotROC(df$extinction, predicted)
# # area under curve 95.52 % --> very good
# 
# #Concordance
# #Ideally, the model-calculated-probability-scores of all actual Positive's, (aka Ones)
# #should be greater than the model-calculated-probability-scores of ALL the Negatives (aka Zeroes).
# #Such a model is said to be perfectly concordant and a highly reliable one. This phenomenon can
# #be measured by Concordance and Discordance.
# #In simpler words, of all combinations of 1-0 pairs (actuals), Concordance is the percentage of
# #pairs, whose scores of actual positive's are greater than the scores of actual negative's. For
# #a perfect model, this will be 100%. So, the higher the concordance, the better is the quality
# #of model.
Concordance(df$extinction, predicted)
# # Concordance > 95 % --> very good
# 
# #Specificity and Sensitivity
# #Sensitivity (or True Positive Rate) is the percentage of 1's (actuals) correctly predicted
# #by the model, while, specificity is the percentage of 0's (actuals) correctly predicted.
# #Specificity can also be calculated as 1???False Positive Rate.
sensitivity(df$extinction, predicted, threshold = optCutOff)
specificity(df$extinction, predicted, threshold = optCutOff)
# # predicts extinctions quite well but non-extinctions even better
# 
# #Confusion Matrix
# #The columns are actuals, while rows are predicteds.
confusionMatrix(df$extinction, predicted, threshold = optCutOff)
# 
# #Summary: good model that explains extinctions well
# Inspect MIC effect

M3 = glm(extinction ~ factor(imigration) * 
         AB *
         log(MIC_sm),
         family = binomial, data = df)

# Let's use AIC to select best model among all interactions

M4 = stepAIC(M3, scope=list(upper= ~factor(imigration)*
                            AB*
                            log(MIC_sm),
                            lower= ~1))

anova(M4, test = "Chi")

# Extinction risk increases with decreasing MIC value and with increasing antibiotic level
# Do same but integrate species and MIC information by grouping strains by MIC resistance
# to experimental streptomycin levels

df$MIC2 = ifelse(df$MIC_sm < 4, "<4", 
                      ifelse(df$MIC_sm >= 4 & df$MIC_sm < 16, "<16", 
                            ifelse(df$MIC_sm >= 16 & df$MIC_sm < 128, "<128", ">128")))
df$MIC2 = factor(df$MIC2, levels = c(">128","<128","<16", "<4"))

M5 = glm(extinction ~ factor(imigration) * 
         AB *
         MIC2,
         family = binomial, data = df)

# Let's use AIC to select best model among all interactions

M6 = stepAIC(M5, scope=list(upper= ~factor(imigration)*
                            AB*
                            MIC2,
                            lower= ~1))

anova(M6, test = "Chi")

# Same result
# Predict extinction risk at increasing MIC values

mean_probs = aggregate(Fitted ~ Immigration-AB-MIC2, 
                      data = MyData2, FUN = function(x){mean(x)})

head(mean_probs)
mean_probs$MIC2 = factor(mean_probs$MIC2, levels = c("<4", "<16", "<128", ">128"),
                        labels = c(expression("MIC"<"4"), 
                                   expression("MIC = 4"-"16"), 
                                   expression("MIC = 16"-"128"), 
                                   expression("MIC">"128")))
mean_probs$AB = as.factor(as.character(gsub("~.*", "", mean_probs$AB)))
mean_probs$AB = factor(mean_probs$AB, levels = c("0", "4", "16", "128"))
p2 = ggplot(mean_probs, aes(AB, Fitted, colour = AB)) +
    stat_summary(fun.y = "mean", fun.ymin = "mean", 
                 fun.ymax= "mean", 
                 size=0.5,
                 width = 0.6,
                 geom = "crossbar") +
    scale_color_manual(values = c("black", "#522639","#8f4364","#cd6090")) + 
    facet_grid(MIC2 ~ Immigration, labeller=label_parsed) +
    ylab("Mean extinction probability across species") +
    xlab("Streptomycin level") +
    theme_classic() +
    labs(colour = "Intrinsic resistance level") +
    theme(strip.background = element_rect(colour = "white", fill = "white"),
          strip.text = element_text(face = "italic"),
         legend.position="none")

suppressWarnings(print(p2))

ggsave("../../../manuscript/figures/extinction_means.pdf")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12989451099588772, 'fraction_numerical': 0.027802610405864472, 'mean_word_length': 3.013993541442411, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 26, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '118675', 'license': 'MIT', 'n_tokens_mistral': 3814, 'n_tokens_neox': 3383, 'n_words': 983, 'path': 'src/notebooks/amplicon_data/.ipynb_checkpoints/extinction-checkpoint.ipynb', 'repo_name': 'johannescairns/repeatable_dynamics', 'size': '11198'}