Poder Calorífico Superior em função da Composição Química da biomassa

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 Química

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 Química da biomassa (teor de lignina, polissacarídeos e extrativos). 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

tabela3 <- read_excel("Dados.xlsx",
                 sheet = "Composicao_quimica")

tabela3 %>%
  kbl(caption = "Tabela dos Dados", align = 'cc') %>%
  kable_classic(full_width = F, html_font = "Times", font_size = 12)

Tabela dos Dados

extrativos	lignina	holocelulose	PCS
2.9434407	29.86403	67.19253	17.52720
2.7417815	29.96502	67.29320	17.64854
2.6133189	30.74463	66.64205	17.67364
1.0326479	30.52474	68.44261	17.48536
1.2695355	29.73810	68.99236	17.49372
0.8540925	29.79510	69.35080	17.58996
1.2518575	29.12709	69.62106	17.62762
1.3169764	29.30707	69.37595	17.78661
1.0532025	28.41579	70.53101	17.59414
1.6052158	29.68110	68.71369	17.07531
1.7641220	30.03399	68.20188	17.32218
1.8730302	29.15542	68.97155	17.25105
2.5166644	27.63171	69.85163	17.34310
1.9856993	27.77445	70.23986	17.29707
1.9593548	28.24153	69.79912	17.36820
1.8144348	29.25707	68.92849	17.43933
1.9256605	29.18541	68.88893	17.38494
1.2597507	30.25092	68.48932	17.42678
1.8393206	28.46292	69.69776	17.31381
1.9124465	28.74425	69.34330	17.46025
1.9287968	28.08158	69.98963	17.54812
2.7028249	29.19124	68.10593	17.80335
3.7316788	31.36373	64.90459	17.76569
3.3813145	30.27092	66.34777	17.74059
1.9874762	30.61163	67.40089	17.48954
2.1662125	31.01380	66.81999	17.41423
2.3980380	30.81384	66.78812	17.58577

2.2 Análise Descritiva

# Poder Calorífico Superior
y <- tabela3$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)

# Extrativos
y <- tabela3$extrativos
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)
extrativos<-data.frame(Med_e_sd, Coef_de_Variacao)

# Lignina
y <- tabela3$lignina
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)
lignina<-data.frame(Med_e_sd, Coef_de_Variacao)

# Holocelulose
y <- tabela3$holocelulose
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)
holocelulose<-data.frame(Med_e_sd, Coef_de_Variacao)

# Junção dos data.frame
dados<-rbind(PCS,extrativos)
dados<-rbind(dados, lignina)
dados<-rbind(dados, holocelulose)

# Criar uma coluna com o nome das variáveis
dados<-mutate(dados, 
      Variaveis=c("Poder Calorífico Superior (MJ/kg)", 
                  "Extrativos", "Lignina", "Holocelulose"))
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)	17.50 \(\pm\) 0.18	1.02
Extrativos	1.99 \(\pm\) 0.71	35.78
Lignina	29.53 \(\pm\) 1.01	3.43
Holocelulose	68.48 \(\pm\) 1.39	2.03

2.2.1 Poder Calorífico Superior - outlier

y <- tabela3$PCS
scale(y) # padronização da variável 

##              [,1]
##  [1,]  0.16212574
##  [2,]  0.84462282
##  [3,]  0.98582911
##  [4,] -0.07321808
##  [5,] -0.02614931
##  [6,]  0.51514147
##  [7,]  0.72695091
##  [8,]  1.62125742
##  [9,]  0.53867585
## [10,] -2.37958750
## [11,] -0.99105897
## [12,] -1.39114346
## [13,] -0.87338706
## [14,] -1.13226526
## [15,] -0.73218077
## [16,] -0.33209628
## [17,] -0.63804324
## [18,] -0.40269942
## [19,] -1.03812773
## [20,] -0.21442437
## [21,]  0.27979765
## [22,]  1.71539494
## [23,]  1.50358551
## [24,]  1.36237922
## [25,] -0.04968370
## [26,] -0.47330257
## [27,]  0.49160709
## attr(,"scaled:center")
## [1] 17.49837
## attr(,"scaled:scale")
## [1] 0.1777867

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      26

2.2.2 Extrativos - outlier

y <- tabela3$extrativos
scale(y) # padronização da variável 

##               [,1]
##  [1,]  1.331445495
##  [2,]  1.048749709
##  [3,]  0.868664554
##  [4,] -1.347198040
##  [5,] -1.015117292
##  [6,] -1.597505782
##  [7,] -1.039899188
##  [8,] -0.948612380
##  [9,] -1.318383529
## [10,] -0.544544085
## [11,] -0.321781574
## [12,] -0.169108733
## [13,]  0.733169486
## [14,] -0.011163553
## [15,] -0.048094631
## [16,] -0.251250556
## [17,] -0.095328870
## [18,] -1.028834087
## [19,] -0.216364485
## [20,] -0.113852965
## [21,] -0.090932378
## [22,]  0.994138401
## [23,]  2.436436619
## [24,]  1.945278655
## [25,] -0.008672669
## [26,]  0.241888773
## [27,]  0.566873104
## attr(,"scaled:center")
## [1] 1.993663
## attr(,"scaled:scale")
## [1] 0.7133434

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      26

2.2.3 Lignina - outlier

y <- tabela3$lignina
scale(y) # padronização da variável 

##             [,1]
##  [1,]  0.3318672
##  [2,]  0.4315091
##  [3,]  1.2007107
##  [4,]  0.9837564
##  [5,]  0.2076259
##  [6,]  0.2638626
##  [7,] -0.3952338
##  [8,] -0.2176548
##  [9,] -1.0970325
## [10,]  0.1513780
## [11,]  0.4995639
## [12,] -0.3672771
## [13,] -1.8706458
## [14,] -1.7298169
## [15,] -1.2689705
## [16,] -0.2669823
## [17,] -0.3376925
## [18,]  0.7135992
## [19,] -1.0505316
## [20,] -0.7729587
## [21,] -1.4267870
## [22,] -0.3319350
## [23,]  1.8115447
## [24,]  0.7333263
## [25,]  1.0694917
## [26,]  1.4662866
## [27,]  1.2689963
## attr(,"scaled:center")
## [1] 29.52767
## attr(,"scaled:scale")
## [1] 1.013532

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
## tipical 
##      27

2.2.4 Holocelulose - outlier

y <- tabela3$holocelulose
scale(y) # padronização da variável 

##               [,1]
##  [1,] -0.927465940
##  [2,] -0.854870991
##  [3,] -1.324430479
##  [4,] -0.025999083
##  [5,]  0.370436514
##  [6,]  0.628920745
##  [7,]  0.823805163
##  [8,]  0.647056570
##  [9,]  1.479993882
## [10,]  0.169479455
## [11,] -0.199595501
## [12,]  0.355428315
## [13,]  0.990075132
## [14,]  1.270038207
## [15,]  0.952210656
## [16,]  0.324379028
## [17,]  0.295852212
## [18,]  0.007684768
## [19,]  0.879116955
## [20,]  0.623510368
## [21,]  1.089592398
## [22,] -0.268789189
## [23,] -2.577357669
## [24,] -1.536648217
## [25,] -0.777213098
## [26,] -1.196115607
## [27,] -1.219094595
## attr(,"scaled:center")
## [1] 68.47867
## attr(,"scaled:scale")
## [1] 1.38672

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      26

2.3 Gráfico Box Plot

par(mfrow=c(1,4))
boxplot(tabela3$extrativos, col="blue", main="Box Plot", 
        ylab= "Extrativos", adj=0.5, ylim=c(0,4))

boxplot(tabela3$lignina, col="purple",main="Box Plot", 
        ylab= "Lignina", adj=0.5, ylim=c(26,32))

boxplot(tabela3$holocelulose,col="green2",main="Box Plot", ylab= "Holocelulose", 
        adj=0.5, ylim=c(64,72)) 

boxplot(tabela3$PCS,col="gray",main="Box Plot", ylab= "PCS (MJ/kg)", 
        adj=0.5, ylim=c(17.0,17.8)) 

2.4 Dispersão, Histograma e Correlação

chart.Correlation(tabela3, histogram =TRUE)

2.5 Regressão Linear Múltipla

# Ajuste do modelo 
modelo3<-lm(formula =PCS~-1+extrativos+lignina+holocelulose, data=tabela3) 
modelo3

## 
## Call:
## lm(formula = PCS ~ -1 + extrativos + lignina + holocelulose, 
##     data = tabela3)
## 
## Coefficients:
##   extrativos       lignina  holocelulose  
##       0.2244        0.2057        0.1603

# Estatísticas das estimativas dos paramentros do modelo
summary(modelo3)

## 
## Call:
## lm(formula = PCS ~ -1 + extrativos + lignina + holocelulose, 
##     data = tabela3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.40514 -0.09816  0.00261  0.10095  0.34158 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## extrativos   0.224352   0.050789   4.417 0.000183 ***
## lignina      0.205672   0.024158   8.514 1.03e-08 ***
## holocelulose 0.160314   0.009798  16.361 1.61e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1687 on 24 degrees of freedom
## Multiple R-squared:  0.9999, Adjusted R-squared:  0.9999 
## F-statistic: 9.678e+04 on 3 and 24 DF,  p-value: < 2.2e-16

# Anova da regressão do modelo
anova(modelo3)

## Analysis of Variance Table
## 
## Response: PCS
##              Df Sum Sq Mean Sq   F value    Pr(>F)    
## extrativos    1 7376.7  7376.7 259062.87 < 2.2e-16 ***
## lignina       1  883.0   883.0  31011.91 < 2.2e-16 ***
## holocelulose  1    7.6     7.6    267.69 1.614e-14 ***
## Residuals    24    0.7     0.0                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

PCS3_est<- predict(modelo3)
resid<-(tabela3$PCS-PCS3_est)
Res2<-resid^2
SQResiduo<- sum(Res2)

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

Syx_Porc<-(Syx/(mean(tabela3$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.17	0.96

2.5.2 Predição do PCS aplicando o modelo

newdatas<-data.frame(extrativos= c(1.25),lignina= c(30), holocelulose= c(68))
PCS_estimado1 = predict(modelo3, newdatas, interval = "confidence");
PCS_estimado1 # Valores do limite de confiança dos valores preditos

##        fit      lwr      upr
## 1 17.35194 17.23403 17.46985

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(modelo3$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  modelo3$residuals
## W = 0.98761, p-value = 0.981

ks.test(modelo3$residuals, "pnorm")

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  modelo3$residuals
## D = 0.37484, p-value = 0.0006472
## 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(modelo3)

##  lag Autocorrelation D-W Statistic p-value
##    1       0.4867611      1.023195       0
##  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(modelo3)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo3
## BP = 1.326, df = 2, p-value = 0.5153

ncvTest(modelo3)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 1.173309, Df = 1, p = 0.27872

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

Valores em vermelho estimados; valores em azul observados.

PCS3_est<-predict.lm(modelo3)
newdata3 = data.frame(tabela3, PCS3_est)
Res<-(((newdata3$PCS-newdata3$PCS3_est)/(newdata3$PCS))*100)
newdata3 = data.frame(tabela3, PCS3_est,Res)

G1<-ggplot(newdata3) + 
  geom_point(size=2, aes(x=PCS, y=PCS3_est), color="red2")+
  geom_abline(intercept = 0, slope = 1, color = "black")+
  geom_point(size=2, aes(x=PCS3_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(17.0, 18.0, 0.20), limits= c(17.0, 18.0))+
  scale_y_continuous(breaks = seq(17.0, 18.0, 0.20), limits= c(17.0, 18.0))+
  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(newdata3) + 
  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(17.0, 18.0, 0.20), limits= c(17.0, 18.0)) + 
  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.