Refresh tutorials with robust variance diagnostics

Author: jaydu1

docs/source/main_function/lfc.rst

@@ -9,5 +9,73 @@ with effect estimates, standard errors, and BH-adjusted p-values.
 For screens with hundreds of perturbations use :func:`gcate_lfc_batch`, which
 runs GCATE and LFC in batches to keep peak memory bounded.
 
+Choosing the variance estimator
+-------------------------------
+
+``LFC`` defaults to ``usevar='unequal'`` (Welch inference), which estimates
+the treatment and control variances separately. Retain this default for
+perturbation screens and for case-control, bulk, and donor-level pseudo-bulk
+analyses. Independence of the rows does not imply equal treatment-arm
+variances: disease severity, biological response, residual composition,
+library size, treatment imbalance, and heterogeneous expression can all make
+pooled inference anti-conservative.
+
+Treat ``usevar='pooled'`` as an opt-in sensitivity or legacy analysis. Use it
+only when equal arm variances have a strong scientific and empirical
+justification; do not select it merely because a case-control study has a
+small number of independent samples. Pooled inference can produce much smaller
+standard errors and substantially more discoveries. For a deliberately
+justified batched sensitivity analysis, pass
+``lfc_kwargs=dict(usevar='pooled')``.
+
+Donor-level independence alone does not establish equal arm variances. For
+example, the SEA-AD tutorial uses ``usevar='unequal'`` because disease severity,
+inter-individual response, residual cell composition, and library-size
+variation can produce different gene-wise variability between disease groups.
+
+Welch inference does not itself model within-subject correlation. Repeated
+cells from the same donor or experimental unit should still be pseudo-bulked
+or handled with cluster-aware inference; ``usevar='unequal'`` only protects
+against unequal arm variances.
+
+Propensity diagnostics
+----------------------
+
+:func:`estimate_propensity_scores` estimates propensity scores without fitting
+the outcome model.  Use ``K=5`` to obtain out-of-fold scores for positivity and
+overfitting diagnostics.  :func:`summarize_propensity_scores` reports overlap,
+tail mass, and inverse-weight effective sample size, while
+:func:`plot_propensity_scores` compares treatment and control distributions.
+
+Both the standalone estimator and ``LFC`` use class-balanced logistic
+propensity fitting by default. This preserves historical causarray behavior and
+ensures that standalone overlap diagnostics describe the same nuisance model
+used for effect estimation. Pass ``class_weight=None`` to
+``estimate_propensity_scores`` or ``ps_class_weight=None`` to ``LFC`` for a
+calibrated-probability sensitivity analysis.
+
+``LFC`` uses the standard AIPW pseudo-outcome, which may be negative for
+individual cells even though its counterfactual mean is positive. Individual
+pseudo-outcomes are never clipped because doing so can bias the arm means,
+particularly when a large shared control group is compared with much smaller
+treatment groups.
+
+For a calibrated-propensity sensitivity analysis, use
+``LFC(..., ps_class_weight=None)`` and diagnose the matching scores with
+``estimate_propensity_scores(..., class_weight=None)``.
+
+The result includes ``mean_control``, ``mean_treated``, and ``estimable``. These
+are computed from the unclipped pseudo-outcomes. For numerical stability, valid
+aggregate arm means are floored at ``thres_diff`` only when constructing the
+log ratio and delta-method denominator. A pair with a nonfinite or nonpositive
+raw aggregate mean remains non-estimable, so the aggregate floor cannot create
+an extreme discovery from an invalid estimate.
+
 .. automodule:: causarray.DR_learner
    :members:
+
+.. automodule:: causarray.DR_estimation
+   :members: estimate_propensity_scores
+
+.. automodule:: causarray.diagnostics
+   :members:

docs/source/tutorial/case_control/sea_ad_case_control.ipynb

@@ -19,18 +19,19 @@ Here, we just use a subset of the data to demonstrate the workflow of the analys
 
 ```{r}
 library(Seurat)
-library(caret)
 
 sc.seurat <- readRDS("perturbseq-exneu.rds")
 
-Y <- data.frame(t(as.matrix(sc.seurat[['RNA']]$counts))) # cell-by-gene matrix
+Y <- data.frame(t(as.matrix(GetAssayData(sc.seurat, assay = "RNA", layer = "counts")))) # cell-by-gene matrix
 metadata <- sc.seurat@meta.data
 
 perturb <- metadata
 colnames(perturb) <- gsub("Perturbation", "trt_", colnames(perturb))
 perturb$trt_ <- relevel(as.factor(perturb$trt_), ref = "GFP")
-dmy <- dummyVars(" ~ trt_", data = perturb)
-A <- data.frame(predict(dmy, newdata = perturb))[,-1] # cell-by-trt matrix
+A <- data.frame(
+  model.matrix(~ trt_ - 1, data = perturb)[, -1, drop = FALSE],
+  check.names = FALSE
+) # cell-by-trt matrix; remove the first (GFP control) column
 ```
 
 For running causarray, we require the following inputs:
@@ -68,11 +69,12 @@ res_gate <- causarray$fit_gcate(Y, X, A, r, verbose=TRUE) # a list of results fr
 U <- res_gate[[2]]$U
 ```
 
-Next, we apply causarray to estimate the causal effects of perturbations on gene expression.
+Next, we apply causarray to estimate the causal effects of perturbations on gene expression. We use unequal (Welch) variance because treatment and control cells can have different pseudo-outcome variability; pooling those variances can be anti-conservative in a perturbation screen.
 
 ```{r}
 offsets <- log(res_gate[[2]][['kwargs_glm']][['size_factor']]) # use the precomputed size factors
-res <- causarray$LFC(Y, cbind(X, U), A, cbind(X_A, U), offset=offsets, verbose=TRUE)
+res <- causarray$LFC(Y, cbind(X, U), A, cbind(X_A, U), offset=offsets,
+                    usevar="unequal", verbose=TRUE)
 ```
 
 ```{r}
@@ -104,4 +106,4 @@ ggplot(discovery_counts, aes(x = Perturbation, y = Count)) +
   ggtitle('Number of Discoveries (padj < 0.1) for Each Perturbation Condition') +
   xlab('Perturbation Condition') +
   ylab('Number of Discoveries')
-```
+```

docs/source/tutorial/perturbseq/perturbseq-py.ipynb

@@ -29,25 +29,19 @@ library(Seurat)
     ## 
     ##     intersect
 
-``` r
-library(caret)
-```
-
-    ## Loading required package: ggplot2
-
-    ## Loading required package: lattice
-
 ``` r
 sc.seurat <- readRDS("perturbseq-exneu.rds")
 
-Y <- data.frame(t(as.matrix(sc.seurat[['RNA']]$counts))) # cell-by-gene matrix
+Y <- data.frame(t(as.matrix(GetAssayData(sc.seurat, assay = "RNA", layer = "counts")))) # cell-by-gene matrix
 metadata <- sc.seurat@meta.data
 
 perturb <- metadata
 colnames(perturb) <- gsub("Perturbation", "trt_", colnames(perturb))
 perturb$trt_ <- relevel(as.factor(perturb$trt_), ref = "GFP")
-dmy <- dummyVars(" ~ trt_", data = perturb)
-A <- data.frame(predict(dmy, newdata = perturb))[,-1] # cell-by-trt matrix
+A <- data.frame(
+  model.matrix(~ trt_ - 1, data = perturb)[, -1, drop = FALSE],
+  check.names = FALSE
+) # cell-by-trt matrix; remove the first (GFP control) column
 ```
 
 For running causarray, we require the following inputs:
@@ -75,7 +69,7 @@ causarray <- import("causarray")
 cat(causarray$`__version__`)
 ```
 
-    ## 0.0.6
+    ## 0.0.7
 
 ``` r
 # (Y, A) should be either data.frame or matrix
@@ -138,11 +132,14 @@ U <- res_gate[[2]]$U
 ```
 
 Next, we apply causarray to estimate the causal effects of perturbations
-on gene expression.
+on gene expression. We use unequal (Welch) variance because treatment and
+control cells can have different pseudo-outcome variability; pooling those
+variances can be anti-conservative in a perturbation screen.
 
 ``` r
 offsets <- log(res_gate[[2]][['kwargs_glm']][['size_factor']]) # use the precomputed size factors
-res <- causarray$LFC(Y, cbind(X, U), A, cbind(X_A, U), offset=offsets, verbose=TRUE)
+res <- causarray$LFC(Y, cbind(X, U), A, cbind(X_A, U), offset=offsets,
+                    usevar="unequal", verbose=TRUE)
 ```
 
     ## 'Estimating LFC...'