ME613 - Análise de Regressão

Modelo de Regressão Polinomial

Introdução

Podemos considerar funções polinomiais como um caso particular do modelo de regressão linear já visto.

Modelo com um preditor - segunda ordem

\[Y=\beta_0+\beta_1X_*+\beta_2X_*^2+\varepsilon\]

em que $X_*=X-\bar{X}$.

Função de resposta quadrática.

$\beta_0$: valor esperado de $Y$ quando $X_*$ é zero, isto é, $X=\bar{X}$.

$\beta_1$: coeficiente de efeito linear.

$\beta_2$: coeficiente de efeito quadrático.

Exemplo

Modelo com um preditor - terceira ordem

\[Y=\beta_0+\beta_1X_*+\beta_2X_*^2+\beta_3X_*^3+\varepsilon\]

em que $X_*=X-\bar{X}$.

Exemplos:

Modelo com dois preditores - segunda ordem

\[Y=\beta_0+\beta_1X_{1*}+\beta_2X_{1*}^2+\beta_3X_{2*}+\beta_4X_{2*}^2+\beta_5X_{1*}X_{2*}+\varepsilon\]

em que $X_{1*}=X_1-\bar{X}_1$ e $X_{2*}=X_2-\bar{X}_2$.

Exemplo

\[E(Y)=1740-4X_{1*}^2-3X_{2*}^2-3X_{1*}X_{2*}\]

Método hierárquico de ajuste de modelo

Pode-se começar com um modelo de segunda ou terceira ordem e ir testando se os coeficientes de ordem maiores são significativos.

Por exemplo:

\[Y=\beta_0+\beta_1X_{*}+\beta_2X_*^2+\beta_3X_*^3+\varepsilon\]

Para testar se $\beta_3=0$ podemos utilizar $SQReg(X_*^3\mid X_*,X_*^2)$. Se quisermos testar se $\beta_2=\beta_3=0$: $SQReg(X_*^2,X_*^3\mid X_*)=SQReg(X_*^2\mid X_*)+SQReg(X_*^3\mid X_*, X_*^2)$

Se um termo de ordem mais alta é mantido no modelo, os de ordem mais baixa devem obrigatoriamente ser mantidos também.

Exemplo

$Y$: número de ciclos

$X_1$: carga, $X_{1*}=(X_1-\bar{X}_1)/0.4$.

$X_2$: temperatura, $X_{2*}=(X_2-\bar{X}_2)/10$.

dados <- read.table("./dados/CH08TA01.txt")
names(dados) <- c("Y","X1","X2")
dados$x1 <- (dados$X1-mean(dados$X1))/0.4 
dados$x2 <- (dados$X2-mean(dados$X2))/10

##      Y  X1 X2 x1 x2
## 1  150 0.6 10 -1 -1
## 2   86 1.0 10  0 -1
## 3   49 1.4 10  1 -1
## 4  288 0.6 20 -1  0
## 5  157 1.0 20  0  0
## 6  131 1.0 20  0  0
## 7  184 1.0 20  0  0
## 8  109 1.4 20  1  0
## 9  279 0.6 30 -1  1
## 10 235 1.0 30  0  1
## 11 224 1.4 30  1  1

Exemplo

Correlação entre $X_1$ e $X_1^2$: 0.99.

Correlação entre $X_{1*}$ e $X_{1*}^2$: 0.

Correlação entre $X_2$ e $X_2^2$: 0.99.

Correlação entre $X_{2*}$ e $X_{2*}^2$: 0.

Exemplo

\[Y=\beta_0+\beta_1X_{1*}+\beta_2X_{2*}+\beta_3X_{1*}^2+\beta_4X_{2*}^2+\beta_5X_{1*}X_{2*}+\varepsilon\]

modelo <- lm(Y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1*x2),data=dados)
summary(modelo)$coef

##              Estimate Std. Error    t value    Pr(>|t|)
## (Intercept) 162.84211   16.60761  9.8052730 0.000187839
## x1          -55.83333   13.21670 -4.2244519 0.008292287
## x2           75.50000   13.21670  5.7124677 0.002297266
## I(x1^2)      27.39474   20.34008  1.3468353 0.235856323
## I(x2^2)     -10.60526   20.34008 -0.5213973 0.624352247
## I(x1 * x2)   11.50000   16.18709  0.7104426 0.509183728

Exemplo

library(alr3)
pureErrorAnova(modelo)

## Analysis of Variance Table
## 
## Response: Y
##              Df Sum Sq Mean Sq F value  Pr(>F)  
## x1            1  18704   18704 26.6315 0.03556 *
## x2            1  34202   34202 48.6970 0.01992 *
## I(x1^2)       1   1646    1646  2.3436 0.26546  
## I(x2^2)       1    285     285  0.4057 0.58935  
## I(x1 * x2)    1    529     529  0.7532 0.47696  
## Residuals     5   5240    1048                  
##  Lack of fit  3   3836    1279  1.8205 0.37378  
##  Pure Error   2   1405     702                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Não rejeitamos $H_0$, isto é, não encontramos evidências para rejeitar que o modelo de segunda ordem é um bom ajuste

Exemplo

Será que um modelo de primeira ordem já seria suficiente?

\[Y=\beta_0+\beta_1X_{1*}+\beta_2X_{2*}+\beta_3X_{1*}^2+\beta_4X_{2*}^2+\beta_5X_{1*}X_{2*}+\varepsilon\]

$H_0$: $\beta_3=\beta_4=\beta_5=0$.

$H_a$: pelo menos um entre $\beta_3, \beta_4$ e $\beta_5$ é diferente de zero.

## Analysis of Variance Table
## 
## Response: Y
##            Df Sum Sq Mean Sq F value   Pr(>F)   
## x1          1  18704   18704 17.8460 0.008292 **
## x2          1  34202   34202 32.6323 0.002297 **
## I(x1^2)     1   1646    1646  1.5704 0.265552   
## I(x2^2)     1    285     285  0.2719 0.624352   
## I(x1 * x2)  1    529     529  0.5047 0.509184   
## Residuals   5   5240    1048                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exemplo

$H_0$: $\beta_q=\beta_{q+1}=\ldots=\beta_{p-1}=0$.
$H_1$: pelo menos um $\beta_q,\ldots,\beta_{p-1}$ não é zero.

(por conveniência, a notação assume que os últimos $p-q$ coeficientes do modelo serão testados)

Estatística do teste:

\[\begin{eqnarray} F^*&=&\frac{SQReg(X_q,\ldots, X_{p-1}\mid X_1,\ldots,X_{q-1})}{p-q}\div\frac{SQE(X_1,\ldots,X_{p-1})}{n-p}\\ &\overset{\mbox{sob $H_0$}}{\sim}&F_{p-q,n-p} \end{eqnarray}\]

Exemplo

$p=6$

$n = 11$

$q=3$

\[F^*=\frac{SQReg(X_{1*}^2,X_{2*}^2,X_{1*}X_{2*}\mid X_{1*},X_{2*})/3}{SQE(X_{1*},X_{2*},X_{1*}^2,X_{2*}^2,X_{1*}X_{2*})/5}\overset{\mbox{sob $H_0$}}{\sim}F_{3,5}\]

\[\begin{eqnarray} SQReg(X_{1*}^2,X_{2*}^2,X_{1*}X_{2*}\mid X_{1*},X_{2*}) &=& SQReg(X_{1*}^2\mid X_{1*},X_{2*})\\ &+& SQReg(X_{2*}^2\mid X_{1*},X_{2*},X_{1*}^2)\\ &+&SQReg(X_{1*}X_{2*}\mid X_{1*},X_{2*},X_{1*}^2,X_{2*}^2)\\ &=&1646 + 284.9 +529\\ &=& 2459.9 \end{eqnarray}\]

\[F_{obs}=\frac{2459.9/3}{1048.1}=0.7823363\]

Comparando com $F(0.95;3,5)=5.41$, não encontramos evidências contra a hipótese nula.

Exemplo

modeloreduz <- lm(Y ~ x1 + x2,data=dados)
anova(modeloreduz,modelo)

## Analysis of Variance Table
## 
## Model 1: Y ~ x1 + x2
## Model 2: Y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1      8 7700.3                           
## 2      5 5240.4  3    2459.9 0.7823 0.5527

Exemplo

Modelo de primeira ordem:

\[Y=\beta_0+\beta_1X_{1*}+\beta_2X_{2*}+\varepsilon\]

modelo1 <- lm(Y ~ x1 + x2,data=dados)
summary(modelo1)$coef

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept) 172.00000   9.354346 18.387175 7.880002e-08
## x1          -55.83333  12.665844 -4.408181 2.261894e-03
## x2           75.50000  12.665844  5.960913 3.378234e-04

Exemplo

Modelo de primeira ordem (variáveis nas escalas originais):

\[Y=\beta_0+\beta_1X_{1}+\beta_2X_{2}+\varepsilon\]

modelo1 <- lm(Y ~ X1 + X2,data=dados)
summary(modelo1)$coef

##              Estimate Std. Error   t value     Pr(>|t|)
## (Intercept)  160.5833  41.615451  3.858743 0.0048174887
## X1          -139.5833  31.664611 -4.408181 0.0022618941
## X2             7.5500   1.266584  5.960913 0.0003378234

Exemplo

Modelo de Regressão com Interação

Efeitos de interação

Um modelo de regressão com $p-1$ variáveis preditoras com efeitos aditivos tem função de regressão da forma:

\[E(Y)=f_1(X_1)+f_2(X_2)+\ldots+f_{p-1}(X_{p-1})\]

em que $f_1,f_2,\ldots,f_{p-1}$ podem ser quaisquer funções.

Por exemplo:

\[E(Y)=\underbrace{\beta_0+\beta_1X_1+\beta_2X_1^2}_{f_1(X_1)}+\underbrace{\beta_3X_2}_{f_2(X_2)}\]

O efeito de $X_1$ e $X_2$ em $Y$ é aditivo.

Efeitos de interação

Já no exemplo a seguir, o efeito não é aditivo, há efeito de interação:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2X_1^2+\beta_3X_2+\beta_3X_1X_2\]

Outro exemplo:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+\beta_4X_1X_2+\beta_5X_1X_3\]

O efeito de uma variável sobre $Y$ irá depender do nível da variável com a qual ela interage.

Interpretação: interação e efeitos lineares

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2+\varepsilon\]

Suponha que $X_1=a$:

\[E(Y)=\beta_0+\beta_1a+\beta_2X_2+\beta_3aX_2\]

Suponha que $X_1=a+1$:

\[E(Y)=\beta_0+\beta_1(a+1)+\beta_2X_2+\beta_3(a+1)X_2\]

Diferença no valor esperado de $Y$ quando aumentamos $X_1$ em 1 unidade:

\[\beta_0+\beta_1(a+1)+\beta_2X_2+\beta_3(a+1)X_2 - (\beta_0+\beta_1a+\beta_2X_2+\beta_3aX_2)\]

\[=\beta_1+\beta_3X_2\]

Interpretação: interação e efeitos lineares

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2+\varepsilon\]

Suponha que $X_2=a$:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2a+\beta_3X_1a\]

Suponha que $X_2=a+1$:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2(a+1)+\beta_3X_1(a+1)\]

Diferença no valor esperado de $Y$ quando aumentamos $X_2$ em 1 unidade:

\[\beta_0+\beta_1X_1+\beta_2(a+1)+\beta_3X_1(a+1) - (\beta_0+\beta_1X_1+\beta_2a+\beta_3X_1a)\]

\[=\beta_2+\beta_3X_1\]

Interpretação: interação e efeitos lineares

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2+\varepsilon\]

Diferença no valor esperado de $Y$ quando aumentamos $X_1$ em 1 unidade:

\[\frac{\partial E(Y)}{\partial X_1}=\beta_1+\beta_3X_2\]

Diferença no valor esperado de $Y$ quando aumentamos $X_2$ em 1 unidade:

\[\frac{\partial E(Y)}{\partial X_2}=\beta_2+\beta_3X_1\]

Interpretação: interação e efeitos lineares

Modelo aditivo:

\[E(Y)=10+2X_1+5X_2\]

$\beta_1$: mudança no valor esperado de $Y$ quando $X_1$ aumenta em 1 unidade, mantendo $X_2$ constante.

Mantendo $X_2$ constante: não importa se $X_2=1$ ou $X_2=3$ o efeito é sempre $\beta_1$ no valor esperado quando $X_1$ aumenta em 1 unidade (retas paralelas).

Interpretação: interação e efeitos lineares

Modelo com interação:

\[E(Y)=10+2X_1+5X_2 +0.5X_1X_2\]

Se $X_2=1$:

\[E(Y)=10+2X_1+5\times1 +0.5X_1\times 1=15+2.5X_1\]

Se $X_2=3$:

\[E(Y)=10+2X_1+5\times3 +0.5X_1\times 3=25+3.5X_1\]

Interpretação: interação e efeitos lineares

Para avaliarmos o efeito de 1 unidade de aumento em $X_1$, devemos considerar o valor de $X_2$ (retas não paralelas).

Interpretação: interação e efeitos curvilinear

Exemplo:

\[E(Y)=65+3X_1+4X_2-10X_1^2-15X_2^2+35X_1X_2\]

Interpretação: interação e efeitos curvilinear

Se $X_1=1$:

\[E(Y)=65+3\times 1+4X_2-10\times(1^2)-15X_2^2+35\times 1\times X_2\]

\[E(Y)=58+39X_2-15X_2^2\]

Se $X_1=-1$:

\[E(Y)=65+3\times (-1)+4X_2-10\times(-1^2)-15X_2^2+35\times (-1)\times X_2\]

\[E(Y)=52-31X_2-15X_2^2\]

Interpretação: interação e efeitos curvilinear

Exemplo

$X_1$: tríceps, $X_{1*}=X_1-\bar{X}_1$.

$X_2$: coxa, $X_{2*}=X_2-\bar{X}_2$.

$X_3$: antebraço, $X_{3*}=X_1-\bar{X}_3$.

$Y$: gordura corporal

##      X1   X2   X3    Y      x1    x2    x3
## 1  19.5 43.1 29.1 11.9  -5.805 -8.07  1.48
## 2  24.7 49.8 28.2 22.8  -0.605 -1.37  0.58
## 3  30.7 51.9 37.0 18.7   5.395  0.73  9.38
## 4  29.8 54.3 31.1 20.1   4.495  3.13  3.48
## 5  19.1 42.2 30.9 12.9  -6.205 -8.97  3.28
## 6  25.6 53.9 23.7 21.7   0.295  2.73 -3.92
## 7  31.4 58.5 27.6 27.1   6.095  7.33 -0.02
## 8  27.9 52.1 30.6 25.4   2.595  0.93  2.98
## 9  22.1 49.9 23.2 21.3  -3.205 -1.27 -4.42
## 10 25.5 53.5 24.8 19.3   0.195  2.33 -2.82
## 11 31.1 56.6 30.0 25.4   5.795  5.43  2.38
## 12 30.4 56.7 28.3 27.2   5.095  5.53  0.68
## 13 18.7 46.5 23.0 11.7  -6.605 -4.67 -4.62
## 14 19.7 44.2 28.6 17.8  -5.605 -6.97  0.98
## 15 14.6 42.7 21.3 12.8 -10.705 -8.47 -6.32
## 16 29.5 54.4 30.1 23.9   4.195  3.23  2.48
## 17 27.7 55.3 25.7 22.6   2.395  4.13 -1.92
## 18 30.2 58.6 24.6 25.4   4.895  7.43 -3.02
## 19 22.7 48.2 27.1 14.8  -2.605 -2.97 -0.52
## 20 25.2 51.0 27.5 21.1  -0.105 -0.17 -0.12

Exemplo

\[E(Y)=\beta_0+\beta_1X_{1*}+\beta_2X_{2*}+\beta_3X_{3*}+\beta_4X_{1*}X_{2*}+\beta_5X_{1*}X_{3*}+\beta_6X_{2*}X_{3*}+\varepsilon\]

modelo <- lm(Y ~ x1 + x2 + x3 + I(x1*x2) + I(x1*x3) + I(x2*x3),data=dat)
summary(modelo)$coef

##                 Estimate Std. Error    t value     Pr(>|t|)
## (Intercept) 20.526893531 1.07362646 19.1192136 6.699796e-11
## x1           3.437808068 3.57866572  0.9606396 3.542612e-01
## x2          -2.094717339 3.03676957 -0.6897848 5.024579e-01
## x3          -1.616337237 1.90721068 -0.8474875 4.120550e-01
## I(x1 * x2)   0.008875562 0.03085046  0.2876963 7.781144e-01
## I(x1 * x3)  -0.084790836 0.07341774 -1.1549093 2.689155e-01
## I(x2 * x3)   0.090415385 0.09200130  0.9827621 3.436619e-01

Exemplo

anova(modelo)

## Analysis of Variance Table
## 
## Response: Y
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## x1          1 352.27  352.27 52.2238 6.682e-06 ***
## x2          1  33.17   33.17  4.9173   0.04503 *  
## x3          1  11.55   11.55  1.7117   0.21343    
## I(x1 * x2)  1   1.50    1.50  0.2217   0.64552    
## I(x1 * x3)  1   2.70    2.70  0.4009   0.53760    
## I(x2 * x3)  1   6.51    6.51  0.9658   0.34366    
## Residuals  13  87.69    6.75                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exemplo

$H_0$: $\beta_4=\beta_5=\beta_6=0$

$H_1$: pelo menos um dentre $\beta_4,\beta_5,\beta_6$ é diferente de 0.

$p=7$

$n = 20$

$q=4$

\[F^*=\frac{SQReg(X_{1*}X_{2*},X_{1*}X_{3*},X_{2*}X_{3*}\mid X_{1*},X_{2*},X_{3*})/3}{SQE(X_{1*},X_{2*},X_{3*},X_{1*}X_{2*},X_{1*}X_{3*},X_{2*}X_{3*})/13}\overset{\mbox{sob $H_0$}}{\sim}F_{3,13}\]

Exemplo

\[\begin{eqnarray} SQReg(X_{1*}X_{2*},X_{1*}X_{3*},X_{2*}X_{3*}\mid X_{1*},X_{2*},X_{3*}) &=& SQReg(X_{1*}X_{2*}\mid X_{1*},X_{2*},X_{3*})\\ &+& SQReg(X_{1*}X_{3*}\mid X_{1*},X_{2*},X_{3*},X_{1*}X_{2*})\\ &+&SQReg(X_{2*}X_{3*}\mid X_{1*},X_{2*},X_{3*},X_{1*}X_{2*},X_{1*}X_{3*})\\ &=&1.5 + 2.7 +6.514836\\ &=& 10.714836 \end{eqnarray}\]

\[F_{obs}=\frac{10.714836/3}{6.7}=0.5330764\]

Comparando com $F(0.95;3,13)=3.41$, não encontramos evidências contra a hipótese nula.

Exemplo

modeloreduz <- lm(Y ~ x1 + x2 + x3,data=dat)
anova(modeloreduz,modelo)

## Analysis of Variance Table
## 
## Model 1: Y ~ x1 + x2 + x3
## Model 2: Y ~ x1 + x2 + x3 + I(x1 * x2) + I(x1 * x3) + I(x2 * x3)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     16 98.405                           
## 2     13 87.690  3    10.715 0.5295 0.6699

Preditores Qualitativos

Exemplo: Seguros

$Y$ = meses até a implementação

$X_1$ = tamanho da firma (em milhões de dólares)

\[\begin{eqnarray} X_2 = \begin{cases} 1, & \mbox{se a firma tem açoes na bolsa} \\ 0, & \mbox{caso contrário} \end{cases} \end{eqnarray}\]

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\varepsilon\]

\[E(Y)=\beta_0+\beta_1X_1+\beta_2X_2\]

Exemplo: Seguros

\[E(Y)=\beta_0+\beta_1X_1+\beta_2X_2\]

Se a firma não tem ações na bolsa, então $X_2=0$:

\[E(Y)=\beta_0+\beta_1X_1\]

Se a firma tem ações na bolsa, então $X_2=1$:

\[E(Y)=(\beta_0+\beta_2)+\beta_1X_1\]

Exemplo: Seguros

##               Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept) 33.8740690 1.813858297  18.675146 9.145269e-13
## X1          -0.1017421 0.008891218 -11.442990 2.074687e-09
## X2           8.0554692 1.459105700   5.520826 3.741874e-05

Exemplo: Seguros

Incluindo termo de interação:

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2+\varepsilon\]

\[E(Y)=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2\]

Se a firma não tem ações na bolsa, então $X_2=0$:

\[E(Y)=\beta_0+\beta_1X_1\]

Se a firma tem ações na bolsa, então $X_2=1$:

\[E(Y)=(\beta_0+\beta_2)+(\beta_1+\beta_3)X_1\]

Exemplo: Seguros

##                  Estimate Std. Error     t value     Pr(>|t|)
## (Intercept) 33.8383694765 2.44064985 13.86449166 2.472768e-10
## X1          -0.1015306249 0.01305254 -7.77861250 7.965637e-07
## X2           8.1312501223 3.65405169  2.22526959 4.079375e-02
## I(X1 * X2)  -0.0004171412 0.01833121 -0.02275578 9.821265e-01

Exemplo: Seguros

Variável preditora com mais de duas classes

Exemplo: Desgaste ($Y$), velocidade ($X_1$) e modelo de uma peça.

Existem 4 tipos de modelos: M1, M2, M3 e M4.

Definimos 3 variáveis "dummy":

\[\begin{eqnarray} X_2 = \begin{cases} 1, & \mbox{se M1} \\ 0, & \mbox{caso contrário} \end{cases} \end{eqnarray}\]

\[\begin{eqnarray} X_3 = \begin{cases} 1, & \mbox{se M2} \\ 0, & \mbox{caso contrário} \end{cases} \end{eqnarray}\]

\[\begin{eqnarray} X_4 = \begin{cases} 1, & \mbox{se M3} \\ 0, & \mbox{caso contrário} \end{cases} \end{eqnarray}\]

Variável preditora com mais de duas classes

\[E(Y)=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+\beta_4X_4\]

Se a peça é do tipo M4:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2\times0+\beta_3\times0+\beta_4\times0=\beta_0+\beta_1X_1\]

Se a peça é do tipo M1:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2\times1+\beta_3\times0+\beta_4\times0=(\beta_0+\beta_2)+\beta_1X_1\]

Se a peça é do tipo M2:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2\times0+\beta_3\times1+\beta_4\times0=(\beta_0+\beta_3)+\beta_1X_1\]

Se a peça é do tipo M3:

\[E(Y)=\beta_0+\beta_1X_1+\beta_2\times0+\beta_3\times0+\beta_4\times1=(\beta_0+\beta_4)+\beta_1X_1\]

Variável preditora com mais de duas classes

O modelo de primeira ordem implica no fato de que o efeito da velocidade é linear e com o mesmo coeficiente angular para todos os modelos de peça. Temos diferentes interceptos para cada modelo.

Variável preditora com mais de duas classes

$\beta_1$: mudança esperada no desgaste da peça ($Y$) para cada unidade de aumento na velocidade ($X_1$), considerando mesmo modelo de peça.
$\beta_2$: diferença esperada do desgaste da peça entre modelos M1 e M4, considerando a mesma velocidade.
$\beta_3$: diferença esperada do desgaste da peça entre modelos M2 e M4, considerando a mesma velocidade.
$\beta_4$: diferença esperada do desgaste da peça entre modelos M3 e M4, considerando a mesma velocidade.

Variável preditora com mais de duas classes

Qual a diferença esperada do desgaste da peça entre modelos M3 e M2, mantendo a mesma velocidade?

Para modelo M3:

\[E(Y)=(\beta_0+\beta_4)+\beta_1X_1\]

Para modelo M2:

\[E(Y)=(\beta_0+\beta_3)+\beta_1X_1\]

A diferença entre M3 e M2, mantendo a mesma velocidade:

\[(\beta_0+\beta_4)+\beta_1X_1 - [(\beta_0+\beta_3)+\beta_1X_1] = \beta_4-\beta_3\]

Variável preditora com mais de duas classes

Após obtermos estimativas: $\hat{\beta}_4-\hat{\beta}_3$ e devemos também fornecer o erro-padrão da estimativa.

Lembre que:

\[Var(\hat{\beta}_4-\hat{\beta}_3)=Var(\hat{\beta}_4)+Var(\hat{\beta}_3)-2Cov(\hat{\beta}_4,\hat{\beta}_3)\]

Exemplo: Fábrica de sabão

$Y$: resíduo de sabão

$X_1$: velocidade

\[\begin{eqnarray} X_2 = \begin{cases} 1, & \mbox{se produção na linha 1} \\ 0, & \mbox{caso contrário} \end{cases} \end{eqnarray}\]

##      Y  X1 X2
## 1  218 100  1
## 2  248 125  1
## 3  360 220  1
## 4  351 205  1
## 5  470 300  1
## 6  394 255  1
## 7  332 225  1
## 8  321 175  1
## 9  410 270  1
## 10 260 170  1
## 11 241 155  1
## 12 331 190  1
## 13 275 140  1
## 14 425 290  1
## 15 367 265  1
## 16 140 105  0
## 17 277 215  0
## 18 384 270  0
## 19 341 255  0
## 20 215 175  0
## 21 180 135  0
## 22 260 200  0
## 23 361 275  0
## 24 252 155  0
## 25 422 320  0
## 26 273 190  0
## 27 410 295  0

Exemplo: Fábrica de sabão

Iremos ajustar um modelo assumindo que:

a relação entre a quantidade de resíduo e velocidade é linear para as duas linhas de produção;
retas diferentes para as duas linhas de produção;
as variâncias dos termos de erros ao redor de cada reta são iguais.

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_1X_2+\varepsilon\]

Para a linha 1: $E(Y)=(\beta_0+\beta_2)+(\beta_1+\beta_3)X_1$.

Para a linha 2: $E(Y)=\beta_0+\beta_1X_1$.

Exemplo: Fábrica de sabão

modelo <- lm(Y ~ X1 + X2 + I(X1*X2), data=dados)
summary(modelo)$coef

##               Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)  7.5744646 20.8696979  0.3629408 7.199633e-01
## X1           1.3220488  0.0926247 14.2731771 6.446165e-13
## X2          90.3908632 28.3457320  3.1888703 4.085851e-03
## I(X1 * X2)  -0.1766614  0.1288377 -1.3711932 1.835463e-01

anova(modelo)

## Analysis of Variance Table
## 
## Response: Y
##            Df Sum Sq Mean Sq  F value    Pr(>F)    
## X1          1 149661  149661 347.5548 2.224e-15 ***
## X2          1  18694   18694  43.4129 1.009e-06 ***
## I(X1 * X2)  1    810     810   1.8802    0.1835    
## Residuals  23   9904     431                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exemplo: Fábrica de sabão

Se quisermos testar a hipótese nula de que temos apenas uma reta para representar as duas linhas:

$H_0$: $\beta_2=\beta_3=0$

$H_a$: pelo menos um entre $\beta_2$ e $\beta_3$ é diferente de zero.

Estatística do teste:

\[\begin{eqnarray} F^*&=&\frac{SQReg(X_q,\ldots, X_{p-1}\mid X_1,\ldots,X_{q-1})}{p-q}\div\frac{SQE(X_1,\ldots,X_{p-1})}{n-p}\\ &\overset{\mbox{sob $H_0$}}{\sim}&F_{p-q,n-p} \end{eqnarray}\]

Exemplo: Fábrica de sabão

$p=4$

$n = 27$

$q=2$

\[F^*=\frac{SQReg(X_2,X_1X_2\mid X_1)/2}{SQE(X_1,X_2,X_1X_2)/23}\overset{\mbox{sob $H_0$}}{\sim}F_{2,23}\]

\[\begin{eqnarray} SQReg(X_2,X_1X_2\mid X_1) &=& SQReg(X_2\mid X_1) + SQReg(X_1X_2\mid X_1,X_2\\ &=&1.86941\times 10^{4} + 809.6\\ &=& 1.95037\times 10^{4} \end{eqnarray}\]

\[F_{obs}=\frac{1.95037\times 10^{4}/2}{430.6}=22.6471203\]

Comparando com $F(0.95;2,23)=3.42$, encontramos evidências contra a hipótese nula.

Exemplo: Fábrica de sabão

modeloreduz <- lm(Y ~ X1, data=dados)
anova(modeloreduz,modelo)

## Analysis of Variance Table
## 
## Model 1: Y ~ X1
## Model 2: Y ~ X1 + X2 + I(X1 * X2)
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     25 29407.8                                  
## 2     23  9904.1  2     19504 22.646 3.669e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exemplo: Fábrica de sabão

Se quisermos testar a hipótese nula de que para as duas linhas de produção o coeficiente angular é o mesmo:

$H_0$: $\beta_3=0$

$H_a$: $\beta_3\neq0$.

$p=4$

$n = 27$

$q=3$

\[F^*=\frac{SQReg(X_1X_2\mid X_1, X_2)/1}{SQE(X_1,X_2,X_1X_2)/23}\overset{\mbox{sob $H_0$}}{\sim}F_{1,23}\]

Exemplo: Fábrica de sabão

\[F_{obs}=\frac{809.6/1}{430.6}=1.8801672\]

Comparando com $F(0.95;1,23)=4.28$, não encontramos evidências contra a hipótese nula.

Leitura

Applied Linear Statistical Models: Seções 8.1-8.3, 8.5-8.7.
Faraway - Linear Models with R: Capítulo 14.
Draper & Smith - Applied Regression Analysis: Capítulo 12.