Paper 3

Holding sex fixed, each additional gram of weight is associated with an estimated $0.03378 ng/dL$ increase in mean testosterone concentration, according to the fitted multiple linear regression.

Solution

At the same weight, male mice have an estimated mean testosterone concentration $7.12162 ng/dL$ higher than female mice; the sex regression coefficient is an adjusted difference between the groups.

Solution

At a conventional 5% significance level, there is evidence for an association with sex after adjusting for weight, because its two-sided p-value is

0.00131 < 0.05

. There is insufficient evidence for a weight association after adjusting for sex, because its p-value is

0.62392 > 0.05

These Student t-tests test whether the corresponding regression coefficients are zero within the additive multiple linear regression. Failure to reject the weight hypothesis does not establish an exactly zero effect; these associations alone also do not justify a causal inference. The later model with a statistical interaction asks a different question, allowing the weight association to differ by sex.

e

Solution

Two assumptions of the usual normal linear model used for the reported coefficient Student t-tests are:

The errors are independent random variables across mice.
Conditional on the predictors, each error has a normal distribution with mean zero and the same variance $σ^{2}$ : $ε_{i} ∣ w_{i}, s_{i} \sim N (0, σ^{2})$ .

The second statement includes homoscedasticity. Normality concerns the conditional errors, rather than requiring the weights themselves to have a normal distribution.

f

Solution

Include both main effects and their interaction term:

fit_interaction <- lm(testosterone ~ weight * sex, data = mice)

Here mice denotes the data frame containing the named variables. In this R formula, weight * sex expands to weight + sex + weight:sex; it does not omit the main effects.

g

Solution

The interaction term estimates that the weight-associated change in mean testosterone concentration is $0.40688 ng/dL$ per gram greater in males than females: the female slope is

- 0.04518

and the male slope is

- 0.04518 + 0.40688 = 0.36170 ng/dL

per gram.

h

Solution

Use Poisson regression for the offspring count, with testosterone as a predictor. It models a nonnegative integer response through its conditional mean. If the conditional variance substantially exceeds that mean, a negative binomial regression is a possible alternative for overdispersion.

i

Solution

Use the log link:

lo g λ_{i} = β_{0} + β_{T} T_{i}, λ_{i} = E [N_{i} ∣ T_{i}],

where

N_{i}

is the offspring count and

T_{i}

is testosterone concentration. This link function keeps the fitted conditional mean positive.

j

Solution

A one-unit increase in testosterone concentration, here

1 ng/dL

, multiplies the expected offspring count by $e^{β_{T}}$ under the Poisson regression with a log link. Thus the exact percentage change in the conditional mean is $100 (e^{β_{T}} - 1) %$ . The regression coefficient

β_{T}

itself is the change in the logarithm of the expected count, not an additive change in offspring number.

2

a

Solution

Each observation has normal distribution probability density function

f (x_{i} ∣ μ) = \frac{1}{2 π} exp [- \frac{1}{2} (x_{i} - μ)^{2}] .

The observations are independent and identically distributed, so their joint probability density function, viewed as a function of the unknown mean, is the product:

L (μ) = i = 1 \prod n f (x_{i} ∣ μ) = (2 π)^{- n /2} exp [- \frac{1}{2} i = 1 \sum n (x_{i} - μ)^{2}] .

Independence justifies multiplication; the common mean and known unit variance justify using this same density factor for every observation.

b

i

Solution

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

The estimate requested is the sample mean, $\overset{x}{ˉ} = n^{- 1} i \sum x_{i}$ . The named electronic file Iris_data.csv is not present in the supplied source archive or this repository, and its measurements are not printed in the paper. Consequently its numerical sample mean cannot be recovered from the paper alone. The familiar built-in iris data must not be silently substituted: the question does not identify those observations as its sample.

The following R code reads the actual file and prints its sample mean. The later code blocks in Question 2 continue this same script and use the resulting x, n and xbar.

dat <- read.csv("Iris_data.csv", check.names = FALSE)
# Recognize common labels, or a single measurement column.
keys <- tolower(gsub("[^[:alnum:]]", "", names(dat)))
j <- which(keys %in% c("sepallength", "sepallengthcm", "x"))
if (length(j) == 0L && ncol(dat) == 1L) j <- 1L
if (length(j) != 1L) {
  stop("Select the sepal-length column explicitly after inspecting names(dat).")
}
x <- dat[[j]]
stopifnot(is.numeric(x), length(x) > 0L, all(is.finite(x)))
n <- length(x)
xbar <- mean(x)
cat("n =", n, " sample mean =", format(xbar, digits = 12), "\n")

This reader assumes a header row and one observation per row. An alternative layout should be identified explicitly rather than dropping measurements or treating an identifier column as a length.

ii

Solution

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

The identity

\sum_{i} (x_{i} - μ)^{2} = \sum_{i} (x_{i} - \overset{x}{ˉ})^{2} + n (μ - \overset{x}{ˉ})^{2}

shows that the likelihood function is a bell-shaped curve in

μ

, maximized at the sample mean, with width

1/ n

. Its actual vertical scale also depends on the within-sample sum of squares.

Continue the R script with:

sse <- sum((x - xbar)^2)
loglik <- function(mu) {
  -n/2 * log(2*pi) - 0.5 * (sse + n * (mu - xbar)^2)
}
grid <- seq(3, 6, length.out = 2001L)
ell <- loglik(grid)
png("q2-likelihood.png", width = 900, height = 600)
par(mar = c(4.5, 4.5, 4, 1))
if (max(ell) > log(.Machine$double.xmin)) {
  plot(grid, exp(ell), type = "l", xlab = expression(mu),
       ylab = expression(L(mu)), main = "Likelihood on [3,6]")
} else {
  # Explicitly label a proportional likelihood if the true scale underflows.
  plot(grid, exp(ell - max(ell)), type = "l", xlab = expression(mu),
       ylab = "Likelihood relative to maximum on plotted grid",
       main = "Likelihood on [3,6] (rescaled)")
}
abline(v = xbar, col = "steelblue", lty = 2)
dev.off()

Working first with the log-likelihood avoids premature numerical underflow. The fallback preserves the shape on the requested interval and labels the altered scale. A data-specific plot requires the missing CSV file; no fabricated observation set is used here.

c

Solution

The exponential distribution prior has rate

λ = 1/3

and support

μ \geq 0

, so its prior density is

(1/3) e^{- μ /3} 1_{{μ \geq 0}}

. By Bayes' theorem,

p (μ ∣ x) = C exp [- \frac{1}{2} i \sum (x_{i} - μ)^{2} - \frac{μ}{3}] 1_{{μ \geq 0}},

where

C

is a normalizing constant independent of

μ

Completing the square gives a useful exact reference for the numerical sampling. Put

m = \overset{x}{ˉ} - \frac{1}{3 n}, s = \frac{1}{n} .

Then this normal mean posterior with an exponential prior is a normal distribution with location

m

and variance

s^{2}

, truncated below at zero. With

ϕ

the standard normal distribution density and

Φ

its cumulative distribution function,

p (μ ∣ x) = \frac{ϕ (( μ - m ) / s )}{s Φ ( m / s )} 1_{{μ \geq 0}} .

(3)

The exponential rate is

1/3

, not three; the prior mean is three. The support restriction must be retained in both the formula and the sampler.

d

Solution

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Use the random-walk Metropolis algorithm with proposal

μ^{'} = μ + ε

, where

ε \sim N (0, 1)

. The proposal is symmetric, so for a nonnegative proposal the acceptance probability is

a (μ, μ^{'}) = min {1, exp [ℓ_{p} (μ^{'}) - ℓ_{p} (μ)]}, ℓ_{p} (u) = - \frac{1}{2} \sum_{i} (x_{i} - u)^{2} - \frac{u}{3} .

Negative proposals are rejected. They are not redrawn until positive, since that would change the proposal kernel and require a different acceptance ratio. On rejection, store the current state again; keeping accepted moves alone would bias the empirical posterior distribution.

This R continuation makes exactly 10,000 proposals, starting at

μ_{0} = 4

set.seed(20240528)
B <- 10000L
mu <- numeric(B + 1L)
mu[1L] <- 4
logpost <- function(u) {
  if (!is.finite(u) || u < 0) return(-Inf)
  loglik(u) - u/3
}
lp <- logpost(mu[1L])
accepted <- 0L
for (i in seq_len(B)) {
  proposal <- rnorm(1L, mean = mu[i], sd = 1)
  lp_new <- logpost(proposal)
  if (log(runif(1L)) < min(0, lp_new - lp)) {
    mu[i + 1L] <- proposal
    lp <- lp_new
    accepted <- accepted + 1L
  } else {
    mu[i + 1L] <- mu[i]
  }
}
cat("Acceptance fraction =", accepted/B, "\n")
# A provisional discard; assess it with the trace in part (e).
burn <- 1000L
post <- mu[(burn + 2L):(B + 1L)]
m <- xbar - 1/(3*n)
s <- 1/sqrt(n)
posterior_density <- function(u) {
  ans <- numeric(length(u))
  ok <- u >= 0
  ans[ok] <- exp(dnorm(u[ok], m, s, log = TRUE) -
                   pnorm(m/s, log.p = TRUE))
  ans
}
png("q2-posterior.png", width = 900, height = 600)
hist(post, breaks = 40, probability = TRUE, xlab = expression(mu),
     main = "Retained Metropolis draws and exact posterior")
u <- seq(0, max(post) + 3*s, length.out = 4001L)
lines(u, posterior_density(u), col = "firebrick", lwd = 2)
dev.off()

The exact truncated normal distribution density overlays the histogram as a check; it does not replace the requested Markov chain Monte Carlo procedure. The 9,000 retained states exclude the initial state and the first 1,000 transitions. The particular histogram and acceptance fraction remain dependent on the missing measurements.

e

Solution

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Plot all 10,001 stored states, including the initial value, and indicate the provisional burn-in boundary. Also display the retained chain's autocorrelation function:

png("q2-trace.png", width = 1000, height = 700)
par(mfrow = c(2, 1), mar = c(4, 4.5, 2.5, 1))
plot(0:B, mu, type = "l", xlab = "Iteration", ylab = expression(mu),
     main = "Complete Metropolis trace")
abline(v = burn, col = "firebrick", lty = 2)
abline(h = m, col = "steelblue", lty = 3)
acf(post, main = "Autocorrelation after provisional burn-in")
dev.off()
# Check sensitivity to the chosen discard and to the retained time window.
for (discard in c(500L, 1000L, 2000L)) {
  kept <- mu[(discard + 2L):(B + 1L)]
  cat("Discard", discard, ": posterior mean", mean(kept), "\n")
}
cat("Means of retained halves:", mean(post[1:4500]),
    mean(post[4501:9000]), "\n")

A satisfactory trace plot should show an initial transient, if any, followed by repeated exploration of one stable range with no continuing drift. Repeated values are expected because rejections are stored. If

n

is large, the prescribed proposal standard deviation one can be much wider than the posterior scale

1/ n

, giving low acceptance and conspicuous plateaus. A flat-looking trace is therefore not by itself evidence of good exploration.

Compare the retained histogram with the exact posterior density, inspect autocorrelation, compare separate retained windows, and assess sensitivity to the discarded count. Agreement supports the adequacy of this run; neither 10,000 iterations nor a chosen burn-in of 1,000 proves convergence. The effective sample size of a Markov chain can be much smaller than the number of stored states.

The actual trace cannot be diagnosed without Iris_data.csv. The code produces the required trace when the measurements are available; the comments here specify what to examine rather than asserting an unseen chain has converged.

f

Solution

TODO: Complete the data-specific numerical results and diagnostics when the original Iris_data.csv is available.

Estimate the posterior mean by averaging the retained states:

E [μ ∣ x] = \frac{1}{9000} t = 1001 \sum 10000 μ_{t} .

The exact truncated normal distribution reference is

E [μ ∣ x] = m + s \frac{ϕ ( m / s )}{Φ ( m / s )}, m = \overset{x}{ˉ} - \frac{1}{3 n}, s = n^{- 1/2} .

For a posterior well separated from zero, the truncation correction is tiny and the posterior mean is approximately

\overset{x}{ˉ} - 1/ (3 n)

, slightly below the sample mean.

The R continuation prints the Monte Carlo estimator, its approximate batch-based Monte Carlo standard error, and the exact value for comparison:

posterior_mean_mc <- mean(post)
posterior_mean_exact <- m + s * exp(dnorm(m/s, log = TRUE) -
                                    pnorm(m/s, log.p = TRUE))
# Nonoverlapping batches account approximately for serial dependence.
batch_length <- floor(sqrt(length(post)))
number_batches <- floor(length(post)/batch_length)
used <- post[seq_len(batch_length * number_batches)]
batch_averages <- colMeans(matrix(used, nrow = batch_length))
mcse <- sd(batch_averages)/sqrt(number_batches)
cat("Posterior mean, MCMC =", posterior_mean_mc,
    " approximate MCSE =", mcse, "\n")
cat("Posterior mean, exact =", posterior_mean_exact, "\n")

The batch length must be long enough relative to the chain's dependence for this uncertainty estimate to be reliable. No numerical data-specific posterior mean is determined by the supplied paper alone.

3

a

Solution

Treat

X \geq 0

as a prescribed input protein concentration, and

Y (t), Z (t) \geq 0

as the downstream gene-product concentrations. Define the Hill activation function

h (u; K, n) = \frac{u ^{n}}{K ^{n} + u ^{n}}, K > 0, n > 0.

A minimal model with no basal expression is

\frac{d Y}{d t} = a_{Y} h (X; K_{X Y}, n_{X Y}) - Y, \frac{d Z}{d t} = a_{Z} h (X; K_{XZ}, n_{XZ}) - Z .

Here

a_{Y}, a_{Z} > 0

are maximal production rates,

K_{X Y}, K_{XZ}

are Hill constants with the same units as

X

, and

n_{X Y}, n_{XZ}

are dimensionless Hill coefficients. The negative terms describe a first-order reaction for loss of each gene product, with both rate constants set to one in the chosen time units. Thus

a_{Y}, a_{Z}

have units of concentration per time and the unit loss rates have units of inverse time.

The equations model regulation of gene expression at the level of its effective protein output, without separate messenger RNA dynamics. The independent branches mean that neither

Y

nor

Z

regulates the other in this initial model. Additive basal production could be included, but is not required by the question and is set to zero throughout this example.

b

Solution

For a constant

X

, set both ordinary differential equation derivatives to zero:

Y^{*} (X) = a_{Y} \frac{X ^{n_{X Y}}}{K _{X Y}^{n_{X Y}} + X ^{n_{X Y}}}, Z_{direct}^{*} (X) = a_{Z} \frac{X ^{n_{XZ}}}{K _{XZ}^{n_{XZ}} + X ^{n_{XZ}}} .

These equilibrium concentrations are unique for each input. For a loss rate

γ

other than one, divide the corresponding production term by

γ

Choose

a_{Y} = 2

K_{X Y} = 1

n_{X Y} = 2

and

a_{Z} = 3

K_{XZ} = 2

n_{XZ} = 4

. The following Python code plots the requested equilibrium concentrations and defines the functions used in part (c):

import numpy as np
import matplotlib.pyplot as plt

def hill(u, K, n):
    u = np.asarray(u)
    return u**n / (K**n + u**n)

X_grid = np.linspace(0, 5, 501)
Y_eq = 2 * hill(X_grid, 1, 2)
Z_eq = 3 * hill(X_grid, 2, 4)
plt.figure(figsize=(6, 4))
plt.plot(X_grid, Y_eq, label="Y equilibrium: K=1, n=2")
plt.plot(X_grid, Z_eq, label="Z equilibrium: K=2, n=4")
plt.xlabel("X concentration")
plt.ylabel("Equilibrium concentration")
plt.legend()
plt.tight_layout()
plt.show()

The Hill constant sets the half-maximum input:

Y^{*} (1) = 1

and

Z_{direct}^{*} (2) = 1.5

. Increasing that constant shifts activation toward larger

X

. The Hill coefficient changes steepness rather than the half-maximum position; the normalized slope at

X = K

n / (4 K)

. Because the two examples also differ in

K

and maximal production, compare normalized responses when isolating the effect of

n

alone. Both curves vanish at zero and saturate toward their respective maxima.

Figure 1.
Independent Hill activation and convergence from two initial states
. The left panel shows different half-activation inputs and Hill coefficients. The right panel integrates the two independent output equations at fixed X=2; dashed curves start above equilibrium and solid curves start at zero. **Details**.

c

Solution

At a fixed input the production terms are constant, so both linear ordinary differential equations can also be solved exactly:

Y (t) = Y^{*} + [Y (0) - Y^{*}] e^{- t}, Z (t) = Z_{direct}^{*} + [Z (0) - Z_{direct}^{*}] e^{- t} .

Each displacement from its equilibrium concentration decays with time constant one. The Jacobian matrix is

- I

, with two eigenvalues equal to

- 1

, so the fixed point is globally asymptotically stable for the nonnegative state space.

Here is a numerical solver using the classical fourth-order Runge-Kutta method; it uses only the repository's existing NumPy and matplotlib dependencies. Continue the definitions from part (b):

def integrate(rhs, initial, times):
    out = np.empty((len(times), len(initial)))
    out[0] = initial
    for j, dt in enumerate(np.diff(times)):
        state = out[j]
        k1 = rhs(state)
        k2 = rhs(state + dt*k1/2)
        k3 = rhs(state + dt*k2/2)
        k4 = rhs(state + dt*k3)
        out[j+1] = state + dt*(k1 + 2*k2 + 2*k3 + k4)/6
    return out

def direct_rhs(state, X):
    Y, Z = state
    return np.array([2*hill(X, 1, 2) - Y,
                     3*hill(X, 2, 4) - Z])

times = np.linspace(0, 8, 801)
plt.figure(figsize=(7, 4))
for initial, style in [([0, 0], "-"), ([3, 4], "--")]:
    trajectory = integrate(lambda s: direct_rhs(s, 2), initial, times)
    exact = np.array([1.6, 1.5]) + (np.array(initial) - [1.6, 1.5]) \
            * np.exp(-times[:, None])
    assert np.max(np.abs(trajectory - exact)) < 1e-8
    plt.plot(times, trajectory[:, 0], style, label=f"Y, initial {initial}")
    plt.plot(times, trajectory[:, 1], style, label=f"Z, initial {initial}")
plt.xlabel("Time")
plt.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()

For

X = 2

, $Y^{*} = 1.6, Z_{direct}^{*} = 1.5$ . The initial state

(0, 0)

approaches those values from below, while

(3, 4)

approaches from above. Over eight time units the initial displacement is reduced by

e^{- 8} \approx 0.000335

, explaining why both trajectories nearly coincide with the equilibria from part (b). The two-panel figure above displays these simulations.

d

Solution

Add the activating link from

Y

Z

and multiply the two Hill activation functions for the target:

\frac{d Y}{d t} = a_{Y} h (X; K_{X Y}, n_{X Y}) - Y, \frac{d Z}{d t} = a_{Z} h (X; K_{XZ}, n_{XZ}) h (Y; K_{Y Z}, n_{Y Z}) - Z .

The extra Hill constant

K_{Y Z} > 0

has units of

Y

concentration and

n_{Y Z} > 0

is its Hill coefficient. Both inputs are required for substantial target production: if either is zero, the production term for

Z

vanishes. This all-activating circuit is a coherent feed-forward loop. A product of graded input functions implements the specified AND rule; it is a phenomenological rule rather than a complete molecular promoter model.

Retain the parameter values of parts (b) and (c), and set

K_{Y Z} = 1

n_{Y Z} = 4

. Continue the same Python script with the altered right-hand side:

def and_rhs(state, X):
    Y, Z = state
    return np.array([2*hill(X, 1, 2) - Y,
                     3*hill(X, 2, 4)*hill(Y, 1, 4) - Z])

times = np.linspace(0, 10, 1001)
plt.figure(figsize=(7, 4))
for X in [1, 3]:
    trajectory = integrate(lambda s: and_rhs(s, X), [0, 0], times)
    plt.plot(times, trajectory[:, 0], "--", label=f"Y, X={X}")
    plt.plot(times, trajectory[:, 1], label=f"Z, X={X}")
plt.xlabel("Time")
plt.ylabel("Concentration")
plt.legend()
plt.tight_layout()
plt.show()

X = 1

Y^{*} = 1

and

Z^{*} = (3/17) (1/2) = 3/34 \approx 0.08824

. At

X = 3

Y^{*} = 1.8

and

Z^{*} \approx 2.2873

. In both simulations

Y

first accumulates and opens the additional gate;

Z

therefore initially grows more slowly than under direct activation alone. The equilibrium concentration of

Y

is unchanged by the added outgoing interaction.

Figure 1.
AND feed-forward dynamics and the effect of activation thresholds
. The left panel compares fixed inputs X=1 and X=3 from the zero state. The right panel compares the full feed-forward equilibrium response with direct activation alone and with individual Hill constants increased. **Details**.

e

Solution

The quasi-steady-state approximation for

Y

sets

Y = a_{Y} h (X; K_{X Y}, n_{X Y})

. Substituting into the

Z

equation gives

\frac{d Z}{d t} = a_{Z} h (X; K_{XZ}, n_{XZ}) h (a_{Y} h (X; K_{X Y}, n_{X Y}); K_{Y Z}, n_{Y Z}) - Z .

Its equilibrium concentration is therefore

Z_{AND}^{*} (X) = a_{Z} \frac{X ^{n_{XZ}}}{K _{XZ}^{n_{XZ}} + X ^{n_{XZ}}} \frac{( a _{Y} X ^{n_{X Y}} / ( K _{X Y}^{n_{X Y}} + X ^{n_{X Y}} ) ) ^{n_{Y Z}}}{K _{Y Z}^{n_{Y Z}} + ( a _{Y} X ^{n_{X Y}} / ( K _{X Y}^{n_{X Y}} + X ^{n_{X Y}} ) ) ^{n_{Y Z}}} .