Identificabilidade em um problema de regressão não linear

Suponha que eu esteja trabalhando com o seguinte modelo

$y_i = \alpha(1-\exp(-\beta t_i))+\gamma(1-\exp(-\delta t_i)) + \varepsilon_i$ .

O é um gaussiano com média zero e estou tentando encontrar os melhores valores de ajuste de . $\varepsilon_i$ $\alpha,\beta,\gamma,\delta$

Para concretização, digamos que este seja um modelo para a quantidade total de algumas espécies bacterianas com duas subespécies que crescem no tempo de acordo com o primeiro e o segundo termos no RHS, mas apenas medimos a população total. Nota: essa não é a configuração real, mas é suficiente para a pergunta.

O modelo não é identificável no sentido usual, pois sempre posso trocar e por exemplo, e obter exatamente a mesma densidade / probabilidade. $\alpha$ $\gamma$

Como você pode esperar, se eu executar um MCMC, acabo com posteriores terrivelmente amplos, e qualquer abordagem de mínimos quadrados não lineares é incrivelmente sensível às suposições iniciais - temos um grande platô na função de probabilidade.

Um melhor projeto experimental não é uma opção nesta fase - obviamente medir a subespécie separadamente seria a melhor opção.

Existe algo que eu possa fazer com esse problema ou o melhor design experimental é a única opção?

regression nonlinear-regression identifiability Vai
fonte

Por que você tem dois termos com a mesma estrutura? Como você disse, eles são identificados incorretamente.

Anna SDTC

Eles representam duas 'sub-quantidades' da quantidade medida. Não formulei esse modelo, recebi a tentativa de ajustar alguns dados e as pessoas que me deram disseram que não querem alterá-lo. Eu queria ter certeza de que não havia nada que eu pudesse fazer antes de dizer a eles que não era possível.

Will

Eu não acho que você pode com base neste modelo. Você pode fazer algum tipo de estimativa em duas etapas, na qual na primeira etapa existe apenas um desses termos exponenciais e na segunda etapa você calcula as "subquantidades"?

Anna SDTC

Isso pode funcionar. Definitivamente, preciso corrigir alguns parâmetros para chegar a algum lugar e isso pode ser uma boa maneira de fazê-lo. E se eu colocar um exponencial anterior (com o parâmetro 'big' rate) em ou (eles não são negativos)? Pensando ao longo das linhas de tentar regularizar alguma forma ...

α

$\alpha$

γ

$\gamma$

Will

Para informações gerais, como métodos de mínimos quadrados não lineares, existe "Levenberg-Marquardt" (LM) e o "reflexo da região de confiança" (TRR), você já tentou um deles? Eu sei que o TRR permite que você defina limites nos parâmetros, enquanto o LM não. Dito isto, para o seu problema, acho que você precisa ter algo mais que permita distinguir os dois termos. Você pode identificar medições com apenas uma população e fazer uma estimativa em duas etapas, como Anna observou, ou alterar seu modelo para dar conta das especificidades de cada subespécie. Caso contrário, não vejo a utilidade desse modelo.

Apr

Não há problema de identificabilidade, exceto no sentido trivial de que qualquer modelo em particular pode ter duas descrições. O problema real parece ser a dificuldade de ajustar o modelo - mas isso se deve à maneira como os modelos são parametrizados, e não à falta de identificabilidade.

Esse problema tem uma solução igualmente trivial: declare, sem qualquer perda de generalidade, que . Se você quer ser realmente exigente, insista também que, se , então . $\beta \ge \delta$ $\beta=\delta$ $\alpha \ge \gamma$

Infelizmente, isso requer qualquer procedimento para ajustar o modelo a respeitar essas restrições. A introdução de uma restrição aqui não é tão ruim, porque o aplicativo é tal que obviamente todos os parâmetros não são negativos de qualquer maneira: o espaço do parâmetro já possui limites nítidos. A inclusão de mais uma restrição não forçará nenhuma alteração na forma como ajustamos o modelo.

Um método conhecido para converter uma otimização restrita em uma irrestrita é remetereter o problema para que, no novo espaço de parâmetro, os limites sejam empurrados para o infinito. Existem muitas maneiras de conseguir isso aqui. Uma consideração do significado dos parâmetros nos guiará. Em particular, é o máximo atingido pela função para . Dado , então necessariamente e $\nu = \alpha + \gamma$

t \to g (t; α, β, γ, δ) = α (1 - e^{- β t}) + γ (1 - e^{- δ t})

$t\to g(t; \alpha,\beta,\gamma,\delta) = \alpha\left(1 - e^{-\beta t}\right) + \gamma\left(1 - e^{-\delta t}\right)$

t \geq 0

$t \ge 0$

ν

$\nu$

0 \leq α \leq ν

$0 \le \alpha\le \nu$

γ = ν - α

$\gamma = \nu - \alpha$ . Quando valores não negativos somam um todo fixo, geralmente funciona para parametrizar suas proporções do todo em termos de ângulos: seja uma proporção o cosseno ao quadrado e a outra o seno ao quadrado. Além disso, uma maneira simples de garantir que , e sejam positivos é torná-los exponenciais - ou seja, use seus logaritmos como parâmetros. Finalmente, para aplicar , defina como o cosseno quadrado de alguns ângulos vezes . Assim, podemos remeterar o problema ajustando a função

ν

$\nu$

β

$\beta$

δ

$\delta$

δ \leq β

$\delta \le \beta$

δ

$\delta$

β

$\beta$

t \to f (t; n, a, b, d) = e^{n} (1 - \cos (a)^{2} \exp (- e^{b} t) - \sin (a)^{2} \exp (- e^{b \cos (d)^{2}} t)) .

$t \to f(t;n,a,b,d) = e^n\left(1 - \cos(a)^2\exp\left(-e^{b} t\right) - \sin(a)^2\exp\left(-e^{b\cos(d)^2} t\right)\right).$

A partir de estimativas desses parâmetros (que, a propósito, não são "identificáveis" devido à ambiguidade nos ângulos e ), é possível recuperar os originais como $a$ $d$

\begin{aligned} α & = e^{n} \cos (a)^{2} \\ β & = e^{b} \\ γ & = e^{n} \sin (a)^{2} \\ δ & = e^{b} \cos (d)^{2} . \end{aligned}

$\eqalign{ \alpha &= e^n\cos(a)^2 \\ \beta &=e^b \\ \gamma &=e^n\sin(a)^2 \\ \delta &= e^b\cos(d)^2. }$

As propriedades das funções exponencial e trigonométrica asseguram todas as restrições: , e . (Como os flutuadores de precisão dupla podem se tornar astronomicamente pequenos, não há distinção prática entre e nessas restrições.) $\alpha \gt 0$ $\beta \ge \delta \gt 0$ $\gamma \gt 0$ $\gt$ $\ge$

Nesse sentido bem definido, o modelo é identificável, mesmo que os parâmetros usados para ajustá-lo não sejam identificáveis.

Embora se possa usar o MCMC, se o objetivo é apenas ajustar a curva, é mais simples usar um solucionador numérico como Newton-Raphson. O truque é encontrar um bom valor inicial . O máximo de seria uma leve superestimação de ; então comece, talvez com . Você pode começar com , supondo que cada componente faça uma contribuição substancial ao todo. Faça algumas suposições razoáveis sobre e base nas taxas esperadas de decaimento. Por exemplo, se o intervalo de for razoável, considere como uma fração da maior $y_i$ $e^n$ $n=\log(\max(y_i)/2)$ $a=\pi/4$ $e^b$ $e^d$ $t$ $b$ $t$ talvez, escolha arbitrariamente ; talvez use um valor inicial menor. (Muitas vezes você vai ter diferentes valores para as estimativas de parâmetros, dependendo estas escolhas, mas normalmente eles não afetará significativamente o função em si .) $d=\pi/4$ $f$

Em muitas circunstâncias, essa abordagem funciona muito bem. Exceto quando a variação dos erros é do mesmo tamanho que ou maior (onde será difícil discernir qualquer sinal sem uma grande quantidade de dados), o ajuste funciona mesmo com pequenas quantidades de dados: todos precisa é de quatro. $\max{y_i}$

Observe que, independentemente de como o modelo é adequado, normalmente haverá uma grande incerteza nos parâmetros: essa família de curvas é essencialmente uma pequena perturbação da família exponencial de dois parâmetros . Em muitas circunstâncias, então, dois dos parâmetros (correspondentes à amplitude e à maior taxa de decaimento ) podem ser identificados com razoável precisão, mas os outros dois parâmetros, que refletem pequenas variações dessa forma exponencial, geralmente são altamente incertos. $t\to Ae^{-Bt}$ $A$ $B$

A figura mostra um exemplo de um ajuste desafiador. A curva subjacente é mostrada em preto. Por fim, atinge um máximo de , muito lentamente. Apenas pontos de dados estão disponíveis, plotados como pontos cinza. O desvio padrão dos erros aleatórios é , uma proporção considerável desse máximo. Muitos dos erros foram positivos, fazendo com que a curva ajustada em vermelho seja um pouco maior. Os dois componentes exponenciais da curva ajustada são mostrados como linhas cinza tracejadas e pontilhadas. Um mostra um rápido aumento para um limiar de no momento ; o outro reflete o outro exponencial subindo ao seu limite de $4/3$ $24$ $1/2$ $1/3$ $t=1$ $1$ . (Você terá pouca esperança de reproduzir esse "ombro" afiado próximo de até ter pontos de dados ou mais: experimente variando no código abaixo.) $t=1$ $1000$ n

Seu sucesso em qualquer problema específico dependerá da magnitude dos erros; a faixa de valores de que são amostrados; como esses valores são espaçados; quantos valores estão disponíveis; e escolha dos valores iniciais. No entanto, este parece ser um problema tratável em geral, com soluções que podem ser obtidas rapidamente. Além disso, qualquer ajustador de probabilidade máxima procederá de maneira semelhante para minimizar a soma dos quadrados dos resíduos - e, além disso, fornecerá regiões de confiança para os parâmetros. $t$

Este é o Rcódigo que usei para testar esta proposta. Ele reproduz a figura e é facilmente modificado - altere os valores das variáveis no início - para estudar dados parecidos com os que você possa ter.

#
# Describe the underlying model
#
set.seed(17)
alpha <- 1
beta <- 2
gamma <- 1/3
delta <- 1/10
sigma <- 1/2  # Error SD.
n <- 24
x.max <- 20   # Largest value of t.
#
# The original parameterization.
#
g <- function(x, alpha, beta, gamma, delta) {
  alpha * (1 - exp(-beta * x)) + gamma * (1 - exp(-delta * x))
}
#
# The re-parameterization.  `f.1` and `f.2` are the two exponential components.
#
f <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  alpha <- n*a
  gamma <- n*(1-a)
  beta <- exp(log.b)
  delta <- cos(t.d)^2 * beta
  n - alpha * exp(-beta * x) - gamma * exp(-delta * x)
}
f.1 <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  alpha <- n*a
  beta <- exp(log.b)
  alpha * (1 - exp(-beta * x))
}
f.2 <- function(x, nu, t.a, log.b, t.d) {
  n <- exp(nu)
  a <- cos(t.a)^2
  gamma <- n*(1-a)
  beta <- exp(log.b)
  delta <- cos(t.d)^2 * beta
  gamma * (1 - exp(-delta * x))
}
#
# The objective to minimize is the mean squared residual.
# This is equivalent to finding the MLE for Gaussian errors.
#
obj <- function(theta, x, y) {
  crossprod(y - f(x, theta[1], theta[2], theta[3], theta[4])) / length(x)
}
#
# Create data and plot them.
#
x <- seq(0, x.max, length.out=n)
y <- g(x, alpha, beta, gamma, delta) + rnorm(length(x), 0, sigma)
plot(x,y, pch=16, col="#00000040", xlab="t")
#
# Fit the curve.
#
theta <- c(nu=log(max(y)/2), t.a=pi/4, log.b=log(max(x)/10), t.d=pi/4)
fit <- nlm(obj, theta, x=x, y=y, gradtol=1e-14)
theta.hat <- fit$estimate
#
# Plot relevant curves.
#
curve(g(x, alpha, beta, gamma, delta), add=TRUE, lwd=2)
curve(f(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Red", lwd=2)
curve(f.1(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Gray", lty=2, lwd=2)
curve(f.2(x, theta.hat[1], theta.hat[2], theta.hat[3], theta.hat[4]),
      add=TRUE, col="Gray", lty=3, lwd=2)

whuber
fonte

Identificabilidade em um problema de regressão não linear

Respostas: