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.

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 Poder Calorífico Superior e a Composição Imediata

O Poder calorífico é uma medida da quantidade de energia térmica liberada, na forma de calor, na queima completa de 1kg de combustível. O objetivo principal dessa modelagem é verificar se o Poder Calorífico (PC) tem correlação com a Composição Imediata da biomassa (teor de carbono fixo, voláteis e cinzas). Caso haja associação entre essas variáveis, interpretar a modelagem estatística realizada e a influência de cada variável no valor do PC, 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

tabela4 <- read_excel("Dados.xlsx",
                 sheet = "Analise_imediata")
tabela4 %>%
  kbl(caption = "Tabela dos Dados") %>%
  kable_classic(full_width = F, html_font = "Times", font_size = 12)%>%
  scroll_box(height = "300px")
Tabela dos Dados
CF TV CZ PCS
20.9 78.7 0.4 20.70
10.1 89.9 0.0 19.76
17.2 81.6 1.3 19.84
20.3 79.6 0.0 20.44
17.8 81.3 0.9 18.98
14.2 85.6 0.2 19.91
16.6 80.6 2.9 19.35
19.5 79.4 1.1 19.16
20.7 77.7 1.6 22.22
18.0 76.1 5.9 16.18
27.1 71.0 1.9 20.28
22.5 75.9 1.7 21.10
15.3 84.6 0.1 20.21
9.1 90.9 0.1 19.91
13.3 86.7 0.1 20.14
21.4 78.4 0.2 22.01
24.7 70.6 4.7 18.67
19.6 76.2 4.2 16.55
18.3 78.5 3.2 18.69
20.0 79.6 0.4 19.41
17.8 82.0 0.2 19.72
17.4 82.4 0.2 20.11
16.0 83.2 0.7 19.70
14.8 84.9 0.3 19.84
5.4 93.9 0.8 19.97
19.7 78.4 1.8 19.85
17.6 76.2 6.1 16.65
15.3 84.1 0.7 19.31
20.0 79.4 0.6 22.22
15.6 82.6 1.8 19.88
18.2 80.9 0.9 20.49
17.1 82.0 1.0 18.92
19.3 79.7 1.0 20.55
19.1 76.9 4.1 19.99
19.6 80.1 0.3 19.85
16.7 82.8 1.7 20.20
22.3 74.8 2.9 19.53
19.8 78.0 2.3 19.59
17.6 81.8 0.6 19.14
16.8 82.6 0.7 19.87

2.2 Análise Descritiva

# Poder Calorífico Superior
y <- tabela4$PCS
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.2f $\\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)
PCS<-data.frame(Med_e_sd, Coef_de_Variacao)

# Carbono Fixo
y <- tabela4$CF
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.2f $\\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)
CF<-data.frame(Med_e_sd, Coef_de_Variacao)

# Teor de Voláteis
y <- tabela4$TV
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.2f $\\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)
TV<-data.frame(Med_e_sd, Coef_de_Variacao)

# Teor de Cinzas
y <- tabela4$CZ
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.2f $\\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)
CZ<-data.frame(Med_e_sd, Coef_de_Variacao)

# Junção dos data.frame
dados<-rbind(PCS,CF)
dados<-rbind(dados, TV)
dados<-rbind(dados, CZ)

# Criar uma coluna com o nome das variáveis
dados<-mutate(dados, 
      Variaveis=c("Poder Calorífico Superior (MJ/kg)", 
                  "Carbono Fixo", "Teor de Voláteis", "Teor de Cinzas"))
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(%)
Poder Calorífico Superior (MJ/kg) 19.72 \(\pm\) 1.24 6.30
Carbono Fixo 17.82 \(\pm\) 3.93 22.03
Teor de Voláteis 80.74 \(\pm\) 4.71 5.84
Teor de Cinzas 1.49 \(\pm\) 1.61 107.78

2.2.1 Poder Calorífico Superior - outlier

y <- tabela4$PCS
scale(y) # padronização da variável 
##               [,1]
##  [1,]  0.786308913
##  [2,]  0.030358641
##  [3,]  0.094694835
##  [4,]  0.577216285
##  [5,] -0.596919244
##  [6,]  0.150989004
##  [7,] -0.299364350
##  [8,] -0.452162809
##  [9,]  2.008696587
## [10,] -2.848686012
## [11,]  0.448543898
## [12,]  1.107989880
## [13,]  0.392249729
## [14,]  0.150989004
## [15,]  0.335955560
## [16,]  1.839814080
## [17,] -0.846221993
## [18,] -2.551131117
## [19,] -0.830137945
## [20,] -0.251112205
## [21,] -0.001809455
## [22,]  0.311829487
## [23,] -0.017893504
## [24,]  0.094694835
## [25,]  0.199241149
## [26,]  0.102736859
## [27,] -2.470710876
## [28,] -0.331532446
## [29,]  2.008696587
## [30,]  0.126862931
## [31,]  0.617426406
## [32,] -0.645171389
## [33,]  0.665678551
## [34,]  0.215325197
## [35,]  0.102736859
## [36,]  0.384207705
## [37,] -0.154607915
## [38,] -0.106355770
## [39,] -0.468246857
## [40,]  0.118820907
## attr(,"scaled:center")
## [1] 19.72225
## attr(,"scaled:scale")
## [1] 1.243468
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 
##       5      35

2.2.2 Carbono Fixo - outlier

y <- tabela4$CF
scale(y) # padronização da variável 
##               [,1]
##  [1,]  0.785246830
##  [2,] -1.965982940
##  [3,] -0.157304110
##  [4,]  0.632400732
##  [5,] -0.004458011
##  [6,] -0.921534601
##  [7,] -0.310150208
##  [8,]  0.428605934
##  [9,]  0.734298131
## [10,]  0.046490688
## [11,]  2.364656514
## [12,]  1.192836426
## [13,] -0.641316754
## [14,] -2.220726438
## [15,] -1.150803749
## [16,]  0.912618579
## [17,]  1.753272120
## [18,]  0.454080284
## [19,]  0.122913737
## [20,]  0.555977683
## [21,] -0.004458011
## [22,] -0.106355410
## [23,] -0.462996306
## [24,] -0.768688503
## [25,] -3.163277378
## [26,]  0.479554634
## [27,] -0.055406711
## [28,] -0.641316754
## [29,]  0.555977683
## [30,] -0.564893705
## [31,]  0.097439388
## [32,] -0.182778459
## [33,]  0.377657235
## [34,]  0.326708535
## [35,]  0.454080284
## [36,] -0.284675858
## [37,]  1.141887727
## [38,]  0.505028983
## [39,] -0.055406711
## [40,] -0.259201508
## attr(,"scaled:center")
## [1] 17.8175
## attr(,"scaled:scale")
## [1] 3.925517
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      37

2.2.3 Teor de Voláteis - outlier

y <- tabela4$TV
scale(y) # padronização da variável 
##              [,1]
##  [1,] -0.43272206
##  [2,]  1.94300688
##  [3,]  0.18242204
##  [4,] -0.24181527
##  [5,]  0.11878645
##  [6,]  1.03089666
##  [7,] -0.02969661
##  [8,] -0.28423900
##  [9,] -0.64484071
## [10,] -0.98423056
## [11,] -2.06603570
## [12,] -1.02665429
## [13,]  0.81877801
## [14,]  2.15512553
## [15,]  1.26422718
## [16,] -0.49635765
## [17,] -2.15088316
## [18,] -0.96301869
## [19,] -0.47514579
## [20,] -0.24181527
## [21,]  0.26726951
## [22,]  0.35211697
## [23,]  0.52181189
## [24,]  0.88241360
## [25,]  2.79148150
## [26,] -0.49635765
## [27,] -0.96301869
## [28,]  0.71271868
## [29,] -0.28423900
## [30,]  0.39454070
## [31,]  0.03393898
## [32,]  0.26726951
## [33,] -0.22060340
## [34,] -0.81453564
## [35,] -0.13575594
## [36,]  0.43696443
## [37,] -1.25998481
## [38,] -0.58120511
## [39,]  0.22484577
## [40,]  0.39454070
## attr(,"scaled:center")
## [1] 80.74
## attr(,"scaled:scale")
## [1] 4.714343
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 
##       4      36

2.2.4 Teor de Cinzas - outlier

y <- tabela4$CZ
scale(y) # padronização da variável 
##              [,1]
##  [1,] -0.67875278
##  [2,] -0.92783637
##  [3,] -0.11831470
##  [4,] -0.92783637
##  [5,] -0.36739829
##  [6,] -0.80329458
##  [7,]  0.87801965
##  [8,] -0.24285650
##  [9,]  0.06849799
## [10,]  2.74614657
## [11,]  0.25531068
## [12,]  0.13076888
## [13,] -0.86556547
## [14,] -0.86556547
## [15,] -0.86556547
## [16,] -0.80329458
## [17,]  1.99889580
## [18,]  1.68754132
## [19,]  1.06483234
## [20,] -0.67875278
## [21,] -0.80329458
## [22,] -0.80329458
## [23,] -0.49194009
## [24,] -0.74102368
## [25,] -0.42966919
## [26,]  0.19303978
## [27,]  2.87068837
## [28,] -0.49194009
## [29,] -0.55421099
## [30,]  0.19303978
## [31,] -0.36739829
## [32,] -0.30512740
## [33,] -0.30512740
## [34,]  1.62527042
## [35,] -0.74102368
## [36,]  0.13076888
## [37,]  0.87801965
## [38,]  0.50439427
## [39,] -0.55421099
## [40,] -0.49194009
## attr(,"scaled:center")
## [1] 1.49
## attr(,"scaled:scale")
## [1] 1.605887
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      38

2.3 Gráfico Box Plot

par(mfrow=c(1,4))
boxplot(tabela4$CF, col="blue", main="Box Plot", 
        ylab= "Carbono Fixo (%)", adj=0.5, ylim=c(5,30))

boxplot(tabela4$TV, col="purple",main="Box Plot", 
        ylab= "Teor de Voláteis (%)", adj=0.5, ylim=c(70, 95))

boxplot(tabela4$CZ,col="green2",main="Box Plot", ylab= "Teor de Cinzas (%)", 
        adj=0.5, ylim=c(0,6)) 

boxplot(tabela4$PCS,col="gray",main="Box Plot", ylab= "PCS (MJ/kg)", 
        adj=0.5) 

2.4 Matriz de correlação (valor de r e o pvalor)

chart.Correlation(tabela4, histogram =TRUE)

2.5 Regressão Linear Múltipla

# Ajuste do modelo 
modelo4<-lm(formula =PCS~-1+CF+TV+CZ, data=tabela4); modelo4
## 
## Call:
## lm(formula = PCS ~ -1 + CF + TV + CZ, data = tabela4)
## 
## Coefficients:
##      CF       TV       CZ  
##  0.2907   0.1876  -0.4051
# Estatísticas das estimativas dos paramentros do modelo
summary(modelo4)
## 
## Call:
## lm(formula = PCS ~ -1 + CF + TV + CZ, data = tabela4)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.74134 -0.45320 -0.05022  0.38873  2.27426 
## 
## Coefficients:
##     Estimate Std. Error t value Pr(>|t|)    
## CF  0.290744   0.031081   9.354 2.74e-11 ***
## TV  0.187586   0.006478  28.958  < 2e-16 ***
## CZ -0.405080   0.092066  -4.400 8.84e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8667 on 37 degrees of freedom
## Multiple R-squared:  0.9982, Adjusted R-squared:  0.9981 
## F-statistic:  6919 on 3 and 37 DF,  p-value: < 2.2e-16
# Anova da regressão do modelo
anova(modelo4)
## Analysis of Variance Table
## 
## Response: PCS
##           Df  Sum Sq Mean Sq   F value    Pr(>F)    
## CF         1 14888.0 14888.0 19821.305 < 2.2e-16 ***
## TV         1   688.7   688.7   916.870 < 2.2e-16 ***
## CZ         1    14.5    14.5    19.359 8.843e-05 ***
## Residuals 37    27.8     0.8                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

PCS4_est<- predict(modelo4)
resid<-(tabela4$PCS-PCS4_est)
Res2<-resid^2

SQResiduo<- sum(Res2)

Syx<-sqrt(SQResiduo/((length(tabela4$PCS))-3))

Syx_Porc<-(Syx/(mean(tabela4$PCS)))*100

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

newdata33 = data.frame(Syx,Syx_Porc)

# Criar da Tabela Descritiva
knitr::kable(newdata33, align = 'cc', 
             caption = "Medidas de Precisão do Modelo",
             col.names = c("S$_{yx}$", "S$_{yx}$%")) %>%
  kable_classic(full_width = F, html_font = "Times", font_size = 16)
Medidas de Precisão do Modelo
S\(_{yx}\) S\(_{yx}\)%
0.87 4.39

2.5.2 Predição do PCS aplicando o modelo

newdatas<-data.frame(CF= c(20),TV= c(78), CZ= c(0.4))
PCS_estimado4 = predict(modelo4, newdatas, interval = "confidence");
PCS_estimado4 # Valores do limite de confiança dos valores preditos
##       fit      lwr      upr
## 1 20.2846 19.87058 20.69861

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(modelo4$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo4$residuals
## W = 0.96481, p-value = 0.2436
ks.test(modelo4$residuals, "pnorm")
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  modelo4$residuals
## D = 0.13132, p-value = 0.4568
## 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(modelo4)
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1595594       2.31727   0.366
##  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(modelo4)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo4
## BP = 1.3741, df = 2, p-value = 0.5031
ncvTest(modelo4)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.2500128, Df = 1, p = 0.61707

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

Valores em vermelho estimados; valores em azul observados.

PCS4_est<-predict.lm(modelo4)
newdata4 = data.frame(tabela4, PCS4_est)
Res<-(((newdata4$PCS-newdata4$PCS4_est)/(newdata4$PCS))*100)
newdata4 = data.frame(tabela4, PCS4_est,Res)

G1<-ggplot(newdata4) + 
  geom_point(size=2, aes(x=PCS, y=PCS4_est), color="red2")+
  geom_abline(intercept = 0, slope = 1, color = "black")+
  geom_point(size=2, aes(x=PCS4_est, y=PCS), color="blue2")+
  theme_classic()+ # aparencia clean do gráfico
  labs(x = "PCS Observada (MJ/kg)",
       y = "PCS Estimada (Mj/kg)",
       title = "Gráfico 1:1")+ # legendas
  scale_x_continuous(breaks = seq(15, 25, 5), limits= c(15, 25))+
  scale_y_continuous(breaks = seq(15, 25, 5), limits= c(15, 25))+
  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(newdata4) + 
  geom_point(size=2, aes(x=PCS, y=Res)
             ,color="black")+ # gráfico de pontos
  theme_classic()+ # aparencia clean do gráfico
  labs(x = "PCS Observada (g/cm³)",
       y = "Resíduo (%)",
       title = "Gráfico de Resíduos")+ # legendas
  scale_x_continuous(breaks = seq(15, 25, 5), limits= c(15, 25)) + 
  scale_y_continuous(breaks = seq(-30, 30, 10), limits= c(-30,30)) +
  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) 

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.