Regression stats question

Question

I am currently focusing on developing a way to predict the onset of sexual maturity in lake sturgeon by using the elemental composition in the fin ray.

Lake sturgeon spawn later in life (females around age 15-25 and males age 12-20).

I have run Wilcoxon signed ranked tests on elemental levels before/after age 20 in females and before/after age 15 in males (Table 1.) and have found that most divalent ion levels are significantly higher after age 20 or 15.

I believe that these ions are being taken up at a greater rate (after onset of sexual maturity) because they are being used in development of gametes (and therefore are needed more) or because they are replacing calcium at a higher rate because Ca is being used for gametes (or a mixture of the two).

I have attached a plot that shows Mg, Mn, and Ba across the fin ray of a female lake sturgeon. The dotted line is age 20 and the blue dots are where I believe she became sexually mature. Is there a way to run some type of regression to indicate/prove that is when those elements began to increase? I would appreciate any suggestions or advice that you have.

Answer 1

This is a type of change point analysis. It's often used to evaluate shifts in time series, but you can also apply it to shifts over space as in your data. (Actually, as there evidently is a mapping from distance along the fin to time, you could think of your data as a time series too.)

The problem that you face is that there doesn't seem to be any simple expression that describes the shapes of these elemental profiles between potential discontinuities. So you need a method that allows for generic shape-fitting along with detection of one (or more, it seems from your data) change points. This page outlines ways to combine LOESS smoothing and change-point detection, with links to further reading and resources. Other threads on change-point methods might also be helpful.

Answer 2

Here's a solution using the R package mcp. Say we simplify your problem as a plateau followed by an abrupt change to a negative slope, followed by a plateau. Let's simulate some data:

data = data.frame(
  x = 1:190,
  distance = c(
    rnorm(100, 30, 4),  # Plateau
    rnorm(40, 60 - (1:40)*0.6, 4),  # Intercept + negative slope
    rnorm(50, 30, 4)  # P
  )
)

Now let's try and recover the onset of segment 2! We list the formulas for the segments:

model = list(
  distance ~ 1,  # Plateau
  ~ 1 + x,  # Intercept change + slope
  ~ 0  # Plateau
)

Let's fit it. This runs MCMC sampling and this case needs longer-than-usual adaptation:

fit = mcp(model, data, adapt = 5000)
plot(fit)

The gray areas on the x-is are posterior distributions of the change point locations. You're looking for the first. See cp_1 (change point 1) in summary(fit) or fixef(fit).

Refining the model

You can expand the model to include more segments, autocorrelation (e.g., an ar(1) term), other formulas, etc. You can also impose restrictions on the parameters. E.g., that Intercept_2 has to be greater than Intercept 1 and that the slope in segment 2 has to be negative (read more here):

prior = list(
  Intercept_2 = "dnorm(30, 10) T(Intercept_1, )",
  x_2 = "dnorm(0, 5) T(, 0)"
)
fit_refined = mcp(model, data, prior)

This was fit on a single series. The upcoming version of mcp (v0.4) will handle multivariate models, so you can infer a common change point across series.