Análise de ponto de mudança usando nls de R ()

Estou tentando implementar uma análise de "ponto de mudança" ou uma regressão multifásica usando nls()em R.

Aqui estão alguns dados falsos que eu criei . A fórmula que eu quero usar para ajustar os dados é:

$y = \beta_0 + \beta_1x + \beta_2\max(0,x-\delta)$

O que isso deve fazer é ajustar os dados até um certo ponto com uma certa interceptação e inclinação ( e ), depois de um certo valor x ( ), aumente a inclinação em . É disso que se trata a coisa toda. Antes do ponto , será igual a 0 e \ beta_2 será zerado.$\beta_0$$\beta_1$$\delta$$\beta_2$$\delta$$\beta_2$

Então, aqui está minha função para fazer isso:

changePoint <- function(x, b0, slope1, slope2, delta){ 
   b0 + (x*slope1) + (max(0, x-delta) * slope2)
}

E eu tento encaixar o modelo dessa maneira

nls(y ~ changePoint(x, b0, slope1, slope2, delta), 
    data = data, 
    start = c(b0 = 50, slope1 = 0, slope2 = 2, delta = 48))

Eu escolhi esses parâmetros de partida, porque eu sei que esses são os parâmetros de partida, porque criei os dados.

No entanto, eu recebo este erro:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Acabei de criar dados infelizes? Tentei ajustar isso com dados reais primeiro e estava recebendo o mesmo erro e concluí que meus parâmetros de inicialização iniciais não eram bons o suficiente.

(No começo eu pensei que poderia ser um problema resultante do fato de que maxnão é vetorizado, mas isso não é verdade It. Não torná-lo uma dor de trabalhar com Changepoint, portanto a seguinte modificação:

changePoint <- function(x, b0, slope1, slope2, delta) { 
   b0 + (x*slope1) + (sapply(x-delta, function (t) max(0, t)) * slope2)
}

Esta postagem da lista de discussão da R-help descreve uma maneira pela qual esse erro pode resultar: o rhs da fórmula é super-parametrizado, de modo que a alteração de dois parâmetros em conjunto dá o mesmo ajuste aos dados. Não vejo como isso é verdade no seu modelo, mas talvez seja.

Em qualquer caso, você pode escrever sua própria função objetivo e minimizá-la. A função a seguir fornece o erro quadrado dos pontos de dados (x, y) e um determinado valor dos parâmetros (a estrutura de argumento estranha da função é responsável por explicar como optimfunciona):

sqerror <- function (par, x, y) {
  sum((y - changePoint(x, par[1], par[2], par[3], par[4]))^2)
}

Então dizemos:

optim(par = c(50, 0, 2, 48), fn = sqerror, x = x, y = data)

E veja:

$par
[1] 54.53436800 -0.09283594  2.07356459 48.00000006

Observe que, para meus dados falsos ( x <- 40:60; data <- changePoint(x, 50, 0, 2, 48) + rnorm(21, 0, 0.5)), existem muitos máximos locais, dependendo dos valores iniciais dos parâmetros que você fornecer. Suponho que, se você quisesse levar isso a sério, ligaria para o otimizador muitas vezes com parâmetros iniciais aleatórios e examinaria a distribuição dos resultados.

Só queria acrescentar que você pode fazer isso com muitos outros pacotes. Se você deseja obter uma estimativa da incerteza em torno do ponto de mudança (algo que o nls não pode fazer), tente o mcppacote.

# Simulate the data
df = data.frame(x = 1:100)
df$y = c(rnorm(20, 50, 5), rnorm(80, 50 + 1.5*(df$x[21:100] - 20), 5))

# Fit the model
model = list(
  y ~ 1,  # Intercept
  ~ 0 + x  # Joined slope
)
library(mcp)
fit = mcp(model, df)

Vamos plotá-lo com um intervalo de previsão (linha verde). A densidade azul é a distribuição posterior para o local do ponto de mudança:

# Plot it
plot(fit, q_predict = T)

Você pode inspecionar parâmetros individuais em mais detalhes usando plot_pars(fit)e summary(fit).