Esta página dedica-se a dar suporte ào livro “Modelo de regressão linear. Teoria e aplicações” da minha autoria. O arquivo contendo o conteúdo desde livro assim como os dois primeiros capítulos poder ser acessado aqui.

Utilizo também este ambiente para oferecer suporte na resolução de alguns exemplos; sempre utilizando a linguagem de programação R, version 4.3.1 (2023-06-16) – “Beagle Scouts”, Copyright (C) 2023 The R Foundation for Statistical Computing. Também o leitor pode encontrar aqui informações adicionais referentes aos exemplos no referido livro assim como os arquivos de dados.



Capítulo 1 Introdução

Exemplo 1.2


Utilizaremos comandos básicos no R. Criamos dois vetores veloc e distan com os dados da velocidade até quando é acionado o freio e a distância percorrida até parar, respectivamente.

veloc = c(4,7,17,14,12,11,20,15,17,13,15,19,10,18,22,18,8,4,12,20,23,18,12,16,18, 
          19,24,14,12,9,10,15,24,25,20,19,13,10,7,16,14,20,24,24,17,13,11,13,14,20)*1.60934
distan = c(2,4,50,36,20,28,48,54,40,34,26,68,26,56,66,84,16,10,14,56,54,76,24,32,42, 
           46,93,26,28,10,34,20,70,85,64,36,26,18,22,40,60,52,120,92,32,34,17,46,80,32)*0.3048
par(mar=c(5,5,1,1), pch=19, cex.axis=0.6)
plot(veloc, distan, xlab="Velocidade em quilômetros por hora \n quando o sinal é dado", 
     ylab=" Distância percorrida após o sinal \n antes de parar, em metros")
grid()

Figura 1.2: Gráfico da relação entre a velocidade do automóvel com a distância que leva para parar, como mostrado por observações individuais.

Exemplo 1.3


Para este exemplo utilizamos somente o comando qplot na biblioteca R qplot, por isso a chamada espefifica ggplot2::. Posteriormente utilizamos somente o comando plot_grid, no pacote cowplot para mostrar a matriz de gráficos apresentada na Figura 1.3.

Tamanho = c(24,89,73,32,48,40,69,44,65,93,28,48,97,65,36,44,89,44,65,32)
Vacas = c(18,0,14,6,1,9,6,12,7,2,17,15,7,0,12,16,2,6,12,15)
Renda = c(960,830,1260,610,590,900,820,880,860,760,1020,1080,960,700,800,1130,760,740,980,800)
par(mar=c(1,1,1,1),pch=19,cex.axis=0.6)
graf01 = ggplot2::qplot(Tamanho, Renda, ylab="Renda por hectare", xlab="Tamanho da fazendo em hectare")
graf02 = ggplot2::qplot(Vacas, Renda, ylab="Renda por hectare", xlab="Número de vacas")
cowplot::plot_grid(graf01, graf02, labels = c("A", "B"), ncol = 2, nrow = 1)

Figura 1.3: Gráfico descritivo da possível relação entre o tamanho em hectares da fazenda e o número de vacas com a renda obtida. Este gráfico é chamado de gráfico de dispersão.

Figura 1.5 Densidade normal 


par(mar=c(2,4,1,1), lwd = 2)
x = seq(-5,5,by=0.01)
plot(x,dnorm(x,sd=1), type = "l", axes = FALSE, xlab = "", ylab = "", ylim = c(0,0.85))
lines(x,dnorm(x,sd=0.5), col="blue")
lines(x,dnorm(x,sd=sqrt(2)), col = "red")
box()
axis(2)
text(0,0.45,expression(paste(sigma^2,"=1")))
text(0.85,0.75,expression(paste(sigma^2,"=1/4")), col="blue")
text(2.5,0.15,expression(paste(sigma^2,"=2")), col="red")
axis(1, at=0, labels = expression(mu))

Figura 1.5: Densidade normal para diversos valores de variância e igual média.

Exemplo 1.4.1 Quanto custa ser deputado?


Leitura dos dados e primeiras linhas de informação:

EPOCA = read.csv("http://leg.ufpr.br/~lucambio/Linear/EPOCA.csv", sep=";") 
head(EPOCA)
##   UF                                 Nome Partido N.votos Arrecadado Eleito
## 1 AC ELISABETH APARECIDA GARCIA RODRIGUES    PSDB    1357      440.0      0
## 2 AC               GLADSON DE LIMA CAMELI      PP   32623  2373076.2      1
## 3 AC      SEBASTIAO SIBA MACHADO OLIVEIRA      PT   25158   525793.2      1
## 4 AC             TAUMATURGO LIMA CORDEIRO      PT   17932   505191.2      1
## 5 AC   ANTONIA LUCILEIA CRUZ RAMOS CAMARA     PSC   15849   430291.9      1
## 6 AC            MARIA PERPETUA DE ALMEIDA PC DO B   33235   413048.2      1


Para mostrar algumas estatísticas descritivas utilizamos o pacote vtable. o qual serve ao propósito de gerar documentação automática de variáveis que pode ser facilmente visualizada enquanto continua a trabalhar com dados. A biblioteca vtable contém quatro funções principais e nos contramos na função sumtable ou st.

vtable::st(EPOCA[c(4,5)])
Summary Statistics
Variable N Mean Std. Dev. Min Pctl. 25 Pctl. 75 Max
N.votos 3775 30707 209900 1 732 18412 7565377
Arrecadado 3775 257836 968807 8.2 3327 143870 38036984


Para criar os gráficos na Figura 1.7 utilizamos somente a função xyplot no pacote lattice. A biblioteca ou pacote de funções lattice representa um sistema elegante e poderoso de visualização de dados de alto nível que visa melhorar os gráficos básicos R.

A chamada específica da função xyplot, do pocote lattice, aconte com a sintaxis mostrada a continuação. Caso queira-se utilizar outras funções também dentro da biblioteca de funções lattice, utilizar-se-á library(lattice).

Eleito = ifelse(EPOCA$Eleito == 1, "Sim", "Não")
par(mar = c(4,4,1,1))
lattice::xyplot(EPOCA$N.votos ~ EPOCA$Arrecadado | Eleito, pch =19, type = c("p", "g"), 
                xlab = "Total arrecadado (R$)", ylab = "Número de votos")


Pode-se perceber que estes gráficos não são muito informativos. Existe um forte agrupamento de dados, ou seja, existem muitos deputados com pouca arrecadação de recursos e com poucos votos recebidos, mas também observamos deputados com relativamente altos valores arrecadados e, mesmo assim, poucos votos recebidos. Tudo isso somente observando o gráficos à esquerda, nos deputados não eleitos. No gráfico à direita na figura acima podemos ter afirmações semelhantes.

Faremos a contiuação uma mudaça da escala nos dados, isso não altera a suposta relação entre elas, somente permite melhor visualização a pretendida relação entre as variáveis consideradas.

Gráfico da Figura 1.7 em escala logaritmica:

Eleito = ifelse(EPOCA$Eleito == 1, "Sim", "Não")
par(mar = c(4,4,1,1))
lattice::xyplot(log(EPOCA$N.votos) ~ log(EPOCA$Arrecadado) | Eleito, pch = 19, type = c("p", "g"),
                xlab = "Logaritmo do total arrecadado (R$)", ylab = "Logaritmo do número de votos")


Percebemos que, todos os deputados eleitos, tiveram uma grande arrecadação de dinheiro para a campanha eletoral e também tiveram um grande número de votos. No caso dos deputados não eleitos, a arrecadação de dinheiro para a campanha, as vezes, foi maior do que no caso dos aleitos e inclussiva, o número e votos foi maior.

Observa-se também que, dentre os deputados não aleitos, houveram deputados com alta arrecadação e quase nulo o número de votos recebidos.

Exemplo 1.4.2 Peso dos mexilhões


Leitura dos dados e primeiras linhas de informação:

mexilhoes = read.csv("http://leg.ufpr.br/~lucambio/Linear/mexilhoes.csv", sep=";") 
head(mexilhoes)
##     W  H   L   S  M
## 1 318 68 158 345 47
## 2 312 56 148 290 52
## 3 265 46 124 167 27
## 4 222 38 104  67 13
## 5 274 51 143 238 31
## 6 216 35  99  68 14


Gráficos descritivos contidos na Fogura 1.8.

library(ggplot2)
grafico1 = qplot(W, M, data = mexilhoes, geom = "point", 
                 xlab = "Largura da concha (mm)", ylab = "Massa muscular (g)")
grafico2 = qplot(H, M, data = mexilhoes, geom = "point",
                 xlab = "Altura da concha (mm)", ylab = "Massa muscular (g)")
grafico3 = qplot(L, M, data = mexilhoes, geom = "point",
                 xlab = "Comprimento da concha (mm)", ylab = "Massa muscular (g)")
grafico4 = qplot(S, M, data = mexilhoes, geom = "point",
                 xlab = "Peso da concha (mm)", ylab = "Massa muscular (g)")
cowplot::plot_grid(grafico1, grafico2, grafico3, grafico4, 
                   labels = c("A", "B", "C", "D"), ncol = 2, nrow = 2)


Percebemos em duas situações um comportamento clao, mas nas outras duas não! As variáveis pesdo da conche (S) e altura da concha (H) apresentam comportamento claramente crescente, ou seja, quanto maior elas maior é a massa muscular (M) do mexilhão. Estes são os gráficos B e D, à esquerda na figura acima. Alguma não linearidade é observada no gráfico D, podendo ser do tipo exponencial a relação entre estas duas variáveis.

As outras duas variáveis não mostram-se comportamento direto com a variável resposta. No gráfico em A, claramente quanto maior a largura da concha (W) mas, existem situações nas quais a massa muscular do mexilhão é muito superior às outras observações obtidas com a mesma largura de conha, por exemplo, observemos que para valores de largura de concha entre 220 e 250 temos as seguintes observaçoes:

mexilhoes[mexilhoes$W>220 & mexilhoes$W<250,]
##      W  H   L   S  M
## 4  222 38 104  67 13
## 20 226 38 104  69 13
## 30 230 40 106  68 16
## 35 221 37 104  58 15
## 37 228 46 129 188 33
## 44 241 39 110 104 23
## 54 224 36 107  69 13
## 57 234 37 114  72 17
## 58 221 37 108  74 15
## 61 227 35 118  76 14
## 63 230 47 112 125 18
## 67 246 37 110  90 17
## 70 226 35 111  64 16
## 77 227 38 112  88 15
## 79 242 45 112  61 12


Na linha 37 vemos que a largura da concha (W) é 228 e a massa mucular (M) observada foi de 33, ponto este que destaca-se no centro do gráfico A. Além deste, observamos muitas outras observações nas quais a massa muscular é muito maior do que outras conchas com similares largura de concha.

No gráfico C percebemos existe um comportamento não linear, possivelmente logaritmico, entre comprimento da concha (L) e a massa muscular (M), além de observarmos diversos mexilhões nos quais a massa muscular observada é muito supeior à esperada, fenômeno que acontece para comprimento da concha igual ou supeior a 120 (mm).

Exemplo 1.4.3 Biomassa do capim-marinho


Leitura dos dados e primeiras linhas de informação:

capim.marinho = read.csv("http://leg.ufpr.br/~lucambio/Linear/capim-marinho.csv", sep=";") 
head(capim.marinho)
##   Loc.  Type  BIO SAL   pH       K      Na      Zn
## 1   OI  DVEG  676  33 5.00 1441.67 35185.5 16.4500
## 2   OI  DVEG  516  35 4.75 1299.19 28170.4 13.9900
## 3   OI  DVEG 1052  32 4.20 1154.27 26455.0 15.3276
## 4   OI  DVEG  868  30 4.40 1045.15 25072.9 17.3128
## 5   OI  DVEG 1008  33 5.55  521.62 31664.2 22.3300
## 6   OI SHIRT  436  33 5.05 1273.02 25491.7 12.2778


Gráficos na parte superior da Figura 1.9. Nesta situação utilizamos novamente a função xyplot, nestas situações a relação descritiva entre as variáveis BIO e SAL e BIO e pH são mostradas segundo os diferentes níveis de dois fatores Loc e Type, por isso é que depois do símbolo “|” aparecem multiplicadas estas duas variáveis.

par(mar=c(0,0,0,0))
ngrafico01 = lattice::xyplot(capim.marinho$BIO ~ capim.marinho$SAL 
                             | (capim.marinho$Loc)*(capim.marinho$Type), 
                             pch =19, type = c("p", "g"), xlab = "SAL: percentual de salinidade", 
                             ylab = expression(paste("BIO: biomassa aérea ",gm^{-2})))
ngrafico02 = lattice::xyplot(capim.marinho$BIO ~ capim.marinho$pH 
                             | (capim.marinho$Loc)*(capim.marinho$Type), 
                             pch =19, type = c("p", "g"), xlab = "pH: acidez da água", 
                             ylab = expression(paste("BIO: biomassa aérea ",gm^{-2})))
print(ngrafico01, pos = c(0,0,0.52,1), more = TRUE)
print(ngrafico02, pos = c(0.48,0,1,1))


Algumas relações podem ser observações, por exemplo, quando o tipo de vegetação (Type) é TALL e a localidade (Loc) é Oak Island ou OI percebemos uma relação linear entre o percentual de salinidade (SAL) e a resposta, a biomassa aérea (BIO). Nas outras situações não percebe-se claramente alguma relação linear. Isso não significa que não existam relações lineares entre as variáveis SAL e BIO.

Quando observamos o gráfico à direita, mostrando a relação entre as variáveis pH e BIO, percebemos linearidade forte quando o tipo de vegetação é TALL na localidade Smith Island (SM); quando o tipo de vegetação é DVEG e a localidade OI apresenta-se uma relação não linear quadrática entre pH e BIO.

par(mar=c(0,0,0,0))
ngrafico03 = lattice::xyplot(capim.marinho$BIO ~ capim.marinho$K 
                             | (capim.marinho$Loc)*(capim.marinho$Type), 
                             pch =19, type = c("p", "g"), xlab = "K: potássio em partes por milhão (ppm)", 
                             ylab = expression(paste("BIO: biomassa aérea ",gm^{-2})))
ngrafico04 = lattice::xyplot(capim.marinho$BIO ~ capim.marinho$Na 
                             | (capim.marinho$Loc)*(capim.marinho$Type), 
                             pch =19, type = c("p", "g"), xlab = "Na: sódio em ppm", 
                             ylab = expression(paste("BIO: biomassa aérea ",gm^{-2})))
ngrafico05 = lattice::xyplot(capim.marinho$BIO ~ capim.marinho$Zn 
                             | (capim.marinho$Loc)*(capim.marinho$Type), 
                             pch =19, type = c("p", "g"), xlab = "Zn: zinco em ppm", 
                             ylab = expression(paste("BIO: biomassa aérea ",gm^{-2})))
print(ngrafico03, pos = c(0,0,0.52,1), more = TRUE)
print(ngrafico04, pos = c(0.48,0,1,1))


ùltimo gráfico descritivo.

ngrafico05



Capítulo 2 Estimação e testes de hipóteses

Exemplo 2.4 Custo 


Custo = c(9.73,9.61,8.15,6.98,5.87,4.98,5.09, 4.79,4.02,4.46,3.82)
Producao = c(100,120,140,160,180,200,220,240, 260,280,300)


qplot(Producao, Custo, geom = "point", xlab = "Produção", ylab = "Custo")


n = length(Custo)
y.media = mean(Custo)
x.media = mean(Producao)
sum.xy = sum(Custo*Producao)
sum.x.quadrado = sum(Producao^2)
beta.1 = (sum.xy-n*y.media*x.media)/(sum.x.quadrado-n*x.media^2)
beta.0 = y.media-x.media*beta.1
beta.0
## [1] 12.29091
beta.1
## [1] -0.03077273


ajuste = lm(Custo ~ Producao)
summary(ajuste)
## 
## Call:
## lm(formula = Custo ~ Producao)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1564 -0.4091 -0.1154  0.6386  1.0118 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.290909   0.749146  16.407 5.17e-08 ***
## Producao    -0.030773   0.003571  -8.616 1.22e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7491 on 9 degrees of freedom
## Multiple R-squared:  0.8919, Adjusted R-squared:  0.8799 
## F-statistic: 74.24 on 1 and 9 DF,  p-value: 1.218e-05


Exemplo 2.5 Health


Leitura dos dados e últimas linhas de informação:

Health = read.csv("http://leg.ufpr.br/~lucambio/Linear/Health.csv", sep=";") 
tail(Health)
##         Peso Pulso   Pernas Treino Tempo
## 25  67.58526    57 78.47147     68   374
## 26  73.02837    57 78.47147     65   309
## 27  89.81128    59 99.79031     62   367
## 28 111.13012    70 98.88313     69   469
## 29  63.95652    63 87.54332     60   252
## 30  80.28584    53 83.00740     75   338


Estudo descritvo

par(mar=c(0,0,0,0),pch=19,cex.axis=0.6)
psych::pairs.panels(Health,smooth=F,cex=1.5,ellipses=F,pch=20)


O gráfico acima mostra o estudo descritivo na parte inferior, os valores dos coeficientes de correlação na parte superior e na diagonal mostra os histogramas e densidades estimadas de cada variável no arquivo de dados.

ajuste.Health = lm(Tempo ~ ., data=Health)
summary(ajuste.Health)
## 
## Call:
## lm(formula = Tempo ~ ., data = Health)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -55.223 -18.821  -5.321  18.928  44.487 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.6186    56.1027  -0.064 0.949086    
## Peso          2.7947     0.6324   4.419 0.000168 ***
## Pulso        -0.5252     0.8628  -0.609 0.548194    
## Pernas       -1.1134     0.5422  -2.054 0.050614 .  
## Treino        3.9030     0.7477   5.220 2.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.67 on 25 degrees of freedom
## Multiple R-squared:  0.8531, Adjusted R-squared:  0.8296 
## F-statistic:  36.3 on 4 and 25 DF,  p-value: 4.51e-10


Figura 2.5 


novos.valores = seq(100,300,by=10)
novos = data.frame(Producao = novos.valores)
preditos.lim = predict(ajuste, newdata=novos, interval="confidence", level=0.95)
grafico.Custo1 = qplot(Producao, Custo, xlab="Produção", ylab="Custo unitário",
                       ylim=c(3,12), main="Intervalo de confiança") + 
                       geom_line(aes(y = preditos.lim[,1], x = novos.valores)) +
                       geom_line(aes(y = preditos.lim[,2], x = novos.valores)) +
                       geom_line(aes(y = preditos.lim[,3], x = novos.valores))
preditos.lim1 = predict(ajuste, newdata=novos, interval="prediction", level=0.95)
grafico.Custo2 = qplot(Producao, Custo, xlab="Produção", ylab="Custo unitário",
                       ylim=c(3,12), main="Intervalo de predição") + 
                       geom_line(aes(y = preditos.lim1[,1], x = novos.valores)) +
                       geom_line(aes(y = preditos.lim1[,2], x = novos.valores)) +
                       geom_line(aes(y = preditos.lim1[,3], x = novos.valores))
library(ggpubr)
ggarrange(grafico.Custo1,grafico.Custo2)

Figura 2.5: Gráficos de dispersão, média estimada e intervalos de confiança dos dados apresentados no Exemplo 2.4. A esquerda mostramos o gráfico de dispersão junto a reta estimada e o intervalo de confiança para a média. A direita mostramos o gráfico de dispersão, a média estimada junto ao intervalo de confiança para novas observações.

Figura 2.7 


summary(ajuste, correlation = TRUE)
## 
## Call:
## lm(formula = Custo ~ Producao)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.1564 -0.4091 -0.1154  0.6386  1.0118 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.290909   0.749146  16.407 5.17e-08 ***
## Producao    -0.030773   0.003571  -8.616 1.22e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7491 on 9 degrees of freedom
## Multiple R-squared:  0.8919, Adjusted R-squared:  0.8799 
## F-statistic: 74.24 on 1 and 9 DF,  p-value: 1.218e-05
## 
## Correlation of Coefficients:
##          (Intercept)
## Producao -0.95
library(car)
confidenceEllipse(ajuste, col = "red")

Figura 2.7: Elipse de confiança no modelo de regressão do Exemplo 2.4. Coeficientes da regressão altamente dependentes negativamente.