Poder Calorífico Superior em função da Composição Imediata
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
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.
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")| 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 |
# 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)| 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 |
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
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
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
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
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) chart.Correlation(tabela4, histogram =TRUE)# 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
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)| S\(_{yx}\) | S\(_{yx}\)% |
|---|---|
| 0.87 | 4.39 |
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
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
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) 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.