Um ponto de quebra, o que fazer?
Como exemplo é utilizado um conjunto de dados.
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} \]
- \(y_{i}\) é o custo de manejo para cada carregamento \(i\);
- \(x_{1i}\) é o tamanho de cada carregamento \(i\);
- \(x_{2i}\) é uma variável dummy (\(0\), se \(x_{1i} \leq 250\), e \(1\), se \(x_{1i} > 250\)) de cada carregamento \(i\);
- \(x_{2i}^{*}\) denota \((x_{1i} - 250) x_{2i}\).
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))
Mais de um ponto de quebra, o que fazer?
Exemplo
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.
Mais um exemplo
O conjunto de dados aqui estudo é do livro do Montgomery, Peck e Vining, 2001, Seção 7.2.2.
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.
