Modelo de Regressão Linear Simples

Nesta aula computacional, aplicaremos os conceitos de modelo de regressão linear simples para explorar a relação entre duas variáveis quantitativas.

1 Dados de resistência ao cisalhamento versus idade do propelente

Inicialmente analisaremos os dados da Tabela 2.1 do livro Montgomery, Peck, e Vining (2021).

Contextualização dos dados:

Motores de foguete são fabricados unindo um propelente de ignição e um propelente sustentador dentro de uma carcaça metálica. Um fator crítico de qualidade nesse processo é a resistência ao cisalhamento da ligação entre os dois tipos de propelente. Suspeita-se que essa resistência ao cisalhamento esteja relacionada à idade do lote de propelente sustentador. Para investigar essa possível relação, foram coletadas 20 observações contendo valores de resistência ao cisalhamento (medida em libras por polegada quadrada - pounds per square inch - psi) e idade do lote correspondente (medida em semanas). Esses dados estão apresentados na tabela a seguir.

Tabela 2.1

**Tabela 2.1**
Observação	Resistência ao Cisalhamento	Idade do Propelente
\(i\)	\(y_i\) (psi)	\(x_i\) (em semanas)
1	2158,70	15,50
2	1678,15	23,75
3	2316,00	8,00
4	2061,30	17,00
5	2207,50	5,50
6	1708,30	19,00
7	1784,70	24,00
8	2575,00	2,50
9	2357,90	7,50
10	2256,70	11,00
11	2165,20	13,00
12	2399,55	3,75
13	1779,80	25,00
14	2336,75	9,75
15	1765,30	22,00
16	2053,50	18,00
17	2414,40	6,00
18	2200,50	12,50
19	2654,20	2,00
20	1753,70	21,50

1.1 Lendo os Dados

Coloque os dados no R.

#Importando os dados
x = c(15.5, 23.75, 8, 17, 5.5, 19, 24, 2.5, 7.5, 11, 13, 3.75, 25, 9.75, 22, 18, 6, 12.5, 2, 21.5)
y = c(2158.7, 1678.15, 2316, 2061.3, 2207.5, 1708.3, 1784.7,
      2575, 2357.9, 2256.7, 2165.2, 2399.55, 1779.8, 2336.75, 1765.3, 2053.5, 2414.4, 2200.5, 2654.2, 1753.7)

1.2 Normalidade

Faça um histograma da variável Resistência ao cisalhamento. Usando o teste do Shapiro-Wilk e do Kolmogorov-Smirnov o que você diria sobre a normalidade para esta variável?

#analisando a densidade da variavel resposta
hist(y, freq=F, main="", cex.lab=1.4, cex.axis=1.4, ylab="densidade", xlab="Resistência ao cisalhamento")

#testando normalidade para a variável resposta
shapiro.test(y)

## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.93181, p-value = 0.1673

ks.test(y, "pnorm", mean(y), sd(y))

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.17719, p-value = 0.5011
## alternative hypothesis: two-sided

1.3 Gráfico de dispersão

Faça um gráfico de dispersão para as variáveis Resistência ao Cisalhamento e Idade do Propelente e calcule o coeficiente de correlação.

#grafico de dispersao
plot(x, y, xlab="Idade do propelente", ylab="Resistência ao cisalhamento", cex.lab=1.4, cex.axis=1.4, bty="n")

#calculando a correlacao linear
cor(x, y)

## [1] -0.9496533

1.4 Ajustando MRLS I

Considere o modelo de regressão linear simples: \[Y_i = \beta_0 + \beta_1 X_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2),\] sendo \(Y_i\) a Resistência ao Cisalhamento da \(i\)-ésima observação e \(X_i\) a Idade do Propelente. Encontre estimativas para \(\beta_0\), \(\beta_1\), \(\sigma^2\), \(\mu_i = \beta_0 + \beta_1 X_i\) e \(e_i\), \(i=1, \ldots, n=20\), usando as fórmulas encontradas na aula anterior.

############################################
#Modelo de Regressao Linear Simples (MRLS)
#Yi = beta0 + beta1 * Xi + ei
#beta0: intercepto (valor do eixo vertical 
#                   quando xi=0)
#beta1: coeficiente angular
#Yi   : variavel resposta
#Xi   : covariavel ou variavel explicativa
#ei   : efeito aleatório
############################################

#####ajustando um Modelo de Regressao Linear Simples

#Estimando os parametros desconhecidos usando as formulas encontradas em sala de aula
beta1chapeu   = sum(y*(x-mean(x))) / 
                sum((x-mean(x))^2)
beta0chapeu   = mean(y) - beta1chapeu * mean(x)
yajustado     = beta0chapeu + beta1chapeu*x
eajustado     = y - yajustado
n             = length(y)
sigma2chapeu  = sum((y-yajustado)^2)/(n-2)

1.5 Ajustando MRLS II

Encontre estimativas para \(\beta_0\), \(\beta_1\), \(\sigma^2\), \(\mu_i = \beta_0 + \beta_1 X_i\) e \(e_i\), \(i=1, \ldots, n=20\), usando a função lm do R.

#Estimando os parametros desconhecidos usando a funcao lm do R
ajuste        = lm(y~x)
names(ajuste)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

ajuste$coefficients

## (Intercept)           x 
##  2627.82236   -37.15359

ajuste$fitted.values

##        1        2        3        4        5        6        7        8 
## 2051.942 1745.425 2330.594 1996.211 2423.478 1921.904 1736.136 2534.938 
##        9       10       11       12       13       14       15       16 
## 2349.170 2219.133 2144.826 2488.496 1698.983 2265.575 1810.443 1959.058 
##       17       18       19       20 
## 2404.901 2163.402 2553.515 1829.020

ajuste$residuals

##           1           2           3           4           5           6 
##  106.758301  -67.274574  -14.593631   65.088687 -215.977609 -213.604131 
##           7           8           9          10          11          12 
##   48.563824   40.061618    8.729573   37.567141   20.374323  -88.946393 
##          13          14          15          16          17          18 
##   80.817415   71.175153  -45.143358   94.442278    9.499187   37.097528 
##          19          20 
##  100.684823  -75.320154

resultados = summary(ajuste)
resultados$sigma^2

## [1] 9236.381

1.6 Comparando ajustes

Compare os resultados encontrados nos itens 4 e 5.

#Comparando os valores encontrados

#para os coeficientes
c(beta0chapeu, ajuste$coefficients[1])

##             (Intercept) 
##    2627.822    2627.822

c(beta1chapeu, ajuste$coefficients[2])

##                   x 
## -37.15359 -37.15359

#para yajustado = ychapeu = beta0chapeu + beta1chapeu * xi
plot(yajustado, ajuste$fitted.values)

summary(round(yajustado - ajuste$fitted.values,8))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0

#para eajustado = echapeu = y - ychapeu
plot(eajustado, ajuste$residuals)

summary(round(eajustado - ajuste$residuals,8))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0

#para sigma2chapeu
resultados = summary(ajuste)
c(sigma2chapeu, resultados$sigma^2)

## [1] 9236.381 9236.381

1.7 Plotando ajuste versus observações

Faça novamente o gráfico de dispersão para as variáveis Resistência ao Cisalhamento e Idade do Propelente e plote em vermelho a reta ajustada.

#grafico de dispersao com a reta ajustada
plot(x, y, xlab="Idade do propelente", ylab="Resistência ao cisalhamento", cex.lab=1.4, cex.axis=1.4, bty="n")
lines(sort(x), beta0chapeu + beta1chapeu*sort(x), col="red", lwd=2)

#usando ggplot2
library(ggplot2)
data = data.frame(Idade_do_propente = x, Resistencia_ao_cisalhamento = y)
data$linha_ajuste <- beta0chapeu + beta1chapeu * data$Idade_do_propente
ggplot(data, aes(x = Idade_do_propente, y = Resistencia_ao_cisalhamento)) +
  geom_point() +  # Adicionar os pontos
  geom_line(aes(y = linha_ajuste), color = "red", linewidth=1) +  # Adicionar a linha de ajuste
  labs(x = "Idade do propelente", y = "Resistência ao cisalhamento") +  # Rótulos dos eixos
  theme(axis.title = element_text(size = 14),  # Tamanho do texto dos rótulos
        axis.text = element_text(size = 14),   # Tamanho do texto dos eixos
        panel.border = element_blank(),         # Remover a borda do gráfico
        panel.background = element_blank())      # Remover o fundo do gráfico

2 Dados Simulados

Considere que o tamanho da amostra é \(n=1000\). Simule uma covariável \(X_i\) da distribuição normal padrão.
Fixe o intercepto em \(\beta_0 = 10\), o coeficiente angular em \(\beta_1=2\) e a variância em \(\sigma^2=2\).
Gere a variável resposta \(Y_1, \ldots, Y_n\) considerando o seguinte modelo: \[Y_i = \beta_0 + \beta_1 X_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2).\]

#Simulando dados
n       = 1000 #tamanho da amostra
beta0   = 10 #intercepto
beta1   = 2 #coeficiente angular
sigma2  = 2 #variancia
x       = rnorm(n, 0, sqrt(1)) #covariavel
e       = rnorm(n, 0, sqrt(sigma2)) #efeito aleatorio
y       = beta0 + beta1 * x + e #variavel resposta

2.1 Análise descritiva

Analise os dados simulados conforme fizemos anteriormente e analise a normalidade dos efeitos aleatórios.

#analisando a densidade da variavel resposta
hist(y, freq=F, main="", cex.lab=1.4, cex.axis=1.4, ylab="densidade", xlab="Y")

#testando normalidade para a variável resposta
shapiro.test(y)

## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.99851, p-value = 0.5613

ks.test(y, "pnorm", mean(y), sd(y))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  y
## D = 0.025401, p-value = 0.5389
## alternative hypothesis: two-sided

#grafico de dispersao
plot(x, y, xlab="X", ylab="Y", cex.lab=1.4, cex.axis=1.4, bty="n")

#calculando a correlacao linear
cor(x, y)

## [1] 0.8286286

#analisando a densidade do efeito aleatório
hist(e, freq=F, main="", cex.lab=1.4, cex.axis=1.4, ylab="densidade", xlab="e")

#testando normalidade para o efeito aleatório
shapiro.test(e)

## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.99857, p-value = 0.6002

ks.test(e, "pnorm", mean(e), sd(e))

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  e
## D = 0.020151, p-value = 0.8115
## alternative hypothesis: two-sided

2.2 Ajuste 1

Suponha que os valores do intercepto, do coeficiente angular e da variância são desconhecidos. Ajuste o modelo de regressão linear simples: \[Y_i = \beta_0 + \beta_1 X_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2).\] Os valores verdadeiros estão próximos dos valores estimados?

#Estimando os parametros desconhecidos
ajuste      = lm(y~x)
resultados  = summary(ajuste)

#Comparando com os valores verdadeiros
c(beta0, ajuste$coefficients[1])

##             (Intercept) 
##   10.000000    9.899355

c(beta1, ajuste$coefficients[2])

##                 x 
## 2.000000 2.020756

c(sigma2, resultados$sigma^2)

## [1] 2.00000 1.94373

#grafico de dispersao com a reta ajustada
plot(x, y, xlab="X", ylab="Y", cex.lab=1.4, cex.axis=1.4, bty="n")
beta0chapeu = ajuste$coefficients[1]
beta1chapeu = ajuste$coefficients[2]
lines(sort(x), beta0chapeu + beta1chapeu*sort(x), col="red", lwd=2)

2.3 Ajuste 2

Gere uma covariável \(W_i \sim Beta(3,2)\), mas não altere a variável resposta (utilize a variável resposta criada acima).

Ajuste o modelo de regressão linear simples: \[Y_i = \beta_0 + \beta_1 W_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2).\]

Os valores verdadeiros estão próximos dos valores estimados?

#simulando uma nova covariavel
w = rbeta(n, 3, 2)

#grafico de dispersao
plot(w, y, xlab="W", ylab="Y", cex.lab=1.4, cex.axis=1.4, bty="n")

#calculando a correlacao linear
cor(w, y)

## [1] -0.004054791

#Estimando os parametros desconhecidos
ajuste = lm(y~w)
resultados = summary(ajuste)

#comparando os valores ajustados com os verdadeiros
#note que a covariavel usada no ajuste nao eh a que
#gerou a variavel resposta e por isso espero que o 
#coeficiente seja zero
c(beta0, ajuste$coefficients[1])

##             (Intercept) 
##   10.000000    9.954161

c(0, ajuste$coefficients[2])

##                       w 
##  0.00000000 -0.05100855

c(sigma2, resultados$sigma^2)

## [1] 2.000000 6.202475

#veja o que acontece com os residuos
mean(resultados$residuals)

## [1] -9.569233e-17

hist(ajuste$residuals, freq=F, cex.lab=1.4, cex.axis=1.4, main="Histograma",ylab="densidade", xlab="resíduos")

#grafico de dispersao com a reta ajustada
plot(w, y, xlab="W", ylab="Y", cex.lab=1.4, cex.axis=1.4, bty="n")
beta0chapeu = ajuste$coefficients[1]
beta1chapeu = ajuste$coefficients[2]
lines(sort(w), beta0chapeu + beta1chapeu*sort(w), col="red", lwd=2)

3 Dados Simulados II

Fixe o tamanho da amostra em \(n=1000\).
Simule 1000 observações de uma covariável, sendo que os 500 primeiros valores devem ser gerados a partir de uma distribuição normal padrão (média 0 e variância 1) e os 500 restantes a partir de uma distribuição normal com média 5 e variância 1.
Fixe o intercepto em \(\beta_0 = 10\), o coeficiente angular em \(\beta_1=2\) e a variância em \(\sigma^2=2\).
Gere a variável resposta \(Y_1, \ldots, Y_n\) considerando o seguinte modelo: \[Y_i = \beta_0 + \beta_1 X_i + e_i, \quad e_i \stackrel{iid}{\sim} N(0,\sigma^2).\]

n       = 1000
beta0   = 10
beta1   = 2
sigma2  = 2
x       = c(rnorm(500), rnorm(500, 5, 1))
e       = rnorm(n, 0, sqrt(sigma2))
y       = beta0 + beta1 * x + e

3.1 Análise descritiva