Poder Calorífico Superior em função da Composição Elementar
Poder Calorífico Superior em função da Composição Elementar
“Elaborado por Profª Drª Isabel Homczinski”
“04 agosto, 2022”
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
## TRUE2 Poder Calorífico Superior e a Composição Elementar
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 Elementar da biomassa (teor de carbono, hidrogênio e oxigênio). 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
tabela2 <- read_excel("Dados.xlsx",
sheet = "Analise_elementar")tabela2 %>%
kbl(caption = "Tabela dos Dados", align = 'cc') %>%
kable_classic(full_width = F, html_font = "Times", font_size = 12)| C | H | O | PCS |
|---|---|---|---|
| 83.67 | 3.56 | 2.84 | 32.856 |
| 82.62 | 3.02 | 3.66 | 33.000 |
| 63.89 | 4.97 | 24.54 | 25.100 |
| 92.04 | 2.45 | 2.96 | 34.388 |
| 89.13 | 0.43 | 0.98 | 31.124 |
| 46.90 | 6.07 | 43.99 | 18.261 |
| 49.14 | 6.34 | 43.52 | 19.423 |
| 54.41 | 4.99 | 39.69 | 21.010 |
| 48.79 | 5.91 | 43.41 | 19.260 |
| 46.58 | 5.87 | 45.46 | 18.770 |
| 48.10 | 5.99 | 45.74 | 19.916 |
| 47.84 | 5.80 | 45.76 | 18.981 |
| 50.64 | 5.98 | 42.88 | 20.720 |
| 48.15 | 5.87 | 44.75 | 19.777 |
| 46.04 | 5.82 | 44.49 | 18.640 |
2.2 Análise Descritiva
# Poder Calorífico Superior
y <- tabela2$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
y <- tabela2$C
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)
Cb<-data.frame(Med_e_sd, Coef_de_Variacao)
#Hidrogenio
y <- tabela2$H
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)
Hd<-data.frame(Med_e_sd, Coef_de_Variacao)
#Oxigenio
y <- tabela2$O
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)
Ox<-data.frame(Med_e_sd, Coef_de_Variacao)
# Junção dos data.frame
dados<-rbind(PCS,Cb)
dados<-rbind(dados, Hd)
dados<-rbind(dados, Ox)
# Criar uma coluna com o nome das variáveis
dados<-mutate(dados,
Variaveis=c("Poder Calorífico Superior (MJ/kg)",
"Carbono (%)", "Hidrogênio (%)", "Oxigênio (%)"))
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) | 23.42 \(\pm\) 6.13 | 26.18 |
| Carbono (%) | 59.86 \(\pm\) 17.53 | 29.28 |
| Hidrogênio (%) | 4.87 \(\pm\) 1.73 | 35.42 |
| Oxigênio (%) | 31.64 \(\pm\) 18.85 | 59.57 |
2.2.1 Poder Calorífico Superior - outlier
y <- tabela2$PCS
scale(y) # padronização da variável ## [,1]
## [1,] 1.5402241
## [2,] 1.5637167
## [3,] 0.2748854
## [4,] 1.7901595
## [5,] 1.2576601
## [6,] -0.8408509
## [7,] -0.6512785
## [8,] -0.3923703
## [9,] -0.6778709
## [10,] -0.7578110
## [11,] -0.5708489
## [12,] -0.7233878
## [13,] -0.4396818
## [14,] -0.5935258
## [15,] -0.7790197
## attr(,"scaled:center")
## [1] 23.41507
## attr(,"scaled:scale")
## [1] 6.129584score_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
## 152.2.2 Carbono - outlier
y <- tabela2$C
scale(y) # padronização da variável ## [,1]
## [1,] 1.3583062
## [2,] 1.2983994
## [3,] 0.2297759
## [4,] 1.8358492
## [5,] 1.6698217
## [6,] -0.7395734
## [7,] -0.6117722
## [8,] -0.3110971
## [9,] -0.6317412
## [10,] -0.7578307
## [11,] -0.6711085
## [12,] -0.6859426
## [13,] -0.5261911
## [14,] -0.6682558
## [15,] -0.7886400
## attr(,"scaled:center")
## [1] 59.86267
## attr(,"scaled:scale")
## [1] 17.52722score_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
## 152.2.3 Hidrogênio - outlier
y <- tabela2$H
scale(y) # padronização da variável ## [,1]
## [1,] -0.76009633
## [2,] -1.07309990
## [3,] 0.05719078
## [4,] -1.40349256
## [5,] -2.57435777
## [6,] 0.69479064
## [7,] 0.85129243
## [8,] 0.06878350
## [9,] 0.60204884
## [10,] 0.57886340
## [11,] 0.64841974
## [12,] 0.53828886
## [13,] 0.64262338
## [14,] 0.57886340
## [15,] 0.54988158
## attr(,"scaled:center")
## [1] 4.871333
## attr(,"scaled:scale")
## [1] 1.72522score_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 142.2.4 Oxigênio - outlier
y <- tabela2$O
scale(y) # padronização da variável ## [,1]
## [1,] -1.5280538
## [2,] -1.4845538
## [3,] -0.3768942
## [4,] -1.5216880
## [5,] -1.6267246
## [6,] 0.6549055
## [7,] 0.6299725
## [8,] 0.4267955
## [9,] 0.6241371
## [10,] 0.7328872
## [11,] 0.7477409
## [12,] 0.7488019
## [13,] 0.5960213
## [14,] 0.6952226
## [15,] 0.6814299
## attr(,"scaled:center")
## [1] 31.64467
## attr(,"scaled:scale")
## [1] 18.85056score_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
## 152.3 Gráfico Box Plot
par(mfrow=c(1,4))
boxplot(tabela2$C, col="blue", main="Box Plot",
ylab= "Carbono (%)", adj=0.5, ylim=c(40,100))
boxplot(tabela2$H, col="purple",main="Box Plot",
ylab= "Hidrogênio (%)", adj=0.5, ylim=c(2,7))
boxplot(tabela2$O,col="green2",main="Box Plot", ylab= "Oxigênio (%)",
adj=0.5, ylim=c(0,50))
boxplot(tabela2$PCS,col="gray",main="Box Plot", ylab= "PCS (MJ/kg)",
adj=0.5, ylim=c(15,35)) 
2.4 Dispersão, Histograma e Correlação
chart.Correlation(tabela2, histogram =TRUE)
2.5 Regressão Linear Múltipla
# Ajuste do modelo
modelo2<-lm(formula =PCS~-1+C+H+O, data=tabela2); modelo2##
## Call:
## lm(formula = PCS ~ -1 + C + H + O, data = tabela2)
##
## Coefficients:
## C H O
## 0.34876 1.14845 -0.09662# Estatísticas das estimativas dos paramentros do modelo
summary(modelo2)##
## Call:
## lm(formula = PCS ~ -1 + C + H + O, data = tabela2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.81639 -0.35435 0.05688 0.26616 1.07100
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## C 0.348757 0.006411 54.401 9.8e-16 ***
## H 1.148451 0.224550 5.114 0.000256 ***
## O -0.096622 0.024853 -3.888 0.002157 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5792 on 12 degrees of freedom
## Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994
## F-statistic: 8690 on 3 and 12 DF, p-value: < 2.2e-16# Anova da regressão do modelo
anova(modelo2)## Analysis of Variance Table
##
## Response: PCS
## Df Sum Sq Mean Sq F value Pr(>F)
## C 1 8729.7 8729.7 26022.698 < 2.2e-16 ***
## H 1 11.2 11.2 33.345 8.818e-05 ***
## O 1 5.1 5.1 15.115 0.002157 **
## Residuals 12 4.0 0.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.5.1 Predição do PCS aplicando o modelo
newdatas<-data.frame(C= c(83),H= c(3), O= c(3))
PCS_estimado2 = predict(modelo2, newdatas, interval = "confidence");
PCS_estimado2 # Valores do limite de confiança dos valores preditos## fit lwr upr
## 1 32.10232 31.43632 32.768332.5.2 Medidas de Precisão para Análise do Modelo
PCS2_est<-predict.lm(modelo2)
newdata1 = data.frame(tabela2, PCS2_est)
Res<-(((newdata1$PCS-newdata1$PCS2_est)/(newdata1$PCS))*100)
newdata1 = data.frame(tabela2, PCS2_est,Res)
resid<-(newdata1$PCS-newdata1$PCS2_est)
Res2<-resid^2
SQResiduo<- sum(Res2)
SQTotal<-var(newdata1$PCS)*((length(newdata1$PCS))-1)
R2<-(1-(SQResiduo/SQTotal))
a<-(length(newdata1$PCS)-1)/(length(newdata1$PCS)-2)
R2adj<-(1-(a*(1-R2)))
Syx<-sqrt(SQResiduo/((length(newdata1$PCS))-2))
Syx_Porc<-(Syx/(mean(newdata1$PCS)))*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)| R\(^{2}\) | R\(_{adj}^{2}\) | S\(_{yx}\) | S\(_{yx}\)% |
|---|---|---|---|
| 0.99 | 0.99 | 0.56 | 2.38 |
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(modelo2$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.97638, p-value = 0.9387ks.test(modelo2$residuals, "pnorm")##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: modelo2$residuals
## D = 0.20714, p-value = 0.4781
## alternative hypothesis: two-sidedINDEPENDÊ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(modelo2)## lag Autocorrelation D-W Statistic p-value
## 1 0.1754137 1.634667 0.256
## Alternative hypothesis: rho != 0HOMOSCEDASTICIDADE 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(modelo2)##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 1.6078, df = 2, p-value = 0.4476ncvTest(modelo2)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.1071283, Df = 1, p = 0.743442.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=PCS, y=PCS2_est), color="red2")+
geom_abline(intercept = 0, slope = 1, color = "black")+
geom_point(size=2, aes(x=PCS2_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(10, 40, 5), limits= c(10,40))+
scale_y_continuous(breaks = seq(10, 40, 5), limits= c(10,40))+
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=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(10, 40, 5), limits= c(10,40)) +
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.