Lyft Data Scientist (Entry Level) - Comprehensive Interview Preparation Guide

Direct answer: Feature selection keeps a subset of the ORIGINAL features (so each surviving feature is still directly interpretable), while dimensionality reduction (like principal component analysis, PCA) projects all features into a new, smaller set of derived axes that are combinations of the originals; the choice mostly comes down to whether you need to explain individual features to a stakeholder or regulator, or whether you only care about downstream predictive performance.

Structured elaboration:
Feature selection preserves interpretability by construction (a surviving feature is still "age" or "income," not some blend of the two), which matters when stakeholders, auditors, or regulators need to understand what's driving a prediction. It's also compatible with any downstream model without extra transformation overhead at serving time.

Dimensionality reduction can capture correlated structure across many original features into fewer derived dimensions, often preserving more of the total predictive signal per retained dimension than selection would, at the direct cost of interpretability (a principal component is a weighted combination of many original features, not a thing a business stakeholder can point to) and at the operational cost of needing to persist and correctly re-apply the exact same transformation (the same fitted components) at both training and serving time, forever, as a new pipeline dependency.

Worked example: A credit-risk model that must explain to a regulator why an application was denied needs feature selection, or at minimum an interpretable model on selected features, since "your application was influenced by principal component 3" is not an acceptable explanation to a regulator or a rejected applicant. A large-scale internal ranking model with no such explainability requirement, and many genuinely correlated raw signals, might get more predictive lift per dimension from PCA, since it doesn't need to explain any individual dimension to anyone outside the modeling team.

Trade-offs and pitfalls: A subtler cost of dimensionality reduction that's easy to miss operationally: the fitted transformation (the principal components themselves) has to be versioned and persisted exactly like a model artifact, and applied identically at training and serving time; if the components are ever refit on a different data slice than the one that trained the downstream model, the model's inputs silently no longer mean what it was trained on, which is the same training-serving-consistency problem that applies to any other feature transformation.

Approach

list(orig) is a shallow copy: it allocates a new outer list, but every element it holds is the exact same dict object as in orig, so mutating something inside one of those shared dicts (like the list stored under 'a') is visible through both orig and sh. copy.deepcopy(orig) recursively duplicates the dicts and the lists inside them, so dp shares nothing with orig at any level and is unaffected by mutating orig's nested state.

Code (Python 3.12)

import copy

orig = [{'a': [1, 2]}, {'b': [3]}]
sh = list(orig)
dp = copy.deepcopy(orig)

orig[0]['a'].append(99)   # mutate the nested list two levels deep

print("orig:", orig)
print("sh:  ", sh)
print("dp:  ", dp)
print("sh[0] is orig[0]:", sh[0] is orig[0])
print("dp[0] is orig[0]:", dp[0] is orig[0])

Output:

orig: [{'a': [1, 2, 99]}, {'b': [3]}]
sh:   [{'a': [1, 2, 99]}, {'b': [3]}]
dp:   [{'a': [1, 2]}, {'b': [3]}]
sh[0] is orig[0]: True
dp[0] is orig[0]: False

Trace through both

sh = list(orig): a new list object is created for sh, but sh[0] is literally orig[0], the same dict, and sh[1] is literally orig[1]. orig[0]['a'] is a list object nested inside that shared dict, so sh[0]['a'] is that same list too, three references deep (orig[0]['a'], sh[0]['a'], and the object itself are all one thing). Calling .append(99) mutates that one shared list in place, so the change is visible through orig, sh, or any other name bound to the same dict, because nothing about list(orig) ever touched anything below the outer list.

dp = copy.deepcopy(orig): deepcopy walks into orig, sees it is a list, creates a new list, then for each element (a dict) creates a new dict, then for each value inside that dict (a list) creates a new list too, copying [1, 2]'s contents into a brand-new list object. By the time orig[0]['a'].append(99) runs, dp[0]['a'] is an entirely separate list from orig[0]['a'], so the mutation has nowhere to propagate to.

Key points

• The depth of sharing in a shallow copy is exactly one level: the container is new, its immediate contents are shared. list(orig), orig.copy(), orig[:], and copy.copy(orig) all produce the same one-level-shallow result for a plain list.
• deepcopy uses an internal memo dictionary keyed by id() to avoid copying the same shared object twice or infinitely recursing on cyclic structures; that memoization is also why two elements of orig that happen to reference the same nested object stay correctly linked to each other in the deep copy, rather than becoming two independent copies that have drifted apart.

Complexity and edge cases

list(orig) is $O(k)$ where $k$ is the number of top-level elements (just copying references). copy.deepcopy(orig) is $O(n)$ where $n$ is the total number of objects in the entire nested structure, since every one of them gets visited and duplicated.

• If orig[0] and orig[1] both referenced the same inner dict object (not the case in this example, but a common real bug), a shallow copy would obviously still share it; a deep copy would still correctly produce two shallow-copy targets that both point to one newly-copied dict, preserving the original aliasing relationship rather than accidentally splitting a shared object into two independent ones.
• Tuples, sets, and custom objects nested inside the structure are handled by deepcopy too (it dispatches on type, and respects a class's own __deepcopy__ if one is defined); list(orig)-style shallow copying has no equivalent generalization beyond "copy this one container".

Direct answer

I treat two weeks as a research spike with a decision at the end, not open-ended learning time. I spend it testing the vendor's own specific claims against a real slice of our workload in an isolated trial that can't touch production, and I decide in advance what result would make me say no, so the verdict isn't a last-minute gut call.

Structured elaboration

• Find the two or three claims that actually gate the decision. I don't try to become an expert in the whole component. I identify the handful of things that, if false, would kill the decision (does it handle our real data volume, is it compatible with what we already depend on, does its failure behavior make sense), and I aim the whole two weeks at testing those.
• Test the claims myself instead of trusting the documentation. Vendor docs and marketing describe the happy path. I build the smallest thing that proves or disproves the specific claim using our own representative data or traffic shape, not the vendor's demo dataset.
• Keep the trial isolated with a clear way back out. The evaluation runs in a sandbox or a feature-flagged path (gated behind a feature flag, a toggle that turns a new component on for only a slice of traffic, without needing a separate deploy to turn it back off) that can't reach real customer data, and I know before I start how quickly we could rip it back out if it doesn't work, so trying it never becomes a one-way door.
• Decide the "say no" triggers before I see the results, not after. Examples: it fails under our expected traffic at even a modest multiple, there's no realistic exit path if we need to remove it later, or its security posture doesn't clear a bar we've already set. Deciding this in advance keeps the deadline from quietly lowering the bar.
• Under a genuinely compressed timeline this same shape compresses further. If instead of two weeks I had days, for instance needing to understand and counter an unfamiliar type of threat quickly, I'd skip the exploratory tour entirely and go straight at the one or two claims that actually gate whether we're safe, using whatever cheap check answers that fastest.
• Write the finding down either way. A short adoption note (what I tested, what passed, what didn't, the verdict) means the next person evaluating something similar doesn't redo this from scratch.
• If we adopt it, the first real use is staged, not a big rollout. A small, reversible slice of production traffic with its own explicit checks, expanded only once that holds up.

Worked example

A team I was on had two weeks to decide whether to adopt a third-party message-queuing service for a path that mattered a lot, with thin documentation and nobody on the team who'd used it in production. Instead of reading everything, I picked out the two claims that actually mattered to us: that it could sustain our peak message rate, and that we could get our data back out cleanly if we ever needed to leave. I spent the first three days building a minimal proof of concept against a sandbox account, fed it a replay of a real day's traffic rather than a toy example, and it held up. I spent a day specifically testing the export path, since a dead end there was one of my pre-agreed reasons to say no, and it worked cleanly. With about five days left I wrote up a one-page recommendation with what I'd tested, what I hadn't had time to test, and the specific evidence behind each claim, and we adopted it behind a feature flag on a low-traffic queue first, with its own success checks, before moving anything critical onto it.

Trade-offs and pitfalls

The biggest trap is spending the whole window reading and exploring instead of testing the load-bearing claims, which leaves you with broad but shallow familiarity and no real evidence at decision time. The opposite trap, trusting the vendor's claims at face value because the deadline is tight, is worse: it just moves the real evaluation to production, after you've already committed. Testing directly against live systems instead of an isolated trial is the other classic mistake, since it turns an evaluation into an incident risk. And skipping the write-up because the deadline already felt tight just guarantees the next evaluator repeats your work.

Direct answer

The underlying finding stays identical, but you change altitude, vocabulary, and level of supporting detail. An executive gets the headline, the business impact, and the recommended decision in one or two lines up front. The implementation team gets the mechanism, the caveats, and enough of the underlying data to act on it correctly.

Structured elaboration

• Altitude: conclusion-first for the executive, versus enough method detail for the team to trust and reproduce the diagnosis.
• Vocabulary: business-impact language (revenue, risk, timeline) for the executive, technical specifics (segments, funnels, thresholds) for the team.
• Format: a one-slide or one-paragraph summary versus a working document with a data appendix.
• What must never change: the number itself and the direction of the conclusion, in both versions.

Worked example

Finding: onboarding drop-off at step 3 is costing an estimated 6% of new signups per month.
Executive version: "we're losing about 6 of every 100 new signups at the step-3 confirmation screen, fixing it could recover meaningful revenue this quarter, recommend prioritizing it."
Team version: "62% of that drop-off happens on mobile between form submit and confirmation render, median time to abandon is 9 seconds, this looks like a loading-state issue on mobile specifically."
Both versions agree on the 6% headline number and the recommendation to prioritize the fix.

Trade-offs and pitfalls

The two versions can quietly drift into different conclusions if you're not careful, always trace both back to the same underlying analysis. Over-simplifying for the executive can also strip out the one caveat that would have changed their decision, so pick what to omit deliberately, not by default.

What the interviewer probes next

Expect a question about what happens when the executive summary gets forwarded on without you in the room, and how you prevent it from being read out of context.

Direct answer

For a single numeric column, start with a histogram (or a KDE, a kernel density estimate) for shape, a boxplot for a compact five-number summary and outlier flags, and the standard descriptive statistics (mean, median, standard deviation, min, max). For a single categorical column, start with a bar chart of value counts and the distinct-value count itself, since cardinality changes what "summarizing" even means.

What each reveals that the others don't

A histogram shows the full shape: modality, skew, and gaps that a single summary number can't convey. A boxplot compresses that shape into quartiles and flags candidate outliers, which is more compact but hides multimodality (two distinct humps can produce a perfectly ordinary-looking boxplot). Descriptive statistics give you exact numbers to reference or compare across columns, but on their own they can be actively misleading on a skewed or multimodal distribution (the classic case: a dataset can have identical mean, variance, and correlation as another while looking completely different when plotted). For categorical data, a bar chart of value counts shows both the shape (is one category dominant, or is it fairly even) and immediately reveals cardinality problems (hundreds of near-unique values where you expected a handful of categories).

Worked example

For a price column: the histogram shows a long right tail with most values clustered under $50 and a scattering up past $500; the boxplot confirms a compact box around $20 to $45 with several points flagged above the upper fence; and the mean ($68) sitting noticeably above the median ($38) confirms the skew numerically. For a payment_method column: value counts show credit_card at 71%, paypal at 24%, and four other methods splitting the remaining 5%, with a distinct count of 6, telling you immediately this is a genuinely low-cardinality categorical column you can treat straightforwardly.

Trade-offs and pitfalls

Don't stop at descriptive statistics alone for a column you haven't looked at before: they're fast to compute but can hide exactly the shape information (multimodality, heavy tails, a spike of a single repeated value) that a chart would reveal in seconds. For categorical columns, always check the distinct count before choosing a visualization, since a bar chart of value counts is useless once cardinality climbs into the hundreds or thousands.

Direct answer

Keep alignment across time zones with three levers: shrink what actually needs real-time overlap by defaulting to async updates on a fixed template, protect a small deliberately scheduled overlap window for anything that truly needs live discussion, and make handoffs explicit in writing so context transfers cleanly across the boundary instead of depending on someone's memory.

Framework

Reduce dependence on overlap. Default to async status updates on a fixed cadence, and use written decision docs rather than requiring a live meeting for every decision. Most updates don't need a room, only genuinely ambiguous or high-stakes calls do.

Protect a deliberate overlap window. Negotiate a recurring block, even a short one, and rotate who takes the inconvenient time so the burden doesn't always fall on the same region.

Make handoffs explicit. When work crosses a time-zone boundary, produce a short written artifact rather than relying on a quick chat message. This matters most in ops-heavy, always-on contexts.

Worked example

Consider an on-call rotation providing 24/7 production coverage across three time zones (for example [Region A], [Region B], and [Region C]), where the two outer regions have little or no live overlap with each other.

• Shadow and overlap periods: the incoming region's on-call shadows the outgoing region's on-call for a short deliberate window at the shift boundary, even 15 to 30 minutes, to ask questions live before the outgoing engineer signs off.
• Written handoff template: a standard document filled at every handoff covering open incidents, any systems in a degraded state, changes deployed in the last shift, and explicit 'known risk' or 'do not touch' notes.
• Escalation expectations: a written policy defining what counts as page-worthy versus a handoff note, who the secondary on-call is in each region, and how long the incoming engineer has to acknowledge before it auto-escalates.
  Result: even with zero live overlap between two of the three regions, the written handoff plus the short shadow window from the middle region means each incoming on-call starts already briefed, instead of reconstructing state from raw logs.
  For non-ops roles the same mechanism applies with a different artifact, for example a design or product handoff might be a written decision log plus a recorded walkthrough rather than an incident handoff, but the principle (explicit written handoff over a live conversation) is the same.

Trade-offs and pitfalls

• Repeatedly scheduling occasional syncs at painful hours burns out whichever time zone draws the short straw. Rotate it deliberately.
• Async-only breaks down for genuinely ambiguous or high-stakes decisions. Some live channel for true emergencies still has to exist.
• A handoff template that's too heavy gets skipped under time pressure. Keep it short enough to fill in within a few minutes.
• Assuming a chat message counts as a handoff is the actual failure mode this whole approach is designed to prevent. The structured artifact is the point, not the tool it's written in.

Direct answer

Split the logic into a fit step and a transform step: fit computes and stores a median per group (only for groups with enough non-missing observations) plus one overall median as the fallback, and transform looks up each row's group median, falls back to the overall median when the group is unseen, too small, or was entirely missing, and fills only the rows where the target is actually NaN (not-a-number, pandas' missing-value marker). Storing the fitted medians as plain attributes means the exact same transform logic can be reapplied to new data without recomputing anything from the original training set.

Structured elaboration

• Fit: group by the region column, compute the median of the target column per group (pandas' median() already ignores NaN within each group), and separately count non-null observations per group.
• Reliability threshold: a group's stored median is only trusted if it had at least min_group_size non-missing observations at fit time. Groups below that threshold, or groups that never appeared at fit time at all, fall back to the single overall median. This directly matches "too few observations, or was never seen at all" from the question, both routes land on the same fallback.
• Transform: map each row's group to its stored median (or NaN if the group wasn't kept), fill the gaps in that mapped Series with the global median, then use it only to fill rows where the target is actually missing, leaving already-present values untouched.
• Packaging this as a small class with fit/transform methods (the scikit-learn estimator convention) is what makes "fit once, reapply consistently" concrete: the fitted medians live on the object, not recomputed inline, so new data goes through transform alone and gets exactly the same group-to-median mapping the training data used.

Worked example

import pandas as pd

class GroupMedianImputer:
    def __init__(self, group_col, target_col, min_group_size=3):
        self.group_col = group_col
        self.target_col = target_col
        self.min_group_size = min_group_size

    def fit(self, X):
        non_null = X.dropna(subset=[self.target_col])
        counts = non_null.groupby(self.group_col)[self.target_col].size()
        medians = non_null.groupby(self.group_col)[self.target_col].median()
        # only trust a group's median if it had enough observed values
        reliable = medians[counts >= self.min_group_size]
        self.group_medians_ = reliable.to_dict()
        self.global_median_ = X[self.target_col].median()
        return self

    def transform(self, X):
        X_out = X.copy()
        fill_values = X_out[self.group_col].map(self.group_medians_).fillna(self.global_median_)
        mask = X_out[self.target_col].isna()
        X_out.loc[mask, self.target_col] = fill_values[mask].values
        return X_out

# fit on historical data with three regions: two well-observed, one tiny, one all-missing
train_df = pd.DataFrame({
    'region': ['east', 'east', 'east', 'west', 'west', 'west', 'tiny', 'north'],
    'income': [50000, 52000, 48000, 60000, 62000, 64000, 70000, float('nan')],
})
imputer = GroupMedianImputer('region', 'income', min_group_size=3)
imputer.fit(train_df)

# reapply the same fitted logic to new data, including an unseen region
incoming_df = pd.DataFrame({
    'region': ['east', 'west', 'tiny', 'south', 'north'],
    'income': [float('nan')] * 5,
})
new_df = imputer.transform(incoming_df)

Verified: imputer.group_medians_ keeps only east (50000.0) and west (62000.0), the two regions with at least min_group_size=3 non-missing observations; tiny (1 observation) and north (all-missing) are correctly excluded from group_medians_. Running transform on incoming_df fills east and west with their own group medians, and correctly falls back to the overall median (60000.0) for tiny, the never-seen south, and the all-missing north, confirming all three fallback triggers named in the question, too few observations, never seen, and entirely missing, resolve to the same value.

Trade-offs and pitfalls

• Complexity: fit is O(n) for the groupby plus counting pass; transform is O(n) for the map and fill. Memory is O(g) for the stored medians, where g is the number of distinct groups, independent of how large future batches are.
• Edge cases: a group with observations but they're all above/below normal ranges still gets a legitimate median, that's expected, not a bug; the target column being entirely NaN at fit time makes the global median itself NaN, so decide explicitly whether that should raise rather than silently leaving rows unfilled; non-numeric or badly-typed group keys (e.g. a mix of int and str region codes) will map inconsistently, so normalize the group key's dtype before fitting.
• A common mistake is computing the fallback as "global median of the fit set" but then never re-checking it against new data: if the new data's distribution shifts meaningfully from what was fit, the stored medians go stale silently, since nothing here re-fits automatically. This is a reasonable trade for reproducibility (the same input always imputes the same way), but it does mean this needs periodic re-fitting in production, the same way any other fitted preprocessing step does.
• min_group_size is itself a judgment call: too low and a two-observation "median" is really just picking one noisy value; too high and legitimately small-but-real segments always fall back to the (potentially very different) global median. There's no universally correct default; pick it based on how much you trust a small sample's median in this specific dataset.

Direct answer: Combining stability selection with permutation importance across cross-validation folds gives a feature-selection procedure that's robust to both a single unlucky data split and to a single noisy importance estimate, aggregating evidence across many folds rather than trusting any one run, as long as the tuning of that procedure itself is kept separate from the fold used to report its final performance.

Structured elaboration:
The algorithm: for each of several cross-validation folds, fit a model on the fold's training portion and compute permutation importance on its held-out portion; separately (or combined), run a stability-selection-style resampling within each fold to get a selection frequency per feature. Aggregate across all folds: a feature that's both consistently selected across resampling AND consistently shows a real permutation-importance drop across folds is a much stronger candidate than one that only clears the bar in a single fold or by a single method.

Pitfalls to guard against: correlated features can distort BOTH stability selection (as discussed elsewhere, Lasso arbitrarily picks one of a correlated group) and permutation importance (shuffling one of a correlated pair understates its importance since the model partially recovers signal from its correlated partner); aggregating across folds helps average out some of this instability but doesn't eliminate the underlying correlated-feature distortion, which needs to be checked for separately (a correlation-clustering step) if it's suspected.

Selection bias from reusing the same folds for tuning and evaluation is a separate, easy-to-miss pitfall: if the stability-selection threshold (what selection frequency counts as "stable") or the permutation-importance cutoff is chosen by looking at how well the resulting feature set performs on the same folds later reported as the final cross-validated score, that score is optimistically biased, since the threshold was implicitly tuned against the very data used to report performance. The fix is to nest the procedure: an outer set of folds is held out purely for final, unbiased evaluation, while an inner loop (further splits of only the remaining training data) is used to tune the stability threshold, the permutation-importance cutoff, and any other selection hyperparameter; only after those choices are fixed does the already-frozen procedure get scored on the outer fold's held-out data, which never influenced any tuning decision.

Worked example: A feature selected in 90% of stability-selection resamples across folds, AND showing a consistently large permutation-importance drop in every fold, is strong, well-corroborated evidence of genuine importance; a feature selected inconsistently across folds (40% in one, 85% in another) or with permutation-importance drops that vary wildly fold-to-fold is a signal of instability worth investigating before trusting it. Separately, if the 90%-selection-frequency threshold itself was chosen by scanning several candidate thresholds and picking whichever gave the best cross-validated accuracy on those same folds, that reported accuracy is inflated; nesting an inner tuning loop inside each outer fold (tune the threshold only on the inner splits, then apply the frozen threshold to the outer fold's held-out data) removes that inflation.

Trade-offs and pitfalls: This combined approach costs meaningfully more compute than either method alone (multiple folds times multiple resamples times multiple model fits), and nesting nested cross-validation to avoid the fold-reuse bias multiplies that cost further; for a very large candidate pool or an expensive-to-fit model, this cost needs to be weighed against a cheaper, single-pass alternative, accepting that the cheaper alternative's reported performance is likely somewhat optimistic if it skips the nested structure.

MECE principle: Mutually Exclusive, Collectively Exhaustive — partition problems into non-overlapping (no double-counting) categories that together cover all possibilities. As a data scientist, MECE helps build clean features, avoid leakage and ensure analyses capture full drivers of churn.

Issue tree (text form):

1. Product/Value
   1.1) Poor product fit (features don't meet needs)
   1.2) Low engagement/usage (inactive users)

2. Price/Packaging
   2.1) Price too high vs. perceived value
   2.2) Inflexible plans / bad trial-to-paid conversion

3. Customer Experience/Support
   3.1) Onboarding friction / setup issues
   3.2) Poor support response / unresolved tickets

4. External/Market
   4.1) Competitive displacement (better alternatives)
   4.2) Macro factors (economic downturn, industry changes)

Why MECE:

• Mutually exclusive: each branch targets a distinct causal domain (product vs. price vs. experience vs. external) so a single reason maps to one branch.
• Collectively exhaustive: together they cover internal (product, price, CX) and external reasons; sub-branches capture typical mechanisms (engagement, pricing sensitivity, support, competition, macro).
  How I'd operationalize: map churn labels to metrics (DAU, NPS, support SLA, plan price), run survival analysis and feature importance (Cox, random forest) to quantify contributions and prioritize interventions.

Approach (brief):

• Simulate many trials of sequential A/B tests with K correlated metrics. At each interim look compute test statistics for each metric (e.g., difference-in-means → z-scores using known/pooled SD). Generate correlated metric outcomes per user/experiment arm via multivariate normal with specified correlation matrix and per-metric effect sizes. Use an alpha-spending function (O'Brien–Fleming) to allocate cumulative alpha across L looks, convert to per-look boundaries (two-sided) then apply them to the simulated z-scores. Record first look crossing for any metric (family-wise decision). Aggregate empirical Type I error (when true effects = 0) and power (when effects > 0) across simulations.

Key choices & reasoning:

• Simulate raw outcomes per subject via multivariate normal to respect correlation structure; this properly captures joint distribution of z-statistics.
• Use O'Brien–Fleming spending: cumulative alpha_t = 2*(1 - Phi(z_alpha / sqrt(t_frac^{-1}))) — implement standard formula (or use numeric root to find boundary per look so cumulative spent equals spending function).
• At each look, test family-wise: reject null if any metric's |z| > boundary. Because correlation is simulated directly, FWER is captured empirically without ad-hoc correction.

Code skeleton (Python):

import numpy as np
from scipy.stats import norm, multivariate_normal
from math import sqrt

def obrien_fleming_alpha_spend(alpha, look_fractions):
    # returns two-sided critical z for each look so cumulative alpha follows OF spending
    # simple conservative approximation: boundary z_t = norm.ppf(1 - alpha/(2* (1.0/len(look_fractions))))
    # Better: use repeated numeric search so that cumulative spent at look t equals spending function.
    def spend(t):
        # OF cumulative spending at fraction t (Lan-DeMets OF approx)
        return 2 * (1 - norm.cdf(norm.ppf(1 - alpha/2) / sqrt(t)))
    z_bounds = []
    for t in look_fractions:
        cum_alpha = spend(t)
        z = norm.ppf(1 - cum_alpha/2)
        z_bounds.append(z)
    return np.array(z_bounds)

def simulate_one_trial(n_per_arm_by_look, K, mu_control, mu_treatment, cov_matrix):
    # simulate incremental additional users between looks
    # returns list of z-scores per look (K,) computed from cumulative data
    cumulative_N = 0
    sum_control = np.zeros(K)
    sum_treat = np.zeros(K)
    for add_n in n_per_arm_by_look:
        # simulate add_n new users per arm with multivariate normal metrics
        Xc = multivariate_normal.rvs(mean=mu_control, cov=cov_matrix, size=add_n)
        Xt = multivariate_normal.rvs(mean=mu_treatment, cov=cov_matrix, size=add_n)
        sum_control += Xc.sum(axis=0)
        sum_treat += Xt.sum(axis=0)
        cumulative_N += add_n
        # compute pooled mean and sd per metric
        mean_c = sum_control / (cumulative_N)
        mean_t = sum_treat / (cumulative_N)
        # sample variances (pooled) — using per-subject simulated variance
        # For speed in large sims, approximate se = sqrt(2*diag(cov)/(cumulative_N))
        se = np.sqrt(2 * np.diag(cov_matrix) / cumulative_N)
        z = (mean_t - mean_c) / se
        yield z  # z is array length K for this look

def run_simulation(n_sims, n_per_arm_by_look, K, effect_sizes, corr, alpha=0.05):
    # prepare parameters
    cov = np.outer(np.sqrt(effect_sizes_var := np.ones(K)), np.sqrt(effect_sizes_var)) * corr
    # here assume unit variances for metrics; effect_sizes are shifts in means for treatment
    mu_control = np.zeros(K)
    mu_treatment = np.array(effect_sizes)
    look_fracs = np.cumsum(n_per_arm_by_look) / sum(n_per_arm_by_look)
    z_bounds = obrien_fleming_alpha_spend(alpha, look_fracs)
    rejections = 0
    for sim in range(n_sims):
        for look_idx, z in enumerate(simulate_one_trial(n_per_arm_by_look, K, mu_control, mu_treatment, cov)):
            if np.any(np.abs(z) > z_bounds[look_idx]):
                rejections += 1
                break
    power = rejections / n_sims
    return power

Key points:

• Use multivariate_normal.rvs to simulate correlated metric vectors per subject; this preserves metric correlation.
• Compute z using cumulative sample sizes; se approximated from known covariances (for large-sample A/B tests this is fine). For small samples compute empirical variances from simulated data instead.
• O'Brien–Fleming spending implemented via Lan-DeMets approximation; for exact boundaries perform numeric solving to match cumulative alpha.

Complexity:

• Time: O(n_sims * total_sample_size * K) — dominated by random draws.
• Memory: O(K^2) for covariance; per-sim incremental memory small.

Edge cases and extensions:

• Binary or non-normal metrics: simulate from appropriate joint copula (e.g., Gaussian copula with marginals transformed).
• Unknown variances: estimate variances from data within each simulation.
• Multiplicity rules: current approach evaluates family-wise crossing of any metric; you can instead control per-metric alpha or use simultaneous test statistics with multivariate crossing thresholds derived analytically (more complex).
• Validate with null (all effects 0) to estimate Type I error, and with specified effects for power.

This skeleton gives a reproducible simulator to estimate empirical Type I error and power under sequential monitoring with correlated metrics and an OF alpha-spending rule.