1 Estatistica Aplicada - AULA 3

1.1 Pacotes

library(descr)
library(kableExtra)
library(car)
library(carData)
library(sciplot)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(GGally)
library(plotly)
library(ExpDes.pt)
library(DT)
library(readxl)
library(xlsx)
library(relaimpo)
library(hrbrthemes)

1.2 Regressao Linear Multipla

Etapas principais para a modelagem de regressao:

  1. Definir a variavel dependente;

  2. Selecionar variaveis independentes;

  3. Coleta dos dados pareados das respectivas variaveis;

  4. Ajustar o modelo obtendo os coeficientes da equacao e valor-p;

  5. Reavaliacao das variaveis que permanecerao no modelo;

  6. Ajustar o modelo novamente com as variaveis definitivas, obtendo os coeficientes da equacao e valor-p;

  7. Calculo dos residuos para avaliar o ajuste do modelo;

  8. Revisao para redefinicao e possivel transformacao de variaveis para melhoria do ajuste do modelo;

  9. Interpretacao do modelo ajustado (equacao) e calculos de previsao;

  10. Utilizacao da equacao para predicao da base que nao apresenta a variavel resposta;

  11. Exportacao dos dados gerados no modelo.

1.3 Exemplo 1:

Uma instituicao de ensino superior deseja avaliar se e possivel identificar fatores associados ao desempenho de seus alunos formandos em uma avaliacao do MEC. As variaveis explicativas testadas foram nota media em disciplinas do curso, idade, sexo e distancia da residencia ao campus (local da prova), para uma amostra aleataria de 127 alunos.

Dados do exemplo1.

bd <- read_excel("C:/LEONARDO/Rpubs/estatistica_aplicada_aula3_leonardo.xlsx")
datatable(
  bd, extensions = 'Buttons', options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
  )
)

1.3.1 a. A variavel resposta Y:

Avaliacao_mec

Como e um exemplo para calculos de previsao, algumas observacoes (alunos) estao sem valores da variavel Avaliacao_mec, portanto, incialmente devemos separar a parte da base de dados com os casos que apresentem essa variavel resposta, formando uma nova base de dados para treino do modelo.

bd_modelo <- bd %>% 
  na.omit()
bd_teste <- setdiff(bd, bd_modelo)

#opcoes para retirar NA
#bd_modelo<- subset.data.frame(bd, Avaliacao_mec>=0)
#bd_modelo<- subset.data.frame(bd, Avaliacao_mec!="NA")
str(bd_modelo)
## tibble [127 x 14] (S3: tbl_df/tbl/data.frame)
##  $ X1_CURSO     : chr [1:127] "Administração" "Administração" "Administração" "Administração" ...
##  $ X2_DISC      : chr [1:127] "Planejamento contábil" "Planejamento contábil" "Planejamento contábil" "Planejamento contábil" ...
##  $ X3_EXIG_MAT  : num [1:127] 2 2 2 2 2 2 2 2 2 2 ...
##  $ X4_MODAL     : chr [1:127] "1" "1" "1" "1" ...
##  $ X5_NOTAS_ANT : num [1:127] 5.8 6.4 9 9.2 7.3 9.3 7.5 7 9.4 6.4 ...
##  $ X6_NUM_EVASAO: num [1:127] 3 2 0 0 0 0 1 2 0 2 ...
##  $ X7_NUM_ATRAS : num [1:127] 0 0 1 0 2 1 0 0 0 0 ...
##  $ X8_SEXO      : num [1:127] 0 0 0 0 1 1 1 0 1 0 ...
##  $ X9_IDADE     : num [1:127] 30 27 23 25 30 28 25 25 23 25 ...
##  $ X10_DIST     : num [1:127] 1 3 2 2 3 1 2 1 1 1 ...
##  $ X11_SEMESTRE : chr [1:127] "1" "1" "1" "1" ...
##  $ X12_TAM_TURMA: num [1:127] 25 30 30 33 30 43 25 35 42 20 ...
##  $ Avaliacao_mec: num [1:127] 62 71 87 93.2 82.4 81 74.6 71.8 68.3 72.2 ...
##  $ Y_SITUACAO   : num [1:127] 0 0 0 0 0 0 0 0 0 0 ...
##  - attr(*, "na.action")= 'omit' Named int [1:243] 13 14 15 16 17 18 19 20 21 22 ...
##   ..- attr(*, "names")= chr [1:243] "13" "14" "15" "16" ...

Visualizacao de um histograma da variavel resposta Avaliacao_mec na base de dados de modelagem.

hist(bd_modelo$Avaliacao_mec, col="gray", las=1, xlab="Avaliacao_mec", ylab="Frequência", 
      main="histograma", cex.lab=1.25, cex.axis=1)
abline(v=mean(bd_modelo$Avaliacao_mec), col="red", lty=3)

1.3.2 b. Selecionar variaveis independentes para o modelo:

1.3.3 c. Coleta dos dados pareados das respectivas variaveis:

X5_NOTAS_ANT; X9_IDADE; X8_SEXO; X10_DIST

par(mfrow=c(2,2))

#X5_NOTAS_ANT
plot(bd_modelo$Avaliacao_mec,bd_modelo$X5_NOTAS_ANT, pch=16, col="black", ylab="NOTAS_ANT", xlab="Avaliacao_mec", las=1)

#X9_IDADE
plot(bd_modelo$Avaliacao_mec,bd_modelo$X9_IDADE, pch=16, col="black", ylab="IDADE", xlab="Avaliacao_mec", las=1)

#X8_SEXO
plot(bd_modelo$Avaliacao_mec,bd_modelo$X8_SEXO, pch=16, col="black", ylab="SEXO", xlab="Avaliacao_mec", las=1)

#X10_DIST
plot(bd_modelo$Avaliacao_mec,bd_modelo$X10_DIST, pch=16, col="black", ylab="DIST", xlab="Avaliacao_mec", las=1)

#X5_NOTAS_ANT
correlacao.X5_NOTAS_ANT <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X5_NOTAS_ANT)
#X9_IDADE
correlacao.X9_IDADE <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X9_IDADE)
#X8_SEXO
correlacao.X8_SEXO <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X8_SEXO)
#X10_DIST
correlacao.X10_DIST <- cor(bd_modelo$Avaliacao_mec,bd_modelo$X10_DIST)

Os valores de correlacao entre a variavel resposta Avaliacao_mec e cada uma das variaveis explicativas sao:

X5_NOTAS_ANT: 0.4062551.

X9_IDADE: 0.5934118.

X8_SEXO: 0.0138624.

X10_DIST: -0.026586.

bd_modelo.ggpairs <- bd_modelo[, c(5, 8, 9, 10, 13)]
ggpairs(bd_modelo.ggpairs)

1.3.4 d. Ajustar o modelo obtendo os coeficientes da equacao e valor-p:

resultado.rlm <- lm(Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO + X10_DIST, 
                    data=bd_modelo)
summary(resultado.rlm)
## 
## Call:
## lm(formula = Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO + 
##     X10_DIST, data = bd_modelo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.3539  -4.4799  -0.8912   3.0237  29.2744 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.47420    6.26408  -0.395    0.694    
## X5_NOTAS_ANT  4.42930    0.58765   7.537 9.36e-12 ***
## X9_IDADE      1.72169    0.16592  10.377  < 2e-16 ***
## X8_SEXO      -2.70870    1.27075  -2.132    0.035 *  
## X10_DIST      0.05424    0.81972   0.066    0.947    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.971 on 122 degrees of freedom
## Multiple R-squared:  0.559,  Adjusted R-squared:  0.5446 
## F-statistic: 38.67 on 4 and 122 DF,  p-value: < 2.2e-16
r2.resultado.rlm <- summary(resultado.rlm)$r.squared
r2.resultado.rlm <- round(r2.resultado.rlm*100, digits = 2)
adj.r2.resultado.rlm <- summary(resultado.rlm)$adj.r.squared

O modelo utilizado com as 4 variaveis explicou (r.squared) 55.9% da variacao dos dados de Avaliacao_mec. O valor de r2 ajustado ficou em 0.5445919.

A seguir vemos a equacao gerada do modelo de regressao linear multiplo com seus coeficientes estao descritos:

equacao.resultado.rlm <- as.formula(
  paste0("Avaliacao_meq ~ ", round(coefficients(resultado.rlm)[1],4), " + ", 
    paste(sprintf("%.2f * %s", 
                  coefficients(resultado.rlm)[-1],  
                  names(coefficients(resultado.rlm)[-1])), 
          collapse=" + ")
  )
)
equacao.resultado.rlm
## Avaliacao_meq ~ -2.4742 + 4.43 * X5_NOTAS_ANT + 1.72 * X9_IDADE + 
##     -2.71 * X8_SEXO + 0.05 * X10_DIST

Os coeficientes de cada uma das variaveis explicativas utilizadas no modelo estao destacadas a seguir.

resultado.rlm$coefficients
##  (Intercept) X5_NOTAS_ANT     X9_IDADE      X8_SEXO     X10_DIST 
##  -2.47419518   4.42930473   1.72169305  -2.70869542   0.05424299
p.value.coef.X5_NOTAS_ANT <- coef(summary(resultado.rlm))["X5_NOTAS_ANT","Pr(>|t|)"]
p.value.coef.X9_IDADE <- coef(summary(resultado.rlm))["X9_IDADE","Pr(>|t|)"]
p.value.coef.X8_SEXO <- coef(summary(resultado.rlm))["X8_SEXO","Pr(>|t|)"]
p.value.coef.X10_DIST <- coef(summary(resultado.rlm))["X10_DIST","Pr(>|t|)"]

O valor-p para o coeficiente de X5_NOTAS_ANT resultou em: 9.364891610^{-12}.

O valor-p para o coeficiente de X9_IDADE resultou em:1.813273610^{-18}.

O valor-p para o coeficiente de X8_SEXO resultou em:0.0350484.

O valor-p para o coeficiente de X10_DIST resultou em:0.9473487.

1.3.5 e. Reavaliacao das variaveis que permanecerao no modelo:

Para avaliacao do modelo utilizamos a verificacao dos valores de p das variaveis explicativas. Observa-se que a variavel X10_DIST apresentou 0.9473487 nao significativo, portanto sera retirada do modelo.

Busca-se tambem fazer uma visualizacao geral do modelo e verificar a importancia de cada uma das variaveis explicativas (X1, X2, X3, … , Xn) sobre a variavel resposta (Y).

avaliacao.resultado.rlm <- calc.relimp(resultado.rlm, type = "lmg")
avaliacao.resultado.rlm
## Response variable: Avaliacao_mec 
## Total response variance: 106.6941 
## Analysis based on 127 observations 
## 
## 4 Regressors: 
## X5_NOTAS_ANT X9_IDADE X8_SEXO X10_DIST 
## Proportion of variance explained by model: 55.9%
## Metrics are not normalized (rela=FALSE). 
## 
## Relative importance metrics: 
## 
##                       lmg
## X5_NOTAS_ANT 0.1833746094
## X9_IDADE     0.3688272024
## X8_SEXO      0.0064981174
## X10_DIST     0.0003493333
## 
## Average coefficients for different model sizes: 
## 
##                      1X        2Xs        3Xs         4Xs
## X5_NOTAS_ANT  3.8883631  4.0268083  4.2065801  4.42930473
## X9_IDADE      1.6255199  1.6486793  1.6806107  1.72169305
## X8_SEXO       0.2859671 -0.6292728 -1.6268172 -2.70869542
## X10_DIST     -0.3619761 -0.2482600 -0.1128168  0.05424299
#dataframe com os resultados da avaliacao.resultado.rlm
impo.variaveis <- data.frame(avaliacao.resultado.rlm$lmg, avaliacao.resultado.rlm$lmg.rank, avaliacao.resultado.rlm$namen[-1])
#rank de $lmg
impo.variaveis <- arrange(impo.variaveis, desc(avaliacao.resultado.rlm$lmg))

impo.variaveis$p.values <- c(round(p.value.coef.X5_NOTAS_ANT, digits = 5), 
                             round(p.value.coef.X9_IDADE, digits = 5), 
                             round(p.value.coef.X8_SEXO, digits = 5), 
                             round(p.value.coef.X10_DIST, digits = 5))

datatable(head(impo.variaveis), colnames = c('lmg', 'rank', 'Variaveis', "valor-p"), 
          caption = 'Tabela 1: Rank de importancia das variedades explicativas no modelo de regressao linear multipla.') %>%
  formatStyle('p.values',
    backgroundColor = styleInterval(0.05, c('green', 'red'))
  )

As variaveis que permanecerao no modelo sao: X5_NOTAS_ANT, X9_IDADE e X8_SEXO.

Retiramos do modelo X10_DIST.

1.3.6 f. Ajustar o modelo novamente com as variaveis definitivas, obtendo os coeficientes da equacao e valor-p:

resultado.rlm.ajust <- lm(Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO, data=bd_modelo)
summary(resultado.rlm.ajust)
## 
## Call:
## lm(formula = Avaliacao_mec ~ X5_NOTAS_ANT + X9_IDADE + X8_SEXO, 
##     data = bd_modelo)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.3283  -4.4596  -0.9163   2.9915  29.2384 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -2.3745     6.0555  -0.392   0.6956    
## X5_NOTAS_ANT   4.4290     0.5852   7.568 7.72e-12 ***
## X9_IDADE       1.7212     0.1651  10.427  < 2e-16 ***
## X8_SEXO       -2.7060     1.2649  -2.139   0.0344 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.942 on 123 degrees of freedom
## Multiple R-squared:  0.559,  Adjusted R-squared:  0.5483 
## F-statistic: 51.98 on 3 and 123 DF,  p-value: < 2.2e-16
r2.resultado.rlm.ajust <- summary(resultado.rlm.ajust)$r.squared
r2.resultado.rlm.ajust <- round(r2.resultado.rlm.ajust*100, digits = 2)
adj.r2.resultado.rlm.ajust <- summary(resultado.rlm.ajust)$adj.r.squared

O modelo utilizado com as 3 variaveis explicou (r.squared) 55.9% da variacao dos dados de Avaliacao_mec. O valor de r2 ajustado ficou em 0.5482782.

A seguir vemos a equacao gerada do modelo de regressao linear multiplo com seus coeficientes estao descritos:

equacao.resultado.rlm.ajust <- as.formula(
  paste0("Avaliacao_meq ~ ", round(coefficients(resultado.rlm.ajust)[1],3), " + ", 
    paste(sprintf("%.2f * %s", 
                  coefficients(resultado.rlm.ajust)[-1],  
                  names(coefficients(resultado.rlm.ajust)[-1])), 
          collapse=" + ")
  )
)
equacao.resultado.rlm.ajust
## Avaliacao_meq ~ -2.374 + 4.43 * X5_NOTAS_ANT + 1.72 * X9_IDADE + 
##     -2.71 * X8_SEXO

Avaliacao da equacao final, com as variaveis explicativas ajustadas.

avaliacao.resultado.rlm.ajust <- calc.relimp(resultado.rlm.ajust, type = "lmg")
avaliacao.resultado.rlm.ajust
## Response variable: Avaliacao_mec 
## Total response variance: 106.6941 
## Analysis based on 127 observations 
## 
## 3 Regressors: 
## X5_NOTAS_ANT X9_IDADE X8_SEXO 
## Proportion of variance explained by model: 55.9%
## Metrics are not normalized (rela=FALSE). 
## 
## Relative importance metrics: 
## 
##                      lmg
## X5_NOTAS_ANT 0.183382480
## X9_IDADE     0.369150704
## X8_SEXO      0.006500253
## 
## Average coefficients for different model sizes: 
## 
##                     1X       2Xs       3Xs
## X5_NOTAS_ANT 3.8883631  4.096111  4.428951
## X9_IDADE     1.6255199  1.660340  1.721192
## X8_SEXO      0.2859671 -1.094494 -2.705966
#dataframe com os resultados da avaliacao.resultado.rlm
impo.variaveis.ajust <- data.frame(avaliacao.resultado.rlm.ajust$lmg, avaliacao.resultado.rlm.ajust$lmg.rank, avaliacao.resultado.rlm.ajust$namen[-1])
#rank de $lmg
impo.variaveis.ajust <- arrange(impo.variaveis.ajust, desc(avaliacao.resultado.rlm.ajust$lmg))


p.value.coef.X5_NOTAS_ANT.ajust <- coef(summary(resultado.rlm.ajust))["X5_NOTAS_ANT","Pr(>|t|)"]

p.value.coef.X9_IDADE.ajust <- coef(summary(resultado.rlm.ajust))["X9_IDADE","Pr(>|t|)"]

p.value.coef.X8_SEXO.ajust <- coef(summary(resultado.rlm.ajust))["X8_SEXO","Pr(>|t|)"]


impo.variaveis.ajust$p.values <- c(round(p.value.coef.X5_NOTAS_ANT.ajust, digits = 5), 
                             round(p.value.coef.X9_IDADE.ajust, digits = 5), 
                             round(p.value.coef.X8_SEXO.ajust, digits = 5))

datatable(head(impo.variaveis.ajust), colnames = c('lmg', 'rank', 'Variaveis', "valor-p"), 
          caption = 'Tabela 2: Rank de importancia das variedades explicativas no modelo de regressao linear multipla ajustado.') %>%
  formatStyle('p.values',
    backgroundColor = styleInterval(0.05, c('green', 'red'))
  )

1.3.7 g. Calculo dos residuos para avaliar o ajuste do modelo:

par(mfrow=c(2,2))

plot(fitted(resultado.rlm.ajust),residuals(resultado.rlm.ajust),xlab="Valores Ajustados",ylab="Residuos", pch=16)
abline(h=0, col="red")

plot(bd_modelo$X5_NOTAS_ANT,residuals(resultado.rlm.ajust),xlab="X5_NOTAS_ANT",ylab="Residuos", pch=16)
abline(h=0, col="red")

plot(bd_modelo$X9_IDADE,residuals(resultado.rlm.ajust),xlab="X9_IDADE",ylab="Residuos", pch=16)
abline(h=0, col="red")

plot(bd_modelo$X8_SEXO,residuals(resultado.rlm.ajust),xlab="X8_SEXO",ylab="Residuos", pch=16)
abline(h=0, col="red")

Avaliacao da suposicao de Normalidade dos erros.

par(mfrow=c(1,2))
qqnorm(resultado.rlm.ajust$residuals, ylab="Residuos",xlab="Quantis teoricos",main="")
qqline(resultado.rlm.ajust$residuals)
hist(resultado.rlm.ajust$residuals, main="histogramas residuos")

shapiro.test(resultado.rlm.ajust$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  resultado.rlm.ajust$residuals
## W = 0.93053, p-value = 6.01e-06

1.3.8 h. Revisao para redefinicao e possivel transformacao de variaveis para melhoria do ajuste do modelo:

1.3.9 i. Interpretacao do modelo ajustado (equacao) e calculos de previsao:

O modelo indicou que 55.9% da variacao da nota dos alunos junto ao MEC foi explicada pela variacao da idade, sexo e media das notas finais de disciplinas conforme a equacao:

equacao.resultado.rlm.ajust
## Avaliacao_meq ~ -2.374 + 4.43 * X5_NOTAS_ANT + 1.72 * X9_IDADE + 
##     -2.71 * X8_SEXO

Intepretacao dos coeficientes B do modelo:

coef.resultado.rlm.ajust.X5_NOTAS_ANT <- resultado.rlm.ajust$coefficients[2]
coef.resultado.rlm.ajust.X9_IDADE <- resultado.rlm.ajust$coefficients[3]
coef.resultado.rlm.ajust.X8_SEXO <- resultado.rlm.ajust$coefficients[4]

X1 (X5_NOTAS_ANT): A cada um ponto a mais na nota media de disciplinas concluidas pelo aluno, observou-se em media, 4.4289509 pontos a mais na avaliacao do MEC.

X2 (X9_IDADE): A cada um ano a mais de idade do aluno, observou-se em media,1.7211919 pontos a mais na avaliacao do MEC.

X3 (X8_SEXO): Em media, as mulheres tem a nota na avaliacao do MEC, -2.7059656 pontos superiores, em relacao aos homens.

1.3.10 j. Utilizacao da equacao para predicao da base que nao apresenta a variavel resposta:

Formacao da base de dados com as variaveis utilizadas na equacao.

bd_new <- subset(bd, select = c('X5_NOTAS_ANT','X9_IDADE','X8_SEXO'))

bd_nova <- bd %>% 
  dplyr::select("X5_NOTAS_ANT" | 'X9_IDADE' | 'X8_SEXO' )

Aplicar o modelo de predicao de valores na base de dados.

predicao <- predict(resultado.rlm.ajust, newdata=bd, se.fit = TRUE, type="response")
bd <- data.frame(bd,predicao)

Avaliacao dos Residuos (Avaliacao_meq - fit) na aplicacao do modelo.

bd$residuos <- bd$Avaliacao_mec - bd$fit

resultado.modelo <-lm(bd$fit~bd$Avaliacao_mec)
summary(resultado.modelo)
## 
## Call:
## lm(formula = bd$fit ~ bd$Avaliacao_mec)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.4559  -3.4715  -0.0137   3.0987  16.1492 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      31.10932    3.16608   9.826   <2e-16 ***
## bd$Avaliacao_mec  0.55903    0.04441  12.588   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.149 on 125 degrees of freedom
##   (243 observations deleted due to missingness)
## Multiple R-squared:  0.559,  Adjusted R-squared:  0.5555 
## F-statistic: 158.5 on 1 and 125 DF,  p-value: < 2.2e-16
plot.resultado.modelo <- ggplot(bd, aes(x=bd$fit, y=bd$Avaliacao_mec)) +
  geom_point() +
  geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
  theme_ipsum()
plot.resultado.modelo + ggtitle("Dispersao de Avaliacao_mec e valores do modelo") +
  xlab("Predicao do modelo") + ylab("Avaliacao_mec")

Erro de estimativa percentual (EEP) de cada observacao da previsao do modelo sobre o valor real da variavel resposta (Avaliacao_mec).

bd$erro_estimativa_relativo <- (bd$fit-bd$Avaliacao_mec)/bd$Avaliacao_mec
bd$erro_estimativa_percentual <- round(((bd$fit-bd$Avaliacao_mec)/bd$Avaliacao_mec)*100, digits = 2)

par(mfrow=c(1,2))
hist(bd$erro_estimativa_percentual , border=F , col=rgb(0.1,0.8,0.3,0.5) , xlab="distribuicao do erro de estimativa percentual" , main="Histograma de EEP")

boxplot(bd$erro_estimativa_percentual , xlab="erro de estimativa percentual" , col=rgb(0.8,0.8,0.3,0.5) , las=2, main="Boxplot de EEP")

Calculo do Erro Relativo Medio das Previsoes (MAPE mean absolute percentage error):

Esta e uma medida de acuracia de predicao do modelo estatistico, quanto menor o valor, melhor este indice.

MAPE <- mean(abs(bd$fit-bd$Avaliacao_mec)/bd$Avaliacao_mec,na.rm = TRUE)*100

O erro relativo medio das previsoes do modelo resultou em 6.967906% (MAPE).

1.3.11 k. Exportacao dos dados gerados no modelo:

bd_final <- write.xlsx(bd, file="estatistica_aplicada_aula3_leonardo_exemplo_finalizado.xls")