Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(brglm2)    # Firth/bias-reduced logistic
library(broom)
library(dplyr)
df <- read.csv("KDM_data_25.csv") 
# Ensure factors with desired baselines
df <- df %>%
  mutate(
    Mutation = relevel(factor(Mutation),ref = "A3"),
    Genotype = relevel(factor(Genotype), ref = "HWT")
  )
mods <- setdiff(levels(df$Mutation), "A3")
one_mod_test <- function(mod_name) {
  sub <- df |> filter(Mutation %in% c("A3", mod_name)) |> droplevels()
  # Make sure LexA is baseline and HWT is baseline
  sub$Mutation     <- relevel(sub$Mutation, "A3")
  sub$Genotype <- relevel(sub$Genotype, "HWT")
  # Full model has the interaction = DID we care about
  m_full <- glm(cbind(Males, Females) ~ Genotype * Mutation,
                family = binomial(), data = sub, method = brglm2::brglmFit)
  # Reduced model = no interaction (no extra K16R drop under this modifier)
  m_red  <- update(m_full, . ~ . - Genotype:Mutation)
  # Likelihood-ratio test for that one interaction
  lrt <- anova(m_red, m_full, test = "LRT")
  p_lrt_two <- lrt$`Pr(>Chi)`[2]
  # Pull the interaction coefficient (effect size)
  rn <- grep("^GenotypeK16R:Mutation", rownames(coef(summary(m_full))), value = TRUE)
  co <- coef(summary(m_full))[rn, , drop = FALSE]
  beta <- co[1, "Estimate"]         # log ratio-of-odds-ratios (log ROR = DID)
  se   <- co[1, "Std. Error"]
  ROR  <- exp(beta)
  # One-sided p (directional: expecting beta < 0)
  p_one <- if (beta < 0) p_lrt_two/2 else 1 - p_lrt_two/2
  tibble(
    Mutation = mod_name,
    log_ROR = beta, SE = se, ROR = ROR,
    p_LRT_two_sided = p_lrt_two,
    p_LRT_one_sided_DIR_is_ROR_lt_1 = p_one
  )
}
res <- bind_rows(lapply(mods, one_mod_test)) |>
  mutate(p_FDR_one_sided = p.adjust(p_LRT_one_sided_DIR_is_ROR_lt_1, method = "BH")) |>
  arrange(p_FDR_one_sided)
library(dplyr)
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(brglm2)
# --- Load data (edit path if needed) ---
df <- read.csv("KDM_data_25.csv", stringsAsFactors = FALSE)
# --- Ensure factors & baselines (HWT, A3) ---
df <- df %>%
  mutate(
    Mutation     = relevel(factor(Mutation),     ref = "A3"),
    Genotype = relevel(factor(Genotype), ref = "HWT")
  )
# --- Firth (bias-reduced) logistic on the full dataset ---
model <- glm(
  cbind(Males, Females) ~ Genotype * Mutation,
  family = binomial(),
  data   = df,
  method = brglm2::brglmFit
)
# --- EMMs on logit scale, simple effect (K16R - HWT) within each RNAi ---
emm  <- emmeans(model, ~ Genotype | Mutation)          # stays on link (logit) scale
simp <- contrast(emm, "revpairwise")               # K16R - HWT because HWT is the ref
# --- DID vs LexA: (simple_effect at modifier) - (simple_effect at LexA) ---
levs <- levels(simp)$Mutation
stopifnot("A3" %in% levs)
did <- contrast(
  simp,
  method = "trt.vs.ctrl",
  ref    = which(levs == "A3"),
  by     = "contrast"
)
# --- Summarize (no multiple-testing adjustment here; we do our own) ---
s <- summary(did, infer = TRUE, adjust = "none")
# Parse modifier name from “X - LexA” and clean any parentheses
mod_name <- gsub("[()]", "", sub(" - A3$", "", s$contrast))
# One-sided p for the direction you expect (K16R < HWT => log_ROR < 0)
z <- s$estimate / s$SE
p_one <- pnorm(z, lower.tail = TRUE)  # one-sided "less"
# BH–FDR across modifiers
p_fdr <- p.adjust(p_one, method = "BH")
# Add interpretable effect sizes (ratio of odds ratios) and 95% CI
res <- data.frame(
  Mutation        = mod_name,
  log_ROR     = s$estimate,
  SE          = s$SE,
  z           = z,
  p_one_sided = p_one,
  p_FDR_one_sided = p_fdr,
  ROR         = exp(s$estimate),
  ROR_LCL95   = exp(s$estimate - 1.96 * s$SE),
  ROR_UCL95   = exp(s$estimate + 1.96 * s$SE)
) %>%
  arrange(p_FDR_one_sided)
res

```