Quando analisamos o relacionamento entre uma resposta, \(y\), e uma covariável, \(x\), pode ser aparente que para diferentes intervalos de \(x\), diferentes relações lineares ocorram. Nesses casos, um único modelo linear pode não fornecer uma adequada descrição e um modelo não linear pode também não ser apropriado.
Piecewise linear regression (ou Regressão Linear por Partes) é uma forma de regressão que permite o ajuste de múltiplos modelos lineares aos dados para diferentes intervalos de \(x\).
Pontos de quebra (breakpoints) são os valores de \(x\) onde a inclinação da função linear muda. O valor do ponto de quebra pode ou não ser conhecido antes da análise, tipicamente sendo desconhecido e consequentemente, estimado.1
Na figura abaixo temos um comportamento que exemplifica a necessidade da utilização deste tipo de modelo.
pkg <- c("segmented",
"latticeExtra",
"RSADBE")
sapply(pkg, require,
character.only = TRUE)
Loading required package: segmented
Loading required package: latticeExtra
Loading required package: RColorBrewer
Loading required package: lattice
Loading required package: RSADBE
segmented latticeExtra RSADBE
TRUE TRUE TRUE
data(down)
xyplot(log(cases/births) ~ age,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
scales = list(x = list(draw = FALSE),
y = list(draw = FALSE)),
xlab = "x",
ylab = list(label = "y",
rot = 0),
down)
Existem diversas formas de escrever o modelo, ao percorrer este documento você irá se deparar com mais algumas.
Usualmente o modelo é escrito da seguinte maneira (considerando um único ponto de quebra, \(\theta\)):
\[ \begin{align} y = \beta_{0} + \beta_{1} x , & \quad x \leq \theta \\ y = \beta_{0} + \beta_{2} x + \theta (\beta_{1}-\beta_{2}), & \quad x > \theta \end{align} \]
Quando a covariável, \(x\), é menor ou igual ao ponto de quebra, \(\theta\), temos uma reta, \(\beta_{0} + \beta_{1} x\). Quando a covariável é maior que \(\theta\), temos uma nova reta, \(\beta_{0} + \beta_{2} x + \theta (\beta_{1}-\beta_{2})\), começando no ponto de quebra e com uma nova inclinação.
Como exemplo é utilizado um conjunto de dados2.
Uma companhia de eletrônicos periodicamente realiza grandes importações de peças utilizadas em seus produtos. O tamanho dos carregamentos varia de acordo com a produção/demanda. O conjunto de dados contém informações do custo (\(y\)) de manejo dos carregamentos até um depósito (em milhares de dólares) e o tamanho (\(x_{1}\)) dos carregamentos (em milhares de peças).
Com base numa gráfico de dispersão foi escolhido como ponto de quebra \(x = 250\).
Quando o ponto de quebra é conhecido, o modelo de regressão linear por partes pode ser formulado da seguinte maneira:
\[ \begin{align} y_{i} & = \beta_{0} + \beta_{1} x_{1i} + \beta_{2} (x_{1i} - 250) x_{2i} + \epsilon_{i}, \quad \text{alternativamente:} \\ y_{i} & = \beta_{0} + \beta_{1} x_{1i} + \beta_{2} x_{2i}^{*} + \epsilon_{i}, \text{onde} \end{align} \]
Foram ajustados dois modelos, um modelo de regressão linear simples é um modelo de regressão linear por partes.
path <- "~/Dropbox/CE092 - 2015-2/Alunos/Henrique/Piecewise Regression/shipment.txt"
data <- read.table(path,
header = TRUE)
m0 <- lm(y ~ x1,
data)
coef.m0 <- coefficients(m0)
attr(coef.m0, "names") <- c("beta0", "beta1")
data$x2 <- with(data,
ifelse(x1 <= 250, 0, x1 - 250))
m1 <- lm(y ~ x1 + x2,
data)
coef.m1 <- coefficients(m1)
attr(coef.m1, "names") <- c("beta0", "beta1", "beta2")
xyplot(y ~ x1,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
xlab = expression(x[1]),
ylab = list(label = "y",
rot = 0),
scales = list(y = list(tick.number = 8)),
xlim = c(100, 425),
ylim = c(7, 16),
key = list(space = "right",
title = "Regressão Linear",
cex.title = 1.15,
text = list(c("Simples", "por Partes")),
lines = list(col = c(3, 6),
lwd = 2)),
data) +
layer(with(as.list(coef.m0),
panel.curve(beta0+beta1*x,
col = 3,
lwd = 2))) +
layer(with(as.list(coef.m1),
panel.curve(beta0+{beta1*x}*(x <= 250)+{beta1*x+beta2*(x-250)}*(x > 250),
col = 6,
lwd = 2)))
Pela retas ajustadas, a impressão que fica é que o ajuste da regressão por partes se saiu muito melhor, dado sua proximidade aos pontos observados.
r.m0 <- residuals(m0,
type = "pearson")
f.m0 <- fitted(m0)
ar.1 <- xyplot(r.m0 ~ f.m0,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = "Resíduos \n de Pearson",
rot = 0),
main = "(a) Ajuste do modelo",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m0, r.m0,
lwd = 2,
col = 2)})
ar.2 <- xyplot(sqrt(abs(r.m0)) ~ f.m0,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = expression(sqrt(abs(atop("Resíduos", "de Pearson")))),
rot = 0),
main = "(b) Relação média - variância",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m0, sqrt(abs(r.m0)),
lwd = 2,
col = 2)})
ar.3 <- qqmath(r.m0,
type = c("p", "g"),
pch = 16,
xlab = "Quantis teóricos",
ylab = list(label = "Quantis \n amostrais",
rot = 0),
main = "(c) Normalidade")
print(ar.1,
position = c(0, 0,
1/3, 1),
more = TRUE)
print(ar.2,
position = c(1/3, 0,
2/3, 1),
more = TRUE)
print(ar.3,
position = c(2/3, 0,
1, 1))
Pelos três gráficos acima mostrados, vemos uma clara não adequação aos pressupostos (a) e (b).
Nos três gráficos abaixo temos a análise de resíduos para o modelo de regressão por partes. Ela também não se mostra muito boa, contudo, está bem melhor do que no caso da regressão linear simples. Também deve ser levado em consideração o tamanho da amostra, que é composta por apenas 10 observações.
r.m1 <- residuals(m1,
type = "pearson")
f.m1 <- fitted(m1)
ar.1 <- xyplot(r.m1 ~ f.m1,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = "Resíduos \n de Pearson",
rot = 0),
main = "(a) Ajuste do modelo",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m1, r.m1,
lwd = 2,
col = 2)})
ar.2 <- xyplot(sqrt(abs(r.m1)) ~ f.m1,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = expression(sqrt(abs(atop("Resíduos", "de Pearson")))),
rot = 0),
main = "(b) Relação média - variância",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m1, sqrt(abs(r.m1)),
lwd = 2,
col = 2)})
ar.3 <- qqmath(r.m1,
type = c("p", "g"),
pch = 16,
xlab = "Quantis teóricos",
ylab = list(label = "Quantis \n amostrais",
rot = 0),
main = "(c) Normalidade")
print(ar.1,
position = c(0, 0,
1/3, 1),
more = TRUE)
print(ar.2,
position = c(1/3, 0,
2/3, 1),
more = TRUE)
print(ar.3,
position = c(2/3, 0,
1, 1))
Quando existe mais de um ponto de quebra, o procedimento empregado é o mesmo.
No exemplo aqui considerado existem dois pontos de quebra, \(15\) e \(30\). O modelo fica expresso da seguinte forma:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} (x_{i} - 15)_{+} + \beta_{3} (x_{i} - 30)_{+} + \epsilon_{i}, \quad \epsilon_{i} \overset{\text{iid}}{\sim} \text{N}(0, \sigma^{2}) \]
data("PW_Illus")
data <- PW_Illus
names(data) <- c("x1", "y")
data$x2 <- with(data,
ifelse(x1 <= 15, 0, x1 - 15))
data$x3 <- with(data,
ifelse(x1 <= 30, 0, x1 - 30))
m0 <- lm(y ~ x1+x2+x3,
data)
coef.m0 <- coefficients(m0)
attr(coef.m0, "names") <- c("beta0", "beta1", "beta2", "beta3")
xyplot(y ~ x1,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
xlab = expression(x[1]),
ylab = list(label = "y",
rot = 0),
scales = list(x = list(tick.number = 11),
y = list(tick.number = 9)),
data) +
layer(with(as.list(coef.m0),
panel.curve(beta0+{beta1*x}*(x <= 15)+
{beta1*x+beta2*(x-15)}*(x > 15 & x <= 30)+
{beta1*x+beta2*(x-15)+beta3*(x-30)}*(x > 30),
col = 6,
lwd = 2)))
No intervalo \([.5; 15]\) temos o modelo:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \epsilon_{i}. \]
No intervalo \((15; 30]\) temos o modelo:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} (x_{i} - 15) + \epsilon_{i}. \]
No intervalo \((30; 50]\) temos o modelo:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} (x_{i} - 15) + \beta_{3} (x_{i} - 30) + \epsilon_{i}. \]
r.m0 <- residuals(m0,
type = "pearson")
f.m0 <- fitted(m0)
ar.1 <- xyplot(r.m0 ~ f.m0,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = "Resíduos \n de Pearson",
rot = 0),
main = "(a) Ajuste do modelo",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m0, r.m0,
lwd = 2,
col = 2)})
ar.2 <- xyplot(sqrt(abs(r.m0)) ~ f.m0,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = expression(sqrt(abs(atop("Resíduos", "de Pearson")))),
rot = 0),
main = "(b) Relação média - variância",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m0, sqrt(abs(r.m0)),
lwd = 2,
col = 2)})
ar.3 <- qqmath(r.m0,
type = c("p", "g"),
pch = 16,
xlab = "Quantis teóricos",
ylab = list(label = "Quantis \n amostrais",
rot = 0),
main = "(c) Normalidade")
print(ar.1,
position = c(0, 0,
1/3, 1),
more = TRUE)
print(ar.2,
position = c(1/3, 0,
2/3, 1),
more = TRUE)
print(ar.3,
position = c(2/3, 0,
1, 1))
Pelos três gráficos de análise de resíduos acima mostrados, vemos uma adequação razoável aos pressupostos.
O conjunto de dados aqui estudo é do livro do Montgomery, Peck e Vining, 2001, Seção 7.2.24.
O dado consiste num experimento em que a queda de tensão da bateria (em voltagem) de um motor de míssil guiado é observada ao longo do tempo (em segundos). A queda de tensão da bateria foi observada em 41 instantes da trajetória do míssil.
Foi considerado dois pontos de quebra, \(6.5\) e \(30\) segundos. O modelo fica expresso da seguinte forma:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} (x_{i} - 6.5)_{+}^{3} + \beta_{5} (x_{i} - 13)_{+}^{3} + \epsilon_{i}, \quad \epsilon_{i} \overset{\text{iid}}{\sim} \text{N}(0, \sigma^{2}) \]
data("VD")
data <- VD
data$x2 <- with(data,
ifelse(Time <= 6.5, 0, Time - 6.5))
data$x3 <- with(data,
ifelse(Time <= 13, 0, Time - 13))
m0 <- lm(Voltage_Drop ~ Time+I(Time**2)+I(Time**3)+I(x2**3)+I(x3**3),
data)
coef.m0 <- coefficients(m0)
attr(coef.m0, "names") <- c("beta0", "beta1", "beta2", "beta3", "beta4", "beta5")
xyplot(Voltage_Drop ~ Time,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
xlab = "Tempo",
ylab = list(label = "Queda \n de \n voltagem",
rot = 0),
scales = list(x = list(tick.number = 11),
y = list(tick.number = 9)),
data) +
layer(with(as.list(coef.m0),
panel.curve(beta0+beta1*x+beta2*x**2+{beta3*x**3}*(x <= 6.5)+
{beta3*x**3+beta4*(x-6.5)**3}*(x > 6.5 & x <= 13)+
{beta3*x**3+beta4*(x-6.5)**3+beta5*(x-13)**3}*(x > 13),
col = 6,
lwd = 2)))
No intervalo \([0; 6.5]\) temos o modelo:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \epsilon_{i}. \]
No intervalo \((6.5; 13]\) temos o modelo:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} (x_{i} - 6.5)^{3} + \epsilon_{i}. \]
No intervalo \((30; 50]\) temos o modelo:
\[ y_{i} = \beta_{0} + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} (x_{i} - 6.5)^{3} + \beta_{5} (x_{i} - 13)^{3} + \epsilon_{i}. \]
r.m0 <- residuals(m0,
type = "pearson")
f.m0 <- fitted(m0)
ar.1 <- xyplot(r.m0 ~ f.m0,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = "Resíduos \n de Pearson",
rot = 0),
main = "(a) Ajuste do modelo",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m0, r.m0,
lwd = 2,
col = 2)})
ar.2 <- xyplot(sqrt(abs(r.m0)) ~ f.m0,
type = c("p", "g"),
pch = 16,
jitter.x = TRUE,
xlab = "Valores ajustados",
ylab = list(label = expression(sqrt(abs(atop("Resíduos", "de Pearson")))),
rot = 0),
main = "(b) Relação média - variância",
panel = function(...){
panel.xyplot(...)
panel.loess(f.m0, sqrt(abs(r.m0)),
lwd = 2,
col = 2)})
ar.3 <- qqmath(r.m0,
type = c("p", "g"),
pch = 16,
xlab = "Quantis teóricos",
ylab = list(label = "Quantis \n amostrais",
rot = 0),
main = "(c) Normalidade")
print(ar.1,
position = c(0, 0,
1/3, 1),
more = TRUE)
print(ar.2,
position = c(1/3, 0,
2/3, 1),
more = TRUE)
print(ar.3,
position = c(2/3, 0,
1, 1))
Pelos três gráficos de análise de resíduos acima mostrados, vemos uma adequação razoável aos pressupostos.
No R, o ajuste de uma Piecewise Regression pode ser feita com o package segmented (relacionamentos segmentados em modelos de regressão com estimação de pontos de quebra/pontos de mudança).
Com este package, somos capazes de estimar modelos de regressão com relacionamentos lineares por partes com um número fixo de pontos de quebra.
O objetivo do package segmented é estimar modelos lineares e modelos lineares generalizados (e praticamente qualquer modelo de regressão) com um ou mais relacionamentos segmentados no preditor linear. Estimação das inclinações e dos possíveis múltiplos pontos de quebra são providas. O package incluí funções de teste/estimação e métodos para retornar, resumir e imprimir gráficos dos resultados.
O algoritmo utilizado pela segmented não é uma busca num grid (grid-search). É um procedimento iterativo (Muggeo, 2003)6 que precisa de valores (chutes) iniciais apenas para o parâmetro correspondente ao ponto de quebra e por este motivo ele é bastante eficiente mesmo com vários pontos de quebra a serem estimados.
Uma outra forma de representar o modelo com uma única covariável é a seguinte:
\[ \beta_{1} x_{i} + \beta_{2} (x_{i} - \theta)_{+}, \quad \text{onde} \\ (x_{i} - \theta)_{+} = (x_{i} - \theta) \cdot I(x_{i} > \theta), \quad \text{e} \]
Da forma como o modelo está acima representado, temos um modelo não linear. Muggeo (2003) mostra uma aproximada representação linear intrínseca que, até certo ponto, nos permite traduzir o problema em uma estrutura linear padrão: dado um chute inicial para o ponto de quebra, \(\tilde{\theta}\), sucessivos segmentos são tentados para estimar o modelo pelo ajuste iterativo do modelo linear com o seguinte preditor linear:
\[ \beta_{1} x_{i} + \beta_{2} (x_{i} - \tilde{\theta})_{+} + \gamma I(x_{i} > \tilde{\theta})^{-}, \quad \text{onde} \]
A cada iteração, um modelo linear padrão é ajustado, e o valor do ponto de quebra é atualizado via \(\hat{\theta} = \tilde{\theta} + \hat{\gamma} / \hat{\beta}_{2}\); perceba que \(\hat{\gamma}\) representa o intervalo, na atual estimativa de \(\theta\), entre as duas retas ajustadas oriundas do modelo descrito. Quando o algoritmo converge, o ‘intervalo’ se torna menor, i.e., \(\hat{\gamma} \approx 0\), e o erro padrão (standard error) de \(\hat{\theta}\) pode ser obtido via método Delta para a razão \(\hat{\gamma} / \hat{\beta}_{2}\) que se reduz em \(\text{SE}(\hat{\gamma})/|\hat{\beta}_{2}|\) se \(\hat{\gamma} = 0\).
A ideia pode ser utilizada para ajustar múltiplos relacionamentos segmentados, apenas pela inclusão no preditor linear de apropriadas variáveis construídas para que os adicionais pontos de quebra possam ser estimados: em cada passo, todo ponto de quebra estimado é atualizado através dos importantes coeficientes de ‘intervalo’ e ‘diferença na inclinação’. Devido a facilidade computacional, o algoritmo é capaz de executar a estimação de múltiplos pontos de quebra de uma maneira bem eficiente.
Um outro package do R que lida com pontos de quebra é o strucchange de Zeileis et al. (2002)7. O package trabalha com modelos de regressão tendo um diferente conjunto de parâmetros para cada ‘pedaço’ da variável segmentada, tipicamente, o tempo; strucchange realiza a estimação dos pontos de quebra via um algoritmo de grid de busca dinâmico e possibilita o teste da instabilidade do parâmetro.
Alguns modelos não lineares são segmentados, tendo em alguns aspectos vantagens e facilidades em sua utilização.
Reparametrizações podem ser feitas com o intuito de facilitar a interpretabilidade dos parâmetros9. O procedimento de reparametrização pode ser sistematizado em três etapas:
- Expressar o parâmetro de interesse \(\vartheta\), como função dos elementos de \(\theta\), ou seja, \(\vartheta = g(\theta)\);
- Escolher um dos elementos \(\theta_{i}\) de \(\theta = (\theta_{i}, \theta_{i-1})\) para ser colocado em função de \(\vartheta\) de tal forma a obter \(\theta_{i} = h(\theta_{-i}, \vartheta)\);
- Substituir \(\theta_{i}\) na equação do modelo não linear pela expressão obtida no passo anterior, \(h(\theta_{-i}, \vartheta)\), fazendo as simplificações convenientes. Assim o modelo \(f(\textbf{x}, \theta)\) pode ser expresso como \(f(\textbf{x}, \theta_{-i}, \vartheta)\).
A função \(h\) é a inversa de \(g\) em \(\theta_{i}\).
Trata-se da junção de dois segmentos de reta no qual um deles é horizontal. O modelo pode ser escrito da seguinte forma:
\[ \left\{\begin{matrix} \theta_{0} + \theta_{1} x, \quad x \leq \theta_{b} \\ \theta_{0} + \theta_{1} \theta_{b}, \quad x > \theta_{b} \end{matrix}\right. \]
\[ \begin{align} \vartheta & = g(\theta) \\ \vartheta_{b} & = \theta_{0} + \theta_{1} \theta_{b} \end{align}\]
\[ \begin{align} \theta_{i} & = g^{-1}(\vartheta, \theta_{-1}) \\ \theta_{0} & = \vartheta_{b} - \theta_{1} \theta_{b} \end{align}\]
Resultando em:
\[ \left\{\begin{matrix} \vartheta_{b} + \theta_{1} (x - \theta_{b}), \quad & x \leq \theta_{b} \\ \vartheta_{b}, \quad & x > \theta_{b} \end{matrix}\right. \]
xyplot(1:50 ~ 1:50,
pch = 16,
jitter.x = TRUE,
type = c("n", "g"),
scales = list(x = list(draw = FALSE),
y = list(draw = FALSE)),
xlab = "x",
ylab = list(label = "y",
rot = 0),
key = list(space = "right",
title = expression(theta[1]),
cex.title = 1.15,
text = list(c("> 0", "< 0")),
lines = list(col = c(3, 6),
lwd = 2))) +
layer(with(list(vthB = 35,
th1 = 1,
thB = 30),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB),
col = 3,
lwd = 2))) +
layer(with(list(vthB = 15,
th1 = -1,
thB = 30),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB),
col = 6,
lwd = 2)))
Trata-se da junção de dois segmentos de reta. O modelo pode ser escrito da seguinte forma:
\[ \left\{\begin{matrix} \theta_{0} + \theta_{1} x, \quad & x \leq \theta_{b} \\ \theta_{0} + \theta_{1} \theta_{b} + \theta_{2} (x - \theta_{b}), \quad & x > \theta_{b} \end{matrix}\right. \]
A função assume diversas formas dependendo dos valores dos parâmetros.
\[ \vartheta_{b} = \theta_{0} + \theta_{1} \theta_{b} \]
\[ \theta_{0} = \vartheta_{b} - \theta_{1} \theta_{b} \]
Resultando em:
\[ \left\{\begin{matrix} \vartheta_{b} + \theta_{1} (x - \theta_{b}), \quad & x \leq \theta_{b} \\ \vartheta_{b} + \theta_{2} (x - \theta_{b}), \quad & x > \theta_{b} \end{matrix}\right. \]
xyplot(1:50 ~ 1:50,
pch = 16,
jitter.x = TRUE,
type = c("n", "g"),
scales = list(x = list(draw = FALSE),
y = list(draw = FALSE)),
xlab = "x",
ylab = list(label = "y",
rot = 0),
key = list(space = "right",
text = list(c(expression(theta[1]!=0),
expression(theta[1]!=theta[2]),
expression(theta[2]!=0),
expression(theta[2]==0),
expression(theta[1]==0),
expression(theta[1]==theta[2]))),
lines = list(col = c("white", "green3", "white", "magenta", "red", "blue"),
lwd = 2))) +
layer(with(list(vthB = 30,
th1 = 1,
thB = 20,
th2 = 2),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB)+th2*(x-thB),
col = 3,
lwd = 2))) +
layer(with(list(vthB = 30,
th1 = 1,
thB = 20,
th2 = 0),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB)+th2*(x-thB),
col = 6,
lwd = 2))) +
layer(with(list(vthB = 30,
th1 = 0,
thB = 20,
th2 = 4),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB)+th2*(x-thB),
col = 2,
lwd = 2))) +
layer(with(list(vthB = 30,
th1 = 1,
thB = 20,
th2 = 1),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB)+th2*(x-thB),
col = 4,
lwd = 2)))
Trata-se da junção de reta horizontal com o ponto crítico de modelo quadrático. O modelo pode ser escrito da seguinte forma:
\[ \left\{\begin{matrix} \theta_{0} + \theta_{1} x + \theta_{2} x^{2}, \quad & x \leq -\theta_{1} / (2 \theta_{2}) \\ \theta_{0} + \theta_{1} \left( \frac{-\theta_{1}}{2 \theta_{2}} \right) + \theta_{2} \left( \frac{-\theta_{1}}{2 \theta_{2}} \right)^{2}, \quad & x > -\theta_{1} / (2 \theta_{2}) \end{matrix}\right. \]
A junção das funções é no ponto crítico do modelo quadrático e é calculada por \(-\theta_{1}/(2 \theta_{2})\).
\[ \begin{align} \vartheta_{x} & = \frac{-\theta_{1}}{2 \theta_{2}} \\ \vartheta_{y} & = \theta_{0} + \theta_{1} \vartheta_{x} + \theta_{2} \vartheta_{x}^{2} \end{align} \]
\[ \begin{align} \theta_{1} & = -2 \theta_{2} \vartheta_{x} \\ \theta_{0} & = \vartheta_{y} - \theta_{1} \vartheta_{x} - \theta_{2} \vartheta_{x}^{2} \end{align} \]
Resultando em:
\[ \left\{\begin{matrix} \vartheta_{y} + \theta_{2} (x - \vartheta_{x})^{2}, \quad & x \leq \vartheta_{x} \\ \vartheta_{y}, \quad & x > \vartheta_{x} \end{matrix}\right. \]
xyplot(1:50 ~ 1:50,
pch = 16,
jitter.x = TRUE,
type = c("n", "g"),
scales = list(x = list(draw = FALSE),
y = list(draw = FALSE)),
xlab = "x",
ylab = list(label = "y",
rot = 0),
key = list(space = "right",
title = expression(theta[2]),
cex.title = 1.15,
text = list(c("> 0", "< 0")),
lines = list(col = c(3, 6),
lwd = 2))) +
layer(with(list(vthY = 30,
th2 = -1,
vthX = 20),
panel.curve(vthY+{th2*(x-vthX)**2}*(x <= vthX),
col = 3,
lwd = 2))) +
layer(with(list(vthY = 20,
th2 = 1,
vthX = 20),
panel.curve(vthY+{th2*(x-vthX)**2}*(x <= vthX),
col = 6,
lwd = 2)))
Quando temos modelos não lineares mas com parâmetros condicionalmente lineares, podemos fazer o ajuste com a função nls e com o algoritmo plinear.
Um modelo é linear baseado na forma como seus parâmetros aparecem. Um modelo pode ser dito linear se toda derivada parcial
\[ \frac{\partial \mu_{y|\textbf{x}}}{\partial \beta_{i}}, \quad i = 0, 1, 2, ..., p \]
não depender de qualquer parâmetro. (\(\mu_{y|\textbf{x}} = \beta_{0} + \beta_{1} x_{1} + \beta_{2} x_{2}\), por exemplo).
No caso dos modelos não lineares, derivamos parcialmente nossa função \(f(\textbf{x}, \theta)\) em relação a cada parâmetro de \(\theta\).
Suponha um modelo com três parâmetros, \(\theta^{T} = [\theta_{1}\ \theta_{2}\ \theta_{3}]\). Se dado o conhecimento do valor do parâmetro \(\theta_{1}\), a derivada parcial de \(f(\textbf{x}, \theta_{-1})\) em relação a \(\theta_{2}\) não depender de parâmetros, então dizemos que \(\theta_{2}\) é um parâmetro condicionalmente linear.
Mostramos dois exemplos utilizando o mesmo conjunto de dados da seção Um ponto de quebra, o que fazer?
Este modelo já foi mostrado na seção Modelo segmentado conhecido como bi-linear segmentado
. Reparametrizado ele é dado da seguinte forma:
\[ \left\{\begin{matrix} \vartheta_{b} + \theta_{1} (x - \theta_{b}), \quad & x \leq \theta_{b} \\ \vartheta_{b} + \theta_{2} (x - \theta_{b}), \quad & x > \theta_{b} \end{matrix}\right. \]
Temos um modelo linear a partir do momento em que conhecemos o ponto de quebra \(\theta_{b}\). Logo, os demais parâmetros do modelo são condicionalmente lineares. Veja as derivadas parciais:
\[ f(\textbf{x}, \theta) = \vartheta_{b} + {\theta_{1} (x - \theta_{b})} \cdot (x \leq \theta_{b}) + \theta_{2} \cdot (x - \theta_{b}) \]
\[ \begin{matrix} \frac{\partial f(\textbf{x}, \theta)}{\partial \vartheta_{b}} = & 1 \\ \frac{\partial f(\textbf{x}, \theta)}{\partial \theta_{1}} = & x - \theta_{b} \\ \frac{\partial f(\textbf{x}, \theta)}{\partial \theta_{2}} = & x - \theta_{b} \\ \frac{\partial f(\textbf{x}, \theta)}{\partial \theta_{b}} = & \left\{\begin{matrix} -\theta_{1}, \quad & x \leq \theta_{b} \\ -\theta_{2}, \quad & x > \theta_{b} \end{matrix}\right. \end{matrix} \]
Basicamente qualquer modelo não linear pode ser ajustado no R com a função nls.
Com o algoritmo plinear os parâmetros condicionalmente lineares não são ajustados. Para que isso aconteça devem ser passadas para a função as derivadas parciais e chutes iniciais são dados apenas para os demais parâmetros, \(\theta_{b}\) neste caso.
path <- "~/Dropbox/CE092 - 2015-2/Alunos/Henrique/Piecewise Regression/shipment.txt"
data <- read.table(path,
header = TRUE)
m2 <- nls(y ~ cbind(1, {x1-thB}*(x1 <= thB), {x1-thB}*(x1 > thB)),
start = list(thB = 250),
algorithm = "plinear",
data)
coef <- coefficients(m2)
attr(coef, "names")[2:4] <- c("vthB", "th1", "th2")
xyplot(y ~ x1,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
xlab = expression(x[1]),
ylab = list(label = "y",
rot = 0),
scales = list(y = list(tick.number = 8)),
xlim = c(100, 425),
ylim = c(7, 16),
data) +
layer(with(as.list(coef),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB)+th2*(x-thB),
col = 3,
lwd = 2)))
m0 <- lm(y ~ x1,
data)
coef.m0 <- coefficients(m0)
attr(coef.m0, "names") <- c("beta0", "beta1")
data$x2 <- with(data,
ifelse(x1 <= 250, 0, x1 - 250))
m1 <- lm(y ~ x1 + x2,
data)
coef.m1 <- coefficients(m1)
attr(coef.m1, "names") <- c("beta0", "beta1", "beta2")
xyplot(y ~ x1,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
xlab = expression(x[1]),
ylab = list(label = "y",
rot = 0),
scales = list(y = list(tick.number = 8)),
xlim = c(100, 425),
ylim = c(7, 16),
key = list(space = "right",
title = "Regressão",
cex.title = 1.15,
text = list(c("Linear Simples", "Linear por Partes", "Não Linear - plinear")),
lines = list(col = c(3, 4, 2),
lwd = 2)),
data) +
layer(with(as.list(coef.m0),
panel.curve(beta0+beta1*x,
col = 3,
lwd = 2))) +
layer(with(as.list(coef.m1),
panel.curve(beta0+{beta1*x}*(x <= 250)+{beta1*x+beta2*(x-250)}*(x > 250),
col = 4,
lwd = 2))) +
layer(with(as.list(coef),
panel.curve(vthB+{th1*(x-thB)}*(x <= thB)+th2*(x-thB),
col = 2,
lwd = 2)))
O modelo Michaelis-Menten reparametrizado é dado da seguinte forma:
\[ f(\textbf{x}, \theta) = \frac{\theta_{a} x}{\vartheta_{q} \frac{1 - q}{q} + x} \]
Aqui consideramos \(q = 0.5\).
Temos um modelo linear a partir do momento em que conhecemos \(\vartheta_{q}\). Logo, \(\theta_{a}\) é condicionalmente linear. Veja as derivadas parciais:
\[ \frac{\partial f(\textbf{x}, \theta)}{\partial \theta_{a}} = \frac{x}{\vartheta_{q} + x} \qquad \frac{\partial f(\textbf{x}, \theta)}{\partial \vartheta_{q}} = - \frac{\theta_{a} x}{(\vartheta_{q} + x)^{2}} \]
m0 <- nls(y ~ x1/(vthQ+x1),
start = list(vthQ = 110),
algorithm = "plinear",
data)
coef <- coefficients(m0)
attr(coef, "names")[2] <- "thA"
xyplot(y ~ x1,
pch = 16,
jitter.x = TRUE,
type = c("p", "g"),
xlab = expression(x[1]),
ylab = list(label = "y",
rot = 0),
scales = list(y = list(tick.number = 8)),
xlim = c(100, 425),
ylim = c(7, 16),
data) +
layer(with(as.list(coef),
panel.curve(thA*x/(vthQ+x),
col = 3,
lwd = 2)))
Informações traduzidas de: A Tutorial on the Piecewise Regression Approach Applied to Bedload Transport Data
↩
O conjunto de dados analisado está disponível na seguinte página da Penn State University
↩
Este exemplo foi retirado do material Exempo Regressão por Partes
do Professor Gilberto A. Paula↩
Montgomery, D. C.; Peck, E. A. e Vining, G. G. (2001). Introduction to Linear Regression Analysis, Third Edition. Hoboken: Wiley.↩
As informações contidas nesta seção foram retiradas do artigo de Muggeo (2003)↩
Muggeo, V.M.R. (2003) Estimating regression models with unknown break-points. Statistics in Medicine 22, 3055-3071.↩
A. Zeileis, F. Leisch, K. Hornik, and C. Kleiber. \(\text{strucchange}\): An R package for testing for struc-tural change in linear regression models. Journal of Statistical Software, 7(2):1-38, 2002.↩
Informações retiradas da Tese de Doutorado de Walmes Marques Zeviani. Parametrizações interpretáveis em modelos não lineares, 2013
↩
Este não é o único intuito e objetivo das reparametrizações, longe disso. Contudo, neste material, este é o intuito abordado↩