1 Instalação dos pacotes

Instalação do pacote e verificação se todos estão instalados. O Comando abaixo ajuda a verificar se os pacotes estão instalados, ou se necessitam ser instalados e habilitados.

pacotes <- c("tidyverse", "readxl", "ggpubr", "kableExtra", 
             "PerformanceAnalytics", "scatterplot3d", "MASS", 
             "gridExtra", "car", "lmtest", "plyr", "knitr", "dplyr")

if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
    break()}
  sapply(pacotes, require, character = T) 
} else {
  sapply(pacotes, require, character = T) 
}
##            tidyverse               readxl               ggpubr 
##                 TRUE                 TRUE                 TRUE 
##           kableExtra PerformanceAnalytics        scatterplot3d 
##                 TRUE                 TRUE                 TRUE 
##                 MASS            gridExtra                  car 
##                 TRUE                 TRUE                 TRUE 
##               lmtest                 plyr                knitr 
##                 TRUE                 TRUE                 TRUE 
##                dplyr 
##                 TRUE

2 Densidade e Propriedades Mecânicas do Combustível

A densidade se refere a quantidade de combustível por unidade de volume. O objetivo principal dessa modelagem é verificar se a Densidade tem correlação com as propriedades mecânicas do combustível (resistência e rigidez). Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência de cada variável no valor da densidade, por meio da comparação da magnitude dos betas da equação de regressão ajustada. Adicionalmente, avaliar a qualidade dos dados experimentais por meio das seguintes estatísticas descritivas: média, desvio padrão, coeficiente de variação e score z. Como referência para a Área de Tecnologia da Madeira, nas disciplinas ministradas pela professora Gilmara, deve-se considerar dados de excelente qualidade se o Coeficiente de variação tiver valor de até 10%, de 10 a 20% dados bons, de 20 a 30 % dados ruins e acima de 30% dados péssimos. Também, no que se refere ao score Z, valores em módulo acima de 2 serão considerados como outliers (pontos atípicos). Os outliers também serão identificados por meio de gráficos box plot das variáveis do modelo.

2.1 Banco de Dados

tabela1 <- read_excel("Dados.xlsx",
                 sheet = "Densidade")
tabela1 %>%
  kbl(caption = "Tabela dos Dados", align = 'cc') %>%
  kable_classic(full_width = F, html_font = "Times", font_size = 12)%>%
  scroll_box(width = "400px", height = "300px")
Tabela dos Dados
Densidade Resistencia Rigidez
688 50.5 12876
1170 79.5 20827
694 59.8 12912
1170 76.7 16694
803 48.1 13481
677 59.1 14098
871 52.0 14613
801 56.0 16224
759 54.8 11105
504 39.0 9839
500 31.5 8058
1090 93.2 23002
838 54.4 13627
1221 83.8 19426
705 47.3 13409
899 48.0 13286
999 62.0 18421
822 51.8 13963
690 48.9 18029
640 40.3 12813
931 63.5 18099
924 48.3 14431
929 54.9 16782
1087 72.7 19881
952 51.6 15561
948 78.5 19360
731 46.8 14933
899 57.7 17198
755 53.9 14617
889 42.7 14577
739 46.0 13166
892 78.4 18359
825 71.4 14624
919 62.4 17212
1068 76.0 18011
1074 93.3 23607
684 56.5 14185
1143 82.9 22733
856 71.4 18971
756 69.9 14719
544 37.8 9067
1106 95.2 21724
940 79.5 19583
580 40.9 15225
579 35.4 8431
537 32.6 7110
535 42.3 9868
560 40.4 11889
538 43.6 10904
645 44.4 13304

2.2 Análise Descritiva

# Densidade
y <- tabela1$Densidade
mean_and_sd <- function(y) {
    m <- mean(y) # média
    s <- sd(y) # desvio padrão
    # sprintf("%0.3f U00B1 %0.3f", m, s)
    sprintf("%0.0f $\\pm$ %0.2f", m, s)}
Med_e_sd<- mean_and_sd(y)
cv <- function(y) {100 * sd(y)/mean(y)}
Coef_de_Variacao<-cv(y) # coeficiente de variação
Coef_de_Variacao<-round(Coef_de_Variacao, 2)
Densidade<-data.frame(Med_e_sd, Coef_de_Variacao)

# Resistencia
y <- tabela1$Resistencia
mean_and_sd <- function(y) {
    m <- mean(y) # média
    s <- sd(y) # desvio padrão
    # sprintf("%0.3f U00B1 %0.3f", m, s)
    sprintf("%0.0f $\\pm$ %0.2f", m, s)}
Med_e_sd<- mean_and_sd(y)
cv <- function(y) {100 * sd(y)/mean(y)}
Coef_de_Variacao<-cv(y) # coeficiente de variação
Coef_de_Variacao<-round(Coef_de_Variacao, 2)
Resistencia<-data.frame(Med_e_sd, Coef_de_Variacao)

#Rigidez
y <- tabela1$Rigidez
mean_and_sd <- function(y) {
    m <- mean(y) # média
    s <- sd(y) # desvio padrão
    # sprintf("%0.3f U00B1 %0.3f", m, s)
    sprintf("%0.0f $\\pm$ %0.2f", m, s)}
Med_e_sd<- mean_and_sd(y)
cv <- function(y) {100 * sd(y)/mean(y)}
Coef_de_Variacao<-cv(y) # coeficiente de variação
Coef_de_Variacao<-round(Coef_de_Variacao, 2)
Rigidez<-data.frame(Med_e_sd, Coef_de_Variacao)

# Junção dos data.frame
dados<-rbind(Densidade,Resistencia)
dados<-rbind(dados, Rigidez)

# Criar uma coluna com o nome das variáveis
dados<-mutate(dados, 
      Variaveis=c("Densidade (g/cm$^{3}$)", "Resistência (MPa)", "Rigidez (MPa)"))
dados<-data.frame(dados)

#Organizar as colunas no data.frame
date<-dplyr::select(dados, Variaveis, Med_e_sd, Coef_de_Variacao) 

# Criar da Tabela Descritiva
knitr::kable(date, align = 'cc', caption = "Tabela Descritiva",
             col.names = c("Variáveis", "Média $\\pm$ sd", "  CV(%)")) %>%
  kable_classic(full_width = F, html_font = "Times", font_size = 16)
Tabela Descritiva
Variáveis Média \(\pm\) sd CV(%)
Densidade (g/cm\(^{3}\)) 822 \(\pm\) 197.37 24.01
Resistência (MPa) 58 \(\pm\) 16.77 28.83
Rigidez (MPa) 15297 \(\pm\) 3930.17 25.69

2.2.1 Densidade - outlier

y <- tabela1$Densidade

scale(y) # padronização da variável 
##                [,1]
##  [1,] -0.6795386825
##  [2,]  1.7625851242
##  [3,] -0.6491388010
##  [4,]  1.7625851242
##  [5,] -0.0968742888
##  [6,] -0.7352717984
##  [7,]  0.2476577006
##  [8,] -0.1070075826
##  [9,] -0.3198067524
## [10,] -1.6118017124
## [11,] -1.6320683000
## [12,]  1.3572533720
## [13,]  0.0804583528
## [14,]  2.0209841161
## [15,] -0.5934056851
## [16,]  0.3895238138
## [17,]  0.8961885040
## [18,] -0.0006079976
## [19,] -0.6694053886
## [20,] -0.9227377337
## [21,]  0.5516565147
## [22,]  0.5161899863
## [23,]  0.5415232208
## [24,]  1.3420534313
## [25,]  0.6580560996
## [26,]  0.6377895120
## [27,] -0.4616728657
## [28,]  0.3895238138
## [29,] -0.3400733400
## [30,]  0.3388573448
## [31,] -0.4211396905
## [32,]  0.3540572855
## [33,]  0.0145919431
## [34,]  0.4908567518
## [35,]  1.2457871402
## [36,]  1.2761870216
## [37,] -0.6998052701
## [38,]  1.6257856578
## [39,]  0.1716579970
## [40,] -0.3350066931
## [41,] -1.4091358363
## [42,]  1.4383197224
## [43,]  0.5972563368
## [44,] -1.2267365478
## [45,] -1.2318031947
## [46,] -1.4446023646
## [47,] -1.4547356584
## [48,] -1.3280694859
## [49,] -1.4395357177
## [50,] -0.8974044992
## attr(,"scaled:center")
## [1] 822.12
## attr(,"scaled:scale")
## [1] 197.3692
score_z <- function(y, media = NULL, despad = NULL) {if (is.null(media)) {
    media <- mean(y)}
  if (is.null(despad)) {despad <- sd(y)}
  (y - media)/despad}
status <- ifelse(abs(score_z(y)) > 2, "outlier", "tipical") # ocorrência de outlier
table(status)
## status
## outlier tipical 
##       1      49

2.2.2 Resistência - outlier

y <- tabela1$Resistencia
scale(y) # padronização da variável 
##              [,1]
##  [1,] -0.45634808
##  [2,]  1.27314673
##  [3,]  0.09828302
##  [4,]  1.10616103
##  [5,] -0.59947868
##  [6,]  0.05653659
##  [7,] -0.36689145
##  [8,] -0.12834044
##  [9,] -0.19990575
## [10,] -1.14218223
## [11,] -1.58946537
## [12,]  2.09018394
## [13,] -0.22376085
## [14,]  1.52958907
## [15,] -0.64718889
## [16,] -0.60544246
## [17,]  0.22948607
## [18,] -0.37881900
## [19,] -0.55176848
## [20,] -1.06465315
## [21,]  0.31894270
## [22,] -0.58755113
## [23,] -0.19394197
## [24,]  0.86761002
## [25,] -0.39074655
## [26,]  1.21350898
## [27,] -0.67700776
## [28,] -0.02695626
## [29,] -0.25357972
## [30,] -0.92152255
## [31,] -0.72471796
## [32,]  1.20754521
## [33,]  0.79008094
## [34,]  0.25334117
## [35,]  1.06441460
## [36,]  2.09614771
## [37,] -0.09852157
## [38,]  1.47591509
## [39,]  0.79008094
## [40,]  0.70062431
## [41,] -1.21374753
## [42,]  2.20945944
## [43,]  1.27314673
## [44,] -1.02887050
## [45,] -1.35687814
## [46,] -1.52386384
## [47,] -0.94537765
## [48,] -1.05868938
## [49,] -0.86784857
## [50,] -0.82013837
## attr(,"scaled:center")
## [1] 58.152
## attr(,"scaled:scale")
## [1] 16.7679
score_z <- function(y, media = NULL, despad = NULL) {if (is.null(media)) {
    media <- mean(y)}
  if (is.null(despad)) {despad <- sd(y)}
  (y - media)/despad}
status <- ifelse(abs(score_z(y)) > 2, "outlier", "tipical") # ocorrência de outlier
table(status)
## status
## outlier tipical 
##       3      47

2.2.3 Rigidez - outlier

y <- tabela1$Rigidez
scale(y) # padronização da variável 
##              [,1]
##  [1,] -0.61592270
##  [2,]  1.40714577
##  [3,] -0.60676279
##  [4,]  0.35553692
##  [5,] -0.46198528
##  [6,] -0.30499456
##  [7,] -0.17395692
##  [8,]  0.23594917
##  [9,] -1.06653951
## [10,] -1.38866311
## [11,] -1.84182434
## [12,]  1.96055716
## [13,] -0.42483675
## [14,]  1.05067251
## [15,] -0.48030511
## [16,] -0.51160147
## [17,]  0.79495828
## [18,] -0.33934423
## [19,]  0.69521701
## [20,] -0.63195255
## [21,]  0.71302795
## [22,] -0.22026537
## [23,]  0.37792782
## [24,]  1.16644363
## [25,]  0.06725411
## [26,]  1.03387934
## [27,] -0.09253547
## [28,]  0.48377570
## [29,] -0.17293915
## [30,] -0.18311683
## [31,] -0.54213452
## [32,]  0.77918288
## [33,] -0.17115806
## [34,]  0.48733788
## [35,]  0.69063705
## [36,]  2.11449458
## [37,] -0.28285810
## [38,]  1.89211226
## [39,]  0.93490139
## [40,] -0.14698606
## [41,] -1.58509234
## [42,]  1.63538026
## [43,]  1.09061990
## [44,] -0.01823840
## [45,] -1.74691747
## [46,] -2.08303537
## [47,] -1.38128429
## [48,] -0.86705697
## [49,] -1.11768236
## [50,] -0.50702152
## attr(,"scaled:center")
## [1] 15296.68
## attr(,"scaled:scale")
## [1] 3930.169
score_z <- function(y, media = NULL, despad = NULL) {if (is.null(media)) {
    media <- mean(y)}
  if (is.null(despad)) {despad <- sd(y)}
  (y - media)/despad}
status <- ifelse(abs(score_z(y)) > 2, "outlier", "tipical") # ocorrência de outlier
table(status)
## status
## outlier tipical 
##       2      48

2.3 Gráfico Box Plot

par(mfrow=c(1,3))
boxplot(tabela1$Densidade, col="blue", main="Box Plot", 
        ylab= "Densidade (g/cm³)", adj=0.5, ylim= c(400,1200))

boxplot(tabela1$Resistencia, col="purple",main="Box Plot", 
        ylab= "Resistência (MPa)", adj=0.5, ylim= c(30,100))

boxplot(tabela1$Rigidez,col="green2",main="Box Plot", ylab= "Rigidez (MPa)", 
        adj=0.5, ylim= c(1000,30000))

2.4 Dispersão, Histograma e Correlação

chart.Correlation(tabela1, histogram =TRUE)

2.5 Regressão Linear Múltipla

# Ajuste do modelo 
modelo1<-lm(formula =Densidade~Resistencia+Rigidez, 
            data=tabela1); modelo1
## 
## Call:
## lm(formula = Densidade ~ Resistencia + Rigidez, data = tabela1)
## 
## Coefficients:
## (Intercept)  Resistencia      Rigidez  
##   164.91421      4.26645      0.02674
# Estatísticas das estimativas dos paramentros do modelo
summary(modelo1)
## 
## Call:
## lm(formula = Densidade ~ Resistencia + Rigidez, data = tabela1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -166.60  -84.37  -22.51   75.49  231.38 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 1.649e+02  5.712e+01   2.887  0.00586 **
## Resistencia 4.266e+00  1.809e+00   2.359  0.02253 * 
## Rigidez     2.675e-02  7.716e-03   3.466  0.00114 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 99.55 on 47 degrees of freedom
## Multiple R-squared:  0.756,  Adjusted R-squared:  0.7456 
## F-statistic:  72.8 on 2 and 47 DF,  p-value: 4.022e-15
# Anova da regressão do modelo
anova(modelo1)
## Analysis of Variance Table
## 
## Response: Densidade
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Resistencia  1 1323934 1323934 133.593 2.442e-15 ***
## Rigidez      1  119063  119063  12.014  0.001139 ** 
## Residuals   47  465778    9910                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.5.1 Medidas de Precisão para Análise do Modelo

Densidade_est<-predict.lm(modelo1)
newdata1 = data.frame(tabela1, Densidade_est)
Res<-(((newdata1$Densidade-newdata1$Densidade_est)/(newdata1$Densidade))*100)
newdata1 = data.frame(tabela1, Densidade_est,Res)

resid<-(newdata1$Densidade-newdata1$Densidade_est)
Res2<-resid^2

SQResiduo<- sum(Res2)

SQTotal<-var(newdata1$Densidade)*((length(newdata1$Densidade))-1)

R2<-(1-(SQResiduo/SQTotal))
a<-(length(newdata1$Densidade)-1)/(length(newdata1$Densidade)-2)

R2adj<-(1-(a*(1-R2)))

Syx<-sqrt(SQResiduo/((length(newdata1$Densidade))-2))

Syx_Porc<-(Syx/(mean(newdata1$Densidade)))*100

R2<-round(R2, 2)
R2adj<-round(R2adj,2)
Syx<-round(Syx,2)
Syx_Porc<-round(Syx_Porc, 2)

newdata11 = data.frame(R2,R2adj,Syx,Syx_Porc)

# Criar da Tabela Descritiva
knitr::kable(newdata11, align = 'cc', 
             caption = "Medidas de Precisão do Modelo",
             col.names = c("R$^{2}$", "R$_{adj}^{2}$", "S$_{yx}$", "S$_{yx}$%")) %>%
  kable_classic(full_width = F, html_font = "Times", font_size = 16)
Medidas de Precisão do Modelo
R\(^{2}\) R\(_{adj}^{2}\) S\(_{yx}\) S\(_{yx}\)%
0.76 0.75 98.51 11.98

2.5.2 Predição da Densidade aplicando o modelo

newdatas<-data.frame(Resistencia= c(50),Rigidez= c(12876))
Densidade_estimado1 = predict(modelo1, newdatas, interval = "confidence");
Densidade_estimado1 # Valores do limite de confiança dos valores preditos
##        fit      lwr      upr
## 1 722.5998 689.0572 756.1425

2.5.3 Testes de hipoteses

NORMALIDADE:
Ho: resíduos seguem distribuição normal;
Hi: resíduos não seguem distribuição normal
Resultado:
Valor de p > 0.05, não rejeita Ho para alfa=5%.

## Normalidade dos Resíduos
shapiro.test(modelo1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo1$residuals
## W = 0.96394, p-value = 0.13
ks.test(modelo1$residuals, "pnorm")
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  modelo1$residuals
## D = 0.54, p-value = 3.642e-14
## alternative hypothesis: two-sided

INDEPENDÊNCIA DOS RESÍDUOS:
Ho: resíduos são independentes;
Hi: resíduos não são independentes
Resultado:
Valor de p > 0.05, não rejeita Ho para alfa=5%.

## Independencia dos Resíduos (Durbin-Watson)
durbinWatsonTest(modelo1)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1013323      1.785325   0.388
##  Alternative hypothesis: rho != 0

HOMOSCEDASTICIDADE DOS RESÍDUOS:
Ho: resíduos são homogênios;
Hi: resíduos não são homogênios;
Resultado:
Valor de p > 0.05, não rejeita Ho para alfa=5%.

## Homocedasticidade (Breusch-Pagan)
bptest(modelo1)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo1
## BP = 1.6976, df = 2, p-value = 0.4279
ncvTest(modelo1)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.9172137, Df = 1, p = 0.33821

2.5.4 Gráfico 1:1 e de Resíduos

Valores em vermelho estimados; valores em azul observados.

G1<-ggplot(newdata1) + 
  geom_point(size=2, aes(x=Densidade, y=Densidade_est), color="red2")+
  geom_abline(intercept = 0, slope = 1, color = "black")+
  geom_point(size=2, aes(x=Densidade_est, y=Densidade), color="blue2")+
  theme_classic()+ # aparencia clean do gráfico
  labs(x = "Densidade Observada (g/cm³)",
       y = "Densidade Estimada (g/cm³)",
       title = "Gráfico 1:1")+ # legendas
  theme( plot.title = element_text(hjust = 0.5),
         legend.title = element_text(hjust = 0.5))+ # centralizar títulos
  theme(plot.title = element_text(size = 13, face = "bold"),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10))  # negrito nos títulos

G2<-ggplot(newdata1) + 
  geom_point(size=2, aes(x=Densidade, y=Res)
             ,color="black")+ # gráfico de pontos
  theme_classic()+ # aparencia clean do gráfico
  labs(x = "Densidade Observada (g/cm³)",
       y = "Resíduo (%)",
       title = "Gráfico de Resíduos")+ # legendas
  scale_x_continuous(breaks = seq(450, 1250, 250), limits= c(450,1250)) + # eixo x
  scale_y_continuous(breaks = seq(-30, 30, 10), limits= c(-30,30)) + # eixo y
  theme( plot.title = element_text(hjust = 0.5),
         legend.title = element_text(hjust = 0.5))+ # centralizar títulos
  theme(plot.title = element_text(size = 13, face = "bold"),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10))+
  # negrito nos títulos
  geom_hline(yintercept = 0, color="black")+ # linha no 0 do eixo y
  theme(axis.title.y = element_text(angle = 90))

grid.arrange(G1,G2, ncol=1)

2.5.5 Gráfico 3D de Dispersão em função do modelo

graph<-scatterplot3d(tabela1$Densidade~tabela1$Resistencia+tabela1$Rigidez,
                     pch=16, angle=45,color="blue", box=FALSE,
                     xlab="Resistência (MPa)", ylab="Rigidez (Mpa)", zlab = "Densidade (g/cm³)")
# previsão do modelo: Quando o modelo é apropriado todos os pontos ficam dentro
#do plano.
graph$plane3d(modelo1, col="black", draw_polygon=TRUE)

3 Citações

De acordo com GROLEMUND (2018), a ciência de dados é uma área em crescimento.

A regressão não linear geralmente possui parâmetros interpretáveis (ROSSE; VENCOVSKY, 2000). Para CHEN et al. (2021), a análise multivariada é a melhor opção.

Referências

CHEN, Xue; LI, Xiaohui; XIE, Jinmei; YANG, Hao; LIU, Aichun. Non-invasive discrimination of multiple myeloma using label-free serum surface-enhanced raman scattering spectroscopy in combination with multivariate analysis. Analytica Chimica Acta, p. 339296, Nov. 2021. DOI 10.1016/j.aca.2021.339296. Available at: https://doi.org/10.1016/j.aca.2021.339296.
GROLEMUND, Garrett. R para data science. S.l: Alta Books, 2018.
ROSSE, Leornado Novaes; VENCOVSKY, Roland. Modelo de regressão não-linear aplicado ao estudo da estabilidade fenotı́pica de genótipos de feijão no estado do paraná. Bragantia, vol. 59, no. 1, p. 99–107, 2000. DOI 10.1590/s0006-87052000000100016. Available at: https://doi.org/10.1590/s0006-87052000000100016.