Incorrect behavior of minimization especially for low values of `p`

[ErdaradunGaztea]

The current implementation doesn't work as advertised. Let's dissect it:

sample(
  c(high_prob_arms, low_prob_arms), 1,
  prob = c(
    rep(
      p / length(high_prob_arms),
      length(high_prob_arms)
    ),
    rep(
      (1 - p) / length(low_prob_arms),
      length(low_prob_arms)
    )
  )
)

What happens if we pass p = 0 (i.e. complete randomness)? Turns out it'll always pick one of low_prob_arms.

What you probably intended on writing is -- give the same probability of picking for p = 0 (50%:50% for 2 arms case), always pick one of the best arms for p = 1 (100%:0% for 2 arms case), and scale linearly for all values inbetween.

How can we express it differently, so that we arrive at the correct solution?

Let's rephrase it. Let there be two components to weights: random and deterministic. If p = 0, we use 100% of the random component for all weights. If p = 1, we use 100% of the deterministic component for high prob weights only (and nothing for low prob). And to cover everything inbetween, we'll use an appropriate proportion of each of these components.

Summing up the previous paragraph, the values should look more-or-less like follows:

# high prob arm weight
hpa_weight_deterministic <- p * 1
hpa_weight_random <- (1 - p) * 1
hpa_weight <- hpa_weight_deterministic + hpa_weight_random

# low prob arm weight
lpa_weight_deterministic <- p * 0
lpa_weight_random <- (1 - p) * 1
lpa_weight <- lpa_weight_deterministic + lpa_weight_random

It would be silly to multiply by 1 or 0, so it'll simplify to the following:

hpa_weight <- 1
lpa_weight <- 1 - p

Now, the best part is that you can plug it into sample() without standardizing. sample() will do that internally, so you can write:

sample(
  c(high_prob_arms, low_prob_arms), 1,
  prob = c(
    rep(
      1,
      length(high_prob_arms)
    ),
    rep(
      1 - p,
      length(low_prob_arms)
    )
  )
)

and not care about it not summing up to 1.

[kamilsi]

Thanks again for raising this, @ErdaradunGaztea — your observation helped uncover a subtle but important edge case in the minimization randomization implementation.

After our internal discussion, we are leaning towards keeping the current logic as-is (i.e., no override or fallback for p = 0), since it does follow the mathematical definition of the minimization method. However, we also realized that the behavior for p ≤ 1/N can be surprising — and in many practical cases, leads to worse imbalance, even though it might look like it introduces "more randomness".

To address this, we're going to:

• ✍️ Improve the documentation to clearly describe what happens at different values of p
• ⚠️ Add a warning if p ≤ 1/N to inform users about the potential for unbalanced allocation
• ✅ Keep the default at p = 0.8, as it's a good trade-off in most common trial settings (2–4 arms), while clarifying how to adjust it for other scenarios

Here’s the proposed new documentation for the p parameter:

---

Parameter p

A numeric value in the interval (0, 1] that defines the probability of assigning a new subject to the treatment arm that minimizes imbalance.

• p = 1 – always assigns to the arm that minimizes imbalance (fully deterministic).
• p ∈ (1/N, 1) – assigns with a bias toward minimizing imbalance. This is the recommended range.
• p = 1/N – all arms have equal probability of selection (similar to randomize_simple()). The minimization mechanism has no effect.
• p < 1/N – the minimizing arm has a lower chance of being selected than in simple randomization. This may reduce balance or even lead to systematic imbalance.
• p = 0 – the minimizing arm is never selected. The patient is randomly assigned among the remaining arms. This typically leads to maximum imbalance.

Default value: p = 0.8
This value works well in most practical scenarios, especially for 2–4 treatment arms. It provides a good balance between minimizing imbalance and preserving unpredictability of allocation.

---

🧭 Recommendations

• Avoid values of p ≤ 1/N, unless you are explicitly testing edge cases or specific allocation behaviors.
• For completely random allocation, use the simple_randomization() function instead.
• Use higher p values when group balance is critical, particularly in small sample trials.
• For trials with more treatment arms, the required strength of p decreases:
  • For 2 arms: use p > 0.5
  • For 3 arms: use p > 0.33
  • For 5 arms: use p > 0.2
  • For 10 arms: use p > 0.1
• In larger trials with 5 or more arms, moderate values of p (e.g., 0.6–0.85) often provide a good compromise between balance and allocation unpredictability.

---

We’re really grateful you brought this up so clearly — your instinct about how the algorithm “should” behave reflected what many users might intuitively expect. Would you be interested in contributing a small PR to update the documentation? Or reviewing it once we have it pushed?

Let us know what you think — and thanks again for your thoughtful contribution!

[ErdaradunGaztea]

Certainly, I can write a PR or review it if you do it first. However, what you wrote is not exactly how your code behaves.

See, p is the probability of choosing one of the minimizing arms. For 2 arms, that's not an issue. For 3 and more, though, the issue arises when there's more than one minimizing arm. If we have 3 arms and we pick p = 0.33, then we have 2 minimizing arms -- what we get is not 0.33 chance of picking any random arm; it's P = 0.167 for those 2 minimizing arms and P = 0.66 for the third arm, the one that increases the imbalance.

If you wish to keep the current behaviour, please consider adding a parameter that allows the user to configure the interpretation of the p parameter, making it possible to switch between these two approaches.

[aleksandraduda2]

@ErdaradunGaztea I understand the situation you're referring to. For example, in the case of N = 3 arms and p = 0.33, it may happen that more than one arm minimizes the imbalance (because several of them have the same minimum value for the "total amount of imbalance"). This can happen, but doesn't have to. According to Sequential Treatment Assignment with Balancing for Prognostic Factors in the Controlled Clinical Trial, if such a situation occurs, the proposed solution is: "In the case of ties a random ordering can be determined." So, in that case, we would need to randomly select one of the "preferred" arms by assigning a random ranking, and for the remaining arms, we would apply the formula (1 - p) / (N - 1), which is currently implemented.

[kamilsi]

Thanks @aleksandraduda2 — your point is factually accurate: Pocock & Simon (1975) do indeed mention tie-breaking by randomly ordering tied arms before applying probabilities. So yes, selecting one minimizing arm at random and assigning the full p to it is a valid interpretation, and it appears in some variants of the method.

That said, the wording “may not behave as intended” is a bit misleading here, because:

• Our implementation explicitly follows the design from the Minirand package, which also splits p equally among all minimizing arms in the case of ties.
• This design choice is consistent with the Pocock & Simon framework, which does not prescribe a single tie-breaking strategy, but allows different ways of resolving ties.
• The equal-split approach ensures that the total probability mass p is fairly distributed across all minimizing arms. For example, with 3 arms and p = 0.8, and a tie between two arms:
  • Equal-split (current): 0.4, 0.4, 0.2 → total for optimal arms = 0.8
  • Random tie-break: 0.8, 0.1, 0.1 → total for optimal arms = 0.9
    This means the total probability assigned to optimal arms increases from 0.8 to 0.9 — which may be unintuitive, especially if the user expects p to represent the total weight given to minimizing arms.

We'll clarify this explicitly in the documentation, including a short explanation of the trade-offs between these approaches. If there's interest, we can also consider adding a tie_method = "equal" | "random" parameter to support both interpretations.

Would you be open to a short meeting next week to go over this together and decide whether we want to stick with one approach, offer both, or adjust the default? We'd be happy to sync up with anyone else from the team who’s interested.

Thanks again for prompting such a productive and thoughtful discussion!

[kamilsi]

Thanks all for the discussion today (@lwalejko, @ErdaradunGaztea, @salatak). Here's a summary of what we agreed and a brief technical review.

✅ Agreements

1. Parameter redesign
   We aim to introduce a more intuitive parameterization where:

   • p = 1 → fully deterministic assignment to minimizing arm(s)
   • p approaching 0 → allocation becomes more random, but not equivalent to simple randomization
   • Note: p = 0 is not allowed — users seeking full randomization should use a different endpoint (e.g., simple_randomization). In such cases, the system will reject the request and point users to the correct alternative.

   Two implementation options are under consideration:

   • v2: replace the current parameter with the new one; document and archive v1
   • extend v1: add a second parameter (e.g., p1) with the new semantics, supporting both options with documentation

2. Tie-breaking behavior
   We do not consider tie-breaking to be a problem.

   • When multiple arms equally minimize imbalance, the current implementation splits the probability mass p equally among them.
   • This effectively acts as a probabilistic tie-breaker and is consistent with approaches described by Pocock & Simon (1975), who allow random ordering in case of ties.
   • It also aligns with implementations like the minirand package.

3. API harmonization
   Based on feedback challenges with the REST API #99 (e.g., covariates vs. strata, unclear method, trailing slash issues, GET/POST inconsistencies), we will audit the API and align it before finalizing the parameter update strategy (v1 extension or v2).

If anything above looks off, please leave a comment or reaction so we can fix it!