Morris Elementary Effects for Screening Mechanistic Parameters

Why Screen Parameters With Elementary Effects?

A mechanistic ODE model may contain dozens of parameters, while each model evaluation may require solving a stiff system, running a population simulation, or repeating calculations across many virtual subjects. Before investing in an expensive global sensitivity analysis, it is useful to identify parameters whose variation appears to have little effect on a specified model output. The Morris method provides a relatively economical way to do this by evaluating the output along short, one-parameter-at-a-time steps through the input space.

The result is a screening assessment, not a final attribution of output variation to each parameter. Morris effects can reveal overall influence and patterns consistent with nonlinearity or interactions, but they do not uniquely separate those causes. The output of interest, the parameter ranges or distributions, and the scale on which the output is measured all shape the result.

Definition: Let $g(\mathbf{u})$ be a scalar model output with inputs $\mathbf{u}=(u_1,\ldots,u_k)$ in the unit hypercube. A Morris elementary effect for input $i$ is the finite change in $g$ across a one-coordinate step, divided by that step’s signed length: $$ EE_i(\mathbf{u})=\frac{g(\mathbf{u}+\Delta u_i\mathbf{e}_i)-g(\mathbf{u})}{\Delta u_i}. $$ Here $\mathbf{e}_i$ is the unit vector in coordinate $i$, and $\Delta u_i$ is either $+\delta$ or $-\delta$. Thus the effect is expressed as output units per unit change in the scaled input.

Scaling Inputs and Building Trajectories

The elementary effect should be calculated on a clearly defined input scale. For bounded parameters, a common choice is to map each parameter range to $[0,1]$. If parameter $i$ ranges from $L_i$ to $U_i$, the affine mapping is $\theta_i=L_i+(U_i-L_i)u_i$. For parameters specified by independent distributions, one can instead use their marginal quantile functions, $\theta_i=F_i^{-1}(u_i)$, and evaluate the ODE model at $\boldsymbol{\theta}(\mathbf{u})$. In either case, the resulting elementary effect describes a change in the chosen unit-coordinate scale, not a change per unit of the original parameter.

A classical Morris design uses a grid with $p$ levels in each unit-coordinate dimension, where $p$ is an integer greater than one. A trajectory begins at a valid grid point and changes each of the $k$ coordinates exactly once by a step of magnitude $\delta=p/[2(p-1)]$. It therefore contains $k+1$ model evaluations. The starting point and the direction of each step are selected so every changed coordinate remains within the grid. Repeating the design with different trajectories produces a collection of effects for each parameter.

This construction is most straightforward for independent inputs. If the parameters are correlated, independently changing one unit coordinate at a time may create combinations that are implausible under the joint input distribution. A standard independent-factor Morris design should not be treated as a sensitivity analysis under such a joint distribution without a suitable dependence-aware design.

Worked Example: An Elementary Effect for ODE Clearance

Consider the hypothetical ODE $dC/dt=-kC$, with initial concentration $C(0)=1$. The target is the fraction remaining at $T=1$, so $g(u)=\exp[-k(u)]$. Suppose the clearance rate is bounded by $0.2$ to $1.0$ per time unit and is mapped linearly as $k(u)=0.2+0.8u$. Use $p=4$, which gives $\delta=4/[2(4-1)]=2/3$.

For a trajectory step from $u=0$ to $u=2/3$, the corresponding rates are $k(0)=0.2$ and $k(2/3)=0.7333$ per time unit. The ODE solution gives $g(u)=e^{-k(u)}$, hence $g(0)=e^{-0.2}=0.8187$ and $g(2/3)=e^{-0.7333}=0.4803$.

EE=\frac{0.4803-0.8187}{2/3} =\frac{-0.3384}{0.6667} \approx -0.5076

The estimated effect is about $-0.5076$ fraction units per unit change in the scaled input $u$. Its negative sign indicates that increasing clearance over this step decreases the fraction remaining. It is not a derivative per unit of clearance rate: the rate changes by $0.8$ per unit $u$, so the two scales must not be confused.

Summarizing Effects Across Trajectories

For each parameter, suppose $r$ trajectories produce effects $EE_i^{(1)},\ldots,EE_i^{(r)}$. Three summaries are commonly useful. The signed mean, $\mu_i$, indicates the average direction and size of the effects. The mean absolute effect, $\mu_i^*$, measures average magnitude without allowing positive and negative effects to cancel. The sample standard deviation, $\sigma_i$, describes how much the effects vary across the sampled locations and directions.

\mu_i=\frac{1}{r}\sum_{j=1}^{r}EE_i^{(j)},\qquad \mu_i^*=\frac{1}{r}\sum_{j=1}^{r}\left|EE_i^{(j)}\right|,\qquad \sigma_i=\sqrt{\frac{1}{r-1}\sum_{j=1}^{r}\left(EE_i^{(j)}-\mu_i\right)^2}.

A large $\mu_i^*$ signals that the output often changes substantially when parameter $i$ is perturbed on the selected scale. A small $\mu_i^*$ suggests the parameter may be a low-priority candidate for a more expensive analysis over the explored input space. A large $\sigma_i$ means the effect varies across the sampled steps. This variation can be consistent with nonlinearity, interactions with other inputs, or both; $\sigma_i$ does not distinguish those explanations. The signed mean alone can be misleading if effects change sign, which is why $\mu_i^*$ is also reported.

Worked Example: Reading Means and Spread Together

A hypothetical two-parameter ODE model predicts a cytokine concentration at 24 hours. A Morris design with $p=4$ gives $\delta=2/3$, and three trajectories are run. The table shows the concentration before and after each step. For simplicity, the listed steps all increase the relevant unit coordinate by $2/3$; concentration is in ng/L.

Parameter	Trajectory	Before (ng/L)	After (ng/L)	Elementary effect (ng/L per unit $u$)
Production rate	1	10.0000	11.3333	$(11.3333-10.0000)/(2/3)=2$
Production rate	2	12.0000	14.6667	$(14.6667-12.0000)/(2/3)=4$
Production rate	3	8.0000	10.0000	$(10.0000-8.0000)/(2/3)=3$
Clearance rate	1	10.0000	9.3333	$(9.3333-10.0000)/(2/3)=-1$
Clearance rate	2	12.0000	12.6667	$(12.6667-12.0000)/(2/3)=1$
Clearance rate	3	9.0000	9.0000	$(9.0000-9.0000)/(2/3)=0$

For production rate, the effects are $2,4,3$, so $\mu=(2+4+3)/3=3$ and $\mu^*=(2+4+3)/3=3$ ng/L per unit $u$. The sample standard deviation is $\sigma=\sqrt{[(-1)^2+1^2+0^2]/2}=1$ ng/L per unit $u$.

For clearance rate, the effects are $-1,1,0$. Thus $\mu=(-1+1+0)/3=0$, while $\mu^*=(1+1+0)/3=0.6667$ ng/L per unit $u$. Its sample standard deviation is $\sigma=\sqrt{[(-1)^2+1^2+0^2]/2}=1$ ng/L per unit $u$. The zero signed mean does not imply no effect: positive and negative changes cancel in $\mu$, but the nonzero $\mu^*$ reveals effects of both signs. In this hypothetical screening, production rate has the larger average magnitude, while clearance rate has sign-changing effects across the sampled trajectories.

Why Effect Spread Occurs: A Short Justification

The finite difference is an average slope over a step, not necessarily the derivative at one fixed point shared by all trajectories. If $g$ is differentiable along a step in coordinate $i$, the fundamental theorem of calculus gives a direct justification:

EE_i =\frac{g(\mathbf{u}+\Delta u_i\mathbf{e}_i)-g(\mathbf{u})}{\Delta u_i} =\frac{1}{\Delta u_i}\int_{0}^{\Delta u_i} \frac{\partial g}{\partial u_i}(\mathbf{u}+t\mathbf{e}_i)\,dt.

This holds for either sign of $\Delta u_i$, with the signed denominator matching the integration direction. Each effect averages the local slope along its own interval. If that slope is constant across the input space, effects will be similar apart from numerical error. If the slope changes with the parameter value or with other coordinates, effects can differ across trajectories. This explains why spread is useful as a warning that influence may be context-dependent, while also explaining why spread alone cannot identify whether nonlinearity or interactions are responsible.

Worked Example: The Input Distribution Changes the Screening Scale

Return to the decay model $dC/dt=-kC$, $C(0)=1$, with target $C(1)$, but now let $k$ range from $0.1$ to $10$ per time unit. Compare two mappings for the same trajectory step $u=0$ to $u=2/3$, using $p=4$ and $\delta=2/3$.

With a uniform rate mapping, $k(u)=0.1+9.9u$. The two rates are $0.1$ and $6.7$, giving outputs $e^{-0.1}=0.9048$ and $e^{-6.7}=0.001231$. Therefore,

EE_{\mathrm{uniform}} =\frac{0.001231-0.9048}{2/3} \approx -1.355

With a log-uniform rate distribution, the quantile mapping is $k(u)=0.1(100)^u$. At $u=2/3$, $k=2.154$, giving $C(1)=e^{-2.154}=0.1160$. The corresponding effect is

EE_{\mathrm{log\ uniform}} =\frac{0.1160-0.9048}{2/3} \approx -1.183

Both calculations use the same ODE and unit-coordinate step, but they explore different rate values. Their effects therefore describe sensitivity under different input mappings; they are not conflicting estimates of a scale-free property. The mapping should match the scientific input uncertainty or range that the screening is intended to represent.

Practical Design, Diagnostics, and Limitations

A practical screening analysis begins by naming one scalar target, such as a concentration at a specified time, an AUC over a specified interval, or a maximum response. It then specifies the input distributions or parameter ranges, maps the inputs to the unit cube, chooses $p$ and the number $r$ of trajectories, and calculates effects from each one-coordinate change. For an ODE model, each evaluation must use consistent initial conditions, solver tolerances, event handling, and output extraction. If several outputs are screened, repeat the calculation for each; their effects have different output units or require an explicitly chosen normalization.

Several checks are important before interpreting a ranking:

Check numerical resolution. Compare the output difference across a Morris step with the change caused by tightening ODE solver tolerances. If solver error is not small relative to the stepwise output change, the elementary effect may mostly reflect numerical noise.
Check the trajectory count. Repeating trajectories and comparing rankings can reveal whether a small design is giving unstable summaries. Morris effects are empirical summaries of sampled steps, not exact global quantities.
Check boundaries and discontinuities. Threshold rules, switching equations, events, and failed simulations can produce abrupt changes. Confirm that these are intended model behavior and that the trajectory stays within valid parameter bounds.
Check the chosen scale. A low $\mu_i^*$ only indicates a small effect over the chosen range or distribution and on the selected input and output scales. It does not establish that the parameter is intrinsically unimportant.

Local derivative methods are cheaper when sensitivity near one nominal parameter set is the only question, but they can miss effects elsewhere in the input space. One-at-a-time changes from a single baseline are easy to interpret, yet do not sample a range of contexts and can conceal changes in effect size or sign. More computationally demanding global methods can provide richer quantitative descriptions of output variation, but Morris screening is useful when the immediate goal is to reduce a candidate list before those analyses.

Morris screening has limits. It does not estimate a causal effect, prove that a parameter can be removed, or quantify a unique share of output variance. A parameter with a small effect for one target may matter for another target, at another time, or under another input distribution. Correlated parameters, discrete choices, and parameters constrained by mechanistic relationships require particular care; naive independent coordinate steps may not represent scientifically feasible combinations. Screening should guide follow-up, not replace model knowledge or validation.

Common Pitfall

Common mistake: “The signed mean elementary effect is zero, so the parameter has no influence.” This is wrong when effects of opposite signs cancel. In the clearance-rate example, effects of $-1$, $1$, and $0$ have $\mu=0$, but $\mu^*=0.6667$ ng/L per unit $u$. Examine $\mu^*$ for average magnitude and $\sigma$ for variation, alongside the signed mean. Neither a small nor a large value should be interpreted without the chosen target, input mapping, and numerical checks.

Define the target and input space.
Specify the model output, its units, and parameter ranges or distributions. Map inputs to a unit coordinate scale that matches the intended analysis.

Generate trajectories and solve the model.
Choose the grid and step size, then evaluate every point with consistent ODE settings. Record solver failures and verify that all parameter combinations are valid.

Calculate and summarize effects.
Use the signed output difference divided by the signed coordinate change. Report $\mu$, $\mu^*$, and $\sigma$ for each parameter, with the input and output scales stated.

Use the results as a screening decision.
Check numerical stability and ranking stability across trajectories. Prioritize follow-up for parameters with substantial effects or strongly varying effects; do not infer irrelevance from one small summary alone.

Check Your Understanding

Answer each question, then compare with the solutions below.

For $p=6$, calculate the Morris step size $\delta$. How many model evaluations are in one trajectory for $k=8$ parameters?
A parameter has effects $3,-3,0$ in output units per unit $u$. Calculate $\mu$ and $\mu^*$. Which better describes its average effect magnitude?
Why can a large $\sigma_i$ be consistent with either nonlinearity or interactions?
A solver’s output changes by $0.002$ units when its tolerance is tightened, while a Morris step changes the output by $0.001$ units. What diagnostic concern should be raised?
What is the main interpretive problem with applying an independent-factor Morris design to strongly correlated parameters?

Solutions. (1) $\delta=6/[2(6-1)]=0.6$. A trajectory changes each of eight parameters once, so it contains $8+1=9$ evaluations. (2) $\mu=(3-3+0)/3=0$, while $\mu^*=(3+3+0)/3=2$ output units per unit $u$. The mean absolute effect describes average magnitude without cancellation. (3) Effects can vary because the slope changes with the value of the parameter itself, because its effect depends on other parameter values, or because both occur; spread alone does not separate these mechanisms. (4) The solver-related change is larger than the stepwise effect, so the estimated elementary effect may be dominated by numerical error; tighten and validate solver settings before interpreting it. (5) Independent coordinate changes can create parameter combinations that are implausible under the joint distribution, so the resulting effects may not describe sensitivity over the scientifically relevant correlated input space.