Comandos R para análises estatísticas

Capítulo 16: Regressão Linear Simples

Figura 16.1

Anteriormente consideramos idade como um fator. Vamos agora considerá-lo como uma variável numérica e ajustar o modelo linear para a variável tempo de reação ao estímulo (Y).

attach(tab15_1)
n_idade_X<-as.numeric(as.character(idade_X))
plot(n_idade_X,tempo_Y, pch=20, xlab="Idade(x)", ylab="Estímulo(y)", col="darkblue")
abline(lm(tempo_Y~n_idade_X), lwd=2, col="red")

Figura 16.1: Gráfico de dispersão de idade e reação ao estímulo, com reta ajustada.

Exemplo 16.1

Para ajustar o modelo linear entre idade(X) e tempo de reação do i-ésimo indivíduo(Y), utilizaremos o comando lm:.

fit16_1<-lm(tempo_Y~n_idade_X)
fit16_1
## 
## Call:
## lm(formula = tempo_Y ~ n_idade_X)
## 
## Coefficients:
## (Intercept)    n_idade_X  
##        80.5          0.9

Para calcular o valor de \(\hat{y}(20)\), podemos utilizar a função predict ou estimar o valor de acordo com as estimativas dos parâmetros:

# fit16_1$coefficients nos dá os coeficientes ajustados, sendo fit16_1$coefficients[1] o intercepto e fit16_1$coefficients [2] a inclinação.

y_20=fit16_1$coefficients[1]+fit16_1$coefficients[2]*20
y_20
## (Intercept) 
##        98.5
predict(object = fit16_1)
##     1     2     3     4     5     6     7     8     9    10    11    12 
##  98.5  98.5  98.5  98.5 103.0 103.0 103.0 103.0 107.5 107.5 107.5 107.5 
##    13    14    15    16    17    18    19    20 
## 112.0 112.0 112.0 112.0 116.5 116.5 116.5 116.5
y_33=fit16_1$coefficients[1]+fit16_1$coefficients[2]*33
y_33
## (Intercept) 
##       110.2

Note que os primeiros 4 valores preditos, que correspondem às observações em que \(x=20\), coincidem com o valor de \(\hat{y}(20)\) calculado manualmente.

Na função predict é possivel também calcular os valores preditos para valores específicos das variáveis explicativas:

newdata=data.frame(n_idade_X=20) # `data.frame` contém os dados das variáveis explicativas utilizados para calcular as previsões.
predict(fit16_1, newdata)
##    1 
## 98.5

Tabela 16.1

Para calcular o resíduo de um modelo linear, utilize a função resid:

tab16_1<-cbind(tab15_1[,1:3],resid_fit16_1=resid(fit16_1))
print(kable(tab16_1,caption="**Tabela 16.1**: Resíduos para o modelo (16.18).", col.names = c("Tempo de Reação", "Sexo", "Idade", "Resíduos")))
Tabela 16.1: Resíduos para o modelo (16.18).
Tempo de Reação Sexo Idade Resíduos
96 H 20 -2.5
92 M 20 -6.5
106 H 20 7.5
100 M 20 1.5
98 M 25 -5.0
104 H 25 1.0
110 H 25 7.0
101 M 25 -2.0
116 M 30 8.5
106 H 30 -1.5
109 H 30 1.5
100 M 30 -7.5
112 M 35 0.0
105 M 35 -7.0
118 H 35 6.0
108 H 35 -4.0
113 M 40 -3.5
112 M 40 -4.5
127 H 40 10.5
117 H 40 0.5

As estatísticas que aparecem nas últimas 4 linhas da tabela são calculadas abaixo:

anova_fit16_1<-anova(fit16_1)
SQRes=sum(tab16_1$resid_fit16_1^2)
S2e=SQRes/(20-2)

SQRes
## [1] 563
S2e
## [1] 31.278
sqrt(S2e)
## [1] 5.5927
2*sqrt(S2e)
## [1] 11.185

Exemplo 16.2

s2<-var(tempo_Y)
s2
## [1] 72.263

Exemplo 16.3

O \(R^2\) é calculado pela função summary aplicada ao ajuste do modelo:

summary(fit16_1)$r.squared 
## [1] 0.58995

Tabela 16.3

A tabela ANOVA é a mesma calculada no capítulo anterior:

anova(fit16_1)
## Analysis of Variance Table
## 
## Response: tempo_Y
##           Df Sum Sq Mean Sq F value  Pr(>F)
## n_idade_X  1    810     810    25.9 7.7e-05
## Residuals 18    563      31

Exemplo 16.4

Para a regressão linear, os intervalos de confiança para os coeficientes são facilmente calculados com a função confint:

confint(fit16_1)
##                2.5 %  97.5 %
## (Intercept) 69.04778 91.9522
## n_idade_X    0.52844  1.2716

Exemplo 16.5

A função summary apresenta detalhes do modelo ajustado, inclusive as estatísticas t para cada coeficiente e seus respectivos valores-p:

summary(fit16_1)
## 
## Call:
## lm(formula = tempo_Y ~ n_idade_X)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -7.50  -4.12  -0.75   2.62  10.50 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   80.500      5.451   14.77  1.7e-11
## n_idade_X      0.900      0.177    5.09  7.7e-05
## 
## Residual standard error: 5.59 on 18 degrees of freedom
## Multiple R-squared:  0.59,   Adjusted R-squared:  0.567 
## F-statistic: 25.9 on 1 and 18 DF,  p-value: 7.66e-05

Exemplo 16.6

A função predict que utilizamos anteriormente também calcula o intervalo de previsão utilizando a opção interval="predict":

newdata=data.frame(n_idade_X=28) # `data.frame` contém os dados das variáveis explicativas utilizados para calcular as previsões.
predict(fit16_1, newdata, interval="predict") 
##     fit    lwr    upr
## 1 105.7 93.637 117.76

Exemplo 16.7

attach(tab16_1)
nui=1/length(n_idade_X)+(n_idade_X-mean(n_idade_X))^2/sum((n_idade_X-mean(n_idade_X))^2)
ri=resid_fit16_1/(sqrt(S2e)*sqrt(1-nui))
tab16_4<-cbind(tab16_1[3:4],zi=resid_fit16_1/sqrt(S2e),ri)
Tabela 16.4: Resíduos para o modelo (16.18)
Idade ei zi ri
20 -2.5 -0.447 -0.485
20 -6.5 -1.162 -1.261
20 7.5 1.341 1.455
20 1.5 0.268 0.291
25 -5.0 -0.894 -0.930
25 1.0 0.179 0.186
25 7.0 1.252 1.301
25 -2.0 -0.358 -0.372
30 8.5 1.520 1.559
30 -1.5 -0.268 -0.275
30 1.5 0.268 0.275
30 -7.5 -1.341 -1.376
35 0.0 0.000 0.000
35 -7.0 -1.252 -1.301
35 6.0 1.073 1.115
35 -4.0 -0.715 -0.744
40 -3.5 -0.626 -0.679
40 -4.5 -0.805 -0.873
40 10.5 1.877 2.036
40 0.5 0.089 0.097 
par(mfrow=c(1,2))
plot(n_idade_X,resid_fit16_1, cex=2, col="darkblue", pch=20, xlab="Idade", ylab="Resíduos",main="(a)")
abline(h=0)
plot(n_idade_X,tab16_4$zi, cex=2, col="darkblue", pch=20, xlab="Idade", ylab="Resíduos padronizados",main="(b)")
abline(h=0)

Figura 16.4: Resíduos para o Exemplo 16.1. (a) \(\hat{e}_i=y_i-\hat{y}_i\); (b) resíduos padronizados.

hist(resid_fit16_1, xlab="Resíduos", ylab="", col="lightblue3", border="white", main="")

Figura 16.8: Histograma dos resíduos do modelo (16.18).

qqnorm(resid_fit16_1, cex=1.5, col="darkblue", xlab="Quantis da normal padrão", ylab="Quantis dos resíduos", pch=20)
qqline(resid_fit16_1, col="darkred")

Figura 16.9: Gráfico q x q (normalidade) para os resíduos do modelo (16.18).

Como antes, é possivel plotar os gráficos de resíduos a partir do comando plot:

par(mfrow=c(2,2))
plot(fit16_1)

Exemplo 16.8

tab16_5<-data.frame(n=1:25,
      xi=c(35.3,29.7,30.8,58.8,61.4,71.3,74.4,76.7,70.7,57.5,46.4,28.9,28.1,39.1,46.8,48.5,59.3,70,70,74.5,72.1,58.1,44.6,33.4,28.6),
      yi=c(10.98,11.13,12.51,8.4,9.27,8.73,6.36,8.5,7.82,9.14,8.24,12.19,11.88,9.57,10.94,9.58,10.09,8.11,6.83,8.88,7.68,8.47,8.86,10.36,11.08))

Tabela 16.5: Temperatura e quantidade de vapor de um processo industrial.

Ajuste do modelo para a Tabela 16.5

fit16_5<-lm(tab16_5$yi~tab16_5$xi)
fit16_5
## 
## Call:
## lm(formula = tab16_5$yi ~ tab16_5$xi)
## 
## Coefficients:
## (Intercept)   tab16_5$xi  
##     13.6230      -0.0798

Figura 16.10

attach(tab16_5)
plot(xi,yi, pch=20, xlab="Temperatura", ylab="Vapor", col="darkblue",main="(a)")
abline(lm(yi~xi), lwd=2, col="red")

Figura 16.10(a): gráfico de dispersão com reta ajustada;

par(mfrow=c(1,2))
plot(xi,resid(fit16_5), cex=2, col="darkblue", pch=20, xlab="Temperatura", ylab="Resíduos",main="(b)")
abline(h=0)

Figura 16.10(b): resíduos vs temperatura;

qqnorm(resid(fit16_5), cex=1.5, col="darkblue", xlab="Quantis da normal padrão", ylab="Quantis dos resíduos", pch=20,main="(c)")
qqline(resid(fit16_5), col="darkred")

Figura 16.10(c): gráfico q x q (normalidade).

Exemplo 16.9

tab16_9A<-data.frame(x=c(119,155,174,190,196,233,272,253,276),
                     y=c(112,152,172,183,192,228,263,239,263))

Ajuste do modelo para os dados do Exemplo 16.9 sem o intercepto é ajustado adicionando -1 à equação:

fit16_9<-lm(tab16_9A$y~tab16_9A$x-1)
fit16_9
## 
## Call:
## lm(formula = tab16_9A$y ~ tab16_9A$x - 1)
## 
## Coefficients:
## tab16_9A$x  
##      0.965
summary(fit16_9)
## 
## Call:
## lm(formula = tab16_9A$y ~ tab16_9A$x - 1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.082 -2.805  0.588  2.909  4.133 
## 
## Coefficients:
##            Estimate Std. Error t value Pr(>|t|)
## tab16_9A$x  0.96475    0.00513     188  7.2e-16
## 
## Residual standard error: 3.29 on 8 degrees of freedom
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 3.53e+04 on 1 and 8 DF,  p-value: 7.2e-16
S2e=sum(resid(fit16_9)^2)/(9-1)
S2e
## [1] 10.847
sqrt(S2e)
## [1] 3.2934
confint(fit16_9)
##              2.5 %  97.5 %
## tab16_9A$x 0.95291 0.97659
plot(tab16_9A$x,tab16_9A$y, pch=20, xlab="x", ylab="y", col="darkblue")
abline(lm(tab16_9A$y~tab16_9A$x), lwd=2, col="red")

Figura 16.11: Dados e reta ajustada para o Exemplo 16.8.

16.6.2 Modelo Não Lineares

Tabela 16.6

ano=seq(1961,1979,2)
inflacao_Y=c(9,24,72,128,192,277,373,613,1236,2639)

tab16_6<-data.frame(
  ano=ano,
  t=ano-mean(ano),
  inflacao_Y=inflacao_Y,
  log_Y=log(inflacao_Y))
Tabela 16.6: Taxa de Inflação no Brasil de 1961 a 1979.
Ano t Inflação(Y) Y*=log Y
1961 -9 9 2.1972
1963 -7 24 3.1780
1965 -5 72 4.2767
1967 -3 128 4.8520
1969 -1 192 5.2575
1971 1 277 5.6240
1973 3 373 5.9216
1975 5 613 6.4184
1977 7 1236 7.1196
1979 9 2639 7.8782

Para construir o gráfico da Figura 16.12, utilizamos o modelo (16.74) que foi ajustado no Exemplo 16.10 seguinte.

plot(t,inflacao_Y, pch=20, col="darkblue", ylim=c(0,max(inflacao_Y)+5),ylab="Inflação",xlab="Ano(t)")
par(new=T)
plot(t,exp(predict(fit16_10)),pch="+", col="darkred", ylim=c(0,max(inflacao_Y)+5),ylab="",xlab="")

Figura 16.12: Dados de inflação no Brasil (pontos) e modelo exponencial ajustado(+).

Exemplo 16.10

Queremos ajustar o modelo(16.68):

\[ y_i=\alpha e^{\beta x_i}\epsilon_i,i=1,...,n\\ \]

Ajustaremos este modelo a partir da equivalente formulação:

\[ log(y_i)=log(\alpha) +\beta x_i+log(\epsilon_i)\\ y_i^* = \alpha^* +\beta x_i +\epsilon_i^*,i=1,...,n \]

fit16_10<-lm(log_Y~t)
fit16_10
## 
## Call:
## lm(formula = log_Y ~ t)
## 
## Coefficients:
## (Intercept)            t  
##       5.272        0.282

O Modelo na forma da equação (16.71) pode ser obtido a partir da exponenciação do parâmetro do intercepto:

exp(fit16_10$coefficients[1])
## (Intercept) 
##      194.87
plot(tab16_6$t,tab16_6$log_Y, pch=20, xlab="x", ylab="y", col="darkblue")
abline(lm(tab16_6$log_Y~tab16_6$t), lwd=2, col="red")

Figura 16.13: Diagrama de dispersão para o logaritmo da inflação com reta ajustada.

resid_exp=inflacao_Y-exp(predict(fit16_10))
tab16_7<-data.frame(t=tab16_6$t,resid_reta=resid(fit16_10),resid_exp)
print(kable(tab16_7,caption="**Tabela 16.7**: Resíduos para os modelos linear e exponencial.", col.names=c("t", "Resíduos Reta", "Resíduos Exponencial")))
Tabela 16.7: Resíduos para os modelos linear e exponencial.
t Resíduos Reta Resíduos Exponencial
-9 -0.53865 -6.4232
-7 -0.12148 -3.0999
-5 0.41348 24.3833
-3 0.42519 44.3334
-1 0.26700 44.9909
1 0.06987 18.6927
3 -0.19623 -80.8674
5 -0.26310 -184.4828
7 -0.12548 -165.2437
9 0.06938 176.8982 
plot(tab16_7$t,tab16_7$resid_reta, xlab="ano", ylab="Resíduos reta", pch=20)
abline(h=0,lty=2, col="gray")

Figura 16.14: Resíduos da reta ajustada ao logarítmo da inflação versus ano.

plot(tab16_7$t,tab16_7$resid_exp, xlab="ano", ylab="Resíduos exponencial", pch=20)
abline(h=0,lty=2, col="gray")

Figura 16.15: Resíduos do modelo exponencial ajustado aos dados originais versus ano.

par(mfrow=c(1,2))
hist(tab16_7$resid_reta,col="lightblue3",border="white", main="(a)",xlab="",ylab="Resíduos reta")
hist(tab16_7$resid_exp,col="lightblue3",border="white", main="(b)",xlab="",ylab="Resíduos exponencial")

Figura 16.16: Histogramas: (a) resíduos reta ajustada ao log(inflação) e (b) resíduos modelo exponencial.

par(mfrow=c(1,2))
qqnorm(tab16_7$resid_reta,col="darkblue", main="(a)",xlab="Quantis da normal padrão",ylab="Resíduos reta", pch=20,cex=1.5)
qqline(tab16_7$resid_reta, col="darkred")
qqnorm(tab16_7$resid_exp,col="darkblue", main="(b)",xlab="Quantis da normal padrão",ylab="Resíduos exponencial", pch=20,cex=1.5)
qqline(tab16_7$resid_exp, col="darkred")

Figura 16.16: Gráficos q x q dos resíduos: (a) reta e (b) exponencial.

Exemplo 16.11

Os grupos E,C e D são dados por :

# Definin
grupo_E=c(2,1,4,3,5,8,6)
grupo_C=c(7,12,10,11,9,14)
grupo_D=c(16,13,15,18,17,20,19)
print(kable(t(tab15_1[grupo_E,c(1,4,2)]),caption = "**Grupo E**"))
Grupo E
2 1 4 3 5 8 6
tempo_Y 92 96 100 106 98 101 104
acuidade_Z 100 90 90 80 100 90 90
sexo_W M H M H M M H 
print(kable(t(tab15_1[grupo_C,c(1,4,2)]),caption = "**Grupo C**"))
Grupo C
7 12 10 11 9 14
tempo_Y 110 100 106 109 116 105
acuidade_Z 80 80 90 90 70 80
sexo_W H M H H M M 
print(kable(t(tab15_1[grupo_D,c(1,4,2)]),caption = "**Grupo D**"))
Grupo D
16 13 15 18 17 20 19
tempo_Y 108 112 118 112 113 117 127
acuidade_Z 90 90 70 90 90 80 60
sexo_W H M H M M H H

Vamos programar uma função para calcular a regressão resistente:

reg_resist<-function(tab,plot=TRUE){
  # tab é um data.frame contendo as colunas: x e y.
  
  tab<-tab[order(tab$x,tab$y),]
  g_E<-tab[1:ceiling(length(tab$x)/3),]
  g_C<-tab[(ceiling(length(tab$x)/3)+1):(floor(length(tab$x)/3*2)),]
  g_D<-tab[(floor(length(tab$x)/3*2)+1):(length(tab$x)),]
  
  md_y_E<-median(g_E[,2])
  md_x_E<-median(g_E[,1])
  
  md_y_C<-median(g_C[,2])
  md_x_C<-median(g_C[,1])
  
  md_y_D<-median(g_D[,2])
  md_x_D<-median(g_D[,1])
  
  tab_resist=data.frame(x=c(md_x_E,md_x_C,md_x_D),y=c(md_y_E,md_y_C,md_y_D))
  fit<-lm(tab_resist$y~tab_resist$x)
  
  if (plot == TRUE){
    plot(tab$x,tab$y, pch=20, xlab="x", ylab="y", col="darkblue")
    abline(lm(tab$y~tab$x), lwd=2, col="red")
    abline(lm(tab_resist$y~tab_resist$x), lwd=2, col="blue", lty=2)
    legend("topleft", legend =c("Min.Quad","Resistente") ,lty=c(1,2),col=c("red","blue"),lwd=c(2,2))
  }
  return(fit)
}
fit16_11<-reg_resist(data.frame(x=n_idade_X,y=tab15_1$tempo_Y))

Figura 16.19: Reta MQ e Resistente para o Exemplo 16.11.

O objeto fit_resist_16_11 armazena o resultado dos coeficientes ajustados para a reta resistente:

fit16_11
## 
## Call:
## lm(formula = tab_resist$y ~ tab_resist$x)
## 
## Coefficients:
##  (Intercept)  tab_resist$x  
##        87.33          0.65

Exemplo 16.12

No Exemplo 4.13 tínhamos:

Leitura dos dados:

cd_mercado <- read.table("cd-mercado.csv",h=T,skip=4,sep=";",dec=",") # Leitura dos dados
attach(cd_mercado)
fit16_12<-lm(telebras[1:39]~indice[1:39])

plot(indice[1:39],telebras[1:39],ylab="Tel",xlab="IBV",pch=16,col="darkblue")
abline(fit16_12)

Figura 16.20: Gráfico de dispersão das variáveis X e Y para o Exemplo 16.12 e reta ajustada.

Os resultados do Quadro 16.1 podem ser impressos utilizando os comandos summary e anova:

summary(fit16_12)
## 
## Call:
## lm(formula = telebras[1:39] ~ indice[1:39])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.2380 -0.6073  0.0599  0.5410  1.4595 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)   -5.2893     1.0803    -4.9  1.9e-05
## indice[1:39]   0.9168     0.0297    30.9  < 2e-16
## 
## Residual standard error: 0.749 on 37 degrees of freedom
## Multiple R-squared:  0.963,  Adjusted R-squared:  0.962 
## F-statistic:  952 on 1 and 37 DF,  p-value: <2e-16
anova(fit16_12)
## Analysis of Variance Table
## 
## Response: telebras[1:39]
##              Df Sum Sq Mean Sq F value Pr(>F)
## indice[1:39]  1    534     534     952 <2e-16
## Residuals    37     21       1
resid_fit16_12<-resid(fit16_12)
par(mfrow=c(2,2))
qqnorm(resid_fit16_12, pch=20, col="darkblue", xlab="Normal Score(a)", ylab="Residual")
qqline(resid_fit16_12, col="darkred")
plot.ts(resid_fit16_12, pch=20, col="darkblue",xlab="Obs. Number (b)", ylab="Residual")
abline(h=c(-0.7932,0.7932), col="red", lty=2)
hist(resid_fit16_12, col="lightblue3", border = "white", xlab="Residual(c)", )
plot(fitted(fit16_12),resid_fit16_12, pch=20, col="darkblue", xlab="Fit(d)", ylab="Residual")
abline(h=0)

Figura 16.21: Gráficos após o ajuste do modelo: análise de resíduos, Exemplo 16.12

Os gráficos da Figura 16.21 são gráficos da análise de resíduos do Minitab. O R imprime difrerentes gráficos em sua análise de resíduos através do comando plot(), a saber:

par(mfrow=c(2,2))
plot(fit16_12)

Exemplo 16.21:

tab16_8<-data.frame(t=1:15, vt=c(22.2,61.1,13,27.8,22.2,7.4,7.4,7.4,20.4,20.4,20.4,11.1,13,7.4,14.8))
print(kable(tab16_8,caption="**Tabela 16.8**: Velocidade do vento no aeroporto de Philadelphia."))
Tabela 16.8: Velocidade do vento no aeroporto de Philadelphia.
t vt
1 22.2
2 61.1
3 13.0
4 27.8
5 22.2
6 7.4
7 7.4
8 7.4
9 20.4
10 20.4
11 20.4
12 11.1
13 13.0
14 7.4
15 14.8 
attach(tab16_8)
fit_16_13<-lm(vt~t)
fit_16_13
## 
## Call:
## lm(formula = vt ~ t)
## 
## Coefficients:
## (Intercept)            t  
##       30.03        -1.45
fit_resist16_13<-reg_resist(data.frame(x=t,y=vt))

Figura 16.22: Reta de MQ e resistente para os dados de velocidade do vento.

fit_resist16_13
## 
## Call:
## lm(formula = tab_resist$y ~ tab_resist$x)
## 
## Coefficients:
##  (Intercept)  tab_resist$x  
##        21.56         -0.92

Capítulo Anterior | Introdução