Interessa estudar a relação entre a resistência de um determinado tipo de plástico (Y) e o tempo que decorre a partir da conclusão do processo de moldagem até ao momento de medição da resistência (x em horas). As observações que se seguem foram efectuadas em 12 peças construídas com este plástico, escolhidas aleatoriamente.

  1. Represente graficamente as observações e desenhe a recta que, no seu entender, melhor se ajusta às observações.
  2. Considere um modelo de regressão linear simples para explicar as observações. Obtenha a estimativa dos mínimos quadrados dos coeficientes da recta de regressão e desenhe-a no gráfico.
  3. Calcule o coeficiente de determinação e comente o valor obtido.
  4. Proceda ao teste da hipótese ``O coeficiente angular é nulo’’. Qual o interesse desta hipótese? Relacione-o com o resultado obtido em c).
  5. Calcule o intervalo de confiança a 95% para o valor esperado da resistência obtida 48 horas depois de concluída a moldagem. Acha legítimo usar o mesmo procedimento tratando-se de um período de 10 horas em vez de 48 horas? Justifique a sua resposta.

Leitura dos dados

rm(list=ls())
x <- c(32,48,72,64,48,16,40,48,48,24,80,56) 
y <- c(230,262,323,298,255,199,248,279,267,214,359,305)

Figura 1

Valores observados da resistência do tipo de plástico (Y) versus o tempo de moldagem até ao momento da sua medição (x).

plot(x,y,col="red")

Modelo de regressão linear simples

\[ Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \ i\!=\!1,\ldots,n \]

onde \(\varepsilon_i, i\!=\!1,\ldots,n, \ \sim N(0,\sigma^2)\) independentes.

Estimativas de mínimos quadrados e de máxima verosimilhança: \[ \begin{array}{l} {\hat{\beta}_1} = \frac{\sum\limits_{i = 1}^{12} {x_i y_i } - n\bar{x}\bar{y}}{\sum\limits_{i = 1}^{12} {x_i^2 } - n\bar{x}^2} = \frac{164752 - 12 \times 48 \times 269.92}{31486 - 12 \times (48)^2} = 2.4167 \\ {\hat {\beta }_0} = \bar{y} - {\hat {\beta }_1} \bar{x} = 269.92 - 2.4167 \times 48 = 153.9167 \end{array} \]

Equação de regressão estimada: \[ {\hat y}_i \equiv {\widehat E(Y_i)} = {\hat \beta_0} + {\hat{\beta}_1} x_i = 153.91 + 2.4167 x_i \]

Resultados 1 com comando lm (linear model):

mrl <- lm(y~x) 
mrl
#> 
#> Call:
#> lm(formula = y ~ x)
#> 
#> Coefficients:
#> (Intercept)            x  
#>     153.917        2.417

Figura 2

Valores ajustados da resistência do tipo de plástico (Y) a partir do tempo de moldagem até ao momento da sua medição (x).

plot(x,y)
lines(x,fitted(mrl),type="l",lty=1,col="blue")

Avaliação do modelo (coeficiente de determinação)

\[ r^2 = \frac{\left( {\sum_{i = 1}^n {x_i y_i } - n\bar{x}\bar{y}} \right)^2}{\left( {\sum_{i = 1}^n {x_i^2 } - n\bar {x}^2} \right)\times \left( {\sum_{i = 1}^n {y_i^2 } - n\bar{y}^2} \right)} = \frac{(9278.08)^2}{3838 \times 23378.92} = 0.9593 \] Isto é, 95.93% da variabilidade total da resistência do plástico é explicada pelo modelo de regressão com o tempo decorrido entre a moldagem e a medição da resistência.

Resultados 2 com comando lm (linear model):

mrl <- lm(y~x) 
summary(mrl)$r.squared 
#> [1] 0.9592689

Teste de hipóteses (hipótese de interesse: coeficiente angular nulo)

Teste de hipóteses:

  1. Hipóteses: \(\ H_0: \beta_1 = 0 \ (\equiv \beta_1^0)\ \) versus \(\ H_1: \beta_1 \not= 0\).
  2. Estatística do teste: \[ T_0 = \frac{\hat {\beta}_1 - \beta_1^0}{\sqrt {\frac{\hat {\sigma }^2}{(\sum_i {x_i^2 } - 10 \times \bar{x}^2)}} } \ \overset{H_0}{\sim} \ t_{\left( {10} \right)} , \] cujo valor observado é dado por \(t_0 = (2.4167)/\sqrt{(9.758)^2/(3838)} = 15.35\).
  3. Valor-P: \(P = 2 \times P(T_0>15.35) = 2.81 \times 10^{-8}\). Note-se que \(P < 0.001 = 2 \times 0.0005\ \) pois \(\ F_{t(10)}^{-1}(0.9995) = 4.587\).
  4. Conclusão: Rejeita-se \(H_0\) para níveis de significância de pelo menos \(2.81 \times 10^{-8}\). Portanto, há evidência contra \(H_0\), nomeadamente, para \(\alpha=0.01\), i.e., o tempo decorrido entre a moldagem e a medição da resistência influencia significativamente a resistência do plástico.

Resultados 3 com comando lm (linear model):

mrl <- lm(y~x) 
summary(mrl)
#> 
#> Call:
#> lm(formula = y ~ x)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -14.917  -5.667  -1.917   7.083  15.750 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 153.9167     8.0666   19.08 3.40e-09 ***
#> x             2.4167     0.1575   15.35 2.81e-08 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 9.758 on 10 degrees of freedom
#> Multiple R-squared:  0.9593, Adjusted R-squared:  0.9552 
#> F-statistic: 235.5 on 1 and 10 DF,  p-value: 2.807e-08

Intervalo de confiança a 95% para \(E(Y|x_0=48)\)

Intervalo de confiança a 95% \(E(Y|x_0)\):

  1. Variável fulcral para \(E(Y|x_0)=\beta_0 + \beta_1 x_0\): \[ W = \frac{\left( {\hat {\beta }_0 + \hat {\beta }_1 x_0 } \right) - \left( {\beta _0 + \beta _1 x_0 } \right)}{\sqrt {\left( {\frac{1}{n} + \frac{\left( {\bar {x} - x_0 } \right)^2}{(\sum_i {x_i^2 } - n \bar{x}^2)}} \right)\hat {\sigma }^2} } \ \sim \ t_{(10)} \]

  2. \(P(-2.228 < W < 2.228 ) = 0.95\)  pois  \(F_{t(10)}^{-1}(0.975) = 2.228\)

  3. Intervalo aleatório de confiança para \(E(Y|x_0)\) a \(95\%\): \[ (\hat{\beta}_0 + \hat{\beta}_1 x_0) \pm 2.228 \sqrt{\left({\frac{1}{n} + \frac{\left( {\bar{x}-x_0} \right)^2}{\sum_i {x_i^2} - n\bar{x}^2} } \right) \hat{\sigma}^2} \]

  4. Estimativa pontual: \[ \hat{y}_0 \equiv {\widehat E(Y|x_0=48)} = 153.91 + 2.4167 \times 48 = 269.91 \]

  5. Intervalo de confiança para \(E(Y|x_0=48)\) a \(95\%\): \[ \begin{array}{l} 269.91 \pm 2.228 \times 9.758 \sqrt{\frac{1}{10} + \frac{(48-48)^2}{3838}} \\ 269.91 \pm 6.875 \\ (263.035; 276.785) \end{array} \]

Resultados 4 com comando lm (linear model):

mrl <- lm(y~x) 
newdata1 = data.frame(x=48) 
predict(mrl, newdata1, interval="confidence") 
#>        fit    lwr      upr
#> 1 269.9167 263.64 276.1933

Intervalo de confiança a 95% para \(E(Y|x_0=10)\)

Note-se que não é aconselhado usarmos \(x_0\) fora do domínio de actuação dos dados observados, visto que não há informação fora desse domínio. O que acontece com \(x_0=10\notin (16,80)\).

Resultados 5 com comando lm (linear model):

mrl <- lm(y~x) 
newdata2 = data.frame(x=10) 
predict(mrl, newdata2, interval="confidence") 
#>        fit      lwr      upr
#> 1 178.0833 163.3466 192.8201

Resíduos padronizados

\[ r_i^S = (y_i-\hat{y}_i)/\sqrt{\hat{\sigma}^2} \]

Resultados 6 com comando lm (linear model):

mrl <- lm(y~x) 
stdres = rstandard(mrl) 
plot(x, stdres, ylab="Resíduos", xlab="tempo de moldagem", main="Resistência do plástico", col="red")


Referências:

Silva, G.L. e Paulino, C.D. (2023). Notas de Probabilidades e Estatística, 2ª edição. Instituto Superior Técnico, Lisboa. https://www.math.tecnico.ulisboa.pt/~gsilva/PE_slides-print.pdf

R Core Team, R: A language and environment for statistical computing, R Foundation for Statistical Computing, 2019, Vienna, Austria. https://www.R-project.org/