Introduction

High-throughput time-lapse microscopy allows us to track thousands of cells across many wells and plates. However, these experiments generate complex data with multiple sources of variability:

1. Biological variability: Cells on different experimental plates (biological replicates) behave differently even under the same conditions.

2. Technical variability: Differences between wells on the same plate.

3. Batch effects: Systematic differences between plates (e.g., imaging days, reagent batches).

Standard statistical tests often ignore this hierarchical structure (cells \(\rightarrow\) wells \(\rightarrow\) plates), which can lead to false positives or missed discoveries.

The Bioconductor package cellmig addresses this by using Bayesian hierarchical models. It separates biological signals from technical noise and provides uncertainty-aware estimates (credible intervals) for treatment effects. This vignette guides you through installing cellmig, formatting your data, fitting a model, and interpreting the results.

Installation

To install this package, start R (version “4.5”) and enter:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("cellmig")

Data Structure

data("d", package = "cellmig")
str(d)

FALSE 'data.frame': 7560 obs. of  6 variables:
FALSE  $ well    : chr  "1" "1" "1" "1" ...
FALSE  $ plate   : chr  "1" "1" "1" "1" ...
FALSE  $ compound: chr  "C1" "C1" "C1" "C1" ...
FALSE  $ dose    : chr  "D1" "D1" "D1" "D1" ...
FALSE  $ v       : num  21.905 0.535 3.348 5.351 1.194 ...
FALSE  $ offset  : num  1 1 1 1 1 1 1 1 1 1 ...

head(d)

FALSE   well plate compound dose          v offset
FALSE 1    1     1       C1   D1 21.9054754      1
FALSE 2    1     1       C1   D1  0.5351645      1
FALSE 3    1     1       C1   D1  3.3484316      1
FALSE 4    1     1       C1   D1  5.3511913      1
FALSE 5    1     1       C1   D1  1.1942105      1
FALSE 6    1     1       C1   D1 26.7272783      1

Exploratory Data Analysis

Before modeling, it is good practice to visualize the raw data to check for obvious batch effects, outliers, multimodal distributions, or other imaging- related effects.

ggplot(data = d)+
  facet_wrap(facets = ~paste0("compound=", compound), 
             scales = "free_y", ncol = 2)+
  geom_sina(aes(x = as.factor(dose), col = plate, y = v, group = well), 
            size = 0.5)+
  theme_bw()+
  theme(legend.position = "top",
        strip.text.x = element_text(margin = margin(0.03,0,0.03,0, "cm")))+
  ylab(label = "migration velocity")+
  xlab(label = '')+
  scale_color_grey()+
  guides(color = guide_legend(override.aes = list(size = 3)))+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  scale_y_log10()+
  annotation_logticks(base = 10, sides = "l")

dm <- aggregate(v~well+plate+compound+dose, data = d, FUN = mean)
ggplot(data = dm)+
  facet_wrap(facets = ~paste0("compound=", compound), 
             scales = "free_y", ncol = 2)+
  geom_sina(aes(x = as.factor(dose), col = plate, y = v, group = well), 
            size = 1.5, alpha = 0.7)+
  theme_bw()+
  theme(legend.position = "top",
        strip.text.x = element_text(margin = margin(0.03,0,0.03,0, "cm")))+
  ylab(label = "migration velocity")+
  xlab(label = '')+
  scale_color_grey()+
  guides(color = guide_legend(override.aes = list(size = 3)))+
  guides(shape = guide_legend(override.aes = list(size = 3)))+
  scale_y_log10()+
  annotation_logticks(base = 10, sides = "l")

Model Fitting

Now we fit the hierarchical Bayesian model. The cellmig() function handles the complex Stan modeling behind the scenes.

o <- cellmig(x = d,
             control = list(mcmc_warmup = 300,  # Warmup iterations
                            mcmc_steps = 1000,  # Sampling iterations
                            mcmc_chains = 2,    # Number of chains
                            mcmc_cores = 2))    # Parallel cores

Interpreting Results

The model output o contains posterior distributions for all parameters. We focus on the overall treatment effects (\(\delta_t\)), which represent the change in velocity relative to the control (offset).

knitr::kable(o$posteriors$delta_t, digits = 2)

ggplot(data = o$posteriors$delta_t) +
  geom_line(aes(x = dose, y = mean, col = compound, group = compound)) +
  geom_point(aes(x = dose, y = mean, col = compound)) +
  geom_errorbar(aes(x = dose, y = mean, ymin = X2.5., ymax = X97.5., 
                    col = compound), width = 0.1) +
  ylab(label = expression("Log-Fold Change ("*delta*")")) +
  xlab("Dose") +
  theme(legend.position = "top")

ggplot(data = o$posteriors$delta_t) +
  geom_line(aes(x = dose, y = exp(mean), col = compound, group = compound)) +
  geom_point(aes(x = dose, y = exp(mean), col = compound)) +
  geom_errorbar(aes(x = dose, y = exp(mean), ymin = exp(X2.5.), 
                    ymax = exp(X97.5.), col = compound), width = 0.1) +
  ylab(label = expression("Fold Change ("*delta*"')")) +
  xlab("Dose") +
  theme(legend.position = "top")

Pairwise Comparisons

Often, you want to compare specific treatments against each other (e.g., Drug A vs. Drug B), not just against the control.

# Get pairwise comparisons (log-scale)
u <- get_pairs(x = o, exponentiate = FALSE)

# vislualize matrix of rhos
u$plot_rho

# visualize matrix of pis
u$plot_pi

# visualize volcano plot
u$plot_volcano

# View available group labels
groups <- get_groups(x = o)
head(groups)

FALSE   group_id group compound dose
FALSE 1        1 C2|D1       C2   D1
FALSE 2        2 C2|D2       C2   D2
FALSE 3        3 C2|D3       C2   D3
FALSE 4        4 C2|D4       C2   D4
FALSE 5        5 C2|D5       C2   D5
FALSE 6        6 C2|D6       C2   D6

# Compare all groups against Compound 2, Dose 1
u_violin <- get_violins(x = o, 
                        from_groups = groups$group,
                        to_group = "C2|D1",
                        exponentiate = FALSE)
u_violin$plot

Model Diagnostics and Calibration

Reliable uncertainty quantification requires verifying that the model fits the data well and that probability metrics (\(\pi\)) are calibrated to your experimental design. cellmig provides tools for both.

g <- get_ppc_violins(x = o, wrap = TRUE, ncol = 3)
g + scale_y_log10()

g <- get_ppc_means(x = o)
g

# run loo
loo_o <- rstan::loo(x = o$fit)

# Print diagnostic table
print(loo_o)

FALSE 
FALSE Computed from 1400 by 7560 log-likelihood matrix.
FALSE 
FALSE          Estimate    SE
FALSE elpd_loo -22300.4 106.7
FALSE p_loo        69.1   3.6
FALSE looic     44600.9 213.4
FALSE ------
FALSE MCSE of elpd_loo is 0.2.
FALSE MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 2.7]).
FALSE 
FALSE All Pareto k estimates are good (k < 0.68).
FALSE See help('pareto-k-diagnostic') for details.

# Plot diagnostic estimates for each cell
plot(loo_o)

# Which cells have k>0.7? Inspect these cells
which(loo_o$pointwise[,"influence_pareto_k"]>0.7)

FALSE integer(0)

rstan::stan_rhat(o$fit)

summary_stats <- summary(o$fit)$summary
hist(summary_stats[,"n_eff"])

rstan::check_hmc_diagnostics(object = o$fit)

FALSE 
FALSE Divergences:

FALSE 
FALSE Tree depth:

FALSE 
FALSE Energy:

Variance Components

You can inspect how much variability comes from different sources (plates, wells, treatments). This helps in planning future experiments.

# Plate-specific baseline effects
g_alpha_p <- ggplot(data = o$posteriors$alpha_p) +
  geom_errorbarh(aes(y = plate, x = mean, xmin = X2.5., xmax = X97.5.),
                 height = 0.2) +
  geom_point(aes(y = plate, x = mean)) +
  xlab("Plate Effect (log-scale)")

# Variance parameters (Biological vs Technical)
g_sigma <- ggplot() +
  geom_errorbarh(data = o$posteriors$sigma_bio,
                 aes(y = "Biological (Plate)",
                     x = mean, xmin = X2.5., xmax = X97.5.), height = 0.2) +
  geom_errorbarh(data = o$posteriors$sigma_tech,
                 aes(y = "Technical (Well)",
                     x = mean, xmin = X2.5., xmax = X97.5.), height = 0.2) +
  geom_errorbarh(data = o$posteriors$sigma_delta,
                 aes(y = "Treatment Variation",
                     x = mean, xmin = X2.5., xmax = X97.5.), height = 0.2) +
  geom_point(data = o$posteriors$sigma_bio,
             aes(y = "Biological (Plate)", x = mean)) +
  geom_point(data = o$posteriors$sigma_tech,
             aes(y = "Technical (Well)", x = mean)) +
  geom_point(data = o$posteriors$sigma_delta,
             aes(y = "Treatment Variation", x = mean)) +
  xlab("Standard Deviation")

g_alpha_p | g_sigma

Advanced: Dose-Response Profiles

For screens with multiple compounds and overlapping doses, cellmig can cluster compounds based on their dose-response similarity, accounting for uncertainty.

• Left: Dendrogram showing similarity between compounds.
• Middle/Right: Dose-response curves for overall (\(\delta\)) and plate-specific (\(\gamma\)) effects.

get_dose_response_profile(x = o, exponentiate = TRUE) +
  patchwork::plot_layout(widths = c(.7, 1, 2))

Session Info

sessionInfo()

FALSE R version 4.6.0 (2026-04-24)
FALSE Platform: x86_64-pc-linux-gnu
FALSE Running under: Ubuntu 24.04.4 LTS
FALSE 
FALSE Matrix products: default
FALSE BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
FALSE LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
FALSE 
FALSE locale:
FALSE  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
FALSE  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
FALSE  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
FALSE  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
FALSE  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
FALSE [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
FALSE 
FALSE time zone: Etc/UTC
FALSE tzcode source: system (glibc)
FALSE 
FALSE attached base packages:
FALSE [1] stats     graphics  grDevices utils     datasets  methods   base     
FALSE 
FALSE other attached packages:
FALSE [1] rstan_2.32.7        StanHeaders_2.32.10 ggforce_0.5.0      
FALSE [4] ggplot2_4.0.3       cellmig_1.3.4       BiocStyle_2.41.0   
FALSE 
FALSE loaded via a namespace (and not attached):
FALSE  [1] ggiraph_0.9.6           tidyselect_1.2.1        dplyr_1.2.1            
FALSE  [4] farver_2.1.2            loo_2.9.0               S7_0.2.2               
FALSE  [7] fastmap_1.2.0           lazyeval_0.2.3          tensorA_0.36.2.1       
FALSE [10] tweenr_2.0.3            fontquiver_0.2.1        digest_0.6.39          
FALSE [13] lifecycle_1.0.5         tidytree_0.4.7          posterior_1.7.0        
FALSE [16] magrittr_2.0.5          compiler_4.6.0          rlang_1.2.0            
FALSE [19] sass_0.4.10             tools_4.6.0             yaml_2.3.12            
FALSE [22] knitr_1.51              labeling_0.4.3          htmlwidgets_1.6.4      
FALSE [25] pkgbuild_1.4.8          curl_7.1.0              plyr_1.8.9             
FALSE [28] RColorBrewer_1.1-3      abind_1.4-8             aplot_0.2.9            
FALSE [31] withr_3.0.2             purrr_1.2.2             sys_3.4.3              
FALSE [34] grid_4.6.0              polyclip_1.10-7         stats4_4.6.0           
FALSE [37] gdtools_0.5.0           inline_0.3.21           scales_1.4.0           
FALSE [40] MASS_7.3-65             cli_3.6.6               rmarkdown_2.31         
FALSE [43] treeio_1.37.0           generics_0.1.4          RcppParallel_5.1.11-2  
FALSE [46] ggtree_4.3.0            reshape2_1.4.5          ape_5.8-1              
FALSE [49] cachem_1.1.0            stringr_1.6.0           parallel_4.6.0         
FALSE [52] ggplotify_0.1.3         BiocManager_1.30.27     matrixStats_1.5.0      
FALSE [55] vctrs_0.7.3             yulab.utils_0.2.4       V8_8.2.0               
FALSE [58] jsonlite_2.0.0          fontBitstreamVera_0.1.1 gridGraphics_0.5-1     
FALSE [61] patchwork_1.3.2         systemfonts_1.3.2       maketools_1.3.2        
FALSE [64] tidyr_1.3.2             jquerylib_0.1.4         glue_1.8.1             
FALSE [67] codetools_0.2-20        distributional_0.7.0    stringi_1.8.7          
FALSE [70] gtable_0.3.6            QuickJSR_1.10.0         tibble_3.3.1           
FALSE [73] pillar_1.11.1           rappdirs_0.3.4          htmltools_0.5.9        
FALSE [76] R6_2.6.1                evaluate_1.0.5          lattice_0.22-9         
FALSE [79] backports_1.5.1         ggfun_0.2.0             fontLiberation_0.1.0   
FALSE [82] bslib_0.11.0            rstantools_2.6.0        Rcpp_1.1.1-1.1         
FALSE [85] checkmate_2.3.4         gridExtra_2.3           nlme_3.1-169           
FALSE [88] xfun_0.57               fs_2.1.0                buildtools_1.0.0       
FALSE [91] pkgconfig_2.0.3